Fri Feb 7 13:32:22 US/Mountain 2003 Diffsrc executing from /home/foster/tiegcm/tgcm15/modsrc_amie Source directory = /home/tgcm/tgcm15 ======================================================================== Diff of /home/tgcm/tgcm15/dynphi.h and dynphi.h: 5c5 < real :: dynpot,phim,phim3d,phih,phihm,zzm,ex,ey,ez --- > real :: dynpot,phim,phim3d,phih,phihm,zzm,ex,ey,ez,phisym,phiant 10c10,11 < 3 EZ(IMAXGP,JMAXG,ZKMXP) --- > 3 EZ(IMAXGP,JMAXG,ZKMXP), > 4 PHISYM(IMAXMP,JMX0),PHIANT(IMAXMP,JMAXM) ======================================================================== WARNING: cannot find file /home/tgcm/tgcm15/extrafields.h ======================================================================== Diff of /home/tgcm/tgcm15/adheel.F and adheel.F: 17c17,99 < integer :: ncc,imx,jmx,n --- > integer :: ncc,imx,jmx,n,jj,j0,j,i > > ! crit(1) = 0.261799387 > ! crit(2) = 0.523598775 > C **** > C **** FILL ARRAY COLATC WITH NORTHERN AURORAL COLATITUDES > C **** CORRESPONDING TO EACH NORTHERN GEOMAGNETIC GRID > C **** POINT. VALUES IN COLATC(IMX0,JMX0). > C **** ALSO CALCULATE FRACTIONAL PRESENCE OF DYNAMO > C **** EQUATION GIVEN CRITICAL COLATITUDES CRIT(2). > C **** VALUES IN P(IMX0,JMX0). > C **** > CALL COLATH > C **** > C **** CALCULATE PHIHM, THE HEELIS POTENTIAL IN > C **** GEOMAGNETIC COORDINATES. > C **** > istar = 0 > CALL POTM > C **** > C **** MODIFY STENCILS > C **** > ! NCC = 1 > ! IMX = IMX0 > ! JMX = JMX0 > ! DO 5 N = 1,5 > ! CALL STENMOD(IMX,JMX,CEE(NCC),PHIHM(1,JMX0),P,COLATC) > ! NCC = NCC+9*IMX*JMX > ! IF(N.EQ.1)NCC = NCC+IMX*JMX > ! IMX = (IMX+1)/2 > ! JMX = (JMX+1)/2 > ! 5 CONTINUE > > C **** > C **** Calculate the symmetric and antisymmetric parts of PHIHM > C **** Put the symmetric part in as if it were the NH part (equ-pole) > C **** (IMX0=IMAXMP=IMAXM+1=80+1=81, JMXO=(JMAXM+1)/2=49, JMAXM=97) > C **** > DO 100 I=1,IMAXMP > PHISYM(I,1)=PHIHM(I,JMX0) > PHIANT(I,JMX0) = 0. > DO 100 J=1,JMX0-1 > JJ = JMAXM + 1 - J > J0 = JMX0 - J + 1 > PHISYM(I,J0) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) > PHIANT(I,J) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) - PHIHM(I,J) > 100 PHIANT(I,JJ) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) - PHIHM(I,JJ) > C Plot of symmetric hi-lat potential > C CALL EZCNTR (PHISYM,IMX0,JMX0) > C Plot of anti-symmetric hi-lat potential > C CALL EZCNTR (PHIANT,IMX0,JMAXM) > C **** > C **** MODIFY STENCILS > C **** > NCC = 1 > IMX = IMX0 > JMX = JMX0 > DO 5 N = 1,5 > CALL STENMOD(IMX,JMX,CEE(NCC),PHISYM,P,COLATC) > NCC = NCC+9*IMX*JMX > IF(N.EQ.1)NCC = NCC+IMX*JMX > IMX = (IMX+1)/2 > JMX = (JMX+1)/2 > 5 CONTINUE > C P=0 from pole to crit1, P=1 from crit2 to equator > DO 200 I=1,IMAXMP > DO 200 J=1,JMX0-1 > JJ = JMAXM + 1 - J > J0 = JMX0 - J + 1 > PHIANT(I,J) = PHIANT(I,J) * (1. - P(I,J0)) > 200 PHIANT(I,JJ) = PHIANT(I,JJ) * (1. - P(I,J0)) > > RETURN > END > C > C-------------------------------------------------------------------------- > #include "dims.h" > SUBROUTINE ADHEELMD > ! > ! am_02/02: adheelmd: adheel for modified mudpack solver modmud > ! (isolve = 2 in dyncal) > ! subroutine stenmd: calculate modified and unmodified > ! coefficient stencil 19,20c101 < crit(1) = 0.261799387 < crit(2) = 0.523598775 --- > implicit none 21a103,120 > C **** MODIFY STENCILS AND RHS SO AS TO INSERT HEELIS > C **** POTENTIAL FOR GEOMAGNETIC LATITUDE .GE. 60.0 DEG. > C **** > C **** > #include "params.h" > #include "ceee.h" > #include "cterp.h" > #include "dynphi.h" > #include "consdyn.h" > #include "ioncr.h" > ! > ! Locals: > integer :: ncc,imx,jmx,n,jj,j0,j,i > ! > ! crit(1) = 0.261799387 > ! crit(2) = 0.523598775 > ! > C **** 38a138,170 > ! NCC = 1 > ! IMX = IMX0 > ! JMX = JMX0 > ! DO 5 N = 1,5 > ! CALL STENMD(IMX,JMX,CEE(NCC),PHIHM(1,JMX0),P,COLATC) > ! NCC = NCC+9*IMX*JMX > ! IF(N.EQ.1)NCC = NCC+IMX*JMX > ! IMX = (IMX+1)/2 > ! JMX = (JMX+1)/2 > ! 5 CONTINUE > > > C **** > C **** Calculate the symmetric and antisymmetric parts of PHIHM > C **** Put the symmetric part in as if it were the NH part (equ-pole) > C **** (IMX0=IMAXMP=IMAXM+1=80+1=81, JMXO=(JMAXM+1)/2=49, JMAXM=97) > C **** > DO 100 I=1,IMAXMP > PHISYM(I,1)=PHIHM(I,JMX0) > PHIANT(I,JMX0) = 0. > DO 100 J=1,JMX0-1 > JJ = JMAXM + 1 - J > J0 = JMX0 - J + 1 > PHISYM(I,J0) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) > PHIANT(I,J) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) - PHIHM(I,J) > 100 PHIANT(I,JJ) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) - PHIHM(I,JJ) > C Plot of symmetric hi-lat potential > C CALL EZCNTR (PHISYM,IMX0,JMX0) > C Plot of anti-symmetric hi-lat potential > C CALL EZCNTR (PHIANT,IMX0,JMAXM) > C **** > C **** MODIFY STENCILS > C **** 43c175 < CALL STENMOD(IMX,JMX,CEE(NCC),PHIHM(1,JMX0),P,COLATC) --- > CALL STENMOD(IMX,JMX,CEE(NCC),PHISYM,P,COLATC) 48a181,187 > C P=0 from pole to crit1, P=1 from crit2 to equator > DO 200 I=1,IMAXMP > DO 200 J=1,JMX0-1 > JJ = JMAXM + 1 - J > J0 = JMX0 - J + 1 > PHIANT(I,J) = PHIANT(I,J) * (1. - P(I,J0)) > 200 PHIANT(I,JJ) = PHIANT(I,JJ) * (1. - P(I,J0)) ======================================================================== Diff of /home/tgcm/tgcm15/advnce.F and advnce.F: 7c7,8 < use input_module,only: calday,start,step,f107,f107a,power,ctpoten --- > use input_module,only: calday,start,step,f107,f107a, > | power,ctpoten,iamie,amie_ibkg 12a14 > use amie_module,only: getamie 25d26 < #include "amie.h" 86d86 < sectgcm = float(iter)*dt 88c88,93 < if (iamie==1) call setamie(sec) --- > > ! Get AMIE data > if (iprint>0) write(6,"('advnce before getamie')") > if (iamie==1) then > call getamie(iyear,iday,int(secs),amie_ibkg,iprint) > endif 365a371 > ! write(6,"('calling mp_snddyn: jx=',i3)") jx 373a380,385 > ! gl - 23/07/2002 > ! Calculate FAC and save extra secondary history fields > ! (qamie,qwind,famie,fwind,wamie,wqind,wtot, etc.) > ! print *,'Calling extraoutput: j = ',j,'ixtimep =',ixtimep > if (mytid==0) call extraoutput > ! 391c403 < if (wrsech) call mp_sechtomaster --- > ! if (wrsech) call mp_sechtomaster 485,522d496 < subroutine setamie(sec) < implicit none < #include "amie.h" < ! < ! Args: < real,intent(in) :: sec < ! < ! Local: < integer :: ja,ia,ka < ! < 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 (iuamie) 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 < end subroutine setamie < !----------------------------------------------------------------------- ======================================================================== WARNING: cannot find file /home/tgcm/tgcm15/amie_mod.F ======================================================================== Diff of /home/tgcm/tgcm15/amiepa.F and amiepa.F: 2c2,3 < SUBROUTINE AMIEPA (VELNS,VELEW,POTEN,CUSP,ALFA,FLUX,DRIZL,WION) --- > SUBROUTINE AMIEPA (VELNS,VELEW,POTEN,CUSP,ALFA1,ALFA2, > + FLUX1,FLUX2,DRIZL,WION) 3a5 > use amie_module,only: ekvg, efxg, crad, phida 4a7,9 > > ! gl 15/07/2002 - modified by taking out the unused AMIE parameters > ! like VELNS, VELEW, POTEN, WION 13d17 < #include "amie.h" 17,18c21,24 < | ALFA(ZIMXP,2),FLUX(ZIMXP,2),DRIZL(ZIMX),CUSP(ZIMX),DLAT(ZIMX), < | DLON(ZIMX) --- > ! | ALFA(ZIMXP,2),FLUX(ZIMXP,2),DRIZL(ZIMX),CUSP(ZIMX),DLAT(ZIMX), > ! | DLON(ZIMX) > | ALFA1(ZIMX),ALFA2(ZIMX),FLUX1(ZIMX),FLUX2(ZIMX), > | DRIZL(ZIMX),CUSP(ZIMX),DLAT(ZIMX),DLON(ZIMX) 21,22c27,31 < DATA S5/0.08726646/, S10/0.174532925/, S20/0.34906585/ < DATA PI/3.14159265358979/ --- > > s5 = 0.08726646 > s10 = 0.174532925 > s20 = 0.34906585 > pi = 3.14159265358979 27,29c36,40 < DO 1000 I=1,IMAX < DLAT(I) = RLATM(I+2,J) < 1000 DLON(I) = RLONM(I+2,J)-DLONS(J) --- > DO I=1,IMAX > DLAT(I) = RLATM(I+2,J) > DLON(I) = RLONM(I+2,J)-DLONS(J) > ENDDO > 33,40c44,53 < DO 2000 I=1,IMAX < SINLAT(I) = SIN(ABS(DLAT(I))) < COSLAT(I) = COS(DLAT(I)) < SINLON(I) = SIN(DLON(I)) < COSLON(I) = COS(DLON(I)) < COLAT(I) = ACOS( SINLAT(I) ) < 2000 ALON(I) = AMOD(ATAN2(+SINLON(I)*COSLAT(I),COSLAT(I)*COSLON(I)) < | +3.*PI,2.*PI) - PI --- > DO I=1,IMAX > SINLAT(I) = SIN(ABS(DLAT(I))) > COSLAT(I) = COS(DLAT(I)) > SINLON(I) = SIN(DLON(I)) > COSLON(I) = COS(DLON(I)) > COLAT(I) = ACOS( SINLAT(I) ) > ALON(I) = AMOD(ATAN2(+SINLON(I)*COSLAT(I),COSLAT(I)*COSLON(I)) > | +3.*PI,2.*PI) - PI > ENDDO > 45,84c58,83 < C EFLXA IN W/M2, NFLUX ASSUMES IT IS IN ERG/CM2-S < SECA = AMAX1(SECUTA(1),SECTGCM) < SECA = AMIN1(SECUTA(2),SECA) < F1 = SECA - SECUTA(1) < F2 = SECUTA(2) - SECA < DO 3000 I=1,IMAX < POTEN(I) = (F2*POTKVA1(I,J)+F1*POTKVA2(I,J))*1000./FLOAT(NSPSEC) < VELEW(I) = (F2*VIXYZA1(I,J,1)+F1*VIXYZA2(I,J,1)) * 100. / < | FLOAT(NSPSEC) < VELNS(I) = (F2*VIXYZA1(I,J,2)+F1*VIXYZA2(I,J,2)) * 100. / < | FLOAT(NSPSEC) < WION(I) = (F2*VIXYZA1(I,J,3)+F1*VIXYZA2(I,J,3)) * 100. / < | FLOAT(NSPSEC) < ALFA(I,1) = (F2*EKEVA1(I,J)+F1*EKEVA2(I,J))*0.5/FLOAT(NSPSEC) < C WARNING --- THIS WORKS ONLY BECAUSE ALFA21=ALFA22=ALFA20 < ALFA(I,2) = ALFA20 < C CUSP AND DRIZL ARE FRACTIONS, NUMBERS BETWEEN 0 AND 1. < C ON RECOMMENDATION OFROD HEELIS, PLACE THE CUSP AT THE REVERSAL < C BOUNDARY WITH A RADIUS OF 5 DEGREES < CUSP(I)=(EXP(-((CRAD(IH)-COLAT(I))/S5)**2)+EXP(-((PI-CRAD(IH)- < |COLAT(I))/S5)**2))*EXP(-(ATAN2(SIN(ALON(I)-PHIDA(IH)),COS(ALON(I) < | - PHIDA(IH)))/S20)**2) < 3000 DRIZL(I) = EXP(-(( COLAT(I)-CRAD(IH) + ABS(COLAT(I)-CRAD(IH)) ) < | / S10)**2) < DO 4000 I=1,IMAX < FLUX(I,1) = (F2*EFLXA1(I,J)+F1*EFLXA2(I,J))*1000. / (2.*ALFA(I,1) < | * 1.602E-9 * FLOAT(NSPSEC) ) < C 7/92 The following flux is much too large, and immediately develops < C large Te's and low Ne's in the auroral region < C FLUX(I,2) = (F2*EFLXA1(I,J)+F1*EFLXA2(I,J))*1000. / (2.*ALFA20 < C | * 1.602E-9 * FLOAT(NSPSEC) ) < C WK1 = COS(LAMDA) < WK1(I) = COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE))) < C WK2 = auroral half-width < WK2(I) = H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH),COS(ALON(I)- < | RROTH))) ) < C WK3 = colat - rad < WK3(I) = COLAT(I) - CRAD(IH) < 4000 FLUX(I,2) = E20*(1.-RE2*WK1(I)) * EXP(-(WK3(I)/WK2(I))**2) < | / (2. * ALFA20 * 1.602E-9) --- > C EFLXA IN mW/M2, NFLUX ASSUMES IT IS IN ERG/CM2-S > DO I=1,IMAX > ! print *,'AMIEPA: i,j,ekvg(i,j),efxg(i,j)=',i,j,ekvg(i,j),efxg(i,j) > if (ekvg(i,j) .lt. 1.0) then > ! write(6,"('amiepa: negative mean energy at i=',i3, > ! | 'j=',i3,e12.4)") i,j,ekvg(i,j) > ekvg(i,j) = 1.0 > endif > ALFA1(I) = ekvg(I,J)/2. > ALFA2(I) = ALFA20 > CUSP(I)=(EXP(-((CRAD(IH)-COLAT(I))/S5)**2)+EXP(-((PI-CRAD(IH)- > | COLAT(I))/S5)**2))*EXP(-(ATAN2(SIN(ALON(I)-PHIDA(IH)), > | COS(ALON(I) - PHIDA(IH)))/S20)**2) > DRIZL(I) = EXP(-(( COLAT(I)-CRAD(IH) + ABS(COLAT(I)-CRAD(IH)) ) > | / S10)**2) > > FLUX1(I) = AMAX1(efxg(I,J)/(2.*ALFA1(I) * 1.602E-9),1.E-20) > WK1(I) = COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE))) > WK2(I) = H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH),COS(ALON(I)- > | RROTH))) ) > WK3(I) = COLAT(I) - CRAD(IH) > FLUX2(I) = E20*(1.-RE2*WK1(I)) * EXP(-(WK3(I)/WK2(I))**2) > | / (2. * ALFA20 * 1.602E-9) > ENDDO > ! print *,'End of amiepa >>>' > ======================================================================== Diff of /home/tgcm/tgcm15/dt.F and dt.F: 396,400c396,399 < ! do k=1,kmaxp1 < ! s4(imax+3,k) = s4(3,k) < ! enddo < ! call addfsech('QJOULE','JOULE HEATING from DT','DEG K', < ! | s4,zimxp,zkmxp,zkmx,j) --- > do k=1,kmaxp1 > s4(imax+3,k) = s4(3,k) > enddo > call addfsech('QJOULE',' ',' ',s4,zimxp,zkmxp,zkmx,j) ======================================================================== Diff of /home/tgcm/tgcm15/dyncal.F and dyncal.F: 145c145,146 < PHIM(I,J) = RIM(I,J,1) --- > ! PHIM(I,J) = RIM(I,J,1) > PHIM(I,J) = RIM(I,J,1) - PHIANT(I,J) ======================================================================== Diff of /home/tgcm/tgcm15/elden.F and elden.F: 4c4 < | boltz,p0 --- > | boltz,p0,imax 18a19 > #include "extrafields.h" 21c22 < integer :: i,k --- > integer :: i,k,nzk 23a25,26 > real :: tec(zimxp1) > 161a165 > if (F(I,NN2PK) < 0.0) F(I,NN2PK) = 0. 163a168 > if (F(I,NO2PK) < 0.0) F(I,NO2PK) = 0. 167a173 > if (F(I,NNOPK) < 0.0) F(I,NNOPK) = 0. 168a175,189 > C len2=kmax*(imax+4) so i should be i+2 in S5 > C Get ht integral of electron density and save in tec of /sight/ > tec_sec(:,j,:) = 0. > tec(:) = 0. > nzk = nj+nz-1 > do k=1,kmax > nzk = nzk+1 > do i=1,imax+1 > tec_sec(i,j,k) = s5(i+2,k) > tec(i) = tec(i) + > + (f(i,nzk+1)-f(i,nzk)) * s5(i+2,k) > enddo > enddo > tec_sec(:,j,zkmxp) = tec(:) > ======================================================================== WARNING: cannot find file /home/tgcm/tgcm15/extraoutput.F ======================================================================== Diff of /home/tgcm/tgcm15/heelis.F and heelis.F: 3c3 < use input_module,only: aurora --- > use input_module,only: aurora, iamie 15c15 < #include "amie.h" --- > !#include "amie.h" 52,54c52,56 < write(6,"('heelis calling amiepa')") < CALL AMIEPA (T2(3),T1(3),T4(3),T5(3),T6(3),T8(3),T3(3), < | WION(3)) --- > ! write(6,"('heelis calling amiepa')") > ! CALL AMIEPA (T2(3),T1(3),T4(3),T5(3),T6(3),T8(3),T3(3), > ! | WION(3)) > CALL AMIEPA (T2(3),T1(3),T4(3),T5(3),T6(3),T7(3), > | T8(3),T9(3),T3(3),WION(3)) 94,99c96,101 < ! 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 --- > ! write(6,"('t3(DRIZL)=',/,(6e12.4))") t3 > ! write(6,"('t5(CUSP)=',/,(6e12.4))") t5 > ! write(6,"('t6(ALFA)=',/,(6e12.4))") t6 > ! write(6,"('t7(ALFA2)=',/,(6e12.4))") t7 > ! write(6,"('t8(FLUX)=',/,(6e12.4))") t8 > ! write(6,"('t9(FLUX2)=',/,(6e12.4))") t9 ======================================================================== Diff of /home/tgcm/tgcm15/init_amie.F and init_amie.F: 4c4,8 < use input_module,only: tempdir,amievol --- > ! gl - 22/07/2002 > ! this routine is not called any more - instead by rdamie_nh and rdamie_sh > ! also amie.h has been taken out of all relavent subroutines > ! > use input_module,only: tempdir,amievol,iamie 6d9 < #include "amie.h" 9c12 < integer :: ios --- > integer :: ios, iuamie ======================================================================== Diff of /home/tgcm/tgcm15/init_mod.F and init_mod.F: 70c70 < | nmc,source_start --- > | nmc,source_start,amienh,amiesh 76a77 > use amie_module,only: rdamie_nh,rdamie_sh 160c161,170 < call init_amie --- > ! call init_amie > if (len_trim(amienh) > 0) then > write(6,"('Reading AMIENH = ',a)") trim(amienh) > call rdamie_nh > endif > if (len_trim(amiesh) > 0) then > write(6,"('Reading AMIESH = ',a)") trim(amiesh) > call rdamie_sh > endif > ======================================================================== Diff of /home/tgcm/tgcm15/input_mod.F and input_mod.F: 2a3,4 > ! am_02/02: added general wave (origin Kelvin wave) > 42c44,46 < | amievol ! file or mss path of amie data file (optional) --- > | amievol, ! file or mss path of amie SH data file (optional) > | amiesh, ! file or mss path of amie SH data file (optional) > | amienh ! file or mss path of amie NH data file (optional) 52c56,58 < | aurora ! 0/1 flag for aurora --- > | aurora, ! 0/1 flag for aurora > | iamie, ! 0/1 flag for AMIE data > | amie_ibkg ! 0/1/2 flag for read real, 1st, or 24-hr averaged data 58a65,67 > ! cg_kw am_1200 > | wave(6), ! general wave: typ,period,amplitude,phase,h_n,zonal wave no > ! cg_kw am_1200 109c118 < | label, tempdir, magvol, amievol, gpi_ncfile --- > | label, tempdir, magvol, amievol, amiesh, amienh, gpi_ncfile 112c121,123 < | dispose, difhor, iuivi, nmc, tideann, aurora --- > | dispose, difhor, iuivi, nmc, tideann, aurora, iamie, > ! 0/1/2 flag for read real, 1st, or 24-hr averaged data > | amie_ibkg 114c125 < | mag(2,2), tide(10),tide2(2),tide3m3(2), --- > | mag(2,2), tide(10),tide2(2),tide3m3(2),wave(6), 148,149c159,160 < | label,tempdir,magvol,amievol,date,calday,step,dispose, < | source,source_start,output,start,stop,hist,save, --- > | label,tempdir,magvol,amievol,amiesh,amienh,date,calday,step, > | dispose,source,source_start,output,start,stop,hist,save, 151,153c162,164 < | mag,difhor,iuivi,nmc,tide,tide2,tide3m3,f107,f107a, < | power,ctpoten,byimf,colfac,tideann,aurora,gpi_ncfile, < | mxhist_prim,mxhist_sech,msreten --- > | mag,difhor,iuivi,nmc,tide,tide2,tide3m3,wave,f107,f107a, > | power,ctpoten,byimf,colfac,tideann,aurora,iamie,gpi_ncfile, > | mxhist_prim,mxhist_sech,msreten,amie_ibkg 200a212,213 > amiesh = ' ' > amienh = ' ' 219a233 > wave = spval 220a235,236 > iamie = 0 ! default is no AMIE data > amie_ibkg = 0 ! default is read in real time AMIE data 253c269 < integer :: i,n --- > integer :: i,n,ictr 420a437,471 > ! cg_kw am_1200 > ! Wave: > n = size(wave)-count(wave==spval) > if (n /= 6) then > write(6,"(/,'>>> INPUT: must have 6 entries for WAVE,', > | ' got ',i3)") n > stop 'WAVE' > endif > ictr = 0 > if(wave(1) < 0.) then > write(6,"(/,'>>> INPUT: typ for WAVE(1) must ', > | 'be > 0.: wave_typ=',e12.4)") wave(1) > ictr = ictr +1 > elseif(wave(1) > 0.) then > if(wave(2) == 0.) then > write(6,"(/,'>>> INPUT: period for WAVE(2) must ', > | 'be > or < 0.: wave_period=',e12.4)") wave(2) > ictr = ictr +1 > elseif (wave(3) <= 0.) then > write(6,"(/,'>>> INPUT: amplitude for WAVE(3) must ', > | 'be > 0.: wave_amp=',e12.4)") wave(3) > ictr = ictr +1 > elseif(wave(5) <= 0.) then > write(6,"(/,'>>> INPUT: equival. depth for WAVE(5) must ', > | 'be > 0.: wave_hequ=',e12.4)") wave(5) > ictr = ictr +1 > elseif(wave(6) <= 0.) then > write(6,"(/,'>>> INPUT: zonal wave no. for WAVE(6) must ', > | 'be > 0.: wave_zono.=',e12.4)") wave(5) > ictr = ictr +1 > endif > endif > if(ictr /= 0) stop 'WAVE' > ! cg_kw am_1200 > ! 570a622,632 > ! > ! Read AMIE data for either NH or SH or Both > ! gl - 07/12/2002 > ! > n = len_trim(amiesh) + len_trim(amienh) > if (n > 0) then > write(6,"('Read AMIE Data as inputs')") > iamie = 1 > endif > ! > ! NMC boundaries for T and Z are not supported in tiegcm: 1353a1416,1417 > inp%amiesh = amiesh > inp%amienh = amienh 1366a1431 > inp%wave(:) = wave(:) 1367a1433,1434 > inp%iamie = iamie > inp%amie_ibkg = amie_ibkg 1420a1488,1495 > if (len_trim(amiesh) > 0) > | write(6,"(' amiesh = ',a,/,4x, > | '(file or mss path containing amie SH data)')") > | trim(inp%amiesh) > if (len_trim(amienh) > 0) > | write(6,"(' amienh = ',a,/,4x, > | '(file or mss path containing amie NH data)')") > | trim(inp%amienh) 1513a1589,1591 > write(6,"(' wave (typ, period, amplitude, phase, ', > | 'equi. depth, zonal wave no.)=', > | /,4x,2f6.2,e8.1,x,f7.2,x,e8.1,x,f7.2))") inp%wave ======================================================================== Diff of /home/tgcm/tgcm15/lamdas.F and lamdas.F: 1a2,3 > ! am_02/02: option to change electron-neutral collision frequency > 20a23 > #include "extrafields.h" 38c41,45 < | ntek,nlxxk,nlyyk,nlxyk,nzk,nuk,nvk,nwk --- > | ntek,nlxxk,nlyyk,nlxyk,nzk,nuk,nvk,nwk,nuik,nvik > real :: tm1(zimxp),tm2(zimxp),tm3(zimxp),sndip(zimxp) > real :: sumfaus, sumfavs, sumfwus, sumfwvs, sumfaun, sumfavn, > | sumfwun, sumfwvn, qwind0, qamie0, wtot0, fwindu0, fwindv0, > | famieu0, famiev0, work0 86a94,96 > if (RM1(I,1,L,M) .lt. 0.) > | print *,'LAMDAS: negative RM1 for J,I,M,L = ', > | J,I,M,L,RM1(I,1,L,M),V6(I,1,1),F(I,NTIK),F(I,NTK) 107a118,119 > if (RM1(I,1,1,1) .lt. 0.) > | print *,'LAMDAS: negative O2+',I,J,F(I,NO2PK) 108a121,122 > if (RM1(I,1,2,1) .lt. 0.) > | print *,'LAMDAS: negative O+',I,J,F(I,NOPK) 109a124,125 > if (RM1(I,1,3,1) .lt. 0.) > | print *,'LAMDAS: negative NO+',I,J,F(I,NNOPK) 142a159,167 > > ! am 2001-10-2 modifying the electron-neutral collision frequency > ! paper: Atmos. Electrodynamics > ! V4(I,1,1)=V6(I,1,1)*(7.2E-15*v4tmp*(F(I,NTEK)/300.)**0.95 > ! | +5.2E-15*F(I,NPS1K)*(F(I,NTEK)/300.)**0.79*rmassinv(1) > ! | +1.9E-15*F(I,NPS2K)*rmassinv(2)*(F(I,NTEK)/300.)**0.85) > ! | *1E6 ! convertion from m^3 to cm^3 > ! > ! original Banks & Kockards 146a172,173 > ! end am 2001-10-2 > 152a180,181 > ! am 2001-10-17 multiply collision frequency by factor of 4. > ! V4(I,K,1)=4.*V4(I,K,1)/T3(I) 190a220,223 > if (V3(I,1,1)*V3(I,2,1) .lt. 0.) print *, > | 'lamdas: I,J,V3(I,1,1),V3(I,2,1) = ',I,J,V3(I,1,1),V3(I,2,1) > if (V3(I,1,1) .lt. 0.) V3(I,1,1) = 0. > if (V3(I,2,1) .lt. 0.) V3(I,2,1) = 0. 191a225,228 > if (V4(I,1,1)*V4(I,2,1) .lt. 0.) print *, > | 'lamdas: I,J,V4(I,1,1),V4(I,2,1) = ',I,J,V4(I,1,1),V4(I,2,1) > if (V4(I,1,1) .lt. 0.) V4(I,1,1) = 0. > if (V4(I,2,1) .lt. 0.) V4(I,2,1) = 0. 192a230,233 > if (V5(I,1,1)*V5(I,2,1) .lt. 0.) print *, > | 'lamdas: I,J,V5(I,1,1),V5(I,2,1) = ',I,J,V5(I,1,1),V5(I,2,1) > if (V5(I,1,1) .lt. 0.) V5(I,1,1) = 0. > if (V5(I,2,1) .lt. 0.) V5(I,2,1) = 0. 255a297,482 > > ! > ! add additional fields calculated here > ! gl - 20/07/2002 > C Get ht integral of Ped, Hall conductance > C write (6,"('LAMDAS: ht integration of sigp and sigh for > C | j = ',i3)") j > do i=1,imax+1 > sigp(i,j) = 0. > sigh(i,j) = 0. > enddo > nzk = nj+nz-1 > do k=1,kmax > nzk = nzk+1 > do i=1,imax+1 > C Times ht integ by 1.e-2 because appears off for siemens (1.e-2cm/m?) > sigp(i,j) = sigp(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * sigma1(i,j,k) * 1.e-2 > sigh(i,j) = sigh(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * sigma2(i,j,k) * 1.e-2 > enddo > enddo > > C calculate the Joule heating due to wind and ion drift > C qwind = Joule Heting due to wind > C qamie = Joule Heating due to ion drift > C work = the mechanical work > C wtot = total electromagnetic power > do i=1,len1 > C sndip = sin(dipmag) > sndip(i) = sin(dipmag(i,j)) > ! if (abs(sndip(i)) .lt. 1.e-3) > ! | print *, 'sndip is zero at i= ',i,sndip(i) > C term1=(bx**2 + bz**2) > tm1(i) = sn2dec(i,j)+(1.-sn2dec(i,j))*sndip(i)**2 > C term2=(by**2 + bz**2) > tm2(i) = (1.-sn2dec(i,j))+sn2dec(i,j)*sndip(i)**2 > C term3=bxby > tm3(i) = sndec(i,j)*csdec(i,j) *(1.-sndip(i)**2) > enddo > > C write (6,"('LAMDAS: ht integration of qwind and qamie')") > > qwind(:,j) = 0. > qamie(:,j) = 0. > wtot(:,j) = 0. > work(:,j) = 0. > fwindu(:,j) = 0. > fwindv(:,j) = 0. > famieu(:,j) = 0. > famiev(:,j) = 0. > > nuik=ndj+nui-1 > nvik=ndj+nvi-1 > nuk=nj+nu-1 > nvk=nj+nv-1 > nzk = nj+nz-1 > do k=1,kmax > nuik = nuik + 1 > nvik = nvik + 1 > nuk = nuk + 1 > nvk = nvk + 1 > nzk = nzk+1 > do i=1,imax+1 > if (abs(sndip(i)) .gt. 0.2) then > C write (6,"('LAMDAS: integration of qwind and qamie at > C | i= ',i3)") i > C qwind is the Joule heating associated with nuetral wind terms > C qamie=sigp{(Bx^2+Bz^2)*Wx^2 + (By^2+Bz^2)*Wy^2-2.*Bx*By*Wx*Wy > C -2.*B^2*Vx*Wx - 2.*B^2*Vy*Wy} > qwind0 = sigma1(i,j,k)*bmod(i,j)**2* > + (tm2(i)*f(i,nuk)**2+tm1(i)* > + f(i,nvk)**2-2.*tm3(i)*f(i,nuk)*f(i,nvk) > + -2.*(f(i,nuk)*f(i,nuik)+f(i,nvk)*f(i,nvik))) > C qamie is the Joule heating withou nuetral wind > C qamie=sigp{(Bx^2+Bz^2)*Vx^2 + (By^2+Bz^2)*Vy^2+2.*Bx*By*Vx*Vy}/sin(dip)^2 > qamie0 = sigma1(i,j,k)*(tm1(i)*f(i,nuik)**2+tm2(i)* > + f(i,nvik)**2+2.*tm3(i)*f(i,nuik)*f(i,nvik)) > qamie0 = qamie0*bmod(i,j)**2 / sndip(i)**2 > C wtot = total electric power > C wtot = sigp (E + UxB).E = Qamie + sigp*{(UxB).E} > wtot0 = - sigma1(i,j,k)*bmod(i,j)**2*(f(i,nuk)*f(i,nuik) > + +f(i,nvk)*f(i,nvik)) + sigma2(i,j,k)*bmod(i,j)**2* > + (tm3(i)*(f(i,nuk)* > + f(i,nuik)-f(i,nvk)*f(i,nvik))-tm1(i)*f(i,nvk)*f(i,nuik) > + +tm2(i)*f(i,nuk)*f(i,nvik))/sndip(i) > wtot0 = qamie0 + wtot0 > C Calculate the mechanical work done by wind > work0 = wtot0 - qwind0 - qamie0 > > qwind_sec(i,j,k) = qwind0 * 1.e-9 > qamie_sec(i,j,k) = qamie0 * 1.e-9 > work_sec(i,j,k) = work0 * 1.e-9 > wtot_sec(i,j,k) = wtot0 * 1.e-9 > > > C To convert to mW/m^2 by mutipling by e-11 (since sigp (mho/m),height(cm) > C B(in Gauss), and velocity (cm/s) > qwind(i,j) = qwind(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * qwind0 * 1.e-11 > qamie(i,j) = qamie(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * qamie0 * 1.e-11 > wtot(i,j) = wtot(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * wtot0 * 1.e-11 > work(i,j) = work(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * work0 * 1.e-11 > > C Calculate the horizontal currents > fwindu0 = -tm1(i)*f(i,nvk)*sigma1(i,j,k) > + + tm3(i)*f(i,nuk)*sigma1(i,j,k) > + + sndip(i)*f(i,nuk)*sigma2(i,j,k) > fwindv0 = tm2(i)*f(i,nuk)*sigma1(i,j,k) > + + sndip(i)*f(i,nvk)*sigma2(i,j,k) > + - tm3(i)*f(i,nvk)*sigma1(i,j,k) > C Convert wind velocity from cm/s to m/s by multipling 1.e-2 and > C Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m > fwindu0 = fwindu0*bmod(i,j)*1.e-8/sndip(i) > fwindv0 = fwindv0*bmod(i,j)*1.e-8/sndip(i) > famieu0 = -tm1(i)*sigma2(i,j,k)* > + f(i,nuik) + sndip(i)*sigma1(i,j,k)*f(i,nvik) > + - tm3(i)*sigma2(i,j,k)*f(i,nvik) > famiev0 = -tm2(i)*sigma2(i,j,k)* > + f(i,nvik) - sndip(i)*sigma1(i,j,k)*f(i,nuik) > + - tm3(i)*sigma2(i,j,k)*f(i,nuik) > C Convert wind velocity from cm/s to m/s by multipling 1.e-2 and > C Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m > famieu0 = famieu0 * bmod(i,j)*1.e-8 / sndip(i)**2 > famiev0 = famiev0 * bmod(i,j)*1.e-8 / sndip(i)**2 > fwindu(i,j) = fwindu(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * fwindu0 > fwindv(i,j) = fwindv(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * fwindv0 > famieu(i,j) = famieu(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * famieu0 > famiev(i,j) = famiev(i,j) + > + (f(i,nzk+1)-f(i,nzk)) * famiev0 > > fwindu_sec(i,j,k) = fwindu0 > fwindv_sec(i,j,k) = fwindv0 > famieu_sec(i,j,k) = famieu0 > famiev_sec(i,j,k) = famiev0 > endif > enddo > enddo > > qwind_sec(:,j,zkmxp) = qwind(:,j) > qamie_sec(:,j,zkmxp) = qamie(:,j) > work_sec(:,j,zkmxp) = work(:,j) > wtot_sec(:,j,zkmxp) = wtot(:,j) > fwindu_sec(:,j,zkmxp) = fwindu(:,j) > fwindv_sec(:,j,zkmxp) = fwindv(:,j) > famieu_sec(:,j,zkmxp) = famieu(:,j) > famiev_sec(:,j,zkmxp) = famiev(:,j) > > ! print *,'LAMDAS: j,famiev(:,j) = ',j, famiev(:,j) > C Set the values at the poles > sumfaus = 0. > sumfavs = 0. > sumfwus = 0. > sumfwvs = 0. > sumfaun= 0. > sumfavn= 0. > sumfwun= 0. > sumfwvn= 0. > do i=1,imax > sumfwus = sumfwus + fwindu(i,1) > sumfwun = sumfwun + fwindu(i,36) > sumfwvs = sumfwvs + fwindv(i,1) > sumfwvn = sumfwvn + fwindv(i,36) > sumfaus = sumfaus + famieu(i,1) > sumfaun = sumfaun + famieu(i,36) > sumfavs = sumfavs + famiev(i,1) > sumfavn = sumfavn + famiev(i,36) > enddo > do i=1,imax+1 > fwindu(i,0) = sumfwus/72. > fwindu(i,37) = sumfwun/72. > fwindv(i,0) = sumfwvs/72. > fwindv(i,37) = sumfwvn/72. > famieu(i,0) = sumfaus/72. > famieu(i,37) = sumfaun/72. > famiev(i,0) = sumfavs/72. > famiev(i,37) = sumfavn/72. > enddo > ! end of calculation of additional fields > ! ======================================================================== Diff of /home/tgcm/tgcm15/mpi_mod.F and mpi_mod.F: 26c26,27 < integer :: msg_dynv(7)=(/1,2,3,4,5,6,7/) --- > integer :: msg_dynv(18)=(/1,2,3,4,5,6,7,8,9,10,11,12, > | 13,14,15,16,17,18/) 213a215 > #include "extrafields.h" 221c223 < real :: ftmp(zimxp1,zkmx),ftmpkp(zimxp1,zkmxp) --- > real :: ftmp(zimxp1,zkmx),ftmpkp(zimxp1,zkmxp),ftmp1(zimxp1) 222a225 > 273a277,364 > ! famiev_sec (amie meridinal current component with height) > ftmpkp(:,:) = famiev_sec(:,lat,:) > call mpi_send(ftmpkp,zimxp1*zkmxp,MPI_REAL8,idest,msg_dynv(8), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for famiev_sec: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! fwindv_sec (wind meridinal current component with height) > ftmpkp(:,:) = fwindv_sec(:,lat,:) > call mpi_send(ftmpkp,zimxp1*zkmxp,MPI_REAL8,idest,msg_dynv(9), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for fwindv_sec: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! famiev (amie meridional current component) > ftmp1(:) = famiev(:,lat) > call mpi_send(ftmp1,zimxp1,MPI_REAL8,idest,msg_dynv(10), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for famiev: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! fwindv (wind meridinal current component) > ftmp1(:) = fwindv(:,lat) > call mpi_send(ftmp1,zimxp1,MPI_REAL8,idest,msg_dynv(11), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for fwindv: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! qamie_sec (amie Joule heating) > ftmpkp(:,:) = qamie_sec(:,lat,:) > call mpi_send(ftmpkp,zimxp1*zkmxp,MPI_REAL8,idest,msg_dynv(12), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for qamie_sec: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! qwind_sec (wind Joule heating) > ftmpkp(:,:) = qwind_sec(:,lat,:) > call mpi_send(ftmpkp,zimxp1*zkmxp,MPI_REAL8,idest,msg_dynv(13), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for qwind_sec: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! work_sec (mechanical work) > ftmpkp(:,:) = work_sec(:,lat,:) > call mpi_send(ftmpkp,zimxp1*zkmxp,MPI_REAL8,idest,msg_dynv(14), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for work_sec: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! wtot_sec (total pmechanical power) > ftmpkp(:,:) = wtot_sec(:,lat,:) > call mpi_send(ftmpkp,zimxp1*zkmxp,MPI_REAL8,idest,msg_dynv(15), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for wtot_sec: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! tec_sec (TEC) > ftmpkp(:,:) = tec_sec(:,lat,:) > call mpi_send(ftmpkp,zimxp1*zkmxp,MPI_REAL8,idest,msg_dynv(16), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for tec_sec: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! sigp (height-integrated Pedersen conductivity) > ftmp1(:) = sigp(:,lat) > call mpi_send(ftmp1,zimxp1,MPI_REAL8,idest,msg_dynv(17), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for sigp: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! > ! sigh (height-integrated Hall conductivity) > ftmp1(:) = sigh(:,lat) > call mpi_send(ftmp1,zimxp1,MPI_REAL8,idest,msg_dynv(18), > | MPI_COMM_WORLD,ier) > if (ier /= 0) write(6,"('>>> mp_snddyn: error from mpi_send', > | ' for sigh: ier=',i3,' lat=',i3,' idest=',i2)") > | ier,lat,idest > ! 285a377 > #include "extrafields.h" 292c384 < real :: ftmp(zimxp1,zkmx),ftmpkp(zimxp1,zkmxp) --- > real :: ftmp(zimxp1,zkmx),ftmpkp(zimxp1,zkmxp),ftmp1(zimxp1) 346a439,518 > ! famiev_sec (amie meridinal current component with height) > call mpi_recv(ftmpkp,zimxp1*zkmxp,MPI_REAL8,isrc,msg_dynv(8), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for famiev_sec: ier=',i3,' lat=',i3)") ier,lat > famiev_sec(:,lat,:) = ftmpkp(:,:) > ! > ! fwindv_sec (wind meridinal current component with height) > call mpi_recv(ftmpkp,zimxp1*zkmxp,MPI_REAL8,isrc,msg_dynv(9), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for fwindv_sec: ier=',i3,' lat=',i3)") ier,lat > fwindv_sec(:,lat,:) = ftmpkp(:,:) > ! > ! famiev (amie meridinal current component) > call mpi_recv(ftmp1,zimxp1,MPI_REAL8,isrc,msg_dynv(10), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for famiev: ier=',i3,' lat=',i3)") ier,lat > famiev(:,lat) = ftmp1(:) > > ! write(6,"('mp_rcvdyn: lat=',i3,' famiev=',/,(6e12.4))") > ! | lat,famiev(:,lat) > ! > ! fwindv (wind meridinal current component) > call mpi_recv(ftmp1,zimxp1,MPI_REAL8,isrc,msg_dynv(11), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for fwindv: ier=',i3,' lat=',i3)") ier,lat > fwindv(:,lat) = ftmp1(:) > ! > ! qamie_sec (amie Joule heating) > call mpi_recv(ftmpkp,zimxp1*zkmxp,MPI_REAL8,isrc,msg_dynv(12), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for qamie_sec: ier=',i3,' lat=',i3)") ier,lat > qamie_sec(:,lat,:) = ftmpkp(:,:) > ! > ! qwind_sec (wind Joule heating) > call mpi_recv(ftmpkp,zimxp1*zkmxp,MPI_REAL8,isrc,msg_dynv(13), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for qwind_sec: ier=',i3,' lat=',i3)") ier,lat > qwind_sec(:,lat,:) = ftmpkp(:,:) > ! > ! work_sec (mechanical work) > call mpi_recv(ftmpkp,zimxp1*zkmxp,MPI_REAL8,isrc,msg_dynv(14), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for work_sec: ier=',i3,' lat=',i3)") ier,lat > work_sec(:,lat,:) = ftmpkp(:,:) > ! > ! wtot_sec (total mechanical work) > call mpi_recv(ftmpkp,zimxp1*zkmxp,MPI_REAL8,isrc,msg_dynv(15), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for wtot_sec: ier=',i3,' lat=',i3)") ier,lat > wtot_sec(:,lat,:) = ftmpkp(:,:) > ! > ! tec_sec (TEC) > call mpi_recv(ftmpkp,zimxp1*zkmxp,MPI_REAL8,isrc,msg_dynv(16), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for tec_sec: ier=',i3,' lat=',i3)") ier,lat > tec_sec(:,lat,:) = ftmpkp(:,:) > ! > ! sigp (height-integrated Pedersen conductivity) > call mpi_recv(ftmp1,zimxp1,MPI_REAL8,isrc,msg_dynv(17), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for sigp: ier=',i3,' lat=',i3)") ier,lat > sigp(:,lat) = ftmp1(:) > ! > ! sigh (height-integrated Hall conductivity) > call mpi_recv(ftmp1,zimxp1,MPI_REAL8,isrc,msg_dynv(18), > | MPI_COMM_WORLD,irstat,ier) > if (ier /= 0) write(6,"('>>> mp_rcvdyn: error from mpi_recv', > | ' for sigh: ier=',i3,' lat=',i3)") ier,lat > sigh(:,lat) = ftmp1(:) > ! ======================================================================== Diff of /home/tgcm/tgcm15/potm.F and potm.F: 3a4,5 > use input_module,only: iamie > use amie_module,only: tiepot 28a31,44 > real,dimension(imaxmp,-2:zkmxp) :: phihm_sech > > IF (IAMIE .EQ. 1) THEN > ! write(6,"('POTM: calling AMIE potential for iamie=',i3)") > ! | iamie > do J=1,JMAXM > do I=1,IMAXM > PHIHM(I,J) = TIEPOT(I,J) > enddo > enddo > > GO TO 1 > ENDIF > 32c48 < DO 1 J = 1,JMAXM --- > do J = 1,JMAXM 36c52 < DO 2 I = 1,IMAXM --- > do I = 1,IMAXM 41c57 < 2 CONTINUE --- > enddo 48c64 < DO 3 I = 1,IMAXM --- > do I = 1,IMAXM 50c66 < 3 CONTINUE --- > enddo 51a68,69 > enddo > 58a77,86 > > C save into secondary history > ! do j=1,jmaxm > ! do i=1,imaxmp > ! phihm_sech(i,:) = phihm(i,j) > ! enddo > ! call addfsech('PHIHM',' ',' ',phihm_sech, > ! | imaxmp,zkmxp+3,zkmxp+3,j) > ! enddo > ======================================================================== WARNING: cannot find file /home/tgcm/tgcm15/rgrd1.F ======================================================================== WARNING: cannot find file /home/tgcm/tgcm15/rgrd2.F ======================================================================== WARNING: cannot find file /home/tgcm/tgcm15/rgrd3.F ======================================================================== Diff of /home/tgcm/tgcm15/settei.F and settei.F: 388,390c388,391 < S10(I,1) = 1.57E-12*exp((2.4E4+0.3*(F(I,NTEK)-1500.)-1.947E-5 < 1 *(F(I,NTEK)-1500.)*(F(I,NTEK)-4000.))*(F(I,NTEK)-3000.)/(3000. < 2 *F(I,NTEK)))*S10(I,1) --- > ! S10(I,1) = 1.57E-12*exp((2.4E4+0.3*(F(I,NTEK)-1500.)-1.947E-5 > ! 1 *(F(I,NTEK)-1500.)*(F(I,NTEK)-4000.))*(F(I,NTEK)-3000.)/(3000. > ! 2 *F(I,NTEK)))*S10(I,1) > S10(I,1) = 0. ======================================================================== Diff of /home/tgcm/tgcm15/tail.F and tail.F: 4c4,14 < use input_module,only: power,ctpoten,byimf,f107,f107a --- > use input_module,only: power,ctpoten,byimf,f107,f107a,iamie > use amie_module,only: hpi_nh_amie,hpi_sh_amie, > | pcp_nh_amie,pcp_sh_amie,crad > > ! am_03/02 the parameterization was changed to fit the parameterization > ! from B.Emery > ! before the integration of the flux didn't give the right > ! hpower index (it was too high eg. POWER=3 resulted into > ! an index of 13) > ! Ray Roble used snoe values which are fitted to the > ! observations 15d24 < #include "amie.h" 32c41 < | plevel,hp,cp,arad(2),alfa21,alfa22 --- > | plevel,hp,cp,arad(2),alfa21,alfa22,c25,c35 42a52,77 > > crit(1) = 0.261799387 > crit(2) = 0.523598775 > > ! Use AMIE data > if (iamie == 1) then > > theta0(1) = crad(1)*180./pi > theta0(2) = crad(2)*180./pi > hp = amax1(hpi_sh_amie,hpi_nh_amie) > cp = amax1(pcp_sh_amie,pcp_nh_amie) > PLEVEL = 2.09 * ALOG(hp) > H1 = AMIN1 (2.35, 0.83 + 0.33 * PLEVEL ) > H2 = 2.87 + 0.15 * PLEVEL > ROTH = (12.18 - 0.89*PLEVEL) * 15. > ROTE = (2.62 - 0.55*PLEVEL) * 15. > OFFC(1) = 1. > OFFC(2) = 1. > DSKOFC(1) = 0. > DSKOFC(2) = 0. > ! RDEG(5) = OFFC(1)*DTR > ! RDEG(6) = OFFC(2)*DTR > ! RDEG(9) = DSKOFC(1)*DTR > ! RDEG(10) = DSKOFC(2)*DTR > > endif 48,51c83,86 < ! ec = 0.01+0.09*power/100. ! kibo < ! alfac = 0.5 ! kibo < ec = 0.5 ! snoe < alfac = 1.0 ! snoe --- > ec = 0.01+0.09*power/100. ! kibo > alfac = 0.5 ! kibo > ! ec = 0.5 ! snoe > ! alfac = 1.0 ! snoe 53,56c88,91 < ! ed = 0.01+0.2*power/100. ! kibo < ! alfad = 0.75 ! kibo < ed = 0.5 ! snoe < alfad = 2.0 ! snoe --- > ed = 0.01+0.2*power/100. ! kibo > alfad = 0.75 ! kibo > ! ed = 0.5 ! snoe > ! alfad = 2.0 ! snoe 65a101,103 > > if (iamie == 1) goto 190 > 67,72c105 < 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) --- > theta0(i) = -3.80 + 8.48 * (ctpoten**0.1875) 125,126c158,160 < C E1 = AMAX1( 0.50, -2.15 + 0.62 * PLEVEL) < C E2 = 0.95 + 0.117 * power --- > PLEVEL = AMAX1(1.0,PLEVEL) > E1 = AMAX1( 0.50, -2.15 + 0.62 * PLEVEL) > E2 = 0.95 + 0.117 * power 133,134c167,168 < E1 = (1.0 + 0.25* power) ! snoe < E2 = (1.0 + 0.25* power) ! snoe --- > ! E1 = (1.0 + 0.25* power) ! snoe > ! E2 = (1.0 + 0.25* power) ! snoe 136,139c170,173 < ! H1 = AMIN1( 2.35, 0.83 + 0.33 * PLEVEL ) ! kibo < ! H2 = 2.87 + 0.15 * PLEVEL ! kibo < H1 = 3.+0.1*power ! snoe < H2 = 3.+0.1*power ! snoe --- > H1 = AMIN1( 2.35, 0.83 + 0.33 * PLEVEL ) ! kibo > H2 = 2.87 + 0.15 * PLEVEL ! kibo > ! H1 = 3.+0.1*power ! snoe > ! H2 = 3.+0.1*power ! snoe 153,154c187,188 < TWA6 = 0. < TWA21 = 0. --- > TWA6 = 0.36 + 0.48 * PLEVEL > TWA21 = AMAX1( 1.00, -1.75 + 0.69 * PLEVEL ) 184,185c218,219 < ! ALFA2 = 3. ! kibo < ALFA2 = 2. ! snoe --- > ALFA2 = 3. ! kibo > ! ALFA2 = 2. ! snoe 203c237,258 < IF (IAMIE .EQ. 1) RETURN --- > > if (iamie == 1) then > C Set critical radii: CRIT(1) = theta0 + 5.0 deg > C for AMIE CRIT(2) = theta0 + 10.0 deg > C CRIT(1) = ((THETA0(1)+THETA0(2))*0.5 + 5.0) * DTR > C CRIT(2) = ((THETA0(1)+THETA0(2))*0.5 + 10.0) * DTR > C Set critical colatitude crit(2) 40 deg -- G. Lu 5/11/98 > C Therefore, crit2 = 40, crit1 = 25-30, depending on theta0 > C35 = AMIN1(30.0,(THETA0(1)+THETA0(2))*0.5 + 5.0) > C25 = AMAX1(25.,C35) > CRIT(1) = C25 * DTR > CRIT(2) = 40. * DTR > theta0(:) = theta0(:) * dtr > offa(:) = offa(:) * dtr > offc(:) = offc(:) * dtr > dskofa(:) = dskofa(:) * dtr > dskofc(:) = dskofc(:) * dtr > > return > > endif >