*/
*/
*/
*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)
*/
*/
*/
*CDK NEWRAT
      COMMON/NEWRAT/HO2OH(KMXP),HOH(KMXP),OHHOX(KMXP)                   
*CDK NEWRATF
      COMMON/NEWRATF/HO2OH(KMXP,3)
*CDK NEWRAT2
      COMMON/NEWRAT2/HO2OH(KMXP,IMXP),HOH(KMXP,IMXP),OHHOX(KMXP,IMXP)
*CDK NEWRAT2F
      COMMON/NEWRAT2/HO2OHF(KMXP,IMXP,3)
*/
*/
*/
*CDK NOZNOZ
      COMMON/NOZNOZ/SNO2NO(KMXP),RNONOZ(KMXP),FNNOZ(KMXP)
*CDK NOZNOZF
      COMMON/NOZNOZF/SNO2NO(KMXP,3)
*CDK NOZNOZ2
      COMMON/NOZNOZ2/SNO2NO(KMXP,IMXP),RNONOZ(KMXP,IMXP),
     1  FNNOZ(KMXP,IMXP)
*CDK NOZNOZ2F
      COMMON/NOZNOZ2/SNO2NOF(KMXP,IMXP,3)
*/
*/
*/
*CDK OPION
      COMMON/OPION/FIOP(KMXP)                                           
*CDK OPION2
      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 
*CDK PSIBF
      COMMON/PSIB/PARB(10)
*CDK PSIB2
      COMMON/PSIB2/PARB(IMXP),PHB(IMXP),PNOB(IMXP),PN4SB(IMXP),
     1  PCH4B(IMXP),PH2B(IMXP),PCO2B(IMXP),PCOB(IMXP),PHOXB(IMXP),
     2  PH2OB(IMXP)
*CDK PSIB2F
      COMMON/PSIB2/PARBF(IMXP,10)
*/
*/
*/
*CDK PHOTO
      COMMON/PHOTO/XNO,XN4S,XCH4,XH2,XCO2,XCO,XH2O                      
*CDK PHOTOF
      COMMON/PHOTO/XNO(7)
*CDK PHOTO2
      COMMON/PHOTO2/XNO(IMXP),XN4S(IMXP),XCH4(IMXP),XH2(IMXP),XCO2(IMXP)
     1  ,XCO(IMXP),XH2O(IMXP)
*CDK PHOTO2F
      COMMON/PHOTO2/XNOF(IMXP,7)
*/
*/
*/
*CDK PHOTOO
      COMMON/PHOTOO/XHE,XARG                                            
*CDK PHOTOOF
      COMMON/PHOTOO/XHE(2)
*CDK PHOTOO2
      COMMON/PHOTOO2/XHE(IMXP),XARG(IMXP)
*CDK PHOTOO2F
      COMMON/PHOTOO2/XHEF(IMXP,2)
*/
*/
*/
*CDK RATEBLK
      COMMON/RATEBLK/RKK(25,KMXP),ALPP(10,KMXP),RBB(20,KMXP),           
     1RKMM(50,KMXP)                                                     
*CDK RATEBLK2
      COMMON/RATEBLK2/RKK(25,KMXP,IMXP),ALPP(10,KMXP,IMXP),
     1RBB(20,KMXP,IMXP),RKMM(50,KMXP,IMXP)
*CDK RATBLK2F
      COMMON/RATEBLK2/RKKF(25,KMXP,IMXP),ALPPF(10,KMXP,IMXP),
     1RBBF(20,KMXP,IMXP),RKMMF(50,KMXP,IMXP)
*/
*/
*/
*CDK RJH2O2O
      COMMON/RJH2O2O/RJH2O2(KMXP)
*CDK RJH2O2O2
      COMMON/RJH2O2O2/RJH2O2(KMXP,IMXP)
*CDK RJH2O22F
      COMMON/RJH2O2O2/RJH2O2F(KMXP,IMXP)
*/
*/
*/
*CDK SOURCE
      COMMON/SOURCE/FS11(KMXP),FS12(KMXP),FS21(KMXP),FS22(KMXP),        
     1FS1(KMXP),FS2(KMXP)                                               
*CDK SOURCEF
      COMMON/SOURCE/FS11(KMXP,6)
*CDK SOURCE2
      COMMON/SOURCE2/FS11(KMXP,IMXP),FS12(KMXP,IMXP),FS21(KMXP,IMXP),
     1FS22(KMXP,IMXP),FS1(KMXP,IMXP),FS2(KMXP,IMXP)
*CDK SOURCE2F
      COMMON/SOURCE2/FS11F(KMXP,IMXP,6)
*/
*/
*/
*CDK VSCR
      COMMON/VSCR/S15(IMXP,KMXP),S14(IMXP,KMXP),S13(IMXP,KMXP),
     1S12(IMXP,KMXP),S11(IMXP,KMXP),S10(IMXP,KMXP),S9(IMXP,KMXP),
     2S8(IMXP,KMXP),S7(IMXP,KMXP),S6(IMXP,KMXP),S5(IMXP,KMXP),
     3S4(IMXP,KMXP),S3(IMXP,KMXP),S2(IMXP,KMXP),S1(IMXP,KMXP),T1(IMXP)
     4,T2(IMXP),T3(IMXP),T4(IMXP),T5(IMXP),T6(IMXP),T7(IMXP)
*/
*/
*/
*CDK VSCR2
      COMMON/VSCR/GAMA(IMXP,KMXP,2,2),Z(IMXP,KMXP,2),EMBAR(IMXP,
     1KMXP),AK(IMXP,2,2,2),EP(IMXP,2,2),FK(IMXP,2),PK(IMXP,
     22,2),QK(IMXP,2,2), RK(IMXP,2,2),WKM1(IMXP,2,2),WKM2(IMXP,
     32,2),WKV1(IMXP,2),WKV2(IMXP,2),WKS1(IMXP),WKS2(IMXP),WKS3(IMXP),
     4EMBAR0(IMXP),PS0(IMXP,2)
*/
*/
*/
*CDK VSCR3
      COMMON/VSCR/WK1(KMXP,IMXP),WK2(KMXP,IMXP),WK3(KMXP,IMXP),
     1  WK4(KMXP,IMXP),WK5(KMXP,IMXP),WK6(KMXP,IMXP),WK7(KMXP,IMXP),
     2  WK8(KMXP,IMXP),WK9(KMXP,IMXP),WK10(KMXP,IMXP),WK11(KMXP,IMXP),
     3  WK12(KMXP,IMXP),WK13(KMXP,IMXP)
*/
*/
*/
*/
*DK BLKDATA
*/
*/
*/
CDIR$ NOLIST
       BLOCK DATA BLOCKS
*CA PARAMS
      COMMON/BLK1/ZERO1
*CA ASAS2
      COMMON/BLK2/ZERO2
*CA CH4H2CO2
      COMMON/BLK3/ZERO3
*CA CONS2
      COMMON/BLK4/ZERO4
*CA CRATES2
      COMMON/BLK5/ZERO5
*CA DISSHOX2
      COMMON/BLK6/ZERO6
*CA FIELDS2
      COMMON/BLK7/ZERO7
*CA FIELDT2
      COMMON/BLK8/ZERO8
*CA FLUX2
      COMMON/BLK9/ZERO9
*CA H2FORG2
      COMMON/BLK10/ZERO10
*CA HEFLUX2
      COMMON/BLK11/ZERO11
*CA HOXSTUF2
      COMMON/BLK12/ZERO12
*CA HOXUPP2
      COMMON/BLK13/ZERO13
*CA INDEX
      COMMON/BLK14/ZERO14
*CA IONHOX2
      COMMON/BLK15/ZERO15
*CA NEWRAT2
      COMMON/BLK16/ZERO16
*CA NOZNOZ2
      COMMON/BLK17/ZERO17
*CA OPION2
      COMMON/BLK18/ZERO18
*CA OXOXOX2
      COMMON/BLK19/ZERO19
*CA PSIB2
      COMMON/BLK20/ZERO20
*CA PHOTO2
      COMMON/BLK21/ZERO21
*CA PHOTOO2
      COMMON/BLK22/ZERO22
*CA RATEBLK2
      COMMON/BLK23/ZERO23
*CA RJH2O2O2
      COMMON/BLK24/ZERO24
*CA SOURCE2
      COMMON/BLK25/ZERO25
*CA VSCR
      COMMON/BLK26/ZERO26
      DATA ZERO1,ZERO2,ZERO3,ZERO4,ZERO5,ZERO6,ZERO7,ZERO8,ZERO9
     1  /1.,2.,3.,4.,5.,6.,7.,8.,9./
      DATA ZERO10,ZERO11,ZERO12,ZERO13,ZERO14,ZERO15,ZERO16,ZERO17,
     1  ZERO18,ZERO19/10.,11.,12.,13.,14.,15.,16.,17.,18.,19./
      DATA ZERO20,ZERO21,ZERO22,ZERO23,ZERO24,ZERO25,ZERO26
     2  /20.,21.,22.,23.,24.,25.,26./
      END
C
C
      SUBROUTINE CHKBLKS
      SAVE

      COMMON/BLK1/ZERO1
      COMMON/BLK2/ZERO2
      COMMON/BLK3/ZERO3
      COMMON/BLK4/ZERO4
      COMMON/BLK5/ZERO5
      COMMON/BLK6/ZERO6
      COMMON/BLK7/ZERO7
      COMMON/BLK8/ZERO8
      COMMON/BLK9/ZERO9
      COMMON/BLK10/ZERO10
      COMMON/BLK11/ZERO11
      COMMON/BLK12/ZERO12
      COMMON/BLK13/ZERO13
      COMMON/BLK14/ZERO14
      COMMON/BLK15/ZERO15
      COMMON/BLK16/ZERO16
      COMMON/BLK17/ZERO17
      COMMON/BLK18/ZERO18
      COMMON/BLK19/ZERO19
      COMMON/BLK20/ZERO20
      COMMON/BLK21/ZERO21
      COMMON/BLK22/ZERO22
      COMMON/BLK23/ZERO23
      COMMON/BLK24/ZERO24
      COMMON/BLK25/ZERO25
      COMMON/BLK26/ZERO26
C     WRITE(6,100)ZERO1,ZERO2,ZERO3,ZERO4,ZERO5,ZERO6,ZERO7,ZERO8,ZERO9,
C    1  ZERO10,ZERO11,ZERO12,ZERO13,ZERO14,ZERO15,ZERO16,ZERO17,
C    2  ZERO18,ZERO19,ZERO20,ZERO21,ZERO22,ZERO23,ZERO24,ZERO25,ZERO26
C 100 FORMAT(10E12.4)
      RETURN
      END
*/
*DK CMPN4S
*/
      SUBROUTINE CMPN4S                                                 
CDIR$ BOUNDS
      SAVE
C     ****                                                              
C     ****     ADVANCE N4S COMPOSITION BY ONE TIME STEP                 
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CONS, FIELDS, FIELDT, FLUX, INDEX, PHOTO,
C     ****          PSIB, RATEBLK, VSCR
C     ****                                                              
*CA PARAMS
*CA ASAS2
*CA CONS2
*CA FIELDS2
*CA FIELDT2
*CA FLUX2
*CA INDEX
*CA PHOTO2
*CA PSIB2
*CA RATEBLK2
*CA VSCR
      DIMENSION PHIN4S(3)                                               
      DATA RMN4S,PHIN4S/14.,0.651,0.731,0.741/                          
