€€ @€ ªCDK ASAS       €  €€       COMMON/ASAS/ F(1,KMXP55)          €  @€ªCDK ASASF      €  €€      COMMON/ASAS/ F(1,KMXP55)          €  @€ªCDK ASAS2      €  €€      COMMON/ASAS2/ F(IMXP,KMXP55)              €€ @€ªCDK ASAS2F     €€ €€      COMMON/ASAS2/ FF(IMXP,KMXP55)             €€ @€ªCDK CH4H2COF   €$  €€      COMMON/CH4H2COF/RJCH4(KMXP),FK9OP(KMXP),SCO2I(KMXP),RSKIONS(KMXP),€  À€     1RJCO2T(KMXP),RKCO2P(KMXP),SCOI(KMXP)              €€ @€ªCDK CH4H2CFF   €€ €€      COMMON/CH4H2COF/RJCH4(KMXP,7)             €€ @€ªCDK CH4H2CO2   €$  €€      COMMON/CH4H2CO2/RJCH4(KMXP,IMXP),FK9OP(KMXP,IMXP),SCO2I(KMXP,IMXP)€   À€     1  ,RSKIONS(KMXP,IMXP),RJCO2T(KMXP,IMXP),RKCO2P(KMXP,IMXP),€€ €     2  SCOI(KMXP,IMXP)         €€ @€ªCDK CH4H2C2F   €€ €€      COMMON/CH4H2CO2/RJCH4F(KMXP,IMXP,7)               €€ @€ªCDK CONS       €#  €€      COMMON/CONS/C(120),RMASS(3),T0(KMXP),LEN1,LEN2,LEN3,KMAX,KMAXP1,  €€ À€     1EXPS(KMXP),DIFK(KMXP)             €  @€	ªCDK CONS2      €#  €€	      COMMON/CONS/C(120),RMASS(3),T0(KMXP),LEN1,LEN2,LEN3,KMAX,KMAXP1,  €  À€	     1EXPS(KMXP),DIFK(IMXP,KMXP),IMAX,IMAXP2,IMAXP4,FB(IMXP,2)          €  @€
ªCDK CONS3      €#€ €€
      COMMON/CONS/CC(120),RMASS(3),T0(KMXP),LEN1,LEN2,LEN3,KMAX,KMAXP1, €  À€
     1EXPS(KMXP),DIFK(IMXP,KMXP),IMAX,IMAXP2,IMAXP4,FB(IMXP,2)          €€ @€ªCDK CRATES     €€ €€      COMMON/CRATES/RK12(KMXP),RK13(KMXP),RK14(KMXP),RK15(KMXP)         €  @€ªCDK CRATESF    €  €€      COMMON/CRATES/RK12(KMXP,4)        €  @€ªCDK CRATES2    €"€ €€      COMMON/CRATES2/RK12(KMXP,IMXP),RK13(KMXP,IMXP),RK14(KMXP,IMXP),   €€ À€     1  RK15(KMXP,IMXP)         €€ @€ªCDK CRATES2F   €€ €€      COMMON/CRATES2/RK12F(KMXP,IMXP,4)         €  @€ªCDK DISSHOX    €!€ €€      COMMON/DISSHOX/SH2OT(KMXP),YIOP(KMXP),YIHP(KMXP),YIOHP(KMXP),     €  À€     1SH2OSRB(KMXP),SH2OLYA(KMXP),SH2OSRC(KMXP),SH2OEUV(KMXP),  €€ €     2SHOXI(KMXP)               €€ @€ªCDK DISSHOXF   €  €€      COMMON/DISSHOX/SH2OT(KMXP,9)              €€ @€ªCDK DISSHOX2   €#€ €€      COMMON/DISSHOX2/SH2OT(KMXP,IMXP),YIOP(KMXP,IMXP),YIHP(KMXP,IMXP), €€ À€     1  YIOHP(KMXP,IMXP),SH2OSRB(KMXP,IMXP),SH2OLYA(KMXP,IMXP), €  €     2  SH2OSRC(KMXP,IMXP),SH2OEUV(KMXP,IMXP),SHOXI(KMXP,IMXP)          €€ @€ªCDK DISSHX2F   €€ €€      COMMON/DISSHOX2/SH2OTF(KMXP,IMXP,9)               €€ @€ªCDK FIELDS     €$  €€      COMMON/FIELDS/FT(KMXP),FNO2(KMXP),FNO(KMXP),FNN2(KMXP),FNAR(KMXP),€#€ À€     1FNHE(KMXP),FNNO(KMXP),FNN2D(KMXP),FNN4S(KMXP),FRJ(KMXP),FW(KMXP), € € €     2FNHOX(KMXP),PPN2D(KMXP),FNE(KMXP),PPN4S(KMXP),FK4O2P(KMXP),       € @€     3FK5O2P(KMXP),F107,FTE(KMXP),FDIFK(KMXP),FU(KMXP)          €  @€ªCDK FIELDSF    €  €€      COMMON/FIELDS/FT(KMXP,17),F107,FTE(KMXP,3)        €  @€ªCDK FIELDS2    €!  €€      COMMON/FIELDS2/FT(KMXP,IMXP),FNO2(KMXP,IMXP),FNO(KMXP,IMXP),      €#  À€     1FNN2(KMXP,IMXP),FNAR(KMXP,IMXP),FNHE(KMXP,IMXP),FNNO(KMXP,IMXP),  €"€ €     2FNN2D(KMXP,IMXP),FNN4S(KMXP,IMXP),FRJ(KMXP,IMXP),FW(KMXP,IMXP),   €$ @€     3FNHOX(KMXP,IMXP),PPN2D(KMXP,IMXP),FNE(KMXP,IMXP),PPN4S(KMXP,IMXP),€ €€     4FK4O2P(KMXP,IMXP),FK5O2P(KMXP,IMXP),F107,FTE(KMXP,IMXP),  €! À€     5FDIFK(KMXP,IMXP),FU(KMXP,IMXP),FWI(KMXP,IMXP),FUI(KMXP,IMXP)      €$  €     6,FNAS(KMXP,IMXP),FNAO(KMXP,IMXP),FNAO2(KMXP,IMXP),FNAOH(KMXP,IMXP)        €€ @€ªCDK FIELDS2F   €  €€      COMMON/FIELDS2/FTF(KMXP,IMXP,17),F107F,FTEF(KMXP,IMXP,9)          €€ @€ªCDK FIELDT     €"  €€      COMMON/FIELDT/FNH2O(KMXP),FNH2O2(KMXP),FNO1D(KMXP),FNH2(KMXP),    €$  À€     1FNOH(KMXP),FNHO2(KMXP),FNO3(KMXP),FNH(KMXP),FNCH4(KMXP),FNCO(KMXP)€  €     2,FNCO2(KMXP),FNNO2(KMXP)          €  @€ªCDK FIELDTF    €  €€      COMMON/FIELDT/FNH2O(KMXP,12)              €  @€ªCDK FIELDT2    €$  €€      COMMON/FIELDT2/FNH2O(KMXP,IMXP),FNH2O2(KMXP,IMXP),FNO1D(KMXP,IMXP)€$  À€     1,FNH2(KMXP,IMXP),FNOH(KMXP,IMXP),FNHO2(KMXP,IMXP),FNO3(KMXP,IMXP),€#€ €     2FNH(KMXP,IMXP),FNCH4(KMXP,IMXP),FNCO(KMXP,IMXP),FNCO2(KMXP,IMXP), € @€     3FNNO2(KMXP,IMXP)          €€ @€ªCDK FIELDT2F   €€ €€      COMMON/FIELDT2/FNH2OF(KMXP,IMXP,12)               €€ @€ªCDK FLUX       €€ €€      COMMON/FLUX/FLUXNO,IBNDNO,FLUXCO,IBNDCO           €  @€ªCDK FLUX2      €  €€      COMMON/FLUX2/FLUXNO(IMXP),IBNDNO,FLUXCO(IMXP),IBNDCO              €€ @€ªCDK FLUX2F     €  €€      COMMON/FLUX2/FLUXNOF(IMXP),IBNDNOF,FLUXCOF(IMXP),IBNDCOF          €€ @€ªCDK H2FORG     €€ €€      COMMON/H2FORG/RJH2OLY(KMXP)               €  @€ªCDK H2FORG2    €€ €€      COMMON/H2FORG2/RJH2OLY(KMXP,IMXP)         €€ @€ ªCDK H2FORG2F   €  €€       COMMON/H2FORG2/RJH2OLYF(KMXP,IMXP)        €€ @€!ªCDK HEFLUX     €  €€!      COMMON/HEFLUX/FLUX        €  @€"ªCDK HEFLUX2    €€ €€"      COMMON/HEFLUX2/FLUX(IMXP)         €€ @€#ªCDK HEFLUX2F   €  €€#      COMMON/HEFLUX2/FLUXF(IMXP)        €  @€$ªCDK HOXSTUF    €!  €€$      COMMON/HOXSTUF/HO2HR(KMXP),OHHR(KMXP),HHOXR(KMXP),ANT(KMXP),      €"€ À€$     1BNT(KMXP),CNT(KMXP),DNT(KMXP),ASNT(KMXP),CONPR(KMXP),DSNT(KMXP)           €€ @€%ªCDK HOXSTUFF   €€ €€%      COMMON/HOXSTUF/HO2HR(KMXP,10)             €€ @€&ªCDK HOXSTUF2   €$  €€&      COMMON/HOXSTUF2/HO2HR(KMXP,IMXP),OHHR(KMXP,IMXP),HHOXR(KMXP,IMXP),€"  À€&     1  ANT(KMXP,IMXP),BNT(KMXP,IMXP),CNT(KMXP,IMXP),DNT(KMXP,IMXP),    €  €&     2  ASNT(KMXP,IMXP),CONPR(KMXP,IMXP),DSNT(KMXP,IMXP)        €€ @€'ªCDK HOXSTF2F   €  €€'      COMMON/HOXSTUF2/HO2HRF(KMXP,IMXP,10)              €€ @€(ªCDK HOXUPP     €  €€(      COMMON/HOXUPP/FHOX        €  @€)ªCDK HOXUPP2    €€ €€)      COMMON/HOXUPP2/FHOX(IMXP)         €€ @€*ªCDK HOXUPP2F   €  €€*      COMMON/HOXUPP2/FHOXF(IMXP)        €  @€+ªCDK INDEX      €"€ €€+      COMMON/INDEX/NJ,NT,NU,NPS,NPS2,NPSA,NPSH,NPNO,NPN4S,NPCH4,NPH2,   €!  À€+     1NPCO2,NPCO,NPHOX,NPH2O,NW,NPSNM,NPS2NM,NPSANM,NPSHNM,NPNONM,      €!€ €+     2NN4SNM,NCH4NM,NPH2NM,NCO2NM,NPCONM,NHOXNM,NH2ONM,NJNP,NDJ,NRJ             €€ @€,ªCDK INDEXDAT   €!€ €€,      DATA NJ,NT,NU,NPS,NPS2,NPSA,NPSH,NPNO,NPN4S,NPCH4,NPH2,NPCO2,     €"€ À€,     1  NPCO,NPHOX,NPH2O,NW,NPSNM,NPS2NM,NPSANM,NPSHNM,NPNONM,NN4SNM,   €$  €,     2  NCH4NM,NPH2NM,NCO2NM,NPCONM,NHOXNM,NH2ONM,NJNP,NDJ,NRJ/1,0,KMXP,€# @€,     3  KMXP2,KMXP3,KMXP4,KMXP5,KMXP6,KMXP7,KMXP8,KMXP9,KMXP10,KMXP11,  €#€€€,     4  KMXP12,KMXP13,KMXP14,KMXP15,KMXP16,KMXP17,KMXP18,KMXP19,KMXP20, €! À€,     5  KMXP21,KMXP22,KMXP23,KMXP24,KMXP25,KMXP26,KMP27P,KMP54P,0/              €€ @€-ªCDK IONHOX     €  €€-      COMMON/IONHOX/PHOXIC(KMXP)        €  @€.ªCDK IONHOX2    €  €€.      COMMON/IONHOX2/PHOXIC(KMXP,IMXP)          €€ @€/ªCDK IONHOX2F   €€ €€/      COMMON/IONHOX2/PHOXICF(KMXP,IMXP)         €€ @€0ªCDK NEWRAT     €€ €€0      COMMON/NEWRAT/HO2OH(KMXP),HOH(KMXP),OHHOX(KMXP)           €  @€1ªCDK NEWRATF    €  €€1      COMMON/NEWRATF/HO2OH(KMXP,3)              €  @€2ªCDK NEWRAT2    €"€ €€2      COMMON/NEWRAT2/HO2OH(KMXP,IMXP),HOH(KMXP,IMXP),OHHOX(KMXP,IMXP)           €€ @€3ªCDK NEWRAT2F   €  €€3      COMMON/NEWRAT2/HO2OHF(KMXP,IMXP,3)        €€ @€4ªCDK NOZNOZ     €€ €€4      COMMON/NOZNOZ/SNO2NO(KMXP),RNONOZ(KMXP),FNNOZ(KMXP)               €  @€5ªCDK NOZNOZF    €€ €€5      COMMON/NOZNOZF/SNO2NO(KMXP,3)             €  @€6ªCDK NOZNOZ2    €€ €€6      COMMON/NOZNOZ2/SNO2NO(KMXP,IMXP),RNONOZ(KMXP,IMXP),       €  À€6     1  FNNOZ(KMXP,IMXP)        €€ @€7ªCDK NOZNOZ2F   €€ €€7      COMMON/NOZNOZ2/SNO2NOF(KMXP,IMXP,3)               €  @€8ªCDK OPION      €€ €€8      COMMON/OPION/FIOP(KMXP)           €€ @€9ªCDK OPION2     €€ €€9      COMMON/OPION2/FIOP(KMXP,IMXP)             €  @€:ªCDK OPION2F    €  €€:      COMMON/OPION2/FIOPF(KMXP,IMXP)            €€ @€;ªCDK OXOXOX     €"€ €€;      COMMON/OXOXOX/O3OR(KMXP),XNOX(KMXP),OOXR(KMXP),FNOX(KMXP),XNOX1   €€ À€;     1,PS2B             €  @€<ªCDK OXOXOXF    €€ €€<      COMMON/OXOXOX/O3OR(KMXP,4),XNOX1(2)               €  @€=ªCDK OXOXOX2    €   €€=      COMMON/OXOXOX2/O3OR(KMXP,IMXP),XNOX(KMXP),OOXR(KMXP,IMXP),€  À€=     1  FNOX(KMXP,IMXP),XNOX1,PS2B(IMXP)        €€ @€>ªCDK OXOXOX2F   €"€ €€>      COMMON/OXOXOX2/O3ORF(KMXP,IMXP),XNOXF(KMXP),OOXRF(KMXP,IMXP,2),   €  À€>     1  XNOX1,PS2B(IMXP)        €€ @€?ªCDK PARAMS     €€ €€?      PARAMETER (KMX=95,IMX=20) €€ À€?      PARAMETER (KMXP=KMX+1,IMXP=IMX+4) €  €?      PARAMETER (KMXP2=2*KMXP,KMXP3=3*KMXP,KMXP4=4*KMXP,€$ @€?     1  KMXP5=5*KMXP,KMXP6=6*KMXP,KMXP7=7*KMXP,KMXP8=8*KMXP,KMXP9=9*KMXP€"€€€?     2  ,KMXP10=10*KMXP,KMXP11=11*KMXP,KMXP12=12*KMXP,KMXP13=13*KMXP,   €" À€?     3  KMXP14=14*KMXP,KMXP15=15*KMXP,KMXP16=16*KMXP,KMXP17=17*KMXP,    €"  €?     4  KMXP18=18*KMXP,KMXP19=19*KMXP,KMXP20=20*KMXP,KMXP21=21*KMXP,    €" @€?     5  KMXP22=22*KMXP,KMXP23=23*KMXP,KMXP24=24*KMXP,KMXP25=25*KMXP,    €$ €€?     6  KMXP26=26*KMXP,KMP27P=27*KMXP+1,KMP54P=54*KMXP+1,KMXP55=55*KMXP)        €€ @€@ªCDK PSIB       €#€ €€@      COMMON/PSIB/PARB,PHB,PNOB,PN4SB,PCH4B,PH2B,PCO2B,PCOB,PHOXB,PH2OB         €  @€AªCDK PSIBF      €  €€A      COMMON/PSIB/PARB(10)              €  @€BªCDK PSIB2      €€ €€B      COMMON/PSIB2/PARB(IMXP),PHB(IMXP),PNOB(IMXP),PN4SB(IMXP), €!  À€B     1  PCH4B(IMXP),PH2B(IMXP),PCO2B(IMXP),PCOB(IMXP),PHOXB(IMXP),      €	€ €B     2  PH2OB(IMXP)             €€ @€CªCDK PSIB2F     €€ €€C      COMMON/PSIB2/PARBF(IMXP,10)               €  @€DªCDK PHOTO      €  €€D      COMMON/PHOTO/XNO,XN4S,XCH4,XH2,XCO2,XCO,XH2O              €€ @€EªCDK PHOTOF     €€ €€E      COMMON/PHOTO/XNO(7)               €€ @€FªCDK PHOTO2     €$  €€F      COMMON/PHOTO2/XNO(IMXP),XN4S(IMXP),XCH4(IMXP),XH2(IMXP),XCO2(IMXP)€€ À€F     1  ,XCO(IMXP),XH2O(IMXP)           €  @€GªCDK PHOTO2F    €  €€G      COMMON/PHOTO2/XNOF(IMXP,7)        €€ @€HªCDK PHOTOO     €  €€H      COMMON/PHOTOO/XHE,XARG            €  @€IªCDK PHOTOOF    €  €€I      COMMON/PHOTOO/XHE(2)              €  @€JªCDK PHOTOO2    €€ €€J      COMMON/PHOTOO2/XHE(IMXP),XARG(IMXP)               €€ @€KªCDK PHOTOO2F   €€ €€K      COMMON/PHOTOO2/XHEF(IMXP,2)               €  @€LªCDK RATEBLK    €€ €€L      COMMON/RATEBLK/RKK(25,KMXP),ALPP(10,KMXP),RBB(20,KMXP),   €	€ À€L     1RKMM(50,KMXP)             €€ @€MªCDK RATEBLK2   €€ €€M      COMMON/RATEBLK2/RKK(25,KMXP,IMXP),ALPP(10,KMXP,IMXP),     €  À€M     1RBB(20,KMXP,IMXP),RKMM(50,KMXP,IMXP)              €€ @€NªCDK RATBLK2F   €€ €€N      COMMON/RATEBLK2/RKKF(25,KMXP,IMXP),ALPPF(10,KMXP,IMXP),   €  À€N     1RBBF(20,KMXP,IMXP),RKMMF(50,KMXP,IMXP)            €  @€OªCDK RJH2O2O    €€ €€O      COMMON/RJH2O2O/RJH2O2(KMXP)               €€ @€PªCDK RJH2O2O2   €€ €€P      COMMON/RJH2O2O2/RJH2O2(KMXP,IMXP)         €€ @€QªCDK RJH2O22F   €  €€Q      COMMON/RJH2O2O2/RJH2O2F(KMXP,IMXP)        €€ @€RªCDK SOURCE     €   €€R      COMMON/SOURCE/FS11(KMXP),FS12(KMXP),FS21(KMXP),FS22(KMXP),€€ À€R     1FS1(KMXP),FS2(KMXP)               €  @€SªCDK SOURCEF    €  €€S      COMMON/SOURCE/FS11(KMXP,6)        €  @€TªCDK SOURCE2    €"€ €€T      COMMON/SOURCE2/FS11(KMXP,IMXP),FS12(KMXP,IMXP),FS21(KMXP,IMXP),   €€ À€T     1FS22(KMXP,IMXP),FS1(KMXP,IMXP),FS2(KMXP,IMXP)             €€ @€UªCDK SOURCE2F   €€ €€U      COMMON/SOURCE2/FS11F(KMXP,IMXP,6)         €€ @€VªCDK VSCR       €€ €€V      COMMON/VSCR/S15(IMXP,KMXP),S14(IMXP,KMXP),S13(IMXP,KMXP), € € À€V     1S12(IMXP,KMXP),S11(IMXP,KMXP),S10(IMXP,KMXP),S9(IMXP,KMXP),       €  €V     2S8(IMXP,KMXP),S7(IMXP,KMXP),S6(IMXP,KMXP),S5(IMXP,KMXP),  €# @€V     3S4(IMXP,KMXP),S3(IMXP,KMXP),S2(IMXP,KMXP),S1(IMXP,KMXP),T1(IMXP)  € €€V     4,T2(IMXP),T3(IMXP),T4(IMXP),T5(IMXP),T6(IMXP),T7(IMXP)            €  @€WªCDK VSCR2      €   €€W      COMMON/VSCR/GAMA(IMXP,KMXP,2,2),Z(IMXP,KMXP,2),EMBAR(IMXP,€€ À€W     1KMXP),AK(IMXP,2,2,2),EP(IMXP,2,2),FK(IMXP,2),PK(IMXP,     €€ €W     22,2),QK(IMXP,2,2), RK(IMXP,2,2),WKM1(IMXP,2,2),WKM2(IMXP, €# @€W     32,2),WKV1(IMXP,2),WKV2(IMXP,2),WKS1(IMXP),WKS2(IMXP),WKS3(IMXP),  € €€W     4EMBAR0(IMXP),PS0(IMXP,2)          €  @€XªCDK VSCR3      €€ €€X      COMMON/VSCR/WK1(KMXP,IMXP),WK2(KMXP,IMXP),WK3(KMXP,IMXP), €"  À€X     1  WK4(KMXP,IMXP),WK5(KMXP,IMXP),WK6(KMXP,IMXP),WK7(KMXP,IMXP),    €#  €X     2  WK8(KMXP,IMXP),WK9(KMXP,IMXP),WK10(KMXP,IMXP),WK11(KMXP,IMXP),  €€@€X     3  WK12(KMXP,IMXP),WK13(KMXP,IMXP)         €€ @€YªDK BLKDATA     €  €€YCDIR$ NOLIST    €  À€Y       BLOCK DATA BLOCKS€  €YªCA PARAMS      €€@€Y      COMMON/BLK1/ZERO1 €€€€YªCA ASAS2       €€À€Y      COMMON/BLK2/ZERO2 €  €YªCA CH4H2CO2    €€@€Y      COMMON/BLK3/ZERO3 €€€€YªCA CONS2       €€À€Y      COMMON/BLK4/ZERO4 €€ €YªCA CRATES2     €€@€Y      COMMON/BLK5/ZERO5 € €€YªCA DISSHOX2    €€À€Y      COMMON/BLK6/ZERO6 €€ €YªCA FIELDS2     €€@€Y      COMMON/BLK7/ZERO7 €€€€YªCA FIELDT2     €€À€Y      COMMON/BLK8/ZERO8 €€ €YªCA FLUX2       €€@€Y      COMMON/BLK9/ZERO9 €€€€YªCA H2FORG2     €€À€Y      COMMON/BLK10/ZERO10       €€ €YªCA HEFLUX2     €€@€Y      COMMON/BLK11/ZERO11       € €€YªCA HOXSTUF2    €€À€Y      COMMON/BLK12/ZERO12       €€ €YªCA HOXUPP2     €€@€Y      COMMON/BLK13/ZERO13       €€€€YªCA INDEX       €€À€Y      COMMON/BLK14/ZERO14       €€ €YªCA IONHOX2     €€@€Y      COMMON/BLK15/ZERO15       €€€€YªCA NEWRAT2     €€À€Y      COMMON/BLK16/ZERO16       €€	 €YªCA NOZNOZ2     €€	@€Y      COMMON/BLK17/ZERO17       € 	€€YªCA OPION2      €€	À€Y      COMMON/BLK18/ZERO18       €€
 €YªCA OXOXOX2     €€
