SUBROUTINE ADVNCE C **** ADVANCES THE MODEL IN TIME include "params.h" include "blnk.h" include "vscr.h" include "buff.h" include "cons.h" include "index.h" include "strt.h" include "unit.h" include "phys.h" include "sechis.h" C **** C **** DIMENSIONS FOR AMIE FIELDS C **** include "amie.h" real :: wall_time0,wall_time1 ! write(6,"('ADVANCE entered at iter=',i6)") iter C **** INITIALIZE BUFFER INDICES MODULO=8 CALL INIT C **** TIME LOOP TIME = SECOND(X) TIMET = TIMEF(X) 100 CONTINUE ITER=ITER+1 call wtimer(wall_time0) write(6,"(/,'advnce begin iter ',i6,' at wall_time=',f10.3, | ' (secs)')") iter,wall_time0/1000. JTOTS=0. JTOTN=0. C **** C **** FIND UT TIME BETWEEN 0-24 HR (MODEL DAY = UT DAY, 11/89) C **** SECS = AMOD(FLOAT(ITER)*C(4),86400.) IF (IADVDA .EQ. 1) THEN IDAYIT = ITER*IFIX(C(4))/86400 IF (IDAYIT*86400 .EQ. ITER*IFIX(C(4))) THEN C **** C **** ADVANCE DAY C **** IDAPR = IIDAY IYRPR = IYEAR IIDAY = IIDAY + 1 C **** C **** LPYR = 1(0) IF IS (NOT) A LEAP YEAR C **** IYR4 = IYEAR/4 IYR100 = IYEAR/100 LPYR = 0 IF (IYR4*4 .EQ. IYEAR .AND. IYR100*100 .NE. IYEAR) LPYR=1 IENDA = 365 + LPYR IF (IIDAY .GT. IENDA) THEN IYEAR = IYEAR + 1 IIDAY = IIDAY - IENDA ENDIF ! FOR PAST YEAR'S END C **** C **** RECALCULATE SUN'S DECLINATION C **** PI = 3.14159265358979 DELTA = ATAN(TAN(23.5*PI/180.)*SIN(2.*PI*FLOAT(IIDAY-80)/365.)) C(95) = SIN(DELTA) C(96) = COS(DELTA) WRITE(6,"(1X,'ADVNCE: ADVANCING DAY (PREVIOUS,PRESENT)=',4I5)") 1 IDAPR,IYRPR,IIDAY,IYEAR ENDIF ! FOR ADVANCING DAY ENDIF ! FOR IADVDA FLAG ON C **** C **** CALCULATE AURORAL PARAMETERS FOR THIS TIME STEP C **** C **** FIND UT TIME (MODEL DAY = UT DAY, 11/89) SEC = FLOAT(ITER)*C(4) SECTGCM = FLOAT(ITER)*C(4) CALL SUN IF (IAMIE.EQ.1) THEN IF (NTST.LT.NTIMS .AND. SEC.GT.SECUTA(2)) THEN NTST = NTST + 1 SECUTA(1) = SECUTA(2) HPNHA(1) = HPNHA(2) HPSHA(1) = HPSHA(2) CPNHA(1) = CPNHA(2) CPSHA(1) = CPSHA(2) CUSPSTA(1) = CUSPSTA(2) CUSPSLA(1) = CUSPSLA(2) CUSPNTA(1) = CUSPNTA(2) CUSPNLA(1) = CUSPNLA(2) DO 120 JA=1,36 DO 120 IA=1,73 POTKVA1(IA,JA) = POTKVA2(IA,JA) EKEVA1(IA,JA) = EKEVA2(IA,JA) EFLXA1(IA,JA) = EFLXA2(IA,JA) 120 CONTINUE DO 121 KA=1,3 DO 121 JA=1,36 DO 121 IA=1,73 VIXYZA1(IA,JA,KA) = VIXYZA2(IA,JA,KA) 121 CONTINUE READ (IUAMI) SECUTA(2),HPSHA(2),HPNHA(2),CPSHA(2),CPNHA(2), 1 CUSPSTA(2),CUSPNTA(2),CUSPSLA(2),CUSPNLA(2),POTKVA2,VIXYZA2, 2 EKEVA2,EFLXA2 ENDIF ! FOR READING AMIE ENDIF ! FOR IAMIE=1 call getgpi CALL TAIL (SEC) IHIS=IHIS-1 ISAV=ISAV-1 if (iter.le.isecstop) then ihissech = ihissech-1 isavsech = isavsech-1 endif C **** FOR THE FIRST FORWARD TIME STEP OF A C **** FORECAST START, HALVE DT. IF(MODEFC)150,130,140 130 C(4) = C(111)/2. C(6) = 2.*C(4) C(7) = 1./C(6) MODEFC=1 GO TO 150 C **** RESTORE DT TO ORIGINAL VALUE 140 C(4) = C(111) C(6) = 2.*C(4) C(7) = 1./C(6) MODEFC=-1 150 CONTINUE C **** BEGIN SWEEP OF MESH FROM SOUTH TO NORTH C **** BRING IN LINE 1 CALL VDRIFT CALL FLIP CALL UNLCRD(IUNIN,NJIN,LBDSKI,NSIN) CALL UNLCCK(IUNIN,NSIN,LEN) IF(LEN.EQ.LBDSKI.AND.NSIN.EQ.0) GO TO 155 NUM=155 GO TO 4000 C **** BRING IN LINES 2 AND 3 C **** AND REARRANGE LINES 1 AND 2 C **** TO FORM LINES 0 AND -1 155 DO 270 JP2=1,2 CALL FLIP CALL UNLCRD(IUNIN,NJIN,LBDSKI,NSIN) C **** ADD DIAGNOSTIC QUANTITIES ONTO LINE JP2 J=JP2-2 CALL ADDIAG C **** REARRANGE LINE JP2 TO GET LINE -JP2+1 NJPUT=NJP1 IF(JP2.EQ.2) NJPUT=NJM1 IF(NHEMI.NE.0) GO TO 200 C **** CONSTRUCT LINES -1,0 FOR GLOBAL CASE ILIM=IMAXH+2 DO 195 N=1,NFLDS C **** CHECK IF FIELD IS NEEDED ACROSS POLE IF(KPOLE(N)) 160,195,160 C **** TRANSFER Z TYPES FIRST ITERATION ONLY 160 IF(IFRST.EQ.1 .AND. KPOLE(N).EQ.2) GO TO 195 C **** DETERMINE LENGTH OF FIELD IF(KFLDS(N)-2) 165,167,170 165 KLIM=1 LEN=LEN1 GO TO 175 167 KLIM=KMAX LEN=LEN2 GO TO 175 170 KLIM=KMAXP1 LEN=LEN3 C **** SET POINTERS TO INPUT AND OUTPUT BUFFERS 175 NXK = NJP2+NDEXA(N+1) NYK = NJPUT+NDEXA(N+1) C **** CHECK IF FIELD IS TO HAVE SIGN CHANGED IF(KPOLE(N)) 180,195,185 C **** CHANGE SIGN 180 DO 184 K=1,KLIM ID=IMAXH+1 DO 182 I=1,ILIM F(ID+2,NYK)=-F(I+2,NXK) F(I ,NYK)=-F(ID ,NXK) ID=ID+1 182 CONTINUE NXK=NXK+1 NYK=NYK+1 184 CONTINUE GO TO 195 C **** DO NOT CHANGE SIGN 185 DO 188 K=1,KLIM ID=IMAXH+1 DO 187 I=1,ILIM F(ID+2,NYK)=F(I+2,NXK) F(I ,NYK)=F(ID ,NXK) ID=ID+1 187 CONTINUE NXK=NXK+1 NYK=NYK+1 188 CONTINUE 195 CONTINUE GO TO 267 C **** CONSTRUCT LINES -1,0 FOR HEMISPHERIC CASE 200 DO 265 N=1,NFLDS C **** CHECK IF FIELD IS NEEDED ACROSS EQUATOR IF(KEQTR(N)) 205,265,205 C **** TRANSFER Z TYPES FIRST ITERATION ONLY 205 IF(IFRST.EQ.1 .AND. KEQTR(N).EQ.2) GO TO 265 C **** DETERMINE LENGTH OF FIELD IF(KFLDS(N)-2) 210,215,220 210 KLIM=1 LEN=LEN1 GO TO 225 215 KLIM=KMAX LEN=LEN2 GO TO 225 220 KLIM=KMAXP1 LEN=LEN3 C **** SET POINTERS TO FIELDS IN SCM BUFFERS 225 NJP2K = NJP2+NDEXA(N+1) NJPUTK= NJPUT+NDEXA(N+1) C **** CHECK IF FIELD IS TO HAVE SIGN CHANGED IF(KEQTR(N)) 230,265,260 C **** CHANGE SIGN AND STORE IN PROPER SCM BUFFER 230 DO 235 I=1,LEN 235 F(I,NJPUTK)=-1.*F(I,NJP2K) GO TO 265 C **** DO NOT CHANGE SIGN AND STORE IN PROPER C **** SCM BUFFER. 260 DO 263 I=1,LEN 263 F(I,NJPUTK)=F(I,NJP2K) 265 CONTINUE 267 CONTINUE C **** CHECK FOR COMPLETION OF BUFFERING CALL UNLCCK(IUNIN,NSIN,LEN) IF(LEN.EQ.LBDSKI.AND.NSIN.EQ.0) GO TO 270 NUM=270 GO TO 4000 270 CONTINUE C **** CHANGE BUFFER POINTERS SO THAT LINE J IS -2 C **** SO THAT INITIALIZATION CALLS TO HDIF C **** CAN BE MADE. DO 280 J=1,5 280 CALL FLIP C **** ITITIALIZE HORIZONTAL DIFFUSION QUANTITIES DO 290 JJ=1,3 CALL FLIP J=JJ-3 290 CALL HDIF C **** GENERAL LINE LOOP (LATITUDINAL VARIATION) DO 3000 J=1,JMAX write(6,"(' iter=',i6,' j=',i3)") iter,j RACS=1./(C(51)*CS(J)) CALL FLIP C **** START IN NEXT LINE NORTH IF(J.GE.JMAXM2) GO TO 310 CALL UNLCRD(IUNIN,NJIN,LBDSKI,NSIN) C **** START OUT PREVIOUS LINE SOUTH 310 IF(J.EQ.1) GO TO 315 CALL UNLCWT(IUNOT,NJOT,LBDSKI,NSOT) C **** ADD DIAGNOSTIC QUANTITIES ONTO LINE J+2 315 CONTINUE IF(J.GT.JMAXM2)GO TO 405 CALL ADDIAG 405 CONTINUE IF(J.EQ.JMAX) IFRST=1 C ****************************************************************** C **** C **** DYNAMICS C **** C ****************************************************************** C **** C **** CALL TO BNDCMP SHOULD LEAD CALLS IN DYNAMICS SECTION OF C **** ADVNCE TO ENSURE THAT MATRICES B(ZIMXP,2,2) AND VECTORS C **** FB(ZIMXP,2) ARE AVAILABLE WHEN NEEDED C **** CALL BNDCMP C **** C **** CALCULATE W AT N*DT C **** CALL SWDOT C **** C **** CALCULATE ION DRIFT VELOCITY C **** CALL VDRIFT2 C **** CALCULATE CP, KT, AND KM CALL CPKTKM C **** CLEAR NQO2P, NQOP, NQN2P, NQNOP, NQNP NQPK=NQO2P DO 1 I=1,5*LEN3 F(I,NQPK)=0. 1 CONTINUE C **** CALCULATE COLUMN NUMBER DENSITIES CALL CHAPMN C **** CALCULATE SOLAR HEATING C **** CALCULATE REACTION RATES CALL ALTV CALL RATES CALL QRJ C **** CALCULATE BACKGROUND (NIGHT-TIME) IONIZATION CALL QINITE C **** ADD IONIZATION DUE TO SOLAR X-RAYS CALL XRAY C **** GENERATE AURORAL FIELDS, UI, VI, WI. FIELDS RETURNED C **** IN T1 THRU T7 CALL HEELIS C **** CALCULATE AURORAL ADDITIONS TO IONIZATION RATES CALL ORORA C **** SOLVE FOR N(O+) C **** SET FLUX AT UPPER BOUNDARY CALL OPFLUX CALL OPLUS C **** COMPUTE ELECTRON DENSITY, N(E), AND ION NUMBER DENSITIES, C **** N(O2+), N(N2+), N(NO+), N(N+) CALL ELDEN C **** ION-DRAG PARAMETERS, LXX, LYY, LXY (IN DIPOLE COORDINATE C **** SYSTEM) CALL LAMDAS C **** ROTATE LAMDAS TO GEOGRAPHIC SYSTEM CALL ROTATE C **** ELECTRON PRODUCTION OF N4S AND N2D CALL QTIEFF C **** ODD NITROGEN CHEMISTRY CALL CMPN2D CALL CMPN4S CALL CMPNO CALL CMPO2O C **** HEATING AND O2 DISSOCIATION DUE TO ODD NITROGEN AND ION C **** CHEMISTRY CALL QJNNO CALL QJION C **** CALCULATE ELECTRON AND ION TEMPERATURES CALL SETTEI C **** CALCULATE IMPLICIT AND EXPLICIT COOLING TERMS IN NWTI C **** AND NWTE CALL NEWTON C **** ADD IN O3P COOLING CALL NEWTO3P C **** CALCULATE FLF, FPH, NQDH, NKMH, NPSDH, NPSDH2 C **** AT (N-1)*DT CALL HDIF C **** ADVANCE T, U, V BY ONE TIME STEP CALL DT CALL DUV C **** COMPOSITION OF MAJOR SPECIES CALL COMP C ****************************************************************** C **** C **** END OF DYNAMICS C **** C ****************************************************************** C **** CHECK IF ON LINE JMAXM2 OR JMAXM1 IF(J.LT.JMAXM2) GO TO 2980 IF(J.EQ.JMAX) GO TO 2990 C **** REARRANGE LINE NJGET TO GET LINE NJIN NJGET=NJP2 IF(J.EQ.JMAXM1) NJGET=NJ ILIM=IMAXH+2 DO 2975 N=1,NFLDS C **** CHECK IF FIELD IS NEEDED ACROSS POLE IF(KPOLE(N)) 2916,2975,2916 C **** TRANSFER Z TYPES FIRST ITERATION ONLY 2916 IF(IFRST.EQ.1.AND.KPOLE(N).EQ.2) GO TO 2975 C **** DETERMINE LENGTH OF FIELD IF(KFLDS(N)-2) 2920,2925,2930 2920 KLIM=1 LEN=LEN1 GO TO 2935 2925 KLIM=KMAX LEN=LEN2 GO TO 2935 2930 KLIM=KMAXP1 LEN=LEN3 C **** BRING IN FIELD FROM SCM BUFFER 2935 NXK = NJGET+NDEXA(N+1) NYK = NJIN+NDEXA(N+1) C **** CHECK IF FIELD IS TO HAVE SIGN CHANGED IF(KPOLE(N)) 2940,2975,2955 C **** CHANGE SIGN 2940 DO 2950 K=1,KLIM ID=IMAXH+1 DO 2945 I=1,ILIM F(ID+2,NYK)=-F(I+2,NXK) F(I ,NYK)=-F(ID ,NXK) ID=ID+1 2945 CONTINUE NXK=NXK+1 NYK=NYK+1 2950 CONTINUE GO TO 2975 C **** DO NOT CHANGE SIGN 2955 DO 2965 K=1,KLIM ID=IMAXH+1 DO 2960 I=1,ILIM F(ID+2,NYK)=F(I+2,NXK) F(I ,NYK)=F(ID ,NXK) ID=ID+1 2960 CONTINUE NXK=NXK+1 NYK=NYK+1 2965 CONTINUE C **** TRANSFER REARRANGED FIELD C **** TO CORRECT LOCATION IN SCM BUFFER 2975 CONTINUE GO TO 2990 C **** CHECK FOR COMPLETION OF BUFFERING 2980 CALL UNLCCK(IUNIN,NSIN,LEN) IF(LEN.EQ.LBDSKI.AND.NSIN.EQ.0) GO TO 2990 NUM=2990 GO TO 4000 2990 CONTINUE IF(J.EQ.1)GO TO 3000 CALL UNLCCK(IUNOT,NSOT,LEN) IF(LEN.EQ.LBDSKI.AND.NSOT.EQ.0) GO TO 3000 NUM=3000 GO TO 4000 3000 CONTINUE C **** TRANSFER LINE JMAX TO DISK CALL FLIP CALL UNLCWT(IUNOT,NJOT,LBDSKI,NSOT) CALL UNLCCK(IUNOT,NSOT,LEN) IF(LEN.EQ.LBDSKI.AND.NSOT.EQ.0) GO TO 3100 NUM=3100 GO TO 4000 3100 ITEMP=IUNIN IUNIN=IUNOT IUNOT=ITEMP REWIND IUNIN REWIND IUNOT c JTOTN=JTOTN*C(81)*C(55)*C(3)*C(53)*C(1)*C(2) c IF(NHEMI.EQ.0)JTOTS=JTOTS*C(81)*C(55)*C(3)*C(53)*C(1)*C(2) C **** C **** SET UP AND SOLVE DYNAMO EQUATIONS C **** CALL DYN C **** CHECK FOR EXIT C **** CHECK FOR SAVE TAPE WRITE call wtimer(wall_time1) write(6,"(' Iteration ',i6,' completed in ',f10.3,' wall clock ', | 'seconds')") iter,wall_time1/1000.-wall_time0/1000. ! c write(6,"('ADVNCE iter',i5,': ihis=',i2, c + ' isav=',i2,' sechis=',i2,' secsav=',i2,' secstrt,stop=', c + 2i5)") iter,ihis,isav,ihissech,isavsech,isecstart,isecstop 3130 if (isav.eq.0.or.isavsech.eq.0) goto 3150 if (ihis.ne.0.and.ihissech.ne.0) goto 3140 if (iter.ge.nstp) then isav = 0 isavsech = 0 endif goto 3150 3140 if (iter.lt.nstp) goto 100 isav = 0 isavsech = 0 3150 TIME = SECOND(X)-TIME TIMET = TIMEF(X)-TIMET WRITE(6,3170) TIMET 3170 format('TIME FROM TIMEF CALL = ',e12.5) WRITE(6,3160) ITER,TIME 3160 format('ADVANCE EXITED AT ITER = ',i6,' TIME IN INNER LOOP = ', + e12.5,' SECS') RETURN 4000 WRITE(6,4005) ITER 4005 FORMAT(*1ERROR EXIT AT ITER =*I6,5X,*RESTART*) WRITE(6,4010)NUM,LEN,NSIN,NSOT,NSIN1,NSIN2,J 4010 FORMAT(*0ADVANCE- I/O PROBLEM AT NUM=*I5,2X,*LEN=*I5,2X,*NSIN=*I5, A2X,*NSOT=*I5,2X,*NSIN1=*I5,2X,*NSIN2=*I5,2X,*J=*I2) CALL EXIT END C !----------------------------------------------------------------------- subroutine wtimer(t) ! ! Use fortran intrinsic date_and_time to return elapsed wall clock ! time t since this routine was last called. ! implicit none real,intent(out) :: t character (len=8) :: date character (len=10) :: time character (len=5) :: zone integer :: val1(8), val2(8), dv(8), its integer(kind=8) :: ms data its / 0 / save ! if ( its == 0 ) then call date_and_time(date,time,zone,val1) t = 0.0 its = 1 else call date_and_time(date,time,zone,val2) dv = val2-val1 ms = dv(8)+1000*dv(7)+60*1000*dv(6)+60*60*1000*dv(5) end if t = real(ms) end subroutine wtimer