C     ****
C     ****     UPPER BOUNDARY: DIFFUSIVE EQUILIBRIUM                    
C     ****     LOWER BOUNDARY: PHOTO-CHEMICAL EQUILIBRIUM               
C     ****                                                              
      DO 1 I=3,LEN1-2
	T4(I)=0.
	T1(I)=0.
	T2(I)=1.
	T3(I)=-XN4S(I)*RMN4S/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+
     1    FNN2(1,I)*RMASS(3))
    1   CONTINUE
C     ****                                                              
C     ****     SOURCES                                                  
C     ****                                                              
      DO 2 K=1,KMAXP1                                                   
	DO 2 I=3,LEN1-2
	  S1(I,K)=-(FK4O2P(K,I)+RBB(1,K,I)*FNO2(K,I)+RBB(3,K,I)*
     1      FNNO(K,I)+RBB(10,K,I)*FNOH(K,I))
	  S2(I,K)=PPN4S(K,I)+RBB(4,K,I)*FNN2D(K,I)*FNO(K,I)+RBB(5,K,I)*
     1      FNN2D(K,I)*FNE(K,I)+RBB(7,K,I)*FNN2D(K,I)+RBB(8,K,I)*
     2      FNNO(K,I)
    2 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
      IBND=0
      IBNDB=0
      ALFA=0.                                                           
      CALL MINOR(NPN4S,NN4SNM,RMN4S,PHIN4S,ALFA,IBND,IBNDB,PN4SB)       
      RETURN                                                            
      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)                          
      DO 112 L=1,2
	NPSK=NPSNM+(L-1)*KMAXP1+K-1
	NPSNMK=NJ+NPSK
	DO 112 I=3,IMAXP2
	  WKV2(I,L)=F(I,NPSNMK)-C(26)*(F(I+2,NPSNMK)+F(I-2,NPSNMK)-4.*
     1       (F(I+1,NPSNMK)+F(I-1,NPSNMK))+6.*F(I,NPSNMK))
  112 CONTINUE
C     ****     COMPLETE CALCULATION OF RHS IN FK                        
      DO 12 L=1,2                                                       
	DO 12 I=3,IMAXP2
	  FK(I,L)=EXPS(K)*(WKV2(I,L)*C(7)-FK(I,L))
   12 CONTINUE                                                          
C     ****     INSERT PERIODIC POINTS
      DO 121 L = 1,2
	DO 121 I = 1,2
	  FK(I,L) = FK(I+IMAX,L)
	  FK(I+IMAXP2,L) = FK(I+2,L)
  121 CONTINUE
      IF(K.EQ.1) GO TO 13                                               
      IF(K.EQ.KMAX) GO TO 14                                            
      GO TO 15                                                          
   13 CONTINUE                                                          
C     ****     LOWER BNDRY                                              
      DO 16 L=1,2                                                       
	DO 16 M=1,2
	  DO 16 KK=1,2
	    DO 16 I=1,LEN1
	      QK(I,L,M)=QK(I,L,M)+PK(I,L,KK)*B(KK,M)
   16 CONTINUE                                                          
      DO 33 L=1,2                                                       
	DO 33 M=1,2
	  DO 33 I=1,LEN1
	    FK(I,L)=FK(I,L)-PK(I,L,M)*FB(I,M)
	    PK(I,L,M)=0.
   33 CONTINUE                                                          
      GO TO 15                                                          
   14 CONTINUE                                                          
C     ****     UPPER BNDRY                                              
      DO 17 L=1,2                                                       
	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                                                       
	DO 22 I=1,LEN1
	  WKV1(I,L)=FK(I,L)
   22 CONTINUE                                                          
      DO 21 M=1,2                                                       
	DO 23 I=1,LEN1
	  GAMA(I,K+1,L,M)=0.
	  WKV1(I,L)=WKV1(I,L)-PK(I,L,M)*Z(I,K,M)
   23   CONTINUE
	DO 21 KK=1,2
	  DO 21 I=1,LEN1
	    GAMA(I,K+1,L,M)=GAMA(I,K+1,L,M)+WKM1(I,L,KK)*RK(I,KK,M)
   21 CONTINUE                                                          
C     ****     Z(K+1)=WKM1*WKV1                                         
      DO 231 L=1,2                                                      
	DO 232 I=1,LEN1
	  Z(I,K+1,L)=0.
  232   CONTINUE
	DO 231 M=1,2
	  DO 231 I=1,LEN1
	    Z(I,K+1,L)=Z(I,K+1,L)+WKM1(I,L,M)*WKV1(I,M)
  231 CONTINUE                                                          
    6 CONTINUE                                                          
C     ****     SET PSNP(KMAXP1) TO ZERO                                 
      DO 24 L=1,2                                                       
	NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX
	DO 24 I=1,LEN1
	  F(I,NPSL)=0.
   24 CONTINUE                                                          
C     ****     DOWNWARD SWEEP                                           
      DO 25 KK=1,KMAX                                                   
	K=KMAX+1-KK
	DO 25 L=1,2
	  NPSL=NJNP+NPS+(L-1)*KMAXP1+(K-1)
	  DO 26 I=1,LEN1
	    F(I,NPSL)=Z(I,K+1,L)
   26     CONTINUE
	  DO 25 M=1,2
	    NPSM=NJNP+NPS+(M-1)*KMAXP1+K
	    DO 25 I=1,LEN1
	      F(I,NPSL)=F(I,NPSL)-GAMA(I,K+1,L,M)*F(I,NPSM)
   25 CONTINUE                                                          
C     ****     INSERT VALUE OF PS(KMAXP1) USING BNDRY CONDITION         
      DO 27 L=1,2                                                       
	NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX
	DO 27 I=1,LEN1
	  F(I,NPSL)=(1.+.5*EP(I,L,KP)*C(3))/(1.-.5*EP(I,L,KP)*C(3))*
     A      F(I,NPSL-1)
   27 CONTINUE                                                          
C     ****     TIME SMOOTHING OF PSI
      NPSNMK=NJ+NPSNM-1
      NPSK=NJ+NPS-1
      NPSNPK=NJNP+NPS-1
      NPSMNK=NJNP+NPSNM-1
      DO 29 K = 1,2*KMAXP1
      DO 29 I=1,LEN1
	F(I,NPSMNK+K)=C(30)*F(I,NPSK+K)+C(31)*(F(I,NPSNMK+K)+
     1  F(I,NPSNPK+K))
   29 CONTINUE
C     ****     INSERT PERIODIC POINTS
      K1=NJNP+NPS
      K2=K1+2*KMAXP1-1
      DO 31 N=1,2
	DO 32 I=1,2
	  DO 32 K=K1,K2
	  F(I,K)=F(I+IMAX,K)
	  F(I+IMAXP2,K)=F(I+2,K)
   32   CONTINUE
	K1=NJNP+NPSNM
	K2=K1+2*KMAXP1-1
   31 CONTINUE
C     ****     INSURE NON-NEGATIVE PSI
      NPS1K=NJNP+NPS-1
      NPS2K=NJNP+NPS2-1
      NPS1MK=NJNP+NPSNM-1
      NPS2MK=NJNP+NPS2NM-1
      DO 41 K=1,KMAXP1
      DO 41 I=1,LEN1
	F(I,NPS1K+K)=CVMGP(F(I,NPS1K+K),SMALL,F(I,NPS1K+K)-SMALL)
	F(I,NPS2K+K)=CVMGP(F(I,NPS2K+K),SMALL,F(I,NPS2K+K)-SMALL)
	F(I,NPS1MK+K)=CVMGP(F(I,NPS1MK+K),SMALL,F(I,NPS1MK+K)-SMALL)
	F(I,NPS2MK+K)=CVMGP(F(I,NPS2MK+K),SMALL,F(I,NPS2MK+K)-SMALL)
	F(I,NPS1K+K)=F(I,NPS1K+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1K+K)+
     1    F(I,NPS2K+K)),1.-SMALL-F(I,NPS1K+K)-F(I,NPS2K+K))
	F(I,NPS2K+K)=F(I,NPS2K+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1K+K)+
     1    F(I,NPS2K+K)),1.-SMALL-F(I,NPS1K+K)-F(I,NPS2K+K))
	F(I,NPS1MK+K)=F(I,NPS1MK+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1MK+K)+
     1  F(I,NPS2MK+K)),1.-SMALL-F(I,NPS1MK+K)-F(I,NPS2MK+K))
	F(I,NPS2MK+K)=F(I,NPS2MK+K)*CVMGP(1.,(1.-SMALL)/(F(I,NPS1MK+K)+
     1  F(I,NPS2MK+K)),1.-SMALL-F(I,NPS1MK+K)-F(I,NPS2MK+K))
   41 CONTINUE
      RETURN                                                            
      END                                                               
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)
    9   CONTINUE
C       ****
C       ****     PERIODIC POINTS
C       ****
	NPSK = NJ+NPS-1
	DO 16 I=1,2
	  DO 16 K=1,12*KMAXP1
	    F(I,NPSK+K) = F(I+IMAX,NPSK+K)
	    F(I+IMAXP2,NPSK+K) = F(I+2,NPSK+K)
   16   CONTINUE
	IFP(1) = NPS
	IFP(2) = NPS2
	IFP(3) = NPSA
	IFP(4) = NPSH
	IFP(5) = NPNO
	IFP(6) = NPN4S
	IFP(7) = NPCH4
	IFP(8) = NPH2
	IFP(9) = NPCO2
	IFP(10) = NPCO
	IFP(11) = NPHOX
	IFP(12) = NPH2O
C       CALL PLOT(IFP,12,NJ)
C       ****
C       ****     COPY CURRENT FIELDS TO NM ARRAYS
C       ****
	NPSK = NJ+NPS-1
	NPSNMK = NJ+NPSNM-1
	DO 14 K = 1,12*KMAXP1
	DO 14 I = 1,LEN1
	  F(I,NPSNMK+K) = F(I,NPSK+K)
   14   CONTINUE
	IFP(1) = NPSNM
	IFP(2) = NPS2NM
	IFP(3) = NPSANM
	IFP(4) = NPSHNM
	IFP(5) = NPNONM
	IFP(6) = NN4SNM
	IFP(7) = NCH4NM
	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)
	PS2B(I+IMAXP2) = PS2B(I+2)
	DO 17 K = 1,KMAXP1
	  DIFK(I,K) = DIFK(I+IMAX,K)
	  DIFK(I+IMAXP2,K) = DIFK(I+2,K)
	  FT(K,I) = FT(K,I+IMAX)
	  FT(K,I+IMAXP2) = FT(K,I+2)
	  FUI(K,I) = FUI(K,I+IMAX)
	  FU(K,I) = FU(K,I+IMAX)
	  FUI(K,I+IMAXP2) = FUI(K,I+2)
	  FU(K,I+IMAXP2) = FU(K,I+2)
	  FW(K,I) = FW(K,I+IMAX)
	  FWI(K,I) = FWI(K,I+IMAX)
	  FW(K,I+IMAXP2) = FW(K,I+2)
	  FWI(K,I+IMAXP2) = FWI(K,I+2)
	  FRJ(K,I) = FRJ(K,I+IMAX)
	  FRJ(K,I+IMAXP2) = FRJ(K,I+2)
   17 CONTINUE