@€Y      COMMON/BLK19/ZERO19       €€
€€YªCA PSIB2       €€
À€Y      COMMON/BLK20/ZERO20       €  €YªCA PHOTO2      €€@€Y      COMMON/BLK21/ZERO21       €€€€YªCA PHOTOO2     €€À€Y      COMMON/BLK22/ZERO22       €  €YªCA RATEBLK2    €€@€Y      COMMON/BLK23/ZERO23       € €€YªCA RJH2O2O2    €€À€Y      COMMON/BLK24/ZERO24       €€ €YªCA SOURCE2     €€@€Y      COMMON/BLK25/ZERO25       € €€YªCA VSCR€€À€Y      COMMON/BLK26/ZERO26       €   €Y      DATA ZERO1,ZERO2,ZERO3,ZERO4,ZERO5,ZERO6,ZERO7,ZERO8,ZERO9€ @€Y     1  /1.,2.,3.,4.,5.,6.,7.,8.,9./    €!€€€Y      DATA ZERO10,ZERO11,ZERO12,ZERO13,ZERO14,ZERO15,ZERO16,ZERO17,     € À€Y     1  ZERO18,ZERO19/10.,11.,12.,13.,14.,15.,16.,17.,18.,19./  €€ €Y      DATA ZERO20,ZERO21,ZERO22,ZERO23,ZERO24,ZERO25,ZERO26     €€@€Y     2  /20.,21.,22.,23.,24.,25.,26./   €€€€Y      END       € €À€YC       € € €YC       € @€Y      SUBROUTINE CHKBLKS€ €€Y      SAVE      €  À€Y€€ €Y      COMMON/BLK1/ZERO1 €€@€Y      COMMON/BLK2/ZERO2 €€€€Y      COMMON/BLK3/ZERO3 €€À€Y      COMMON/BLK4/ZERO4 €€ €Y      COMMON/BLK5/ZERO5 €€@€Y      COMMON/BLK6/ZERO6 €€€€Y      COMMON/BLK7/ZERO7 €€À€Y      COMMON/BLK8/ZERO8 €€ €Y      COMMON/BLK9/ZERO9 €€@€Y      COMMON/BLK10/ZERO10       €€€€Y      COMMON/BLK11/ZERO11       €€À€Y      COMMON/BLK12/ZERO12       €€ €Y      COMMON/BLK13/ZERO13       €€@€Y      COMMON/BLK14/ZERO14       €€€€Y      COMMON/BLK15/ZERO15       €€À€Y      COMMON/BLK16/ZERO16       €€ €Y      COMMON/BLK17/ZERO17       €€@€Y      COMMON/BLK18/ZERO18       €€€€Y      COMMON/BLK19/ZERO19       €€À€Y      COMMON/BLK20/ZERO20       €€ €Y      COMMON/BLK21/ZERO21       €€@€Y      COMMON/BLK22/ZERO22       €€€€Y      COMMON/BLK23/ZERO23       €€À€Y      COMMON/BLK24/ZERO24       €€ €Y      COMMON/BLK25/ZERO25       €€@€Y      COMMON/BLK26/ZERO26       €$ €€YC     WRITE(6,100)ZERO1,ZERO2,ZERO3,ZERO4,ZERO5,ZERO6,ZERO7,ZERO8,ZERO9,€  À€YC    1  ZERO10,ZERO11,ZERO12,ZERO13,ZERO14,ZERO15,ZERO16,ZERO17,€#  €YC    2  ZERO18,ZERO19,ZERO20,ZERO21,ZERO22,ZERO23,ZERO24,ZERO25,ZERO26  €
€@€YC 100 FORMAT(10E12.4)   € €€Y      RETURN    €€À€Y      END               €  @€ZªDK CMPN4S      €€ €€Z      SUBROUTINE CMPN4S €  À€ZCDIR$ BOUNDS    €  €Z      SAVE      € @€ZC     ****      €€€€ZC     ****     ADVANCE N4S COMPOSITION BY ONE TIME STEP € À€ZC     ****      €  €ZC     ****     COMMON DECKS:    €$ @€ZC     ****       PARAMS, ASAS, CONS, FIELDS, FIELDT, FLUX, INDEX, PHOTO,€€€€ZC     ****          PSIB, RATEBLK, VSCR € À€ZC     ****      €  €ZªCA PARAMS      €€@€ZªCA ASAS2       €€€€ZªCA CONS2       €€À€ZªCA FIELDS2     €€ €ZªCA FIELDT2     €€@€ZªCA FLUX2       €€€€ZªCA INDEX       € À€ZªCA PHOTO2      €€ €ZªCA PSIB2       € @€ZªCA RATEBLK2    € €€ZªCA VSCR€€À€Z      DIMENSION PHIN4S(3)       €  €Z      DATA RMN4S,PHIN4S/14.,0.651,0.731,0.741/  € @€ZC     ****      € €€ZC     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM    €€À€ZC     ****     LOWER BOUNDARY: PHOTO-CHEMICAL EQUILIBRIUM       €  €ZC     ****      €
€@€Z      DO 1 I=3,LEN1-2   €€€€Z	T4(I)=0.       €€À€Z	T1(I)=0.       €€ €Z	T2(I)=1.       € @€Z	T3(I)=-XN4S(I)*RMN4S/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+    €€€€Z     1    FNN2(1,I)*RMASS(3))   € À€Z    1   CONTINUE€ 	 €ZC     ****      € 	@€ZC     ****     SOURCES  € 	€€ZC     ****      €
€	À€Z      DO 2 K=1,KMAXP1   € 
 €Z	DO 2 I=3,LEN1-2€€
@€Z	  S1(I,K)=-(FK4O2P(K,I)+RBB(1,K,I)*FNO2(K,I)+RBB(3,K,I)*       € 
€€Z     1      FNNO(K,I)+RBB(10,K,I)*FNOH(K,I))    €  
À€Z	  S2(I,K)=PPN4S(K,I)+RBB(4,K,I)*FNN2D(K,I)*FNO(K,I)+RBB(5,K,I)*€ € €Z     1      FNN2D(K,I)*FNE(K,I)+RBB(7,K,I)*FNN2D(K,I)+RBB(8,K,I)*       €
€@€Z     2      FNNO(K,I)   € €€Z    2 CONTINUE  € À€ZC     ****      €  €ZC     ****     PERIODIC POINTS  € @€ZC     ****      €	 €€Z      DO 3 I = 1,2      €	€À€Z	T1(I) = T1(I+IMAX)     €€ €Z	T1(I+IMAXP2) = T1(I+2) €	€@€Z	T2(I) = T2(I+IMAX)     €€€€Z	T2(I+IMAXP2) = T2(I+2) €	€À€Z	T3(I) = T3(I+IMAX)     €€ €Z	T3(I+IMAXP2) = T3(I+2) €	€@€Z	T4(I) = T4(I+IMAX)     €€€€Z	T4(I+IMAXP2) = T4(I+2) €
 À€Z	  DO 3 K = 1,KMAXP1    €€ €Z	    S1(I,K) = S1(I+IMAX,K)     €€@€Z	    S1(I+IMAXP2,K) = S1(I+2,K) €€€€Z	    S2(I,K) = S2(I+IMAX,K)     €€À€Z	    S2(I+IMAXP2,K) = S2(I+2,K) €  €Z    3 CONTINUE  € @€Z      IBND=0    €€€€Z      IBNDB=0   €€À€Z      ALFA=0.   € € €Z      CALL MINOR(NPN4S,NN4SNM,RMN4S,PHIN4S,ALFA,IBND,IBNDB,PN4SB)       € @€Z      RETURN    €€€€Z      END               €  @€[ªDK PPHOTO      €  €€[CDIR$ NOLIST    €  À€[      SUBROUTINE PHOTO  €  €[      SAVE      € @€[C     ****      €€€€[C     ****     COMMON DECKS NEEDED:     €# À€[C     ****       PARAMETER, CONS, FIELDS, FIELDT, FLUX, PHOTO, PHOTOO,  €  €[C     ****         RATEBLK      € @€[C     ****      € €€[ªCA PARAMS      €€À€[ªCA CONS3       €€ €[ªCA FIELDS2     €€@€[ªCA FIELDT2     €€€€[ªCA FLUX2       € À€[ªCA PHOTO2      €€ €[ªCA PHOTOO2     € @€[ªCA RATEBLK2    €$ €€[      DIMENSION A(IMXP),B(IMXP),C(IMXP),D(IMXP),E(IMXP),R(IMXP),S(IMXP),€ À€[     1  T(IMXP),FM(IMXP)€€ €[      DO 1 I = 3,LEN1-2 €€@€[	A(I) = PPN4S(1,I)+RBB(4,1,I)*FNN2D(1,I)*FNO(1,I)+RBB(5,1,I)*   €€€€[     1    FNN2D(1,I)*FNE(1,I)+RBB(7,1,I)*FNN2D(1,I)     €€À€[	B(I)=RBB(2,1,I)*FNO2(1,I)*FNN2D(1,I)+RBB(13,1,I)*FNO(1,I)*     €  €[     1    FNNO2(1,I)+RBB(14,1,I)**FNNO2(1,I)    € @€[	C(I)=FK4O2P(1,I)+RBB(1,1,I)*FNO2(1,I)+RBB(10,1,I)*FNOH(1,I)    € €€[	D(I)=FK5O2P(1,I)+RBB(6,1,I)*FNN2D(1,I)+RBB(8,1,I)+RBB(9,1,I)+  € À€[     1    RBB(11,1,I)*FNO3(1,I)+RBB(12,1,I)*FNHO2(1,I)  €  €[	E(I)=RBB(1,1,I)*FNO2(1,I)+RBB(10,1,I)*FNOH(1,I)€€@€[	R(I)=RBB(3,1,I)*E(I)+RBB(3,1,I)*C(I)   €€€€[	S(I)=-(RBB(8,1,I)*E(I)-RBB(3,1,I)*B(I)-D(I)*C(I)+RBB(3,1,I)*   €€À€[     1    A(I)) €  €[	T(I)=-(RBB(8,1,I)*B(I)+D(I)*A(I))      €€@€[	XN4S(I)=(-S(I)+SQRT(S(I)*S(I)-4.*R(I)*T(I)))/(2.*R(I)) €€€€[	XNO(I)=(C(I)*XN4S(I)-A(I))/(RBB(8,1,I)-RBB(3,1,I)*XN4S(I))     €€À€[	FM(I)=FNO(1,I)+FNO2(1,I)+FNN2(1,I)     €
 	 €[	XNO(I)=1.1E-8*FM(I)    € 	@€[	XCH4(I)=0.16E-6*FM(I)  €
 	€€[	XH2(I)=0.5E-6*FM(I)    € 	À€[	XCO2(I)=350.E-6*FM(I)  €
 
 €[	XCO(I)=3.6E-8*FM(I)    €
€
@€[	XH2O(I)=6.0E-6*FM(I)   €
€
€€[	XHE(I)=5.24E-6*FM(I)   € 
À€[	XARG(I)=9.34E-3*FM(I)  €  €[    1 CONTINUE  € @€[C     WRITE(6,1234)(XN4S(I),XNO(I),XCH4(I),XH2(I),XCO2(I),      €€€€[C    1  XCO(I),XH2O(I),XHE(I),XARG(I),I=1,LEN1) €€À€[C     *****     €€ €[C     *****     PERIODIC POINTS €€@€[C     *****     €	 €€[      DO 2 I = 1,2      €€À€[	XN4S(I) = XN4S(I+IMAX) €€ €[	XN4S(I+IMAXP2) = XN4S(I+2)     €
€@€[	XNO(I) = XNO(I+IMAX)   €€€€[	XNO(I+IMAXP2) = XNO(I+2)       €€À€[	XCH4(I) = XCH4(I+IMAX) €€ €[	XCH4(I+IMAXP2) = XCH4(I+2)     €
€@€[	XH2(I) = XH2(I+IMAX)   €€€€[	XH2(I+IMAXP2) = XH2(I+2)       €€À€[	XCO2(I) = XCO2(I+IMAX) €€ €[	XCO2(I+IMAXP2) = XCO2(I+2)     € @€[    2 CONTINUE  €€€€[      DO 3 I =1,2       €
€À€[	XCO(I) = XCO(I+IMAX)   €€ €[	XCO(I+IMAXP2) = XCO(I+2)       €€@€[	XH2O(I) = XH2O(I+IMAX) €€€€[	XH2O(I+IMAXP2) = XH2O(I+2)     €
€À€[	XHE(I) = XHE(I+IMAX)   €€ €[	XHE(I+IMAXP2) = XHE(I+2)       €€@€[	XARG(I) = XARG(I+IMAX) €€€€[	XARG(I+IMAXP2) = XARG(I+2)     € À€[    3 CONTINUE  €  €[C     WRITE(6,1234)(XN4S(I),XNO(I),XCH4(I),XH2(I),XCO2(I),      €€@€[C    1  XCO(I),XH2O(I),XHE(I),XARG(I),I=1,LEN1) €
 €€[C1234 FORMAT(9E12.4)    € À€[      RETURN    €€ €[      END       € @€[CDIR$ NOLIST            €  @€\ªDK COMP€  €€\CDIR$ NOLIST    €
€ À€\      SUBROUTINE COMP   €  €\      SAVE      €€@€\C     ****    ADVANCE PSI IN TIME       € €€\C     ****      €€À€\C     ****    COMMON DECKS:     €  €\C     ****      €  @€\C     ****      PARAMS, ASAS, CONS, CRATES, INDEX, SOURCE, VSCR2€ €€\C     ****      € À€\ªCA PARAMS      €€ €\ªCA ASAS2       €€@€\ªCA CONS2       €€€€\ªCA CRATES2     €€À€\ªCA INDEX       €€ €\ªCA SOURCE2     €€@€\ªCA VSCR2       € €€\      DIMENSION AK0(2,2),S(2,2),PHI(2,3),DELTA(2,2),B(2,2)      €€À€\      DIMENSION SS(2,2) €"  €\      DIMENSION ENMBAR(IMXP),FS11ST(IMXP),FS12ST(IMXP),FS21ST(IMXP),    €€@€\     1  FS22ST(IMXP),ST(IMXP)   €€€€\      EQUIVALENCE (C(40),B)     €€À€\      DATA PHI/0.,0.673,1.35,0.,1.11,0.769/,TAU/1.86E+3/,       €€ €\     AS/-1.,1.,0.,0./,DELTA/1.,0.,0.,1./,       €€@€\     BT00/273./ €€€€\      DATA SS/0.,0.,1.,-1./     €€À€\      DATA SMALL/1.E-6/ €  €\      DFACT=1.  € @€\C     ****     CALCULATE EMBAR  €
 €€\      NPS1K=NJ+NPS-1    € À€\      NPS2K=NPS1K+KMAXP1€€ €\      DO 1 K = 1,KMAXP1 €	€@€\      DO 1 I=1,LEN1     €  €€\	EMBAR(I,K)=1./(F(I,NPS1K+K)/RMASS(1)+F(I,NPS2K+K)/RMASS(2)+(1.-€ À€\     A  F(I,NPS1K+K)-F(I,NPS2K+K))/RMASS(3))    € 	 €\    1 CONTINUE  € €	@€\C     CALCULATE PS0, EMBAR0 AND .5*((DMBAR/DZ)/MBAR) AT LEVEL 1/2       €	 	€€\      NPS1K=NJ+NPS      € 	À€\      NPS2K=NPS1K+KMAXP1€
 
 €\      DO 30 I=1,LEN1    €€
@€\	PS0(I,1)=B(1,1)*F(I,NPS1K)+B(1,2)*F(I,NPS2K)+FB(I,1)   €€
€€\	PS0(I,2)=B(2,1)*F(I,NPS1K)+B(2,2)*F(I,NPS2K)+FB(I,2)   €€
À€\	EMBAR0(I)=1./(PS0(I,1)/RMASS(1)+PS0(I,2)/RMASS(2)+(1.-PS0(I,1) €  €\     A  -PS0(I,2))/RMASS(3))    €€@€\	WKS3(I)=(EMBAR(I,1)-EMBAR0(I))/(C(3)*(EMBAR0(I)+EMBAR(I,1)))   € €€\   30 CONTINUE  €  À€\C     ****     CALCULTE EP AND AK AT LEVEL 1/2, SET Z(0) TO ZERO€  €\      KM=1      € @€\      KP=2      € €€\C     ****     EP(1/2)  € À€\      DO 2 L=1,2€  €\	DO 2 I=1,LEN1  €€@€\	  EP(I,L,KP)=1.-(2./(EMBAR0(I)+EMBAR(I,1)))*(RMASS(L)+ € €€\     A    (EMBAR(I,1)-EMBAR0(I))/C(3))  € À€\	  Z(I,1,L)=0.  €  €\    2 CONTINUE  € @€\      DO 3 L=1,2€€€€\	LL=3-L €€À€\	NPSL=NJ+NPS+(L-1)*KMAXP1       €€ €\	DO 3 M=1,2     € @€\	  DO 3 I=1,LEN1€€€€\	    AK(I,L,M,KP)=-DELTA(L,M)*(PHI(LL,3)+(PHI(LL,L)-PHI(LL,3))* € €À€\     A        .5*(PS0(I,L)+F(I,NPSL)))-(1.-DELTA(L,M))*(PHI(L,M)-       €€ €\     B        PHI(L,3))*.5*(PS0(I,L)+F(I,NPSL)) € @€\    3 CONTINUE  €" €€\C     ****     WKS1=MBAR/M3*(T00/(T0+T))*0.25/(TAU*DET) AT LEVEL 1/2    €
 À€\      NTK=NJ+NT+KMAX    €	€ €\      DO 4 I=1,LEN1     €€@€\	WKS1(I)=0.5*(EMBAR0(I)+EMBAR(I,1))/RMASS(3)*(T00/(T0(1)+       €# €€\     A  F(I,NTK)))**0.25/(TAU*(AK(I,1,1,KP)*AK(I,2,2,KP)-AK(I,1,2,KP)*  € À€\     B  AK(I,2,1,KP)))  €  €\    4 CONTINUE  € @€\C     ****     COMPLETE CALCULATION OF AK(1/2)  € €€\      DO 5 L=1,2€€À€\	DO 5 M=1,2     €  €\	  DO 5 I=1,LEN1€ @€\	    AK(I,L,M,KP)=AK(I,L,M,KP)*WKS1(I)  € €€\C           ****     SET GAMA(0) TO 0.  €
€À€\	    GAMA(I,1,L,M)=0.   €  €\    5 CONTINUE  € @€\C     ****     MAIN HEIGHT DO LOOP      €	€€€\      DO 6 K=1,KMAX     €€À€\C       ****     CYCLE K LEVELS €  €\	KB=KM  € @€\	KM=KP  € €€\	KP=KB  € À€\C       **** EP(K+1/2)  €€ €\	DO 7 L=1,2     €€@€\	 DO 7 I=1,LEN1 €€€€\	 EP(I,L,KP)=1.-(2./(EMBAR(I,K)+EMBAR(I,K+1)))*(RMASS(L)+       € À€\     A   (EMBAR(I,K+1)-EMBAR(I,K))/C(3))€  €\    7 CONTINUE  € @€\C     ****     AK(K+1/2)€ €€\      DO 8 L=1,2€€À€\	LL=3-L €€ €\	NPSL=NJ+NPS+(L-1)*KMAXP1+K     €€@€\	DO 8 M=1,2     € €€\	  DO 8 I=1,LEN1€€À€\	    AK(I,L,M,KP)=-DELTA(L,M)*(PHI(LL,3)+(PHI(LL,L)-PHI(LL,3))* €"  €\     A        .5*(F(I,NPSL-1)+F(I,NPSL)))-(1.-DELTA(L,M))*(PHI(L,M)-    € @€\     B        PHI(L,3))*.5*(F(I,NPSL-1)+F(I,NPSL))      € €€\    8 CONTINUE  € À€\C     ****     WKS1=MBAR/M3*(T00/(T0+T))**0.25/(TAU*DET(ALFA))  €€ €\      NTK=NJ+NT+K       €	€@€\      DO 9 I=1,LEN1     € €€€\	WKS1(I)=0.5*(EMBAR(I,K)+EMBAR(I,K+1))/RMASS(3)*(T00/(T0(K+1)+.5*       €$ À€\     A    (F(I,NTK-1)+F(I,NTK))))**0.25/(TAU*(AK(I,1,1,KP)*AK(I,2,2,KP)-€€ €\     B    AK(I,1,2,KP)*AK(I,2,1,KP)))   € @€\C       ****     EDDY DIFFUSION TERMS IN WKS1 AND WKS2  € €€\	WKS2(I)=WKS3(I)€€À€\	WKS3(I)=(EMBAR(I,K+1)-EMBAR(I,K))/(C(3)*(EMBAR(I,K)+   €  €\     C    EMBAR(I,K+1)))€ @€\    9 CONTINUE  € €€\C     ****      €€À€\C     ****     CALCULATE AUGMENTED SOURCE TERMS €  €\C     ****      €	€@€\      NTK=NJ+NT+K-1     € €€\      NPS1K=NJ+NPS+K-1  € À€\      NPS2K=NPS1K+KMAXP1€  €\C     ****     ENMBAR=P0*EXPS*MBAR/(K*T)€ @€\      DO 91 I = 1,LEN1  € €€\	ENMBAR(I)=C(81)*EXPS(K)*EMBAR(I,K)/(C(84)*F(I,NTK))    € À€\	FS11ST(I)=ENMBAR(I)**2*(-.5/.3*(RK13(K,I)+RK13(K+1,I))*€ € €\     1    (F(I,NPS2K)/RMASS(2))**2-1./3.*(RK14(K,I)+RK14(K+1,I))*       €  @€\     2    F(I,NPS1K)/RMASS(1)*F(I,NPS2K)/RMASS(2)+.5*(RK15(K,I)+€!€€€\     3    RK15(K+1,I))*(-.5*F(I,NPS2K)/(RMASS(2)*EMBAR(I,K))+1./3.*     €€À€\     4    (F(I,NPS2K)/RMASS(2))**2+2./3.*F(I,NPS1K)/RMASS(1)*   €€  €\     5    F(I,NPS2K)/RMASS(2))) €€ @€\	FS21ST(I)=.5*(FS21(K,I)+FS21(K+1,I))+FS11ST(I) €€ €€\	FS11ST(I)=.5*(FS11(K,I)+FS11(K+1,I))+FS11ST(I) €  À€\	FS12ST(I)=-1./3.*(RK13(K,I)+RK13(K+1,I))*F(I,NPS1K)/RMASS(1)*  €$ ! €\     1    F(I,NPS2K)/RMASS(2)-1./6.*(RK14(K,I)+RK14(K+1,I))*(F(I,NPS1K)/€# !@€\     2    RMASS(1))**2+.5*(RK15(K,I)+RK15(K+1,I))*F(I,NPS1K)/RMASS(1)*  €"€!€€\     3    (-.5/EMBAR(I,K)+1./3.*F(I,NPS1K)/RMASS(1)+2./3.*F(I,NPS2K)/   €	€!À€\     4    RMASS(2))     €  " €\	FS22ST(I)=.5*(FS22(K,I)+FS22(K+1,I))+ENMBAR(I)**2*(-(RK12(K,I)+€! "@€\     1    RK12(K+1,I))*F(I,NPS2K)/(RMASS(2)*EMBAR(I,K))+FS12ST(I))      € €"€€\	FS12ST(I)=.5*(FS12(K,I)+FS12(K+1,I))+ENMBAR(I)**2*(.5*(RK12(K,I)       €!€"À€\     1    +RK12(K+1,I))*F(I,NPS2K)/(RMASS(2)*EMBAR(I,K))+FS12ST(I))     € # €\C       ****    €€#@€\C       ****     TOTAL SOURCE   € #€€\C       ****    €!€#À€\C       ST(I)=.5*(ENMBAR(I)*((2.*(FS11(K,I)+FS11(K+1,I))+FS21(K,I)+     €#€$ €\C    1    FS21(K+1,I))*F(I,NPS1K)/RMASS(1)+(2.*(FS12(K,I)+FS12(K+1,I))+ €!€$@€\C    2    FS22(K,I)+FS22(K+1,I))*F(I,NPS2K)/RMASS(2))+2.*(FS1(K,I)+     €# $€€\C    3    FS1(K+1,I))+FS2(K,I)+FS2(K+1,I)-3.*ENMBAR(I)**3*((RK13(K,I)+  €" $À€\C    4    RK13(K+1,I))*F(I,NPS1K)/RMASS(1)*(F(I,NPS2K)/RMASS(2))**2+    €# % €\C    5    (RK14(K,I)+RK14(K+1,I))*(F(I,NPS1K)/RMASS(1))**2*F(I,NPS2K)/  €€%@€\C    6    RMASS(2)+(RK15(K,I)+RK15(K+1,I))*F(I,NPS1K)/RMASS(1)* €" %€€\C    7    F(I,NPS2K)/RMASS(2)*(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)))*    €€%À€\C    8    RMASS(2)/ENMBAR       €€& €\	ST(I)=0.       € &@€\C       ****    €€&€€\C       ****     SET UP SOURCE MATRIX IN WKM2   € &À€\C       ****    € ' €\	WKM2(I,1,1)=FS11ST(I)-ST(I)    €€'@€\	WKM2(I,1,2)=2.*FS12ST(I)       €€'€€\	WKM2(I,2,1)=.5*FS21ST(I)       € 'À€\	WKM2(I,2,2)=FS22ST(I)-ST(I)    € ( €\   91 CONTINUE  €!€(@€\C     ****     FINISH CALCULATING AK(K+1/2) AND GENERATE PK, QK, RK     €€(€€\      NWK=NJ+NW+K       €€(À€\      DO 10 L=1,2       € ) €\	DO 10 M=1,2    €€)@€\	  DO 10 I=1,LEN1       € )€€\	    AK(I,L,M,KP)=AK(I,L,M,KP)*WKS1(I)  € )À€\	    PK(I,L,M)=(AK(I,L,M,KM)*(1./C(3)+EP(I,M,KM)/2.)-EXPS(K)    €€* €\     A        *(C(87)*DIFK(I,K)*DFACT*(1./C(3)-WKS2(I))+0.25*   €€*@€\     B        (F(I,NWK-1)+F(I,NWK)))*DELTA(L,M))/C(3)   € *€€\	    RK(I,L,M)=(AK(I,L,M,KP)*(1./C(3)-EP(I,M,KP)/2.)-EXPS(K)    €€*À€\     A        *(C(86)*DIFK(I,K+1)*DFACT*(1./C(3)+WKS3(I))-0.25* €€+ €\     B        (F(I,NWK-1)+F(I,NWK)))*DELTA(L,M))/C(3)   € +@€\	    QK(I,L,M)=-(AK(I,L,M,KM)*(1./C(3)-EP(I,M,KM)/2.)+  € €+€€\     A        AK(I,L,M,KP)*(1./C(3)+EP(I,M,KP)/2.))/C(3)+EXPS(K)*       €"€+À€\     B        (((C(86)*DIFK(I,K+1)*(1./C(3)-WKS3(I))+C(87)*DIFK(I,K)*   €€, €\     C        (1./C(3)+WKS2(I)))*DFACT/C(3)+C(7))*DELTA(L,M)-   € ,@€\     D        WKM2(I,L,M))      € ,€€\   10 CONTINUE  €€,À€\C     ****     CALCULATE HORIZONTAL ADVECTION IN FK ARRAY       €€- €\      DO 11 L=1,2       € -@€\	NPSL=NPS+(L-1)*KMAXP1  € -€€\	CALL ADVECL(NPSL,FK(1,L),K)    € -À€\   11 CONTINUE  € . €\C     ****      €€.@€\C     ****     ADD EXPLICIT SOURCE TERMS FK     € .€€\C     ****      €€.À€\      DO 111 I = 1,LEN1 €€/ €\	FK(I,1)=FK(I,1)-(FS1(K,I)+FS1(K+1,I))*RMASS(2)/ENMBAR(I)       € /@€\	FK(I,2)=FK(I,2)-.5*(FS2(K,I)+FS2(K+1,I))*RMASS(2)/ENMBAR(I)    € /€€\  111 CONTINUE  € /À€\C     ****     SHAPIRO SMOOTHER   (ELIMINATED)  €	 0 €\      DO 112 L=1,2      € 0@€\	NPSK=NPSNM+(L-1)*KMAXP1+K-1    €€0€€\	NPSNMK=NJ+NPSK €	 0À€\	DO 112 I=3,IMAXP2      €€1 €\	  WKV2(I,L)=F(I,NPSNMK)-C(26)*(F(I+2,NPSNMK)+F(I-2,NPSNMK)-4.* € 1@€\     1       (F(I+1,NPSNMK)+F(I-1,NPSNMK))+6.*F(I,NPSNMK))      € 1€€\  112 CONTINUE  € 1À€\C     ****     COMPLETE CALCULATION OF RHS IN FK€€2 €\      DO 12 L=1,2       €€2@€\	DO 12 I=3,IMAXP2       €€2€€\	  FK(I,L)=EXPS(K)*(WKV2(I,L)*C(7)-FK(I,L))     € 2À€\   12 CONTINUE  €€3 €\C     ****     INSERT PERIODIC POINTS   €
 3@€\      DO 121 L = 1,2    €€3€€\	DO 121 I = 1,2 €€3À€\	  FK(I,L) = FK(I+IMAX,L)       €€4 €\	  FK(I+IMAXP2,L) = FK(I+2,L)   € 4@€\  121 CONTINUE  €€4€€\      IF(K.EQ.1) GO TO 13       € 4À€\      IF(K.EQ.KMAX) GO TO 14    € 5 €\      GO TO 15  € 5@€\   13 CONTINUE  € 5€€\C     ****     LOWER BNDRY      €€5À€\      DO 16 L=1,2       € 6 €\	DO 16 M=1,2    €€6@€\	  DO 16 KK=1,2 €	€6€€\	    DO 16 I=1,LEN1     €€6À€\	      QK(I,L,M)=QK(I,L,M)+PK(I,L,KK)*B(KK,M)   € 7 €\   16 CONTINUE  €€7@€\      DO 33 L=1,2       € 7€€\	DO 33 M=1,2    €€7À€\	  DO 33 I=1,LEN1       € 8 €\	    FK(I,L)=FK(I,L)-PK(I,L,M)*FB(I,M)  €€8@€\	    PK(I,L,M)=0.       € 8€€\   33 CONTINUE  € 8À€\      GO TO 15  € 9 €\   14 CONTINUE  € 9@€\C     ****     UPPER BNDRY      €€9€€\      DO 17 L=1,2       € 9À€\	DO 17 M=1,2    €€: €\	  DO 17 I=1,LEN1       € :@€\	    QK(I,L,M)=QK(I,L,M)+(1.+.5*EP(I,M,KP)*C(3))/(1.-.5*€ :€€\     A        EP(I,M,KP)*C(3))*RK(I,L,M)€€:À€\	    RK(I,L,M)=0.       € ; €\   17 CONTINUE  € ;@€\   15 CONTINUE  €€;€€\C     ****     QK=ALFAK=QK-PK*GAMA(K-1) €€;À€\      DO 18 L=1,2       € < €\	DO 18 M=1,2    €€<@€\	  DO 18 KK=1,2 €	€<€€\	    DO 18 I=1,LEN1     € <À€\	      QK(I,L,M)=QK(I,L,M)-PK(I,L,KK)*GAMA(I,K,KK,M)    € = €\   18 CONTINUE  €€=@€\C     ****     WKS1=DET(ALFA)   €
 =€€\      DO 19 I=1,LEN1    € =À€\	WKS1(I)=QK(I,1,1)*QK(I,2,2)-QK(I,1,2)*QK(I,2,1)€ > €\   19 CONTINUE  €€>@€\C     ****     WKM1=ALFAI       €€>€€\      DO 20 L=1,2       €€>À€\	LL=3-L € ? €\	DO 20 M=1,2    €€?@€\	  DO 20 I=1,LEN1       €€?€€\	    WKM1(I,L,M)=(DELTA(L,M)*QK(I,LL,LL)-(1.-DELTA(L,M))*       € ?À€\     A        QK(I,L,M))/WKS1(I)€ @ €\   20 CONTINUE  €€@@€\C     ****     GAMA(K+1)=ALFAI*RK       € @€€\C     ****     WKV1=FK-PK*Z(K)  €€@À€\      DO 21 L=1,2       €€A €\	DO 22 I=1,LEN1 €
 A@€\	  WKV1(I,L)=FK(I,L)    € A€€\   22 CONTINUE  €€AÀ€\      DO 21 M=1,2       €€B €\	DO 23 I=1,LEN1 €
€B@€\	  GAMA(I,K+1,L,M)=0.   €€B€€\	  WKV1(I,L)=WKV1(I,L)-PK(I,L,M)*Z(I,K,M)       € BÀ€\   23   CONTINUE€€C €\	DO 21 KK=1,2   €€C@€\	  DO 21 I=1,LEN1       € C€€\	    GAMA(I,K+1,L,M)=GAMA(I,K+1,L,M)+WKM1(I,L,KK)*RK(I,KK,M)    € CÀ€\   21 CONTINUE  €€D €\C     ****     Z(K+1)=WKM1*WKV1 €	 D@€\      DO 231 L=1,2      € D€€\	DO 232 I=1,LEN1€ DÀ€\	  Z(I,K+1,L)=0.€ E €\  232   CONTINUE€€E@€\	DO 231 M=1,2   €	 E€€\	  DO 231 I=1,LEN1      € EÀ€\	    Z(I,K+1,L)=Z(I,K+1,L)+WKM1(I,L,M)*WKV1(I,M)€ F €\  231 CONTINUE  € F@€\    6 CONTINUE  €€F€€\C     ****     SET PSNP(KMAXP1) TO ZERO €€FÀ€\      DO 24 L=1,2       € G €\	NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX€€G@€\	DO 24 I=1,LEN1 €€G€€\	  F(I,NPSL)=0. € GÀ€\   24 CONTINUE  €€H €\C     ****     DOWNWARD SWEEP   €
€H@€\      DO 25 KK=1,KMAX   € H€€\	K=KMAX+1-KK    € HÀ€\	DO 25 L=1,2    €€I €\	  NPSL=NJNP+NPS+(L-1)*KMAXP1+(K-1)     €€I@€\	  DO 26 I=1,LEN1       €€I€€\	    F(I,NPSL)=Z(I,K+1,L)       €	 IÀ€\   26     CONTINUE      € J €\	  DO 25 M=1,2  €€J@€\	    NPSM=NJNP+NPS+(M-1)*KMAXP1+K       €	€J€€\	    DO 25 I=1,LEN1     € JÀ€\	      F(I,NPSL)=F(I,NPSL)-GAMA(I,K+1,L,M)*F(I,NPSM)    € K €\   25 CONTINUE  €€K@€\C     ****     INSERT VALUE OF PS(KMAXP1) USING BNDRY CONDITION €€K€€\      DO 27 L=1,2       € KÀ€\	NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX€€L €\	DO 27 I=1,LEN1 €€L@€\	  F(I,NPSL)=(1.+.5*EP(I,L,KP)*C(3))/(1.-.5*EP(I,L,KP)*C(3))*   €€L€€\     A      F(I,NPSL-1) € LÀ€\   27 CONTINUE  € M €\C     ****     TIME SMOOTHING OF PSI    €€M@€\      NPSNMK=NJ+NPSNM-1 €	€M€€\      NPSK=NJ+NPS-1     €€MÀ€\      NPSNPK=NJNP+NPS-1 €€N €\      NPSMNK=NJNP+NPSNM-1       € N@€\      DO 29 K = 1,2*KMAXP1      €
 N€€\      DO 29 I=1,LEN1    € NÀ€\	F(I,NPSMNK+K)=C(30)*F(I,NPSK+K)+C(31)*(F(I,NPSNMK+K)+  € O €\     1  F(I,NPSNPK+K))  € O@€\   29 CONTINUE  €€O€€\C     ****     INSERT PERIODIC POINTS   €€OÀ€\      K1=NJNP+NPS       € P €\      K2=K1+2*KMAXP1-1  €€P@€\      DO 31 N=1,2       € P€€\	DO 32 I=1,2    € PÀ€\	  DO 32 K=K1,K2€
€Q €\	  F(I,K)=F(I+IMAX,K)   €€Q@€\	  F(I+IMAXP2,K)=F(I+2,K)       € Q€€\   32   CONTINUE€ QÀ€\	K1=NJNP+NPSNM  €€R €\	K2=K1+2*KMAXP1-1       € R@€\   31 CONTINUE  € R€€\C     ****     INSURE NON-NEGATIVE PSI  € RÀ€\      NPS1K=NJNP+NPS-1  €€S €\      NPS2K=NJNP+NPS2-1 €€S@€\      NPS1MK=NJNP+NPSNM-1       € S€€\      NPS2MK=NJNP+NPS2NM-1      € SÀ€\      DO 41 K=1,KMAXP1  €
 T €\      DO 41 I=1,LEN1    € T@€\	F(I,NPS1K+K)=CVMGP(F(I,NPS1K+K),SMALL,F(I,NPS1K+K)-SMALL)      € T€€\	F(I,NPS2K+K)=CVMGP(F(I,NPS2K+K),SMALL,F(I,NPS2K+K)-SMALL)      €€TÀ€\	F(I,NPS1MK+K)=CVMGP(F(I,NPS1MK+K),SMALL,F(I,NPS1MK+K)-SMALL)   €€U €\	F(I,NPS2MK+K)=CVMGP(F(I,NPS2MK+K),SMALL,F(I,NPS2MK+K)-SMALL)   €€U@€\	F(I,NPS1K+K)=F(I,NPS1K+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1K+K)+   €€U€€\     1    F(I,NPS2K+K)),1.-SMALL-F(I,NPS1K+K)-F(I,NPS2K+K))     €€UÀ€\	F(I,NPS2K+K)=F(I,NPS2K+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1K+K)+   €€V €\     1    F(I,NPS2K+K)),1.-SMALL-F(I,NPS1K+K)-F(I,NPS2K+K))     €  V@€\	F(I,NPS1MK+K)=F(I,NPS1MK+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1MK+K)+€ V€€\     1  F(I,NPS2MK+K)),1.-SMALL-F(I,NPS1MK+K)-F(I,NPS2MK+K))    €  VÀ€\	F(I,NPS2MK+K)=F(I,NPS2MK+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1MK+K)+€ W €\     1  F(I,NPS2MK+K)),1.-SMALL-F(I,NPS1MK+K)-F(I,NPS2MK+K))    € W@€\   41 CONTINUE  € W€€\      RETURN    €€WÀ€\      END       € X €\CDIR$ NOLIST            €  @€]ªDK FACE€  €€]CDIR$ NOLIST    €€ À€]      SUBROUTINE FACE2(STEP,DS,SB,DX)   €  €]      SAVE      € @€]C     ****      € €€]C     ****     INTERFACE€ À€]C     ****      €  €]C     ****     COMMON DECKS:    €! @€]C     ****       PARAMS, ASAS, CONS, CRATES, FIELDS, FIELDT, FLUX,      €$ €€]C     ****         HOXSTUF, INDEX, INDEXDAT, NEWRAT, PSIB, SOURCE, VSCR3€ À€]C     ****      €  €]ªCA PARAMS      €€@€]ªCA ASAS2       €€€€]ªCA CONS2       €€À€]ªCA CRATES2     €  €]ªCA DISSHOX2    €€@€]ªCA FIELDS2     €€€€]ªCA FIELDT2     €€À€]ªCA FLUX2       €  €]ªCA HOXSTUF2    €€@€]ªCA INDEX       € €€]ªCA INDEXDAT    €€À€]ªCA NEWRAT2     €€ €]ªCA NOZNOZ2     €€@€]ªCA OXOXOX2     € €€]ªCA PHOTO2      €€À€]ªCA PHOTOO2     €€ €]ªCA PSIB2       € @€]ªCA RATEBLK2    € €€]ªCA RJH2O2O2    €€À€]ªCA SOURCE2     €€ €]ªCA VSCR3       €€@€]      COMMON/SODIUMP/PNA(KMXP,IMXP),XLNA(KMXP,IMXP),XLBNA(IMXP) € €€]     1,PNAI(KMXP,IMXP),XLNAI(KMXP,IMXP),XLBNAI(IMXP)    € À€]     2,RNAO(KMXP,IMXP),RNAO2(KMXP,IMXP),RNAOH(KMXP,IMXP)€€	 €]      DIMENSION IFP(12) €€	@€]      DIMENSION FM(IMXP),DH2OT(KMXP,IMXP)       €€	€€]      EQUIVALENCE (C(46),PS1B),(C(47),PS3B)     €€	À€]      DATA PS1B,PS3B/.22,.78/   €€
 €]      DATA ISTAR/0/,PI/3.141592654/     € 
