Tue Mar 14 13:16:01 US/Mountain 2000 Diffsrc executing from /home/foster/tiegcm/tgcm13mt/modsrc.snoe Source directory = /home/tgcm/tgcm14a ======================================================================== Diff of /home/tgcm/tgcm14a/cons.h and cons.h: 4,8c4,7 < | ncols,nprog,nprogp,lflds,kflds,kpole,keqtr,kutsp,kut,kutnp, < | inpt < real :: c,brn2d,cssp,cs,csnp,cor,sn,tn,t0,exps,a,rmass, < | difk,tbound,zbound,dift,decayf,decayh < ! --- > | ncols,nprog,nprogp,kflds,kpole,keqtr,kutsp,kut, > | kutnp,inpt,iaur,idisp,lflds > real :: c,cssp,cs,csnp,cor,sn,tn,t0,exps,a,rmass, > | difk,tbound,dift,decayf,decayh,xmue,alfalp,efluxlp,sfeps 10,20c9,36 < | KMAX,KMAXM1,KMAXP1,LAST,LONG,LEN1,LEN2,LEN3,NCBUF,NCPHYS,NCOLS, < | NPROG,NPROGP,C(120),BRN2D, < | LFLDS(nflds),KFLDS(nflds),KPOLE(nflds),KEQTR(nflds), < | CSSP(2),CS(ZJMX),CSNP(2),COR(ZJMX),KUTSP(2),KUT(ZJMX),KUTNP(2), < | SN(ZJMX),TN(ZJMX),T0(ZKMXP),EXPS(ZKMXP),A(ZKMX), < | RMASS(3),DIFK(ZKMXP),INPT(170), < | TBOUND,ZBOUND,DIFT(ZKMXP),DECAYF(ZJMX,ZKMXP),DECAYH(ZJMX,ZKMX) < ! < ! Field information written to history summary for post-model processing. < ! < character*8 nflds_lab(nflds), nphys_lab(nphys) --- > 1KMAX,KMAXM1,KMAXP1,LAST,LONG,LEN1,LEN2,LEN3,NCBUF,NCPHYS,NCOLS, > 2NPROG,NPROGP,C(120),LFLDS(32),KFLDS(32),KPOLE(32),KEQTR(32), > 3CSSP(2),CS(ZJMX),CSNP(2),COR(ZJMX),KUTSP(2),KUT(ZJMX),KUTNP(2), > 4SN(ZJMX),TN(ZJMX),T0(ZKMXP),EXPS(ZKMX),A(ZKMX),ZB > 5(ZJMX),TB(ZJMX),UB(ZJMX),VB(ZJMX),BND(ZIMXP),CI, > 6RMASS(3),DIFK(ZKMXP),INPT(165),IAUR,IDISP, > 6ALFALP,EFLUXLP,TBOUND,DIFT(ZKMXP),XMUE(ZKMXP), > 7ZB2(ZJMX),TB2(ZJMX),UB2(ZJMX),VB2(ZJMX),BND2(ZIMXP), > 8ZBA(ZJMX),TBA(ZJMX),UBA(ZJMX),VBA(ZJMX),BNDA(ZIMXP),SFEPS > COMPLEX ZB,TB,UB,VB,BND,CI > COMPLEX ZB2,TB2,UB2,VB2,BND2 > COMPLEX ZBA,TBA,UBA,VBA,BNDA > c > c Field information written to history summary for post-model processing. > c > c nfproc is the number of fields processors should read from the history. > c nfproc is usually the number of fields up to and including NPHI > c (see index.h) > c > c The names are floats so they can be put in summary (f(1-lsumm)), but > c are transferred to a character array in the processor. > c The names are defined with data statements in types.h. > c > integer,parameter :: nfproc=16 > real fproc_names(nfproc) > common/fproc/ fproc_names > C > character*8 nflds_lab(nzflds), nphys_lab(nzphys) 22d37 < ! ======================================================================== Diff of /home/tgcm/tgcm14a/crates.h and crates.h: 39d38 < !$OMP THREADPRIVATE (/rates_priv/,/metas/,/fsaray/,/rates_priv1/) ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/othend.h ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/advnce.f ======================================================================== Diff of /home/tgcm/tgcm14a/altv.f and altv.f: 10,12d9 < ! 2/29/00: added to tgcm14: if used, this is called from < ! dynamics, before rates. tvib is referenced in rates.f. < ! 14c11 < include 'fcom.h' --- > include 'blnk.h' 144,145c141 < ! call addfsech('TVIB','VIBRATIONAL TEMPERATURE', < ! | ' ',tvib,zimxp,zkmxp,zkmx,j) --- > c call addfsech('TVIB',tvib,zimxp,zkmxp,kmax,j) ======================================================================== Diff of /home/tgcm/tgcm14a/aurht.f and aurht.f: 1c1 < SUBROUTINE AURHT(ALFA,FLUX,DRIZL,t7,t9,IM,IHEM) --- > SUBROUTINE AURHT(ALFA,FLUX,DRIZL,IM) 5,8d4 < ! < ! Defines ALFA, FLUX, DRIZL. These are passed back through flowv2, flowv1, < ! flowx0, and flowxx to heelis < ! 12,15c8,24 < include "ioncr.h" < include "ovalr.h" < include "aurlp.h" ! low-protons in aurora < real :: alfa3,flux3 --- > > ! include "phys.h" ! for debug only > > integer :: istar,ihem > real :: theta0,dduumm > COMMON /IONCR/ ISTAR,IHEM,THETA0(2), > 1 DDUUMM(IMAXMP*JMX0+IMX0*JMX0+30) > real :: rrad,h0,rh,rroth,e0,ree,rrote,fc,alfac,fd,alfad, > | alfk,alf6,alf21,rrot6,rrot21,rd6,rd6v,rd21,rh6,rh21,rt6, > | rt21,alfa0,ralfa,alfa20,ralfa2,e20,re2 > COMMON /OVALR/ RRAD(2),H0,RH,RROTH,E0,REE,RROTE,FC,ALFAC,FD,ALFAD > 1, ALFK,ALF6,ALF21,RROT6,RROT21,RD6,RD6V,RD21,RH6,RH21,RT6,RT21 > 2, ALFA0,RALFA,ALFA20,RALFA2,E20,RE2 > c > c Low energy auroral proton input (see ALFALP and EFLUXLP from input): > real :: alfa_lp,eflux_lp,qteaur > common/aurlp/ alfa_lp(zimxp),eflux_lp(zimxp),qteaur(zimxp) 18,19c27,28 < integer,intent(in) :: IHEM,im < real,intent(out) :: ALFA(IM),FLUX(IM),DRIZL(IM),t7(IM),t9(IM) --- > real,intent(out) :: ALFA(ZIMXP,2),FLUX(ZIMXP,2),DRIZL(1) > integer,intent(in) :: im 21a31 > real :: pi,alfap1,alfap2,pe1,pe2,alfap0,e0p 23,26c33,34 < real :: alfap1,alfap2,pe1,pe2,alfap0,e0p < ! < ! ALFA0 was defined by tail.f: < ! --- > C > PI = 3.1415926535898 29c37 < DO 1 I=1,IM --- > DO 1 I=1,IM 33c41 < WK1(I)=COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE))) --- > WK1(I)=COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE))) 37,38c45,46 < WK2(I)=H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH), < 1 COS(ALON(I)-RROTH))) ) --- > WK2(I)=H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH), > 1 COS(ALON(I)-RROTH))) ) 42c50 < WK3(I) = COLAT(I) - RRAD(IHEM) --- > WK3(I) = COLAT(I) - RRAD(IHEM) 45,57c53,54 < ALFA(I) = ALFA0*(1.-RALFA*WK1(I)) < < ! write(6,"(/'aurht: i=',i2,' ihem=',i2)") i,ihem < ! write(6,"(' alfa0=',e12.4,' alon(i)=',e12.4)") alfa0,alon(i) < ! write(6,"(' rrote=',e12.4,' h0=',e12.4,' rh=',e12.4)") < ! | rrote,h0,rh < ! write(6,"(' rroth=',e12.4,' rrad(ihem)=',e12.4)") < ! | rroth,rrad(ihem) < ! write(6,"(' ralfa=',e12.4)") ralfa < < 1 CONTINUE < ! < ! ALFA0 <= 0.01: --- > ALFA(I,1) = ALFA0*(1.-RALFA*WK1(I)) > 1 CONTINUE 62c59 < DO 2 I=1,IM --- > DO 2 I=1,IM 66c63 < WK1(I)=COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE))) --- > WK1(I)=COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE))) 70,71c67,68 < WK2(I)=H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH), < 1 COS(ALON(I)-RROTH))) ) --- > WK2(I)=H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH), > 1 COS(ALON(I)-RROTH))) ) 75c72 < WK3(I) = COLAT(I) - RRAD(IHEM) --- > WK3(I) = COLAT(I) - RRAD(IHEM) 79,83c76,81 < C ALFA(I,1) = ALFK+ALF1*EXP(-((WK3(I)-(RD1A+RD1V*COS(ALON(I)))) < C 1 / RAH1)**2)*EXP(-( (ABS(ALON(I)-RAROT1)-IFIX(ABS(ALON(I)- < C 2 RAROT1)/180.)*360.)/RAHT1)**2)+ALF2 * EXP(-( (WK3(I)-RD2) / < C 3 RAH2)**2)*EXP(-( (ABS(ALON(I)-RAROT2)-IFIX(ABS(ALON(I)-RAROT2) < C 4 /180.)*360.) / RAHT2)**2) --- > C ALFA(I,1) = ALFK+ALF1*EXP(-((WK3(I)-(RD1A+RD1V*COS(ALON(I)))) > C 1 / RAH1)**2) > C 2 *EXP(-( (ABS(ALON(I)-RAROT1)-IFIX(ABS(ALON(I)-RAROT1)/180.)*360.) > C 3 /RAHT1)**2) + ALF2 * EXP(-( (WK3(I)-RD2) / RAH2)**2) > C 4 *EXP(-( (ABS(ALON(I)-RAROT2)-IFIX(ABS(ALON(I)-RAROT2)/180.)*360.) > C 5 / RAHT2)**2) 89,95c87,92 < ALFA(I) = ALFK < 1 + ALF6 * EXP(-( (WK3(I)-RD6-RD6V*COS(ALON(I)) )/RH6)**2) * < 2 EXP(-( ATAN2(SIN(ALON(I)-RROT6),COS(ALON(I)-RROT6))/RT6 )**2) < 3 + ALF21 * EXP(-( (WK3(I)-RD21)/RH21)**2) * < 4 EXP(-( ATAN2(SIN(ALON(I)-RROT21),COS(ALON(I)-RROT21))/ < 5 RT21)**2) < 2 CONTINUE --- > ALFA(I,1) = ALFK > 1 + ALF6 * EXP(-( (WK3(I)-RD6-RD6V*COS(ALON(I)) )/RH6)**2) * > 2 EXP(-( ATAN2(SIN(ALON(I)-RROT6),COS(ALON(I)-RROT6))/RT6 )**2 ) > 3 + ALF21 * EXP(-( (WK3(I)-RD21)/RH21)**2) * > 4 EXP(-( ATAN2(SIN(ALON(I)-RROT21),COS(ALON(I)-RROT21))/RT21 )**2 ) > 2 CONTINUE 97,104c94,104 < ! < DO I=1,IM < FLUX(I)=E0*(1.-REE*WK1(I))*EXP(-( WK3(I) /WK2(I))**2) < | / (2. * ALFA(I) * 1.602E-9) < DRIZL(I)=EXP(-(( WK3(I) +ABS( WK3(I) ))/(2.*H0))**2) < t7(i) = ALFA20*(1.-RALFA2*WK1(I)) < t9(i) = E20*(1.-RE2*WK1(I))*EXP(-( WK3(I) /WK2(I))**2) < | / (2. * t7(i) * 1.602E-9) --- > > ! write(6,"('aurht: j=',i3,' e0=',e12.4)") j,e0 > > DO 3 I=1,IM > FLUX(I,1)=E0*(1.-REE*WK1(I))*EXP(-( WK3(I) /WK2(I))**2) > 1 / (2. * ALFA(I,1) * 1.602E-9) > DRIZL(I)=EXP(-(( WK3(I) +ABS( WK3(I) ))/(2.*H0))**2) > ALFA(I,2) = ALFA20*(1.-RALFA2*WK1(I)) > FLUX(I,2)=E20*(1.-RE2*WK1(I))*EXP(-( WK3(I) /WK2(I))**2) > | / (2. * ALFA(I,2) * 1.602E-9) > C The above line was commented out until 8/92, effectively making FLUX2 107,118c107,116 < ALFAP1 = 5. < ALFAP2 = 15. < C PE1 = 0.15 < C PE2 = 0.4 < PE1 = 1.E-20 < PE2 = 1.E-20 < ALFAP0 = 0.5*(ALFAP1+ALFAP2) < alfa_lp(i+2) = ALFAP0*(1-RALFA2*WK1(I)) < E0P = 0.5*(PE1+PE2) < eflux_lp(i+2) = E0P*(1.-REE*WK1(I))*exp(-(wk3(i)/wk2(i))**2) < qteaur(i+2) = -7.e+8*exp(-(wk3(i)/wk2(i))**2) < enddo --- > ALFAP1 = 5. > ALFAP2 = 15. > C PE1 = 0.15 > C PE2 = 0.4 > PE1 = 1.E-20 > PE2 = 1.E-20 > ALFAP0 = 0.5*(ALFAP1+ALFAP2) > alfa_lp(i+2) = ALFAP0*(1-RALFA2*WK1(I)) > E0P = 0.5*(PE1+PE2) > eflux_lp(i+2) = E0P*(1.-REE*WK1(I))*exp(-(wk3(i)/wk2(i))**2) 119a118,120 > qteaur(i+2) = -7.e+8*exp(-(wk3(i)/wk2(i))**2) > 3 CONTINUE > c ======================================================================== Diff of /home/tgcm/tgcm14a/cmpn2d.f and cmpn2d.f: 8c8 < include "fcom.h" --- > include "blnk.h" 14c14,15 < include "mwt.h" --- > real :: rmn4s,rmn2d,rmno,brn2d,cee > COMMON/MASS/ RMN4S,RMN2D,RMNO,BRN2D,CEE 45,46c46,47 < 1(RA1(I,1)*F(I,NNOPK)*0.85+RA3(I,1)*F(I,NN2PK)*0.9)*SQRT(F(I,NEK)* < 2F(I,NEK+1)) --- > 1(RA1(I,1)*F(I,NNOPK)*0.85+RA3(I,1)*F(I,NN2PK)*0.9)*SQRT(F(I,NEK) > 2*F(I,NEK+1)) ======================================================================== Diff of /home/tgcm/tgcm14a/cmpn4s.f and cmpn4s.f: 7c7 < include "fcom.h" --- > include "blnk.h" 13c13,14 < include "mwt.h" --- > real :: rmn4s,rmn2d,rmno,brn2d,cee > COMMON/MASS/ RMN4S,RMN2D,RMNO,BRN2D,CEE 19a21 > DATA RMN4S,RMN2D,RMNO/14.,14.,30./ ======================================================================== Diff of /home/tgcm/tgcm14a/cmpo2o.f and cmpo2o.f: 7c7 < include "fcom.h" --- > include "blnk.h" 13c13,14 < include "mwt.h" --- > real :: rmn4s,rmn2d,rmno,brn2d,cee > COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,CEE 71c72 < 1 *F(I,NPN2DK)/RMN2D)+RK1(I,1,1)*F(I,NOPK)+RK7(I,1)* --- > 1 *F(I,NPN2DK)/RMN2D)+RK1(I,1)*F(I,NOPK)+RK7(I,1)* 77,79d77 < ! < ! 1/7/00: fixed bug in following line (added "I" to F(NN2PK)) < ! S4(I,1) = RK3(I,1)*F(NN2PK)+RK8(I,1)*F(I,NNPK) 93c91 < 2 /RMN2D)+RK1(I,1,1)*F(I,NOPK)+(RK6(I,1)+RK7(I,1))* --- > 2 /RMN2D)+RK1(I,1)*F(I,NOPK)+(RK6(I,1)+RK7(I,1))* ======================================================================== Diff of /home/tgcm/tgcm14a/comp.f and comp.f: 2d1 < use input_module,only: difhor 6c5,6 < include "fcom.h" --- > include "strt.h" > include "blnk.h" 21,30c21,24 < PHI(:,1)=(/0.,0.673/) < PHI(:,2)=(/1.35,0./) < PHI(:,3)=(/1.11,0.769/) < TAU=1.86e+3 < DELTA(:,1)=(/1.,0./) < DELTA(:,2)=(/0.,1./) < T00=273. < ! SMALL = 1.e-20 < SMALL = 1.e-6 < ! --- > DATA PHI/0.,0.673,1.35,0.,1.11,0.769/,TAU/1.86E+3/, > ADELTA/1.,0.,0.,1./, > BT00/273./ > DATA SMALL/1.E-20/ 34c28 < IF (difhor.EQ.1) CALL DFACT(T1,J,IMAXP4) --- > IF(ICOMP.EQ.1)CALL DFACT(T1,J,IMAXP4) ======================================================================== Diff of /home/tgcm/tgcm14a/con.f and con.f: 2,3d1 < use input_module,only: step < use init_module,only: p0,iday 7c5 < include "fcom.h" --- > include "blnk.h" 13,16c11,13 < ! include "strt.h" ! iday < include "amie.h" < real :: xmue < COMMON/EDYVISC/XMUE(ZKMXP) --- > include "strt.h" > include "unit.h" > include "sechis.h" 23,24c20 < ! Moved to types.f: < ! DATA ISAVE/3/ --- > DATA ISAVE/3/ 36c32,37 < iuamie = 16 ! amie.h --- > IHIST = 12 > ISORC = 11 > ! IUNIN=8 > ! IUNOT=9 > IMAG=17 > IUAMI = 16 37a39,40 > MODEFC=-1 > IF(MODEST.EQ.3) MODEFC=0 39c42,43 < C --- > IFRST=0 > NERR=0 48,49c52 < ! 1/00: removing nhemi < ! IF(NHEMI.NE.0) DPHI=PI/(2*JMAX) --- > IF(NHEMI.NE.0) DPHI=PI/(2*JMAX) 89a93,104 > c > c lbsech = j-slice buffer length for secondary histories > c (dependent on fields selected in secflds input) > c > lbsech = 1 ! secondary hist buffer length per lat > if (lusech.gt.0) then > do n=1,mxfsech > if (len_trim(secflds(n)).gt.0) lbsech = lbsech + imaxp2*kmaxp1 > enddo > write(6,"('Buffer length for secondary ', > + 'histories (lbsech) = ',i6,' (per latitude)')") lbsech > endif 155,161c170,173 < ! WRITE(6,1010)MAXFLD,MAXPHY,NCBUF,NCPHYS < !1010 FORMAT(//* F ARRAY DIMENSIONED FOR FLDX=*I5,2X,*PHYX=*I5, < ! A*.*/* MODEL NEEDS NCBUF=*I5,2X,*NCPHYS=*I5, < ! B*.*/) < write(6,"(/,'F-array dimensioned for fldx=',i5,' phyx=',i5, < | /,' Model needs ncbuf=',i5,' ncphy=',i5,/)") < | maxfld,maxphy,ncbuf,ncphys --- > WRITE(6,1010)MAXFLD,MAXPHY,NCBUF,NCPHYS > 1010 FORMAT(//* F ARRAY DIMENSIONED FOR FLDX=*I5,2X,*PHYX=*I5, > A*.*/* MODEL NEEDS NCBUF=*I5,2X,*NCPHYS=*I5, > B*.*/) 164,168c176,178 < ! WRITE(6,1020)NCBUF,NCPHYS < !1020 FORMAT(//* EXIT CALLED BY CON*/* NEED FLDX.GE.*I5,2X,*PHYX.GE. < ! A*I5) < write(6,"(/,'Exit called by con: need fldx.ge.',i5,' phyx.ge.', < | i5)") ncbuf,ncphys --- > WRITE(6,1020)NCBUF,NCPHYS > 1020 FORMAT(//* EXIT CALLED BY CON*/* NEED FLDX.GE.*I5,2X,*PHYX.GE. > A*I5) 172c182 < CALL SETfft(TRIGS,IFAX,IMAX) --- > CALL SET99(TRIGS,IFAX,IMAX) 178c188 < ! IF(NHEMI.NE.0) JS=0 --- > IF(NHEMI.NE.0) JS=0 187a198,200 > C **** IF HEMI CASE, COS HAS SAME SIGN ACROSS EQUA > IF(NHEMI.NE.0) CSSP(1)=CS(2) > IF(NHEMI.NE.0) CSSP(2)=CS(1) 196c209 < DELTA=ATAN(TAN(23.5*PI/180.)*SIN(2.*PI*FLOAT(IDAY-80)/365.)) --- > DELTA=ATAN(TAN(23.5*PI/180.)*SIN(2.*PI*FLOAT(IIDAY-80)/365.)) 201d213 < c(4) = float(step) 254,256d265 < ! < ! 2/00: f107 and f107a are in input_mod.f. < ! c(61),c(62) are no longer used. 264d272 < p0 = c(81) ! init_mod 289a298,304 > ! > ! Intitialize SFEPS: Default solar flux variation due to orbital > ! eccentricity assumes no calendar day advance, thus no variation. > ! If calendar day *is* being advanced, SFEPS will be calculated > ! for each new day in advnce.f. SFEPS is in common in cons.h. > ! > SFEPS = 1. 292,293c307,308 < 70 FORMAT('IMAX=',I3,2X,'JMAX=',I2,2X,'KMAX=',I2,2X,'NFLDS=',I3,2X, < | 'NPHYS=',I3,2X,'DT=',E10.3,2X,'DS=',E10.3) --- > 70 FORMAT(*1IMAX=*I3,2X,*JMAX=*I2,2X,*KMAX=*I2,2X,*NFLDS=*I3,2X, > A*NPHYS=*I3,2X,*DT=*E10.3,2X,*DS=*E10.3) 295c310 < 80 FORMAT(' LBDSK=',I6,2X,'LBSCM=',I6,2X,'LBPHY=',I6,2X,'LBTAP=',I6) --- > 80 FORMAT(* LBDSK=*I6,2X,*LBSCM=*I6,2X,*LBPHY=*I6,2X,*LBTAP=*I6) 297c312 < 95 FORMAT('KUT = ',20I5,/,(7X,20I5)) --- > 95 FORMAT(1X,*KUT = *20I5/(7X,20I5)) 299,300c314,315 < 90 FORMAT('0 J',4X,'CS(J)',7X,'F (J)',7X,'SN(J)',7X,'TN(J)',/,/, < | (3X,E12.3)) --- > 90 FORMAT(*0 J*,4X,*CS(J)*,7X,*F (J)*7X,*SN(J)*,7X,*TN(J)*// > A(3X,E12.3)) 302c317 < 100 FORMAT(' ',I2,4E12.3) --- > 100 FORMAT(1H ,I2,4E12.3) 306c321 < 130 FORMAT('C=',10E12.3,/,(3X,10E12.3)) --- > 130 FORMAT(*1C=*10E12.3/(3X,10E12.3)) 327,328c342,343 < 135 FORMAT('0 K',5X,'EXPS(K)',8X,'A(K)',5X,'DIFK(K)',5X,'DIFT(K)', < | /,/,(' ',I2,4E12.3)) --- > 135 FORMAT(*0 K*,5X,*EXPS(K)*,8X,*A(K)*,5X,*DIFK(K)*,5X,*DIFT(K)*// > A(1H ,I2,4E12.3)) ======================================================================== Diff of /home/tgcm/tgcm14a/dfact.f and dfact.f: 6,7c6,8 < real :: ttbound,nob,avto,hor < COMMON/RUNTIM/TTBOUND,NOB(ZJMX),AVTO,HOR(ZJMX) --- > integer :: nerd > real :: hor > COMMON/RUNTIM/NERD(ZJMX+2),HOR(ZJMX) ======================================================================== Diff of /home/tgcm/tgcm14a/dt.f and dt.f: 2a3 > implicit none 8,11d8 < use input_module,only: iuivi < use init_module,only: iter < use bndry_module,only: tb,tb2,tba,bnd,bnd2,bnda,ci < implicit none 13c10 < include "fcom.h" --- > include "blnk.h" 19c16 < ! include "strt.h" ! for iter --- > include "strt.h" 202c199 < ! write(6,"('dt: t1 (lower boundary, to be stored at kmaxp1):', --- > ! write(6,"('dt: t1 (lower boundary, to be stored at kmaxp1:', 240c237 < S2(I,1) = F(I,NTNMK)-C(26)*(F(I,NTJP2K)+F(I,NTJM2K)-4.* --- > S2(I) = F(I,NTNMK)-C(26)*(F(I,NTJP2K)+F(I,NTJM2K)-4.* 248,249c245,246 < S3(I,1) = S2(I,1)-C(26)*(S2(I+2,1)+S2(I-2,1)-4.*(S2(I+1,1)+ < | S2(I-1,1))+6.*S2(I,1)) --- > S3(I) = S2(I)-C(26)*(S2(I+2)+S2(I-2)-4.*(S2(I+1)+S2(I-1))+ > 1 6.*S2(I)) 375,378c372,375 < S5(I,1) = S5(I,1)+S4(I,1)*1.5 < C S5(I,1) = S5(I,1)+S4(I,1)*2. < C S5(I,1) = S5(I,1)+S4(I,1)*3. < C S5(I,1) = S5(I,1)+S4(I,1)*1. --- > S5(I,1) = S5(I,1)+S4(I,1)*1.5 > C S5(I,1) = S5(I,1)+S4(I,1)*2. > C S5(I,1) = S5(I,1)+S4(I,1)*3. > C S5(I,1) = S5(I,1)+S4(I,1)*1. 382,384c379,380 < ! subroutine addfsech(name,long_name,units,f2d,londim,levdim, < ! | nlev,lat) < ! --- > c subroutine addfsech(fname,f,idim1,idim2,ilat) > c 388,390c384 < ! call addfsech('QSOLAR','SOLAR HEATING from DT','DEG K', < ! | s10,zimxp,zkmxp,zkmx,j) < ! --- > ! call addfsech('QSOLAR',s10,zimxp,kmaxp1,kmaxp1,j) 394,396c388,389 < ! call addfsech('QJOULE','JOULE HEATING from DT','DEG K', < ! | s4,zimxp,zkmxp,zkmx,j) < ! --- > ! call addfsech('QJOULE',s4,zimxp,kmaxp1,kmaxp1,j) > c 670c663 < F(I,NTK) = merge(F(I,NTK),100.,F(I,NTK)-100.>=0.) --- > F(I,NTK) = merge(F(I,NTK),100.,F(I,NTK)-100.>=0.) ======================================================================== Diff of /home/tgcm/tgcm14a/dynamics.f and dynamics.f: 3d2 < use bndry_module,only: bndcmp 8,12c7,10 < ! (fcom.h is task common). < ! These routines and those below them use f(:,njm2-njm1-nj-njp1-nmp2) < ! from the previous time iteration and output new latitude calculations < ! in f(:,njnp). nj is integer pointer to current latitude fields in < ! the f-array. (input arg lat is current latitude index). --- > ! (blnk.h is task common). > ! These routines use f(:,njm2-njm1-nj-njp1-nmp2) from previous > ! iteration and output f(:,njnp). nj is pointer to current > ! latitude fields (lat is current latitude index). 21c19 < include "fcom.h" ! debug only --- > include "blnk.h" ! debug only 32d29 < ! call mkfsech(lat,ixtime) 40a38,41 > ! NQPK=NQO2P > ! DO 1 I=1,5*LEN3 > ! F(I,NQPK)=0. > ! 1 CONTINUE 46c47 < ! call altv ! calc tvib: see tvib in crates.h, rates.f --- > C call altv 104d104 < implicit none 109c109 < include "fcom.h" --- > include "blnk.h" 112d111 < integer :: nqpk,i ======================================================================== Diff of /home/tgcm/tgcm14a/elden.f and elden.f: 7c7 < include "fcom.h" --- > include "blnk.h" 14c14,15 < include "mwt.h" --- > real :: rmn4s,rmn2d,rmno,brn2d,cee > COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,CEE 74c75 < 1S12(I,1)*(RK1(I,1,1)*F(I,NOPK)+RK6(I,1)* --- > 1S12(I,1)*(RK1(I,1)*F(I,NOPK)+RK6(I,1)* 88,100d88 < < ! < ! f(i,nqn2pk) (from orora) bad on irix here (all others ok): < ! < ! write(6,"('i=',i5,' f(i,nqn2pk)=',e12.4,' f(i,nqn2pk+1)=', < ! | e12.4)") i,f(i,nqn2pk),f(i,nqn2pk+1) < ! write(6,"(' rk16=',e12.4,' xiop2p=',e12.4)") rk16(i,1), < ! | xiop2p(i,1) < ! write(6,"(' rk23=',e12.4,' rmass(3)=',e12.4)") rk23(i,1), < ! | rmass(3) < ! write(6,"(' f(i,nps1k)=',e12.4,' f(i,nps2k)=',e12.4)") < ! | f(i,nps1k),f(i,nps2k) < ! if (i > 5) stop 'elden' 145,147d132 < ! < ! 2/00: changed args to vquart (s.a. vquart.f): < ! old call: CALL VQUART(S15,S5,S3,S2,S1,LEN3,1,LEN2,j) 149,153d133 < ! < ! s15,s14,s13,s12,s11 are coefficients for the quartic solver. < ! Positive roots of quartic are returned in s5: < ! < ! call vquart(s15,s14,s13,s12,s11, s5,imaxp4,kmax,kmaxp1,j) 165,175c145,154 < C S5(I,1) = merge(S5(I,1),3.3E+3,S5(I,1)-3.3E+3>=0.) < S5(I,1) = merge(S5(I,1),3.1E+3,S5(I,1)-3.1E+3>=0.) < C S5(I,1) = merge(S5(I,1),2.8E+3,S5(I,1)-2.8E+3>=0.) < C S5(I,1) = merge(S5(I,1),2.5E+3,S5(I,1)-2.5E+3>=0.) < F(I,NN2PK)=S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1)) < F(I,NO2PK)=(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1)))/ < | (S8(I,1)+RA2(I,1)*S5(I,1)) < F(I,NNOPK)=(S10(I,1)+S7(I,1)*(S6(I,1)-S4(I,1))/(S6(I,1)+ < | RA3(I,1)*S5(I,1))+S8(I,1)*(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+ < | RA3(I,1)*S5(I,1)))/(S8(I,1)+RA2(I,1)*S5(I,1)))/(RA1(I,1)* < | S5(I,1)) --- > C S5(I,1) = merge(S5(I,1),3.3E+3,S5(I,1)-3.3E+3>=0.) > S5(I,1) = merge(S5(I,1),3.1E+3,S5(I,1)-3.1E+3>=0.) > C S5(I,1) = merge(S5(I,1),2.8E+3,S5(I,1)-2.8E+3>=0.) > C S5(I,1) = merge(S5(I,1),2.5E+3,S5(I,1)-2.5E+3>=0.) > F(I,NN2PK)=S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1)) > F(I,NO2PK)=(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1)))/ > 1(S8(I,1)+RA2(I,1)*S5(I,1)) > F(I,NNOPK)=(S10(I,1)+S7(I,1)*(S6(I,1)-S4(I,1))/(S6(I,1)+RA3(I,1)* > 1S5(I,1))+S8(I,1)*(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+RA3(I,1)* > 2S5(I,1)))/(S8(I,1)+RA2(I,1)*S5(I,1)))/(RA1(I,1)*S5(I,1)) 183c162 < F(I,NEK+1)=SQRT(S5(I,1)*S5(I,2)) --- > F(I,NEK+1)=SQRT(S5(I,1)*S5(I,2)) 186,189c165,168 < F(I,NEK)=SQRT(S5(I,1)**3/S5(I,2)) < F(I,NEK+KMAX)=SQRT(S5(I,KMAX)**3/S5(I,KMAXM1)) < C F(I,NEK)= 1.5*S5(I,1)-0.5*S5(I,2) < C F(I,NEK+KMAX)=1.5*S5(I,KMAX)-0.5*S5(I,KMAXM1) --- > F(I,NEK)=SQRT(S5(I,1)**3/S5(I,2)) > F(I,NEK+KMAX)=SQRT(S5(I,KMAX)**3/S5(I,KMAXM1)) > C F(I,NEK)= 1.5*S5(I,1)-0.5*S5(I,2) > C F(I,NEK+KMAX)=1.5*S5(I,KMAX)-0.5*S5(I,KMAXM1) ======================================================================== Diff of /home/tgcm/tgcm14a/euvac.f and euvac.f: /home/tgcm/tgcm14a/euvac.f and euvac.f are identical ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/getgpi.f ======================================================================== Diff of /home/tgcm/tgcm14a/heelis.f and heelis.f: 2d1 < use input_module,only: aurora 8c7,8 < include "fcom.h" --- > include "dynphi.h" > include "blnk.h" 11a12 > include "strt.h" 14c15,18 < include "amie.h" --- > real :: bxm,bx,bxp,by,byp,bz,bzp,bmod,bmodp > COMMON/MAGFLD/BXM(ZIMXP,2),BX(ZIMXP,ZJMX),BXP(ZIMXP,4), > 1 BY(ZIMXP,ZJMX),BYP(ZIMXP,4),BZ(ZIMXP,ZJMX),BZP(ZIMXP,4), > 2 BMOD(ZIMXP,ZJMX),BMODP(ZIMXP,2) 19a24 > include "amie.h" 23c28 < integer :: nuik,nvik,nwik,i,jj,ishunk,ihem --- > integer :: nuik,nvik,nwik,i,jj,ishunk 38,42c43 < ! < ! Aurora routines flowv3, polht, and aurht are called only when < ! abs(latitude) <= 32.5 (also orora, called after heelis) < ! (ISHUNK and USERLA are local) < ! --- > ! userla and ishunk are local (s.a. orora): 45,74c46,54 < IF(aurora.EQ.0) ISHUNK=0 < IF(ABS(USERLA).LT.C(110)/6.) ISHUNK=0 < if (ishunk==1) then < ! < ! Amie not yet tested w/ multi-tasking. < IF (IAMIE .EQ. 1) THEN < write(6,"('heelis calling amiepa')") < CALL AMIEPA (T2(3),T1(3),T4(3),T5(3),T6(3),T8(3),T3(3), < | WION(3)) < GO TO 499 < ENDIF < ! < ! Call flowv3, polht, and aurht (flowxx, flowx0, flowv1, and flowv2 < ! have been eliminated). < ! Flowv3 defines taskcommon /cflowv3/, which is needed by < ! polht and aurht. < ! Flowv3 also defines ihem (local), which is passed to polht and aurht. < ! < call flowv3(IHEM) < ! < ! T5=CUSP, T6=ALFA, T8=FLUX, T3=DRIZL T7=ALFA2 T9=FLUX2 < ! < ! polht defines cusp in t5: < ! cusp < call polht(t5(3),imax,IHEM) < ! < ! aurht defines alfa,flux,drizl,alfa2,flux2: < ! alfa flux drizl alfa2 flux2 < call aurht(t6(3),t8(3),t3(3),t7(3),t9(3),imax,IHEM) < ! --- > IF(IAUR.EQ.0)ISHUNK=0 > IF(ABS(USERLA).LT.C(110)/6.)ISHUNK=0 > IF(ISHUNK.EQ.0)GO TO 24 > IF (IAMIE .EQ. 1) THEN > CALL AMIEPA (T2(3),T1(3),T4(3),T5(3),T6(3),T8(3),T3(3),WION(3)) > GO TO 499 > ENDIF > JJ=((2*J-JMAX-1)/IABS(2*J-JMAX-1)+3)/2 > CALL FLOWXX(T5(3),T6(3),T8(3),T3(3)) 84d63 < ! 92,105c71,91 < ! < ! ishunk==0 (abs(lat) < 32.5) < else < DO 26 I=1,LEN1 < T3(I)=0. < T5(I)=0. < T6(I)=0. < T7(I)=0. < T8(I)=0. < T9(I)=0. < 26 CONTINUE < ! < ! End latitude >= 32.5 < endif ! ishunk --- > write(6,"(/,'heelis after aurht: j=',i2)") j > write(6,"('t3=',/,(6e12.4))") t3 > write(6,"('t5=',/,(6e12.4))") t5 > write(6,"('t6=',/,(6e12.4))") t6 > write(6,"('t7=',/,(6e12.4))") t7 > write(6,"('t8=',/,(6e12.4))") t8 > write(6,"('t9=',/,(6e12.4))") t9 > GO TO 25 > 24 CONTINUE > DO 26 I=1,LEN1 > T3(I)=0. > T5(I)=0. > T6(I)=0. > T7(I)=0. > T8(I)=0. > T9(I)=0. > 26 CONTINUE > 25 CONTINUE > C **** > C **** SAVE AURORAL PARAMETERS > C **** ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/input.f ======================================================================== Diff of /home/tgcm/tgcm14a/lamdas.f and lamdas.f: 12c12 < include "fcom.h" --- > include "blnk.h" 29d28 < 35,36d33 < < !$OMP THREADPRIVATE (/vscr/) 57d53 < ! T4(I)=SIN(DIPMAG(I,J)) 58a55 > C T4(I)=SIN(DIPMAG(I,J)) ======================================================================== Diff of /home/tgcm/tgcm14a/magdyn.f and magdyn.f: 9c9 < include "fieldz.h" --- > include "field.h" 33c33 < DIPMAG(I+2,J) = ATAN(ZZB(I,J)/SQRT(XB(I,J)**2+YB(I,J)**2)) --- > DIPMAG(I+2,J) = ATAN(ZB(I,J)/SQRT(XB(I,J)**2+YB(I,J)**2)) 41c41 < BZ(I,J) = -ZZB(I,J)/BMOD(I,J) --- > BZ(I,J) = -ZB(I,J)/BMOD(I,J) ======================================================================== Diff of /home/tgcm/tgcm14a/new.f and new.f: 1a2 > implicit none 8,9d8 < use input_module,only: colfac < implicit none 11c10 < include "fcom.h" --- > include "blnk.h" 15a15,16 > real :: rmn4s,rmn2d,rmno,brn2d,colfac > COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,COLFAC ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/nextlu.f ======================================================================== Diff of /home/tgcm/tgcm14a/opflux.f and opflux.f: 2d1 < use init_module,only: secs 10a10 > include "strt.h" 17d16 < real t7opflux(zimxp,zkmxp) 49,51c48 < ! 1/00: remove nhemi: < ! RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*C(2) < RLAT=-.5*C(110)+(FLOAT(J-1)+.5)*C(2) --- > RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*C(2) ======================================================================== Diff of /home/tgcm/tgcm14a/oplus.f and oplus.f: 9c9 < include "fcom.h" --- > include "blnk.h" 16d15 < include "mwt.h" 22a22,23 > real :: rmn4s,rmn2d,rmno,brn2d,cee > COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,CEE 32d32 < character(len=80) :: title ======================================================================== Diff of /home/tgcm/tgcm14a/orora.f and orora.f: 2d1 < use input_module,only: aurora 13c12 < include "fcom.h" --- > include "blnk.h" 19,20d17 < include "ovalr.h" < include "aurlp.h" 21a19,28 > real :: rrad,h0,rh,rroth,e0,ree,rrote,fc,alfac,fd,alfad, > | alfk,alf6,alf21,rrot6,rrot21,rd6,rd6v,rd21,rh6,rh21,rt6, > | rt21,alfa0,ralfa,alfa20,ralfa2,e20,re2 > COMMON /OVALR/ RRAD(2),H0,RH,RROTH,E0,REE,RROTE,FC,ALFAC,FD,ALFAD > 1, ALFK,ALF6,ALF21,RROT6,RROT21,RD6,RD6V,RD21,RH6,RH21,RT6,RT21 > 2, ALFA0,RALFA,ALFA20,RALFA2,E20,RE2 > c > c Low energy auroral proton input (see ALFALP and EFLUXLP from input): > real :: alfa_lp,eflux_lp,qteaur > common/aurlp/ alfa_lp(zimxp),eflux_lp(zimxp),qteaur(zimxp) 32d38 < character(len=80) :: title 34c40 < ! Ion production from low energy proton source --- > ! Ion production from low energy proton source 56c62 < IF (aurora.EQ.0) ISHUNK=0 --- > IF (IAUR.EQ.0) ISHUNK=0 58c64 < IF (ISHUNK.EQ.0) then --- > IF(ISHUNK.EQ.0) then 61,79c67,72 < ! < ! Put qia on secondary histories (see also below, when ishunk > 0): < ! < ! call addfsech('QIA1_N2P', < ! | 'N2+ production from low-energy proton source', < ! | ' ',qia(:,:,1),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_O2P', < ! | 'O2+ production from low-energy proton source', < ! | ' ',qia(:,:,2),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_OP', < ! | 'O+ production from low-energy proton source', < ! | ' ',qia(:,:,3),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_NOP', < ! | 'NO+ production from low-energy proton source', < ! | ' ',qia(:,:,4),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_NP', < ! | 'N+ production from low-energy proton source', < ! | ' ',qia(:,:,5),zimxp,zkmxp,zkmx,j) < return --- > ! call addfsech('QIA1_N2+',qia(:,:,1),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA2_O2+',qia(:,:,2),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA3_O+ ',qia(:,:,3),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA4_NO+',qia(:,:,4),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA5_N+ ',qia(:,:,5),zimxp,zkmxp,zkmx,j) > RETURN 104,106d96 < ! 9/99: aion is a cpu consumer because of **real < ! and exp. It is located in inline.f. < ! 114c104 < ! (top level zkmxp not defined) (s.a., aurht.f) --- > ! (top level zkmxp not defined) 120,121d109 < ! subroutine addfsech(name,long_name,units,f2d,londim,levdim, < ! | nlev,lat) 123,137c111,115 < ! call addfsech('QIA1_N2P', < ! | 'N2+ production from low-energy proton source', < ! | ' ',qia(:,:,1),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_O2P', < ! | 'O2+ production from low-energy proton source', < ! | ' ',qia(:,:,2),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_OP', < ! | 'O+ production from low-energy proton source', < ! | ' ',qia(:,:,3),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_NOP', < ! | 'NO+ production from low-energy proton source', < ! | ' ',qia(:,:,4),zimxp,zkmxp,zkmx,j) < ! call addfsech('QIA1_NP', < ! | 'N+ production from low-energy proton source', < ! | ' ',qia(:,:,5),zimxp,zkmxp,zkmx,j) --- > ! call addfsech('QIA1_N2+',qia(:,:,1),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA2_O2+',qia(:,:,2),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA3_O+ ',qia(:,:,3),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA4_NO+',qia(:,:,4),zimxp,zkmxp,zkmx,j) > ! call addfsech('QIA5_N+ ',qia(:,:,5),zimxp,zkmxp,zkmx,j) 178c156 < c (called above) for N2+, O2+, O+, NO+, N+ respectively --- > c (called above) for N2+, O2+, O+, NO+, N+ respectively 180d157 < ! 2/28/00: added to tgcm14 188,189c165,167 < S15(I,1)=S1(I,1)*0.7*merge(1.-F(I,NPS1K)-F(I,NPS2K),0.,1.- < 1 F(I,NPS1K)-F(I,NPS2K)>=0.)/(RMASS(3)*S3(I,1))+qia(i,1,1) --- > S15(I,1)=S1(I,1)*0.7*merge(1.-F(I,NPS1K)-F(I,NPS2K),1.E-8, > 1 (1.-F(I,NPS1K)-F(I,NPS2K))-1.E-8>0.)/(RMASS(3)*S3(I,1)) > 2 +qia(i,1,1) 193,194c171,172 < K=(3-M)*KMAXP1+1 < NQPK=NQO2P+(M-1)*KMAXP1 --- > K=(3-M)*KMAXP1+1 > NQPK=NQO2P+(M-1)*KMAXP1 196,198c174,176 < DO 36 I=LEN1+1,LEN2 < F(I,NQPK)=F(I,NQPK)+SQRT(S15(I,K)*S15(I,K-1)) < 36 CONTINUE --- > DO 36 I=LEN1+1,LEN2 > F(I,NQPK)=F(I,NQPK)+SQRT(S15(I,K)*S15(I,K-1)) > 36 CONTINUE 200,204c178,182 < DO 37 I=1,LEN1 < F(I,NQPK) = F(I,NQPK)+1.5*S15(I,K)-0.5*S15(I,K+1) < F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+1.5*S15(I,K+KMAXM1)-0.5* < 1 S15(I,K+KMAX-2) < 37 CONTINUE --- > DO 37 I=1,LEN1 > F(I,NQPK) = F(I,NQPK)+1.5*S15(I,K)-0.5*S15(I,K+1) > F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+1.5*S15(I,K+KMAXM1)-0.5* > 1 S15(I,K+KMAX-2) > 37 CONTINUE 210c188 < F(I,NQPK)=F(I,NQPK)+.22/.7*SQRT(S15(I,1)*S15(I,0)) --- > F(I,NQPK)=F(I,NQPK)+.22/.7*SQRT(S15(I,1)*S15(I,0)) 214,216c192,194 < F(I,NQPK) = F(I,NQPK)+.22/.7*(1.5*S15(I,1)-0.5*S15(I,2)) < F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+.22/.7*(1.5*S15(I,KMAX) < 1 -0.5*S15(I,KMAXM1)) --- > F(I,NQPK) = F(I,NQPK)+.22/.7*(1.5*S15(I,1)-0.5*S15(I,2)) > F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+.22/.7*(1.5*S15(I,KMAX) > 1 -0.5*S15(I,KMAXM1)) 218,219d195 < ! < ! 2/28/00: disn2p (see crates.h, qrj.f, qtieff.f, xray.f) 222c198 < DISN2P(I,K)=DISN2P(I,K)+SQRT(S15(I,K-1)*S15(I,K)) --- > DISN2P(I,K)=DISN2P(I,K)+SQRT(S15(I,K-1)*S15(I,K)) 225,226c201,202 < DISN2P(I,1) = DISN2P(I,1)+1.5*S15(I,1)-0.5*S15(I,2) < DISN2P(I,KMAXP1) = DISN2P(I,KMAXP1)+1.5*S15(I,KMAX)-0.5* --- > DISN2P(I,1) = DISN2P(I,1)+1.5*S15(I,1)-0.5*S15(I,2) > DISN2P(I,KMAXP1) = DISN2P(I,KMAXP1)+1.5*S15(I,KMAX)-0.5* ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/othend.f ======================================================================== Diff of /home/tgcm/tgcm14a/proton.f and proton.f: 17c17 < include 'fcom.h' --- > include 'blnk.h' ======================================================================== Diff of /home/tgcm/tgcm14a/qinite.f and qinite.f: 7c7 < include "fcom.h" --- > include "blnk.h" 14c14,15 < include "mwt.h" --- > real :: rmn4s,rmn2d,rmno,brn2d,cee > COMMON/MASS/ RMN4S,RMN2D,RMNO,BRN2D,CEE 132,134d132 < ! F(I,NQOPK)=F(I,NQOPK)+(0.33*S6(I,1)+S7(I,1))*S5(I,1) < ! F(I,NQN2PK)=F(I,NQN2PK)+0.86*S8(I,1)*S5(I,1) < ! F(I,NQNPK)=F(I,NQNPK)+0.14*S8(I,1)*S5(I,1) ======================================================================== Diff of /home/tgcm/tgcm14a/qrj.f and qrj.f: 1d0 < ! 3,5c2 < use input_module,only: f107 < use init_module,only: sfeps < ! CDIR$ VFUNCTION EXPHF --- > CDIR$ VFUNCTION EXPHF 11c8 < include "fcom.h" --- > include "blnk.h" 19c16,17 < include "mwt.h" --- > real :: rmn4s,rmn2d,rmno,brn2d,cee > COMMON/MASS/ RMN4S,RMN2D,RMNO,BRN2D,CEE 31c29,30 < real :: do2,ad,bd,aband,do22,e3,factor,bband,cband --- > real :: do2,ad,bd,aband,do22,e3,f107,f107a,factor,bband,cband > real :: exphf 73,75c72,74 < S3(I,1) = exp(-S3(I,1)) < S4(I,1) = exp(-S4(I,1)) < S5(I,1) = exp(-S5(I,1)) --- > S3(I,1) = exphf(-S3(I,1)) > S4(I,1) = exphf(-S4(I,1)) > S5(I,1) = exphf(-S5(I,1)) 79,80c78,79 < S6(I,1) = exp(-S1(I,1)) < S7(I,1) = exp(-S2(I,1)) --- > S6(I,1) = exphf(-S1(I,1)) > S7(I,1) = exphf(-S2(I,1)) 179,181d177 < ! This prevents underflow on babyblue, but gets wrong results: < ! S5(I,1)=merge((S5(I,1)+SIGEUV(M,N)*F(I,NNK)),0., < ! | (S5(I,1)+SIGEUV(M,N)*F(I,NNK)) < 1.e-200) 193c189 < S5(I,1) = FEUV(N)*exp(-S5(I,1)) --- > S5(I,1) = FEUV(N)*exphf(-S5(I,1)) 197,201d192 < ! S5 goes wacko here on babyblue: < ! write(6,"('qrj: j=',i2,' i=',i4,' n=',i2,' rlmeuv(n)=', < ! | e12.4,' s5=',e12.4,' s6=',e12.4)") j,i,n,rlmeuv(n), < ! | s5(i,1),s6(i,1) < 310c301 < DISN2P(I,1) = S10(I,1)*S14(I,1)*S5(I,1) --- > DISN2P(I,1) = S10(I,1)*S14(I,1)*S5(I,1) 366c357 < S9(I,1)=SIGSRC(N)*FSRC(N)*exp(-SIGSRC(N)*F(I,NNO2K)) --- > S9(I,1)=SIGSRC(N)*FSRC(N)*exphf(-SIGSRC(N)*F(I,NNO2K)) 386c377 < + F(I,NNO2K)-1.E18>=0.0)*(1.+0.11*(f107-65.)/165.)*sfeps --- > + F(I,NNO2K)-1.E18>=0.0)*(1.+0.11*(C(61)-65.)/165.)*SFEPS 391c382 < end subroutine qrj --- > END 394,395c385 < use input_module,only: f107,f107a < use init_module,only: sfeps --- > implicit none 402a393,394 > ! > ! Local: 405a398 > real :: EUVFLX(37),wave1(lmax),wave2(lmax) ! output from euvac 408,409c401 < real :: EUVFLX(37),wave1(lmax),wave2(lmax) ! output from euvac < real :: flya,hlybr,fexvir,hlya,heiew,xuvfac --- > real :: flya,f107,f107a,hlybr,fexvir,hlya,heiew,xuvfac 411,418c403,406 < ! < ! 2/00: < ! Use f107,f107a from input_mod (these are either provided by < ! the user, or obtained from GPI database, see input_mod.f and < ! gpi_mod.f). < ! F107 = C(61) < ! F107A = C(62) < ! --- > > F107 = C(61) > F107A = C(62) > ISCALE = 0 451c439 < FEUV(N) = SFLUX(N)*SFEPS --- > FEUV(N) = SFLUX(N)*SFEPS 454c442 < FSRC(N) = SFLUX(N)*SFEPS --- > FSRC(N) = SFLUX(N)*SFEPS 458c446 < FEUV(NN) = EUVFLX(N)*SFEPS --- > FEUV(NN) = EUVFLX(N)*SFEPS 599a588 > implicit none 608c597 < parameter (lmax=59) --- > integer,parameter :: lmax=59 622a612,614 > ! Local: > integer :: m,n,nn > ! ======================================================================== Diff of /home/tgcm/tgcm14a/qtieff.f and qtieff.f: 7c7 < include "fcom.h" --- > include "blnk.h" 11a12 > include "crates.h" ======================================================================== Diff of /home/tgcm/tgcm14a/rates.f and rates.f: 3,5c3 < use input_module,only: f107 < use init_module,only: sfeps < ! CDIR$ VFUNCTION EXPHF --- > CDIR$ VFUNCTION EXPHF 13c11 < include "fcom.h" --- > include "blnk.h" 22a21 > real :: exphf 39c38 < ! RK1(I,1,1)=(((9.65E-16*S15(I,1)-5.17E-14)*S15(I,1)+1.073E-12)* --- > ! RK1(I,1)=(((9.65E-16*S15(I,1)-5.17E-14)*S15(I,1)+1.073E-12)* 44c43 < C RK1(I,1,1)=1.7E-11*(300./F(I,NTK))**0.77+8.54E-11*exp(-3467./ --- > C RK1(I,1)=1.7E-11*(300./F(I,NTK))**0.77+8.54E-11*exphf(-3467./ 48c47 < RK1(I,1,1)=1.6E-11*S15(I,1)**(-0.52)+5.5E-11*exp(-22.85/S15(I,1)) --- > RK1(I,1)=1.6E-11*S15(I,1)**(-0.52)+5.5E-11*exphf(-22.85/S15(I,1)) 56,58c55 < ! See also sub altv (altv.f). Tvib is output of sub altv, if it < ! is called (see dynamics.f). < ! --- > C **** 60c57 < S5(I,1) = exp(-3353./TVIB(I,1)) --- > S5(I,1) = exphf(-3353./TVIB(I,1)) 75c72 < RKM12(I,1) = 9.59E-34*exp(480./F(I,NTK)) --- > RKM12(I,1) = 9.59E-34*exphf(480./F(I,NTK)) 84,86c81,83 < BETA1(I,1) = 1.5E-11*exp(-3600./F(I,NTK)) < C BETA3(I,1) = 1.6E-10*exp(-460./F(I,NTK)) < BETA3(I,1) = 2.5E-10*SQRT(F(I,NTK)/300.)*exp(-600./ --- > BETA1(I,1) = 1.5E-11*exphf(-3600./F(I,NTK)) > C BETA3(I,1) = 1.6E-10*exphf(-460./F(I,NTK)) > BETA3(I,1) = 2.5E-10*SQRT(F(I,NTK)/300.)*exphf(-600./ 89,96c86,91 < BETA8(I,1)=4.5E-6*(1.+0.11*(f107-65.)/165.)* < | exp(-1.E-8*F(I,NNO2K)**0.38)*sfeps < < BETA9(I,1)=2.91E11*(1.+0.2*(f107-65.)/100.)*2.E-18* < | exp(-8.E-21*F(I,NNO2K))*sfeps < < BETA9N(I,1)=5.E9*(1.+0.2*(f107-65.)/100.)*2.E-18* < | exp(-8.E-21*F(I,NNVO2K))*sfeps --- > BETA8(I,1)=4.5E-6*(1.+0.11*(C(61)-65.)/165.)* > | exphf(-1.E-8*F(I,NNO2K)**0.38)*SFEPS > BETA9(I,1)=2.91E11*(1.+0.2*(C(61)-65.)/100.)*2.E-18* > | exphf(-8.E-21*F(I,NNO2K))*SFEPS > BETA9N(I,1)=5.E9*(1.+0.2*(C(61)-65.)/100.)*2.E-18* > | exphf(-8.E-21*F(I,NNVO2K))*SFEPS 100,105c95,100 < BETA8(I,KMAXP1)=4.5E-6*(1.+0.11*(f107-65.)/165.)* < | exp(-1.E-8*F(I,NNO2K+KMAX)**0.38) < BETA9(I,KMAXP1)=2.91E11*(1.+0.2*(f107-65.)/100.)*2.E-18* < | exp(-8.E-21*F(I,NNO2K+KMAX)) < BETA9N(I,KMAXP1)=5.E9*(1.+0.2*(f107-65.)/100.)*2.E-18* < | exp(-8.E-21*F(I,NNVO2K+KMAX)) --- > BETA8(I,KMAXP1)=4.5E-6*(1.+0.11*(C(61)-65.)/165.)* > | exphf(-1.E-8*F(I,NNO2K+KMAX)**0.38) > BETA9(I,KMAXP1)=2.91E11*(1.+0.2*(C(61)-65.)/100.)*2.E-18* > | exphf(-8.E-21*F(I,NNO2K+KMAX)) > BETA9N(I,KMAXP1)=5.E9*(1.+0.2*(C(61)-65.)/100.)*2.E-18* > | exphf(-8.E-21*F(I,NNVO2K+KMAX)) 143c138 < BETA4(I,1) = 5.0E-13 --- > BETA4(I,1) = 5.0E-13 ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/runmds.f ======================================================================== Diff of /home/tgcm/tgcm14a/settei.f and settei.f: 2c2,3 < use input_module,only: f107 --- > CDIR$ VFUNCTION EXPHF > implicit none 6,7d6 < use init_module,only: secs < implicit none 9c8 < include "fcom.h" --- > include "blnk.h" 12a12 > include "strt.h" 15d14 < include "aurlp.h" 21a21,22 > real :: alfa_lp,eflux_lp,qteaur > common/aurlp/ alfa_lp(zimxp),eflux_lp(zimxp),qteaur(zimxp) 27c28 < | coslat,sinlat,expsk --- > | coslat,sinlat,expsk,exphf 54c55 < F107TE = f107 --- > F107TE = C(61) 59,61c60 < ! 1/00: removing nhemi < ! RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*C(2) < RLAT=-.5*C(110)+(FLOAT(J-1)+.5)*C(2) --- > RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*C(2) 207c206 < T2(I) = merge(T2(I),.10,T2(I)-.10>=0.) --- > T2(I) = merge(T2(I),.10,T2(I)-.10>=0.) 234,237d232 < < ! write(6,"('settei: i=',i2,' s4=',e12.4,' s1=',e12.4, < ! | ' s10=',e12.4)") i,s4(i,1),s1(i,1),s10(i,1) < 324c319 < S10(I,1) = exp(-((((0.001996*S10(I,1)+0.08034)*S10(I,1)+1.166) --- > S10(I,1) = exphf(-((((0.001996*S10(I,1)+0.08034)*S10(I,1)+1.166) 342,343c337,338 < S9(I,1) = merge(2.E-7*exp(-4605.2/F(I,NTEK)),5.71E-8* < | exp(-3352.6/F(I,NTEK)),F(I,NTEK)-1000.>=0.) --- > S9(I,1) = merge(2.E-7*exphf(-4605.2/F(I,NTEK)),5.71E-8* > | exphf(-3352.6/F(I,NTEK)),F(I,NTEK)-1000.>=0.) 345c340 < | exp(-17620./F(I,NTEK)),2000.-F(I,NTEK)>=0.) --- > | exphf(-17620./F(I,NTEK)),2000.-F(I,NTEK)>=0.) 353c348 < | (1.-exp(S10(I,1)))/S10(I,1) --- > | (1.-exphf(S10(I,1)))/S10(I,1) 376,377c371,372 < | (1.-exp(S10(I,1)))/S10(I,1) < S10(I,1) = 1.57E-12*exp((2.4E4+0.3*(F(I,NTEK)-1500.)-1.947E-5 --- > | (1.-exphf(S10(I,1)))/S10(I,1) > S10(I,1) = 1.57E-12*exphf((2.4E4+0.3*(F(I,NTEK)-1500.)-1.947E-5 ======================================================================== Diff of /home/tgcm/tgcm14a/ssflux.f and ssflux.f: /home/tgcm/tgcm14a/ssflux.f and ssflux.f are identical ======================================================================== WARNING: cannot find file /home/tgcm/tgcm14a/start.f ======================================================================== Diff of /home/tgcm/tgcm14a/tail.f and tail.f: 1,11c1 < ! < subroutine tail(iprint) < use input_module,only: power,ctpoten,byimf,f107,f107a < ! < ! Set auroral parameters. Most of these are optionally read by < ! input (time-dependent input of these is not yet available) < ! < ! This version of tail assumes naurp=3 (only hp,cp,by are input), < ! and constant hp,cp,by in time. It also assumes JSWOLD=1 (OLDALFA) < ! from the old tail. < ! --- > SUBROUTINE TAIL (SECS) 12a3,13 > C **** > C **** ROUTINE TO CALCULATE AURORAL CONSTANTS AS FUNCTION OF > C **** TIME FOR THE NORMAL CASE USING THE REVISED ALFA (10/86) > C **** AND THE NEW POLAR CAP CONVECTION (5/87) > C **** SECS=UNIVERSAL TIME IN SECONDS WHERE > C **** 0.12 (MDAY.IHR) = 0.0 (UTMDAY.UTHR) > C **** 1.0 (MDAY.IHR) = 0.12 (UTMDAY.UTHR) > C **** 1.12 (MDAY.IHR) = 1.0 (UTMDAY.UTHR) > C **** > ! > ! Commons: 14c15,21 < include "amie.h" --- > include "cons.h" > integer,parameter :: NAURPX=55,NLEX=600 > real :: aurps,paramv > integer :: npts,ipr,naurp,jswolda,jdith,jdidk > COMMON/AURDAT/ AURPS(NLEX,2,NAURPX), NPTS(NAURPX),IPR, NAURP, > | PARAMV(NAURPX), JSWOLDA, JDITH, JDIDK > EQUIVALENCE (PARAMV(1),HP), (PARAMV(2),CP), (PARAMV(3),BY) 16,17c23,25 < ! Parameters from /ovalr/ that are set here (or received from input), < ! and referenced by orora: fc,alfac,fd,alfad, e30,fd2,alfad2 --- > ! Geophys indices: > real :: f107,f107a,ctpoten,hpower,byimf > common/ingpi/ f107,f107a,ctpoten,hpower,byimf 19c27,36 < include "ovalr.h" --- > ! AURORAL PARAMETERS IN DEGREES ETC. > real :: theta0,offa,offc,dskofa,dskofc,phid,phin,phidp0, > | phidm0,phinp0,phinm0,psim,psie,pcen,arad,h1,h2,roth,e1, > | e2,rote,ec,ed,alfa1,alfa2,twak,twa6,twa21,rot6,rot21, > | d6,d21,h6,h21,t6,t21 > COMMON /PARAMD/ THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),DSKOFC(2), > | PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),PHINM0(2), > | PSIM(2),PSIE(2),PCEN(2),ARAD(2),H1,H2,ROTH,E1,E2,ROTE, > | EC,ED,ALFA1,ALFA2,TWAK,TWA6,TWA21,ROT6,ROT21,D6,D21,H6,H21, > | T6,T21 21,22c38,48 < ! The first 15 parameters in ioncr.h are optional inputs (see input): < include "ioncr.h" --- > ! AURORAL PARAMETERS IN RADIANS ETC FOR FLOWV3 > integer :: istar,ihem > real :: rdeg,psiv,r1,dduumm > COMMON /IONCR/ ISTAR,IHEM,RDEG(22),PSIV(6),R1(2), > 1 DDUUMM(IMAXMP*JMX0+IMX0*JMX0+2) > real :: rrad,h0,rh,rroth,e0,ree,rrote,fc,alfac,fd,alfad, > | alfk,alf6,alf21,rrot6,rrot21,rd6,rd6v,rd21,rh6,rh21,rt6, > | rt21,alfa0,ralfa,alfa20,ralfa2,e20,re2 > COMMON /OVALR/ RRAD(2),H0,RH,RROTH,E0,REE,RROTE,FC,ALFAC,FD, > | ALFAD,ALFK,ALF6,ALF21,RROT6,RROT21,RD6,RD6V,RD21,RH6,RH21, > | RT6,RT21, ALFA0,RALFA,ALFA20,RALFA2,E20,RE2 23a50,52 > ! Amie common: > include "amie.h" > ! 25c54 < integer,intent(in) :: iprint --- > real,intent(in) :: secs 28,34c57,62 < real :: ec,ed,pi,dtr,e1,e2,rote,h1,h2,roth,rhp,rcp,raur,twak, < | twa6,twa21,rot6,rot21,t6,t21,d6,d6v,d21,h6,h21,e21,e22, < | raursh,raurnh,thetsh,alfa1,alfa2,e,ut,p0,hpgwsh,hpgwnh, < | plevel,hp,cp,arad(2),alfa21,alfa22 < real,parameter :: convrt=3.1211e+8 < integer,parameter :: ish=1,inh=2 ! south,north hemispheres < integer :: i,iut,n --- > integer :: ifrst,i,ihp,ihn,inh,ish,iut,mday,n,ih > real :: hp,cp,by,convrt,dtr,e,rhp,rcp,raursh,raurnh,hpgwsh, > | hpgwnh,pi,ut,seca,f1,f2,clat,cmlt,plevel,bygd,hp35,hpm35, > | e21,e22,alfa21,alfa22,ut24,dispc,disp,raur,d6v,thetsh,p0 > CHARACTER(len=8) :: NAMEHA(38) > real,external :: terp 36,37c64,65 < ! write(6,"('enter tail: /ingpi/ ctpoten, power, byimf=', < ! | 3e12.4)") ctpoten,power,byimf --- > DATA CONVRT/3.1211E8/,DTR/1.7453293E-2/,E/1.E-10/,IFRST/0/ > COMMON /RADII/ RHP,RCP,RAURSH,RAURNH,HPGWSH,HPGWNH 39,57c67,188 < pi = 4.*atan(1.) < dtr = pi/180. < e = 1.e-10 < ! < ! Set cusp and drizzle parameters alfac, alfad, fc, fd (ovalr.h): < ! (power, from input, is in ingpi.h) < ! < ! ec = 0.1+0.9*power/100. < ! ec = 0.01+0.09*power/100. < ! alfac = 0.5 < ec = 0.5 < alfac = 1.0 < ! ed = 0.1+2.*power/100. < ! ed = 0.01+0.2*power/100. < ! alfad = 0.75 < ed = 0.5 < alfad = 2.0 < fc = convrt * ec / alfac < fd = convrt * ed / alfad --- > ! HEELIS CASE PLUS AURORA (ASSYMMETRICAL) > DATA NAMEHA/ ' H POWER', ' C POTEN', ' IMF BY', ' R1', > | ' PHIDM', ' PHIDP', ' PHINM', ' PHINP', ' ARADP,N', > | ' ARADN,S', ' OFFAP,N', ' OFFAN,S', ' DKAP,N', ' DKAN,S', > | ' PHIDP,N', ' PHIDN,S', ' PHINP,N', ' PHINN,S', ' PCEPP,N', > | ' PCEPN,S', ' PCENP,N', ' PCENN,S', 'DPC/DKPN', 'DPC/DKNS', > | ' DP/THPN', ' DP/THNS', ' OFFCP,N', ' OFFCN,S', > | ' E1 ', ' E2 ', ' H1 ', ' H2 ', ' ROTE ', > | ' ROTH ', ' TWA6 ', ' TWA21 ', ' ALFA1 ', ' ALFA2 '/ > DATA PI/3.14159265358979/ > C > IFRST = IFRST + 1 > C IF (IFRST .GT. 1 .AND. IPR .EQ. 0) RETURN > C > C **** > C **** UNIVERSAL TIME IN HOURS > C **** > UT = SECS/3600. > C > C **** INTERPOLATE NAURP PARAMETERS > if (npts(54).ne.0) then ! f107 was not defined by getgpi > I = 54 > F107 = TERP( AURPS(1,1,I), AURPS(1,2,I), UT, NPTS(I) ) > IF (NPTS(I) .GT. 1) THEN > IF (UT .LT. AURPS(1,1,I) .OR. UT .GT. AURPS(NPTS(I),1,I)) THEN > WRITE(6,"(1X,'TAIL: STOP---UT OUT OF BOUNDS, UT', > + ' UTMIN UTMAX =',3F7.2,' F107')") > + UT,AURPS(1,1,I),AURPS(NPTS(I),1,I) > STOP > ENDIF > ENDIF > endif > if (npts(55).ne.0) then ! f107a was not defined by getgpi > I = 55 > F107A = TERP( AURPS(1,1,I), AURPS(1,2,I), UT, NPTS(I) ) > IF (NPTS(I) .GT. 1) THEN > IF (UT .LT. AURPS(1,1,I) .OR. UT .GT. AURPS(NPTS(I),1,I)) THEN > WRITE(6,"(1X,'TAIL: STOP---UT OUT OF BOUNDS, UT', > + ' UTMIN UTMAX =',3F7.2,' F107A')") > + UT,AURPS(1,1,I),AURPS(NPTS(I),1,I) > STOP > ENDIF > ENDIF > endif > DO 100 I=1,NAURP > if (npts(i).eq.0) goto 100 ! parameter was defined by getgpi > IF (NPTS(I).EQ.1) GO TO 99 > IF (UT .LT. AURPS(1,1,I) .OR. UT .GT. AURPS(NPTS(I),1,I)) THEN > WRITE(6,"(1X,'TAIL: STOP---UT OUT OF BOUNDS, UT UTMIN ', > + 'UTMAX =',3F7.2,2X,A8)") UT,AURPS(1,1,I),AURPS(NPTS(I),1,I), > + NAMEHA(I) > STOP > ENDIF > 99 PARAMV(I) = TERP(AURPS(1,1,I),AURPS(1,2,I),UT,NPTS(I)) > 100 continue > C > C REPLACE HP AND CP WITH N. HEM. VALUES FROM THE AMIE VOLUME > IF (IAMIE .EQ. 1) THEN > SECA = AMAX1(SECUTA(1),SECS) > SECA = AMIN1(SECUTA(2),SECA) > F1 = SECA - SECUTA(1) > F2 = SECUTA(2) - SECA > HP = (F2*HPNHA(1) + F1*HPNHA(2)) / FLOAT(NSPSEC) > CP = (F2*CPNHA(1) + F1*CPNHA(2)) / FLOAT(NSPSEC) > CLAT = (F2*CUSPSLA(1) + F1*CUSPSLA(2)) / FLOAT(NSPSEC) > CMLT = (F2*CUSPSTA(1) + F1*CUSPSTA(2)) / FLOAT(NSPSEC) > CRAD(1) = (90.-CLAT)*PI/180. > PHIDA(1) = (CMLT - 12.) * PI / 12. > CLAT = (F2*CUSPNLA(1) + F1*CUSPNLA(2)) / FLOAT(NSPSEC) > CMLT = (F2*CUSPNTA(1) + F1*CUSPNTA(2)) / FLOAT(NSPSEC) > CRAD(2) = (90.-CLAT)*PI/180. > PHIDA(2) = (CMLT - 12.) * PI / 12. > C Compute variables for eflx2 > PLEVEL = 2.09 * ALOG(HP) > C H1 = AMIN1 (2.35, 0.83 + 0.33 * PLEVEL ) > C H2 = 2.87 + 0.15 * PLEVEL > H1 = 3.+0.1*HP > H2 = 3.0+0.1*HP > ROTH = (12.18 - 0.89*PLEVEL) * 15. > ROTE = (2.62 - 0.55*PLEVEL) * 15. > ENDIF > C > PLEVEL = 0. > IF (HP .GE. 0.01) PLEVEL = 2.09 * ALOG(HP) > C **** PUT F107,F107A INTO CONSTANT ARRAY > ! (only if they were not defined by getgpi) > if (npts(54).ne.0) c(61) = f107 > if (npts(55).ne.0) c(62) = f107a > C > C FIND LIMITS OF BY WHICH ARE VALID FOR MODEL OF (3/89) > BYGD = AMIN1(7.,BY) > BYGD = AMAX1(-11.,BYGD) > IHP = 1 > IHN = 2 > IF (BY .GT. 0.) THEN > IHP = 2 > IHN = 1 > ENDIF > C ASSUME THAT WE ARE ALWAYS DEALING WITH S AND N HEM IN PARAMV > INH = 2 > ISH = 1 > C > HP35 = AMIN1(HP,35.) > HPM35 = AMAX1(0.,HP-35.) > C > C SET CUSP AND DRIZZLE PARAMETERS > C **** > C **** CHANGE CUSP AND DRIZZLE PARAMETERS > C **** (IN TGCM3 HAD EC=0.40, ALFAC=0.09, ED=0.05, ALFAD=0.20) > C **** (WAS TOO MUCH ED ESPECIALLY.) > C EC = 0.1+0.9*HP/100. > C EC = 0.01+0.09*HP/100. > C ALFAC = 0.5 > EC = 0.5 > ALFAC = 1.0 > C ED = 0.1+2.*HP/100. > C ED = 0.01+0.2*HP/100. > C ALFAD = 0.75 > ED = 0.5 > ALFAD = 2.0 > FC = CONVRT * EC / ALFAC > FD = CONVRT * ED / ALFAD 59,60c190,191 < E21 = 1.E-80 < E22 = 1.E-80 --- > E21 = 1.E-80 > E22 = 1.E-80 65,123c196,342 < do i=ish,inh < if (iamie==0) then < theta0(i) = -3.80 + 8.48 * (ctpoten**0.1875) < else < theta0(i) = -1.92 + 8.10 * (ctpoten**0.1875) < endif < theta0(i) = theta0(i) < offc(i) = 1.1 < phidp0(i) = 85. < phidm0(i) = 50. < phinp0(i) = 57.5 < phinm0(i) = 100. < psim(i) = 0.44 * ctpoten < psie(i) = -0.56 * ctpoten < enddo < offa(ish) = 4.3 < offa(inh) = 3.7 < dskofa(ish) = (-1.26 + 0.15 * byimf) < dskofa(inh) = (-1.26 - 0.15 * byimf) < dskofc(ish) = (-0.08 + 0.15 * byimf) < dskofc(inh) = (-0.08 - 0.15 * byimf) < phid(ish) = ((9.39 + 0.21 * byimf - 12.) *15.) < phid(inh) = ((9.39 - 0.21 * byimf - 12.) *15.) < phin(ish) = ((23.5 + 0.15 * byimf - 12.) *15.) < phin(inh) = ((23.5 - 0.15 * byimf - 12.) *15.) < pcen(ish) = (-0.168 - 0.027 * byimf) * ctpoten < pcen(inh) = (-0.168 + 0.027 * byimf) * ctpoten < rr1(:) = -2.6 < ! < ! Report auroral parameters in degrees: < ! < if (iprint > 0) then < ! < write(6,"(/'Tail: input variables:')") < write(6,"(' f107 =',e12.4,' f107a=',e12.4)") f107,f107a < write(6,"(' ctpoten=',e12.4)") ctpoten < write(6,"(' power =',e12.4)") power < write(6,"(' byimf =',e12.4)") byimf < ! write(6,"(' fkp =',e12.4)") fkp < ! < write(6,"(/'Tail: /ioncr/ variables (deg south,north):')") < write(6,"(' theta0 = ',2e12.4)") theta0 < write(6,"(' offa = ',2e12.4)") offa < write(6,"(' offc = ',2e12.4)") offc < write(6,"(' dskofa = ',2e12.4)") dskofa < write(6,"(' dskofc = ',2e12.4)") dskofc < write(6,"(' phid = ',2e12.4)") phid < write(6,"(' phin = ',2e12.4)") phin < write(6,"(' phidp0 = ',2e12.4)") phidp0 < write(6,"(' phidm0 = ',2e12.4)") phidm0 < write(6,"(' phinp0 = ',2e12.4)") phinp0 < write(6,"(' phinm0 = ',2e12.4)") phinm0 < write(6,"(' psim = ',2e12.4)") psim < write(6,"(' psie = ',2e12.4)") psie < write(6,"(' pcen = ',2e12.4)") pcen < write(6,"(' ')") < endif < plevel = 0. < if (power >= 0.01) plevel = 2.09 * alog(power) --- > C***************** > IUT = UT + E > UT24 = AMOD(UT,24.) > MDAY = UT/24. > IF (ABS(UT-IUT) .LT. 0.001) WRITE (6,"(1X,'TAIL: UT HP CP BY =', > | F8.2,I3,F8.2,3F7.1/1X,'EC ALFAC FC ED ALFAD FD ALFA21,22 ALFA20 > |RALFA2=',2F7.3,E10.3,2F7.3,E10.3,4F7.3)") UT,MDAY,UT24,HP,CP,BYIMF > |,EC,ALFAC,FC,ED,ALFAD,FD,ALFA21,ALFA22,ALFA20,RALFA2 > IF (IAMIE .EQ. 1) GO TO 190 > C***************** > C OLD ALFA > IF (JSWOLDA .EQ. 1) THEN > C ALFA1 = 0.91 > C ALFA2 = 0.63 + 0.039*HP35 + 0.0022*HPM35 > C ALFA1 = 2.0 > C ALFA2 = 3.0 > C **** SNOE AURORA PARTICLE ALPHA > ALFA1 = 2.0 > ALFA2 = 2.0 > C > TWA6 = 0. > TWA21 = 0. > ELSE > ALFA1 = 0. > ALFA2 = 0. > ENDIF > C > IF (NAURP .GT. 3) THEN > DO 110 N=1,2 > R1(N) = PARAMV(4) > PHIDM0(N) = PARAMV(5) > PHIDP0(N) = PARAMV(6) > PHINM0(N) = PARAMV(7) > 110 PHINP0(N) = PARAMV(8) > ARAD(INH) = PARAMV(9) > ARAD(ISH) = PARAMV(10) > OFFA(INH) = PARAMV(11) > OFFA(ISH) = PARAMV(12) > DSKOFA(INH) = PARAMV(13) > DSKOFA(ISH) = PARAMV(14) > C THIS PART IS NEW (11/10/88) > PHID(INH) = PARAMV(15) > PHID(ISH) = PARAMV(16) > PHIN(INH) = PARAMV(17) > PHIN(ISH) = PARAMV(18) > PSIE(INH) = -CP * PARAMV(19) > PSIM(INH) = CP * (1.-PARAMV(19)) > PSIE(ISH) = -CP * PARAMV(20) > PSIM(ISH) = CP * (1.-PARAMV(20)) > PCEN(INH) = PARAMV(21) * CP > PCEN(ISH) = PARAMV(22) * CP > C CONVERT BACK TO OLDER METHOD TEMPORARILY > C PHID(IHP) = PARAMV(15) > C PHID(IHN) = PARAMV(16) > C PHIN(IHP) = PARAMV(17) > C PHIN(IHN) = PARAMV(18) > C PSIE(IHP) = -CP * PARAMV(19) > C PSIM(IHP) = CP * (1.-PARAMV(19)) > C PSIE(IHN) = -CP * PARAMV(20) > C PSIM(IHN) = CP * (1.-PARAMV(20)) > C PCEN(IHP) = PARAMV(21) * PSIE(IHP) > C PCEN(IHN) = PARAMV(22) * PSIM(IHN) > C END OF CONVERSION > IF (JDIDK .EQ. 0) THEN > DISPC = PARAMV(23) > DSKOFC(INH) = DSKOFA(INH) + DISPC > DSKOFC(ISH) = DSKOFA(ISH) + DISPC > ELSE > DSKOFC(INH) = PARAMV(23) > DSKOFC(ISH) = PARAMV(24) > DISPC = DSKOFC(INH) - DSKOFA(INH) > ENDIF > IF (JDITH .EQ. 0) THEN > DISP = PARAMV(25) > THETA0(INH) = ARAD(INH) - DISP > THETA0(ISH) = ARAD(ISH) - DISP > ELSE > THETA0(INH) = PARAMV(25) > THETA0(ISH) = PARAMV(26) > DISP = ARAD(INH) - PARAMV(25) > ENDIF > OFFC(INH) = PARAMV(27) > OFFC(ISH) = PARAMV(28) > IF(NAURP .GT. 28) THEN > E1 = PARAMV(29) > E2 = PARAMV(30) > H1 = PARAMV(31) > H2 = PARAMV(32) > ROTE = PARAMV(33) > ROTH = PARAMV(34) > IF (JSWOLDA .EQ. 0) THEN > TWA6 = PARAMV(35) > TWA21 = PARAMV(36) > ELSE > C ALFA1 = 2.0 > C ALFA2 = 3.0 > C **** SNOE AURORA PARTICLE ALPHA > ALFA1 = 2.0 > ALFA2 = 2.0 > ENDIF > GO TO 135 > ENDIF ! FOR NAURP .GT. 3 > ELSE > C MADE NEW MODEL IN 3/89 > C ION CONVECTION PARAMETERS: > DSKOFC(INH) = -0.08 - 0.15*BYGD > DSKOFC(ISH) = -0.08 + 0.15*BYGD > PHID(INH) = (9.39 - 0.21 * BYGD - 12.) * 15. > PHID(ISH) = (9.39 + 0.21 * BYGD - 12.) * 15. > PHIN(INH) = (23.50 - 0.15 * BYGD - 12.) * 15. > PHIN(ISH) = (23.50 + 0.15 * BYGD - 12.) * 15. > PCEN(INH) = (-0.168 - 0.027 * BYGD) * CP > PCEN(ISH) = (-0.168 + 0.027 * BYGD) * CP > DO 3010 IH=1,2 > R1(IH) = -2.6 > IF (IAMIE .EQ. 1) THEN > C FORMULA FOR AMIE POTENTIALS > THETA0(IH) = -1.92 + 8.10 * (CP**0.1875) > ELSE > C FORMULA FOR IMF POTENTIALS > THETA0(IH) = -3.80 + 8.48 * (CP**0.1875) > ENDIF > OFFC(IH) = 1.1 > PSIM(IH) = 0.44 * CP > PSIE(IH) = -0.56 * CP > PHIDP0(IH) = 85. > PHIDM0(IH) = 50. > PHINP0(IH) = 57.5 > PHINM0(IH) = 100. > C AURORAL PRECIPITATION > RHP = 14.20 + 0.96*PLEVEL > IF (IAMIE .EQ. 1) THEN > C FORMULA FOR AMIE POTENTIALS > RCP = 3.06 + 8.49 * (CP**0.1875) > ELSE > C FORMULA FOR IMF POTENTIALS > RCP = -0.43 + 9.69 * (CP**0.1875) > ENDIF > ARAD(IH) = AMAX1(RHP,RCP) > 3010 CONTINUE > OFFA(INH) = 3.7 > OFFA(ISH) = 4.3 > DSKOFA(INH) = -1.26 - 0.15 * BYGD > DSKOFA(ISH) = -1.26 + 0.15 * BYGD > DISP = ARAD(1) - THETA0(1) > DISPC = DSKOFC(1) - DSKOFA(1) > ENDIF 125c344 < C E2 = 0.95 + 0.117 * power --- > C E2 = 0.95 + 0.117 * HP 128,140c347,364 < C E1 = (0.95 + 0.117 * power) < C E2 = (0.95 + 0.117 * power) < C E1 = (1.0 + 0.15* power) < C E2 = (1.0 + 0.15* power) < E1 = (1.0 + 0.25* power) < E2 = (1.0 + 0.25* power) < ROTE = (2.62 - 0.55 * PLEVEL) * 15. < C H1 = AMIN1( 2.35, 0.83 + 0.33 * PLEVEL ) < C H2 = 2.87 + 0.15 * PLEVEL < H1 = 3.+0.1*power < H2 = 3.0+0.1*power < < ROTH = (12.18 - 0.89 * PLEVEL) * 15. --- > C E1 = (0.95 + 0.117 * HP) > C E2 = (0.95 + 0.117 * HP) > C E1 = (1.0 + 0.15* HP) > C E2 = (1.0 + 0.15* HP) > E1 = (1.0 + 0.25* HP) > E2 = (1.0 + 0.25* HP) > C H1 = AMIN1( 2.35, 0.83 + 0.33 * PLEVEL ) > C H2 = 2.87 + 0.15 * PLEVEL > H1 = 3.+0.1*HP > H2 = 3.0+0.1*HP > C > write(6,"('tail: ut=',f8.3,' hp=',e12.4,' plevel=',e12.4, > | ' cp=',e12.4)") ut,hp,plevel,cp > write(6,"(' e1=',e12.4,' e2=',e12.4,' h1=',e12.4,' h2=',e12.4)") > | e1,e2,h1,h2 > C > ROTE = (2.62 - 0.55 * PLEVEL) * 15. > ROTH = (12.18 - 0.89 * PLEVEL) * 15. 142c366 < 135 RHP = 14.20 + 0.96*PLEVEL --- > 135 RHP = 14.20 + 0.96*PLEVEL 144c368 < RCP = 3.06 + 8.49 * (ctpoten**0.1875) --- > RCP = 3.06 + 8.49 * (CP**0.1875) 146,149c370,371 < RCP = -0.43 + 9.69 * (ctpoten**0.1875) < RAUR = max(RCP,RHP) < arad(1) = max(rcp,rhp) < arad(2) = max(rcp,rhp) --- > RCP = -0.43 + 9.69 * (CP**0.1875) > RAUR = AMAX1(RCP,RHP) 151,173c373,376 < TWAK = 0.50 < TWA6 = 0. < TWA21 = 0. < ROT6 = (6.00 - 12.) * 15. < ROT21 = (21.00 - 12.) * 15. < T6 = 7.00 < T21 = 4.00 < D6 = -4.0 < D6V = 0. < D21 = 4.0 < H6 = 7.0 < H21 = 10.0 < E21 = 1.e-80 < E22 = 1.e-80 < IF (power .LT. 0.01) THEN < E1 = 1.E-20 < E2 = 1.E-20 < EC = 1.E-20 < ED = 1.E-20 < E21 = 1.E-20 < E22 = 1.E-20 < FD = 1.E-20 < FC = 1.E-20 --- > TWAK = 0.50 > IF (NAURP .LE. 28 .AND. JSWOLDA .EQ. 0) THEN > TWA6 = 0.36 + 0.48 * PLEVEL > TWA21 = AMAX1( 1.00, -1.75 + 0.69 * PLEVEL ) 174a378,397 > ROT6 = (6.00 - 12.) * 15. > ROT21 = (21.00 - 12.) * 15. > T6 = 7.00 > T21 = 4.00 > D6 = -4.0 > D6V = 0. > D21 = 4.0 > H6 = 7.0 > H21 = 10.0 > 140 CONTINUE > IF (HP .LT. 0.01) THEN > E1 = 1.E-20 > E2 = 1.E-20 > EC = 1.E-20 > ED = 1.E-20 > E21 = 1.E-20 > E22 = 1.E-20 > FD = 1.E-20 > FC = 1.E-20 > ENDIF 182,184d404 < ALFA1 = 2. < ! ALFA2 = 3. < ALFA2 = 2. 187a408,420 > IUT = UT + E > C IF (ABS(UT-IUT) .LT. 0.001) WRITE (6,603) UT, HP, CP, BYIMF, > C | RHP, RCP, RAUR, (THETA0(I),I=1,51), DISPC,ARAD(1),ARAD(2), > C | TWA6,TWA21,JSWOLDA,JDIDK,JDITH > C 603 FORMAT(1X,'TAIL: ',4X,'UT HP CP BYIMF RHP RCP RAUR =', 7F6.1/ > C |1X,'PARAMD=THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),DSKOFC(2),PHID(2)'/ > C |1X,' PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),PHINM0(2),PSIM(2)'/ > C |1X,' PSIE(2),PCEN(2),ARAD(2),H1,H2,ROTH,E1,E2,ROTE,EC,ED,ALFA1'/ > C |1X,' ALFA2,TWAK,TWA6,TWA21,ROT6,ROT21,D6,D21,H6,H21,T6,T21'/ > C | 5X, 'PARAMD(1-10) =', 10E10.3/ 5X, 'PARAMD(11-20)=',10E10.3/ > C | 5X, 'PARAMD(21-30)=', 10E10.3/ 5X, 'PARAMD(31-40)=',10E10.3/ > C | 5X, 'PARAMD(41-50)=', 10E10.3/ 5X, 'PARAMD(51)=',E10.3/ > C | 5X, 'DISPC ARAD(1,2) TWA6 TWA21 JSWOLDA,DIDK,DITH =', 5F6.3,3I3) 211,227d443 < ! < ! Convert to radians: < theta0(:) = theta0(:) * dtr < offa(:) = offa(:) * dtr < offc(:) = offc(:) * dtr < dskofa(:) = dskofa(:) * dtr < dskofc(:) = dskofc(:) * dtr < phid(:) = phid(:) * dtr < phin(:) = phin(:) * dtr < phidp0(:) = phidp0(:) * dtr < phidm0(:) = phidm0(:) * dtr < phinp0(:) = phinp0(:) * dtr < phinm0(:) = phinm0(:) * dtr < ! < psim(:) = psim(:) * 1000. < psie(:) = psie(:) * 1000. < pcen(:) = pcen(:) * 1000. 228a445,449 > DO 200 N=1,22 > 200 RDEG(N) = THETA0(N) * DTR > DO 210 N=1,6 > 210 PSIV(N) = PSIM(N) * 1000. > C 242,271c463,465 < ! < ! COMMON /OVALR/ RRAD(2),H0,RH,RROTH,E0,REE,RROTE,FC,ALFAC,FD,ALFAD, < ! | ALFK,ALF6,ALF21,RROT6,RROT21,RD6,RD6V,RD21,RH6,RH21,RT6,RT21, < ! | ALFA0,RALFA,ALFA20,RALFA2,E20,RE2, ALFA30,E30,FD2,ALFAD2,ED2 < ! < if (iprint > 0) then < write(6,"('Tail: /ovalr/ variables:')") < write(6,"(' rrad =',2e12.4)") rrad < write(6,"(' h0 =',e12.4)") h0 < write(6,"(' rh =',e12.4)") rh < write(6,"(' rroth =',e12.4)") rroth < write(6,"(' e0 =',e12.4,' (e1,e2=',2e12.4,')')") e0,e1,e2 < write(6,"(' ree =',e12.4)") ree < write(6,"(' rrote =',e12.4)") rrote < write(6,"(' fc =',e12.4,' fd =',e12.4)") fc,fd < write(6,"(' alfac =',e12.4,' alfad =',e12.4,' alfk=',e12.4)") < | alfac,alfad,alfk < write(6,"(' alf6 =',e12.4,' alf21 =',e12.4)") alf6,alf21 < write(6,"(' rrot6 =',e12.4,' rrot21=',e12.4)") rrot6,rrot21 < write(6,"(' rd6 =',e12.4,' rd6v =',e12.4)") rd6,rd6v < write(6,"(' rd21 =',e12.4)") rd21 < write(6,"(' rh6 =',e12.4,' rh21 =',e12.4)") rh6,rh21 < write(6,"(' rt6 =',e12.4,' rt21 =',e12.4)") rt6,rt21 < write(6,"(' alfa0 =',e12.4,' alfa20=',e12.4)") alfa0,alfa20 < write(6,"(' ralfa =',e12.4,' ralfa2=',e12.4)") ralfa,ralfa2 < write(6,"(' e20 =',e12.4)") e20 < write(6,"(' re2 =',e12.4)") re2 < write(6,"(' ')") < endif < end subroutine tail --- > RETURN > END > C ======================================================================== Diff of /home/tgcm/tgcm14a/vdrift2.f and vdrift2.f: 10c10 < include "fcom.h" --- > include "blnk.h" 21d20 < !$OMP THREADPRIVATE (/vscr/) 26c25 < character(len=80) :: title --- > ! 30,33d28 < < ! write(6,"('vdrift2: j=',i3,' k=',i3,' f(:,nzk)=',/,(6e12.4))") < ! | j,k,f(:,nzk) < 56c51 < EEZ(I,K) = EZ(I,J,K)/(F(I+2,NZK+1)-F(I+2,NZK-1)) --- > EEZ(I,K) = EZ(I,J,K)/(F(I+2,NZK+1)-F(I+2,NZK-1)) 63,64c58,59 < EEZ(I,1) = 2.*EEZ(I,2)-EEZ(I,3) < EEZ(I,KMAXP1) = 2.*EEZ(I,KMAX)-EEZ(I,KMAX-1) --- > EEZ(I,1) = 2.*EEZ(I,2)-EEZ(I,3) > EEZ(I,KMAXP1) = 2.*EEZ(I,KMAX)-EEZ(I,KMAX-1) ======================================================================== Diff of /home/tgcm/tgcm14a/xray.f and xray.f: 2,3d1 < use input_module,only: f107 < use init_module,only: sfeps 10c8 < include "fcom.h" --- > include "blnk.h" 21a20,23 > ! 11/30/99: add SFEPS (cons.h): > ! SFEPS is initialized to 1 in con.f, and recalculated in > ! advnce.f only if calendar day is being advanced and at > ! new day boundaries. 23,26c25,27 < ! < ! **** SNOE X-RAYS < ! EX=0.3+0.5*(f107-67.)/(176.*SF) < EX=(0.3+0.5*(f107-70.)/240.)*SFEPS --- > C **** SNOE X-RAYS > EX=(0.3+0.5*(C(61)-70.)/240.)*SFEPS > C EX=0.6+0.5*(C(61)-67.)/(176.*SF)