C     ****                                                              
C     ****     PERIODIC POINTS FOR SOURCE TERMS
C     ****                                                              
C     CALL CHKBLKS
      DO 18 I = 1,2
	DO 18 K =1,KMAXP1
	  FS11(K,I) = FS11(K,I+IMAX)
	  FS11(K,I+IMAXP2) = FS11(K,I+2)
	  FS12(K,I) = FS12(K,I+IMAX)
	  FS12(K,I+IMAXP2) = FS12(K,I+2)
	  FS21(K,I) = FS21(K,I+IMAX)
	  FS21(K,I+IMAXP2) = FS21(K,I+2)
	  FS22(K,I) = FS22(K,I+IMAX)
	  FS22(K,I+IMAXP2) = FS22(K,I+2)
	  FS1(K,I) = FS1(K,I+IMAX)
	  FS1(K,I+IMAXP2) = FS1(K,I+2)
	  FS2(K,I) = FS2(K,I+IMAX)
	  FS2(K,I+IMAXP2) = FS2(K,I+2)
   18 CONTINUE
C     CALL CHKBLKS
C     ****
C     ****     PERIODIC POINTS RATE COEFFICIENTS
C     ****
      DO 19 I = 1,2
	DO 20 N = 1,25
	  DO 20 K = 1,KMAXP1
	    RKK(N,K,I) = RKK(N,K,I+IMAX)
	    RKK(N,K,I+IMAXP2) = RKK(N,K,I+2)
   20   CONTINUE
	DO 21 N = 1,10
	  DO 21 K = 1,KMAXP1
	    ALPP(N,K,I) = ALPP(N,K,I+IMAX)
	    ALPP(N,K,I+IMAXP2) = ALPP(N,K,I+2)
   21   CONTINUE
	DO 22 N = 1,20
	  DO 22 K = 1,KMAXP1
	    RBB(N,K,I) = RBB(N,K,I+IMAX)
	    RBB(N,K,I+IMAXP2) = RBB(N,K,I+2)
   22   CONTINUE
	DO 23 N = 1,50
	  DO 23 K = 1,KMAXP1
	    RKMM(N,K,I) = RKMM(N,K,I+IMAX)
	    RKMM(N,K,I+IMAXP2) = RKMM(N,K,I+2)
   23   CONTINUE
	DO 24 K = 1,KMAXP1
	  RK12(K,I) = RK12(K,I+IMAX)
	  RK12(K,I+IMAXP2) = RK12(K,I+2)
	  RK13(K,I) = RK13(K,I+IMAX)
	  RK13(K,I+IMAXP2) = RK13(K,I+2)
	  RK14(K,I) = RK14(K,I+IMAX)
	  RK14(K,I+IMAXP2) = RK14(K,I+2)
	  RK15(K,I) = RK15(K,I+IMAX)
	  RK15(K,I+IMAXP2) = RK15(K,I+2)
   24   CONTINUE
   19 CONTINUE
C     ****
C     ****     NTRANSFER T AND U
C     ****     N
      NTK=NJ+NT-1
      NUK=NJ+NU-1
      DO 5 K=1,KMAX                                                     
	DO 5 I=1,LEN1
	  F(I,NTK+K)=.5*(FT(K,I)+FT(K+1,I))
	  F(I,NUK+K)=.5*(FU(K,I)+FU(K+1,I))
    5 CONTINUE                                                          
C     ****     BOUNDARY                                                 
      DO 10 I = 1,LEN1
	F(I,NTK+KMAXP1)=FT(1,I)
	F(I,NUK+KMAXP1)=FU(1,I)
   10 CONTINUE
C     ****                                                              
C     ****     TRANSFER W AND RJ                                        
C     ****                                                              
      NWK=NJ+NW-1                                                       
      NRJK=NDJ+NRJ-1                                                    
      DO 6 K=1,KMAXP1                                                   
	DO 6 I=1,LEN1
	  F(I,NWK+K)=FW(K,I)
	  F(I,NRJK+K)=FRJ(K,I)
    6 CONTINUE                                                          
C     ****                                                              
C     ****     UPDATE PSI'S                                             
C     ****                                                              
C     CALL CHKBLKS
      CALL RATES(FT)                                                    
C     CALL CHKBLKS
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                                                       
C     CALL CHKBLKS
      CALL CMPN4S                                                       
C     CALL CHKBLKS
      CALL COMPNO                                                       
C     CALL CHKBLKS
      CALL COMPARG                                                      
C     CALL CHKBLKS
C     CALL COMPHE
C     CALL CHKBLKS
      CALL COMPCH4                                                      
C     CALL CHKBLKS
      CALL COMPH2                                                       
C     CALL CHKBLKS
      CALL COMPCO2                                                      
C     CALL CHKBLKS
      CALL COMPCO                                                       
C     CALL CHKBLKS
      CALL COMPHOX
C     CALL CHKBLKS
      CALL COMPH2O
C     CALL CHKBLKS
C     ****
C     ****     NTRANSFER UI
C     ****     N
      NUK=NJ+NU-1
      DO 55 K=1,KMAX
	DO 55 I=1,LEN1
	  F(I,NUK+K)=.5*(FUI(K,I)+FUI(K+1,I))
   55 CONTINUE
C     ****     BOUNDARY                                                 
      DO 67 I = 1,LEN1
	F(I,NUK+KMAXP1)=FUI(1,I)
   67 CONTINUE
C     ****                                                              
C     ****     TRANSFER WI
C     ****                                                              
      NWK=NJ+NW-1                                                       
      DO 66 K=1,KMAXP1
	DO 66 I=1,LEN1
	  F(I,NWK+K)=FWI(K,I)
   66 CONTINUE
      CALL COMPHE                                                       
	IFP(1) = NPS
	IFP(2) = NPS2
	IFP(3) = NPSA
	IFP(4) = NPSH
	IFP(5) = NPNO
	IFP(6) = NPN4S
	IFP(7) = NPCH4
	IFP(8) = NPH2
	IFP(9) = NPCO2
	IFP(10) = NPCO
	IFP(11) = NPHOX
	IFP(12) = NPH2O
C       CALL PLOT(IFP,12,NJNP)
C     WRITE(6,1000)PNOB,PN4SB
C1000 FORMAT(* PNOB, PN4SB *,2E13.4)                                    
C     ****                                                              
C     ****     CONVERT PSI'S TO NUMBER DENSITIES                        
C     ****                                                              
C     ****     WK1=INTERPOLATED PSI1,   WK2=INTERPOLATED PSI2           
C     ****     WK3=INTERPOLATED PSAR,   WK4=INTERPOLATED PSHE           
C     ****     WK5=INTERPOLATED PSNO,   WK6=INTERPOLATED PSN4S          
C     ****     WK7=INTERPOLATED PSCH4,  WK8=INTERPOLATED PSH2           
C     ****     WK9=INTERPOLATED PSCO2,  WK10=INTERPOLATED PSCO          
C     ****     WK11=INTERPOLATED PSHOX, WK12=INTERPOLATED PSH2O         
C     ****                                                              
      NPS1K=NJNP+NPS-1
      NPS2K=NPS1K+KMAXP1                                                
      NPSAK=NJNP+NPSA-1
      NPSHK=NJNP+NPSH-1
      NPNOK=NJNP+NPNO-1
      NPN4SK=NJNP+NPN4S-1
      NPCH4K=NJNP+NPCH4-1
      NPH2K=NJNP+NPH2-1
      NPCO2K=NJNP+NPCO2-1
      NPCOK=NJNP+NPCO-1
      NPHOXK=NJNP+NPHOX-1
      NPH2OK=NJNP+NPH2O-1
      DO 7 K=1,KMAX                                                     
	DO 7 I=1,LEN1
	  WK1(K+1,I)=.5*(F(I,NPS1K+K)+F(I,NPS1K+K+1))
	  WK2(K+1,I)=.5*(F(I,NPS2K+K)+F(I,NPS2K+K+1))
	  WK3(K+1,I)=.5*(F(I,NPSAK+K)+F(I,NPSAK+K+1))
	  WK4(K+1,I)=.5*(F(I,NPSHK+K)+F(I,NPSHK+K+1))
	  WK5(K+1,I)=.5*(F(I,NPNOK+K)+F(I,NPNOK+K+1))
	  WK6(K+1,I)=.5*(F(I,NPN4SK+K)+F(I,NPN4SK+K+1))
	  WK7(K+1,I)=.5*(F(I,NPCH4K+K)+F(I,NPCH4K+K+1))
	  WK8(K+1,I)=.5*(F(I,NPH2K+K)+F(I,NPH2K+K+1))
	  WK9(K+1,I)=.5*(F(I,NPCO2K+K)+F(I,NPCO2K+K+1))
	  WK10(K+1,I)=.5*(F(I,NPCOK+K)+F(I,NPCOK+K+1))
	  WK11(K+1,I)=.5*(F(I,NPHOXK+K)+F(I,NPHOXK+K+1))
	  WK12(K+1,I)=.5*(F(I,NPH2OK+K)+F(I,NPH2OK+K+1))
    7 CONTINUE                                                          
C     ****     LOWER BOUNDARY                                           
      DO 13 I = 1,LEN1
	WK1(1,I)=PS1B
	WK2(1,I)=PS2B(I)
	WK3(1,I)=.5*(F(1,NPSAK+1)+PARB(I))
	WK4(1,I)=.5*(F(1,NPSHK+1)+PHB(I))
	WK5(1,I)=.5*(F(1,NPNOK+1)+PNOB(I))
	WK6(1,I)=.5*(F(1,NPN4SK+1)+PN4SB(I))
	WK7(1,I)=.5*(F(1,NPCH4K+1)+PCH4B(I))
	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                                                               
*/
*DK COMPCO
*/
      SUBROUTINE COMPCO                                                 
      SAVE
C     ****                                                              
C     ****     ADVANCE CO COMPOSITION BY ONE TIME STEP                  
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CH4H2COF, CONS, FIELDS, FIELDT, FLUX,
C     ****         INDEX, PHOTO, PSIB, RATEBLK, 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 PHICO(3)                                                
      DATA RMCO,PHICO/28.,0.833,1.427,0.852/                            
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)=-(RKMM(40,K,I)*FNO(K,I)*(FNO(K,I)+FNO2(K,I)+FNN2(K,I))
     1      +RKMM(41,K,I)*FNOH(K,I))
	  S2(I,K)=RKCO2P(K,I)*FNCO2(K,I)+(RKMM(35,K,I)*FNOH(K,I)
     1      +RKMM(36,K,I)*FNO(K,I)+RKMM(37,K,I)*FNO1D(K,I))*FNCH4(K,I)+
     2      SCOI(K,I)
    2 CONTINUE                                                          
C     ****    IF IBNDCO = 0, FLUX GIVEN AT LOWER BOUNDARY               
      IF(IBNDCO.EQ.0)THEN                                               
	DO 3 I=3,LEN1-2
          T1(I)=0.                                                      
          T2(I)=0.                                                      
	  T3(I)=FLUXCO(I)
    3   CONTINUE
	IBNDB=1
C       ****    IF IBNDCO = 1, PHOTOCHEMICAL EQUILIBRIUM AT LOWER
C       ****      BOUNDARY
      ELSE IF(IBNDCO.EQ.1)THEN                                          
	DO 4 I=3,LEN1-2
          T1(I)=0.                                                      
          T2(I)=1.                                                      
	  T3(I)=-XCO(I)*RMCO/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+
     1      FNN2(1,I)*RMASS(3))
    4   CONTINUE