@€]      IF(ISTAR.EQ.0)THEN€ 
€€]C       ****    € 
À€]C       ****     SET INITIAL CONSTANTS  €  €]C       ****     CONVERT NUMBER DENSITIES TO PSI'S      € @€]C       ****    € €€]	ISTAR=1€€À€]	KMAX=KMX       €  €]	KMAXP1=KMAX+1  €€@€]	IMAX=IMX       € €€]	IMAXP2=IMAX+2  € À€]	IMAXP4=IMAX+4  €  €]	LEN1=IMXP      €€@€]	LEN2=LEN1*KMAX €€€€]	LEN3=LEN1*KMAXP1       €€À€]C       ****     SET DLAMDA     €  €]	DLAMDA=8.*ATAN(1.)/IMAX€ @€]	C(1)=DLAMDA    €€€€]C       ****     FINITE DIFFERENCE OPERATOR CONSTANTS   €
€À€]	C(10)=2./(3.*DLAMDA)   €  €]	C(11)=1./(12.*DLAMDA)  € @€]C       ****     SHAPIRO SMOOTHER CONSTANT      €
 €€]C       C(26)=3.0E-2    €€À€]	C(26)=0.       €  €]C       ****     TIME SMOOTHING CONSTANT€ @€]	ALPHA=.95      € €€]	C(30)=ALPHA    €
 À€]	C(31)=.5*(1.-ALPHA)    €€ €]C       ****       B MATRIX     € @€]	C(40)=-1.      €€€€]	C(41)=0.       €€À€]	C(42)=0.       €  €]	C(43)=-1.      € @€]C       ****        ATOMIC OXYGEN RECOMBINATION RATE    € €€]	C(44)=3.8E-30  €€À€]C       ****        FB VECTOR   €€ €]	DO 12 I = 1,LEN1       €	 @€]	  FB(I,1)=2*C(46)      € €€]   12   CONTINUE€€À€]C       ****    EARTHS RADIUS   €€ €]C       C(51)=6.3 122E8 €
€@€]	C(51)=IMX*DX/(2.*PI)   €€€€]	C(52)=1./C(51) €	 À€]	C(53)=C(51)*C(51)      €  €]C       ****    ACCELERATION DU TO GRAVITY      €€@€]	C(54)=870.     € €€]C       ****    GAS CONSTANT    € À€]	C(57)=8.314E7  €€ €]C       ****         P0 € @€]	C(81)=5.E-4    € €€]C       ****        BOLTZMANN CONSTANT  €€À€]	C(84)=1.38E-16 €  €]C       ****     AVOGRADO'S NUMBER      €€@€]	C(85)=6.023E23 €€€€]	RMASS(1)=32.   €€À€]	RMASS(2)=16.   €€ €]	RMASS(3)=28.   € @€]	RMA=23.€ €€]	RMH=23.€€À€]	RMNO=30.       €  €]	RMN4S=14.      € @€]	RMCH4=16.      € €€]	RMH2=2.€ À€]	RMCO2=44.      €€ €]	RMCO=28.       € @€]	RMH2O=18.      €€€€]	RMHOX=1.       €	€À€]C       ****     T0     €  €]	DO 2 K=1,KMAXP1€€@€]	  T0(K)=0.     € €€]    2   CONTINUE€ À€]C       ****    €  €]C       ****     SET INITIAL PSI'S      € @€]C       ****    €€€€]C       ****     WK1=PSI1, WK2=PSI2, WK3=PSAR,  WK4=PSHE,       € À€]C       ****     WK5=PNO,  WK6=PN4S, WK7=PCH4,  WK8=PH2,€€ €]C       ****     WK9=PCO2, WK10=PCO, WK11=PHOX, WK12=PH2O       € @€]C       ****    € €€]	DO 3 K=1,KMAXP1€	 À€]	  DO 3 I=3,LEN1-1      €  €]	    FNNOZ(K,I)=FNNO(K,I)+FNNO2(K,I)    €€@€]	    WK1(K,I)=FNO2(K,I)*RMASS(1)+FNO(K,I)*RMASS(2)+FNN2(K,I)*   € €€]     1        RMASS(3)  €€À€]	    WK2(K,I)=FNOX(K,I)*RMASS(2)/WK1(K,I)       €  €]	    WK3(K,I)=FNAR(K,I)*RMA/WK1(K,I)    € @€]	    WK4(K,I)=FNHE(K,I)*RMH/WK1(K,I)    € €€]	    WK5(K,I)=FNNOZ(K,I)*RMNO/WK1(K,I)  €€À€]	    WK6(K,I)=FNN4S(K,I)*RMN4S/WK1(K,I) €€  €]	    WK7(K,I)=FNCH4(K,I)*RMCH4/WK1(K,I) €€ @€]	    WK8(K,I)=FNH2(K,I)*RMH2/WK1(K,I)   €€ €€]	    WK9(K,I)=FNCO2(K,I)*RMCO2/WK1(K,I) €  À€]	    WK10(K,I)=FNCO(K,I)*RMCO/WK1(K,I)  € ! €]	    WK11(K,I)=FNHOX(K,I)*RMHOX/WK1(K,I)€ !@€]	    WK12(K,I)=FNH2O(K,I)*RMH2O/WK1(K,I)€€!€€]	    WK1(K,I)=FNO2(K,I)*RMASS(1)/WK1(K,I)       € !À€]    3   CONTINUE€ " €]C       ****    €€"@€]C       ****    TRANSFER PSI'S TO F-ARRAY       € "€€]C       ****    €€"À€]	NPS1K=NJ+NPS-1 €	€# €]	NPS2K=NPS1K+KMAXP1     € #@€]	NPSAK=NJ+NPSA-1€ #€€]	NPSHK=NJ+NPSH-1€ #À€]	NPNOK=NJ+NPNO-1€	 $ €]	NPN4SK=NJ+NPN4S-1      €	 $@€]	NPCH4K=NJ+NPCH4-1      € $€€]	NPH2K=NJ+NPH2-1€	 $À€]	NPCO2K=NJ+NPCO2-1      € % €]	NPCOK=NJ+NPCO-1€	 %@€]	NPHOXK=NJ+NPHOX-1      €	 %€€]	NPH2OK=NJ+NPH2O-1      € %À€]	DO 4 K=1,KMAX  €	 & €]	  DO 4 I=3,LEN1-1      € &@€]	    F(I,NPS1K+K)=.5*(WK1(K,I)+WK1(K+1,I))      € &€€]	    F(I,NPS2K+K)=.5*(WK2(K,I)+WK2(K+1,I))      € &À€]	    F(I,NPSAK+K)=.5*(WK3(K,I)+WK3(K+1,I))      € ' €]	    F(I,NPSHK+K)=.5*(WK4(K,I)+WK4(K+1,I))      € '@€]	    F(I,NPNOK+K)=.5*(WK5(K,I)+WK5(K+1,I))      €€'€€]	    F(I,NPN4SK+K)=.5*(WK6(K,I)+WK6(K+1,I))     €€'À€]	    F(I,NPCH4K+K)=.5*(WK7(K,I)+WK7(K+1,I))     € ( €]	    F(I,NPH2K+K)=.5*(WK8(K,I)+WK8(K+1,I))      €€(@€]	    F(I,NPCO2K+K)=.5*(WK9(K,I)+WK9(K+1,I))     € (€€]	    F(I,NPCOK+K)=.5*(WK10(K,I)+WK10(K+1,I))    €€(À€]	    F(I,NPHOXK+K)=.5*(WK11(K,I)+WK11(K+1,I))   €€) €]	    F(I,NPH2OK+K)=.5*(WK12(K,I)+WK12(K+1,I))   € )@€]    4   CONTINUE€€)€€]C       ****     BOUNDARY       €	 )À€]	DO 9 I = 3,LEN1-1      €€* €]	  F(I,NPS1K+KMAXP1)=1.5*WK1(KMAXP1,I)-.5*WK1(KMAX,I)   €€*@€]	  F(I,NPS2K+KMAXP1)=1.5*WK2(KMAXP1,I)-.5*WK2(KMAX,I)   €€*€€]	  F(I,NPSAK+KMAXP1)=1.5*WK3(KMAXP1,I)-.5*WK3(KMAX,I)   €€*À€]	  F(I,NPSHK+KMAXP1)=1.5*WK4(KMAXP1,I)-.5*WK4(KMAX,I)   €€+ €]	  F(I,NPNOK+KMAXP1)=1.5*WK5(KMAXP1,I)-.5*WK5(KMAX,I)   € +@€]	  F(I,NPN4SK+KMAXP1)=1.5*WK6(KMAXP1,I)-.5*WK6(KMAX,I)  € +€€]	  F(I,NPCH4K+KMAXP1)=1.5*WK7(KMAXP1,I)-.5*WK7(KMAX,I)  €€+À€]	  F(I,NPH2K+KMAXP1)=1.5*WK8(KMAXP1,I)-.5*WK8(KMAX,I)   € , €]	  F(I,NPCO2K+KMAXP1)=1.5*WK9(KMAXP1,I)-.5*WK9(KMAX,I)  €€,@€]	  F(I,NPCOK+KMAXP1)=1.5*WK10(KMAXP1,I)-.5*WK10(KMAX,I) € ,€€]	  F(I,NPHOXK+KMAXP1)=1.5*WK11(KMAXP1,I)-.5*WK11(KMAX,I)€ ,À€]	  F(I,NPH2OK+KMAXP1)=1.5*WK12(KMAXP1,I)-.5*WK12(KMAX,I)€€- €]	  F(I,NPS1K+KMAXP1)=1.5*WK1(KMAXP1,I)-.5*WK1(KMAX,I)   €€-@€]	  F(I,NPS2K+KMAXP1)=1.5*WK2(KMAXP1,I)-.5*WK2(KMAX,I)   €€-€€]	  F(I,NPSAK+KMAXP1)=1.5*WK3(KMAXP1,I)-.5*WK3(KMAX,I)   €€-À€]	  F(I,NPSHK+KMAXP1)=1.5*WK4(KMAXP1,I)-.5*WK4(KMAX,I)   €€. €]	  F(I,NPNOK+KMAXP1)=1.5*WK5(KMAXP1,I)-.5*WK5(KMAX,I)   € .@€]	  F(I,NPN4SK+KMAXP1)=1.5*WK6(KMAXP1,I)-.5*WK6(KMAX,I)  € .€€]	  F(I,NPCH4K+KMAXP1)=1.5*WK7(KMAXP1,I)-.5*WK7(KMAX,I)  €€.À€]	  F(I,NPH2K+KMAXP1)=1.5*WK8(KMAXP1,I)-.5*WK8(KMAX,I)   € / €]	  F(I,NPCO2K+KMAXP1)=1.5*WK9(KMAXP1,I)-.5*WK9(KMAX,I)  €€/@€]	  F(I,NPCOK+KMAXP1)=1.5*WK10(KMAXP1,I)-.5*WK10(KMAX,I) € /€€]	  F(I,NPHOXK+KMAXP1)=1.5*WK11(KMAXP1,I)-.5*WK11(KMAX,I)€ /À€]	  F(I,NPH2OK+KMAXP1)=1.5*WK12(KMAXP1,I)-.5*WK12(KMAX,I)€ 0 €]    9   CONTINUE€ 0@€]C       ****    € 0€€]C       ****     PERIODIC POINTS€ 0À€]C       ****    € 1 €]	NPSK = NJ+NPS-1€ 1@€]	DO 16 I=1,2    € 1€€]	  DO 16 K=1,12*KMAXP1  €€1À€]	    F(I,NPSK+K) = F(I+IMAX,NPSK+K)     €€2 €]	    F(I+IMAXP2,NPSK+K) = F(I+2,NPSK+K) € 2@€]   16   CONTINUE€€2€€]	IFP(1) = NPS   € 2À€]	IFP(2) = NPS2  € 3 €]	IFP(3) = NPSA  € 3@€]	IFP(4) = NPSH  € 3€€]	IFP(5) = NPNO  €€3À€]	IFP(6) = NPN4S €€4 €]	IFP(7) = NPCH4 € 4@€]	IFP(8) = NPH2  €€4€€]	IFP(9) = NPCO2 €€4À€]	IFP(10) = NPCO € 5 €]	IFP(11) = NPHOX€ 5@€]	IFP(12) = NPH2O€ 5€€]C       CALL PLOT(IFP,12,NJ)    € 5À€]C       ****    €€6 €]C       ****     COPY CURRENT FIELDS TO NM ARRAYS       € 6@€]C       ****    € 6€€]	NPSK = NJ+NPS-1€
 6À€]	NPSNMK = NJ+NPSNM-1    € 7 €]	DO 14 K = 1,12*KMAXP1  €€7@€]	DO 14 I = 1,LEN1       € 7€€]	  F(I,NPSNMK+K) = F(I,NPSK+K)  € 7À€]   14   CONTINUE€€8 €]	IFP(1) = NPSNM € 8@€]	IFP(2) = NPS2NM€ 8€€]	IFP(3) = NPSANM€ 8À€]	IFP(4) = NPSHNM€ 9 €]	IFP(5) = NPNONM€ 9@€]	IFP(6) = NN4SNM€ 9€€]	IFP(7) = NCH4NM€ 9À€]	IFP(8) = NPH2NM€ : €]	IFP(9) = NCO2NM€€:@€]	IFP(10) = NPCONM       €€:€€]	IFP(11) = NHOXNM       €€:À€]	IFP(12) = NH2ONM       € ; €]C       CALL PLOT(IFP,12,NJ)    €€;@€]      ENDIF     €€;€€]      C(3)=DS   €	 ;À€]      C(7)=1./STEP      € < €]C     ****     EXP(-.5*DS)      €€<@€]      C(86)=EXP(-.5*DS) €€<€€]C     ****     EXP(.5*DS)       €
 <À€]      C(87)=1./C(86)    €€= €]C     ****     EXPS, DIFK       € =@€]      S=SB+.5*DS€
