PROGRAM THRDM C C ALAN BURNS C JULY 1987 C Ion trajectories included by Alan Burns and Dick Cannata, May 1988 C Minor species calculations included by Alan Burns, July 1988 C Updated to run using UNIX August 1990 Alan Burns. C Modified for the TIEGCM by Alan Burns Oct 17 1994. Note C that changes have not yet been made to make the ion drag C calculations self-consistent. C Modified in 2002 to enable it to run with Wenbin's histories on a work C station. Alan Burns Jan 31 2003 C Checking vertical winds C Currently the TIEGCM histories are being obtained by linking to Ben Foster's C getgcm library. At the moment he outputs height in cm (but with the 97km C added) and neutral winds in cm/s. The diagnostic code has been C modified to account for this. However I have no control over his C modifying his libraries, so the user should be ready to modify this C code in response to Ben's potential modifications. Therefore, neutral C winds and the altitudes of the pressure surfaces should always be C included in any files brought back or any output graphics made. C C C THIS PROGRAM IS MODIFIED FROM THE TGCM ISP PROCESSOR C AND THE DIAGNOSTIC PROCESSOR TO ALLOW TERMS TO BE CALCULATED C ALONG THE PATH OF A THREE DIMENSIONAL TRAJECTORY. THESE C THREE DIMENSIONAL TRAJECTORIES ARE TO BE PLOTTED IN A C 3 D PLOT SET. This package is divided into 6 parts which are C run by typing EXEC next to THRDR2 EXEC. The other component parts C are THRDM JOB, THRS2 JOB,THRDS JOB, THRDJ2 JOB and TREAD DAT. C THRDM JOB is the main program that calls subroutines from THRDJ2 JOB, C THRDS, THRS2 and itself. THRDJ2 JOB is the part that contains the C job control language (JCL) for this set of programs and amalgamates C the ensemble when exec is typed next to the EXEC file. TREAD DAT C has the data for the lexical reads that supply inputs for the auroral C calculations. The following need to be set. C 1) THE LEXICAL READS for THE DATA VOLUME BEING USED C C 2) IBCK TO 1 FOR BACKWARDS TRAJECTORIES, 0 for forward C C 3) THE TRAJECTORY STARTING POINTS (LAT,LON,Z,UTSRT) C C 4) The parameters NTJ and NTM determine the number of C trajectories to be calculated and the number of hours C for which the trajectories are to be run respectively. C NTM equals NUMTMS in THRS2.SH C C 5) The data and history statements in THRS2 must be set to C the data and history files that you want to use. C C 6) REECH can be set to adjust the size of the polar C dial, but 50 degrees is the most aesthetically pleasing C C 8) Set JEOMGN=1 for ion trajectories (else 0). Set KEOMGN C =1 for neutral trajectories in geomagnetic coordinates C C 9) Set ISMIN = 1 for minor species calculations(TIGCM) C C 10) Set OLD,SUBAU,NEW,NEWH and AMIE to the values C that are appropriate for the run (DATA statement) C PARAMETER(NFHIST=17) PARAMETER(IMX=73,JMX=36,KMX=25,NF=7,IZPLT=3,IDIF=0) PARAMETER (NLAT=36,NLONG=73,NZK=25,NTHM=36,NLMB=73) PARAMETER (ITRU=73,JTRU=36,KTRU=25,MXFLDS=17) PARAMETER (MXUT=25,MXDAYS=20,MXVOLS=10,MXTMS=73) PARAMETER(NV=2,NTJ=8,NTM= 12 ) C C LXDICT NEEDS TO BE SET WHEN YOU HAVE TOO MANY LEXICAL READS COMMON /LXDICT/KDICT,NDX,MXD,NID,IDICT(3,64) COMMON /DATAVOL/ PQN(73,25,36),UVIDYN(73,36,2),DISO2(IMX,KMX,JMX), 1 DVNE(IMX,KMX,JMX),WIDYN(73,36) CHARACTER*8 DATAVOM(3),HISTVOL(3,MXVOLS) COMMON/CHINP/DATAVOM,HISTVOL COMMON/INPUT1/MTGCM,MD(MXTMS),MH(MXTMS),MM(MXTMS),MDUT(MXTMS), + MHUT(MXTMS),MMUT(MXTMS),UTG(MXTMS),UTINC(MXTMS),IXGCM(MXFLDS), + LUHIST(MXVOLS),NHVOLS,NFLDS,NTMS,INDEN,IPRNT COMMON /CHAR20/ADATA(3),SOURCE(3),OUTPUT(3,5),ALABEL(40) C *********** New inputs for TGCMRD 4/91 to correct for C *********** for errors that occurred with the new compiler COMMON /RDINPT/IV1,IFRST,UTP,MDUTP,NCALLS,IEOF(MXVOLS) ,MDP,MHP 1 ,MMP DATA IV1/1/,IFRST/0/,UTP/-1./,MDUTP/-1./,NCALLS/0/ DATA IEOF/MXVOLS*0/ C ********** End of changes C COMMON /PSUTS/ UT1(MXDAYS),UT1PL(MXDAYS) C COMMON /BLK20/ HNT(33), REL(17), GLON(73), ZP(25), GLAT(36), 1 IACT,JACT,KACT,INDX(73),JNDX(36),KNDX(25),TLON(ITRU),TLAT(JTRU), 2 ZSRF(25), PNT(ITRU,KTRU,JTRU,NFHIST), T0(27), B22(2,2) , 3 Q0(25),JSWSRB(1),JSWEUV,JSW2FL ,IYRDAY,F107D,F107AV,ZUR, 4 IGMSIS,IGIRIR DIMENSION IND1(73),JND1(36),KND1(25) DATA IND1/73*1/, JND1/36*1/, KND1/25*1/ COMMON /IONPR/ DAYY, IYD, F107, F107A, AP(7), STL, UTSEC, 1 GEOLAT, GEOLONG, GMLATI, GMLONGI COMMON /TRGM/ RLATM(76,36),RLONM(76,36),DIPMAG(76,36), 1 DECMAG(76,36),SNDEC(76,36),CSDEC(76,36),SN2DEC(76,36), 2 SNCSDC(76,36),RLATMP(36),RLONMP(36),DLONS(36) COMMON /JSLICE/ 3 DTOR,RTOD,PI,XNO(KMX),XNN2(KMX),XNO2(KMX),TN(KMX),UN(KMX),VN(KMX) 4 ,WN(KMX),RHO(KMX),SHT(KMX),AMAS(KMX), POTT(IMX,KMX), 5 VIU(IMX,KMX),VIE(IMX,KMX),RE,RZUR, 6 TMO, XNE(IMX,KMX),BOLTZ,GZ(KMX),ALAMXX(KMX),ALAMYY(KMX), 7 ALAMXY(KMX),ALAMYX(KMX),BGAUS,UI(IMX,JMX),VI(IMX,JMX),WI(IMX,JMX) 8 ,UIZ(KMX),VIZ(KMX), 9 GNLAT(IMX,JMX),UM(KMX),VM(KMX),SIGPED(IMX,KMX), 1 SIGHAL(IMX,KMX),TJOULE(IMX,JMX),TAUROR(IMX,JMX),TJ0(IMX,JMX) C C BUFFER ARRAY FOR ABSTRACT VOLUME FOR TERM ANALYSIS PROGRAM COMMON/HISFLD/BUFDAT(6,4,6,37),TOO(27) COMMON/HISTTO/DUFDAT(6,4,6,91) C Extra info on O+ DIMENSION EUFDAT(IMX,KMX,JMX,4) DIMENSION CUFDAT(6,KMX,6), Z2PKT(NTM,NTJ,25), QI3PKT(NTM,NTJ,25) DIMENSION Z3PKT(NTM,NTJ) DIMENSION QPPHT(NTM,NTJ,25) C Array for the total ionization DIMENSION QTI(3) C Dimension the auroral efficiency DIMENSION AUREFF(KMX) C Include the auroral efficiency to calculate the auroral heating terms. DATA AUREFF 1/0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55, 2 0.50,0.45,0.37,0.30,0.20,0.16,0.13,0.10,0.08,0.07,0.06,0.05/ C Dimension the arrays for the TGCM inputs DIMENSION QINH(IMX,JMX),RO2(IMX,JMX) C C IKPARMS CONTOUR LONGITUDE VS PRESSURE AT SELECTED LATITUDES COMMON /IKPARR/ 1 PRHO(IMX,KMX),PSHT(IMX,KMX),PXNO(IMX,KMX),PXNN2(IMX,KMX), 2 PXNO2(IMX,KMX), 6 PLXX(IMX,KMX),PLXY(IMX,KMX),PLYY(IMX,KMX),PLYX(IMX,KMX) C Common block to output ALAM1 from inputs COMMON /ALAHOL/ALAM1HLD(IMX,KMX) DIMENSION ALAMPKT(NTM,KMX,NTJ) DIMENSION ALAMDAT(6,4,6,KMX) C COMMON/CSW/SW(25), ISW, SWC(25) DIMENSION SV(25) DATA SV/25*1/ COMMON/XYTRANS/XX(500), YY(500), IBX, IBY, IREV, IDEF C COMMON/LEVELO/ILEVEL,IAUR COMMON/IONSTR/XIOP(KMX),XINP(KMX),XINOP(KMX),XIO2P(KMX), 1XIN2P(KMX),XNEE(KMX),QNA(KMX) C C COMMON/GLHVI/ AF(73,36), VELEW(73,36), VELNS(73,36), FLUX(73,36), C 1 ALFA(73,36), POTEN(73,36), CUSP(73,36), DRIZL(73,36), RLON(72), C 2 RLAT(72), DPHI, GVIMX C**************************************************************** C GLHVI has been replaced to make it compatitble with flowx12 rather C than flowx8. 6 September 1991 - Alan Burns COMMON/GLHVI/ AF(73,36), VELEW(73,36), VELNS(73,36), FLUX(73,36), 1 ALFA(73,36), POTEN(73,36), CUSP(73,36), DRIZL(73,36), RLON(72), 2 RLAT(72), DPHI, GVIMX, AF2(73,36),ALFA2(73,36),FLUX2(73,36) REAL MAG(4) COMMON/CONS/IMAX,JMAX,WWW(120),INPT(30),DIPOLE(2),PMLONG(2),IYEAR, 1IDAY,SEC C COMMON/OVALR/RROT,RDISP,H0,RH,FC,ALFAC,FD,ALFAD,E0,REE,ALFK, C 1 ALF1,ALF2,RROT1,RAROT2,RD1A,RD1V,RD2,RAH1,RAH2,RAHT1, C 2 RAHT2,FALFC,FALF6,FALF22,ALFA0,RALFA C**************************************************************** C OVALR has been replaced to make it compatitble with flowx12 rather C than flowx8. 6 September 1991 - Alan Burns C Move data to and from the subroutine DIAGPKG3 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 COMMON/TIMIT/UTOO,TIMSTP,JUT COMMON/LIMTOO/IZERO,JZERO,KZERO,IFLAG,LFLAG,IMTGCM COMMON/DATFLD/IZ3 COMMON/FLAGIN/KFLAG(6,4,6),ISMIN,KION COMMON /TRAJ/ UPKT(NTM,NTJ),VPKT(NTM,NTJ),WPKT(NTM,NTJ), 1 UXPKT(NTM,6,NTJ),VXPKT(NTM,6,NTJ),QPKT(NTM,13,NTJ), 2 TPKT(NTM,NTJ),COMPKT(NTM,6,NTJ),UTPKT(NTM),FPKT(NTM,3,NTJ), 3 XOPKT(NTM,NTJ),XO2PKT(NTM,NTJ),XNPKT(NTM,NTJ),XLT(NTM,NTJ), 4 AZPKT(NTM,NTJ),UTSRT(20),AKSTRT(20), 5 FORNOR(NTM,NTJ),FORPAR(NTM,NTJ),PKTMAG(NTM,NTJ),FORMAG(NTM,NTJ), 6 QPKT2(NTM,NTJ),UIPKT(NTM,NTJ),VIPKT(NTM,NTJ), 7 TIPKT(NTM,NTJ),TEPKT(NTM,NTJ) COMMON/TRAJMISC/COMPKT2(NTM,7,NTJ) COMMON /POSIT/ ALON(NTM,NTJ),ALAT(NTM,NTJ),ALATS(20),ALONS(20) COMMON /WGHTS/W1(NTJ),W2(NTJ),W3(NTJ),W4(NTJ),W5(NTJ), 1 W6(NTJ),W7(NTJ),W8(NTJ),KK C Output arrays for the minor species COMMON /TRAJ2/AN4PKT(NTM,16,NTJ),ANOPKT(NTM,10,NTJ), 1 AN2DPKT(NTM,11,NTJ), 2 BN4PKT(NTM,6,NTJ),BNOPKT(NTM,6,NTJ), 3 BN2DPKT(NTM,3,NTJ), 4 PSN4PK(NTM,NTJ),PSNOPK(NTM,NTJ), 5 PSN2DPK(NTM,NTJ),AOPPKT(NTM,10,NTJ), 6 ANPPKT(NTM,7,NTJ),AN2PPKT(NTM,6,NTJ), 7 AO2PPKT(NTM,9,NTJ),ANOPPKT(NTM,11,NTJ), 8 DOPPKT(NTM,NTJ),DNPPKT(NTM,NTJ),DN2PPKT(NTM,NTJ), 9 DO2PPKT(NTM,NTJ),DNOPPKT(NTM,NTJ) DIMENSION KUHLP2(NTJ),U(IMX,KMX,JMX),V(IMX,KMX,JMX), 1 W(IMX,KMX,JMX),POO2(IMX,KMX,JMX),POO(IMX,KMX,JMX),GCNLON(73), , 2 GCNLAT(JMX),IPKT(NTJ),JPKT(NTJ),KPKT(NTJ),IPKU(NTJ),JPKU(NTJ), 3 KPKU(NTJ),ZKP1(NTJ),ZN(NTM),XN(NTM), 4 YN(24),AGLAT(NTJ),AGLON(NTJ) C Inputs for CHAPMN COMMON/DENHT/OXNO2(KMX),OXNO(KMX),OXNN2(KMX),ZPHT(KMX),TNNX(KMX), 1 XNEG(KMX) COMMON/AHTES/ZPP(KMX),PSIXX(KMX,3),AHMBAR(KMX),ZXX(KMX) COMMON/SOPTS/III,JJJ,IIDAY,UTX C This is used in RATECOF and elsewhere COMMON/PSIS/PSI(3),F107G,XNA(3) C Input to SOLHEAT COMMON/FTEN/F107AVE C Output from SOLHEAT COMMON/SOLOUT/QHEAT,QION,QCONT,QBAND,RJ C Input for IONMOD COMMON/XTRAN/BTN,BTI,BTE,XNN2D,XNN4S,XNNO C Output from IONOUT COMMON/IONOUT/QIC,RJIC C Output from NOHEAT COMMON/OUTNO/QNMIN,RJNO C Output from COLHEAT COMMON/COLOUT/QENEI,QEN,QEI C Output from IRCOOL COMMON/OUTIR/CO2COOL,ANOCOOL,O3PCOOL,QRAD C Minor species COMMON/MININP/AEMBAR(3),DAMBAR(3),APSIO(3),APSIO2(3),DPSIO(3), 1 DPSIO2(3),TEMIT(3) COMMON/MINMIX/PSMIN(3,3) C N(4S) sources COMMON/SRCN4S/SN4S4,SN4S5,SN4S7,SN4S8,SN4SK2,SN4SK6,SN4SK8, 1 SN4SA1,SN4SA3,SN4SE C N(4S) loss terms COMMON/LSSN4S/ALN4S1,ALN4S3,ALN4K4 C NO sources COMMON/SRCNOS/SNOS1,SNOS2 C NO loss terms COMMON/LSSNOS/ALNOS3,ALNOS6,ALNOS8,ALNOS9,ALNOK5 C N(2D) sources COMMON/SRCN2D/SN2DK3,SN2DA1,SN2DA3,SN2DE C N(2D) loss terms COMMON/LSSN2D/ALN2D2,ALN2D4,ALN2D5,ALN2D6,ALN2D7,ALN2DK10 C New N(2D) COMMON/PASSN2D/PN2D C Electromagnetic radiative production of ions COMMON/QSIONS/QSINP,QSOP,QSO2P,QSN2P,QSNOP C Ion production by auroral precipitation COMMON/QAIONS/QAINP,QAOP,QAO2P,QAN2P,QANOP C Ion production and loss COMMON/IONPRO/SNP1,ALNP1,ALNP2,ALNP3,ALN2P1,ALN2P2,ALN2P3, 1 SO2P1,SO2P2,SO2P3,ALO2P1,ALO2P2,ALO2P3, 2 SNOP1,SNOP2,SNOP3,SNOP4,SNOP5 3 ,SNOP6,SNOP7,ALNOP1,JION C Ion densities COMMON /DENION/XINPJ,XIN2PJ,XIOPJ,XIO2PJ,XINOPJ C O+ input COMMON/OPIN/ZFAC(3),UNEU(3,3,3),VNEU(3,3,3),UII,VII,WII, 1 CXNOP(3,3),FLAT(3),WNEU(3,3,3) C O+ terms COMMON/OUTOP/DIFOP,SOP1,ALOP1,ALOP2,ALOP3,PHIAV,OPV,ADVOPN,ADVOPI DIMENSION XIOJP(3),CTN(3),CTE(3),CTI(3) C Common block from MAGFIEL giving magnetic field components COMMON /MAGFLD/BX(3,3,3),BY(3,3,3),BZ(3,3,3),BMOD(3,3,3) C Common block to permit solar heating terms to be calculated COMMON /SOLH/IFACT1,IFACT2 C Pass information about whether the program is running on the TIEGCM COMMON /TIEPASS/NTIE C PARAMETER (MXDAY=10) COMMON/OUTOPEX/UPADVOP,UPNEU DIMENSION EYE(3),ZPKT(NTM,NTJ),BLAT(NTM) DIMENSION CHMAG1(NTJ),CHMAG2(NTJ) C Initialize the common block for REECH the variable that defines the C latitudinal extent of the polar plots. COMMON/PLLAT/REECH COMMON/MGBIT/JEOMGN,IHTSW,QHEX(6,2,6),QFRIC(6,2,6),TION(6,2,6), 1 WNSHR(6,2,6),Q2FRIC(6,25,6),QDENS(6,25,6) DIMENSION QIPKT(NTM,3,NTJ),QI2PKT(NTM,3,NTJ), 1ENDENS(NTM,NTJ), VIVN(NTM,3,NTJ) C Character declaration for plotting CHARACTER*40 ICHARA C New UTPKT to allow for plotting over day boundaries DIMENSION UTPKT2(NTM) C Inputs for equatorial trajectory plots DIMENSION ALATLAB(5),ALONLAB(5),ALSTLAB(5) DIMENSION ALONA(NTM,NTJ),ALATA(NTM,NTJ) DATA ALATLAB/-60.,-30.,0.,30.,60./ DATA ALONLAB/-180.,-90.,0.,90.,180./ DATA ALSTLAB/0.,6.,12.,18.,24./ DATA AMINLON,AMAXLON,AMINLAT,AMAXLAT/-180.,180.,-60.,60./ DATA AMINLST,AMAXLST/0.,24./ C COORDS OF S (LAT,LON) AND N OR CENTERED (LAT,LON) DIPOLE DATA MAG/-74.5,127.,79.,-70./ C COORDS FOR EQUAL GEOGRAPHIC AND GEOMAGNETIC POLES C DATA MAG/-89.9,127.,89.9,-70./ DATA SPVAL/-9999./ DATA IBX,IBY,IREV,IDEF/1,1,0,0/ C ******************************************************************** C NEWHA and NEWM are additional elements to the call to TAIL that have C been added (6 September 1991) to make this call compatible with C FLOWX12 - Alan Burns LOGICAL OLD,NEW,NEWH,NEWHA,NEWM,SUBAU,AMIE DATA AMIE/.FALSE./ DATA NEWHA,NEWM/.FALSE.,.FALSE./ DATA OLD,NEW,NEWH,SUBAU/.FALSE.,.TRUE.,.TRUE.,.FALSE./ DATA TIMSTP/3600./ DATA ASTP/1./ C SET IBCK=1 FOR BACK FILES DATA IBCK/1/ IFUNFI=1 C Set the latitude range for polar plots DATA REECH/40./ C this calculates ion trajectories DATA JEOMGN/0/ C This is set to 1 when you wish to calculate neutral trajectories C in geomagnetic coordinates. DATA KEOMGN/1/ C If ISMIN = 1 then include the minor species calculations. DATA ISMIN/1/ C Work out the ion production terms DATA JION/1/ C Switch to make height calculations ie ped and hall when needed C ihtsw=1 means all heights are used DATA IHTSW/0/ C Switch for the new TIGCM arrangement so that new model hour = UT hour C New when NEWTIME = 1 DATA NEWTIME/1/ C C Trajectory starting values C Pressure height surface c DATA AKSTRT/15.0,15.0,15.0,9.0,9.0,9.0,11.0,11.0,11.0,11*0./ C DATA AKSTRT/7.0,7.0,7.0,9.0,9.0,9.0,11.0,11.0,11.0,11*0./ C Sharp study 2 DATA AKSTRT/5.0,7.0,9.0,11.0,5.0,7.0,9.0,11.0,11.0,11*0./ C Sharp study C DATA AKSTRT/4.0,8.0,9.0,9.0,9.0,9.0,11.0,11.0,11.0,11*0./ C IUGG95 c DATA AKSTRT/19.0,8.0,9.0,9.0,9.0,9.0,11.0,11.0,11.0,11*0./ c DATA AKSTRT/19.0,19.0,19.0,19.0,19.0,19.0,19.0,19.0,19.0,19.0, c 1 17.0,15.0,15.0,15.0,15.0,15.0,15.0,15.0,15.0,15.0/ C Geographic latitude C Thule Sonde c DATA ALATS/79.,64.,69.5,60.6,65.1,69.5,60.6,65.1,69.5,11*0./ c DATA ALATS/60.6,65.1,69.5,60.6,65.1,69.5,60.6,65.1,69.5,11*0./ C DATA ALATS/60.6,65.,70.0,60.,65.,70.,60.,65.,70.0,11*0./ C Sharp study 2 May 1996 DATA ALATS/65.1,65.1,65.1,65.1,70.1,70.1,70.1,70.1,70.0,11*0./ C Sharp study May 95 C DATA ALATS/65.1,65.1,65.1,60.6,65.1,69.5,60.6,65.1,69.5,11*0./ C IUGG95 c DATA ALATS/65,65.1,65.1,60.6,65.1,69.5,60.6,65.1,69.5,11*0./ c DATA ALATS/-85.,-62.5,-65.0,-45.,-40.,30.0,85.0,80.0,75.0,70.0, c 1 -55.5,74.0,66.0,60.0,54.0,82.0,74.0,66.0,60.0,54.0/ C Geographic longitude C DATA ALONS/-52.0,-47.5,-41.1,-52.0,-47.5,-41.1,-52.0, C 1 -47.5,-41.1,11*0./ C Sharp Study 2 May 1996 DATA ALONS/8*-147.5, 1 -41.1,11*0./ C Sharp study May 95 c DATA ALONS/-47.5,-47.5,-47.5,-52.0,-147.5,-141.1,-152.0, C 1 -147.5,-141.1,11*0./ c DATA ALONS/90.,-147.5,-147.5,-152.0,-147.5,-141.1,-152.0, c 1 -147.5,-141.1,11*0./ C Thule Sonde c DATA ALONS/-70.0,-50.,-141.1,-152.0,-147.5,-141.1,-152.0, c 1 -147.5,-141.1,11*0./ c DATA ALONS/48.0,48.0,-130.0,48.0,-130.0,120.0,160.0,160.0,160.0, c 1 160.0,155.0,-30.0,-30.0,-30.0,-30.0,180.0,180.0,180.0, c 1 180.0,180.0/ C This is not used until I find a better way to do it DATA UTSRT/14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,14.0,11*0./ C DATA UTSRT/22.0,22.0,22.0,14.0,14.0,14.0,14.0,14.0,14.0,11*0./ C IUGG 95 c DATA UTSRT/22.0,22.0,22.0,14.0,14.0,14.0,14.0,14.0,14.0,11*0./ c DATA UTSRT/5.00,5.00,5.00,5.00,5.00,19.00,16.00, c 1 16.00,16.00,16.00,21.00,17.00,14.40, c 2 14.40,14.4,14.4,14.4,14.4,14.4,14.40/ WRITE (6,'(1X,''SETUP3 --- LOAD IMAGE'')') KION = JION IMAX = IMX JMAX = JMX IFACT1 = 0 IFACT2 = 0 IFRST = 0 IFRST2 = 0 open(17,file='lex1.dat',status='old') read(17,*)F107 read(17,*)F107a read(17,*)year read(17,*)day read(17,*)mday close(17) MTGCM=4 DO z = 1,25,1 ZP(z) = -7. + (z-1.)*0.5 ENDDO do i=1,73 glon(i)=-180.+(i-1)*5. enddo do j=1,36 glat(j)=-87.5+(j-1)*5. enddo do i = 1,16,1 ixgcm(i) = i enddo write(6,*)'f107 f107a',f107,f107a,day C C INUMXY = 0 C Initial time of the trajectories UTO = UTSRT(1) - 1./ASTP F107AVE = F107EAV ISTART = 1 IF(IBCK .EQ. 1) 1UTO=UTO+2./ASTP ISESAV=0 IUTAV=0 DTOR=ATAN(1.)/45. RTOD=1./DTOR BOLTZ=1.38E-16 RE=6371.E+5 PI=180./RTOD DIPOLE(1)=90.+MAG(1) DIPOLE(2)=90.-MAG(3) DLON = 2.*PI/(NLONG-1) DLAT = PI/FLOAT(NLAT) DTR = 0.0174532925 DO 1052 JL = 1,NTJ,1 CHMAG1(JL) = MAG(3) CHMAG2(JL) = MAG(4) IF (ALATS(JL)/ABS(ALATS(JL)) .LT. 0.)THEN CHMAG1(JL) = MAG(1) CHMAG2(JL) = MAG(2) ENDIF PMLONG(1)=CHMAG2(JL)-180 PMLONG(2)=CHMAG2(JL) 1052 CONTINUE DO 1051 J=1,36,1 DO 1051 I=1,72,1 CUSP(I,J)=0. DRIZL(I,J)=0. 1051 CONTINUE C SET UP THE MAGNETIC COORDINATES FOR THE GRID CALL TRGTM C C SET UP LONGITUDE OF TGCM GRID IN RADIANS FOR GLOBVI DPHI = PI/36. RLON(1) = -PI DO 1053 I=2,72 1053 RLON(I) = RLON(I-1) + DPHI RLAT(1) = -PI*87.5/180. DO 1010 J = 2,36,1 1010 RLAT(J) = RLAT(J-1) + DPHI C Geographical latitude and longitude for trajectories GCNLAT(1) = -87.5 GCNLON(1) = -180. DO 1056 I=2,36 1056 GCNLAT(I) = GCNLAT(I-1) + 5. C DO 1059 I=2,73 1059 GCNLON(I) = GCNLON(I-1) + 5.0 C Make lexical reads for data vol and history reads C CALL GCMIN2 WRITE(6,*)'HISTVOL3,MXVOLS**********!!!',HISTVOL(3,MXVOLS) C If you are using the TIEGCM get new type of name for history volumes C IF(NTIE .EQ. 1)HISTCHAR = ' ' c IF(NTIE .EQ. 1)CALL CHGINP(HISTVOL,HISTCHAR,NUMCHARS) C IF(NTIE .EQ. 1)WRITE(6,*)'!!HISTCHAR ',HISTCHAR C IF(NTIE .EQ. 1)WRITE(6,*)'!!NUMCHARS ',NUMCHARS C IF(NTIE .EQ. 1)NUMCHARS = NUMCHARS - 2 C IF(NTIE .EQ. 1)WRITE(6,*)'!!NUMCHARS ',NUMCHARS C Read in the data file. 21 is the input file C CALL TGCM0 (21,MTGCM) C IMTGCM is used to put MTGCM into DIAGPK3 IMTGCM=MTGCM C CALL GCMIND (IND1, JND1, KND1) CALL DATRD C Reset lexical read parameter NID = 1 C Make the lexical reads for the call to TAIL C IPR3P = 0 C IPL3P = 0 C CALL LEXRD(IPR3P,IPL3P) F107X = F107E IYD=YEAR -1900+DAY IYDD=IYD IYEAR=YEAR DAYY=DAY IDAY=DAY IIDAY=DAY MINDAY = IDAY TMO=(DAY+16.)/30. C1=23.5*DTOR DGLON=5. DGLAT=5. DZP=0.5 WRITE(6,*)'***** ZP *******' WRITE(6,*)ZP SUD=RE*RE*DGLON*DGLAT*DTOR*DTOR PS=1./73. C IXDAY=MDAY C IUT1 should be set to 0. The program loops through NTM times IUT2=IUT1 + NTM -1 IDG=1 IZ1 = 19 IT=1 C C UT LOOP:************************************************ DO 2000 IU=IUT1,IUT2 C Set ITT for the units called in TGCMH ITT= 1 C C CALL FINDAP (AKP,1,1,ONLY,UTO,AP) DO 2003 J=1,JMX DO 2003 I=1,IMX 2003 GNLAT(I,J) = RLATM(I+2,J)*RTOD DO 2006 I=1,IMX 2006 XX(I)=GLON(I) DO 2009 K=1,KMX 2009 YY(K)=ZP(K) SUM=0. DO 1300 I=1,IMX DO 1300 J=1,JMX 1300 SUM=SUM+SUD*COS(GLAT(J)*DTOR) DO 1301 K=1,KMX QNA(K)=1.E-20 1301 CONTINUE C CALL tgcmrd(iu) UTO = iu UTO=AMOD(UTO,24.) WRITE(6,5005)IU,UTO 5005 FORMAT(' UNIVERSAL TIME LOOP NUMBER',I4' UT:',F5.1,////) MHR=UTO XMINUT=(UTO-FLOAT(MHR))*60. MIN=XMINUT UT=UTO UTZ=UTO SEC= (FLOAT(MD(IT))*24. + FLOAT(MH(IT)))*3600. 1 -43200. + FLOAT(MM(IT))*60. IF(NEWTIME .EQ. 1) 1 SEC= (FLOAT(MD(IT))*24. + FLOAT(MH(IT)))*3600. 2 + FLOAT(MM(IT))*60. WRITE(6,5010)SEC,MD(IT),MH(IT),MM(IT) 5010 FORMAT(' SEC',E10.3,' MD MH MM',3I5) IT=IT+1 CALL readtail(iu) C Write the characteristics of the history(ies) used to a file for C the first time through the UT loop IF(IU .EQ. IUT1)THEN WRITE(6,5015)(ADATA(IDM),IDM=1,3,1),((OUTPUT(IDM,JDM), 1 IDM=1,3,1),JDM=1,5,1),HNT(5),HNT(6),REL(15),REL(16),REL(17), 2 (REL(IDM),IDM=4,14,1),(REL(IDM),IDM=1,4,1) 5015 FORMAT(1X,'DESCRIPTION OF TGCM INPUTS',//,5X,'THE DATA ', 1 'VOLUME USED IS ',/,10X,3(6A1),//,5X,'THE HISTORY VOLUMES USED ' 2 ,'ARE',/,5(10X,3(6A1),2X,/),//,5X,'THE DATE FOR THESE HISTORIES=' 3 ,2(F6.1,1X),//,5X,'HEMISPHERIC POWER = ',F6.2,//,5X, 4 'CROSS CAP POTENTIAL = ',F6.2,//,5X,'BYIMF = ',F6.2,//,5X, 5 'TIDES = ',10(F7.2,','),//,5X,'MAG POLES = ',4(F7.2,',')) WRITE(6,5016)GCMUT,IYEAR,PNT(1,IZ1,2,IXGCM(15)), 1 PNT(36,IZ1,2,IXGCM(15)),PNT(1,IZ1,35,IXGCM(15)), 2 PNT(36,IZ1,35,IXGCM(15)),PNT(1,IZ1,2,IXGCM(1)), 3 PNT(36,IZ1,2,IXGCM(1)),PNT(1,IZ1,35,IXGCM(1)), 4 PNT(36,IZ1,35,IXGCM(1)) 5016 FORMAT(1X,'GCMUT',F6.1,' IYEAR',I5,'HT SH NH',4E10.3,/, 1 ' TEMPERATURE',4E10.3///) ENDIF DO J = 1,36,1 DO K = 1,25,1 DO I = 1,73,1 PNT(I,K,J,IXGCM(15)) = PNT(I,K,J,IXGCM(15))*1.E-2 + 97.*1.E3 ENDDO ENDDO ENDDO C Tides are given as the amplitudes of the (2,2),(2,3),(2,4),(2,5), C (2,6) tides and the phases of the same. Cross cap potential is in kV C By is in nT C C+++++++++++++++++++++++++++++++++++++++++++++ C CALCULATE TRAJECTORY LOCI C+++++++++++++++++++++++++++++++++++++++++++++ JUT=IU+1 UTOO=UTO DO 2012 KK=1,NTJ KUHLP2(KK)=0 2012 CONTINUE C LOOP TO CALCULATE up to TWENTY TRAJECTORIES DO 2018 KK=1,NTJ C C IF FIRST POINT INITIALISE WRITE(6,5020)UT,UTSRT(KK),IBCK 5020 FORMAT(' UT',F6.2,'UTSRT',F6.2,'IBCK',I2) IF(JUT .EQ. 1)GO TO 2021 C IF UT=SPECIAL SET IN DATA STATEMENT, INITIALISE C IF(UTO .EQ. UTSRT(KK))GO TO 2021 C OTHERWISE SKIP INITIALISATION: GO TO 2024 2021 CONTINUE WRITE(6,5025)UTSRT(KK),UT 5025 FORMAT(' UTSRT =',F6.2,' UT =',F6.2) KUHLP2(KK)=1 C Read in the initial conditions for pressure surface ZPKT(JUT,KK)=AKSTRT(KK) AZPKT(JUT,KK)=ZPKT(JUT,KK) C IF(KK .EQ. 10)WRITE(6,5030)ZPKT(JUT,KK) WRITE(6,5030)ZPKT(JUT,KK) 5030 FORMAT(' ZPKT',F8.3) C FIND NEAREST NEIGHBOUR POINTS using the initial conditions for lat lon IPK1 = (ALONS(KK)+180.)/5. + 1 IPK1 = MOD(IPK1+72,72) IF(IPK1 .EQ. 0)IPK1=IPK1 + 72 IPKHOLD = IPK1 IPK2 = MOD(IPK1+1,72) IF(IPK2 .EQ. 0)IPK2 = IPK2 + 72 JPK1 = (ALATS(KK)+92.5)/5. JPKHOLD = JPK1 JPK2 = JPK1+1 KPK1 = AKSTRT(KK) KPK2 = KPK1 + 1 C IPK3 and IPK4 are used to make transpolar calculations easy IPK3 = IPK1 IPK4 = IPK2 MFLAG = 0 C The next 2 if blocks are invoked for points close to the pole IF(JPK1 .EQ. 0) THEN JPK1 = 1 MFLAG = 1 IPK3 = MOD(IPK1+37,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 ENDIF IF(JPK2 .GT. 36) THEN JPK2 = 36 MFLAG = 2 IPK3 = MOD(IPK1+37,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 ENDIF C Find U,V,W near the trajectory path and allow for transpolar points DO 2027 JJJ=1,6 DO 2027 KKK=1,4 DO 2027 III=1,6 I=IPKHOLD + III - 3 IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 J=JPKHOLD + JJJ - 3 AJPT = 1. IF(J .GT. 36)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(J .LT. 1)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 IF(J .GT. 36)J=73-J IF(J .LT. 1)J=1-J K=KPK1+KKK-2 U(I,K,J)=PNT(I,K,J,2) * AJPT WRITE(6,5035)U(I,K,J) 5035 FORMAT(//1X,'U1 ',F8.2/) V(I,K,J)=PNT(I,K,J,3) *AJPT W(I,K,J)=PNT(I,K,J,IXGCM(14)) POO2(I,K,J)=PNT(I,K,J,4) POO(I,K,J)=PNT(I,K,J,5) C ADJUST VERTICAL WINDS using the mean molecular mass and scale height AMBAR=1./((POO2(I,K,J)/32.)+(POO(I,K,J)/16.)+((1.-POO2(I,K,J)- 1 POO(I,K,J))/28.)) SHGT=8314.*PNT(I,K,J,1)/(AMBAR*8.7) W(I,K,J)=W(I,K,J)/SHGT C Ion velocities must be used for ion trajectories IF(JEOMGN .EQ. 1)THEN C 2CND FIELD IN DATAVOL = UIM*CDEC + VIM*SDEC --- IS CONST IN HT C 3RD FIELD IN DATAVOL = VIM*CDEC - UIM*SDEC --- IS CONST IN HT C 0 LT IS IN I=3, 12 LT IN I=39, ETC. ILT = (UTO + GLON(I)/15. + 0.01)*3. + 72 ILT = MOD(ILT,72) + 3 IF(ILT .GT. 72)ILT = ILT - 72 U(I,K,J)= ( UVIDYN(ILT,J,1)+VELEW(I,J) )*AJPT*100. V(I,K,J)= ( UVIDYN(ILT,J,2)+VELNS(I,J) )*AJPT*100. ENDIF 2027 CONTINUE c calculate weighting terms for nearest neighbour points: X1 = GCNLON(IPK2) - ALONS(KK) Y1 = GCNLAT(JPK2) - ALATS(KK) IF (MFLAG .GT. 0) THEN Y1 = 5. - GCNLAT(JPK1) + ALATS(KK) X1 = 5. - ALONS(KK) + GCNLON(IPK1) ENDIF WRITE(6,5040)X1,Y1 5040 FORMAT('X1 X2',2F6.2) S1 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2 + 0.0000001) S2 = SQRT( (5.0-X1)**2 + Y1**2 + 0.000001) S3 = SQRT( X1**2 + (5.0-Y1)**2 +0.000001) S4 = SQRT( X1**2 + Y1**2+0.000001) S5 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2 +0.000001) S6 = SQRT( (5.0-X1)**2 + Y1**2 + 0.000001) S7 = SQRT( X1**2 + (5.0-Y1)**2 + 0.000001) S8 = SQRT( X1**2 + Y1**2 + 0.000001) AAA = 1.0 / ( (1./S1) + (1./S2) + (1./S3) + (1./S4)) WA = AAA/S1 WB = AAA/S2 WC = AAA/S3 WD = AAA/S4 C Exponentials better describe the variation of mixing ratio with C pressure surface, but thiss is not true for temperature. Other fields C are not strongly dependent on height so the problem does not arise Z1 = FLOAT(KPK2)-ZPKT(JUT,KK) W1(KK) = AAA/S1 *(EXP(Z1)-1.)/(EXP(1.)-1.) W2(KK) = AAA/S2 *(EXP(Z1)-1.)/(EXP(1.)-1.) W3(KK) = AAA/S3 *(EXP(Z1)-1.)/(EXP(1.)-1.) W4(KK) = AAA/S4 *(EXP(Z1)-1.)/(EXP(1.)-1.) W5(KK) = AAA/S5 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) W6(KK) = AAA/S6 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) W7(KK) = AAA/S7 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) W8(KK) = AAA/S8 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) WRITE(6,5045)W1(KK),W2(KK),W3(KK),W4(KK),W5(KK),W6(KK),W7(KK), 1 W8(KK) 5045 FORMAT(8F8.4) UPKT(JUT,KK) = (W1(KK)*U(IPK1,KPK1,JPK1)+W2(KK)* 1 U(IPK3,KPK1,JPK2) + 2 W3(KK)*U(IPK2,KPK1,JPK1)+W4(KK)* 3 U(IPK4,KPK1,JPK2) + 4 W5(KK)*U(IPK1,KPK2,JPK1)+W6(KK)* 5 U(IPK3,KPK2,JPK2) + 6 W7(KK)*U(IPK2,KPK2,JPK1)+W8(KK)* 7 U(IPK4,KPK2,JPK2)) VPKT(JUT,KK) = (W1(KK)*V(IPK1,KPK1,JPK1)+W2(KK)* 1 V(IPK3,KPK1,JPK2) + 2 W3(KK)*V(IPK2,KPK1,JPK1)+W4(KK)* 3 V(IPK4,KPK1,JPK2) + 4 W5(KK)*V(IPK1,KPK2,JPK1)+W6(KK)* 5 V(IPK3,KPK2,JPK2) + 6 W7(KK)*V(IPK2,KPK2,JPK1)+W8(KK)* 7 V(IPK4,KPK2,JPK2)) WPKT(JUT,KK) = (W1(KK)*W(IPK1,KPK1,JPK1)+W2(KK)*W(IPK3,KPK1,JPK2)+ 1 W3(KK)*W(IPK2,KPK1,JPK1)+W4(KK)*W(IPK4,KPK1,JPK2)+ 2 W5(KK)*W(IPK1,KPK2,JPK1)+W6(KK)*W(IPK3,KPK2,JPK2)+ 3 W7(KK)*W(IPK2,KPK2,JPK1)+W8(KK)*W(IPK4,KPK2,JPK2)) IF(JEOMGN .EQ. 1)THEN WPKT(JUT,KK) = 0. ENDIF IF(KK .EQ. 1) 1 WRITE(6,5050)JUT,UPKT(JUT,KK),VPKT(JUT,KK) 5050 FORMAT(' JUT',I4,'UPKT',E10.3,'VPKT',E10.3) IPKT(KK)=IPK1 IPKU(KK)=IPK1 JPKT(KK)=(ALATS(KK) +92.5)/5. JPKU(KK)=JPK1 KPKT(KK)=KPK1 KPKU(KK)=KPK1 WRITE(6,5055)IPKT(KK),JPKT(KK),KPKT(KK) 5055 FORMAT(' STARTING LOCATION. I=',I2,' J=',I2,' K=',I2) C UTPKT(JUT)=UT AUSTP = ASTP IF(IBCK .EQ. 1)AUSTP = -1. UTPKT2(JUT)=UT ALAT(JUT,KK) = ALATS(KK) ALON(JUT,KK) = ALONS(KK) AGLAT(KK) = ALATS(KK) AGLON(KK) = ALONS(KK) ALONA(JUT,KK) = ALONS(KK) ALATA(JUT,KK) = ALATS(KK) C For ion trajectories or neutral trajectories in geomagnetic coords IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN CALL GTM (ALAT(JUT,KK)*DTOR,ALON(JUT,KK)*DTOR,CHMAG1(KK)*DTOR, 1 CHMAG2(KK)*DTOR,A1,A2,DIP1,DEC1,W,1) ALON(JUT,KK) = A2/DTOR ALAT(JUT,KK) = A1/DTOR ENDIF C Dhour is set to allow for local time displacement between geomagnetic C and geographic poles ZKP1(KK) = AKSTRT(KK) DHOUR=0. CHGHR = 4.667 IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN DHOUR = CHGHR IF (ALATS(KK)/ABS(ALATS(KK)) .LT. 0.) THEN CHGHR = -8.467 DHOUR = CHGHR ENDIF ENDIF XLT(JUT,KK) = UT + ALON(JUT,KK)/15. - DHOUR IF(XLT(JUT,KK) .GT. 24.)XLT(JUT,KK)=XLT(JUT,KK)-24. IF(XLT(JUT,KK) .LT. 0. )XLT(JUT,KK)=XLT(JUT,KK)+24. IHEMI=1 IF(ALAT(JUT,KK) .LT. 0.)IHEMI=-1 2024 CONTINUE WRITE(6,5790) 5790 FORMAT(' HIT 2024') IF(KK .EQ. NTJ)GO TO 2030 2018 CONTINUE C GO TO 2034 2030 CONTINUE C find next trajectory point REI = RE/1.E2 + PNT(IPK1,KPK1,JPK1,IXGCM(15))*1.E3 DO 2033 KK=1,NTJ IF(KUHLP2(KK) .EQ. 1)GO TO 2036 C SOLVE SPHERICAL TRIANGLE TO OBTAIN COORDS OF NEW TRAJ POINT ALNGX=AGLON(KK) ALATGX=AGLAT(KK) UPOLD=UPKT(JUT-1,KK) VPOLD=VPKT(JUT-1,KK) WPOLD=WPKT(JUT-1,KK) UPNEW=UPOLD VPNEW=VPOLD WPNEW=WPOLD CZKP1=ZKP1(KK) IHEM=AGLAT(KK)/ABS(AGLAT(KK)) AHEM=FLOAT(IHEM) DO 6103 JNJK=1,2,1 ALNGX=AGLON(KK) ALATGX=AGLAT(KK) ZKP1(KK) = CZKP1 C When the terms causing the changes in mixing ratio were originally C compared with the changes in mixing ratio it was found that the C agreement was poor in some instances. This occurred because the C the position of the parcel was not sufficiently well defined by the C large time step used. As a consequence it was decided to define the C position of the parcel by a linear interpolation in four steps between C the velocity at the final position of the parcel and the velocity at C the initial position of the parcel. DO 6100 INJK=1,4,1 UPLEN = UPOLD*(5.-INJK)/4. + UPNEW*(INJK-1.)/4. VPLEN = VPOLD*(5.-INJK)/4. + VPNEW*(INJK-1.)/4. WPLEN = WPOLD*(5.-INJK)/4. + WPNEW*(INJK-1.)/4. ATMSTP = TIMSTP/4. C SOLVE SPHERICAL TRIANGLE TO OBTAIN COORDS OF NEW TRAJ POINT VMAG = SQRT(VPLEN**2 + UPLEN**2) C ANGULAR DISTANCE TRAVELLED: BETA = VMAG * ATMSTP / REI ALMONE = -(ALNGX - AGLON(KK)) * DTR UCLEN = UPLEN * COS(ALMONE) - AHEM*VPLEN * SIN(ALMONE) VCLEN = AHEM*UPLEN * SIN(ALMONE) + VPLEN * COS(ALMONE) C If backward trajectories are being calculated wind velocities affect C the parcel in the opposite sense ABCK=1. IF(IBCK .EQ. 1)ABCK=-1. ALPHA = ATAN2(ABCK*UCLEN,ABCK*VCLEN) THNEW = ASIN( SIN(ALATGX*DTR) * COS(BETA) 1 + COS(ALATGX*DTR) *SIN(BETA) * COS(ALPHA) ) C XXX = SIN(BETA) * SIN(ALPHA) YYY = COS(ALATGX*DTR) * COS(BETA) 1 - ( SIN(ALATGX*DTR) * SIN(BETA) * COS(ALPHA) ) GLNEW = ALNGX + ATAN2(XXX,YYY)/DTR ALATGX = THNEW/DTR ALNGX = GLNEW ZKP1(KK)=ZKP1(KK)+ABCK*(TIMSTP/2.0)*(WPLEN) C C IF(ALATGX .GT. 90.)ALNGX=180.+ALNGX IF(ALATGX .GT. 90.)ALATGX=180.-ALATGX IF(ALATGX .LT. -90.)ALNGX=180.+ALNGX IF(ALATGX .LT.-90.)ALATGX=-180.-ALATGX IF(ALNGX .LT. -180.)ALNGX=ALNGX+360. IF(ALNGX .GT. 180.)ALNGX=ALNGX-360. 6100 CONTINUE C Define the parcel location in i,j (longitude and latitude) grids IPK1 = (ALNGX+180.)/5. + 1 IPK1 = MOD(IPK1+72,72) IPKHOLD = IPK1 IF(IPK1 .EQ. 0)IPK1=IPK1 + 72 IPK2 = MOD(IPK1+1,72) IF(IPK2 .EQ. 0)IPK2 = IPK2 + 72 JPK1 = ((ALATGX)+92.5)/5. JPKHOLD = JPK1 JPK2 = JPK1+1 KPK1 = ZKP1(KK) KPK2 = KPK1+1 IPK3 = IPK1 IPK4 = IPK2 MFLAG = 0 C This section deals with transpolar points IF(JPK1 .EQ. 0) THEN JPK1 = 1 MFLAG = 1 IPK3 = MOD(IPK1+37,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 ENDIF IF(JPK2 .GT. 36) THEN JPK2 = 36 MFLAG = 2 IPK3 = MOD(IPK1+37,72) IF(IPK3 .EQ. 0)IPK3=IPK3 + 72 IPK4 = MOD(IPK1+36,72) IF(IPK4 .EQ. 0)IPK4 = IPK4+72 ENDIF C Find U,V,W near the trajectory path DO 2039 JJJ=1,6 DO 2039 KKK=1,4 DO 2039 III=1,6 I=IPKHOLD + III - 3 IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 J=JPKHOLD + JJJ - 3 C AJPT inverts the sense of the velocity for transpolar points AJPT = 1. C These sections deal with transpolar points IF(J .GT. 36)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(J .LT. 1)THEN I = MOD(IPK1+37,72) + 3 - III AJPT = -1. ENDIF IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 IF(J .GT. 36)J=73-J IF(J .LT. 1)J=1-J K=KPK1+KKK-2 U(I,K,J)=PNT(I,K,J,2) * AJPT V(I,K,J)=PNT(I,K,J,3) *AJPT W(I,K,J)=PNT(I,K,J,IXGCM(14)) POO2(I,K,J)=PNT(I,K,J,4) POO(I,K,J)=PNT(I,K,J,5) C ADJUST VERTICAL WINDS to w/H, where H is the scale height AMBAR=1./((POO2(I,K,J)/32.)+(POO(I,K,J)/16.)+((1.-POO2(I,K,J)- 1 POO(I,K,J))/28.)) SHGT=8314.*PNT(I,K,J,1)/(AMBAR*8.7) W(I,K,J)=W(I,K,J)/SHGT C If ion trajectories are to be calculated use ion winds IF(JEOMGN .EQ. 1)THEN C 2CND FIELD IN DATAVOL = UIM*CDEC + VIM*SDEC --- IS CONST IN HT C 3RD FIELD IN DATAVOL = VIM*CDEC - UIM*SDEC --- IS CONST IN HT C 0 LT IS IN I=3, 12 LT IN I=39, ETC. ILT = (UTO + GLON(I)/15. + 0.01)*3. + 72 ILT = MOD(ILT,72) + 3 U(I,K,J)= ( UVIDYN(ILT,J,1)+VELEW(I,J) )*AJPT*100. V(I,K,J)= ( UVIDYN(ILT,J,2)+VELNS(I,J) )*AJPT*100. ENDIF 2039 CONTINUE c calculate weighting terms for nearest neighbour points: X1 = GCNLON(IPK2) - ALNGX Y1 = GCNLAT(JPK2) - ALATGX IF (MFLAG .GT. 0) THEN Y1 = 5. - GCNLAT(JPK1) + ALATGX X1 = 5. - ALNGX + GCNLON(IPK1) ENDIF IF(MFLAG .NE. 0)WRITE(6,5060)X1,Y1,ALNGX,ALATGX,MFLAG,IPK1,JPK1 5060 FORMAT(1X,4F8.3,'MFLAG',I3,'IPK1 JPK1',2I3) S1 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2 + 0.0000001) S2 = SQRT( (5.0-X1)**2 + Y1**2 + 0.000001) S3 = SQRT( X1**2 + (5.0-Y1)**2 +0.000001) S4 = SQRT( X1**2 + Y1**2+0.000001) S5 = SQRT( (5.0-X1)**2 + (5.0-Y1)**2 +0.000001) S6 = SQRT( (5.0-X1)**2 + Y1**2 + 0.000001) S7 = SQRT( X1**2 + (5.0-Y1)**2 + 0.000001) S8 = SQRT( X1**2 + Y1**2 + 0.000001) AAA = 1.0 / ( (1./S1) + (1./S2) + (1./S3) + (1./S4)) C See earlier comment Z1 = FLOAT(KPK2)-ZKP1(KK) IF(MFLAG .NE. 0)WRITE(6,5065)S1,S2,S3,S4,S5,S6,S7,S8,Z1,KPK2 5065 FORMAT(1X,'SSSSS',9F7.3,I3) W1(KK) = AAA/S1 *(EXP(Z1)-1.)/(EXP(1.)-1.) W2(KK) = AAA/S2 *(EXP(Z1)-1.)/(EXP(1.)-1.) W3(KK) = AAA/S3 *(EXP(Z1)-1.)/(EXP(1.)-1.) W4(KK) = AAA/S4 *(EXP(Z1)-1.)/(EXP(1.)-1.) W5(KK) = AAA/S5 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) W6(KK) = AAA/S6 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) W7(KK) = AAA/S7 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) W8(KK) = AAA/S8 *(EXP(1.)-EXP(Z1))/(EXP(1.)-1.) IF(MFLAG .NE. 0)WRITE(6,5070)W1(KK),W2(KK),W3(KK),W4(KK) 5070 FORMAT(1X,'WWS',4F8.3) C Find the velocities acting on the parcel UPLEN = (W1(KK)*U(IPK1,KPK1,JPK1)+W2(KK)* 1 U(IPK3,KPK1,JPK2) + 2 W3(KK)*U(IPK2,KPK1,JPK1)+W4(KK)* 3 U(IPK4,KPK1,JPK2) + 4 W5(KK)*U(IPK1,KPK2,JPK1)+W6(KK)* 5 U(IPK3,KPK2,JPK2) + 6 W7(KK)*U(IPK2,KPK2,JPK1)+W8(KK)* 7 U(IPK4,KPK2,JPK2)) VPLEN = (W1(KK)*V(IPK1,KPK1,JPK1)+W2(KK)* 1 V(IPK3,KPK1,JPK2) + 2 W3(KK)*V(IPK2,KPK1,JPK1)+W4(KK)* 3 V(IPK4,KPK1,JPK2) + 4 W5(KK)*V(IPK1,KPK2,JPK1)+W6(KK)* 5 V(IPK3,KPK2,JPK2) + 6 W7(KK)*V(IPK2,KPK2,JPK1)+W8(KK)* 7 V(IPK4,KPK2,JPK2)) WPLEN = (W1(KK)*W(IPK1,KPK1,JPK1)+W2(KK)*W(IPK3,KPK1,JPK2)+ 1 W3(KK)*W(IPK2,KPK1,JPK1)+W4(KK)*W(IPK4,KPK1,JPK2)+ 2 W5(KK)*W(IPK1,KPK2,JPK1)+W6(KK)*W(IPK3,KPK2,JPK2)+ 3 W7(KK)*W(IPK2,KPK2,JPK1)+W8(KK)*W(IPK4,KPK2,JPK2)) C No vertical drifts for ion trajectories IF(JEOMGN .EQ. 1)THEN WPLEN = 0. ENDIF ALMONE = (ALNGX - AGLON(KK)) * DTR C Adjust so that the newly calculated velocitiesare arranged in the C zonal and meridional directions at the point of calculation. UPNEW = UPLEN * COS(ALMONE) - AHEM*VPLEN * SIN(ALMONE) VPNEW = AHEM*UPLEN * SIN(ALMONE) + VPLEN * COS(ALMONE) WPNEW = WPLEN 6103 CONTINUE UPKT(JUT,KK) = UPLEN VPKT(JUT,KK) = VPLEN WPKT(JUT,KK) = WPLEN ALAT(JUT,KK) = ALATGX ALON(JUT,KK) = ALNGX AGLAT(KK) = ALATGX AGLON(KK) = ALNGX ALONA(JUT,KK) = AGLON(KK) ALATA(JUT,KK) = AGLAT(KK) IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN CALL GTM (ALAT(JUT,KK)*DTOR,ALON(JUT,KK)*DTOR,CHMAG1(KK)*DTOR, 1 CHMAG2(KK)*DTOR,A1,A2,DIP1,DEC1,W,1) ALON(JUT,KK) = A2/DTOR ALAT(JUT,KK) = A1/DTOR ENDIF C ANGULAR DISTANCE TRAVELLED: WRITE(6,5075)WPKT(JUT-1,KK) 5075 FORMAT(' WPKT',E10.3) ZPKT(JUT,KK)=ZKP1(KK) AZPKT(JUT,KK)=ZPKT(JUT,KK) C C FIND NEAREST NEIGHBOUR POINTS C UTPKT(JUT)=UT UTPKT2(JUT) = UTPKT2(JUT-1) + AUSTP*TIMSTP/3600. IPKT(KK)=IPK1 JPKT(KK)=((ALATGX) + 92.5)/5. KPKT(KK)=KPK1 C DHOUR fixes local time differences between the geomagnetic and C geographic poles DHOUR=0. CHGHR = 4.667 IF(JEOMGN .EQ. 1 .OR. KEOMGN .EQ. 1)THEN DHOUR = CHGHR IF (ALATS(KK)/ABS(ALATS(KK)) .LT. 0.) THEN CHGHR = -8.467 DHOUR = CHGHR ENDIF ENDIF XLT(JUT,KK) = UT + ALON(JUT,KK)/15. - DHOUR IF(XLT(JUT,KK) .GT. 24.)XLT(JUT,KK)=XLT(JUT,KK)-24. IF(XLT(JUT,KK) .LT. 0. )XLT(JUT,KK)=XLT(JUT,KK)+24. IF(KK .EQ. 1) 1 WRITE (6,5080) ALON(JUT,KK),ALAT(JUT,KK),GCNLON(IPK1), 1 GCNLON(IPK2), GCNLAT(JPK1), GCNLAT(JPK2), 1 W1(KK),W2(KK),W3(KK),W4(KK) 5080 FORMAT(' DEBUG:',2F6.1,2F6.1,6F6.2) C C C 2036 CONTINUE 2033 CONTINUE 2034 CONTINUE GCMUT = UT C C Set the data boundaries around the three d trajectory point C kk is incremented up to the number of trajectories DO 2500 KK=1,NTJ,1 IZERO=IPKT(KK) JZERO=JPKT(KK) KZERO=KPKT(KK) WRITE(6,5085)IZERO,JZERO,KZERO 5085 FORMAT(' IZERO JZERO KZERO',3I4) IZ3=KZERO-1 C **** START OF LATITUDE LOOP **** C JPAL and J are coordinates in latitude space while JJ goes from 1 C to 6 in this dimension DO 2510 JPAL=JZERO-2,JZERO+3,1 J=JPAL JJ = J - JZERO + 3 JFK = J IFLAG=0 LFLAG=0 c Look out if there is a pole IF(J .GT. 36)IFLAG=1 IF(J .LT. 0)LFLAG=1 IF(IFLAG .EQ. 1)J=73-J IF(LFLAG .EQ. 1 .OR. J .EQ. 0)J=1-J GLATJ = GLAT(J) LAT=J C ONLY CALL FLOWXX FOR HIGH LATITUDES: IAUR = IIAUR IF(J.GT.12.AND.J.LT.25) IAUR=0 C IPAL and I are longitude coordinates DO 2520 IPAL=IZERO-4,IZERO+5,1 I=IPAL C Modify I if the parcel is near the pole IF(IFLAG .EQ. 1 .OR. LFLAG .EQ. 1)I=MOD(IPAL+36,72) IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 C 2CND FIELD IN DATAVOL = UIM*CDEC + VIM*SDEC --- IS CONST IN HT C 3RD FIELD IN DATAVOL = VIM*CDEC - UIM*SDEC --- IS CONST IN HT C 0 LT IS IN I=3, 12 LT IN I=39, ETC. ILT = (UTO + GLON(I)/15. + 0.01)*3. + 72 ILT = MOD(ILT,72) + 3 UI(I,J)= ( UVIDYN(ILT,J,1)+VELEW(I,J) ) * 100. VI(I,J)= ( UVIDYN(ILT,J,2)+VELNS(I,J) ) * 100. 2520 CONTINUE C I and IPAL are in longitude space ii goes from 1-6 in this direction DO 2530 IPAL=IZERO-2,IZERO+3,1 I=IPAL II=I-IZERO+3 AJPK = 1. IF(IFLAG .EQ. 1 .OR. LFLAG .EQ. 1)THEN I=MOD(IZERO+37,72) +IZERO-IPAL AJPK = -1. ENDIF C Modify I if the parcel is near the pole IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 IF(MTGCM .LT. 4)THEN DO 3522 K = KZERO-1,KZERO+2,1 PN2 = 1.00 - PNT(I,K,J,4) - PNT(I,K,J,5) PN2 = AMAX1(0.00001, PN2) BARM = 1.00 / (PNT(I,K,J,4)/32. + PNT(I,K,J,5)/16. + PN2/28.) ) RATN = (5.E-4/1.3805E-16) * BARM * EXP(-ZP(K)) / PNT(I,K,J,1) PXNO(I,K) = RATN * PNT(I,K,J,5) / 16. PXNO2(I,K) = RATN * PNT(I,K,J,4) / 32. PXNN2(I,K) = RATN * PN2 / 28. 3522 CONTINUE ENDIF IF(MTGCM .EQ. 4)THEN DO 3250 K = 1,KMX C Inputs for CHPMN ZPP(K) = ZP(K) C Height ZXX(K) = PNT(I,K,J,IXGCM(15)) + 97.*1.E5 ZPHT(K) = ZXX(K) C Number densities of O2, O and N2 PN2 = 1.00 - PNT(I,K,J,4) - PNT(I,K,J,5) PN2 = AMAX1(0.00001, PN2) BARM = 1.00 / (PNT(I,K,J,4)/32. + PNT(I,K,J,5)/16. + PN2/28.) ) RATN = (5.E-4/1.3805E-16) * BARM * EXP(-ZP(K)) / PNT(I,K,J,1) PXNO(I,K) = RATN * PNT(I,K,J,5) / 16. PXNO2(I,K) = RATN * PNT(I,K,J,4) / 32. PXNN2(I,K) = RATN * PN2 / 28. OXNO2(K) = PXNO2(I,K) OXNO(K) = PXNO(I,K) OXNN2(K) = PXNN2(I,K) C times and positions III = I JJJ = J UTX = UTO C Neutral temperature PTN = PNT(I,K,J,IXGCM(1)) TNNX(K) = PTN C Electron density XNEG(K) = PNT(I,K,J,IXGCM(12)) C Mixing ratios of O2, O and N2 PSIXX(K,1) = PNT(I,K,J,4) PSIXX(K,2) = PNT(I,K,J,5) PSIXX(K,3) = 1. - PNT(I,K,J,4) - PNT(I,K,J,5) C mean molecular mass AHMBAR(K) = (PSIXX(K,1)*32. + PSIXX(K,2)*16. + PSIXX(K,3)*28.)/ 1 (PSIXX(K,1) + PSIXX(K,2) + PSIXX(K,3)) 3250 CONTINUE ENDIF C Find vertical and line of sight to the sun column densities IF(MTGCM .EQ. 4) 1CALL CHPMN C Loop in pressure surface space DO 2540 K=KZERO-1,KZERO+2,1 KLK=K-KZERO+2 IF(I .GT. 72)I=I-72 IF(I .LT. 1)I=I+72 LON=I ILT = (UTO + GLON(I)/15. + 0.01)*3. + 72 ILT = MOD(ILT,72) + 3 IF(II .EQ. 1 .AND. JJ .EQ. 1 .AND. KLK .EQ. 1) 1 WRITE(6,5090)I,J,K,II,KK,JJ,KLK,KZERO 5090 FORMAT(' I J K II KK JJ KLK KZERO',8I3) QINH(I,J)= PQN(ILT,K,J) RO2(I,J)= DISO2(ILT,K,J) C C C Calculate the cusp and auroral heating elements in the high latitudes C This section calculates the cusp and auroral heating terms C Find electron densities and time constants for ion-neutral coupling ILM = I JLM = J IF(MTGCM .LT. 4)THEN CALL INPUTS(UT,ILM,JLM,ZPHT,CUSP(1,J),ALFA(1,J),FLUX(1,J), 1 DRIZL(1,J),IGCMCOM,IMSIS,MSIS,AF(1,J),JFK,F107E,MTGCM) IF(IAUR .EQ. 1 .AND.(AF(I,J) .GT. 0.01 .OR. CUSP(I,J) .GT. 0.01 1 .OR. DRIZL(I,J) .GT. 0.01))THEN ACONVR=3.1211E8 AALFAC=0.09 AALFAD=0.20 AFC=ACONVR*0.40/AALFAC AFD=ACONVR*0.05/AALFAD DELTE=35.E-3 DENCOL=PRHO(I,K)*PSHT(I,K) G1X=((DENCOL/4.E-6)**0.606)/ALFA(I,J) G1FX=2.11685*(G1X**2.97035)*EXP(-2.09710*G1X**0.740544)+0.587954* 1 (G1X**1.72746)*EXP(-1.37459*G1X**0.932956) QTI(1)=FLUX(I,J)*ALFA(I,J)*G1FX/(DELTE*PSHT(I,K)) G1X=((DENCOL/4.E-6)**0.606)/AALFAC G1FX=2.11685*(G1X**2.97035)*EXP(-2.09710*G1X**0.740544)+0.587954* 1 (G1X**1.72746)*EXP(-1.37459*G1X**0.932956) QTI(2)=CUSP(I,J)*AFC*AALFAC*G1FX/(DELTE*PSHT(I,K)) G1X=((DENCOL/4.E-6)**0.606)/AALFAD G1FX=2.11685*(G1X**2.97035)*EXP(-2.09710*G1X**0.740544)+0.587954* 1 (G1X**1.72746)*EXP(-1.37459*G1X**0.932956) QTI(3)=DRIZL(I,J)*AFD*AALFAD*G1FX/(DELTE*PSHT(I,K)) QCIRC=35.*1.602E-12*QTI(1)*AUREFF(K)/PRHO(I,K) QCUSP=35.*1.602E-12*QTI(2)*AUREFF(K)/PRHO(I,K) QDRIZ=35.*1.602E-12*QTI(3)*AUREFF(K)/PRHO(I,K) ELSE QCIRC=0. QCUSP=0. QDRIZ=0. ENDIF ENDIF AMBARP = 1./((PNT(I,K+1,J,4)/32.)+(PNT(I,K+1,J,5)/16.)+((1. - 1 PNT(I,K+1,J,4) - PNT(I,K+1,J,5))/28.)) AMBARM = 1./((PNT(I,K-1,J,4)/32.)+(PNT(I,K-1,J,5)/16.)+((1. - 1 PNT(I,K-1,J,4) - PNT(I,K-1,J,5))/28.)) C CALCULATE VISCOSITY, THERMAL CONDUCTION AND HEAT CAPACITY TT = PNT(I,K,J,1)**0.69 UR = 8.3144E+7 XNT = PXNN2(I,K) + PXNO2(I,K) + PXNO(I,K) PN2=PXNN2(I,K)/XNT PO=PXNO(I,K)/XNT PO2=1.-PO-PN2 C UNITS CP=ERGS/(K-CM**3) C CK=ERGS/(CM-K-SEC) C CMU=GMS/(CM-SEC) CK = ((PO2+PN2)*56.+PO*75.9)*TT CP = 0.5 * UR * (PO2*7./32.+PN2*7./28.+PO*5./16.) CMU = (PO2*4.03+PN2*3.42+PO*3.9)*1.E-6*TT C Find the extra heating terms needed for the TIGCM. Note that these C are of the form Q where dT/dZ = Q/Cp QPHOTO = 0. QIC = 0. QNMIN = 0. QRAD = 0. QION = 0. QBAND = 0. QCONT = 0. CO2COOL = 0. ANOCOOL = 0. O3PCOOL = 0. IF(MTGCM .EQ. 4)THEN CALL INPUTS(UT,ILM,JLM,ZPHT,CUSP(1,J),ALFA(1,J),FLUX(1,J), 1 DRIZL(1,J),IGCMCOM,IMSIS,MSIS,AF(1,J),JFK,F107E,MTGCM) QINH(I,J) = 0. QCIRC = 0. QCUSP = 0. QDRIZ = 0. IF(ISMIN .EQ. 1 .OR. JION .EQ. 1)THEN C Initialize all variables to 0. DIFN4S = 0. EDDN4S = 0. SN4S4 = 0. SN4S5 = 0. SN4S7 = 0. SN4S8 = 0. SN4SK2 = 0. SN4SK6 = 0. SN4SK8 = 0. SN4SA1 = 0. SN4SA3 = 0. SN4SE = 0. ALN4S1 = 0. ALN4S3 = 0. ALN4K4 = 0. DIFNO = 0. EDDNO = 0. SNOS1 = 0. SNOS2 = 0. ALNOS3 = 0. ALNOS6 = 0. ALNOS8 = 0. ALNOS9 = 0. ALNOK5 = 0. DIFOP = 0. SOP1 = 0. ALOP1 = 0. ALOP2 = 0. ALOP3 = 0. QSINP = 0. QSOP = 0. QSO2P = 0. QSN2P = 0. QSNOP = 0. QAINP = 0. QAOP = 0. QAO2P = 0. QAN2P = 0. QANOP = 0. SNP1 = 0. ALNP1 = 0. ALNP2 = 0. ALNP3 = 0. ALN2P1 = 0. ALN2P2 = 0. ALN2P3 = 0. SO2P1 = 0. SO2P2 = 0. SO2P3 = 0. ALO2P1 = 0. ALO2P2 = 0. ALO2P3 = 0. SNOP1 = 0. SNOP2 = 0. SNOP3 = 0. SNOP4 = 0. ALNOP1 = 0. PHIAV = 0. OPV = 0. ADVOPN = 0. ADVOPI = 0. DO 3520 KL = 1,3,1 C Mean molecular mass AEMBAR(KL) = AHMBAR(K-2+KL) DAMBAR(KL) = AHMBAR(K-1+KL) - AHMBAR(K-3+KL) C Mixing ratios APSIO(KL) = PSIXX(K-2+KL,2) APSIO2(KL) = PSIXX(K-2+KL,1) DPSIO(KL) = PSIXX(K-1+KL,2) - PSIXX(K-3+KL,2) DPSIO2(KL) = PSIXX(K-1+KL,1) - PSIXX(K-3+KL,1) C Temperatures TEMIT(KL) = PNT(I,K-2+KL,J,IXGCM(1)) C Mixing ratios N(4S) NO and N(2D) PSMIN(KL,1) = PNT(I,K-2+KL,J,IXGCM(6)) PSMIN(KL,2) = PNT(I,K-2+KL,J,IXGCM(7)) PSMIN(KL,3) = PNT(I,K-2+KL,J,IXGCM(9)) C Pressure surface ZFAC(KL) = ZP(K-2+KL) IF(I .EQ. 36 .AND. J .EQ. 2)THEN WRITE(6,5095)(PSMIN(KL,JL),JL=1,3,1),ZFAC(KL) 5095 FORMAT(1X,' PSMIN ZFAC',4E10.3) ENDIF 3520 CONTINUE ENDIF C Set psis for the pressure surface PSI(1) = PSIXX(K,1) PSI(2) = PSIXX(K,2) PSI(3) = PSIXX(K,3) F107G = F107E C mean molecular mass EMBAR = (PSIXX(K,1)*32. + PSIXX(K,2)*16. + PSIXX(K,3)*28.)/ 1 (PSIXX(K,1) + PSIXX(K,2) + PSIXX(K,3)) C Number densities of O2, O and N2 XNA(1) = PXNO2(I,K) XNA(2) = PXNO(I,K) XNA(3) = PXNN2(I,K) C Ion and electron temperatures TIII = PNT(I,K,J,IXGCM(10)) TEEE = PNT(I,K,J,IXGCM(11)) C Use these temperatures to find various reaction rates for ions and C neutrals. PTN = PNT(I,K,J,IXGCM(1)) CALL RATECOF(TIII,TEEE,PTN,K) C Calculate solar UV and EUV heating PSIN4 = PNT(I,K,J,IXGCM(6)) C WRITE(6,5111)PSIN4,PTN,K,PSI,F107G,F107AVE,XNA,III,JJJ,IIDAY,UTX 5111 FORMAT(1X,2E10.3,I4,8E10.3,3I4,F6.1) CALL SOLHEAT(PTN,PSIN4,K) C Calculate the mixing ratio of NO PSINO = PNT(I,K,J,IXGCM(7)) C Find the minor species number densities for XTRAN BARM = 1./(PSIXX(K,1)/32. + PSIXX(K,2)/16. + PSIXX(K,3)/28.) RATN = (5.E-4/1.3805E-16)*BARM*EXP(-ZP(K))/PTN XNN4S = RATN*PNT(I,K,J,IXGCM(6))/14. XNNO = RATN*PNT(I,K,J,IXGCM(7))/30. XNN2D = RATN*PNT(I,K,J,IXGCM(9))/14. C WRITE(6,5100)XNN2D,XNNO,XNN4S 5100 FORMAT(' XN2D XNNO XN4S',3E10.3) C Mixing ratios of the minor species PSIN4 = PNT(I,K,J,IXGCM(6)) PSINO = PNT(I,K,J,IXGCM(7)) PSIN2 = PNT(I,K,J,IXGCM(9)) BTN = PTN BTI = TIII BTE = TEEE YIOP = PNT(I,K,J,IXGCM(8)) C Find the ion densities -1 production by em radiation CALL IONDENS(K,PTN,PSINO) C Find the ion densities -2 production by chemistry and precipitating C protons and electrons CALL IONMOD(IAUR,AF(I,J),CUSP(I,J),DRIZL(I,J),FLUX(I,J), 1 ALFA(I,J),YIOP,K) C Now calculate ion heating CALL IONHEAT(PTN,PSIN4,PSIN2,PSINO,XNE(I,K),EMBAR,K) C Photoelectric electron heating CALL PHOTELEC(PTN,QPHOTO,EMBAR,K) C Minor species neutral heating CALL NOHEAT(PTN,PSIN4,PSIN2,PSINO,TEEE,XNE(I,K),K) C Heating as a result of collisions between electrons and ions C and electrons and neutrals CALL COLHEAT(BTN,BTI,BTE,XNE(I,K),K) C Infrared cooling CALL IRCOOL(K,PTN,PSINO,CP) CO2COOL = -CO2COOL ANOCOOL = -ANOCOOL O3PCOOL = -O3PCOOL QRAD = -QRAD IF(II .EQ. 1 .AND. JJ .EQ. 1 .AND. KLK .EQ. 1) 1 WRITE(6,5105)CP,QPHOTO, QHEAT,QION,QCONT,QBAND,QIC,QNMIN, 1 CO2COOL,ANOCOOL,O3PCOOL,QRAD 5105 FORMAT(1X,'CP',E10.3,' QPHOTO',E10.3,' QHEAT',E10.3,/,' QION', 1 E10.3,'QCONT',E10.3,' QBAND',E10.3,'QIC',E10.3,' QNMIN', 2 E10.3,/,'CO2COOL',E10.3,' NOCOOL',E10.3,' O3PCOOL',E10.3, 3 ' QRAD',E10.3) C WRITE(6,5110)RJ,RJIC,RJNO,RO2(I,J) 5110 FORMAT(1X,'RJ',E10.3,' RJIC',E10.3,'RJNO',E10.3,'RO2',E10.3) RO2(I,J) = RJ + RJIC + RJNO QINH(I,J) = QHEAT C End of block C Make the minor species calculations if required IF(ISMIN .EQ. 1)THEN C N(4S) production and loss C Molecular diffusion CALL MINDIF(1,DIFN4S,ZFAC) C Eddy diffusion CALL MINEDD(1,EDDN4S,ZFAC) C Chemistry CALL CHEMN4(ZFAC,K) C NO production and loss C Molecular diffusion CALL MINDIF(2,DIFNO,ZFAC) C Eddy diffusion CALL MINEDD(2,EDDNO,ZFAC) C Chemistry CALL CHEMNO(ZFAC,K) C N(2D) production and loss. C Chemistry CALL CHEMN2D(ZFAC,K) IF((II .EQ. 3) .AND. (JJ .EQ. 3) .AND. (KLK .EQ. 2))THEN WRITE(6,5095)(PSMIN(2,JL),JL=1,3,1),ZFAC(2) WRITE(6,5115)DIFN4S,EDDN4S,SN4S4,SN4S5,SN4S7, 1 SN4S8,SN4SK2,SN4SK6,SN4SK8,SN4SA1,SN4SA3,SN4SE, 2 ALN4S1,ALN4S3,ALN4K4 5115 FORMAT(1X,/,' *****N4S TERMS*****',/,' DIFN4S',E10.3,'EDDN4S', 1 E10.3,'SB4',E10.3,'SB5',E10.3,/,' SB7',E10.3,'SB8',E10.3, 2 'SK2',E10.3,'SK6',E10.3,'SK8',E10.3,/,' SA1',E10.3,'SA3', 3 E10.3,'SE',E10.3,/,' LB1',E10.3,'LB3',E10.3,'LK4',E10.3) WRITE(6,5120)DIFNO,EDDNO,SNOS1,SNOS2,ALNOS3,ALNOS6,ALNOS8, 1 ALNOS9,ALNOK5 5120 FORMAT(1X,/,' ***NO TERMS***',/,' DIFNO',E10.3,'EDDNO',E10.3, 1 'SB1',E10.3,'SB2',E10.3,/,' LB3',E10.3,'LB6',E10.3,'LB8', 2 E10.3,'LB9',E10.3,'LK5',E10.3) WRITE(6,5121)SN2DK3,SN2DA1,SN2DA3,SN2DE,ALN2D2,ALN2D4, 1 ALN2D5,ALN2D6, 1 ALN2D7,ALN2DK10 5121 FORMAT(1X,/,' ***N(2D) TERMS***',/,' sk3',E10.3,' sa1',E10.3, 1 ' sa3',E10.3,' sion',E10.3,/,' lb2',E10.3, 2 ' lb4',E10.3,' lb5', 2 E10.3,' lb6',E10.3,' lb7',E10.3,' lk10',E10.3) ENDIF ENDIF IF(JION .EQ. 1)THEN XIOJP(1) = PNT(I,K-1,J,IXGCM(8)) XIOJP(2) = PNT(I,K,J,IXGCM(8)) XIOJP(3) = PNT(I,K+1,J,IXGCM(8)) DO 3525 KL = 1,3,1 AEMBAR(KL) = AHMBAR(K-2+KL) CTN(KL) = PNT(I,K-2+KL,J,IXGCM(1)) CTE(KL) = PNT(I,K-2+KL,J,IXGCM(11)) CTI(KL) = PNT(I,K-2+KL,J,IXGCM(10)) 3525 CONTINUE DO 3527 LG = 1,3,1 JGG = J - 2 + LH IF(JGG .LT. 1)JGG = 1 - JGG IF(JGG .GT. 36)JGG = 73 - JGG FLAT(LG) = RLAT(JGG) DO 3527 LH = 1,3,1 IGG = I - 2 + LG IF(IGG .LT. 1)IGG = IGG + 72 IF(IGG .GT. 72)IGG = IGG - 72 JGG = J - 2 + LH IF((JGG .LT. 1) .OR. (JGG .GT. 36))IGG = IGG + 36 IF(IGG .LT. 1)IGG = IGG + 72 IF(IGG .GT. 72)IGG = IGG - 72 IF(JGG .LT. 1)JGG = 1 - JGG IF(JGG .GT. 36)JGG = 73 - JGG C IF(J .EQ. 3)WRITE(6,5917)IGG,JGG,VELEW(I,J),VELNS(I,J), C 1 UI(I,J),VI(I,J),I,J 5917 FORMAT(1X,'IGG JGG',2I4,'UEW VNS UI VI',4E10.3,'I J',2I4) CXNOP(LG,LH) = PNT(IGG,K,JGG,IXGCM(8)) DO 3527 LI = 1,3,1 UNEU(LG,LI,LH) = PNT(IGG,K-2+LI,JGG,IXGCM(2))*AJPK VNEU(LG,LI,LH) = PNT(IGG,K-2+LI,JGG,IXGCM(3))*AJPK WNEU(LG,LI,LH) = PNT(IGG,K-2+LI,JGG,IXGCM(14)) 3527 CONTINUE UII = UI(I,J)*1.E-2*AJPK VII = VI(I,J)*1.E-2*AJPK CALL MAGFIEL(I,J) WII = WIDYN(ILT,J) - (BX(2,2,2)*VELEW(I,J)*AJPK + BY(2,2,2)* 1 VELNS(I,J)*AJPK)/BZ(2,2,2) CALL OPTERM(XIOJP,CTE,CTI,CTN,XNN2D,K) IF((II .EQ. 3) .AND. (JJ .EQ. 3) .AND. (KLK .EQ. 2))THEN WRITE(6,5125)XINPJ,XIN2PJ,XIOPJ,XIO2PJ,XINOPJ,DIFOP,SOP1, 1 ALOP1,ALOP2,ALOP3,ADVOPN,ADVOPI 5125 FORMAT(1X,/,' *****ION DENSITIES O+ TERMS****',/,' XINP',E10.3, 1 'XIN2P',E10.3,'XIOP',E10.3,'XIO2P',E10.3,'XINOP',E10.3,/, 2 ' DIFOP',E10.3,'SOP1',E10.3,'ALOP1',E10.3,'ALOP2',E10.3, 3 'ALOP3',E10.3,/,' ADVOPN',E10.3,'ADVOPI',E10.3) WRITE(6,5130)QSINP,QSOP,QSO2P,QSN2P,QSNOP,QAINP,QAOP,QAO2P, 1 QAN2P,QANOP 5130 FORMAT(1X,/,' ***EM AURORAL IONIZATION***',/,' QSINP',E10.3, 1 'QSOP',E10.3,'QSO2P',E10.3,'QSN2P',E10.3,'QSNOP',E10.3,/, 2 ' QAINP',E10.3,'QAOP',E10.3,'QAO2P',E10.3,'QAN2P',E10.3, 3 'QANOP',E10.3) WRITE(6,5135)SNP1,ALNP1,ALNP2,ALNP3,ALN2P1,ALN2P2,ALN2P3, 1 SO2P1,SO2P2,SO2P3,ALO2P1,ALO2P2,ALO2P3,SNOP1, 2 SNOP2,SNOP3,SNOP4,ALNOP1,PHIAV,OPV 5135 FORMAT(1X,'*****ION TERMS*****',/,' SNP1',E10.3,'ALNP1',E10.3, 1 'ALNP2',E10.3,'ALNP3',E10.3,'ALN2P1',E10.3,'ALN2P2', 2 E10.3,'ALN2P3',E10.3,/,' SO2P1',E10.3,'SO2P2',E10.3, 3 'SO2P3',E10.3,'ALO2P1',E10.3,'ALO2P2',E10.3,'ALO2P3', 4 E10.3,/,'SNOP1',E10.3,'SNOP2',E10.3,'SNOP3',E10.3, 5 'SNOP4',E10.3,'ALNOP1',E10.3,/,' PHIAV',E10.3,'OPV', 6 E10.3) ENDIF ENDIF ENDIF IF(II .EQ. 1 .AND. JJ .EQ. 1 .AND. KLK .EQ. 1) 1 WRITE(6,*)'T0',T0(K) TOO(K)=T0(K) IF(MTGCM .EQ. 4)TOO(K)=0. C Fill Lamda array with values DO K1K = 1,KMX ALAMDAT(II,KLK,JJ,K1K) = ALAM1HLD(II,K1K) ENDDO C C Fill the array with values for term calculation in diags KFLAG(II,KLK,JJ)=0 IF(IFLAG .EQ. 1 .OR. LFLAG .EQ. 1)KFLAG(II,KLK,JJ)=1 C The first 3 terms are the ion-neutral coupling coefficients BUFDAT(II,KLK,JJ,1) = PLXX(I,K) BUFDAT(II,KLK,JJ,2) = PLYY(I,K) BUFDAT(II,KLK,JJ,3) = PLXY(I,K) C Viscosity coefficient BUFDAT(II,KLK,JJ,4) = CMU C Heat capacity at constant pressure BUFDAT(II,KLK,JJ,5) = CP C Solar heat input BUFDAT(II,KLK,JJ,6) = QINH(I,J) C Zonal and meridional ion velocities BUFDAT(II,KLK,JJ,7) = UI(I,J)*1.E-2*AJPK BUFDAT(II,KLK,JJ,8) = VI(I,J)*1.E-2*AJPK C Thermal conductivity coefficient BUFDAT(II,KLK,JJ,9) = CK C PNT(...1) is temperature, 2 is zonal wind, 3 is meridional wind, C 4 is the mixing ratio of O2, 5 is the mixing ratio of O BUFDAT(II,KLK,JJ,10) = PNT(I,K,J,1) BUFDAT(II,KLK,JJ,11) = PNT(I,K,J,2)*AJPK BUFDAT(II,KLK,JJ,12) = PNT(I,K,J,3)*AJPK BUFDAT(II,KLK,JJ,13) = PNT(I,K,J,4) BUFDAT(II,KLK,JJ,14) = PNT(I,K,J,5) C Vertical velocity divided by the scale height BUFDAT(II,KLK,JJ,15) = W(I,K,J) C The next PNT value refers to the height of the pressure surface BUFDAT(II,KLK,JJ,16) = PNT(I,K,J,IXGCM(15)) C Circle, cusp and drizzle heating BUFDAT(II,KLK,JJ,17) = QCIRC BUFDAT(II,KLK,JJ,18) = QCUSP BUFDAT(II,KLK,JJ,19) = QDRIZ C Rate coefficient for the photodissociation of molecular oxygen BUFDAT(II,KLK,JJ,20) = RO2(I,J) C electron density and potential BUFDAT(II,KLK,JJ,21) = XNE(I,K) BUFDAT(II,KLK,JJ,22) = POTT(I,J) C Heating by photoelectrons BUFDAT(II,KLK,JJ,23) = QPHOTO C Heating due to ion chemistry BUFDAT(II,KLK,JJ,24) = QIC C Heating due to minor species chemistry BUFDAT(II,KLK,JJ,25) = QNMIN C Radiative cooling BUFDAT(II,KLK,JJ,26) = QRAD C Heating from ionization by solar EUV radiation BUFDAT(II,KLK,JJ,27) = QION C Heating in the Schumann-Runge continuum BUFDAT(II,KLK,JJ,28) = QCONT C Heating in the Schumann-Runge bands BUFDAT(II,KLK,JJ,29) = QBAND C CO2 cooling BUFDAT(II,KLK,JJ,30) = CO2COOL C NO cooling BUFDAT(II,KLK,JJ,31) = ANOCOOL C O(3P) cooling BUFDAT(II,KLK,JJ,32) = O3PCOOL C Electron ion and electron neutral collisional heating BUFDAT(II,KLK,JJ,33) = QENEI BUFDAT(II,KLK,JJ,34) = QEN BUFDAT(II,KLK,JJ,35) = QEI C It would be nice to have ion and electron temperatures, so BUFDAT(II,KLK,JJ,36) = TIII BUFDAT(II,KLK,JJ,37) = TEEE C write out minor species and ion production terms on the iz1 surface IF ((MTGCM .EQ. 4) .AND. ((ISMIN .EQ. 1) .OR. (JION .EQ. 1)))THEN DO 3540 KNK = 1,91,1 DUFDAT(II,KLK,JJ,KNK) = 0. 3540 CONTINUE C Temperature DUFDAT(II,KLK,JJ,1) = BUFDAT(II,KLK,JJ,10) C Neutral winds DUFDAT(II,KLK,JJ,2) = BUFDAT(II,KLK,JJ,11) DUFDAT(II,KLK,JJ,3) = BUFDAT(II,KLK,JJ,12) C Ion winds (zonal,meridional) DUFDAT(II,KLK,JJ,4) = BUFDAT(II,KLK,JJ,7) DUFDAT(II,KLK,JJ,5) = BUFDAT(II,KLK,JJ,8) C O2 Psi, O Psi DUFDAT(II,KLK,JJ,6) = BUFDAT(II,KLK,JJ,13) DUFDAT(II,KLK,JJ,7) = BUFDAT(II,KLK,JJ,14) C Vertical Winds DUFDAT(II,KLK,JJ,8) = BUFDAT(II,KLK,JJ,15) C Electron Densities DUFDAT(II,KLK,JJ,9) = BUFDAT(II,KLK,JJ,15) C N(4s) DUFDAT(II,KLK,JJ,10) = PNT(I,K,J,IXGCM(6)) C NO DUFDAT(II,KLK,JJ,11) = PNT(I,K,J,IXGCM(7)) C N(2D) C DUFDAT(II,KLK,JJ,12) = PNT(I,K,J,IXGCM(9)) DUFDAT(II,KLK,JJ,12) = PN2D C N+ DUFDAT(II,KLK,JJ,13) = XINPJ C N2+ DUFDAT(II,KLK,JJ,14) = XIN2PJ C O+ DUFDAT(II,KLK,JJ,15) = XIOPJ C O2+ DUFDAT(II,KLK,JJ,16) = XIO2PJ C NO+ DUFDAT(II,KLK,JJ,17) = XINOPJ C N(4S) diffusion DUFDAT(II,KLK,JJ,18) = DIFN4S C N(4S) eddy diffusion DUFDAT(II,KLK,JJ,19) = EDDN4S C Minor species chemistry DUFDAT(II,KLK,JJ,20) = SN4S4 DUFDAT(II,KLK,JJ,21) = SN4S5 DUFDAT(II,KLK,JJ,22) = SN4S7 DUFDAT(II,KLK,JJ,23) = SN4S8 DUFDAT(II,KLK,JJ,24) = SN4SK2 DUFDAT(II,KLK,JJ,25) = SN4SK6 DUFDAT(II,KLK,JJ,26) = SN4SK8 DUFDAT(II,KLK,JJ,27) = SN4SA1 DUFDAT(II,KLK,JJ,28) = SN4SA3 DUFDAT(II,KLK,JJ,29) = SN4SE DUFDAT(II,KLK,JJ,30) = ALN4S1 DUFDAT(II,KLK,JJ,31) = ALN4S3 DUFDAT(II,KLK,JJ,32) = ALN4K4 C NO diffusion DUFDAT(II,KLK,JJ,33) = DIFNO C NO eddy diffusion DUFDAT(II,KLK,JJ,34) = EDDNO C NO chemistry DUFDAT(II,KLK,JJ,35) = SNOS1 DUFDAT(II,KLK,JJ,36) = SNOS2 DUFDAT(II,KLK,JJ,37) = ALNOS3 DUFDAT(II,KLK,JJ,38) = ALNOS6 DUFDAT(II,KLK,JJ,39) = ALNOS8 DUFDAT(II,KLK,JJ,40) = ALNOS9 DUFDAT(II,KLK,JJ,41) = ALNOK5 C N(2D) chemistry DUFDAT(II,KLK,JJ,42) = SN2DK3 DUFDAT(II,KLK,JJ,43) = SN2DA1 DUFDAT(II,KLK,JJ,44) = SN2DA3 DUFDAT(II,KLK,JJ,45) = SN2DE DUFDAT(II,KLK,JJ,46) = ALN2D2 DUFDAT(II,KLK,JJ,47) = ALN2D4 DUFDAT(II,KLK,JJ,48) = ALN2D5 DUFDAT(II,KLK,JJ,49) = ALN2D6 DUFDAT(II,KLK,JJ,50) = ALN2D7 DUFDAT(II,KLK,JJ,51) = ALN2DK10 C O+ diffusion DUFDAT(II,KLK,JJ,52) = DIFOP DUFDAT(II,KLK,JJ,53) = SOP1 DUFDAT(II,KLK,JJ,54) = ALOP1 DUFDAT(II,KLK,JJ,55) = ALOP2 DUFDAT(II,KLK,JJ,56) = ALOP3 C Ionizations DUFDAT(II,KLK,JJ,57) = QSINP DUFDAT(II,KLK,JJ,58) = QSOP DUFDAT(II,KLK,JJ,59) = QSO2P DUFDAT(II,KLK,JJ,60) = QSN2P DUFDAT(II,KLK,JJ,61) = QSNOP DUFDAT(II,KLK,JJ,62) = QAINP DUFDAT(II,KLK,JJ,63) = QAOP DUFDAT(II,KLK,JJ,64) = QAO2P DUFDAT(II,KLK,JJ,65) = QAN2P DUFDAT(II,KLK,JJ,66) = QANOP C Ion chemistry DUFDAT(II,KLK,JJ,67) = SNP1 DUFDAT(II,KLK,JJ,68) = ALNP1 DUFDAT(II,KLK,JJ,69) = ALNP2 DUFDAT(II,KLK,JJ,70) = ALNP3 DUFDAT(II,KLK,JJ,71) = ALN2P1 DUFDAT(II,KLK,JJ,72) = ALN2P2 DUFDAT(II,KLK,JJ,73) = ALN2P3 DUFDAT(II,KLK,JJ,74) = SO2P1 DUFDAT(II,KLK,JJ,75) = SO2P2 DUFDAT(II,KLK,JJ,76) = SO2P3 DUFDAT(II,KLK,JJ,77) = ALO2P1 DUFDAT(II,KLK,JJ,78) = ALO2P2 DUFDAT(II,KLK,JJ,79) = ALO2P3 DUFDAT(II,KLK,JJ,80) = SNOP1 DUFDAT(II,KLK,JJ,81) = SNOP2 DUFDAT(II,KLK,JJ,82) = SNOP3 DUFDAT(II,KLK,JJ,83) = SNOP4 DUFDAT(II,KLK,JJ,84) = SNOP5 DUFDAT(II,KLK,JJ,85) = SNOP6 DUFDAT(II,KLK,JJ,86) = SNOP7 DUFDAT(II,KLK,JJ,87) = ALNOP1 DUFDAT(II,KLK,JJ,88) = PHIAV DUFDAT(II,KLK,JJ,89) = OPV DUFDAT(II,KLK,JJ,90) = ADVOPN DUFDAT(II,KLK,JJ,91) = ADVOPI EUFDAT(II,KLK,JJ,1) = UPADVOP EUFDAT(II,KLK,JJ,2) = PNT(I,K,J,IXGCM(10)) EUFDAT(II,KLK,JJ,3) = PNT(I,K,J,IXGCM(11)) EUFDAT(II,KLK,JJ,4) = UPNEU ENDIF C 2540 CONTINUE IF(JEOMGN .EQ. 1)THEN DO 2544 K = 1,25,1 QDENS(II,K,JJ) = XNE(I,K) 2544 CONTINUE C Calculate the height of the pressure surface DO 2548 KMK = 1,25,1 CUFDAT(II,KMK,JJ) = PNT(I,KMK,J,IXGCM(15)) 2548 CONTINUE ENDIF 2530 CONTINUE 2510 CONTINUE Z3PKT(JUT,KK) = (W1(KK)*BUFDAT(3,2,3,16)+ 1 W2(KK)*BUFDAT(3,2,4,16)+ 2 W3(KK)*BUFDAT(4,2,3,16)+W4(KK)*BUFDAT(4,2,4,16)) C If ion trajectories are being calculated the find the ion heating C terms that were calculated in inputs along the trajectory IF(JEOMGN .EQ. 1)THEN C INSERT SPECIFIC HEAT COEFFICIENT FOR ATOMIC OXYGEN ATMOSPHERE C TO BE USED IN COMPUTING HEATING RATES FOR ION PARCEL CPI = 3.451E-16 C Ion temperature TIPKT(JUT,KK) = (W1(KK)*TION(3,1,3)+ 1 W2(KK)*TION(3,1,4)+ 2 W3(KK)*TION(4,1,3)+W4(KK)*TION(4,1,4)+ 3 W5(KK)*TION(3,2,3)+W6(KK)*TION(3,2,4)+ 4 W7(KK)*TION(4,2,3)+W8(KK)*TION(4,2,4)) C Loss as a result of ion-neutral heat exchange QIPKT(JUT,2,KK) = (W1(KK)*QHEX(3,1,3)/BUFDAT(3,2,3,21)+ 1 W2(KK)*QHEX(3,1,4)/BUFDAT(3,2,4,21)+ 2 W3(KK)*QHEX(4,1,3)/BUFDAT(4,2,3,21)+ 3 W4(KK)*QHEX(4,1,4)/BUFDAT(4,2,4,21)+ 4 W5(KK)*QHEX(3,2,3)/BUFDAT(3,3,3,21)+ 5 W6(KK)*QHEX(3,2,4)/BUFDAT(3,3,4,21)+ 6 W7(KK)*QHEX(4,2,3)/BUFDAT(4,3,3,21)+ 7 W8(KK)*QHEX(4,2,4)/BUFDAT(4,3,4,21)) C CALCULATE ION HEAT RATE (ERGS CM-3 SEC-1)FOR ALL PRESS SFCS DO 2515 K=1,25,1 QI3PKT(JUT,KK,K) = (W1(KK)*Q2FRIC(3,K,3)+ 1 W2(KK)*Q2FRIC(3,K,4)+ 2 W3(KK)*Q2FRIC(4,K,3)+W4(KK)*Q2FRIC(4,K,4)) QPPHT(JUT,KK,K) = W1(KK)*Q2FRIC(3,K,3)/QDENS(3,K,3) + 1 W2(KK)*Q2FRIC(3,K,4)/QDENS(3,K,4) + 2 W3(KK)*Q2FRIC(4,K,3)/QDENS(4,K,3) + 3 W4(KK)*Q2FRIC(4,K,4)/QDENS(4,K,4) C CALCULATE ALT FOR EACH STD PRESS SFC AT EACH TRAJECTORY POINT Z2PKT(JUT,KK,K) = (W1(KK)*CUFDAT(3,K,3)+ 1 W2(KK)*CUFDAT(3,K,4)+ 2 W3(KK)*CUFDAT(4,K,3)+W4(KK)*CUFDAT(4,K,4)) 2515 CONTINUE C Frictional heating of the ions per particle QIPKT(JUT,1,KK) = (W1(KK)*QFRIC(3,1,3)/BUFDAT(3,2,3,21)+ 1 W2(KK)*QFRIC(3,1,4)/BUFDAT(3,2,4,21)+ 2 W3(KK)*QFRIC(4,1,3)/BUFDAT(4,2,3,21)+ 3 W4(KK)*QFRIC(4,1,4)/BUFDAT(4,2,4,21)+ 4 W5(KK)*QFRIC(3,2,3)/BUFDAT(3,3,3,21)+ 5 W6(KK)*QFRIC(3,2,4)/BUFDAT(3,3,4,21)+ 6 W7(KK)*QFRIC(4,2,3)/BUFDAT(4,3,3,21)+ 7 W8(KK)*QFRIC(4,2,4)/BUFDAT(4,3,4,21)) QI2PKT(JUT,2,KK) = (W1(KK)*QHEX(3,1,3)/BUFDAT(3,2,3,5)+ 1 W2(KK)*QHEX(3,1,4)/BUFDAT(3,2,4,5)+ 2 W3(KK)*QHEX(4,1,3)/BUFDAT(4,2,3,5)+ 3 W4(KK)*QHEX(4,1,4)/BUFDAT(4,2,4,5)+ 4 W5(KK)*QHEX(3,2,3)/BUFDAT(3,3,3,5)+ 5 W6(KK)*QHEX(3,2,4)/BUFDAT(3,3,4,5)+ 6 W7(KK)*QHEX(4,2,3)/BUFDAT(4,3,3,5)+ 7 W8(KK)*QHEX(4,2,4)/BUFDAT(4,3,4,5)) C Electron densities ENDENS(JUT,KK) = (W1(KK)*BUFDAT(3,2,3,21)+ 1 W2(KK)*BUFDAT(3,2,4,21)+ 2 W3(KK)*BUFDAT(4,2,3,21)+ 3 W4(KK)*BUFDAT(4,2,4,21)+ 4 W5(KK)*BUFDAT(3,3,3,21)+ 5 W6(KK)*BUFDAT(3,3,4,21)+ 6 W7(KK)*BUFDAT(4,3,3,21)+ 7 W8(KK)*BUFDAT(4,3,4,21)) C Ion-neutral velocity differences VIVN(JUT,2,KK) = (W1(KK)*WNSHR(3,1,3)+ 1 W2(KK)*WNSHR(3,1,4)+ 2 W3(KK)*WNSHR(4,1,3)+ 3 W4(KK)*WNSHR(4,1,4)+ 4 W5(KK)*WNSHR(3,2,3)+ 5 W6(KK)*WNSHR(3,2,4)+ 6 W7(KK)*WNSHR(4,2,3)+ 7 W8(KK)*WNSHR(4,2,4)) C Frictional heating in degrees/second QI2PKT(JUT,1,KK) = QIPKT (JUT,1,KK)/CPI C Total ion heating QIPKT(JUT,3,KK) = QIPKT(JUT,1,KK) + QIPKT(JUT,2,KK) QI2PKT(JUT,3,KK) = QI2PKT(JUT,1,KK) + QI2PKT(JUT,2,KK) ENDIF WRITE(6,*)'KK',KK,'pre5140' DO K1K = 1,KMX ALAMPKT(JUT,K1K,KK) = (W1(KK)*ALAMDAT(3,1,3,K1K)+ 1 W2(KK)*ALAMDAT(3,1,4,K1K)+ 2 W3(KK)*ALAMDAT(4,1,3,K1K)+ 3 W4(KK)*ALAMDAT(4,1,4,K1K)+ 4 W5(KK)*ALAMDAT(3,2,3,K1K)+ 5 W6(KK)*ALAMDAT(3,2,4,K1K)+ 6 W7(KK)*ALAMDAT(4,2,3,K1K)+ 7 W8(KK)*ALAMDAT(4,2,4,K1K)) ENDDO WRITE(6,*)'KK',KK,'Still pre5140' WRITE(6,5140)(BUFDAT(3,2,3,L),L=1,32) 5140 FORMAT(1X,'BUFDAT',7E10.3) CALL DIAGPKG3 WRITE(6,*)'Post Diag' 2500 CONTINUE 2000 CONTINUE C C END UT LOOP: WRITE(6,*)'OUT of UT loop' write(43,5145)tpkt write(43,5145)upkt write(43,5145)vpkt write(43,5145)wpkt write(43,5145)utpkt write(43,5145)alat write(43,5145)alon write(43,5145)alata write(43,5145)alona write(43,5145)zpkt write(43,5145)z3pkt write(43,5145)qpkt write(43,5145)alampkt write(43,5145)uxpkt write(43,5145)vxpkt WRITE(43,5145)TEPKT WRITE(43,5145)COMPKT WRITE(43,5145)COMPKT2 WRITE(43,5145)DN2PPKT if(ismin .eq. 1)then write(43,5145)psn4pk write(43,5145)psnopk write(43,5145)psn2dpk write(43,5145)an4pkt write(43,5145)bn4pkt write(43,5145)anopkt write(43,5145)bnopkt write(43,5145)an2dpkt write(43,5145)bn2dpkt endif WRITE(6,*)'MANaged to write files?' C C If ion trajectories are used then write ion terms to a file IF(JEOMGN .EQ. 1)THEN DO 8540 KK=1,NTJ,1 WRITE(41,5145)(UTPKT(JUT),JUT=1,NTM,1) WRITE(41,5145)(ALATA(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(ALONA(JUT,KK),JUT=1,NTM,1) KLIP = 0 IF(KLIP .EQ. 1)THEN WRITE(41,5145)(ALAT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(ALON(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(TPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(TIPKT(JUT,KK),JUT=1,NTM,1) C WRITE(41,5145)(QIPKT(JUT,1,KK),JUT=1,NTM,1) WRITE(41,5145)(ENDENS(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(VIVN(JUT,2,KK),JUT=1,NTM,1) WRITE(41,5145)(TEPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(QPKT2(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(UPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(VPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(UIPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(VIPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(XNPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(XOPKT(JUT,KK),JUT=1,NTM,1) IF(JION .EQ. 1)THEN WRITE(41,5145)(DOPPKT(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)((AOPPKT(JUT,LUG,KK),JUT=1,NTM,1),LUG=1,10,1) ENDIF C DO 8541 JUT=1,NTM,1 C WRITE(41,5145)(QI3PKT(JUT,KK,K),K = 1,KMX,1) C WRITE(41,5145)(QPPHT(JUT,KK,K),K = 1,KMX,1) 8541 CONTINUE DO 8542 JUT =1,NTM,1 WRITE(41,5145)(Z2PKT(JUT,KK,K),K = 1,KMX,1) 8542 CONTINUE IX3 = 6 C IF(MTGCM .EQ. 4)IX3 = 9 C WRITE(41,5145)(QPKT(JUT,9,KK),JUT=1,NTM,1) ENDIF 8540 CONTINUE ENDIF DO 8640 KK=1,NTJ,1 WRITE(41,5145)(UTPKT(JUT),JUT=1,NTM,1) WRITE(41,5145)(ALATA(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(ALONA(JUT,KK),JUT=1,NTM,1) WRITE(41,5145)(ZPKT(JUT,KK),JUT=1,NTM,1) C WRITE(41,5145)(TPKT(JUT,KK),JUT=1,NTM,1) C WRITE(41,5145)(XNPKT(JUT,KK),JUT=1,NTM,1) C WRITE(41,5145)((COMPKT(JUT,LUG,KK),JUT=1,NTM,1),LUG=1,6,1) 8640 CONTINUE 5145 FORMAT(5E11.4) C CALL PLOTPART STOP END