C     ****
C     ****     PERIODIC POINTS
C     ****
      DO 5 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 5 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)
    5 CONTINUE
      IBNDB=0
      ENDIF                                                             
      IBND=0                                                            
      ALFA=0.                                                           
      CALL MINOR(NPCO,NPCONM,RMCO,PHICO,ALFA,IBND,IBNDB,PCOB)           
      RETURN                                                            
      END                                                               
*/
*DK COMPCO2
*/
      SUBROUTINE COMPCO2                                                
      SAVE
C     ****                                                              
C     ****     ADVANCE CO2 COMPOSITION BY ONE TIME STEP                 
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CH4H2COF, CONS, FIELDS, FIELDT, FLUX,
C     ****         INDEX, PHOTO, PSIB, RATEBLK, 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 PHICO2(3)                                               
      DATA RMCO2,PHICO2/44.,1.199,3.91,1.201/                           
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,IMAXP4-2
	  S1(I,K)=-(RJCO2T(K,I)+RSKIONS(K,I))
	  S2(I,K)=RKMM(40,K,I)*FNCO(K,I)*FNO(K,I)*(FNO(K,I)+FNO2(K,I)+
     1      FNN2(K,I))+RKMM(41,K,I)*FNCO(K,I)*FNOH(K,I)+SCO2I(K,I)
    2 CONTINUE                                                          