€=€€]      EXPS(1)=EXP(-S)   €
 =À€]      EXPDS=EXP(-DS)    €
€> €]      DO 1 K=2,KMAXP1   € >@€]	EXPS(K)=EXPS(K-1)*EXPDS€ >€€]    1 CONTINUE  € >À€]      DO 15 K = 1,KMAXP1€	€? €]	DO 15 I = 3,LEN1-2     €€?@€]	  DIFK(I,K) = FDIFK(K,I)       € ?€€]   15 CONTINUE  € ?À€]C     ****      € @ €]C     ****     PERIODIC POINTS FOR T, U, W, RJ, PS2B, DIFK      € @@€]C     ****      €	€@€€]      DO 17 I = 1,2     €€@À€]	PS2B(I) = PS2B(I+IMAX) €€A €]	PS2B(I+IMAXP2) = PS2B(I+2)     €	€A@€]	DO 17 K = 1,KMAXP1     €€A€€]	  DIFK(I,K) = DIFK(I+IMAX,K)   €€AÀ€]	  DIFK(I+IMAXP2,K) = DIFK(I+2,K)       €€B €]	  FT(K,I) = FT(K,I+IMAX)       €€B@€]	  FT(K,I+IMAXP2) = FT(K,I+2)   €€B€€]	  FUI(K,I) = FUI(K,I+IMAX)     €€BÀ€]	  FU(K,I) = FU(K,I+IMAX)       €€C €]	  FUI(K,I+IMAXP2) = FUI(K,I+2) €€C@€]	  FU(K,I+IMAXP2) = FU(K,I+2)   €€C€€]	  FW(K,I) = FW(K,I+IMAX)       €€CÀ€]	  FWI(K,I) = FWI(K,I+IMAX)     €€D €]	  FW(K,I+IMAXP2) = FW(K,I+2)   €€D@€]	  FWI(K,I+IMAXP2) = FWI(K,I+2) €€D€€]	  FRJ(K,I) = FRJ(K,I+IMAX)     €€DÀ€]	  FRJ(K,I+IMAXP2) = FRJ(K,I+2) € E €]   17 CONTINUE  € E@€]C     ****      €€E€€]C     ****     PERIODIC POINTS FOR SOURCE TERMS € EÀ€]C     ****      €	 F €]C     CALL CHKBLKS      €	€F@€]      DO 18 I = 1,2     €	 F€€]	DO 18 K =1,KMAXP1      €€FÀ€]	  FS11(K,I) = FS11(K,I+IMAX)   €€G €]	  FS11(K,I+IMAXP2) = FS11(K,I+2)       €€G@€]	  FS12(K,I) = FS12(K,I+IMAX)   €€G€€]	  FS12(K,I+IMAXP2) = FS12(K,I+2)       €€GÀ€]	  FS21(K,I) = FS21(K,I+IMAX)   €€H €]	  FS21(K,I+IMAXP2) = FS21(K,I+2)       €€H@€]	  FS22(K,I) = FS22(K,I+IMAX)   €€H€€]	  FS22(K,I+IMAXP2) = FS22(K,I+2)       €€HÀ€]	  FS1(K,I) = FS1(K,I+IMAX)     €€I €]	  FS1(K,I+IMAXP2) = FS1(K,I+2) €€I@€]	  FS2(K,I) = FS2(K,I+IMAX)     €€I€€]	  FS2(K,I+IMAXP2) = FS2(K,I+2) € IÀ€]   18 CONTINUE  €	 J €]C     CALL CHKBLKS      € J@€]C     ****      € J€€]C     ****     PERIODIC POINTS RATE COEFFICIENTS€ JÀ€]C     ****      €	€K €]      DO 19 I = 1,2     €€K@€]	DO 20 N = 1,25 €
€K€€]	  DO 20 K = 1,KMAXP1   €€KÀ€]	    RKK(N,K,I) = RKK(N,K,I+IMAX)       €€L €]	    RKK(N,K,I+IMAXP2) = RKK(N,K,I+2)   € L@€]   20   CONTINUE€€L€€]	DO 21 N = 1,10 €
€LÀ€]	  DO 21 K = 1,KMAXP1   €€M €]	    ALPP(N,K,I) = ALPP(N,K,I+IMAX)     €€M@€]	    ALPP(N,K,I+IMAXP2) = ALPP(N,K,I+2) € M€€]   21   CONTINUE€€MÀ€]	DO 22 N = 1,20 €
€N €]	  DO 22 K = 1,KMAXP1   €€N@€]	    RBB(N,K,I) = RBB(N,K,I+IMAX)       €€N€€]	    RBB(N,K,I+IMAXP2) = RBB(N,K,I+2)   € NÀ€]   22   CONTINUE€€O €]	DO 23 N = 1,50 €
€O@€]	  DO 23 K = 1,KMAXP1   €€O€€]	    RKMM(N,K,I) = RKMM(N,K,I+IMAX)     €€OÀ€]	    RKMM(N,K,I+IMAXP2) = RKMM(N,K,I+2) € P €]   23   CONTINUE€	€P@€]	DO 24 K = 1,KMAXP1     €€P€€]	  RK12(K,I) = RK12(K,I+IMAX)   €€PÀ€]	  RK12(K,I+IMAXP2) = RK12(K,I+2)       €€Q €]	  RK13(K,I) = RK13(K,I+IMAX)   €€Q@€]	  RK13(K,I+IMAXP2) = RK13(K,I+2)       €€Q€€]	  RK14(K,I) = RK14(K,I+IMAX)   €€QÀ€]	  RK14(K,I+IMAXP2) = RK14(K,I+2)       €€R €]	  RK15(K,I) = RK15(K,I+IMAX)   €€R@€]	  RK15(K,I+IMAXP2) = RK15(K,I+2)       € R€€]   24   CONTINUE€ RÀ€]   19 CONTINUE  € S €]C     ****      € S@€]C     ****     NTRANSFER T AND U€ S€€]C     ****     N€€SÀ€]      NTK=NJ+NT-1       €€T €]      NUK=NJ+NU-1       €	€T@€]      DO 5 K=1,KMAX     € T€€]	DO 5 I=1,LEN1  € TÀ€]	  F(I,NTK+K)=.5*(FT(K,I)+FT(K+1,I))    € U €]	  F(I,NUK+K)=.5*(FU(K,I)+FU(K+1,I))    € U@€]    5 CONTINUE  €€U€€]C     ****     BOUNDARY € UÀ€]      DO 10 I = 1,LEN1  € V €]	F(I,NTK+KMAXP1)=FT(1,I)€ V@€]	F(I,NUK+KMAXP1)=FU(1,I)€ V€€]   10 CONTINUE  € VÀ€]C     ****      € W €]C     ****     TRANSFER W AND RJ€ W@€]C     ****      €€W€€]      NWK=NJ+NW-1       €
 WÀ€]      NRJK=NDJ+NRJ-1    €
€X €]      DO 6 K=1,KMAXP1   € X@€]	DO 6 I=1,LEN1  €
€X€€]	  F(I,NWK+K)=FW(K,I)   €€XÀ€]	  F(I,NRJK+K)=FRJ(K,I) € Y €]    6 CONTINUE  € Y@€]C     ****      €€Y€€]C     ****     UPDATE PSI'S     € YÀ€]C     ****      €	 Z €]C     CALL CHKBLKS      €
 Z@€]      CALL RATES(FT)    €	 Z€€]C     CALL CHKBLKS      € ZÀ€]C     ****      €€[ €]C     ****     LOWER BOUNDARY CONDITION € [@€]C     ****      € [€€]      DO 11 I = 1,LEN1  €"€[À€]C       ENMBAR(I)=C(81)*EXPS(1)*C(87)/((PS1B/RMASS(1)+PS3B/RMASS(3))*   € \ €]C    1    C(84)*FT(1,I))€#€\@€]C       AA(I)=ENMBAR(I)**3*(2.*RK12(1,I)*(PS1B/RMASS(1)+PS3B/RMASS(3))+ €" \€€]C    1    RK13(1,I)*PS1B/RMASS(1)-RK15(1,I)*RMASS(2)*PS1B/(RMASS(3)*    €
 \À€]C    2    RMASS(1)))    €" ] €]C       BB(I)=ENMBAR(I)*(ENMBAR(I)**2*(RK14(1,I)*(PS1B/RMASS(1))**2+    €! ]@€]C    1    RK15(1,I)*(PS1B/RMASS(1))*(1.-PS1B)/RMASS(3))-FS22(1,I))      €€]€€]C       CC(I)=-FS21(1,I)*ENMBAR(I)*PS1B/RMASS(1)-FS2(1,I)       €" ]À€]C       PS2B(I)=RMASS(2)*(-BB(I)+SQRT(BB(I)**2-4.*AA(I)*CC(I)))/(2.*    € ^ €]C    1    AA(I))€ ^@€]	FB(I,2)=PS2B(I)€ ^€€]   11 CONTINUE  € ^À€]C     II=1      €
 _ €]C     WRITE(6,101)II    €
€_@€]C 101 FORMAT(* HI*I1)   €	 _€€]C     CALL CHKBLKS      €€_À€]      CALL COMP €	 ` €]C     CALL CHKBLKS      € `@€]      CALL PHOTO€	 `€€]C     CALL CHKBLKS      €€`À€]      CALL CMPN2D       €	 a €]C     CALL CHKBLKS      €€a@€]      CALL CMPN4S       €	 a€€]C     CALL CHKBLKS      €€aÀ€]      CALL COMPNO       €	 b €]C     CALL CHKBLKS      €	 b@€]      CALL COMPARG      €	 b€€]C     CALL CHKBLKS      €€bÀ€]C     CALL COMPHE       €	 c €]C     CALL CHKBLKS      €	 c@€]      CALL COMPCH4      €	 c€€]C     CALL CHKBLKS      €€cÀ€]      CALL COMPH2       €	 d €]C     CALL CHKBLKS      €	 d@€]      CALL COMPCO2      €	 d€€]C     CALL CHKBLKS      €€dÀ€]      CALL COMPCO       €	 e €]C     CALL CHKBLKS      €	 e@€]      CALL COMPHOX      €	 e€€]C     CALL CHKBLKS      €	 eÀ€]      CALL COMPH2O      €	 f €]C     CALL CHKBLKS      € f@€]C     ****      €€f€€]C     ****     NTRANSFER UI     € fÀ€]C     ****     N€€g €]      NUK=NJ+NU-1       €
 g@€]      DO 55 K=1,KMAX    €€g€€]	DO 55 I=1,LEN1 € gÀ€]	  F(I,NUK+K)=.5*(FUI(K,I)+FUI(K+1,I))  € h €]   55 CONTINUE  €€h@€]C     ****     BOUNDARY € h€€]      DO 67 I = 1,LEN1  €€hÀ€]	F(I,NUK+KMAXP1)=FUI(1,I)       € i €]   67 CONTINUE  € i@€]C     ****      € i€€]C     ****     TRANSFER WI      € iÀ€]C     ****      €€j €]      NWK=NJ+NW-1       € j@€]      DO 66 K=1,KMAXP1  €€j€€]	DO 66 I=1,LEN1 € jÀ€]	  F(I,NWK+K)=FWI(K,I)  € k €]   66 CONTINUE  €€k@€]      CALL COMPHE       €€k€€]	IFP(1) = NPS   € kÀ€]	IFP(2) = NPS2  € l €]	IFP(3) = NPSA  € l@€]	IFP(4) = NPSH  € l€€]	IFP(5) = NPNO  €€lÀ€]	IFP(6) = NPN4S €€m €]	IFP(7) = NPCH4 € m@€]	IFP(8) = NPH2  €€m€€]	IFP(9) = NPCO2 €€mÀ€]	IFP(10) = NPCO € n €]	IFP(11) = NPHOX€ n@€]	IFP(12) = NPH2O€ n€€]C       CALL PLOT(IFP,12,NJNP)  €€nÀ€]C     WRITE(6,1000)PNOB,PN4SB   € o €]C1000 FORMAT(* PNOB, PN4SB *,2E13.4)    € o@€]C     ****      € o€€]C     ****     CONVERT PSI'S TO NUMBER DENSITIES€ oÀ€]C     ****      €€p €]C     ****     WK1=INTERPOLATED PSI1,   WK2=INTERPOLATED PSI2   €€p@€]C     ****     WK3=INTERPOLATED PSAR,   WK4=INTERPOLATED PSHE   € p€€]C     ****     WK5=INTERPOLATED PSNO,   WK6=INTERPOLATED PSN4S  €€pÀ€]C     ****     WK7=INTERPOLATED PSCH4,  WK8=INTERPOLATED PSH2   € q €]C     ****     WK9=INTERPOLATED PSCO2,  WK10=INTERPOLATED PSCO  €€q@€]C     ****     WK11=INTERPOLATED PSHOX, WK12=INTERPOLATED PSH2O € q€€]C     ****      € qÀ€]      NPS1K=NJNP+NPS-1  € r €]      NPS2K=NPS1K+KMAXP1€€r@€]      NPSAK=NJNP+NPSA-1 €€r€€]      NPSHK=NJNP+NPSH-1 €€rÀ€]      NPNOK=NJNP+NPNO-1 €€s €]      NPN4SK=NJNP+NPN4S-1       €€s@€]      NPCH4K=NJNP+NPCH4-1       €€s€€]      NPH2K=NJNP+NPH2-1 €€sÀ€]      NPCO2K=NJNP+NPCO2-1       €€t €]      NPCOK=NJNP+NPCO-1 €€t@€]      NPHOXK=NJNP+NPHOX-1       €€t€€]      NPH2OK=NJNP+NPH2O-1       €	€tÀ€]      DO 7 K=1,KMAX     € u €]	DO 7 I=1,LEN1  € u@€]	  WK1(K+1,I)=.5*(F(I,NPS1K+K)+F(I,NPS1K+K+1))  € u€€]	  WK2(K+1,I)=.5*(F(I,NPS2K+K)+F(I,NPS2K+K+1))  € uÀ€]	  WK3(K+1,I)=.5*(F(I,NPSAK+K)+F(I,NPSAK+K+1))  € v €]	  WK4(K+1,I)=.5*(F(I,NPSHK+K)+F(I,NPSHK+K+1))  € v@€]	  WK5(K+1,I)=.5*(F(I,NPNOK+K)+F(I,NPNOK+K+1))  € v€€]	  WK6(K+1,I)=.5*(F(I,NPN4SK+K)+F(I,NPN4SK+K+1))€ vÀ€]	  WK7(K+1,I)=.5*(F(I,NPCH4K+K)+F(I,NPCH4K+K+1))€ w €]	  WK8(K+1,I)=.5*(F(I,NPH2K+K)+F(I,NPH2K+K+1))  € w@€]	  WK9(K+1,I)=.5*(F(I,NPCO2K+K)+F(I,NPCO2K+K+1))€€w€€]	  WK10(K+1,I)=.5*(F(I,NPCOK+K)+F(I,NPCOK+K+1)) €€wÀ€]	  WK11(K+1,I)=.5*(F(I,NPHOXK+K)+F(I,NPHOXK+K+1))       €€x €]	  WK12(K+1,I)=.5*(F(I,NPH2OK+K)+F(I,NPH2OK+K+1))       € x@€]    7 CONTINUE  €€x€€]C     ****     LOWER BOUNDARY   € xÀ€]      DO 13 I = 1,LEN1  € y €]	WK1(1,I)=PS1B  €€y@€]	WK2(1,I)=PS2B(I)       €€y€€]	WK3(1,I)=.5*(F(1,NPSAK+1)+PARB(I))     € yÀ€]	WK4(1,I)=.5*(F(1,NPSHK+1)+PHB(I))      €€z €]	WK5(1,I)=.5*(F(1,NPNOK+1)+PNOB(I))     €€z@€]	WK6(1,I)=.5*(F(1,NPN4SK+1)+PN4SB(I))   €€z€€]	WK7(1,I)=.5*(F(1,NPCH4K+1)+PCH4B(I))   €€zÀ€]	WK8(1,I)=.5*(F(1,NPH2K+1)+PH2B(I))     €€{ €]	WK9(1,I)=.5*(F(1,NPCO2K+1)+PCO2B(I))   € {@€]	WK10(1,I)=.5*(F(1,NPCOK+1)+PCOB(I))    € {€€]	WK11(1,I)=.5*(F(1,NPHOXK+1)+PHOXB(I))  € {À€]	WK12(1,I)=.5*(F(1,NPH2OK+1)+PH2OB(I))  € | €]   13 CONTINUE  € |@€]C     ****      € |€€]C     ****     CONVERT TO NUMBER DENSITIES      € |À€]C     ****      €€} €]C     ****     WK5=MBAR*P0*EXPS/(R*T)   €