C     ****    DENSITY SPECIFIED AT LOWER BOUNDARY                       
      DO 4 I=3,LEN1-2
	T1(I)=0.
	T2(I)=1.
	T3(I)=-XCO2(I)*RMCO2/(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(NPCO2,NCO2NM,RMCO2,PHICO2,ALFA,IBND,IBNDB,PCO2B)       
      RETURN                                                            
      END                                                               
*/
*DK COMPH2
*/
      SUBROUTINE COMPH2                                                 
      SAVE
C     ****                                                              
C     ****     ADVANCE H2 COMPOSITION BY ONE TIME STEP                  
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CH4H2COF, CONS, FIELDS, FIELDT, FLUX,
C     ****         H2FORG, INDEX, PHOTO, PSIB, RATEBLK, VSCR
C     ****                                                              
*CA PARAMS
*CA ASAS2
*CA CH4H2CO2
*CA CONS2
*CA FIELDS2
*CA FIELDT2
*CA FLUX2
*CA H2FORG2
*CA INDEX
*CA PHOTO2
*CA PSIB2
*CA RATEBLK2
*CA VSCR
      DIMENSION PHIH2(3)                                                
      DATA RMH2,PHIH2/2.,0.226,0.321,0.282/                             
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,IMAXP4-2
	  S1(I,K)=-(RKMM(11,K,I)*FNO1D(K,I)+RKMM(20,K,I)*FNO(K,I)+
     1      RKMM(25,K,I)*FNOH(K,I)+FK9OP(K,I))
	  S2(I,K)=0.12*RJH2OLY(K,I)*FNH2O(K,I)+RKMM(30,K,I)*FNH(K,I)*
     1      FNHO2(K,I)+RKMM(33,K,I)*
     2      FNH(K,I)*FNH(K,I)*(FNO(K,I)+FNO2(K,I)+FNN2(K,I))+RJCH4(K,I)*
     3      FNCH4(K,I)
    2 CONTINUE                                                          
C     ****    DENSITY SPECIFIED AT LOWER BOUNDARY                       
      DO 4 I=3,LEN1-2
	T1(I)=0.
	T2(I)=1.
	T3(I)=-XH2(I)*RMH2/(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(NPH2,NPH2NM,RMH2,PHIH2,ALFA,IBND,IBNDB,PH2B)           
      RETURN                                                            
      END                                                               
*/
*DK COMPH2O
*/
      SUBROUTINE COMPH2O                                                
      SAVE
C     ****                                                              
C     ****     ADVANCE H2O COMPOSITION BY ONE TIME STEP                 
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CONS, DISSHOX, FIELDS, FIELDT, FLUX,
C     ****         INDEX, IONHOX, PHOTO, PSIB, RATEBLK, VSCR
C     ****                                                              
C     ****                                                              
*CA PARAMS
*CA ASAS2
*CA CONS2
*CA DISSHOX2
*CA FIELDS2
*CA FIELDT2
*CA FLUX2
*CA INDEX
*CA IONHOX2
*CA PHOTO2
*CA PSIB2
*CA RATEBLK2
*CA VSCR
      DIMENSION PHIH2O(3)                                               
      DATA RMH2O,PHIH2O/18.,0.817,0.922,0.920/                          
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-1
	  S1(I,K)=-(SH2OSRB(K,I)+SH2OSRC(K,I)+SH2OLYA(K,I)+SH2OEUV(K,I)
     1      +RKMM(10,K,I)*FNO1D(K,I)+RKK(10,K,I)*YIOP(K,I)+RKMM(42,K,I)*
     2      FNO(K,I))
	  S2(I,K)=RKMM(22,K,I)*FNOH(K,I)*FNOH(K,I)+RKMM(23,K,I)*
     1      FNOH(K,I)*FNHO2(K,I)+RKMM(24,K,I)*FNOH(K,I)*FNH2O2(K,I)+
     2      RKMM(25,K,I)*FNOH(K,I)*FNH2(K,I)+RKMM(32,K,I)*FNH(K,I)*
     3      FNHO2(K,I)+(2.*RKMM(35,K,I)*FNOH(K,I)+RKMM(36,K,I)*FNO(K,I)+
     4      RKMM(37,K,I)*FNO1D(K,I))*FNCH4(K,I)-PHOXIC(K,I)*0.5
    2 CONTINUE                                                          
C     ****    DENSITY SPECIFIED AT LOWER BOUNDAR                        
      DO 4 I=3,LEN1-2
	T1(I)=0.
	T2(I)=1.
	T3(I)=-XH2O(I)*RMH2O/(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
C     CALL CONREC(SH2OSRB(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(SH2OSRC(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(SH2OLYA(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(SH2OEUV(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(FNO1D(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(YIOP(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(FNOH(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(FNH2O2(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(FNH(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(PHOXIC(1,3),KMAXP1,KMAXP1,IMAX+1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(FNO2(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(S1,LEN1,LEN1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S2,LEN1,LEN1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
      DO 5 I = 3,IMAX+3
	DO 5 K = 1,KMAXP1
	  S3(I,K) = RKK(10,K,I)
	  S4(I,K) = RKMM(10,K,I)
	  S5(I,K) = RKMM(22,K,I)
	  S6(I,K) = RKMM(23,K,I)
	  S7(I,K) = RKMM(24,K,I)
	  S8(I,K) = RKMM(25,K,I)
	  S9(I,K) = RKMM(32,K,I)
	  S10(I,K) = RKMM(35,K,I)
	  S11(I,K) = RKMM(36,K,I)
	  S12(I,K) = RKMM(37,K,I)
	  S13(I,K) = RKMM(42,K,I)
    5 CONTINUE
C     CALL CONREC(S3(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S4(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S5(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S6(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S7(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S8(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S9(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S10(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S11(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S12(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     CALL CONREC(S13(3,1),LEN1,IMAX+1,KMAXP1,0.,0.,0.,0,0,-1430B)
C     CALL FRAME
C     WRITE(6,2345)(T3(I),XH2O(I),I=1,LEN1)
C2345 FORMAT(2E15.4)
      IBNDB=0
      IBND=0                                                            
      ALFA=0.                                                           
      CALL MINOR(NPH2O,NH2ONM,RMH2O,PHIH2O,ALFA,IBND,IBNDB,PH2OB)       
      RETURN                                                            
      END                                                               
*/
*DK COMPHE
*/
      SUBROUTINE COMPHE                                                 
      SAVE
C     ****
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CONS, FIELDS, HEFLUX, INDEX, PHOTOO,
C     ****         PSIB, VSCR
C     ****
*CA PARAMS
*CA ASAS2
*CA CONS2
*CA FIELDS2
*CA HEFLUX2
*CA INDEX
*CA PHOTOO2
*CA PSIB2
*CA VSCR
      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 PHIH(3)                                                 
      DIMENSION FMM(IMXP)
C     DATA RMH,PHIH,ALFAH,PSHB/23.,.270,.404,.322,-0.38,0.73E-6/
      DATA RMH,PHIH,ALFAH,PSHB/23.,1.042,1.509,1.176,0.17,0.131E-1/
C     ****     UPPER BOUNDARY CONDITION                                 
C     ****     T4=TOTAL FLUX AT UPPER BOUNDARY                          
      DO 10 I=3,LEN1-2
      FMM(I)=FNO(1,I)*16.+FNO2(1,I)*32.+FNN2(1,I)*28.
      T4(I) =0.
      T1(I)=0.
      T2(I)=1.
      T3(I)=-XLBNAI(I)*RMH/FMM(I)
   10 CONTINUE
C     DO 10 I=3,LEN1-2
C     FMM(I)=FNO(1,I)*16.+FNO2(1,I)*32.+FNN2(1,I)*28.
C     T4(I) = FLUX(I)
C  10 CONTINUE
C     ****     LOWER BOUNDARY CONDITION - CONSTANT VALUE                
C     DO 4 I=3,LEN1-2
C     T1(I)=0.
C     T2(I)=1.
C     T3(I)=-XHE(I)*RMH/FMM(I)
C   4 CONTINUE
C     ****     SOURCES ARE ZERO                                         
      DO 5 K=1,KMAXP1
      DO 5 I=3,IMAXP4-2
      S1(I,K)=-XLNAI(K,I)
      S2(I,K)=PNAI(K,I)
    5 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)
    3 CONTINUE
      IBND=0
      IBNDB=0                                                           
      CALL MINOR(NPSH,NPSHNM,RMH,PHIH,ALFAH,IBND,IBNDB,PHB)             
      RETURN                                                            
      END                                                               
*/
*DK COMPARG
*/
      SUBROUTINE COMPARG                                                
      SAVE
C     ****     ADVANCE AR COMPOSITION BY ONE TIME STEP                  
C     ****
C     ****   COMMON DECKS:
C     ****   PARAMS,ASAS,CONS,FIELDS,INDEX,PHOTOO,PSIB,VSCR
C     ****
*CA PARAMS
*CA ASAS2
*CA CONS2
*CA FIELDS2
*CA INDEX
*CA PHOTOO2
*CA PSIB2
*CA VSCR
      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 PHIA(3)                                                 
      DIMENSION FMM(IMXP)
      DATA RMA,PHIA,ALFAA,PSAB/23.,1.042,1.509,1.176,0.17,0.131E-1/
C     ****     UPPER BOUNDARY - ZERO DIFFUSIVE FLUX                     
C     ****     LOWER BOUNDARY - CONSTANT VALUE                          
      DO 1 I=3,LEN1-2
      FMM(I)=FNO(1,I)*16.+FNO2(1,I)*32.+FNN2(1,I)*28.
      T4(I)=0.
      T1(I)=0
      T2(I)=1.
      T3(I)=-XLBNA(I)*RMA/FMM(I)
    1 CONTINUE
      IBNDB=0
C     LOWER BOUNDARY FLUX CONDITION
C     FLUXNA=0.
C     DO 1 I=3,LEN1-2
C     T4(I)=0.
C     T1(I)=0
C     T2(I)=0.
C     T3(I)=FLUXNA
C   1 CONTINUE
C     IBNDB=1
C     ****     SOURCES
      DO 2 K=1,KMAXP1                                                   
	DO 2 I=3,IMAXP4-2
	  S1(I,K)=-XLNA(K,I)
	  S2(I,K)=PNA(K,I)
    2 CONTINUE                                                          
      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
      IBND=0                                                            
      CALL MINOR (NPSA,NPSANM,RMA,PHIA,ALFAA,IBND,IBNDB,PARB)           
      RETURN                                                            
      END                                                               
*/
*DK COMPHOX
*/
      SUBROUTINE COMPHOX
      SAVE
C     ****                                                              
C     ****     ADVANCE HOX COMPOSITION BY ONE TIME STEP                 
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CONS, DISSHOX, FIELDS, FIELDT, FLUX,
C     ****         HOXSTUF, HOXUPP, INDEX, IONHOX, NEWRAT, PHOTO, PSIB,
C     ****         RATEBLK, RJH2O2O, VSCR
C     ****                                                              
*CA PARAMS
*CA ASAS2
*CA CONS2
*CA DISSHOX2
*CA FIELDS2
*CA FIELDT2
*CA FLUX2
*CA HOXSTUF2
*CA HOXUPP2
*CA INDEX
*CA IONHOX2
*CA NEWRAT2
*CA PHOTO2
*CA PSIB2
*CA RATEBLK2
*CA RJH2O2O2
*CA VSCR
      DIMENSION FM(IMXP),SS1(IMXP),SS2(IMXP),RR(IMXP),XHOX(IMXP)
      DIMENSION PHIHOX(3)
      DATA RMHOX,PHIHOX/1.,0.146,0.243,0.162/                           
C     ****                                                              
C     ****     UPPER BOUNDARY: TOTAL FLUX=FHOX                          
C     ****                                                              
      DO 1 I=3,LEN1-2
	T4(I)=FHOX(I)
    1 CONTINUE                                                          
      IBND=1                                                            
C     ****                                                              
C     ****     SOURCES                                                  
C     ****                                                              
C     U=9.
C     V=17.
C     W=25.
      U=1.                                                              
      V=1.                                                              
      W=1.                                                              
      DO 2 K=1,KMAXP1                                                   
	DO 2 I=3,LEN1-2
	  FM(I)=FNO(K,I)+FNO2(K,I)+FNN2(K,I)
	  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))
	  SS1(I)=HO2HR(K,I)
	  SS2(I)=OHHR(K,I)
	  RR(I)=HHOXR(K,I)
	  S2(I,K)=(2.*(SH2OSRB(K,I)+SH2OSRC(K,I)+SH2OEUV(K,I)+0.88*
     1      SH2OLYA(K,I))*FNH2O(K,I)*U+2.*RJH2O2(K,I)*FNH2O2(K,I)*V+
     2      PHOXIC(K,I)+(2.*RKMM(11,K,I)*FNH2(K,I)*U+2.*RKMM(10,K,I)*
     3      FNH2O(K,I)*V+2.*RKMM(37,K,I)*FNCH4(K,I))*FNO1D(K,I)*V+
     4      RKK(11,K,I)*YIHP(K,I)*FNO(K,I)+2.*RKMM(20,K,I)* FNH2(K,I)*
     5      FNO(K,I)*U+2.*RKMM(19,K,I)*FNH2O2(K,I)*FNO(K,I)*W+2.*
     6      RKMM(42,K,I)*FNH2O(K,I)*FNO(K,I)*V+2.*RKMM(36,K,I)*
     7      FNCH4(K,I)*FNO(K,I)*V+RKK(9,K,I)*FNH2(K,I)*YIOP(K,I)+
     8      ALPP(4,K,I)*YIOHP(K,I)*FNE(K,I))
	  S1(I,K)=-(2.*RKMM(23,K,I)*SS1(I)*RR(I)*SS2(I)*RR(I)+
     1      2.*RKMM(27,K,I)*SS1(I)*SS1(I)*RR(I)*RR(I)+2.*(RKMM(30,K,I)+
     2      RKMM(32,K,I))*SS1(I)*RR(I)*RR(I)+2.*RKMM(22,K,I)*SS2(I)*
     3      RR(I)*SS2(I)*RR(I)+2.*RKMM(33,K,I)*FM(I)*RR(I)*RR(I))*
     4      FNHOX(K,I)-RKK(12,K,I)*YIOP(K,I)*RR(I)
    2 CONTINUE                                                          
C     ****    PHOTOCHEMICAL EQUILIBRIUM AT LOWER BOUNDARY               
      DO 3 I = 3,LEN1-2
	FM(I)=FNO(1,I)+FNO2(1,I)+FNN2(1,I)
	SS1(I)=HO2HR(1,I)
	SS2(I)=OHHR(1,I)
	RR(I)=HHOXR(1,I)
	XHOX(I)=SQRT((2.*(SH2OSRB(1,I)+SH2OSRC(1,I)+SH2OEUV(1,I)+0.88*
     1    SH2OLYA(1,I))*FNH2O(1,I)*U+2.*RJH2O2(1,I)*FNH2O2(1,I)*V+
     1    PHOXIC(1,I)
     2    +(2.*RKMM(11,1,I)*FNH2(1,I)*U+2.*RKMM(10,1,I)*FNH2O(1,I)*V+2.*
     3    RKMM(37,1,I)*FNCH4(1,I))*FNO1D(1,I)*V+RKK(11,1,I)*YIHP(1,I)*
     4    FNO(1,I)+2.*RKMM(20,1,I)*FNH2(1,I)*FNO(1,I)*U+2.*RKMM(19,1,I)*
     5    FNH2O2(1,I)*FNO(1,I)*W+2.*RKMM(42,1,I)*FNH2O(1,I)*FNO(1,I)*V+
     6    2.*RKMM(36,1,I)*FNCH4(1,I)*FNO(1,I)*V+RKK(9,1,I)*FNH2(1,I)*
     7    YIOP(1,I)+ALPP(4,1,I)*YIOHP(1,I)*FNE(1,I))/(2.*RKMM(23,1,I)*
     8    SS1(I)*RR(I)*SS2(I)*RR(I)+2.*RKMM(27,1,I)*SS1(I)*SS1(I)*RR(I)*
     9    RR(I)+2.*(RKMM(30,1,I)+RKMM(32,1,I))*SS1(I)*RR(I)*RR(I)+2.*
     1    RKMM(22,1,I)*RR(I)*RR(I)*SS2(I)*SS2(I)+2.*RKMM(33,1,I)*FM(I)*
     2    RR(I)*RR(I)))
    3 CONTINUE
      DO 4 I=3,LEN1-2
	T1(I)=0.
	T2(I)=1.
	T3(I)=-XHOX(I)/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+
     1    FNN2(1,I)*RMASS(3))
    4 CONTINUE
C     ****
C     ****     PERIODIC POINTS
C     ****
      DO 5 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 5 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)
    5 CONTINUE
      IBNDB=0
      ALFA=-0.38                                                        
      CALL MINOR(NPHOX,NHOXNM,RMHOX,PHIHOX,ALFA,IBND,IBNDB,PHOXB)       
      RETURN                                                            
      END                                                               
*/
*DK COMPNO
*/
      SUBROUTINE COMPNO
      SAVE
C     ****                                                              
C     ****     ADVANCE NO COMPOSITION BY ONE TIME STEP                  
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CONS, FIELDS, FIELDT, FLUX, INDEX, PHOTO,
C     ****          PSIB, RATEBLK, VSCR
C     ****                                                              
*CA PARAMS
*CA ASAS2
*CA CONS2
*CA FIELDS2
*CA FIELDT2
*CA FLUX2
*CA INDEX
*CA NOZNOZ2
*CA PHOTO2
*CA PSIB2
*CA RATEBLK2
*CA VSCR
      DIMENSION PHINO(3)                                                
      DATA RMNO,PHINO/30.,0.814,0.866,0.926/                            
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)=-(FK5O2P(K,I)+RBB(3,K,I)*FNN4S(K,I)+RBB(6,K,I)*
     1      FNN2D(K,I)+RBB(8,K,I)+RBB(9,K,I))*RNONOZ(K,I)
	  S2(I,K)=RBB(1,K,I)*FNO2(K,I)*FNN4S(K,I)+RBB(2,K,I)*FNO2(K,I)*
     1      FNN2D(K,I)+RBB(10,K,I)*FNOH(K,I)*FNN4S(K,I)
    2 CONTINUE
C     ****    IF IBNDNO = 0, FLUX GIVEN AT LOWER BOUNDARY               
      IF(IBNDNO.EQ.0)THEN                                               
	DO 3 I=3,LEN1-2
          T1(I)=0.                                                      
          T2(I)=0.                                                      
	  T3(I)=FLUXNO(I)
    3   CONTINUE
	IBNDB=1
C       ****    IF IBNDNO = 1, PHOTOCHEMICAL EQUILIBRIUM AT LOWER
C       ****      BOUNDARY
      ELSE IF(IBNDNO.EQ.1)THEN                                          
	DO 4 I=3,LEN1-2
          T4(I)=0.                                                      
          T1(I)=0.                                                      
          T2(I)=1.                                                      
	  T3(I)=-XNO(I)*RMNO/(FNO2(1,I)*RMASS(1)+FNO(1,I)*RMASS(2)+
     1    FNN2(1,I)*RMASS(3))
    4   CONTINUE
	IBNDB=0
      ENDIF                                                             
      DO 5 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 5 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)
    5 CONTINUE
      IBND=0                                                            
      ALFA=0.                                                           
      CALL MINOR(NPNO,NPNONM,RMNO,PHINO,ALFA,IBND,IBNDB,PNOB)           
      RETURN                                                            
      END                                                               
*/
*DK FFACE
*/
      SUBROUTINE FFACE(II)
      SAVE
C     ****
C     ****     INTERFACES ONE- AND TWO-DIMENSIONAL COMMON BLOCKS
C     ****
*CA PARAMS
*CA CONS2
*CA ASASF
*CA ASAS2F
*CA CH4H2CFF
*CA CH4H2C2F
*CA CRATESF
*CA CRATES2F
*CA DISSHOXF
*CA DISSHX2F
*CA FIELDSF
*CA FIELDS2F
*CA FIELDTF
*CA FIELDT2F
*CA FLUX
*CA FLUX2F
*CA H2FORG
*CA H2FORG2F
*CA HEFLUX
*CA HEFLUX2F
*CA HOXSTUFF
*CA HOXSTF2F
*CA HOXUPP
*CA HOXUPP2F
*CA IONHOX
*CA IONHOX2F
*CA NEWRATF
*CA NEWRAT2F
*CA OPION
*CA OPION2F
*CA PSIBF
*CA PSIB2F
*CA PHOTOF
*CA PHOTO2F
*CA PHOTOOF
*CA PHOTOO2F
*CA RATEBLK
*CA RATBLK2F
*CA RJH2O2O
*CA RJH2O22F
*CA SOURCEF
*CA SOURCE2F
C     ****
C     ****     ALLOW FOR PERIODIC WRAP AROUND POINTS (TWO AT EACH END)
C     ****
      I = II+2
C     ****
C     ****     /ASAS/ AND /ASAS2/ ***** NEED TO WORK OUT LEAPFROG
C     ****       POINTERS
C     ****
C     ****     /CH4H2COF/ AND /CH4H2CO2/
C     ****
      DO 2 J = 1,7
	DO 2 K = 1,KMAXP1
	  RJCH4F(K,I,J) = RJCH4(K,J)
    2 CONTINUE
C     ****
C     ****     /CRATES/ AND /CRATES2/
C     ****
      DO 3 J = 1,4
	DO 3 K = 1,KMAXP1
	  RK12F(K,I,J) = RK12(K,J)
    3 CONTINUE
C     ****
C     ****     /DISSHOX/ AND /DISSHOX2/
C     ****
      DO 4 J = 1,9
	DO 4 K = 1,KMAXP1
	  SH2OTF(K,I,J) = SH2OT(K,J)
    4 CONTINUE
C     ****
C     ****     /FIELDS/ AND /FIELDS2/
C     ****
      DO 5 J = 1,17
	DO 5 K = 1,KMAXP1
	  FTF(K,I,J) = FT(K,J)
    5 CONTINUE
      F107F = F107
      DO 6 J = 1,3
	DO 6 K = 1,KMAXP1
	  FTEF(K,I,J) = FTE(K,J)
    6 CONTINUE
C     ****
C     ****     /FIELDT/ AND /FIELDT2/
C     ****
      DO 7 J = 1,12
	DO 7 K = 1,KMAXP1
	  FNH2OF(K,I,J) = FNH2O(K,J)
    7 CONTINUE
C     ****
C     ****     /FLUX/ AND /FLUX2/
C     ****
      FLUXNOF(I) = FLUXNO
      IBNDNOF = IBNDNO
      FLUXCOF(I) = FLUXCO
      IBNDCOF = IBNDCO
C     ****
C     ****     /H2FORG/ AND /H2FORG2/
C     ****
      DO 8 K =1,KMAXP1
	RJH2OLYF(K,I) = RJH2OLY(K)
    8 CONTINUE
C     ****
C     ****     /HEFLUX/ AND /HEFLUX2/
C     ****
      FLUXF(I) = FLUX
C     ****
C     ****     /HOXSTUF/ AND /HOXSTUF2/
C     ****
      DO 9 J = 1,10
	DO 9 K = 1,KMAXP1
	  HO2HRF(K,I,J) = HO2HR(K,J)
    9 CONTINUE
C     ****
C     ****     /HOXUPP/ AND /HOXUPP2/
C     ****
      FHOXF(I) = FHOX
C     ****
C     ****     /IONHOX/ AND /IONHOX2/
C     ****
      DO 10 K = 1,KMAXP1
	PHOXICF(K,I) = PHOXIC(K)
   10 CONTINUE
C     ****
C     ****     /NEWRATF/ AND /NEWRAT2F/
C     ****
      DO 11 J = 1,3
	DO 11 K = 1,KMAXP1
	  HO2OHF(K,I,J) = HO2OH(K,J)
   11 CONTINUE
C     ****
C     ****     /OPION/ AND /OPION2/
C     ****
      DO 12 K = 1,KMAXP1
	FIOPF(K,I) = FIOP(K)
   12 CONTINUE
C     ****
C     ****     /PSIB/ AND /PSIB2/
C     ****
      DO 13 J = 1,10
	PARBF(I,J) = PARB(J)
   13 CONTINUE
C     ****
C     ****     /PHOTO/ AND /PHOTO2/
C     ****
      DO 14 J =1,7
	XNOF(I,J) = XNO(J)
   14 CONTINUE
C     ****
C     ****     /PHOTOO/ AND /PHOTOO2/
C     ****
      DO 15 J = 1,2
	XHEF(I,J) = XHE(J)
   15 CONTINUE
C     ****
C     ****     /RATEBLK/ AND /RATEBLK2/
C     ****
      DO 16 J = 1,25
	DO 16 K = 1,KMAXP1
	  RKKF(J,K,I) = RKK(J,K)
   16 CONTINUE
      DO 17 J = 1,10
	DO 17 K = 1,KMAXP1
	  ALPPF(J,K,I) = ALPP(J,K)
   17 CONTINUE
      DO 18 J = 1,20
	DO 18 K = 1,KMAXP1
	  RBBF(J,K,I) = RBB(J,K)
   18 CONTINUE
      DO 19 J = 1,50
	DO 19 K = 1,KMAXP1
	  RKMMF(J,K,I) = RKMM(J,K)
   19 CONTINUE
C     ****
C     ****     /RJH2O2O/ AND /RJH2O2O2/
C     ****
      DO 21 K = 1,KMAXP1
	RJH2O2F(K,I) = RJH2O2(K)
   21 CONTINUE
C     ****
C     ****     /SOURCE/ AND /SOURCE2/
C     ****
      DO 20 J = 1,6
	DO 20 K = 1,KMAXP1
	  FS11F(K,I,J) = FS11(K,J)
   20 CONTINUE
      RETURN
      END
*/
*DK MINOR
*/
      SUBROUTINE MINOR(NX,NXNM,RMX,PHIX,ALFAX,IBND,IBNDB,PXB)           
      SAVE
C     ****
C     ****     COMMON DECKS:
C     ****       PARAMS, ASAS, CONS, INDEX, VSCR,
C     ****
C     ****                                                              
C     ****     ADVANCES MINOR SPECIES BY ONE TIME STEP                  
C     ****     NX = INDEX OF MINOR SPECIES FOR TIME T(N)                
C     ****     NXNM = INDEX OF MINOR SPECIES FOR TIME T(N-1)            
C     ****     NPSDH = INDEX FOR HORIZONTAL DIFFUSION FIELD             
C     ****     RMX = MOLECULAR WEIGHT OF MINOR SPECIES                  
C     ****     PHIX = DIFFUSION VECTOR, PHI(X,N),N=1,3                  
C     ****     ALFAX = THERMAL DIFFUSION COEFFICIENT                    
C     ****     IBND=0 IF DIFFUSIVE FLUX GIVEN AT UPPER BOUNDARY (IN T4) 
C     ****         =1 IF TOTAL FLUX GIVEN AT UPPER BOUNDARY (IN T4)     
C     ****         NOTE: FLUXES ARE IN 1/(CM**2*SEC)                    
C     ****                                                              
C     ****     ON ENTRY:                                                
C     ****     LOWER BOUNDARY:                                          
C     ****        IF IBNDB = 0, THEN                                    
C     ****        T1 = A, T2 = B, T3 = C DEFINE LOWER BOUNDARY CONDITION
C     ****        WHERE       A*DPSX/DZ + B*PSX + C = 0.                
C     ****        IF IBNDB = 1, THEN                                    
C     ****        T3 = UPWARD FLUX (1/(SEC*CM**2)) AT LOWER BOUNDARY    
C     ****     UPPER BOUNDARY:                                          
C     ****        T4 = PHI, GIVES UPWARD FLUX OF SPECIES X AT UPPER
C     ****             BOUNDARY (1/(SEC*CM**2))
C     ****     SOURCES:                                                 
C     ****        S1 = SX/N(X), WHERE SX IS PORTION OF NUMBER DENSITY   
C     ****        SOURCE PROPORTIONAL TO N(X), THE MINOR SPECIES NUMBER 
C     ****        DENSITY                                               
C     ****        S2 = S0, PORTION OF NUMBER DENSITY SOURCE INDEPENDENT 
C     ****        OF N(X).                                              
*CA PARAMS
*CA ASAS2
*CA CONS2
*CA INDEX
*CA VSCR
      DIMENSION PHIX(3),B(2,2),PHI(2,3)
      DIMENSION PXB(IMXP)
      DATA PHI/0.,0.673,1.35,0.,1.11,0.769/,TAU/1.86E+3/,T00/273./,     
     1SMALL/1.E-15/
      EQUIVALENCE (C(40),B)
C     ****     S13 = HORIZONAL ADVECTION                                
      DO 24 K=1,KMAX
      CALL ADVECL(NX,S13(1,K),K)
   24 CONTINUE                                                          
C     ****     PERIODIC POINTS
      DO 42 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 42 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)
      S13(I,K)=S13(I+IMAX,K)
      S13(I+IMAXP2,K)=S13(I+2,K)
   42 CONTINUE
C     ****     S14=SX/N(X)(K+1/2),  S15=S0(K+1/2)                       
      DO 1 K=1,KMAX
      DO 1 I=1,LEN1
      S15(I,K)=.5*(S2(I,K)+S2(I,K+1))
      S14(I,K)=.5*(S1(I,K)+S1(I,K+1))
    1 CONTINUE                                                          
C     ****     T4 = PS1(-1/2), T5 = PS2(-1/2), T6 = MBAR(-1/2)          
      NPS1K=NJ+NPS                                                      
      NPS2K=NPS1K+KMAXP1                                                
      DO 2 I=1,LEN1                                                     
      T7(I)=T4(I)                                                       
      T4(I)=B(1,1)*F(I,NPS1K)+B(1,2)*F(I,NPS2K)+FB(I,1)
      T5(I)=B(2,1)*F(I,NPS1K)+B(2,2)*F(I,NPS2K)+FB(I,2)
      T6(I)=1./(T4(I)/RMASS(1)+T5(I)/RMASS(2)+(1.-T4(I)-T5(I))/RMASS(3))
    2 CONTINUE                                                          
C     ****     S12 = MBAR(K+1/2)                                        
      NPS1K=NJ+NPS-1
      NPS2K=NPS1K+KMAXP1                                                
      DO 3 K = 1,KMAXP1
      DO 3 I=1,LEN1
      S12(I,K)=1./(F(I,NPS1K+K)/RMASS(1)+F(I,NPS2K+K)/RMASS(2)+(1.-
     1  F(I,NPS1K+K)-F(I,NPS2K+K))/RMASS(3))
    3 CONTINUE                                                          
C     ****     S11 = MBAR(K), S10 = DM/DZ(K), S9 = PS2(K),              
C     ****     S8 = PS1(K), S7 = DPS2/DZ(K), S6 = DPS1/DZ(K)            
C     ****     LOWER BOUNDARY                                           
      DO 4 I=1,LEN1                                                     
      S11(I,1)=.5*(T6(I)+S12(I,1))
      S10(I,1)=(S12(I,1)-T6(I))/C(3)
      S9(I,1)=.5*(T5(I)+F(I,NPS2K+1))
      S8(I,1)=.5*(T4(I)+F(I,NPS1K+1))
      S7(I,1)=(F(I,NPS2K+1)-T5(I))/C(3)
      S6(I,1)=(F(I,NPS1K+1)-T4(I))/C(3)
      T4(I)=T7(I)                                                       
    4 CONTINUE                                                          
C     ****     LEVELS K=2,KMAXP1                                        
      DO 5 K = 2,KMAXP1
      DO 5 I=1,LEN1
      S11(I,K)=.5*(S12(I,K)+S12(I,K-1))
      S10(I,K)=(S12(I,K)-S12(I,K-1))/C(3)
      S9(I,K)=.5*(F(I,NPS2K+K)+F(I,NPS2K+K-1))
      S8(I,K)=.5*(F(I,NPS1K+K)+F(I,NPS1K+K-1))
      S7(I,K)=(F(I,NPS2K+K)-F(I,NPS2K+K-1))/C(3)
      S6(I,K)=(F(I,NPS1K+K)-F(I,NPS1K+K-1))/C(3)
    5 CONTINUE                                                          
C     ****     S5=T(TOT,K)                                              
C     ****     LOWER AND UPPER BOUNDARIES                               
      NTK=NJ+NT+KMAX                                                    
      DO 6 I=1,LEN1                                                     
      S5(I,1)=T0(1)+F(I,NTK)                                            
      S5(I,KMAXP1)=T0(KMAXP1)+F(I,NTK-1)                                
    6 CONTINUE                                                          
C     ****     LEVELS K=2,KMAX                                          
      NTK=NJ+NT                                                         
      DO 7 K=2,KMAX                                                     
      NTK=NTK+1                                                         
      DO 7 I=1,LEN1                                                     
      S5(I,K)=T0(K)+.5*(F(I,NTK)+F(I,NTK-1))                            
    7 CONTINUE                                                          
C     ****     EVALUATE ALFA12, ALFA21, ALFAM1, ALFAM2                  
      ALFA12=PHI(1,2)-PHI(1,3)                                          
      ALFA21=PHI(2,1)-PHI(2,3)                                          
      ALFAX1=PHIX(1)-PHIX(3)                                            
      ALFAX2=PHIX(2)-PHIX(3)                                            
C     ****     S15=(S0*MX/NMBAR)(K+1/2)                                 
      NTK=NJ+NT-1                                                       
      DO 8 K=1,KMAX                                                     
      NTK=NTK+1                                                         
      DO 8 I=1,LEN1                                                     
      S15(I,K)=S15(I,K)*RMX*C(84)*(.5*(T0(K)+T0(K+1))+F(I,NTK))/(C(81)* 
     1EXPS(K)*S12(I,K))                                                 
    8 CONTINUE                                                          
C     ****     IF IBNDB = 1, SET UP FLUX BOUNDARY CONDITION AT BOTTOM   
      IF(IBNDB.EQ.1)THEN                                                
          DO 26 I=1,LEN1                                                
          T1(I)=1.                                                      
          T2(I)=S10(I,1)/S11(I,1)                                       
	  T3(I)=T3(I)*C(54)*RMX/(DIFK(I,1)*C(81)*EXPS(1)*C(87)*C(85))
   26     CONTINUE                                                      
      ENDIF                                                             
C     ****     S1 = ALFA(1,1,K), S2 = ALFA(1,2,K),                      
C     ****     S2 = ALFA(2,1,K), S4 = ALFA(2,2,K)                       
      DO 9 K=1,KMAXP1
      DO 9 I=1,LEN1
      S1(I,K)=-(PHI(1,3)+ALFA12*S9(I,K))
      S2(I,K)=ALFA12*S8(I,K)
      S3(I,K)=ALFA21*S9(I,K)
      S4(I,K)=-(PHI(2,3)+ALFA21*S8(I,K))
C     ****     S12=EX(K)                                                
      S12(I,K)=((ALFAX1*S4(I,K)-ALFAX2*S3(I,K))*(S6(I,K)-(1.-(RMASS(1)+
     1  S10(I,K))/S11(I,K))*S8(I,K))+(ALFAX2*S1(I,K)-ALFAX1*S2(I,K))*
     2  (S7(I,K)-(1.-(RMASS(2)+S10(I,K))/S11(I,K))*S9(I,K)))/(S1(I,K)*
     3  S4(I,K)-S2(I,K)*S3(I,K))+1.-(RMX+S10(I,K))/S11(I,K)
C     ****     S10 = (DM/DZ)/MBAR                                       
      S10(I,K)=S10(I,K)/S11(I,K)
C     ****     S11 = AX(K)                                              
      S11(I,K)=-S11(I,K)/(TAU*RMASS(3))*(T00/S5(I,K))**0.25/(PHIX(3)+
     1  ALFAX1*S8(I,K)+ALFAX2*S9(I,K))
    9 CONTINUE                                                          
C     ****     S12=EX-ALFAX*D/DS(LN(T(TOT))  (THERMAL DIFFUSION TERM)   
      DO 21 K=2,KMAX
      DO 21 I=1,LEN1
      S12(I,K)=S12(I,K)-ALFAX*(S5(I,K+1)-S5(I,K-1))/(2.*C(3)*
     1  S5(I,K))
   21 CONTINUE                                                          
C     ****     LOWER BOUNDARY AND UPPER BOUNDARY                        
      DO 22 I=1,LEN1                                                    
      S12(I,1)=S12(I,1)-ALFAX*(S5(I,2)-S5(I,1))/(C(3)*S5(I,1))
   22 CONTINUE                                                          
      K=KMAXP1
      DO 23 I=1,LEN1
      S12(I,K)=S12(I,K)-ALFAX*(S5(I,K)-S5(I,K-1)
     1  )/(C(3)*S5(I,K))
   23 CONTINUE                                                          
C     ****     T7=DFACT                                                 
      DO 10 I=1,LEN1                                                    
      T7(I)=1.                                                          
   10 CONTINUE                                                          
C     ****     SHAPIRO SMOOTHER,   S10=NXNM  (ELIMINATED)               
      NXNMK=NJ+NXNM-1
      DO 12 K=1,KMAXP1
      DO 12 I=3,LEN1-2
      S9(I,K)=S10(I,K)
       S10(I,K)=F(I,NXNMK+K)-C(26)*(F(I+2,NXNMK+K)+F(I-2,NXNMK+K)-4.*
     1  (F(I+2,NXNMK+K)+F(I-2,NXNMK+K))+6.*F(I,NXNMK+K))
   12 CONTINUE                                                          
C     ****     S1 = P, S2 = Q, S3 = R, S4= F                            
      NWK=NJ+NW-1                                                       
      DO 13 K=1,KMAX                                                    
      NWK=NWK+1                                                         
      DO 13 I=1,LEN1                                                    
      S1(I,K)=S11(I,K)/C(3)*(1./C(3)+.5*S12(I,K))-EXPS(K)*(C(87)*
     1  DIFK(I,K)*T7(I)*(1./C(3)-.5*S9(I,K))+0.25*(F(I,NWK)+
     2  F(I,NWK+1)))/C(3)
      S3(I,K)=S11(I,K+1)/C(3)*(1./C(3)-.5*S12(I,K+1))-EXPS(K)*(C(86)*   
     1  DIFK(I,K+1)*T7(I)*(1./C(3)+.5*S9(I,K+1))-0.25*(F(I,NWK)+
     2  F(I,NWK+1)))/C(3)
      S2(I,K)=-(S11(I,K)/C(3)*(1./C(3)-.5*S12(I,K))+S11(I,K+1)/C(3)*(1./
     1  C(3)+.5*S12(I,K+1)))+EXPS(K)*((C(87)*DIFK(I,K)*(1./C(3)+.5*
     2  S9(I,K))+C(86)*DIFK(I,K+1)*(1./C(3)-.5*S9(I,K+1)))*T7(I)/C(3)-
     3  S14(I,K)+C(7))
      S4(I,K)=EXPS(K)*(S10(I,K)*C(7)-S13(I,K)+S15(I,K))                 
   13 CONTINUE                                                          
C     ****     BOUNDARIES                                               
C     ****     MODIFY EX(KMAXP1) IF IBND=1                              
      IF(IBND.EQ.1)THEN                                                 
           NWK=NJ+NW+KMAX                                               
           DO 25 I=1,LEN1                                               
           T7(I)=EXPS(KMAX)*C(86)*F(I,NWK)/S11(I,KMAXP1)                
           S12(I,KMAXP1)=S12(I,KMAXP1)-T7(I)                            
   25      CONTINUE                                                     
      ENDIF                                                             
      DO 14 I=1,LEN1                                                    
C     ****     LOWER BOUNDARY                                           
      S2(I,1)=S2(I,1)+S1(I,1)*(T1(I)+.5*T2(I)*C(3))/(T1(I)-.5*T2(I)*C(3)
     1)                                                                 
      S4(I,1)=S4(I,1)-S1(I,1)*T3(I)*C(3)/(T1(I)-.5*T2(I)*C(3))          
      S1(I,1)=0.                                                        
C     ****     UPPER BOUNDARY  (UPWARD FLUX GIVEN)                      
      S2(I,KMAX)=S2(I,KMAX)+S3(I,KMAX)*(1.+.5*S12(I,KMAXP1)*C(3))/(1.-  
     1.5*S12(I,KMAXP1)*C(3))                                            
      S4(I,KMAX)=S4(I,KMAX)-S3(I,KMAX)*C(54)*RMX/(C(81)*S11(I,KMAXP1)*  
     1C(85))*T4(I)*C(3)/(1.-.5*S12(I,KMAXP1)*C(3))                      
      S3(I,KMAX)=0.                                                     
   14 CONTINUE                                                          
C     ****     SOLVE TRIDIAGONAL SYSTEM                                 
      NXNPK=NJNP+NX                                                     
      CALL TRSOLV(S1,S2,S3,S4,S7,F(1,NXNPK),LEN1,3,LEN1-2,KMAXP1,1,KMAX
     1,1)                                                               
C     ****     INSERT VALUE OF NXNPK(KMAXP1) USING UPPER BOUNDARY       
C     ****                     CONDITION                                
      NXNPK=NJNP+NX+KMAX                                                
      DO 15 I=3,LEN1-2
      F(I,NXNPK)=(C(54)*RMX/(C(81)*S11(I,KMAXP1)*C(85))*T4(I)*C(3)+     
     1(1.+.5*S12(I,KMAXP1)*C(3))*F(I,NXNPK-1))/(1.-.5*S12(I,KMAXP1)*    
     2C(3))                                                             
   15 CONTINUE                                                          
C     ****     INSERT PERIODIC POINTS
      K1=NJNP+NX
      K2=K1+KMAX
      DO 18 N=1,2
      DO 19 I=1,2
      DO 19 K=K1,K2
      F(I,K)=F(I+IMAX,K)
      F(I+IMAXP2,K)=F(I+2,K)
   19 CONTINUE
      K1=NJNP+NXNM
      K2=K1+KMAX
   18 CONTINUE
C     ****     TIME SMOOTHING
      NXNMK=NJ+NXNM-1
      NXK=NJ+NX-1
      NXNPK=NJNP+NX-1
      NXMNK=NJNP+NXNM-1
      DO 17 K=1,KMAXP1
      DO 17 I=1,LEN1
      F(I,NXMNK+K)=C(30)*F(I,NXK+K)+C(31)*(F(I,NXNMK+K)+F(I,NXNPK+K))
   17 CONTINUE
C     ****     PXB = PSX(-1/2)                                          
      NXNPK=NJNP+NX                                                     
      DO 16 I = 1,LEN1
      PXB(I)=(F(I,NXNPK)*(T1(I)+.5*T2(I)*C(3))+T3(I)*C(3))/
     1(T1(I)-.5*T2(I)*C(3))
   16 CONTINUE
C     ****     INSURE NON-NEGATIVE NX                                   
      NXK=NJNP+NX-1
      NXNMK=NJNP+NXNM-1
      DO 20 K=1,KMAXP1
      DO 20 I=1,LEN1
      F(I,NXK+K)=CVMGP(F(I,NXK+K),SMALL,F(I,NXK+K)-SMALL)
      F(I,NXNMK+K)=CVMGP(F(I,NXNMK+K),SMALL,F(I,NXNMK+K)-SMALL)
   20 CONTINUE                                                          
      RETURN                                                            
      END                                                               
*/
*DK PLOT
*/
      SUBROUTINE PLOT(IFP,N,NJ)
      SAVE
C     ****                                                              
C     ****     PROCEDURE TO PLOT OUT MASS MIXING RATIOS
C     ****                                                              
C     ****     ARRAY IFP CONTAINS INDICES OF FIELDS TO BE PLOTTED       
C     ****                                                              
C     ****     FOR EXAMPLE:                                             
C     ****              NPS    FOR    O2                                
C     ****              NPS2   FOR    O                                 
C     ****              NPSA   FOR    AR                                
C     ****              NPSH   FOR    HE                                
C     ****              NPNO   FOR    NO                                
C     ****              NPN4S  FOR    N4S
C     ****              NPCH4  FOR    CH4                               
C     ****              NPH2   FOR    H2                                
C     ****              NPCO2  FOR    CO2                               
C     ****              NPCO   FOR    CO                                
C     ****              NPHOX  FOR    HOX                               
C     ****              NPH2O  FOR    H2O                               
C     ****                                                              
C     ****     N IS NUMBER OF INDICES PASSED (MAX IS 28)                
C     ****                                                              
C     ****     NJ IS BUFFER POINTER
C     ****                                                              
*CA PARAMS
*CA ASAS2
*CA CONS2
      DIMENSION IFP(1),LABEL(28),LAB2(28)                               
      DATA LABEL/3H O2,3H  O,3H AR,3H HE,3H NO,3HN4S,3HCH4,3H H2,3HCO2, 
     13H CO,3HHOX,3HH2O/                                                
C     ****                                                              
C     ****     SHIFT REQUIRED LABELS TO ARRAY LAB2                      
C     ****                                                              
      DO 1 I=1,N                                                        
       LAB2(I)=LABEL(MOD(IFP(I)/KMAXP1-1,13))
    1 CONTINUE                                                          
C     ****                                                              
C     ****     PLOT FIELDS
C     ****                                                              
      CALL GSCLIP(0)
      DO 2 I = 1,N
	CALL CONREC(F(3,NJ+IFP(I)),IMAXP4,IMAX+1,KMAXP1,0.,0.,0.,0,0,
     1    -1430B)
C       ****
C       ****     WRITE TITLE
C       ****
	CALL GETSET(FXA,FXB,FYA,FYB,XC,XD,YC,YD,LTYPE)
	MXA = KFPX(FXA)
	MXB = KFPX(FXB)
	MYA = KFPY(FYA)
	MYB = KFPY(FYB)
	MX = (MXA+MXB)/2
	MY = MYB+48
	CALL PWRIT(CPUX(MX),CPUY(MY),LAB2(I),3,2,0,0)
	CALL FRAME
    2 CONTINUE                                                          
      CALL GSCLIP(1)
      RETURN
      END                                                               
*/
*DK RATES
*/
      SUBROUTINE RATES(FT)                                              
      SAVE
C     ****                                                              
C     ****     CALCULATES K12, K13, K14, K15 AT ALL LEVELS              
C     ****                                                              
C     ****     COMMON DECKS:
C     ****       PARAMS, CONS, CRATES
C     ****                                                              
*CA PARAMS
*CA CONS2
*CA CRATES2
*CA OXOXOX2
      DIMENSION FT(KMXP,1)
C     ****
C     ****     PERIODIC POINTS FOR OOXR
C     ****
      DO 2 I = 1,2
	DO 2 K = 1,KMAXP1
	  OOXR(K,I) = OOXR(K,I+IMX)
	  OOXR(K,I+IMX+2) = OOXR(K,I+2)
    2 CONTINUE
      DO 1 K=1,KMAXP1                                                   
	DO 1 I=1,LEN1
C       RK13(K,I)=2.15E-34*EXP(345./FT(K,I))
C       RK14(K,I)=2.15E-34*EXP(345./FT(K,I))
C       RK15(K,I)=8.82E-35*EXP(575./FT(K,I))
C NEW   JPL-85 RATES
	RK12(K,I)=OOXR(K,I)*OOXR(K,I)*9.59E-34*EXP(480./FT(K,I))
C       RK12(K,I)=0.
C       RK13(K,I)=OOXR(K,I)*6.0E-34*(300./FT(K,I))**2.8
C       RK14(K,I)=OOXR(K,I)*6.0E-34*(300./FT(K,I))**2.8
C       RK15(K,I)=OOXR(K,I)*6.0E-34*(300./FT(K,I))**2.8
	RK13(K,I)=0.
	RK14(K,I)=0.
	RK15(K,I)=0.
    1 CONTINUE                                                          
      RETURN                                                            
      END                                                               
*/
*DK UTILITY
*/
CDIR$ NOLIST
      SUBROUTINE TRSOLV(A,B,C,F,W,X,IF,I1,I2,KF,K1,K2,NF)               
      SAVE
C     ****     SOLVES TRIDIAGONAL SYSTEM                                
C     ****     A(I,K)*X(I,K-1)+B(I,K)*X(I,K)+C(I,K)*X(I,K+1)=F(I,K)     
C     ****     WHERE I=I1,I2,1 AND K=K1,K2,1 AND A(I,K1)=C(I,K2)=0.     
C     ****     ARRAYS A,B,C,F,X ARE DIMENSIONED (IF,KF)                 
C     ****     WHERE 1.LE.I1.LT.I2.LE.IF AND 1.LE.K1.LT.K2.LE.KF        
C     ****     W IS WORK SPACE DIMENSIONED (IF,3*KF)                    
C     ****     THE ROUTINE USES VECTOR OPERATIONS ON ARRAYS OF LENGTH   
C     ****     IMAX=I2-I1+1                                             
      DIMENSION A(IF,KF),B(IF,KF),C(IF,KF),F(IF,KF),W(IF,3*KF),X(IF,KF) 
      IMAX=I2-I1+1                                                      
      K1P=K1+1                                                          
      GO TO(1,2),NF                                                     
    1 CONTINUE                                                          
C     ****          LOWER BOUNDARY                                      
C     ****          W(K1)=B(K1)                                         
      DO 3 I=I1,I2                                                      
      W(I,K1)=B(I,K1)                                                   
    3 CONTINUE                                                          
C     ****          SWEEP THROUGH K                                     
      DO 4 K=K1P,K2                                                     
C     ****          W(KF+K-1)=C(K-1)/W(K-1)                             
      CALL QVDIV0(W(I1,KF+K-1),C(I1,K-1),W(I1,K-1),IMAX)                
C     ****          W(K)=A(K)*W(KF+K-1)                                 
      CALL QVMPY0(W(I1,K),A(I1,K),W(I1,KF+K-1),IMAX)                    
C     ****          W(K)=B(K)-W(K)                                      
      CALL QVSUB0(W(I1,K),B(I1,K),W(I1,K),IMAX)                         
    4 CONTINUE                                                          
    2 CONTINUE                                                          
C     ****          LOWER BOUNDARY                                      
C     ****          W(2*KF+K1)=F(K1)/W(K1)                              
      CALL QVDIV0(W(I1,2*KF+K1),F(I1,K1),W(I1,K1),IMAX)                 
      DO 5 K=K1P,K2                                                     
C     ****          W(2*KF+K)=A(K)*W(2*KF+K-1)                          
      CALL QVMPY0(W(I1,2*KF+K),A(I1,K),W(I1,2*KF+K-1),IMAX)             
C     ****          W(2*KF+K)=F(K)-W(2*KF+K)                            
      CALL QVSUB0(W(I1,2*KF+K),F(I1,K),W(I1,2*KF+K),IMAX)               
C     ****          W(2*KF+K)=W(2*KF+K)/W(K)                            
      CALL QVDIV0(W(I1,2*KF+K),W(I1,2*KF+K),W(I1,K),IMAX)               
    5 CONTINUE                                                          
C     ****          BACK SUBSTITUTION                                   
C     ****          X(K2)=W(2*KF+K2)                                    
      DO 6 I=I1,I2                                                      
      X(I,K2)=W(I,2*KF+K2)                                              
    6 CONTINUE                                                          
      DO 7 KK=K1P,K2                                                    
      K=K1+K2-KK                                                        
C     ****          X(K)=W(KF+K)*X(K+1)                                 
      CALL QVMPY0(X(I1,K),W(I1,KF+K),X(I1,K+1),IMAX)                    
C     ****          X(K)=W(2*KF+K)-X(K)                                 
      CALL QVSUB0(X(I1,K),W(I1,2*KF+K),X(I1,K),IMAX)                    
    7 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QAA0(W,X,Y,Z,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1),Z(1)                                     
      DO 100 I=1,N                                                      
        W(I) = X(I) + Y(I) + Z(I)                                       
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QAM0(W,X,Y,Z,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1),Z(1)                                     
      DO 100 I=1,N                                                      
        W(I) = ( X(I) + Y(I) ) * Z(I)                                   
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QAM1(W,X,Y,C,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = ( X(I) + Y(I) ) * C                                      
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QAM2(W,X,C,Y,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = ( X(I) + C ) * Y(I)                                      
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QMA0(W,X,Y,Z,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1),Z(1)                                     
      DO 100 I=1,N                                                      
        W(I) = X(I) * Y(I) + Z(I)                                       
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QMA1(W,X,Y,C,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = X(I) * Y(I) + C                                          
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QMA2(W,X,C,Y,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = X(I) * C + Y(I)                                          
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QMA5(W,C,X,S,N)                                        
      SAVE
      DIMENSION W(1),X(1)                                               
      DO 100 I=1,N                                                      
        W(I) = C * X(I) + S                                             
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QMM0(W,X,Y,Z,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1),Z(1)                                     
      DO 100 I=1,N                                                      
        W(I) = X(I) * Y(I) * Z(I)                                       
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QMM1(W,X,Y,C,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = X(I) * Y(I) * C                                          
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QMS0(W,X,Y,Z,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1),Z(1)                                     
      DO 100 I=1,N                                                      
        W(I) = X(I) * Y(I) - Z(I)                                       
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QSA0(W,X,Y,Z,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1),Z(1)                                     
      DO 100 I=1,N                                                      
        W(I) = X(I) - Y(I) + Z(I)                                       
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QSM0(W,X,Y,Z,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1),Z(1)                                     
      DO 100 I=1,N                                                      
        W(I) = ( X(I) - Y(I) ) * Z(I)                                   
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QSM1(W,X,Y,C,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = ( X(I) - Y(I) ) * C                                      
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVADD0(W,X,Y,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = X(I) + Y(I)                                              
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVADD1(W,X,C,N)                                        
      SAVE
      DIMENSION W(1),X(1)                                               
      DO 100 I=1,N                                                      
        W(I) = X(I) + C                                                 
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVDIV0(W,X,Y,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = X(I) / Y(I)                                              
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVDIV2(W,C,X,N)                                        
      SAVE
      DIMENSION W(1),X(1)                                               
      DO 100 I=1,N                                                      
        W(I) = C / X(I)                                                 
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVMPY0(W,X,Y,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = X(I) * Y(I)                                              
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVMPY2(W,C,X,N)                                        
      SAVE
      DIMENSION W(1),X(1)                                               
      DO 100 I=1,N                                                      
        W(I) = C * X(I)                                                 
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVSUB0(W,X,Y,N)                                        
      SAVE
      DIMENSION W(1),X(1),Y(1)                                          
      DO 100 I=1,N                                                      
        W(I) = X(I) - Y(I)                                              
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE QVSUB1(W,X,C,N)                                        
      SAVE
      DIMENSION W(1),X(1)                                               
      DO 100 I=1,N                                                      
        W(I) = X(I) - C                                                 
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               