€}@€]      DO 8 K=1,KMAXP1   € }€€]	DO 8 I=3,LEN1-1€ }À€]	  WK13(K,I)=C(81)*EXPS(K)*C(87)/(C(84)*FT(K,I)*(WK1(K,I)/      €$ ~ €]     1      RMASS(1)+WK2(K,I)/RMASS(2)+(1.-WK1(K,I)-WK2(K,I))/RMASS(3)))€ ~@€]	  FNO2(K,I)=WK1(K,I)*WK13(K,I)/RMASS(1)€ ~€€]	  FNOX(K,I)=WK2(K,I)*WK13(K,I)/RMASS(2)€€~À€]	  FNO(K,I)=OOXR(K,I)*FNOX(K,I) €€ €]	  FNO3(K,I)=O3OR(K,I)*FNO(K,I) € @€]	  FNN2(K,I)=(1.-WK1(K,I)-WK2(K,I))*WK13(K,I)/RMASS(3)  €€€€]	  FM(I)=FNO(K,I)+FNO2(K,I)+FNN2(K,I)   €€À€]	  FNAR(K,I)=WK3(K,I)*WK13(K,I)/RMA     €! € €]C         FNAS(K,I)=FNAR(K,I)/(1.+RNAO(K,I)+RNAO2(K,I)+RNAOH(K,I))      €€€@€]C         FNAO(K,I)=FNAS(K,I)*RNAO(K,I) €€€€€]C         FNAO2(K,I)=FNAS(K,I)*RNAO2(K,I)       €€€À€]C         FNAOH(K,I)=FNAS(K,I)*RNAOH(K,I)       €€ €]	  FNHE(K,I)=WK4(K,I)*WK13(K,I)/RMH     €€@€]	  FNNOZ(K,I)=WK5(K,I)*WK13(K,I)/RMNO   €€€€]	  FNNO(K,I)=RNONOZ(K,I)*FNNOZ(K,I)     €€À€]	  FNNO2(K,I)=SNO2NO(K,I)*FNNO(K,I)     € ‚ €]	  FNN4S(K,I)=WK6(K,I)*WK13(K,I)/RMN4S  € ‚@€]	  FNCH4(K,I)=WK7(K,I)*WK13(K,I)/RMCH4  € ‚€€]	  FNH2(K,I)=WK8(K,I)*WK13(K,I)/RMH2    € ‚À€]	  FNCO2(K,I)=WK9(K,I)*WK13(K,I)/RMCO2  €€ƒ €]	  FNCO(K,I)=WK10(K,I)*WK13(K,I)/RMCO   €€ƒ@€]	  FNHOX(K,I)=WK11(K,I)*WK13(K,I)/RMHOX €€ƒ€€]	  FNH2O(K,I)=WK12(K,I)*WK13(K,I)/RMH2O €€ƒÀ€]	  HO2HR(K,I)=RKMM(28,K,I)*FM(I)*FNO2(K,I)/(RKMM(18,K,I)*       €! „ €]     1      FNO(K,I)+RBB(12,K,I)*FNNO(K,I)+RKMM(26,K,I)*FNO3(K,I))      € „@€]	  OHHR(K,I)=((RKMM(18,K,I)*FNO(K,I)*RKMM(28,K,I)*FNO2(K,I)*    € €„€€]     1      FM(I))/(RKMM(18,K,I)*FNO(K,I)+RKMM(26,K,I)*FNO3(K,I)+       €€„À€]     2      RBB(12,K,I)*FNNO(K,I))+RKMM(29,K,I)*FNO3(K,I))/     € … €]     3      (RKMM(17,K,I)*FNO(K,I)+RKMM(21,K,I)*FNO3(K,I))      € …@€]	  HHOXR(K,I)=1./(1.+HO2HR(K,I)+OHHR(K,I))      €€…€€]	  FNH(K,I)=HHOXR(K,I)*FNHOX(K,I)       €€…À€]	  FNHO2(K,I)=HO2HR(K,I)*FNH(K,I)       €€† €]	  FNOH(K,I)=OHHR(K,I)*FNH(K,I) €€†@€]	  DH2OT(K,I)=SH2OSRC(K,I)+SH2OSRB(K,I)+.88*SH2OLYA(K,I)+       € †€€]     1      SH2OEUV(K,I)€"€†À€]C         FNH2O2(K,I)=RKMM(27,K,I)*FNHO2(K,I)*FNHO2(K,I)/(RJH2O2(K,I)   € ‡ €]C    1    +RKMM(19,K,I)*FNO(K,I)+RKMM(24,K,I)*FNOH(K,I))€
€‡@€]	  FNH2O2(K,I)=1.E-20   € ‡€€]    8 CONTINUE  € ‡À€]C*****************************  €"€ˆ €]C     CALL CONREC(FNO2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € ˆ@€]C     CALL FRAME€"€ˆ€€]C     CALL CONREC(FNOX(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € ˆÀ€]C     CALL FRAME€" ‰ €]C     CALL CONREC(FNO(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)    € ‰@€]C     CALL FRAME€"€‰€€]C     CALL CONREC(FNO3(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € ‰À€]C     CALL FRAME€"€Š €]C     CALL CONREC(FNN2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € Š@€]C     CALL FRAME€"€Š€€]C     CALL CONREC(FNAR(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € ŠÀ€]C     CALL FRAME€"€‹ €]C     CALL CONREC(FNHE(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € ‹@€]C     CALL FRAME€# ‹€€]C     CALL CONREC(FNNOZ(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € ‹À€]C     CALL FRAME€"€Œ €]C     CALL CONREC(FNNO(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € Œ@€]C     CALL FRAME€# Œ€€]C     CALL CONREC(FNNO2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € ŒÀ€]C     CALL FRAME€#  €]C     CALL CONREC(FNN4S(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € @€]C     CALL FRAME€# €€]C     CALL CONREC(FNCH4(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € À€]C     CALL FRAME€"€Ž €]C     CALL CONREC(FNH2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € Ž@€]C     CALL FRAME€# Ž€€]C     CALL CONREC(FNCO2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € ŽÀ€]C     CALL FRAME€"€ €]C     CALL CONREC(FNCO(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € @€]C     CALL FRAME€# €€]C     CALL CONREC(FNHOX(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € À€]C     CALL FRAME€#  €]C     CALL CONREC(FNH2O(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € @€]C     CALL FRAME€# €€]C     CALL CONREC(HO2HR(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € À€]C     CALL FRAME€"€‘ €]C     CALL CONREC(OHHR(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € ‘@€]C     CALL FRAME€# ‘€€]C     CALL CONREC(HHOXR(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € ‘À€]C     CALL FRAME€" ’ €]C     CALL CONREC(FNH(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)    € ’@€]C     CALL FRAME€# ’€€]C     CALL CONREC(FNHO2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € ’À€]C     CALL FRAME€"€“ €]C     CALL CONREC(FNOH(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € “@€]C     CALL FRAME€# “€€]C     CALL CONREC(DH2OT(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € “À€]C     CALL FRAME€#€” €]C     CALL CONREC(FNH2O2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B) € ”@€]C     CALL FRAME€ ”€€]C*****************************  € ”À€]C     ****      € • €]C     ****     FLIP BUFFER INDICES      € •@€]C     ****      € •€€]      ITEMP = NJ€€•À€]      NJ = NJNP €	 – €]      NJNP = ITEMP      € –@€]      RETURN    €€–€€]      END       € –À€]CDIR$ NOLIST            €  @€^ªDK ADVECL      €  €€^CDIR$ NOLIST    €€ À€^      SUBROUTINE ADVECL(NX,S,K) €  €^      SAVE      € @€^ªCA PARAMS      €€€€^ªCA ASAS2       €€À€^ªCA CONS2       €€ €^ªCA INDEX       €€@€^      DIMENSION S(IMXP) €
€€€^      NUK=NJ+NU+(K-1)   €
€À€^      NXK=NJ+NX+(K-1)   €
€ €^      DO 1 I=3,IMAXP2   €"€@€^      S(I)=.5*C(52)*(C(10)*(F(I+1,NXK)-F(I-1,NXK))*(F(I+1,NUK)+F(I-1,   €! €€^     ANUK))-C(11)*(F(I+2,NXK)-F(I-2,NXK))*(F(I+2,NUK)+F(I-2,NUK)))      € À€^    1 CONTINUE  €  €^      RETURN    €€@€^      END               €  @€_ªDK CMPN2D      €€ €€_      SUBROUTINE CMPN2D €  À€_      SAVE      €  €_C     ****      €! @€_C     ****     CALCULATE N(N2D) ASSUMING PHOTOCHEMICAL EQUILIBRIUM      € €€_C     ****      € À€_C     ****     COMMON BLOCKS NEEDED:    €!€ €_C     ****       PARAMS, CONS, FIELDS, FIELDT, FLUX, OPION, RATEBLK     € @€_C     ****      € €€_ªCA PARAMS      €€À€_ªCA CONS2       €€ €_ªCA FIELDS2     €€@€_ªCA FIELDT2     €€€€_ªCA FLUX2       €€À€_ªCA NOZNOZ2     €  €_ªCA OPION2      € @€_ªCA RATEBLK2    €
€€€_      DO 1 K=1,KMAXP1   € À€_	DO 1 I=3,LEN1-2€  €_	  FNN2D(K,I)=PPN2D(K,I)/(RBB(2,K,I)*FNO2(K,I)+RBB(4,K,I)*      €$ @€_     1      FNO(K,I)+RBB(5,K,I)*FNE(K,I)+RBB(6,K,I)*FNNO(K,I)+RBB(7,K,I)€€€€_     2      +FIOP(K,I)) € À€_	  FNNOZ(K,I)=FNNO(K,I)+FNNO2(K,I)      €€ €_	  SNO2NO(K,I)=(RBB(12,K,I)*FNHO2(K,I)+RBB(11,K,I)*FNO3(K,I))   €"€@€_     1      /(RBB(14,K,I)+RBB(13,K,I)*FNO(K,I)+RBB(15,K,I)*FNO3(K,I))   € €€_	  RNONOZ(K,I)=1./(1.+SNO2NO(K,I))      € À€_    1 CONTINUE  €  €_C     ****      € @€_C     ****     PERIODIC POINTS  € €€_C     ****      €	 À€_      DO 2 I = 1,2      €	  €_	DO 2 K = 1,KMAXP1      €€@€_	  FNN2D(K,I) = FNN2D(K,I+IMAX) €€€€_	  FNN2D(K,I+IMAXP2) = FNN2D(K,I+2)     €€À€_	  FNNOZ(K,I) = FNNOZ(K,I+IMAX) €€	 €_	  FNNOZ(K,I+IMAXP2) = FNNOZ(K,I+2)     €€	@€_	  SNO2NO(K,I) = SNO2NO(K,I+IMAX)       €€	€€_	  SNO2NO(K,I+IMAXP2) = SNO2NO(K,I+2)   €€	À€_	  RNONOZ(K,I) = RNONOZ(K,I+IMAX)       €€
 €_	  RNONOZ(K,I+IMAXP2) = RNONOZ(K,I+2)   € 
@€_    2 CONTINUE  € 
€€_      RETURN    €€
À€_      END               €€ @€`ªDK COMPCH4     €  €€`      SUBROUTINE COMPCH4€  À€`      SAVE      €  €`C     ****      €€@€`C     ****     ADVANCE CH4 COMPOSITION BY ONE TIME STEP € €€`C     ****      € À€`C     ****     COMMON DECKS:    €"  €`C     ****       PARAMS, ASAS, CH4H2COF, CONS, FIELDS, FIELDT, FLUX,    €€@€`C     ****         INDEX, PHOTO, PSIB, RATEBLKL, VSCR   € €€`C     ****      € À€`ªCA PARAMS      €€ €`ªCA ASAS2       € @€`ªCA CH4H2CO2    €€€€`ªCA CONS2       €€À€`ªCA FIELDS2     €€ €`ªCA FIELDT2     €€@€`ªCA FLUX2       €€€€`ªCA INDEX       € À€`ªCA PHOTO2      €€ €`ªCA PSIB2       € @€`ªCA RATEBLK2    € €€`ªCA VSCR€€À€`      DIMENSION PHICH4(3)       €  €`      DATA RMCH4,PHICH4/16.,0.921,0.846,1.077/  € @€`C     ****      € €€`C     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM    € À€`C     ****      €	€ €`      DO 1 I=1,LEN1     €€@€`	T4(I)=0.       € €€`    1 CONTINUE  € À€`C     ****      €  €`C     ****     SOURCES  € @€`C     ****      €
€€€`      DO 2 K=1,KMAXP1   € À€`	DO 2 I=3,LEN1-2€€	 €`	 S1(I,K)=-(RJCH4(K,I)+RKMM(35,K,I)*FNOH(K,I)+RKMM(36,K,I)*     € 	@€`     1     FNO(K,I)+RKMM(37,K,I)*FNO1D(K,I))    € 	€€`	 S2(I,K)=0.    € 	À€`    2 CONTINUE  €€
 €`C     ****    DENSITY SPECIFIED AT LOWER BOUNDARY       €
€
@€`      DO 4 I=3,LEN1-2   €€
€€`	T1(I)=0.       €€
À€`	T2(I)=1.       €  €`	T3(I)=-XCH4(I)*RMCH4/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+    €€@€`     1    FNN2(1,I)*RMASS(3))   € €€`    4 CONTINUE  € À€`C     ****      €  €`C     ****     PERIODIC POINTS  € @€`C     ****      €	 €€`      DO 3 I = 1,2      €	€À€`	T1(I) = T1(I+IMAX)     €€ €`	T1(I+IMAXP2) = T1(I+2) €	€@€`	T2(I) = T2(I+IMAX)     €€€€`	T2(I+IMAXP2) = T2(I+2) €	€À€`	T3(I) = T3(I+IMAX)     €€ €`	T3(I+IMAXP2) = T3(I+2) €	€@€`	T4(I) = T4(I+IMAX)     €€€€`	T4(I+IMAXP2) = T4(I+2) €
 À€`	  DO 3 K = 1,KMAXP1    €€ €`	    S1(I,K) = S1(I+IMAX,K)     €€@€`	    S1(I+IMAXP2,K) = S1(I+2,K) €€€€`	    S2(I,K) = S2(I+IMAX,K)     €€À€`	    S2(I+IMAXP2,K) = S2(I+2,K) €  €`    3 CONTINUE  €€@€`      IBNDB=0   € €€`      IBND=0    €€À€`      ALFA=0.   € € €`      CALL MINOR(NPCH4,NCH4NM,RMCH4,PHICH4,ALFA,IBND,IBNDB,PCH4B)       € @€`      RETURN    €€€€`      END               €  @€aªDK COMPCO      €€ €€a      SUBROUTINE COMPCO €  À€a      SAVE      €  €aC     ****      € @€aC     ****     ADVANCE CO COMPOSITION BY ONE TIME STEP  € €€aC     ****      € À€aC     ****     COMMON DECKS:    €"  €aC     ****       PARAMS, ASAS, CH4H2COF, CONS, FIELDS, FIELDT, FLUX,    € @€aC     ****         INDEX, PHOTO, PSIB, RATEBLK, VSCR    € €€aC     ****      € À€aªCA PARAMS      €€ €aªCA ASAS2       € @€aªCA CH4H2CO2    €€€€aªCA CONS2       €€À€aªCA FIELDS2     €€ €aªCA FIELDT2     €€@€aªCA FLUX2       €€€€aªCA INDEX       € À€aªCA PHOTO2      €€ €aªCA PSIB2       € @€aªCA RATEBLK2    € €€aªCA VSCR€ À€a      DIMENSION PHICO(3)€  €a      DATA RMCO,PHICO/28.,0.833,1.427,0.852/    € @€aC     ****      € €€aC     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM    € À€aC     ****      €	€ €a      DO 1 I=1,LEN1     €€@€a	T4(I)=0.       € €€a    1 CONTINUE  € À€aC     ****      €  €aC     ****     SOURCES  € @€aC     ****      €
€€€a      DO 2 K=1,KMAXP1   € À€a	DO 2 I=3,LEN1-2€ €	 €a	  S1(I,K)=-(RKMM(40,K,I)*FNO(K,I)*(FNO(K,I)+FNO2(K,I)+FNN2(K,I))       € 	@€a     1      +RKMM(41,K,I)*FNOH(K,I))    €€	€€a	  S2(I,K)=RKCO2P(K,I)*FNCO2(K,I)+(RKMM(35,K,I)*FNOH(K,I)       €#€	À€a     1      +RKMM(36,K,I)*FNO(K,I)+RKMM(37,K,I)*FNO1D(K,I))*FNCH4(K,I)+ €
€
 €a     2      SCOI(K,I)   € 
@€a    2 CONTINUE  €€
€€aC     ****    IF IBNDCO = 0, FLUX GIVEN AT LOWER BOUNDARY       €€
À€a      IF(IBNDCO.EQ.0)THEN       €  €a	DO 3 I=3,LEN1-2€	 @€a          T1(I)=0.      €	 €€a          T2(I)=0.      €	 À€a	  T3(I)=FLUXCO(I)      €  €a    3   CONTINUE€ @€a	IBNDB=1€ €€€aC       ****    IF IBNDCO = 1, PHOTOCHEMICAL EQUILIBRIUM AT LOWER       € À€aC       ****      BOUNDARY      €  €a      ELSE IF(IBNDCO.EQ.1)THEN  € @€a	DO 4 I=3,LEN1-2€	 €€a          T1(I)=0.      €	 À€a          T2(I)=1.      €  €a	  T3(I)=-XCO(I)*RMCO/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+    €€@€a     1      FNN2(1,I)*RMASS(3)) € €€a    4   CONTINUE€ À€aC     ****      €  €aC     ****     PERIODIC POINTS  € @€aC     ****      €	 €€a      DO 5 I = 1,2      €	€À€a	T1(I) = T1(I+IMAX)     €€ €a	T1(I+IMAXP2) = T1(I+2) €	€@€a	T2(I) = T2(I+IMAX)     €€€€a	T2(I+IMAXP2) = T2(I+2) €	€À€a	T3(I) = T3(I+IMAX)     €€ €a	T3(I+IMAXP2) = T3(I+2) €	€@€a	T4(I) = T4(I+IMAX)     €€€€a	T4(I+IMAXP2) = T4(I+2) €
 À€a	  DO 5 K = 1,KMAXP1    €€ €a	    S1(I,K) = S1(I+IMAX,K)     €€@€a	    S1(I+IMAXP2,K) = S1(I+2,K) €€€€a	    S2(I,K) = S2(I+IMAX,K)     €€À€a	    S2(I+IMAXP2,K) = S2(I+2,K) €  €a    5 CONTINUE  €€@€a      IBNDB=0   €€€€a      ENDIF     € À€a      IBND=0    €€ €a      ALFA=0.   €€@€a      CALL MINOR(NPCO,NPCONM,RMCO,PHICO,ALFA,IBND,IBNDB,PCOB)   € €€a      RETURN    €€À€a      END               €€ @€bªDK COMPCO2     €  €€b      SUBROUTINE COMPCO2€  À€b      SAVE      €  €bC     ****      €€@€bC     ****     ADVANCE CO2 COMPOSITION BY ONE TIME STEP € €€bC     ****      € À€bC     ****     COMMON DECKS:    €"  €bC     ****       PARAMS, ASAS, CH4H2COF, CONS, FIELDS, FIELDT, FLUX,    € @€bC     ****         INDEX, PHOTO, PSIB, RATEBLK, VSCR    € €€bC     ****      € À€bªCA PARAMS      €€ €bªCA ASAS2       € @€bªCA CH4H2CO2    €€€€bªCA CONS2       €€À€bªCA FIELDS2     €€ €bªCA FIELDT2     €€@€bªCA FLUX2       €€€€bªCA INDEX       € À€bªCA PHOTO2      €€ €bªCA PSIB2       € @€bªCA RATEBLK2    € €€bªCA VSCR€€À€b      DIMENSION PHICO2(3)       €€ €b      DATA RMCO2,PHICO2/44.,1.199,3.91,1.201/   € @€bC     ****      € €€bC     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM    € À€bC     ****      €	€ €b      DO 1 I=1,LEN1     €€@€b	T4(I)=0.       € €€b    1 CONTINUE  € À€bC     ****      €  €bC     ****     SOURCES  € @€bC     ****      €
€€€b      DO 2 K=1,KMAXP1   €	 À€b	DO 2 I=3,IMAXP4-2      € 	 €b	  S1(I,K)=-(RJCO2T(K,I)+RSKIONS(K,I))  €€	@€b	  S2(I,K)=RKMM(40,K,I)*FNCO(K,I)*FNO(K,I)*(FNO(K,I)+FNO2(K,I)+ €! 	€€b     1      FNN2(K,I))+RKMM(41,K,I)*FNCO(K,I)*FNOH(K,I)+SCO2I(K,I)      € 	À€b    2 CONTINUE  €€
 €bC     ****    DENSITY SPECIFIED AT LOWER BOUNDARY       €
€
@€b      DO 4 I=3,LEN1-2   €€
€€b	T1(I)=0.       €€
À€b	T2(I)=1.       €  €b	T3(I)=-XCO2(I)*RMCO2/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+    €€@€b     1    FNN2(1,I)*RMASS(3))   € €€b    4 CONTINUE  € À€bC     ****      €  €bC     ****     PERIODIC POINTS  € @€bC     ****      €	 €€b      DO 3 I = 1,2      €	€À€b	T1(I) = T1(I+IMAX)     €€ €b	T1(I+IMAXP2) = T1(I+2) €	€@€b	T2(I) = T2(I+IMAX)     €€€€b	T2(I+IMAXP2) = T2(I+2) €	€À€b	T3(I) = T3(I+IMAX)     €€ €b	T3(I+IMAXP2) = T3(I+2) €	€@€b	T4(I) = T4(I+IMAX)     €€€€b	T4(I+IMAXP2) = T4(I+2) €
 À€b	  DO 3 K = 1,KMAXP1    €€ €b	    S1(I,K) = S1(I+IMAX,K)     €€@€b	    S1(I+IMAXP2,K) = S1(I+2,K) €€€€b	    S2(I,K) = S2(I+IMAX,K)     €€À€b	    S2(I+IMAXP2,K) = S2(I+2,K) €  €b    3 CONTINUE  €€@€b      IBNDB=0   € €€b      IBND=0    €€À€b      ALFA=0.   € € €b      CALL MINOR(NPCO2,NCO2NM,RMCO2,PHICO2,ALFA,IBND,IBNDB,PCO2B)       € @€b      RETURN    €€€€b      END               €  @€cªDK COMPH2      €€ €€c      SUBROUTINE COMPH2 €  À€c      SAVE      €  €cC     ****      € @€cC     ****     ADVANCE H2 COMPOSITION BY ONE TIME STEP  € €€cC     ****      € À€cC     ****     COMMON DECKS:    €"  €cC     ****       PARAMS, ASAS, CH4H2COF, CONS, FIELDS, FIELDT, FLUX,    € @€cC     ****         H2FORG, INDEX, PHOTO, PSIB, RATEBLK, VSCR    € €€cC     ****      € À€cªCA PARAMS      €€ €cªCA ASAS2       € @€cªCA CH4H2CO2    €€€€cªCA CONS2       €€À€cªCA FIELDS2     €€ €cªCA FIELDT2     €€@€cªCA FLUX2       €€€€cªCA H2FORG2     €€À€cªCA INDEX       €  €cªCA PHOTO2      €€@€cªCA PSIB2       € €€cªCA RATEBLK2    € À€cªCA VSCR€  €c      DIMENSION PHIH2(3)€€@€c      DATA RMH2,PHIH2/2.,0.226,0.321,0.282/     € €€cC     ****      € À€cC     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM    €  €cC     ****      €	€@€c      DO 1 I=1,LEN1     €€€€c	T4(I)=0.       € À€c    1 CONTINUE  €  €cC     ****      € @€cC     ****     SOURCES  € €€cC     ****      €
€À€c      DO 2 K=1,KMAXP1   €	 	 €c	DO 2 I=3,IMAXP4-2      €€	@€c	  S1(I,K)=-(RKMM(11,K,I)*FNO1D(K,I)+RKMM(20,K,I)*FNO(K,I)+     € 	€€c     1      RKMM(25,K,I)*FNOH(K,I)+FK9OP(K,I))  € 	À€c	  S2(I,K)=0.12*RJH2OLY(K,I)*FNH2O(K,I)+RKMM(30,K,I)*FNH(K,I)*  € 
 €c     1      FNHO2(K,I)+RKMM(33,K,I)*    €$ 
@€c     2      FNH(K,I)*FNH(K,I)*(FNO(K,I)+FNO2(K,I)+FNN2(K,I))+RJCH4(K,I)*€ 
€€c     3      FNCH4(K,I)  € 
À€c    2 CONTINUE  €€ €cC     ****    DENSITY SPECIFIED AT LOWER BOUNDARY       €
€@€c      DO 4 I=3,LEN1-2   €€€€c	T1(I)=0.       €€À€c	T2(I)=1.       €  €c	T3(I)=-XH2(I)*RMH2/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+      €€@€c     1    FNN2(1,I)*RMASS(3))   € €€c    4 CONTINUE  € À€cC     ****      €  €cC     ****     PERIODIC POINTS  € @€cC     ****      €	 €€c      DO 3 I = 1,2      €	€À€c	T1(I) = T1(I+IMAX)     €€ €c	T1(I+IMAXP2) = T1(I+2) €	€@€c	T2(I) = T2(I+IMAX)     €€€€c	T2(I+IMAXP2) = T2(I+2) €	€À€c	T3(I) = T3(I+IMAX)     €€ €c	T3(I+IMAXP2) = T3(I+2) €	€@€c	T4(I) = T4(I+IMAX)     €€€€c	T4(I+IMAXP2) = T4(I+2) €
 À€c	  DO 3 K = 1,KMAXP1    €€ €c	    S1(I,K) = S1(I+IMAX,K)     €€@€c	    S1(I+IMAXP2,K) = S1(I+2,K) €€€€c	    S2(I,K) = S2(I+IMAX,K)     €€À€c	    S2(I+IMAXP2,K) = S2(I+2,K) €  €c    3 CONTINUE  €€@€c      IBNDB=0   € €€c      IBND=0    €€À€c      ALFA=0.   €€ €c      CALL MINOR(NPH2,NPH2NM,RMH2,PHIH2,ALFA,IBND,IBNDB,PH2B)   € @€c      RETURN    €€€€c      END               €€ @€dªDK COMPH2O     €  €€d      SUBROUTINE COMPH2O€  À€d      SAVE      €  €dC     ****      €€@€dC     ****     ADVANCE H2O COMPOSITION BY ONE TIME STEP € €€dC     ****      € À€dC     ****     COMMON DECKS:    €!€ €dC     ****       PARAMS, ASAS, CONS, DISSHOX, FIELDS, FIELDT, FLUX,     € @€dC     ****         INDEX, IONHOX, PHOTO, PSIB, RATEBLK, VSCR    € €€dC     ****      € À€dC     ****      €  €dªCA PARAMS      €€@€dªCA ASAS2       €€€€dªCA CONS2       € À€dªCA DISSHOX2    €€ €dªCA FIELDS2     €€@€dªCA FIELDT2     €€€€dªCA FLUX2       €€À€dªCA INDEX       €€ €dªCA IONHOX2     € @€dªCA PHOTO2      €€€€dªCA PSIB2       € À€dªCA RATEBLK2    €  €dªCA VSCR€€@€d      DIMENSION PHIH2O(3)       € €€d      DATA RMH2O,PHIH2O/18.,0.817,0.922,0.920/  € À€dC     ****      €  €dC     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM    € @€dC     ****      €	€€€d      DO 1 I=1,LEN1     €€À€d	T4(I)=0.       €  €d    1 CONTINUE  € @€dC     ****      € €€dC     ****     SOURCES  € À€dC     ****      €
€	 €d      DO 2 K=1,KMAXP1   € 	@€d	DO 2 I=3,LEN1-1€  	€€d	  S1(I,K)=-(SH2OSRB(K,I)+SH2OSRC(K,I)+SH2OLYA(K,I)+SH2OEUV(K,I)€$ 	À€d     1      +RKMM(10,K,I)*FNO1D(K,I)+RKK(10,K,I)*YIOP(K,I)+RKMM(42,K,I)*€
€
 €d     2      FNO(K,I))   €€
@€d	  S2(I,K)=RKMM(22,K,I)*FNOH(K,I)*FNOH(K,I)+RKMM(23,K,I)*       €" 
€€d     1      FNOH(K,I)*FNHO2(K,I)+RKMM(24,K,I)*FNOH(K,I)*FNH2O2(K,I)+    €!€
À€d     2      RKMM(25,K,I)*FNOH(K,I)*FNH2(K,I)+RKMM(32,K,I)*FNH(K,I)*     €$  €d     3      FNHO2(K,I)+(2.*RKMM(35,K,I)*FNOH(K,I)+RKMM(36,K,I)*FNO(K,I)+€€@€d     4      RKMM(37,K,I)*FNO1D(K,I))*FNCH4(K,I)-PHOXIC(K,I)*0.5 € €€d    2 CONTINUE  € À€dC     ****    DENSITY SPECIFIED AT LOWER BOUNDAR€
€ €d      DO 4 I=3,LEN1-2   €€@€d	T1(I)=0.       €€€€d	T2(I)=1.       € À€d	T3(I)=-XH2O(I)*RMH2O/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+    €€ €d     1    FNN2(1,I)*RMASS(3))   € @€d    4 CONTINUE  € €€dC     ****      € À€dC     ****     PERIODIC POINTS  €  €dC     ****      €	 @€d      DO 3 I = 1,2      €	€€€d	T1(I) = T1(I+IMAX)     €€À€d	T1(I+IMAXP2) = T1(I+2) €	€ €d	T2(I) = T2(I+IMAX)     €€@€d	T2(I+IMAXP2) = T2(I+2) €	€€€d	T3(I) = T3(I+IMAX)     €€À€d	T3(I+IMAXP2) = T3(I+2) €	€ €d	T4(I) = T4(I+IMAX)     €€@€d	T4(I+IMAXP2) = T4(I+2) €
 €€d	  DO 3 K = 1,KMAXP1    €€À€d	    S1(I,K) = S1(I+IMAX,K)     €€ €d	    S1(I+IMAXP2,K) = S1(I+2,K) €€@€d	    S2(I,K) = S2(I+IMAX,K)     €€€€d	    S2(I+IMAXP2,K) = S2(I+2,K) € À€d    3 CONTINUE  €$  €dC     CALL CONREC(SH2OSRB(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)€ @€dC     CALL FRAME€$ €€dC     CALL CONREC(SH2OSRC(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)€ À€dC     CALL FRAME€$  €dC     CALL CONREC(SH2OLYA(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)€ @€dC     CALL FRAME€$ €€dC     CALL CONREC(SH2OEUV(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)€ À€dC     CALL FRAME€#  €dC     CALL CONREC(FNO1D(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € @€dC     CALL FRAME€"€€€dC     CALL CONREC(YIOP(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € À€dC     CALL FRAME€"  €dC     CALL CONREC(FNO(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)    € @€dC     CALL FRAME€"€€€dC     CALL CONREC(FNOH(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € À€dC     CALL FRAME€#  €dC     CALL CONREC(FNHO2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € @€dC     CALL FRAME€#€€€dC     CALL CONREC(FNH2O2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B) € À€dC     CALL FRAME€"€ €dC     CALL CONREC(FNH2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € @€dC     CALL FRAME€" €€dC     CALL CONREC(FNH(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)    € À€dC     CALL FRAME€#  €dC     CALL CONREC(FNCH4(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)  € @€dC     CALL FRAME€#€€€dC     CALL CONREC(PHOXIC(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B) € À€dC     CALL FRAME€"€ €dC     CALL CONREC(FNO2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € @€dC     CALL FRAME€"€€€dC     CALL CONREC(FNN2(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)   € À€dC     CALL FRAME€  €dC     CALL CONREC(S1,LEN1,LEN1,KMAXP1,0.,0.,0.,0,0,-1430B)      € @€dC     CALL FRAME€ €€dC     CALL CONREC(S2,LEN1,LEN1,KMAXP1,0.,0.,0.,0,0,-1430B)      € À€dC     CALL FRAME€€ €d      DO 5 I = 3,IMAX+3 €	 @€d	DO 5 K = 1,KMAXP1      € €€d	  S3(I,K) = RKK(10,K,I)€€À€d	  S4(I,K) = RKMM(10,K,I)       €€ €d	  S5(I,K) = RKMM(22,K,I)       €€@€d	  S6(I,K) = RKMM(23,K,I)       €€€€d	  S7(I,K) = RKMM(24,K,I)       €€À€d	  S8(I,K) = RKMM(25,K,I)       €€ €d	  S9(I,K) = RKMM(32,K,I)       € @€d	  S10(I,K) = RKMM(35,K,I)      € €€d	  S11(I,K) = RKMM(36,K,I)      € À€d	  S12(I,K) = RKMM(37,K,I)      €  €d	  S13(I,K) = RKMM(42,K,I)      € @€d    5 CONTINUE  € €€€dC     CALL CONREC(S3(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)       € À€dC     CALL FRAME€ € €dC     CALL CONREC(S4(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)       € @€dC     CALL FRAME€ €€€dC     CALL CONREC(S5(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)       € À€dC     CALL FRAME€ €  €dC     CALL CONREC(S6(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)       €  @€dC     CALL FRAME€ € €€dC     CALL CONREC(S7(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)       €  À€dC     CALL FRAME€ €! €dC     CALL CONREC(S8(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)       € !@€dC     CALL FRAME€ €!€€dC     CALL CONREC(S9(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)       € !À€dC     CALL FRAME€! " €dC     CALL CONREC(S10(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)      € "@€dC     CALL FRAME€! "€€dC     CALL CONREC(S11(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)      € "À€dC     CALL FRAME€! # €dC     CALL CONREC(S12(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)      € #@€dC     CALL FRAME€! #€€dC     CALL CONREC(S13(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)      € #À€dC     CALL FRAME€€$ €dC     WRITE(6,2345)(T3(I),XH2O(I),I=1,LEN1)     €
 $@€dC2345 FORMAT(2E15.4)    €€$€€d      IBNDB=0   € $À€d      IBND=0    €€% €d      ALFA=0.   € €%@€d      CALL MINOR(NPH2O,NH2ONM,RMH2O,PHIH2O,ALFA,IBND,IBNDB,PH2OB)       € %€€d      RETURN    €€%À€d      END               €  @€eªDK COMPHE      €€ €€e      SUBROUTINE COMPHE €  À€e      SAVE      €  €eC     ****      € @€eC     ****     COMMON DECKS:    €!€€€eC     ****       PARAMS, ASAS, CONS, FIELDS, HEFLUX, INDEX, PHOTOO,     €€À€eC     ****         PSIB, VSCR   €  €eC     ****      € @€eªCA PARAMS      €€€€eªCA ASAS2       €€À€eªCA CONS2       €€ €eªCA FIELDS2     €€@€eªCA HEFLUX2     €€€€eªCA INDEX       €€À€eªCA PHOTOO2     €€ €eªCA PSIB2       € @€eªCA VSCR€€€€e      COMMON/SODIUMP/PNA(KMXP,IMXP),XLNA(KMXP,IMXP),XLBNA(IMXP) € À€e     1,PNAI(KMXP,IMXP),XLNAI(KMXP,IMXP),XLBNAI(IMXP)    €  €e     2,RNAO(KMXP,IMXP),RNAO2(KMXP,IMXP),RNAOH(KMXP,IMXP)€€@€e      DIMENSION PHIH(3) €€€€e      DIMENSION FMM(IMXP)       €  À€eC     DATA RMH,PHIH,ALFAH,PSHB/23.,.270,.404,.322,-0.38,0.73E-6/€!€ €e      DATA RMH,PHIH,ALFAH,PSHB/23.,1.042,1.509,1.176,0.17,0.131E-1/     €€@€eC     ****     UPPER BOUNDARY CONDITION € €€eC     ****     T4=TOTAL FLUX AT UPPER BOUNDARY  € À€e      DO 10 I=3,LEN1-2  €€ €e      FMM(I)=FNO(1,I)*16.+FNO2(1,I)*32.+FNN2(1,I)*28.   €€@€e      T4(I) =0. € €€e      T1(I)=0.  € À€e      T2(I)=1.  €€ €e      T3(I)=-XLBNAI(I)*RMH/FMM(I)       € @€e   10 CONTINUE  € €€eC     DO 10 I=3,LEN1-2  €€À€eC     FMM(I)=FNO(1,I)*16.+FNO2(1,I)*32.+FNN2(1,I)*28.   €
€	 €eC     T4(I) = FLUX(I)   € 	@€eC  10 CONTINUE  € 	€€eC     ****     LOWER BOUNDARY CONDITION - CONSTANT VALUE€
€	À€eC     DO 4 I=3,LEN1-2   € 
 €eC     T1(I)=0.  € 
@€eC     T2(I)=1.  € 
€€eC     T3(I)=-XHE(I)*RMH/FMM(I)  € 
À€eC   4 CONTINUE  €€ €eC     ****     SOURCES ARE ZERO €
€@€e      DO 5 K=1,KMAXP1   €€€€e      DO 5 I=3,IMAXP4-2 €€À€e      S1(I,K)=-XLNAI(K,I)       €€ €e      S2(I,K)=PNAI(K,I) € @€e    5 CONTINUE  € €€eC     ****      € À€eC     ****     PERIODIC POINTS  €  €eC     ****      €	 @€e      DO 3 I = 1,2      €	€€€e	T1(I) = T1(I+IMAX)     €€À€e	T1(I+IMAXP2) = T1(I+2) €	€ €e	T2(I) = T2(I+IMAX)     €€@€e	T2(I+IMAXP2) = T2(I+2) €	€€€e	T3(I) = T3(I+IMAX)     €€À€e	T3(I+IMAXP2) = T3(I+2) €	€ €e	T4(I) = T4(I+IMAX)     €€@€e	T4(I+IMAXP2) = T4(I+2) € €€e    3 CONTINUE  € À€e      IBND=0    €€ €e      IBNDB=0   €€@€e      CALL MINOR(NPSH,NPSHNM,RMH,PHIH,ALFAH,IBND,IBNDB,PHB)     € €€e      RETURN    €€À€e      END               €€ @€fªDK COMPARG     €  €€f      SUBROUTINE COMPARG€  À€f      SAVE      €  €fC     ****     ADVANCE AR COMPOSITION BY ONE TIME STEP  € @€fC     ****      € €€fC     ****   COMMON DECKS:      €€À€fC     ****   PARAMS,ASAS,CONS,FIELDS,INDEX,PHOTOO,PSIB,VSCR     €  €fC     ****      € @€fªCA PARAMS      €€€€fªCA ASAS2       €€À€fªCA CONS2       €€ €fªCA FIELDS2     €€@€fªCA INDEX       €€€€fªCA PHOTOO2     €€À€fªCA PSIB2       €  €fªCA VSCR€€@€f      COMMON/SODIUMP/PNA(KMXP,IMXP),XLNA(KMXP,IMXP),XLBNA(IMXP) € €€f     1,PNAI(KMXP,IMXP),XLNAI(KMXP,IMXP),XLBNAI(IMXP)    € À€f     2,RNAO(KMXP,IMXP),RNAO2(KMXP,IMXP),RNAOH(KMXP,IMXP)€€ €f      DIMENSION PHIA(3) €€@€f      DIMENSION FMM(IMXP)       €!€€€f      DATA RMA,PHIA,ALFAA,PSAB/23.,1.042,1.509,1.176,0.17,0.131E-1/     €€À€fC     ****     UPPER BOUNDARY - ZERO DIFFUSIVE FLUX     €  €fC     ****     LOWER BOUNDARY - CONSTANT VALUE  €
€@€f      DO 1 I=3,LEN1-2   €€€€f      FMM(I)=FNO(1,I)*16.+FNO2(1,I)*32.+FNN2(1,I)*28.   € À€f      T4(I)=0.  €€ €f      T1(I)=0   € @€f      T2(I)=1.  € €€f      T3(I)=-XLBNA(I)*RMA/FMM(I)€ À€f    1 CONTINUE  €€ €f      IBNDB=0   €€@€fC     LOWER BOUNDARY FLUX CONDITION     €€€€fC     FLUXNA=0. €
€À€fC     DO 1 I=3,LEN1-2   € 	 €fC     T4(I)=0.  €€	@€fC     T1(I)=0   € 	€€fC     T2(I)=0.  €	 	À€fC     T3(I)=FLUXNA      € 
 €fC   1 CONTINUE  €€
@€fC     IBNDB=1   € 
€€fC     ****     SOURCES  €
€
À€f      DO 2 K=1,KMAXP1   €	  €f	DO 2 I=3,IMAXP4-2      €
€@€f	  S1(I,K)=-XLNA(K,I)   €	€€€f	  S2(I,K)=PNA(K,I)     € À€f    2 CONTINUE  €	  €f      DO 3 I = 1,2      €	€@€f	T1(I) = T1(I+IMAX)     €€€€f	T1(I+IMAXP2) = T1(I+2) €	€À€f	T2(I) = T2(I+IMAX)     €€ €f	T2(I+IMAXP2) = T2(I+2) €	€@€f	T3(I) = T3(I+IMAX)     €€€€f	T3(I+IMAXP2) = T3(I+2) €	€À€f	T4(I) = T4(I+IMAX)     €€ €f	T4(I+IMAXP2) = T4(I+2) €
 @€f	  DO 3 K = 1,KMAXP1    €€€€f	    S1(I,K) = S1(I+IMAX,K)     €€À€f	    S1(I+IMAXP2,K) = S1(I+2,K) €€ €f	    S2(I,K) = S2(I+IMAX,K)     €€@€f	    S2(I+IMAXP2,K) = S2(I+2,K) € €€f    3 CONTINUE  € À€f      IBND=0    €€ €f      CALL MINOR (NPSA,NPSANM,RMA,PHIA,ALFAA,IBND,IBNDB,PARB)   € @€f      RETURN    €€€€f      END               €€ @€gªDK COMPHOX     €  €€g      SUBROUTINE COMPHOX€  À€g      SAVE      €  €gC     ****      €€@€gC     ****     ADVANCE HOX COMPOSITION BY ONE TIME STEP € €€gC     ****      € À€gC     ****     COMMON DECKS:    €!€ €gC     ****       PARAMS, ASAS, CONS, DISSHOX, FIELDS, FIELDT, FLUX,     €#€@€gC     ****         HOXSTUF, HOXUPP, INDEX, IONHOX, NEWRAT, PHOTO, PSIB, €€€€gC     ****         RATEBLK, RJH2O2O, VSCR       € À€gC     ****      €  €gªCA PARAMS      €€@€gªCA ASAS2       €€€€gªCA CONS2       € À€gªCA DISSHOX2    €€ €gªCA FIELDS2     €€@€gªCA FIELDT2     €€€€gªCA FLUX2       € À€gªCA HOXSTUF2    €€ €gªCA HOXUPP2     €€@€gªCA INDEX       €€€€gªCA IONHOX2     €€À€gªCA NEWRAT2     €  €gªCA PHOTO2      €€@€gªCA PSIB2       € €€gªCA RATEBLK2    € À€gªCA RJH2O2O2    €  €gªCA VSCR€  @€g      DIMENSION FM(IMXP),SS1(IMXP),SS2(IMXP),RR(IMXP),XHOX(IMXP)€€€€g      DIMENSION PHIHOX(3)       €€À€g      DATA RMHOX,PHIHOX/1.,0.146,0.243,0.162/   €  €gC     ****      € @€gC     ****     UPPER BOUNDARY: TOTAL FLUX=FHOX  € €€gC     ****      €
€À€g      DO 1 I=3,LEN1-2   € 	 €g	T4(I)=FHOX(I)  € 	@€g    1 CONTINUE  € 	€€g      IBND=1    € 	À€gC     ****      € 
 €gC     ****     SOURCES  € 
@€gC     ****      € 
€€gC     U=9.      €€
À€gC     V=17.     €€ €gC     W=25.     € @€g      U=1.      € €€g      V=1.      € À€g      W=1.      €
€ €g      DO 2 K=1,KMAXP1   € @€g	DO 2 I=3,LEN1-2€€€€g	  FM(I)=FNO(K,I)+FNO2(K,I)+FNN2(K,I)   €€À€g	  HO2HR(K,I)=RKMM(28,K,I)*FM(I)*FNO2(K,I)/(RKMM(18,K,I)*       €!  €g     1      FNO(K,I)+RBB(12,K,I)*FNNO(K,I)+RKMM(26,K,I)*FNO3(K,I))      € @€g	  OHHR(K,I)=((RKMM(18,K,I)*FNO(K,I)*RKMM(28,K,I)*FNO2(K,I)*    € €€€g     1      FM(I))/(RKMM(18,K,I)*FNO(K,I)+RKMM(26,K,I)*FNO3(K,I)+       €€À€g     2      RBB(12,K,I)*FNNO(K,I))+RKMM(29,K,I)*FNO3(K,I))/     €  €g     3      (RKMM(17,K,I)*FNO(K,I)+RKMM(21,K,I)*FNO3(K,I))      € @€g	  HHOXR(K,I)=1./(1.+HO2HR(K,I)+OHHR(K,I))      €
 €€g	  SS1(I)=HO2HR(K,I)    €	€À€g	  SS2(I)=OHHR(K,I)     €	€ €g	  RR(I)=HHOXR(K,I)     € @€g	  S2(I,K)=(2.*(SH2OSRB(K,I)+SH2OSRC(K,I)+SH2OEUV(K,I)+0.88*    €" €€g     1      SH2OLYA(K,I))*FNH2O(K,I)*U+2.*RJH2O2(K,I)*FNH2O2(K,I)*V+    €"€À€g     2      PHOXIC(K,I)+(2.*RKMM(11,K,I)*FNH2(K,I)*U+2.*RKMM(10,K,I)*   €!  €g     3      FNH2O(K,I)*V+2.*RKMM(37,K,I)*FNCH4(K,I))*FNO1D(K,I)*V+      €# @€g     4      RKK(11,K,I)*YIHP(K,I)*FNO(K,I)+2.*RKMM(20,K,I)* FNH2(K,I)*  € €€€g     5      FNO(K,I)*U+2.*RKMM(19,K,I)*FNH2O2(K,I)*FNO(K,I)*W+2.*       €€À€g     6      RKMM(42,K,I)*FNH2O(K,I)*FNO(K,I)*V+2.*RKMM(36,K,I)* € € €g     7      FNCH4(K,I)*FNO(K,I)*V+RKK(9,K,I)*FNH2(K,I)*YIOP(K,I)+       € @€g     8      ALPP(4,K,I)*YIOHP(K,I)*FNE(K,I))    €€€€g	  S1(I,K)=-(2.*RKMM(23,K,I)*SS1(I)*RR(I)*SS2(I)*RR(I)+ €#€À€g     1      2.*RKMM(27,K,I)*SS1(I)*SS1(I)*RR(I)*RR(I)+2.*(RKMM(30,K,I)+ €"  €g     2      RKMM(32,K,I))*SS1(I)*RR(I)*RR(I)+2.*RKMM(22,K,I)*SS2(I)*    €! @€g     3      RR(I)*SS2(I)*RR(I)+2.*RKMM(33,K,I)*FM(I)*RR(I)*RR(I))*      € €€g     4      FNHOX(K,I)-RKK(12,K,I)*YIOP(K,I)*RR(I)      € À€g    2 CONTINUE  €€ €gC     ****    PHOTOCHEMICAL EQUILIBRIUM AT LOWER BOUNDARY       €€@€g      DO 3 I = 3,LEN1-2 €€€€g	FM(I)=FNO(1,I)+FNO2(1,I)+FNN2(1,I)     €	 À€g	SS1(I)=HO2HR(1,I)      €€ €g	SS2(I)=OHHR(1,I)       €€@€g	RR(I)=HHOXR(1,I)       €€€€g	XHOX(I)=SQRT((2.*(SH2OSRB(1,I)+SH2OSRC(1,I)+SH2OEUV(1,I)+0.88* €! À€g     1    SH2OLYA(1,I))*FNH2O(1,I)*U+2.*RJH2O2(1,I)*FNH2O2(1,I)*V+      €
€ €g     1    PHOXIC(1,I)   €$ @€g     2    +(2.*RKMM(11,1,I)*FNH2(1,I)*U+2.*RKMM(10,1,I)*FNH2O(1,I)*V+2.*€# €€g     3    RKMM(37,1,I)*FNCH4(1,I))*FNO1D(1,I)*V+RKK(11,1,I)*YIHP(1,I)*  €$ À€g     4    FNO(1,I)+2.*RKMM(20,1,I)*FNH2(1,I)*FNO(1,I)*U+2.*RKMM(19,1,I)*€#€ €g     5    FNH2O2(1,I)*FNO(1,I)*W+2.*RKMM(42,1,I)*FNH2O(1,I)*FNO(1,I)*V+ €"€@€g     6    2.*RKMM(36,1,I)*FNCH4(1,I)*FNO(1,I)*V+RKK(9,1,I)*FNH2(1,I)*   €# €€g     7    YIOP(1,I)+ALPP(4,1,I)*YIOHP(1,I)*FNE(1,I))/(2.*RKMM(23,1,I)*  €$ À€g     8    SS1(I)*RR(I)*SS2(I)*RR(I)+2.*RKMM(27,1,I)*SS1(I)*SS1(I)*RR(I)*€"€ €g     9    RR(I)+2.*(RKMM(30,1,I)+RKMM(32,1,I))*SS1(I)*RR(I)*RR(I)+2.*   €#€@€g     1    RKMM(22,1,I)*RR(I)*RR(I)*SS2(I)*SS2(I)+2.*RKMM(33,1,I)*FM(I)* €€€€g     2    RR(I)*RR(I))) € À€g    3 CONTINUE  €
€ €g      DO 4 I=3,LEN1-2   €€@€g	T1(I)=0.       €€€€g	T2(I)=1.       € À€g	T3(I)=-XHOX(I)/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+  €€ €g     1    FNN2(1,I)*RMASS(3))   € @€g    4 CONTINUE  € €€gC     ****      € À€gC     ****     PERIODIC POINTS  €  €gC     ****      €	 @€g      DO 5 I = 1,2      €	€€€g	T1(I) = T1(I+IMAX)     €€À€g	T1(I+IMAXP2) = T1(I+2) €	€ €g	T2(I) = T2(I+IMAX)     €€@€g	T2(I+IMAXP2) = T2(I+2) €	€€€g	T3(I) = T3(I+IMAX)     €€À€g	T3(I+IMAXP2) = T3(I+2) €	€ €g	T4(I) = T4(I+IMAX)     €€@€g	T4(I+IMAXP2) = T4(I+2) €
 €€g	  DO 5 K = 1,KMAXP1    €€À€g	    S1(I,K) = S1(I+IMAX,K)     €€ €g	    S1(I+IMAXP2,K) = S1(I+2,K) €€@€g	    S2(I,K) = S2(I+IMAX,K)     €€€€g	    S2(I+IMAXP2,K) = S2(I+2,K) € À€g    5 CONTINUE  €€ €g      IBNDB=0   € @€g      ALFA=-0.38€ €€€g      CALL MINOR(NPHOX,NHOXNM,RMHOX,PHIHOX,ALFA,IBND,IBNDB,PHOXB)       € À€g      RETURN    €€ €g      END               €  @€hªDK COMPNO      €€ €€h      SUBROUTINE COMPNO €  À€h      SAVE      €  €hC     ****      € @€hC     ****     ADVANCE NO COMPOSITION BY ONE TIME STEP  € €€hC     ****      € À€hC     ****     COMMON DECKS:    €$  €hC     ****       PARAMS, ASAS, CONS, FIELDS, FIELDT, FLUX, INDEX, PHOTO,€€@€hC     ****          PSIB, RATEBLK, VSCR € €€hC     ****      € À€hªCA PARAMS      €€ €hªCA ASAS2       €€@€hªCA CONS2       €€€€hªCA FIELDS2     €€À€hªCA FIELDT2     €€ €hªCA FLUX2       €€@€hªCA INDEX       €€€€hªCA NOZNOZ2     € À€hªCA PHOTO2      €€ €hªCA PSIB2       € @€hªCA RATEBLK2    € €€hªCA VSCR€ À€h      DIMENSION PHINO(3)€  €h      DATA RMNO,PHINO/30.,0.814,0.866,0.926/    € @€hC     ****      € €€hC     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM    € À€hC     ****      €	€ €h      DO 1 I=1,LEN1     €€@€h	T4(I)=0.       € €€h    1 CONTINUE  € À€hC     ****      €  €hC     ****     SOURCES  € @€hC     ****      €
€€€h      DO 2 K=1,KMAXP1   € À€h	DO 2 I=3,LEN1-2€ 	 €h	  S1(I,K)=-(FK5O2P(K,I)+RBB(3,K,I)*FNN4S(K,I)+RBB(6,K,I)*      €€	@€h     1      FNN2D(K,I)+RBB(8,K,I)+RBB(9,K,I))*RNONOZ(K,I)       €  	€€h	  S2(I,K)=RBB(1,K,I)*FNO2(K,I)*FNN4S(K,I)+RBB(2,K,I)*FNO2(K,I)*€€	À€h     1      FNN2D(K,I)+RBB(10,K,I)*FNOH(K,I)*FNN4S(K,I) € 
 €h    2 CONTINUE  €€
@€hC     ****    IF IBNDNO = 0, FLUX GIVEN AT LOWER BOUNDARY       €€
€€h      IF(IBNDNO.EQ.0)THEN       € 
À€h	DO 3 I=3,LEN1-2€	  €h          T1(I)=0.      €	 @€h          T2(I)=0.      €	 €€h	  T3(I)=FLUXNO(I)      € À€h    3   CONTINUE€  €h	IBNDB=1€ €@€hC       ****    IF IBNDNO = 1, PHOTOCHEMICAL EQUILIBRIUM AT LOWER       € €€hC       ****      BOUNDARY      € À€h      ELSE IF(IBNDNO.EQ.1)THEN  €  €h	DO 4 I=3,LEN1-2€	 @€h          T4(I)=0.      €	 €€h          T1(I)=0.      €	 À€h          T2(I)=1.      €  €h	  T3(I)=-XNO(I)*RMNO/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+    €€@€h     1    FNN2(1,I)*RMASS(3))   € €€h    4   CONTINUE€ À€h	IBNDB=0€€ €h      ENDIF     €	 @€h      DO 5 I = 1,2      €	€€€h	T1(I) = T1(I+IMAX)     €€À€h	T1(I+IMAXP2) = T1(I+2) €	€ €h	T2(I) = T2(I+IMAX)     €€@€h	T2(I+IMAXP2) = T2(I+2) €	€€€h	T3(I) = T3(I+IMAX)     €€À€h	T3(I+IMAXP2) = T3(I+2) €	€ €h	T4(I) = T4(I+IMAX)     €€@€h	T4(I+IMAXP2) = T4(I+2) €	 €€h	DO 5 K = 1,KMAXP1      €€À€h	  S1(I,K) = S1(I+IMAX,K)       €€ €h	  S1(I+IMAXP2,K) = S1(I+2,K)   €€@€h	  S2(I,K) = S2(I+IMAX,K)       €€€€h	  S2(I+IMAXP2,K) = S2(I+2,K)   € À€h    5 CONTINUE  €  €h      IBND=0    €€@€h      ALFA=0.   €€€€h      CALL MINOR(NPNO,NPNONM,RMNO,PHINO,ALFA,IBND,IBNDB,PNOB)   € À€h      RETURN    €€ €h      END               €€ @€iªDK FFACE       €  €€i      SUBROUTINE FFACE(II)      €  À€i      SAVE      €  €iC     ****      €  @€iC     ****     INTERFACES ONE- AND TWO-DIMENSIONAL COMMON BLOCKS€ €€iC     ****      € À€iªCA PARAMS      €€ €iªCA CONS2       €€@€iªCA ASASF       € €€iªCA ASAS2F      € À€iªCA CH4H2CFF    €  €iªCA CH4H2C2F    €€@€iªCA CRATESF     € €€iªCA CRATES2F    € À€iªCA DISSHOXF    €  €iªCA DISSHX2F    €€@€iªCA FIELDSF     € €€iªCA FIELDS2F    €€À€iªCA FIELDTF     €  €iªCA FIELDT2F    € @€iªCA FLUX€ €€iªCA FLUX2F      € À€iªCA H2FORG      €  €iªCA H2FORG2F    € @€iªCA HEFLUX      € €€iªCA HEFLUX2F    € À€iªCA HOXSTUFF    €  €iªCA HOXSTF2F    € @€iªCA HOXUPP      € €€iªCA HOXUPP2F    € À€iªCA IONHOX      €  €iªCA IONHOX2F    €€@€iªCA NEWRATF     € €€iªCA NEWRAT2F    €€À€iªCA OPION       €€	 €iªCA OPION2F     €€	@€iªCA PSIBF       € 	€€iªCA PSIB2F      € 	À€iªCA PHOTOF      €€
 €iªCA PHOTO2F     €€
@€iªCA PHOTOOF     € 
€€iªCA PHOTOO2F    €€
À€iªCA RATEBLK     €  €iªCA RATBLK2F    €€@€iªCA RJH2O2O     € €€iªCA RJH2O22F    €€À€iªCA SOURCEF     €  €iªCA SOURCE2F    € @€iC     ****      €# €€iC     ****     ALLOW FOR PERIODIC WRAP AROUND POINTS (TWO AT EACH END)  € À€iC     ****      €  €i      I = II+2  € @€iC     ****      € €€€iC     ****     /ASAS/ AND /ASAS2/ ***** NEED TO WORK OUT LEAPFROG       €€À€iC     ****       POINTERS       €  €iC     ****      € @€iC     ****     /CH4H2COF/ AND /CH4H2CO2/€ €€iC     ****      €	 À€i      DO 2 J = 1,7      €	  €i	DO 2 K = 1,KMAXP1      €€@€i	  RJCH4F(K,I,J) = RJCH4(K,J)   € €€i    2 CONTINUE  € À€iC     ****      €€ €iC     ****     /CRATES/ AND /CRATES2/   € @€iC     ****      €	 €€i      DO 3 J = 1,4      €	 À€i	DO 3 K = 1,KMAXP1      €€ €i	  RK12F(K,I,J) = RK12(K,J)     € @€i    3 CONTINUE  € €€iC     ****      €€À€iC     ****     /DISSHOX/ AND /DISSHOX2/ €  €iC     ****      €	 @€i      DO 4 J = 1,9      €	 €€i	DO 4 K = 1,KMAXP1      €€À€i	  SH2OTF(K,I,J) = SH2OT(K,J)   €  €i    4 CONTINUE  € @€iC     ****      €€€€iC     ****     /FIELDS/ AND /FIELDS2/   € À€iC     ****      €	€ €i      DO 5 J = 1,17     €	 @€i	DO 5 K = 1,KMAXP1      €€€€i	  FTF(K,I,J) = FT(K,J) € À€i    5 CONTINUE  €	  €i      F107F = F107      €	 @€i      DO 6 J = 1,3      €	 €€i	DO 6 K = 1,KMAXP1      €€À€i	  FTEF(K,I,J) = FTE(K,J)       €  €i    6 CONTINUE  € @€iC     ****      €€€€iC     ****     /FIELDT/ AND /FIELDT2/   € À€iC     ****      €	€ €i      DO 7 J = 1,12     €	 @€i	DO 7 K = 1,KMAXP1      €€€€i	  FNH2OF(K,I,J) = FNH2O(K,J)   € À€i    7 CONTINUE  €  €iC     ****      €€@€iC     ****     /FLUX/ AND /FLUX2/       € €€iC     ****      €€À€i      FLUXNOF(I) = FLUXNO       €  €i      IBNDNOF = IBNDNO  €€@€i      FLUXCOF(I) = FLUXCO       € €€i      IBNDCOF = IBNDCO  € À€iC     ****      €€ €iC     ****     /H2FORG/ AND /H2FORG2/   € @€iC     ****      € €€i      DO 8 K =1,KMAXP1  €€À€i	RJH2OLYF(K,I) = RJH2OLY(K)     €  €i    8 CONTINUE  € @€iC     ****      €€€€iC     ****     /HEFLUX/ AND /HEFLUX2/   € À€iC     ****      €
€ €i      FLUXF(I) = FLUX   € @€iC     ****      €€€€iC     ****     /HOXSTUF/ AND /HOXSTUF2/ € À€iC     ****      €	€ €i      DO 9 J = 1,10     €	 @€i	DO 9 K = 1,KMAXP1      €€€€i	  HO2HRF(K,I,J) = HO2HR(K,J)   € À€i    9 CONTINUE  €  €iC     ****      €€@€iC     ****     /HOXUPP/ AND /HOXUPP2/   € €€iC     ****      €
€À€i      FHOXF(I) = FHOX   €  €iC     ****      €€@€iC     ****     /IONHOX/ AND /IONHOX2/   € €€iC     ****      € À€i      DO 10 K = 1,KMAXP1€€  €i	PHOXICF(K,I) = PHOXIC(K)       €  @€i   10 CONTINUE  €  €€iC     ****      €€ À€iC     ****     /NEWRATF/ AND /NEWRAT2F/ € ! €iC     ****      €	€!@€i      DO 11 J = 1,3     €	€!€€i	DO 11 K = 1,KMAXP1     €€!À€i	  HO2OHF(K,I,J) = HO2OH(K,J)   € " €i   11 CONTINUE  € "@€iC     ****      €€"€€iC     ****     /OPION/ AND /OPION2/     € "À€iC     ****      € # €i      DO 12 K = 1,KMAXP1€
€#@€i	FIOPF(K,I) = FIOP(K)   € #€€i   12 CONTINUE  € #À€iC     ****      €€$ €iC     ****     /PSIB/ AND /PSIB2/       € $@€iC     ****      €
 $€€i      DO 13 J = 1,10    €
€$À€i	PARBF(I,J) = PARB(J)   € % €i   13 CONTINUE  € %@€iC     ****      €€%€€iC     ****     /PHOTO/ AND /PHOTO2/     € %À€iC     ****      €	 & €i      DO 14 J =1,7      €	€&@€i	XNOF(I,J) = XNO(J)     € &€€i   14 CONTINUE  € &À€iC     ****      €€' €iC     ****     /PHOTOO/ AND /PHOTOO2/   € '@€iC     ****      €	€'€€i      DO 15 J = 1,2     €	€'À€i	XHEF(I,J) = XHE(J)     € ( €i   15 CONTINUE  € (@€iC     ****      €€(€€iC     ****     /RATEBLK/ AND /RATEBLK2/ € (À€iC     ****      €
 ) €i      DO 16 J = 1,25    €	€)@€i	DO 16 K = 1,KMAXP1     €€)€€i	  RKKF(J,K,I) = RKK(J,K)       € )À€i   16 CONTINUE  €
 * €i      DO 17 J = 1,10    €	€*@€i	DO 17 K = 1,KMAXP1     €€*€€i	  ALPPF(J,K,I) = ALPP(J,K)     € *À€i   17 CONTINUE  €
 + €i      DO 18 J = 1,20    €	€+@€i	DO 18 K = 1,KMAXP1     €€+€€i	  RBBF(J,K,I) = RBB(J,K)       € +À€i   18 CONTINUE  €
 , €i      DO 19 J = 1,50    €	€,@€i	DO 19 K = 1,KMAXP1     €€,€€i	  RKMMF(J,K,I) = RKMM(J,K)     € ,À€i   19 CONTINUE  € - €iC     ****      €€-@€iC     ****     /RJH2O2O/ AND /RJH2O2O2/ € -€€iC     ****      € -À€i      DO 21 K = 1,KMAXP1€€. €i	RJH2O2F(K,I) = RJH2O2(K)       € .@€i   21 CONTINUE  € .€€iC     ****      €€.À€iC     ****     /SOURCE/ AND /SOURCE2/   € / €iC     ****      €	€/@€i      DO 20 J = 1,6     €	€/€€i	DO 20 K = 1,KMAXP1     €€/À€i	  FS11F(K,I,J) = FS11(K,J)     € 0 €i   20 CONTINUE  € 0@€i      RETURN    €€0€€i      END               €€ @€jªDK MINOR       €€ €€j      SUBROUTINE MINOR(NX,NXNM,RMX,PHIX,ALFAX,IBND,IBNDB,PXB)   €  À€j      SAVE      €  €jC     ****      € @€jC     ****     COMMON DECKS:    €€€€jC     ****       PARAMS, ASAS, CONS, INDEX, VSCR,       € À€jC     ****      €  €jC     ****      € @€jC     ****     ADVANCES MINOR SPECIES BY ONE TIME STEP  € €€jC     ****     NX = INDEX OF MINOR SPECIES FOR TIME T(N)€ À€jC     ****     NXNM = INDEX OF MINOR SPECIES FOR TIME T(N-1)    €€ €jC     ****     NPSDH = INDEX FOR HORIZONTAL DIFFUSION FIELD     € @€jC     ****     RMX = MOLECULAR WEIGHT OF MINOR SPECIES  € €€jC     ****     PHIX = DIFFUSION VECTOR, PHI(X,N),N=1,3  € À€jC     ****     ALFAX = THERMAL DIFFUSION COEFFICIENT    €#€ €jC     ****     IBND=0 IF DIFFUSIVE FLUX GIVEN AT UPPER BOUNDARY (IN T4) €!€@€jC     ****         =1 IF TOTAL FLUX GIVEN AT UPPER BOUNDARY (IN T4)     € €€jC     ****         NOTE: FLUXES ARE IN 1/(CM**2*SEC)    € À€jC     ****      €  €jC     ****     ON ENTRY:€ @€jC     ****     LOWER BOUNDARY:  € €€jC     ****        IF IBNDB = 0, THEN    €$ À€jC     ****        T1 = A, T2 = B, T3 = C DEFINE LOWER BOUNDARY CONDITION€  €jC     ****        WHERE       A*DPSX/DZ + B*PSX + C = 0.€ @€jC     ****        IF IBNDB = 1, THEN    €" €€jC     ****        T3 = UPWARD FLUX (1/(SEC*CM**2)) AT LOWER BOUNDARY    € À€jC     ****     UPPER BOUNDARY:  €!€ €jC     ****        T4 = PHI, GIVES UPWARD FLUX OF SPECIES X AT UPPER     €€@€jC     ****             BOUNDARY (1/(SEC*CM**2)) €€€€jC     ****     SOURCES: €"€À€jC     ****        S1 = SX/N(X), WHERE SX IS PORTION OF NUMBER DENSITY   €#€ €jC     ****        SOURCE PROPORTIONAL TO N(X), THE MINOR SPECIES NUMBER €€@€jC     ****        DENSITY       €#€€€jC     ****        S2 = S0, PORTION OF NUMBER DENSITY SOURCE INDEPENDENT € À€jC     ****        OF N(X).      € 	 €jªCA PARAMS      €€	@€jªCA ASAS2       €€	€€jªCA CONS2       €€	À€jªCA INDEX       € 
 €jªCA VSCR€€
@€j      DIMENSION PHIX(3),B(2,2),PHI(2,3) €€
€€j      DIMENSION PXB(IMXP)       €!€
À€j      DATA PHI/0.,0.673,1.35,0.,1.11,0.769/,TAU/1.86E+3/,T00/273./,     €	€ €j     1SMALL/1.E-15/     €€@€j      EQUIVALENCE (C(40),B)     € €€jC     ****     S13 = HORIZONAL ADVECTION€
 À€j      DO 24 K=1,KMAX    €  €j      CALL ADVECL(NX,S13(1,K),K)€ @€j   24 CONTINUE  € €€jC     ****     PERIODIC POINTS  €€À€j      DO 42 I=1,2       €  €j      T1(I)=T1(I+IMAX)  € @€j      T1(I+IMAXP2)=T1(I+2)      € €€j      T2(I)=T2(I+IMAX)  € À€j      T2(I+IMAXP2)=T2(I+2)      €  €j      T3(I)=T3(I+IMAX)  € @€j      T3(I+IMAXP2)=T3(I+2)      € €€j      T4(I)=T4(I+IMAX)  € À€j      T4(I+IMAXP2)=T4(I+2)      €  €j      DO 42 K=1,KMAXP1  € @€j      S1(I,K)=S1(I+IMAX,K)      € €€j      S1(I+IMAXP2,K)=S1(I+2,K)  € À€j      S2(I,K)=S2(I+IMAX,K)      €  €j      S2(I+IMAXP2,K)=S2(I+2,K)  € @€j      S13(I,K)=S13(I+IMAX,K)    € €€j      S13(I+IMAXP2,K)=S13(I+2,K)€ À€j   42 CONTINUE  €€ €jC     ****     S14=SX/N(X)(K+1/2),  S15=S0(K+1/2)       €	€@€j      DO 1 K=1,KMAX     €	€€€j      DO 1 I=1,LEN1     €€À€j      S15(I,K)=.5*(S2(I,K)+S2(I,K+1))   €€ €j      S14(I,K)=.5*(S1(I,K)+S1(I,K+1))   € @€j    1 CONTINUE  € €€jC     ****     T4 = PS1(-1/2), T5 = PS2(-1/2), T6 = MBAR(-1/2)  €	 À€j      NPS1K=NJ+NPS      €  €j      NPS2K=NPS1K+KMAXP1€	€@€j      DO 2 I=1,LEN1     €€€€j      T7(I)=T4(I)       €€À€j      T4(I)=B(1,1)*F(I,NPS1K)+B(1,2)*F(I,NPS2K)+FB(I,1) €€ €j      T5(I)=B(2,1)*F(I,NPS1K)+B(2,2)*F(I,NPS2K)+FB(I,2) €$ @€j      T6(I)=1./(T4(I)/RMASS(1)+T5(I)/RMASS(2)+(1.-T4(I)-T5(I))/RMASS(3))€ €€j    2 CONTINUE  € À€jC     ****     S12 = MBAR(K+1/2)€
  €j      NPS1K=NJ+NPS-1    € @€j      NPS2K=NPS1K+KMAXP1€€€€j      DO 3 K = 1,KMAXP1 €	€À€j      DO 3 I=1,LEN1     €!€ €j      S12(I,K)=1./(F(I,NPS1K+K)/RMASS(1)+F(I,NPS2K+K)/RMASS(2)+(1.-     € @€j     1  F(I,NPS1K+K)-F(I,NPS2K+K))/RMASS(3))    € €€j    3 CONTINUE  € À€jC     ****     S11 = MBAR(K), S10 = DM/DZ(K), S9 = PS2(K),      €  €jC     ****     S8 = PS1(K), S7 = DPS2/DZ(K), S6 = DPS1/DZ(K)    €€@€jC     ****     LOWER BOUNDARY   €	€€€j      DO 4 I=1,LEN1     € À€j      S11(I,1)=.5*(T6(I)+S12(I,1))      €  €j      S10(I,1)=(S12(I,1)-T6(I))/C(3)    €€@€j      S9(I,1)=.5*(T5(I)+F(I,NPS2K+1))   €€€€j      S8(I,1)=.5*(T4(I)+F(I,NPS1K+1))   €€À€j      S7(I,1)=(F(I,NPS2K+1)-T5(I))/C(3) €€ €j      S6(I,1)=(F(I,NPS1K+1)-T4(I))/C(3) €€@€j      T4(I)=T7(I)       € €€j    4 CONTINUE  € À€jC     ****     LEVELS K=2,KMAXP1€€ €j      DO 5 K = 2,KMAXP1 €	€@€j      DO 5 I=1,LEN1     €€€€j      S11(I,K)=.5*(S12(I,K)+S12(I,K-1)) €€À€j      S10(I,K)=(S12(I,K)-S12(I,K-1))/C(3)       €  €j      S9(I,K)=.5*(F(I,NPS2K+K)+F(I,NPS2K+K-1))  € @€j      S8(I,K)=.5*(F(I,NPS1K+K)+F(I,NPS1K+K-1))  € €€j      S7(I,K)=(F(I,NPS2K+K)-F(I,NPS2K+K-1))/C(3)€ À€j      S6(I,K)=(F(I,NPS1K+K)-F(I,NPS1K+K-1))/C(3)€  €j    5 CONTINUE  € @€jC     ****     S5=T(TOT,K)      €€€€jC     ****     LOWER AND UPPER BOUNDARIES       €
 À€j      NTK=NJ+NT+KMAX    €	€ €j      DO 6 I=1,LEN1     € @€j      S5(I,1)=T0(1)+F(I,NTK)    € €€j      S5(I,KMAXP1)=T0(KMAXP1)+F(I,NTK-1)€ À€j    6 CONTINUE  €  €jC     ****     LEVELS K=2,KMAX  €€@€j      NTK=NJ+NT €	€€€j      DO 7 K=2,KMAX     €€À€j      NTK=NTK+1 €	€ €j      DO 7 I=1,LEN1     € @€j      S5(I,K)=T0(K)+.5*(F(I,NTK)+F(I,NTK-1))    € €€j    7 CONTINUE  € À€jC     ****     EVALUATE ALFA12, ALFA21, ALFAM1, ALFAM2  €   €j      ALFA12=PHI(1,2)-PHI(1,3)  €  @€j      ALFA21=PHI(2,1)-PHI(2,3)  €  €€j      ALFAX1=PHIX(1)-PHIX(3)    €  À€j      ALFAX2=PHIX(2)-PHIX(3)    €€! €jC     ****     S15=(S0*MX/NMBAR)(K+1/2) €€!@€j      NTK=NJ+NT-1       €	€!€€j      DO 8 K=1,KMAX     €€!À€j      NTK=NTK+1 €	€" €j      DO 8 I=1,LEN1     €#€"@€j      S15(I,K)=S15(I,K)*RMX*C(84)*(.5*(T0(K)+T0(K+1))+F(I,NTK))/(C(81)* €€"€€j     1EXPS(K)*S12(I,K)) € "À€j    8 CONTINUE  €"€# €jC     ****     IF IBNDB = 1, SET UP FLUX BOUNDARY CONDITION AT BOTTOM   € #@€j      IF(IBNDB.EQ.1)THEN€ #€€j          DO 26 I=1,LEN1€	 #À€j          T1(I)=1.      €€$ €j          T2(I)=S10(I,1)/S11(I,1)       € $@€j	  T3(I)=T3(I)*C(54)*RMX/(DIFK(I,1)*C(81)*EXPS(1)*C(87)*C(85))  €	 $€€j   26     CONTINUE      €€$À€j      ENDIF     € % €jC     ****     S1 = ALFA(1,1,K), S2 = ALFA(1,2,K),      €€%@€jC     ****     S2 = ALFA(2,1,K), S4 = ALFA(2,2,K)       €
€%€€j      DO 9 K=1,KMAXP1   €	€%À€j      DO 9 I=1,LEN1     € & €j      S1(I,K)=-(PHI(1,3)+ALFA12*S9(I,K))€ &@€j      S2(I,K)=ALFA12*S8(I,K)    € &€€j      S3(I,K)=ALFA21*S9(I,K)    € &À€j      S4(I,K)=-(PHI(2,3)+ALFA21*S8(I,K))€ ' €jC     ****     S12=EX(K)€#€'@€j      S12(I,K)=((ALFAX1*S4(I,K)-ALFAX2*S3(I,K))*(S6(I,K)-(1.-(RMASS(1)+ €"€'€€j     1  S10(I,K))/S11(I,K))*S8(I,K))+(ALFAX2*S1(I,K)-ALFAX1*S2(I,K))*   €# 'À€j     2  (S7(I,K)-(1.-(RMASS(2)+S10(I,K))/S11(I,K))*S9(I,K)))/(S1(I,K)*  €€( €j     3  S4(I,K)-S2(I,K)*S3(I,K))+1.-(RMX+S10(I,K))/S11(I,K)     €€(@€jC     ****     S10 = (DM/DZ)/MBAR       € (€€j      S10(I,K)=S10(I,K)/S11(I,K)€ (À€jC     ****     S11 = AX(K)      €"€) €j      S11(I,K)=-S11(I,K)/(TAU*RMASS(3))*(T00/S5(I,K))**0.25/(PHIX(3)+   € )@€j     1  ALFAX1*S8(I,K)+ALFAX2*S9(I,K))  € )€€j    9 CONTINUE  €"€)À€jC     ****     S12=EX-ALFAX*D/DS(LN(T(TOT))  (THERMAL DIFFUSION TERM)   €
 * €j      DO 21 K=2,KMAX    €
 *@€j      DO 21 I=1,LEN1    €€*€€j      S12(I,K)=S12(I,K)-ALFAX*(S5(I,K+1)-S5(I,K-1))/(2.*C(3)*   € *À€j     1  S5(I,K))€ + €j   21 CONTINUE  € +@€jC     ****     LOWER BOUNDARY AND UPPER BOUNDARY€
 +€€j      DO 22 I=1,LEN1    € +À€j      S12(I,1)=S12(I,1)-ALFAX*(S5(I,2)-S5(I,1))/(C(3)*S5(I,1))  € , €j   22 CONTINUE  € ,@€j      K=KMAXP1  €
 ,€€j      DO 23 I=1,LEN1    € ,À€j      S12(I,K)=S12(I,K)-ALFAX*(S5(I,K)-S5(I,K-1)€ - €j     1  )/(C(3)*S5(I,K))€ -@€j   23 CONTINUE  €€-€€jC     ****     T7=DFACT €
 -À€j      DO 10 I=1,LEN1    € . €j      T7(I)=1.  € .@€j   10 CONTINUE  €€.€€jC     ****     SHAPIRO SMOOTHER,   S10=NXNM  (ELIMINATED)       €
€.À€j      NXNMK=NJ+NXNM-1   € / €j      DO 12 K=1,KMAXP1  € /@€j      DO 12 I=3,LEN1-2  € /€€j      S9(I,K)=S10(I,K)  €"€/À€j       S10(I,K)=F(I,NXNMK+K)-C(26)*(F(I+2,NXNMK+K)+F(I-2,NXNMK+K)-4.*   € 0 €j     1  (F(I+2,NXNMK+K)+F(I-2,NXNMK+K))+6.*F(I,NXNMK+K))€ 0@€j   12 CONTINUE  € 0€€jC     ****     S1 = P, S2 = Q, S3 = R, S4= F    €€0À€j      NWK=NJ+NW-1       €
 1 €j      DO 13 K=1,KMAX    €€1@€j      NWK=NWK+1 €
 1€€j      DO 13 I=1,LEN1    € €1À€j      S1(I,K)=S11(I,K)/C(3)*(1./C(3)+.5*S12(I,K))-EXPS(K)*(C(87)*       € 2 €j     1  DIFK(I,K)*T7(I)*(1./C(3)-.5*S9(I,K))+0.25*(F(I,NWK)+    €€2@€j     2  F(I,NWK+1)))/C(3)       €"€2€€j      S3(I,K)=S11(I,K+1)/C(3)*(1./C(3)-.5*S12(I,K+1))-EXPS(K)*(C(86)*   €  2À€j     1  DIFK(I,K+1)*T7(I)*(1./C(3)+.5*S9(I,K+1))-0.25*(F(I,NWK)+€€3 €j     2  F(I,NWK+1)))/C(3)       €$ 3@€j      S2(I,K)=-(S11(I,K)/C(3)*(1./C(3)-.5*S12(I,K))+S11(I,K+1)/C(3)*(1./€!€3€€j     1  C(3)+.5*S12(I,K+1)))+EXPS(K)*((C(87)*DIFK(I,K)*(1./C(3)+.5*     €# 3À€j     2  S9(I,K))+C(86)*DIFK(I,K+1)*(1./C(3)-.5*S9(I,K+1)))*T7(I)/C(3)-  € 4 €j     3  S14(I,K)+C(7))  €€4@€j      S4(I,K)=EXPS(K)*(S10(I,K)*C(7)-S13(I,K)+S15(I,K)) € 4€€j   13 CONTINUE  €€4À€jC     ****     BOUNDARIES       € 5 €jC     ****     MODIFY EX(KMAXP1) IF IBND=1      €€5@€j      IF(IBND.EQ.1)THEN €€5€€j           NWK=NJ+NW+KMAX       €€5À€j           DO 25 I=1,LEN1       € 6 €j           T7(I)=EXPS(KMAX)*C(86)*F(I,NWK)/S11(I,KMAXP1)€ 6@€j           S12(I,KMAXP1)=S12(I,KMAXP1)-T7(I)    €	€6€€j   25      CONTINUE     €€6À€j      ENDIF     €
 7 €j      DO 14 I=1,LEN1    €€7@€jC     ****     LOWER BOUNDARY   €$ 7€€j      S2(I,1)=S2(I,1)+S1(I,1)*(T1(I)+.5*T2(I)*C(3))/(T1(I)-.5*T2(I)*C(3)€€7À€j     1) € 8 €j      S4(I,1)=S4(I,1)-S1(I,1)*T3(I)*C(3)/(T1(I)-.5*T2(I)*C(3))  € 8@€j      S1(I,1)=0.€ 8€€jC     ****     UPPER BOUNDARY  (UPWARD FLUX GIVEN)      €# 8À€j      S2(I,KMAX)=S2(I,KMAX)+S3(I,KMAX)*(1.+.5*S12(I,KMAXP1)*C(3))/(1.-  € 9 €j     1.5*S12(I,KMAXP1)*C(3))    €# 9@€j      S4(I,KMAX)=S4(I,KMAX)-S3(I,KMAX)*C(54)*RMX/(C(81)*S11(I,KMAXP1)*  € 9€€j     1C(85))*T4(I)*C(3)/(1.-.5*S12(I,KMAXP1)*C(3))      €	€9À€j      S3(I,KMAX)=0.     € : €j   14 CONTINUE  €€:@€jC     ****     SOLVE TRIDIAGONAL SYSTEM €	€:€€j      NXNPK=NJNP+NX     €#€:À€j      CALL TRSOLV(S1,S2,S3,S4,S7,F(1,NXNPK),LEN1,3,LEN1-2,KMAXP1,1,KMAX €€; €j     1,1)       € €;@€jC     ****     INSERT VALUE OF NXNPK(KMAXP1) USING UPPER BOUNDARY       € ;€€jC     ****                     CONDITION€ ;À€j      NXNPK=NJNP+NX+KMAX€ < €j      DO 15 I=3,LEN1-2  €!€<@€j      F(I,NXNPK)=(C(54)*RMX/(C(81)*S11(I,KMAXP1)*C(85))*T4(I)*C(3)+     €" <€€j     1(1.+.5*S12(I,KMAXP1)*C(3))*F(I,NXNPK-1))/(1.-.5*S12(I,KMAXP1)*    €€<À€j     2C(3))     € = €j   15 CONTINUE  €€=@€jC     ****     INSERT PERIODIC POINTS   € =€€j      K1=NJNP+NX€ =À€j      K2=K1+KMAX€€> €j      DO 18 N=1,2       €€>@€j      DO 19 I=1,2       €	€>€€j      DO 19 K=K1,K2     € >À€j      F(I,K)=F(I+IMAX,K)€ ? €j      F(I+IMAXP2,K)=F(I+2,K)    € ?@€j   19 CONTINUE  €	 ?€€j      K1=NJNP+NXNM      € ?À€j      K2=K1+KMAX€ @ €j   18 CONTINUE  €€@@€jC     ****     TIME SMOOTHING   €
€@€€j      NXNMK=NJ+NXNM-1   €€@À€j      NXK=NJ+NX-1       €
€A €j      NXNPK=NJNP+NX-1   €€A@€j      NXMNK=NJNP+NXNM-1 € A€€j      DO 17 K=1,KMAXP1  €
 AÀ€j      DO 17 I=1,LEN1    €"€B €j      F(I,NXMNK+K)=C(30)*F(I,NXK+K)+C(31)*(F(I,NXNMK+K)+F(I,NXNPK+K))   € B@€j   17 CONTINUE  € B€€jC     ****     PXB = PSX(-1/2)  €	€BÀ€j      NXNPK=NJNP+NX     € C €j      DO 16 I = 1,LEN1  €€C@€j      PXB(I)=(F(I,NXNPK)*(T1(I)+.5*T2(I)*C(3))+T3(I)*C(3))/     €€C€€j     1(T1(I)-.5*T2(I)*C(3))     € CÀ€j   16 CONTINUE  €€D €jC     ****     INSURE NON-NEGATIVE NX   €	€D@€j      NXK=NJNP+NX-1     €€D€€j      NXNMK=NJNP+NXNM-1 € DÀ€j      DO 20 K=1,KMAXP1  €
 E €j      DO 20 I=1,LEN1    €€E@€j      F(I,NXK+K)=CVMGP(F(I,NXK+K),SMALL,F(I,NXK+K)-SMALL)       €€E€€j      F(I,NXNMK+K)=CVMGP(F(I,NXNMK+K),SMALL,F(I,NXNMK+K)-SMALL) € EÀ€j   20 CONTINUE  € F €j      RETURN    €€F@€j      END               €  @€kªDK PLOT€€ €€k      SUBROUTINE PLOT(IFP,N,NJ) €  À€k      SAVE      €  €kC     ****      €€@€kC     ****     PROCEDURE TO PLOT OUT MASS MIXING RATIOS € €€kC     ****      € €À€kC     ****     ARRAY IFP CONTAINS INDICES OF FIELDS TO BE PLOTTED       €  €kC     ****      €€@€kC     ****     FOR EXAMPLE:     € €€kC     ****              NPS    FOR    O2€€À€kC     ****              NPS2   FOR    O €  €kC     ****              NPSA   FOR    AR€ @€kC     ****              NPSH   FOR    HE€ €€kC     ****              NPNO   FOR    NO€€À€kC     ****              NPN4S  FOR    N4S       €€ €kC     ****              NPCH4  FOR    CH4       € @€kC     ****              NPH2   FOR    H2€€€€kC     ****              NPCO2  FOR    CO2       € À€kC     ****              NPCO   FOR    CO€€ €kC     ****              NPHOX  FOR    HOX       €€@€kC     ****              NPH2O  FOR    H2O       € €€kC     ****      € À€kC     ****     N IS NUMBER OF INDICES PASSED (MAX IS 28)€  €kC     ****      €€@€kC     ****     NJ IS BUFFER POINTER     € €€kC     ****      € À€kªCA PARAMS      €€ €kªCA ASAS2       €€@€kªCA CONS2       €€€€k      DIMENSION IFP(1),LABEL(28),LAB2(28)       €#€À€k      DATA LABEL/3H O2,3H  O,3H AR,3H HE,3H NO,3HN4S,3HCH4,3H H2,3HCO2, €  €k     13H CO,3HHOX,3HH2O/€ @€kC     ****      € €€kC     ****     SHIFT REQUIRED LABELS TO ARRAY LAB2      € À€kC     ****      € 	 €k      DO 1 I=1,N€€	@€k       LAB2(I)=LABEL(MOD(IFP(I)/KMAXP1-1,13))   € 	€€k    1 CONTINUE  € 	À€kC     ****      € 
 €kC     ****     PLOT FIELDS      € 
@€kC     ****      €
 
€€k      CALL GSCLIP(0)    €	 
À€k      DO 2 I = 1,N      €  €k	CALL CONREC(F(3,NJ+IFP(I)),IMAXP4,IMAX+1,KMAXP1,0.,0.,0.,0,0,  €€@€k     1    -1430B)       € €€kC       ****    € À€kC       ****     WRITE TITLE    €  €kC       ****    €€@€k	CALL GETSET(FXA,FXB,FYA,FYB,XC,XD,YC,YD,LTYPE) € €€k	MXA = KFPX(FXA)€ À€k	MXB = KFPX(FXB)€  €k	MYA = KFPY(FYA)€ @€k	MYB = KFPY(FYB)€€€€k	MX = (MXA+MXB)/2       € À€k	MY = MYB+48    €  €k	CALL PWRIT(CPUX(MX),CPUY(MY),LAB2(I),3,2,0,0)  €€@€k	CALL FRAME     € €€k    2 CONTINUE  €
 À€k      CALL GSCLIP(1)    €  €k      RETURN    €€@€k      END               €€ @€lªDK RATES       €  €€l      SUBROUTINE RATES(FT)      €  À€l      SAVE      €  €lC     ****      € @€lC     ****     CALCULATES K12, K13, K14, K15 AT ALL LEVELS      € €€lC     ****      € À€lC     ****     COMMON DECKS:    €€ €lC     ****       PARAMS, CONS, CRATES   € @€lC     ****      € €€lªCA PARAMS      €€À€lªCA CONS2       €€ €lªCA CRATES2     €€@€lªCA OXOXOX2     € €€l      DIMENSION FT(KMXP,1)      € À€lC     ****      €€ €lC     ****     PERIODIC POINTS FOR OOXR € @€lC     ****      €	 €€l      DO 2 I = 1,2      €	 À€l	DO 2 K = 1,KMAXP1      €  €l	  OOXR(K,I) = OOXR(K,I+IMX)    € @€l	  OOXR(K,I+IMX+2) = OOXR(K,I+2)€ €€l    2 CONTINUE  €
€À€l      DO 1 K=1,KMAXP1   €  €l	DO 1 I=1,LEN1  € @€lC       RK13(K,I)=2.15E-34*EXP(345./FT(K,I))    € €€lC       RK14(K,I)=2.15E-34*EXP(345./FT(K,I))    € À€lC       RK15(K,I)=8.82E-35*EXP(575./FT(K,I))    €
  €lC NEW   JPL-85 RATES    €€@€l	RK12(K,I)=OOXR(K,I)*OOXR(K,I)*9.59E-34*EXP(480./FT(K,I))       €
 €€lC       RK12(K,I)=0.    €€À€lC       RK13(K,I)=OOXR(K,I)*6.0E-34*(300./FT(K,I))**2.8 €€ €lC       RK14(K,I)=OOXR(K,I)*6.0E-34*(300./FT(K,I))**2.8 €€@€lC       RK15(K,I)=OOXR(K,I)*6.0E-34*(300./FT(K,I))**2.8 €€€€l	RK13(K,I)=0.   €€À€l	RK14(K,I)=0.   €€	 €l	RK15(K,I)=0.   € 	@€l    1 CONTINUE  € 	€€l      RETURN    €€	À€l      END               €€ @€mªDK UTILITY     €  €€mCDIR$ NOLIST    €€ À€m      SUBROUTINE TRSOLV(A,B,C,F,W,X,IF,I1,I2,KF,K1,K2,NF)       €  €m      SAVE      € @€mC     ****     SOLVES TRIDIAGONAL SYSTEM€!€€€mC     ****     A(I,K)*X(I,K-1)+B(I,K)*X(I,K)+C(I,K)*X(I,K+1)=F(I,K)     €!€À€mC     ****     WHERE I=I1,I2,1 AND K=K1,K2,1 AND A(I,K1)=C(I,K2)=0.     €€ €mC     ****     ARRAYS A,B,C,F,X ARE DIMENSIONED (IF,KF) €  @€mC     ****     WHERE 1.LE.I1.LT.I2.LE.IF AND 1.LE.K1.LT.K2.LE.KF€ €€mC     ****     W IS WORK SPACE DIMENSIONED (IF,3*KF)    €"€À€mC     ****     THE ROUTINE USES VECTOR OPERATIONS ON ARRAYS OF LENGTH   €€ €mC     ****     IMAX=I2-I1+1     €#€@€m      DIMENSION A(IF,KF),B(IF,KF),C(IF,KF),F(IF,KF),W(IF,3*KF),X(IF,KF) €	 €€m      IMAX=I2-I1+1      € À€m      K1P=K1+1  €	€ €m      GO TO(1,2),NF     € @€m    1 CONTINUE  € €€mC     ****          LOWER BOUNDARY      €€À€mC     ****          W(K1)=B(K1) €	  €m      DO 3 I=I1,I2      €
€@€m      W(I,K1)=B(I,K1)   € €€m    3 CONTINUE  €€À€mC     ****          SWEEP THROUGH K     €	€ €m      DO 4 K=K1P,K2     €€@€mC     ****          W(KF+K-1)=C(K-1)/W(K-1)     € €€m      CALL QVDIV0(W(I1,KF+K-1),C(I1,K-1),W(I1,K-1),IMAX)€€À€mC     ****          W(K)=A(K)*W(KF+K-1) €  €m      CALL QVMPY0(W(I1,K),A(I1,K),W(I1,KF+K-1),IMAX)    € @€mC     ****          W(K)=B(K)-W(K)      €€€€m      CALL QVSUB0(W(I1,K),B(I1,K),W(I1,K),IMAX) € À€m    4 CONTINUE  €  €m    2 CONTINUE  € @€mC     ****          LOWER BOUNDARY      € €€mC     ****          W(2*KF+K1)=F(K1)/W(K1)      €€À€m      CALL QVDIV0(W(I1,2*KF+K1),F(I1,K1),W(I1,K1),IMAX) €	€	 €m      DO 5 K=K1P,K2     € 	@€mC     ****          W(2*KF+K)=A(K)*W(2*KF+K-1)  €€	€€m      CALL QVMPY0(W(I1,2*KF+K),A(I1,K),W(I1,2*KF+K-1),IMAX)     € 	À€mC     ****          W(2*KF+K)=F(K)-W(2*KF+K)    €€
 €m      CALL QVSUB0(W(I1,2*KF+K),F(I1,K),W(I1,2*KF+K),IMAX)       € 
@€mC     ****          W(2*KF+K)=W(2*KF+K)/W(K)    €€
€€m      CALL QVDIV0(W(I1,2*KF+K),W(I1,2*KF+K),W(I1,K),IMAX)       € 
À€m    5 CONTINUE  €€ €mC     ****          BACK SUBSTITUTION   € @€mC     ****          X(K2)=W(2*KF+K2)    €	 €€m      DO 6 I=I1,I2      € À€m      X(I,K2)=W(I,2*KF+K2)      €  €m    6 CONTINUE  €
 @€m      DO 7 KK=K1P,K2    € €€m      K=K1+K2-KK€€À€mC     ****          X(K)=W(KF+K)*X(K+1) €  €m      CALL QVMPY0(X(I1,K),W(I1,KF+K),X(I1,K+1),IMAX)    €€@€mC     ****          X(K)=W(2*KF+K)-X(K) € €€m      CALL QVSUB0(X(I1,K),W(I1,2*KF+K),X(I1,K),IMAX)    € À€m    7 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QAA0(W,X,Y,Z,N)€ À€m      SAVE      €€ €m      DIMENSION W(1),X(1),Y(1),Z(1)     €	 @€m      DO 100 I=1,N      €€€€m        W(I) = X(I) + Y(I) + Z(I)       € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QAM0(W,X,Y,Z,N)€ À€m      SAVE      €€ €m      DIMENSION W(1),X(1),Y(1),Z(1)     €	 @€m      DO 100 I=1,N      €€€€m        W(I) = ( X(I) + Y(I) ) * Z(I)   € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QAM1(W,X,Y,C,N)€ À€m      SAVE      €  €m      DIMENSION W(1),X(1),Y(1)  €	 @€m      DO 100 I=1,N      € €€m        W(I) = ( X(I) + Y(I) ) * C      € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QAM2(W,X,C,Y,N)€ À€m      SAVE      €  €m      DIMENSION W(1),X(1),Y(1)  €	 @€m      DO 100 I=1,N      € €€m        W(I) = ( X(I) + C ) * Y(I)      € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QMA0(W,X,Y,Z,N)€ À€m      SAVE      €€ €m      DIMENSION W(1),X(1),Y(1),Z(1)     €	 @€m      DO 100 I=1,N      €€€€m        W(I) = X(I) * Y(I) + Z(I)       € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QMA1(W,X,Y,C,N)€ À€m      SAVE      €  €m      DIMENSION W(1),X(1),Y(1)  €	 @€m      DO 100 I=1,N      € €€m        W(I) = X(I) * Y(I) + C  € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QMA2(W,X,C,Y,N)€ À€m      SAVE      €  €m      DIMENSION W(1),X(1),Y(1)  €	 @€m      DO 100 I=1,N      € €€m        W(I) = X(I) * C + Y(I)  € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QMA5(W,C,X,S,N)€ À€m      SAVE      €€ €m      DIMENSION W(1),X(1)       €	 @€m      DO 100 I=1,N      €€€€m        W(I) = C * X(I) + S     € À€m  100 CONTINUE  €  €m      RETURN    €€@€m      END       € €€m      SUBROUTINE QMM0(W,X,Y,Z,N)€ À€m      SAVE      €€ €m      DIMENSION W(1),X(1),Y(1),Z(1)     €	 @€m      DO 100 I=1,N      €€€€m        W(I) = X(I) * Y(I) * Z(I)       € À€m  100 CONTINUE  €   €m      RETURN    €€ @€m      END       €  €€m      SUBROUTINE QMM1(W,X,Y,C,N)€  À€m      SAVE      € ! €m      DIMENSION W(1),X(1),Y(1)  €	 !@€m      DO 100 I=1,N      € !€€m        W(I) = X(I) * Y(I) * C  € !À€m  100 CONTINUE  € " €m      RETURN    €€"@€m      END       € "€€m      SUBROUTINE QMS0(W,X,Y,Z,N)€ "À€m      SAVE      €€# €m      DIMENSION W(1),X(1),Y(1),Z(1)     €	 #@€m      DO 100 I=1,N      €€#€€m        W(I) = X(I) * Y(I) - Z(I)       € #À€m  100 CONTINUE  € $ €m      RETURN    €€$@€m      END       € $€€m      SUBROUTINE QSA0(W,X,Y,Z,N)€ $À€m      SAVE      €€% €m      DIMENSION W(1),X(1),Y(1),Z(1)     €	 %@€m      DO 100 I=1,N      €€%€€m        W(I) = X(I) - Y(I) + Z(I)       € %À€m  100 CONTINUE  € & €m      RETURN    €€&@€m      END       € &€€m      SUBROUTINE QSM0(W,X,Y,Z,N)€ &À€m      SAVE      €€' €m      DIMENSION W(1),X(1),Y(1),Z(1)     €	 '@€m      DO 100 I=1,N      €€'€€m        W(I) = ( X(I) - Y(I) ) * Z(I)   € 'À€m  100 CONTINUE  € ( €m      RETURN    €€(@€m      END       € (€€m      SUBROUTINE QSM1(W,X,Y,C,N)€ (À€m      SAVE      € ) €m      DIMENSION W(1),X(1),Y(1)  €	 )@€m      DO 100 I=1,N      € )€€m        W(I) = ( X(I) - Y(I) ) * C      € )À€m  100 CONTINUE  € * €m      RETURN    €€*@€m      END       € *€€m      SUBROUTINE QVADD0(W,X,Y,N)€ *À€m      SAVE      € + €m      DIMENSION W(1),X(1),Y(1)  €	 +@€m      DO 100 I=1,N      € +€€m        W(I) = X(I) + Y(I)      € +À€m  100 CONTINUE  € , €m      RETURN    €€,@€m      END       € ,€€m      SUBROUTINE QVADD1(W,X,C,N)€ ,À€m      SAVE      €€- €m      DIMENSION W(1),X(1)       €	 -@€m      DO 100 I=1,N      €€-€€m        W(I) = X(I) + C € -À€m  100 CONTINUE  € . €m      RETURN    €€.@€m      END       € .€€m      SUBROUTINE QVDIV0(W,X,Y,N)€ .À€m      SAVE      € / €m      DIMENSION W(1),X(1),Y(1)  €	 /@€m      DO 100 I=1,N      € /€€m        W(I) = X(I) / Y(I)      € /À€m  100 CONTINUE  € 0 €m      RETURN    €€0@€m      END       € 0€€m      SUBROUTINE QVDIV2(W,C,X,N)€ 0À€m      SAVE      €€1 €m      DIMENSION W(1),X(1)       €	 1@€m      DO 100 I=1,N      €€1€€m        W(I) = C / X(I) € 1À€m  100 CONTINUE  € 2 €m      RETURN    €€2@€m      END       € 2€€m      SUBROUTINE QVMPY0(W,X,Y,N)€ 2À€m      SAVE      € 3 €m      DIMENSION W(1),X(1),Y(1)  €	 3@€m      DO 100 I=1,N      € 3€€m        W(I) = X(I) * Y(I)      € 3À€m  100 CONTINUE  € 4 €m      RETURN    €€4@€m      END       € 4€€m      SUBROUTINE QVMPY2(W,C,X,N)€ 4À€m      SAVE      €€5 €m      DIMENSION W(1),X(1)       €	 5@€m      DO 100 I=1,N      €€5€€m        W(I) = C * X(I) € 5À€m  100 CONTINUE  € 6 €m      RETURN    €€6@€m      END       € 6€€m      SUBROUTINE QVSUB0(W,X,Y,N)€ 6À€m      SAVE      € 7 €m      DIMENSION W(1),X(1),Y(1)  €	 7@€m      DO 100 I=1,N      € 7€€m        W(I) = X(I) - Y(I)      € 7À€m  100 CONTINUE  € 8 €m      RETURN    €€8@€m      END       € 8€€m      SUBROUTINE QVSUB1(W,X,C,N)€ 8À€m      SAVE      €€9 €m      DIMENSION W(1),X(1)       €	 9@€m      DO 100 I=1,N      €€9€€m        W(I) = X(I) - C € 9À€m  100 CONTINUE  € : €m      RETURN    €€:@€m      END       €  :€€m        (     @        1­üún     @   	   1­üúZ     @      1­üú      @      1­üú     @   &   1­üú     @   ;   1­üúý     @   E   1­üúô     @   `   1­üú3    	 @   k   1­üúe    
 @	   ~   1­üúf     @
   •   1­üúÂ     @   ¬   1­üú     @   ¹   1­üúô     @   Â   1­üú:     @   Ô   1­üú"     @   Þ   1­üúh     @   ù   1­üúT     @     1­üúK     @  #   1­üú·     @  .   1­üúý     @  X   1­üúé     @  c   1­üú/     @  ¬   1­üú¸     @  ¹   1­üúþ     @  Ö   1­üúê     @  à   1­üú0     @     1­üú?     @     1­üúq     @     1­üú·     @  )   1­üú¨     @  6   1­üúÚ      @  @   1­üú     ! @   J   1­üúÌ    " @!  T   1­üúþ    # @"  \   1­üúD    $ @#  e   1­üú1    % @$  n   1­üúw    & @%  †   1­üúc    ' @&     1­üúT    ( @'  °   1­üúä    ) @(  »   1­üú    * @)  Ã   1­üú\    + @*  Ì   1­üúx    , @+  Õ   1­üúQ    - @,  ÷   1­üúÕ    . @-  7   1­üú    / @.  @   1­üúM    0 @/  J   1­üúÑ    1 @0  T   1­üú    2 @1  `   1­üú    3 @2  j   1­üúI    4 @3  x   1­üúî    5 @4  ‚   1­üú4    6 @5     1­üú     7 @6  ™   1­üúf    8 @7  ª   1­üú…    9 @8  µ   1­üú·    : @9  ¾   1­üúý    ; @:  È   1­üúõ    < @;  Ò   1­üú;    = @<  ã   1­üú'    > @=  î   1­üúm    ? @>     1­üúÄ    @ @?     1­üú.    A @@  f   1­üút    B @A  t   1­üú`    C @B  }   1­üú¦    D @C  ˜   1­üúŠ    E @D  ¢   1­üúÐ    F @E  ®   1­üú¼    G @F  ·   1­üú    H @G  Ê   1­üúÙ    I @H  Ó   1­üú    J @I  Ü   1­üú    K @J  å   1­üúQ    L @K  ð   1­üú    M @L  ú   1­üú7    N @M     1­üú8    O @N     1­üúæ    P @O  3   1­üú    Q @P  =   1­üú    R @Q  G   1­üúÑ    S @R  Q   1­üú    T @S  c   1­üú    U @T  l   1­üúI    V @U  ‚   1­üú>    W @V  Œ   1­üúp    X @W  ¿   1­üúq    Y @X  í   1­üúó€  Z @Y     1­üúµ€  [ @Z  Å   1­üúÚ€  \ @[  ø   1­üú/€  ] @\  
v   1­üú€  ^ @]  «   1­üú¯€  _ @^  9   1­üú¤€  ` @_     1­üúî€  a @`  Q   1­üúÁ€  b @a  x   1­üúó€  c @b   ã   1­üú©€  d @c  "   1­üúø€  e @d  #U   1­üú¼€  f @e  &y   1­üú	€  g @f  '¥   1­üú€  h @g  (Æ   1­üúÌ€  i @h  +n   1­üúU€  j @i  ,Ë   1­üú…€  k @j  /¯   1­üú?€  l @k  5Ä   1­üú€  m @l  6è   1­üú4€  n @m  7š   1­üú  <é   nASAS    ASASF   ASAS2   ASAS2F  CH4H2COFCH4H2CFFCH4H2CO2CH4H2C2FCONS    CONS2   CONS3   CRATES  CRATESF CRATES2 CRATES2FDISSHOX DISSHOXFDISSHOX2DISSHX2FFIELDS  FIELDSF FIELDS2 FIELDS2FFIELDT  FIELDTF FIELDT2 FIELDT2FFLUX    FLUX2   FLUX2F  H2FORG  H2FORG2 H2FORG2FHEFLUX  HEFLUX2 HEFLUX2FHOXSTUF HOXSTUFFHOXSTUF2HOXSTF2FHOXUPP  HOXUPP2 HOXUPP2FINDEX   INDEXDATIONHOX  IONHOX2 IONHOX2FNEWRAT  NEWRATF NEWRAT2 NEWRAT2FNOZNOZ  NOZNOZF NOZNOZ2 NOZNOZ2FOPION   OPION2  OPION2F OXOXOX  OXOXOXF OXOXOX2 OXOXOX2FPARAMS  PSIB    PSIBF   PSIB2   PSIB2F  PHOTO   PHOTOF  PHOTO2  PHOTO2F PHOTOO  PHOTOOF PHOTOO2 PHOTOO2FRATEBLK RATEBLK2RATBLK2FRJH2O2O RJH2O2O2RJH2O22FSOURCE  SOURCEF SOURCE2 SOURCE2FVSCR    VSCR2   VSCR3   BLKDATA CMPN4S  PPHOTO  COMP    FACE    ADVECL  CMPN2D  COMPCH4 COMPCO  COMPCO2 COMPH2  COMPH2O COMPHE  COMPARG COMPHOX COMPNO  FFACE   MINOR   PLOT    RATES   UTILITY d*€n  ;Ÿ05/30/96       m       P        HIST OK 