Received: from SGCCP.GSFC.NASA.GOV by spica.npac.syr.edu (4.1/I-1.98K)
	id AA02943; Tue, 2 Feb 93 08:56:47 EST
Date:    Tue, 2 Feb 1993 8:57:31 EST
From: ALLEN@SGCCP.GSFC.NASA.GOV
Message-Id: <930202085731.22c00515@SGCCP.GSFC.NASA.GOV>
Subject: transport.subs
To: kim@npac.syr.edu
X-Vmsmail-To: @KIMM



      SUBROUTINE CHEMDO (CONST,PROD,LOSS,TDT)
      IMPLICIT REAL (A-H,O-Z)
C
C     BASIC CHEMISTRY
C
      PARAMETER (N=144,JGMAX=88,NCON=1)
      PARAMETER (kstrat=11,ktrop=8,KTOP=KTROP+KSTRAT)
C
      REAL PROD(0:JGMAX,KTOP,NCON),LOSS(0:JGMAX,KTOP,NCON)
      REAL CONST(N,0:JGMAX,KTOP,NCON)
 
      DO 120 IC=1,NCON
C        PRODUCTION MUST BE GIVEN IN PARTS/SEC.
         DO 80 K=1,KTOP
         DO 80 J=0,JGMAX
         DO 80 I=1,N
            CONST(I,J,K,IC) = CONST(I,J,K,IC) + TDT *
     1       PROD(J,K,IC)
 80      CONTINUE
 
C        COMPLETE CHEMISTRY
         DO 100 K=1,KTOP
         DO 100 J=0,JGMAX
         DO 100 I=1,N
            CONST(I,J,K,IC)=CONST(I,J,K,IC)/(1.+ TDT*LOSS(J,K,IC))
 100     CONTINUE
 120  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE CHEMIN(PROD,LOSS,LUNCHM,ICHEM2,LUNCHM2,THETA)
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (JG2=78,KTOP=19,JGG=90,JGK=JG2*KTOP,JGMAX=88)
 
      REAL DGLAT3(JG2)
 
      REAL PROD(0:JGMAX,KTOP), LOSS(0:JGMAX,KTOP),THETA(0:JGMAX,KTOP)
      REAL AAA(JG2,KTOP)
 
      COMMON // AA2(JG2,KTOP),DLATR(0:JGG  ),
     1 AA3(0:JGG,KTOP),X1(JG2),X2(JGG+1)
 
      DATA DGLAT3/ -88.245, -85.971, -83.684, -81.394, -79.102, -76.810,
     !    -74.518, -72.226, -69.933, -67.640, -65.348, -63.055, -60.762,
     !    -58.469, -56.177, -53.884, -51.591, -49.298, -47.005, -44.712,
     !    -42.419, -40.126, -37.833, -35.541, -33.248, -30.955, -28.662,
     !    -26.369, -24.076, -21.783, -19.490, -17.197, -14.904, -12.611,
     !    -10.318,  -8.025,  -5.732,  -3.439,  -1.146,   1.146,   3.439,
     !      5.732,   8.025,  10.318,  12.611,  14.904,  17.197,  19.490,
     !     21.783,  24.076,  26.369,  28.662,  30.955,  33.248,  35.541,
     !     37.833,  40.126,  42.419,  44.712,  47.005,  49.298,  51.591,
     !     53.884,  56.177,  58.469,  60.762,  63.055,  65.348,  67.640,
     !     69.933,  72.226,  74.518,  76.810,  79.102,  81.394,  83.684,
     !     85.971,  88.245/
 
      DY = 180. / JGG
      DO 10 J=0,JGG
         DLATR(J) = -90.+J*DY
 10   CONTINUE
C     SKIP OVER PRODUCTION AND LOSS
      MRECS = 2
      REWIND(LUNCHM)
      DO 20 ISKIP = 1,ICHEM2-1
      WRITE (6,*) ' SKIPPING OVER CHEMISTRY FILE',ISKIP
      DO 20 IREC = 1, MRECS
         READ (LUNCHM) AAAF
 20   CONTINUE
 
C     READ IN PRODUCTION IN PARTS / SEC IS READ IN.
      READ (LUNCHM) AA2
 
C     INTERPOLATE ONTO TRANSPORT GRID
      CALL INTERX2D(AA2,JG2,KTOP,AA3,JGG+1,DGLAT3,DLATR,X1,X2)
 
      DO 30 K=1,KTOP
         PROD(0,K) = AA3(0,K)
         DO 25 J=1,JGMAX-1
            JP1=J+1
            PROD(J,K) = AA3(JP1,K)
 25      CONTINUE
         PROD(JGMAX,K) = AA3(JGG,K)
 30   CONTINUE
 
      WRITE (6,*) ' CHEMICAL PRODUCTION IN PARTS / SECOND'
      DO 40 J=0,JGMAX
         WRITE (6,2)  (PROD(J,K),K=1,KTOP)
 40   CONTINUE
 
C     LOSS RATES READ IN.
      READ (LUNCHM) AA2
 
C     INTERPOLATE ONTO TRANSPORT GRID
      CALL INTERX2D(AA2,JG2,KTOP,AA3,JGG+1,DGLAT3,DLATR,X1,X2)
 
      DO 50 K=1,KTOP
         LOSS(0,K) = AA3(0,K)
         DO 45 J=1,JGMAX-1
            JP1=J+1
            LOSS(J,K) = AA3(JP1,K)
 45      CONTINUE
         LOSS(JGMAX,K) = AA3(JGG,K)
 50   CONTINUE
 
      WRITE (6,*) ' CHEMICAL LOSS PER SECOND'
      DO 60 J=0,JGMAX
         WRITE (6,2) (LOSS(J,K),K=1,KTOP)
 60   CONTINUE
C***  ORIGINAL O3 CHEMISTRY SUBROUTINE ENDED HERE.****
 
C     SKIP OVER HEADER AND THETA.
      MRECS = 2
      REWIND(LUNCHM2)
      DO 120 ISKIP = 1,ICHEM2-1
      DO 120 IREC = 1, MRECS
         READ (LUNCHM2)
 120  CONTINUE
      WRITE (6,*) ' SKIPPING OVER ', ICHEM2-1, ' THETA OBSERVATIONS'
 
C     READ IN THETA IN INVERSE (SECONDS * KELVIN)
      READ (LUNCHM2) DOBS
      WRITE (6,*) ' THETA COEFFICIENTS FOR MONTH ', DOBS
      READ (LUNCHM2) THETA
 
      WRITE (6,*) ' THETA IN INVERSE (SEC KELVIN)'
      DO 140 J=0,JGMAX
         WRITE (6,2) (THETA(J,K),K=1,KTOP)
 140  CONTINUE
 
      RETURN
 2    FORMAT(10E10.3,/,9E10.3)
      END
C
C
C
      SUBROUTINE CHEMIN2(PROD,LOSS,LUNCHM,ICHEM2)
      IMPLICIT REAL (A-H,O-Z)
C
 
      PARAMETER (JGG=90,KTOP=19,JGK=JGG*KTOP,JGMAX=88,NCON=1)
 
      REAL PROD(0:JGMAX,KTOP,NCON),LOSS(0:JGMAX,KTOP,NCON)
      INTEGER LUNCHM(NCON)
 
      COMMON // AAA(0:JGG,KTOP),AA2(0:JGG,KTOP)
 
      MRECS = 2
      DO 100 IC=1,NCON
         REWIND(LUNCHM(IC))
         DO 20 ISKIP = 1,ICHEM2-1
         DO 20 IREC = 1, MRECS
            WRITE (6,*) ' SKIPPING OVER CHEMISTRY FILES'
            READ (LUNCHM(IC)) AAAF
 20      CONTINUE
 
C        READ IN PRODUCTION IN PARTS / SEC IS READ IN.
c        READ(LUNCHM(IC)) AAA
c        CALL Q9IC32(AAA,AA2,JGK,IERR)
         read (lunchm(ic)) aa2
         DO 30 K=1,KTOP
            PROD(0,K,IC) = AA2(0,K)
            DO 25 J=1,JGMAX-1
               JP1=J+1
               PROD(J,K,IC) = AA2(JP1,K)
 25         CONTINUE
            PROD(JGMAX,K,IC) = AA2(JGG,K)
 30      CONTINUE
 
         WRITE (6,*) ' CHEMICAL PRODUCTION IN PARTS / SECOND'
         DO 40 J=0,JGMAX
            WRITE (6,2)    (PROD(J,K,IC),K=1,KTOP)
 40      CONTINUE
 
C        LOSS RATES READ IN.
c        READ(LUNCHM(IC)) AAA
c        CALL Q9IC32(AAA,AA2,JGK,IERR)
         read (lunchm(ic)) aa2
 
         DO 50 K=1,KTOP
            LOSS(0,K,IC) = AA2(0,K)
            DO 45 J=1,JGMAX-1
               JP1=J+1
               LOSS(J,K,IC) = AA2(JP1,K)
 45         CONTINUE
            LOSS(JGMAX,K,IC) = AA2(JGG,K)
 50      CONTINUE
 
         WRITE (6,*) ' CHEMICAL LOSS PER SECOND'
         DO 60 J=0,JGMAX
            WRITE (6,2) (LOSS(J,K,IC),K=1,KTOP)
 60      CONTINUE
 100  CONTINUE
 
      RETURN
 2    FORMAT(10E10.3,/,9E10.3)
      END
C
C
C
      SUBROUTINE AJRATE(LJRAT,AJ,PROD)
      IMPLICIT REAL (A-H,O-Z)
 
C     ZONALLY AVERAGED JRATES FOR N2O PHOTOLYSIS ARE READ IN FROM
 
 
      PARAMETER (JMP=18,ktrop=8,kstrat=11,JGG=90,JG1=JGG+1)
      PARAMETER (KTOP=KTROP+KSTRAT,JGMAX=88)
      PARAMETER (KJG1=(KTOP-1)*(JGMAX+1))
 
      REAL AJ(0:JGMAX,KTOP), PROD(0:JGMAX,KTOP)
 
      COMMON // DLAT(JMP), AJ1(JMP,KTOP), P1(KTOP),DGLAT(0:JGG),
     1  AJ2(0:JGG,KTOP),X1(JMP),X2(JG1)
 
      COMMON / VERT/ DELTAS(KTOP), DSIG(0:KTOP), DELZ(KSTRAT),
     2 SIG(KTOP),P(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)
 
      COMMON / THREE/ HTSCALE, HINV, PI
 
C     STEP1. READ IN PHOTOLYSIS RATES (NOTE: FLIPPING NOT NEC JAN)
      WRITE (6,*) ' INPUT N2O PHOTOLYSIS RATES'
 
      DO 10 K=1,KTOP
*        K2 = KTOP-K+1
         READ(LJRAT,11) P1(K), (AJ1(J,K),J=1,JMP)
         WRITE (6  ,11) P1(K), (AJ1(J,K),J=1,JMP)
 10   CONTINUE
 11   FORMAT(F7.3,(9E12.4))
C
C     STEP4: N-S INTERPOLATION ONTO TRANSPORT LATITUDES.
      DY = 180. / JGG
      DO 15 IJ=0,JGG
         DGLAT(IJ) = -90. +IJ*DY
 15   CONTINUE
      DO 20 IJ=1,JMP
         DLAT(IJ) = -85. + (IJ-1)*10.
 20   CONTINUE
 
      CALL INTERX2D(AJ1,JMP,KTOP,AJ2,JG1,DLAT,DGLAT,X1,X2)
 
      DO 45 K=1,KTOP
         AJ(0,K) = AJ2(0,K)
         DO 42 J=1,JGMAX-1
            AJ(J,K) = AJ2(J+1,K)
 42      CONTINUE
         AJ(JGMAX,K) = AJ2(JGG,K)
 45   CONTINUE
 
C     WRITE (6,*) ' LOSS COEFFICIENTS AFTER INTERPOLATION'
      DO 60 K=1,KTOP
         WRITE (6,12) K,(AJ(J,K),J=0,JGMAX)
 60   CONTINUE
 12   FORMAT(I5,(10E11.3))
C
C     PRODUCTION COEFFICIENTS ARE SET FOR N2O.
      DO 80 J=0,JGMAX
         PROD(J,1) = 7.17E-16
 80   CONTINUE
      DO 100 K=2,KTOP
      DO 100 J=0,JGMAX
         PROD(J,K) = 0.
 100  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE HCLCHEM(CONST,KEL1,PROD,LOSS,CL,CLRAT,TDT)
      IMPLICIT REAL (A-H,O-Z)
C
C      CHEMISTRY FOR HCL-DHCL-DHCL1 RUN OF 3-D MODEL
C      1) HCL GAS PHASE CHEMISTRY ONLY
C      2) HCL WITH ADDED HETEROGENEOUS LOSS
C      3) HCL WITH ADDED HETEROGENEOUS LOSS AND DEPLETION-DEPENDENT PROD
C
C
      PARAMETER (N=144,JGMAX=88,NCON=1)
      PARAMETER (kstrat=11,ktrop=8,KTOP=KTROP+KSTRAT)
C
      REAL PROD(0:JGMAX,KTOP,NCON),LOSS(0:JGMAX,KTOP,NCON)
      REAL CONST(N,0:JGMAX,KTOP,NCON), KEL1(N,0:JGMAX,KTOP)
      REAL CL(0:JGMAX,KTOP),CLRAT(0:JGMAX,KTOP)
 
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),
     2 SIGMA(KTOP),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),
     3 SGINT(0:KTOP)
 
      REAL PROD2, LOSS2
 
      COMMON // PROD2(N,0:JGMAX,KTOP), LOSS2(N,0:JGMAX,KTOP)
 
C     PRODUCTION IN MIXING RATIO PER SECOND
 
********************SPECIAL SKIP********************
C     DO 100 IC=NCON,1,-1
C        IF ((IC.EQ.1).OR.(IC.EQ.2)) THEN
C           DO 20 K=1,KTOP
C           DO 20 J=0,JGMAX
C              DO 10 IL=1,N
C                 PROD2(IL,J,K)  = PROD(J,K,IC) * CL(J,K)
C                 LOSS2(IL,J,K)  = LOSS(J,K,IC)
C10            CONTINUE
C20         CONTINUE
C        ELSE
C           DO 40 K=1,KTOP
C           DO 40 J=0,JGMAX
C           DO 40 IL=1,N
C              PROD2(IL,J,K) = PROD(J,K,IC)*(CL(J,K)+CLRAT(J,K)
C    1                      *(CONST(IL,J,K,1)-CONST(IL,J,K,3)))
C              LOSS2(IL,J,K)  = LOSS(J,K,IC)
C40         CONTINUE
C        END IF
C
C        DO 80 K=1,KTOP
C        DO 80 IJ=0,JGMAX
C        DO 80 IL=1,N
C           CONST(IL,IJ,K,IC)=(CONST(IL,IJ,K,IC)+
C    1       TDT*PROD2(IL,IJ,K))/(1.+TDT*LOSS2(IL,IJ,K))
C80      CONTINUE
C100  CONTINUE
*************END SPECIAL SKIP*****************
 
C     HETEROGENEOUS LOSS DHCL AND DHCL1 FOR COLD T'S NORTH OF 50N
C     HARDWIRED FOR 2 DEGREE N-S RESOLUTION
C     LATITUDE J=69 CORRESPONDS TO 50N
      DO 160 IC=2,NCON
         JLAT=69
         JN=(JGMAX-J+1)*N
         DO 110 IJ=JLAT,JGMAX
         DO 110 IL=1,N
            IF (KEL1(IL,IJ,8).LT.199.) CONST(IL,IJ,8,IC) =
     1        CONST(IL,IJ,8,IC)  * EXP(-4.E-6*TDT)
 
            IF (KEL1(IL,IJ,9).LT.197.) CONST(IL,IJ,9,IC) =
     1         CONST(IL,IJ,9,IC) * EXP(-4.E-6*TDT)
 
            IF (KEL1(IL,IJ,10).LT.195.) CONST(IL,IJ,10,IC)=
     1         CONST(IL,IJ,10,IC) * EXP(-4.E-6*TDT)
 
            IF (KEL1(IL,IJ,11).LT.192.) CONST(IL,IJ,11,IC)=
     1          CONST(IL,IJ,11,IC) * EXP(-4.E-6*TDT)
 110     CONTINUE
 160  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE HCLIN(PROD,LOSS,CL,CLRAT,LUNCHM,ICHEM2)
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (JGMAX=88,KTOP=19,JGG=91,JGK=JGG*KTOP,NCON=1)
      PARAMETER (KJG=(JGMAX+1)*KTOP)
 
      INTEGER LUNCHM(NCON)
 
      REAL PROD(0:JGMAX,KTOP,NCON), LOSS(0:JGMAX,KTOP,NCON),
     1 CL(0:JGMAX,KTOP),CLRAT(0:JGMAX,KTOP)
 
      REAL DIOBS, AAA(JGG,KTOP)
 
      COMMON //  BBB(JGG,KTOP),
     1 ALF(0:JGMAX,KTOP),BET(0:JGMAX,KTOP),GAM(0:JGMAX,KTOP)
 
C     HDR,PROD,LOSS,CHLORINE,ALF,BET, AND GAM ARE SKIPPED OVER
      MRECS = 7
      REWIND(LUNCHM(1))
 
      DO 20 ISKIP = 1,ICHEM2
      DO 20 IREC = 1, MRECS
         READ (LUNCHM(1)) AAAF
 20   CONTINUE
      WRITE (6,*) ICHEM2, ' CHEMISTRY FILES SKIPPED OVER '
 
      READ (LUNCHM(1)) DIOBS2
      WRITE (6,22) DIOBS2
 22   FORMAT (' CHEMICAL COEFFICIENT SET ',F4.0, ' ACCESSED')
 
C     PRODUCTION COEFFICIENTS READ IN
      READ (LUNCHM(1)) BBB
      DO 30 K=1,KTOP
         PROD(0,K,1 ) = BBB(1,K)
         DO 25 J=1,JGMAX-1
            JP2=J+2
            PROD(J,K,1 ) = BBB(JP2,K)
 25      CONTINUE
         PROD(JGMAX,K,1 ) = BBB(JGG,K)
 30   CONTINUE
 
        DO 32 K=1,KTOP
       DO 32 IJ=0,JGMAX
*************WARNING 79 EXPERIMENT MULTIPLY BY 2. / 2.8
          PROD(IJ,K,1) = PROD(IJ,K,1) * 2. / 2.8
*******************************************************
C         PROD(IJ,K,2) = PROD(IJ,K,1)
C         PROD(IJ,K,3) = PROD(IJ,K,1)
  32   CONTINUE
 
      WRITE (6,*) ' PRODUCTION / CHLORINE '
      DO 33 J=0,JGMAX
         WRITE (6,2) (PROD(J,K,1),K=1,KTOP)
 33   CONTINUE
 
C     LOSS COEFFICIENTS READ IN
      READ (LUNCHM(1)) BBB
      DO 40 K=1,KTOP
         LOSS(0,K,1 ) = BBB(1,K)
         DO 35 J=1,JGMAX-1
            JP2=J+2
            LOSS(J,K,1 ) = BBB(JP2,K)
 35      CONTINUE
         LOSS(JGMAX,K,1 ) = BBB(JGG,K)
 40   CONTINUE
      WRITE (6,*) ' LOSS COEFFS '
      DO 42 J=0,JGMAX
         WRITE (6,2) (LOSS(J,K,1),K=1,KTOP)
 42   CONTINUE
C      DO 43 K=1,KTOP
C      DO 43 IJ=0,JGMAX
C         LOSS(IJ,K,2) = LOSS(IJ,K,1)
C         LOSS(IJ,K,3) = LOSS(IJ,K,1)
C 43   CONTINUE
 
 
C     CHLORINE READ IN.
      read (lunchm(1)) bbb
      DO 50 K=1,KTOP
         CL(0,K) = BBB(1,K)
         DO 45 J=1,JGMAX-1
            JP2=J+2
            CL(J,K) = BBB(JP2,K)
 45      CONTINUE
         CL(JGMAX,K) = BBB(JGG,K)
 50   CONTINUE
      WRITE (6,*) ' CHLORINE '
      DO 52 J=0,JGMAX
         WRITE (6,2) (CL(J,K),K=1,KTOP)
 52   CONTINUE
 
C     ALFA IS READ IN
      READ (LUNCHM(1)) BBB
      DO 60 K=1,KTOP
         ALF(0,K) = BBB(1,K)
         DO 55 J=1,JGMAX-1
            JP2=J+2
            ALF(J,K) = BBB(JP2,K)
 55      CONTINUE
         ALF(JGMAX,K) = BBB(JGG,K)
 60   CONTINUE
      WRITE (6,*) ' ALF '
      DO 62 J=0,JGMAX
         WRITE (6,2) (ALF(J,K),K=1,KTOP)
 62   CONTINUE
 
C     BETA IS READ IN
      read (lunchm(1)) bbb
      DO 70 K=1,KTOP
         BET(0,K) = BBB(1,K)
         DO 65 J=1,JGMAX-1
            JP2=J+2
            BET(J,K) = BBB(JP2,K)
 65      CONTINUE
         BET(JGMAX,K) = BBB(JGG,K)
 70   CONTINUE
      WRITE (6,*) 'BET  '
      DO 72 J=0,JGMAX
         WRITE (6,2) (BET(J,K),K=1,KTOP)
 72   CONTINUE
 
C     GAMMA IS READ IN
      read (lunchm(1)) bbb
      DO 80 K=1,KTOP
         GAM(0,K) = BBB(1,K)
         DO 75 J=1,JGMAX-1
            JP2=J+2
            GAM(J,K) = BBB(JP2,K)
 75      CONTINUE
         GAM(JGMAX,K) = BBB(JGG,K)
 80   CONTINUE
      WRITE (6,*) ' GAM '
      DO 82 J=0,JGMAX
         WRITE (6,2) (GAM(J,K),K=1,KTOP)
 82   CONTINUE
 
      DO 85 K=1,KTOP
      DO 85 IJ=0,JGMAX
         CLRAT(IJ,K) = ALF(IJ,K) / (1.+ALF(IJ,K)+BET(IJ,K)+
     1    GAM(IJ,K))
 85   CONTINUE
 
      RETURN
 2    FORMAT(10E10.3,/,9E10.3)
      END
C
C
C
      SUBROUTINE N2OIN(PROD,LOSS,LUNCHM,ICHEM2)
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (JGMAX=88,KTOP=19,JGG=91,JGK=JGG*KTOP)
 
      REAL  PROD(0:JGMAX,KTOP),LOSS(0:JGMAX,KTOP)
      COMMON //  BBB(JGG,KTOP)
 
C     N2O PRODUCTION IS SPECIFIED.
      DO 5 J=0,JGMAX
         PROD(J,1) = 7.17E-16
 5    CONTINUE
      DO 10 K=2,KTOP
      DO 10 J=0,JGMAX
         PROD(J,K) = 0.
 10   CONTINUE
 
C     HDR AND LOSS SKIPPED OVER.
      MRECS = 2
      REWIND(LUNCHM)
      DO 20 ISKIP = 1,ICHEM2
      DO 20 IREC = 1, MRECS
         READ (LUNCHM) AAAF
 20   CONTINUE
      WRITE (6,*) ICHEM2, ' CHEMISTRY FILES SKIPPED OVER '
      READ (LUNCHM) DIOBS
      WRITE (6,*) 'CHEMISTRY FILE ',DIOBS, ' ACCESSED '
 
C     LOSS COEFFICIENTS READ IN
      READ (LUNCHM) BBB
      DO 40 K=1,KTOP
         LOSS(0,K) = BBB(1,K)
         DO 35 J=1,JGMAX-1
            JP2=J+2
            LOSS(J,K ) = BBB(JP2,K)
 35      CONTINUE
         LOSS(JGMAX,K ) = BBB(JGG,K)
 40   CONTINUE
      WRITE (6,*) ' LOSS COEFFS '
      DO 42 J=0,JGMAX
         WRITE (6,2) (LOSS(J,K),K=1,KTOP)
 42   CONTINUE
 
      RETURN
 2    FORMAT(10E10.3,/,9E10.3)
      END
C
C
C
      SUBROUTINE NOYIN(PROD,LOSS,LUNCHM,ICHEM2)
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (JGMAX=88,KTOP=19,JGG=91,JGK=JGG*KTOP)
 
      REAL  PROD(0:JGMAX,KTOP),LOSS(0:JGMAX,KTOP)
      COMMON //  BBB(JGG,KTOP)
 
 
C     HDR, PROD,  AND LOSS SKIPPED OVER.
      MRECS = 3
      REWIND(LUNCHM)
      DO 20 ISKIP = 1,ICHEM2
      DO 20 IREC = 1, MRECS
         READ (LUNCHM) AAAF
 20   CONTINUE
      WRITE (6,*) ICHEM2, ' CHEMISTRY FILES SKIPPED OVER '
      READ (LUNCHM) DIOBS
      WRITE (6,*) 'CHEMISTRY FILE ',DIOBS, ' ACCESSED '
 
C     PROD COEFFICIENTS READ IN
      READ (LUNCHM) BBB
      DO 40 K=1,KTOP
         PROD(0,K) = BBB(1,K)
         DO 35 J=1,JGMAX-1
            JP2=J+2
            PROD(J,K ) = BBB(JP2,K)
 35      CONTINUE
         PROD(JGMAX,K ) = BBB(JGG,K)
 40   CONTINUE
      WRITE (6,*) ' PROD COEFFS '
      DO 42 J=0,JGMAX
         WRITE (6,2) (PROD(J,K),K=1,KTOP)
 42   CONTINUE
 
C     LOSS COEFFICIENTS READ IN
      READ (LUNCHM) BBB
      DO 60 K=1,KTOP
         LOSS(0,K) = BBB(1,K)
         DO 55 J=1,JGMAX-1
            JP2=J+2
            LOSS(J,K ) = BBB(JP2,K)
 55      CONTINUE
         LOSS(JGMAX,K ) = BBB(JGG,K)
 60   CONTINUE
      WRITE (6,*) ' LOSS COEFFS '
      DO 70 J=0,JGMAX
         WRITE (6,2) (LOSS(J,K),K=1,KTOP)
 70   CONTINUE
 
      RETURN
 2    FORMAT(10E10.3,/,9E10.3)
      END
C
C
C
      SUBROUTINE CHEMO3 (CONST,PROD,LOSS,TDT,CHEMPROD,CHEMLOSS,NDIAG,
     1 DIAG,KEL,THETA,DIURNAL,IDARK,TIME)
 
      IMPLICIT REAL (A-H,O-Z)
C
C     UNITS OF PRODUCTION ARE PARTS/SEC.
C
      PARAMETER (N=144,JGMAX=88,NCON=1)
      PARAMETER (kstrat=11,ktrop=8,KTOP=KTROP+KSTRAT)
      PARAMETER (JDIAG=5,JDIAG2=3,KDIAG2=3,NVARS=10)
C
      REAL PROD(0:JGMAX,KTOP,NCON),LOSS(0:JGMAX,KTOP,NCON)
      REAL CONST(N,0:JGMAX,KTOP,NCON),F1(N,0:JGMAX,KTOP)
      REAL KEL(N,0:JGMAX,KTOP),THETA(0:JGMAX,KTOP),DARK(N,0:JGMAX)
      REAL CHEMPROD(N,JDIAG,KTOP,NCON), CHEMLOSS(N,JDIAG,KTOP,NCON)
      INTEGER NDIAG(JDIAG)
      LOGICAL DIAG,IDARK(N,0:JGMAX), DIURNAL
 
      COMMON / DIAGNOSTIC / DD(N,JDIAG2,KDIAG2,NVARS), NNLAT(JDIAG2),
     1 NNLEV(KDIAG2)
 
      real KELMEAN(0:JGMAX,KTOP),PROD2(N,0:JGMAX,KTOP)

      COMMON / VERT/ DELTAS(KTOP), DSIG(0:KTOP), DELZ(KSTRAT),
     2 SIG(KTOP),P(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)
 
      IF (DIURNAL) THEN
         DO 10 J=0,JGMAX
         DO 10 I=1,N
            IF (IDARK(I,J)) THEN
               DARK(I,J) = 0.
            ELSE
               DARK(I,J) = 1.
            END IF
 10      CONTINUE
      END IF
 
C     ZONALLY AVERAGED TEMPERATURE IS DETERMINED.
      DN = 1./N
      DO 20 K=1,KTOP
      DO 20 J=0,JGMAX
         KELMEAN(J,K) = SSUM(N,KEL(1,J,K),1) * DN
 20   CONTINUE
 
      DO 140 IC=1,NCON
         IF (DIAG) THEN
            DO 30 K=1,KTOP
            DO 30 J=0,JGMAX
            DO 30 I=1,N
               F1(I,J,K) = CONST(I,J,K,IC)
 30         CONTINUE
         END IF
 
         IF (DIURNAL) THEN
            DO 40 K=1,KTOP
            DO 40 J=0,JGMAX
            DO 40 I=1,N
               CONST(I,J,K,IC) = CONST(I,J,K,IC) + TDT *
     1          (PROD(J,K,IC)*DARK(I,J) - THETA(J,K)*DARK(I,J)
     1       *  (KEL(I,J,K)-KELMEAN(J,K)))
 40         CONTINUE
         ELSE
            do 50 k=1,ktrop
            do 50 j=0,jgmax
            do 50 i=1,n
               prod2(i,j,k) = prod(j,k,ic)
 50         continue 
                     
C           PROD IS NOT ADJUSTED WHERE CHEM NOT DOMINANT.
             DO 90 K=ktrop+1,ktop
             do 90 j=0,jgmax
C             1 < CL/P < 5
C             THE LOSS EXCEEDS 2.31E-6 WHERE CHEMISTRY IS IMPORTANT.
               if ((p(k-ktrop).lt.0.65).or.(p(k-ktrop).gt.31.0).or.
     1          (LOSS(J,K,IC).lT.2.31E-6)) THEN
                  DO 70 I=1,N
                     PROD2(I,J,K) = PROD(J,K,IC)
 70               CONTINUE
               else 
                  DO 80 I=1,N
                     T2 = AMAX1(AMIN1(5.,CONST(I,J,K,IC)*(LOSS(J,K,IC)/
     1                PROD(J,K,IC))),1.)
                     PROD2(I,J,K) = PROD(J,K,IC) * (1.+(T2-1) *
     1                EXP(-1./(LOSS(J,K,IC)*345600.)))
 80               CONTINUE
               end if
 90        CONTINUE
 
            IF (TIME.LT.2000.) THEN
               WRITE (6,*) ' PROD2 COEFS AT 0 EAST '
               DO 95 J=0,JGMAX
                  WRITE (6,2) (PROD2(1,J,K),K=1,KTOP)
 95            CONTINUE
            END IF
 
            DO 100 K=1,KTOP
            DO 100 J=0,JGMAX
            DO 100 I=1,N
               CONST(I,J,K,IC) = CONST(I,J,K,IC) + TDT *
     1          (PROD2(I,J,K) - THETA(J,K)
     1       *  (KEL(I,J,K)-KELMEAN(J,K)))
 100        CONTINUE
         END IF
 
         IF (DIAG) THEN
            CALL DIAG3(F1,CONST(1,0,1,IC),CHEMPROD(1,1,1,IC),
     1       CHEMLOSS(1,1,1,IC),NDIAG,1)
            DO 110 K=1,KTOP
            DO 110 J=0,JGMAX
            DO 110 I=1,N
               F1(I,J,K) = CONST(I,J,K,IC)
 110         CONTINUE
         END IF
 
C        COMPLETE CHEMISTRY
         IF (DIURNAL) THEN
            DO 120 K=1,KTOP
            DO 120 J=0,JGMAX
            DO 120 I=1,N
               CONST(I,J,K,IC)=CONST(I,J,K,IC)/(1.+ TDT*LOSS(J,K,IC)
     1          * DARK(I,J))
 120        CONTINUE
         ELSE
            DO 130 K=1,KTOP
            DO 130 J=0,JGMAX
            DO 130 I=1,N
               CONST(I,J,K,IC)=CONST(I,J,K,IC)/(1.+ TDT*LOSS(J,K,IC))
 130        CONTINUE
         END IF
 
         IF (DIAG) CALL DIAG3(F1,CONST(1,0,1,IC),CHEMPROD(1,1,1,IC),
     1    CHEMLOSS(1,1,1,IC),NDIAG,2)
 140  CONTINUE
 
 2    FORMAT(10E10.3,/,9E10.3)
      RETURN
      END
C
C
C
      SUBROUTINE O3IN(PROD1,PROD2,LOSS1,LOSS2,THETA1,THETA2,ICHEM2,
     1                LUNCHM,LUNCHM2,DIURNAL)
 
      IMPLICIT REAL (A-H,O-Z)
 
C     CHEMISTRY COEFFICIENTS ARE READ IN FOR E368.
 
      PARAMETER (KTOP=19,JGMAX=88,N=144)
 
      REAL PROD1(0:JGMAX,KTOP),LOSS1(0:JGMAX,KTOP),THETA1(0:JGMAX,KTOP)
      REAL PROD2(0:JGMAX,KTOP),LOSS2(0:JGMAX,KTOP),THETA2(0:JGMAX,KTOP)
      REAL FRAC(N,0:JGMAX)
      LOGICAL DIURNAL
 
      IF (DIURNAL) THEN
         REWIND(48)
         MRECS = 1
         DO 15 ISKIP=1,ICHEM2
         DO 15 IREC=1,MRECS
            READ (48)
 15      CONTINUE
         READ (48) FRAC
         WRITE (6,*) ' FRACTIONAL DAYLIGHT AS A F(LAT) '
         WRITE (6,99) (FRAC(1,J),J=0,JGMAX)
 99      FORMAT(16F7.2)
      END IF
 
      REWIND(LUNCHM)
C     SKIP OVER PRODUCTION AND LOSS
      MRECS = 3
      DO 20 ISKIP = 1,ICHEM2
      DO 20 IREC = 1, MRECS
         READ (LUNCHM)
 20   CONTINUE
      WRITE (6,*) ICHEM2, ' CHEMISTRY OBSERVATIONS SKIPPED OVER '
C     READ IN PRODUCTION IN PARTS / SEC IS READ IN.
      READ (LUNCHM) AAA
      WRITE (6,*) ' P + L FOR SEMI-MONTH ', INT(AAA+0.01),' READ IN '
      READ (LUNCHM) PROD1
      READ (LUNCHM) LOSS1
      READ (LUNCHM) AAA1
      WRITE (6,*) ' P + L FOR SEMI-MONTH ', INT(AAA1+0.01),' READ IN '
      READ (LUNCHM) PROD2
      READ (LUNCHM) LOSS2
 
      IF (DIURNAL) THEN
         DO 35 K=1,KTOP
         DO 35 J=0,JGMAX
            IF (FRAC(1,J).GT.0.01) THEN
               PROD1(J,K)  = PROD1(J,K) / FRAC(1,J)
               LOSS1(J,K)  = LOSS1(J,K) / FRAC(1,J)
               PROD2(J,K)  = PROD2(J,K) / FRAC(1,J)
               LOSS2(J,K)  = LOSS2(J,K) / FRAC(1,J)
            END IF
 35      CONTINUE
      END IF
 
      WRITE (6,*) ' CHEMICAL PRODUCTION IN PARTS / SECOND'
      DO 40 J=0,JGMAX
         WRITE (6,2)  (PROD1(J,K),K=1,KTOP)
 40   CONTINUE
 
      WRITE (6,*) ' CHEMICAL LOSS PER SECOND'
      WRITE (6,*) ' DEPOSITION LOSS IS INCLUDED (LAYER 1) '
C     DEPOSITION VELOCITY OF 0.1 CM/S IS ASSUMED
      DO 50 J=0,JGMAX
         LOSS1(J,1) = LOSS1(J,1) + 1.2E-6
         LOSS2(J,1) = LOSS2(J,1) + 1.2E-6
 50   CONTINUE
C
      DO 60 J=0,JGMAX
         WRITE (6,2) (LOSS1(J,K),K=1,KTOP)
 60   CONTINUE
 
C     SKIP OVER HEADER AND THETA.
      MRECS = 2
      REWIND(LUNCHM2)
      DO 120 ISKIP = 1,ICHEM2/2
      DO 120 IREC = 1, MRECS
         READ (LUNCHM2)
 120  CONTINUE
      WRITE (6,*) ICHEM2/2, ' THETA OBS  SKIPPED OVER '
 
C     READ IN THETA IN INVERSE (SECONDS * KELVIN)
      READ (LUNCHM2) DOBS
      WRITE (6,*) ' THETA COEFFICIENTS FOR MONTH ', INT(DOBS+0.01)
      READ (LUNCHM2) THETA1
C     READ IN THETA IN INVERSE (SECONDS * KELVIN)
      READ (LUNCHM2) DOBS1
      WRITE (6,*) ' THETA COEFFICIENTS FOR MONTH ', INT(DOBS1+0.01)
      READ (LUNCHM2) THETA2
      IF (DIURNAL) THEN
         DO 135 K=1,KTOP
         DO 135 J=0,JGMAX
            IF (FRAC(1,J).GT.0.01) THETA1(J,K)=THETA1(J,K)/FRAC(1,J)
            IF (FRAC(1,J).GT.0.01) THETA2(J,K)=THETA2(J,K)/FRAC(1,J)
 135     CONTINUE
      END IF
 
      WRITE (6,*) ' THETA IN INVERSE (SEC KELVIN)'
      DO 140 J=0,JGMAX
         WRITE (6,2) (THETA1(J,K),K=1,KTOP)
 140  CONTINUE
 
      RETURN
 2    FORMAT(10E10.3,/,9E10.3)
      END
C
C
C
      SUBROUTINE JULDAY(JYR,JMNTH,JDAY,JHR,JMN,JD,FJD)
C-----------------------------------------------------------------------
C
C     COMPUTES JULIAN DAY AND FRACTION FROM YEAR, MONTH, DAY AND TIME UT
C     ACCURATE ONLY BETWEEN MARCH 1, 1900 AND FEBRUARY 28, 2100
C     BASED ON JULIAN CALENDAR CORRECTED TO CORRESPOND TO GREGORIAN
C        CALENDAR DURING THIS PERIOD
C
CJEN--------------------------------------------------------------------
      SAVE NDM,JDOR,JYR19
CJEN--------------------------------------------------------------------
      DIMENSION NDM(12)
C
C    JDOR=JD OF DECEMBER 30, 1899 AT 12 HOURS UT
      DATA JDOR/2415019/,JYR19/1900/
      DATA NDM/0,31,59,90,120,151,181,212,243,273,304,334/
C
      JD = JDOR
      JYRM9 = JYR - JYR19
      LP = JYRM9/4
      IF (LP.GT.0) THEN
         JD = JD + 1461*LP
      ENDIF
      NY = JYRM9 - 4*LP
      IC = 0
      IF (NY.LE.0) THEN
         IF (JMNTH.GT.2) IC = 1
      ELSE
         JD = JD + 365*NY + 1
      ENDIF
      JD = JD + NDM(JMNTH) + JDAY + IC
      IF (JHR.LT.12) THEN
         JD = JD - 1
         FJD = .5 + .041666667*FLOAT(JHR) + .00069444444*FLOAT(JMN)
      ELSE
         FJD = .041666667*FLOAT(JHR-12) + .00069444444*FLOAT(JMN)
      ENDIF
      RETURN
      END
C
C
C
      SUBROUTINE RADSLR(JD,FJD,R,DLT,ALP,SLAG,ALAT)
C-----------------------------------------------------------------------
C
CL            COMPUTE RADIUS VECTOR, DECLINATION AND RIGHT ASCENSION OF
CL            SUN, EQUATION OF TIME, AND HOUR ANGLE OF SUN AT SUNSET FOR
CL            EQUALLY SPACED LATITUDES GIVEN JULIAN DAY AND FRACTION.
C
C-----------------------------------------------------------------------
      IMPLICIT REAL(A-H)
      IMPLICIT REAL(O-Z)
CJEN--------------------------------------------------------------------
      SAVE PI,TPI,HPI,RAD,CYEAR,CCR,TPP,SVT6,JDOR
CJEN--------------------------------------------------------------------
      PARAMETER(JGMAX=88)
      DIMENSION HANG(0:JGMAX),COSZ(0:JGMAX),TAUDA(0:JGMAX)
      DIMENSION ALAT(0:JGMAX)
      DATA PI/3.1415927/,TPI/6.2831853/,HPI/1.5707963/,RAD/57.29578/
      DATA CYEAR/365.25/,CCR/1.E-7/
C     TPP=DAYS BETWEEN EPOCH AND PERIHELION PASSAGE OF 1900
C     SVT6=DAYS BETWEEN PERIHELION PASSAGE AND MARCH EQUINOX OF 1900
      DATA TPP/1.55/,SVT6/78.035/
      DATA JDOR/2415020/
C
C     JDOR=JD OF EPOCH WHICH IS JANUARY 0, 1900 AT 12 HOURS UT
C
      DAT=FLOAT(JD-JDOR)-TPP+FJD
C    COMPUTES TIME IN JULIAN CENTURIES AFTER EPOCH
      T=FLOAT(JD-JDOR)/36525.
C    COMPUTES LENGTH OF ANOMALISTIC AND TROPICAL YEARS (MINUS 365 DAYS)
      YEAR=.25964134+.304E-5*T
      TYEAR=.24219879-.614E-5*T
C    COMPUTES ORBIT ECCENTRICITY AND ANGLE OF EARTH"S INCLINATION FROM T
      EC=.01675104-(.418E-4+.126E-6*T)*T
      ANGIN=23.452294-(.0130125+.164E-5*T)*T
      ADOR=JDOR
      JDOE=ADOR+(SVT6*CYEAR)/(YEAR-TYEAR)
C    DELEQN=UPDATED SVT6 FOR CURRENT DATE
      DELEQN=FLOAT(JDOE-JD)*(YEAR-TYEAR)/CYEAR
      YEAR=YEAR+365.
      SNI=SIN(ANGIN/RAD)
      TINI=1./TAN(ANGIN/RAD)
      ER=SQRT((1.+EC)/(1.-EC))
      QQ=DELEQN*TPI/YEAR
C    DETERMINE TRUE ANOMALY AT EQUINOX
      E=1.0
 32   EP=E-(E-EC*SIN(E)-QQ)/(1.-EC*COS(E))
      CD=ABS(E-EP)
      E=EP
      IF(CD.GT.CCR) GO TO 32
      HE=.5*E
      EQ=2.*ATAN(ER*TAN(HE))
C    DATE=DAYS SINCE LAST PERIHELION PASSAGE
      DATE=AMOD(DAT,YEAR)
C    SOLVE ORBIT EQUATIONS BY NEWTON"S METHOD
      EM=TPI*DATE/YEAR
      E=1.0
 31   EP=E-(E-EC*SIN(E)-EM)/(1.-EC*COS(E))
      CR=ABS(E-EP)
      E=EP
      IF(CR.GT.CCR) GO TO 31
      R=1.-EC*COS(E)
      HE=.5*E
      W=2.*ATAN(ER*TAN(HE))
      SIND=SNI*SIN(W-EQ)
      DLT=ASIN(SIND)
      ALP=ASIN(TAN(DLT)*TINI)
      TST=COS(W-EQ)
      IF(TST.LT.0.0) ALP=PI-ALP
      IF(ALP.LT.0.) ALP=ALP+TPI
      SUN=TPI*(DATE-DELEQN)/YEAR
      IF(SUN.LT.0.0) SUN=SUN+TPI
      SLAG=SUN-ALP-.03255
C    COMPUTE HOUR ANGLE OF SUNSET AT ALL LATITUDES
      DO 1 I=0,JGMAX
      PH=ALAT(I)
      SS=SIN(PH)*SIN(DLT)
      CC=COS(PH)*COS(DLT)
      IF(DLT.EQ.0.0) GO TO 16
      AP=ABS(PH)
      EPS=ABS(AP-HPI)
      IF(EPS.GT.CCR) GO TO 14
      H=HPI*ABS(AP/PH+ABS(DLT)/DLT)
      GO TO 5
   14 AR=-SS/CC
      AC=ABS(AR)
      IF(AC-1.+CCR) 4,2,3
  2   H=(AC-AR)*HPI
      GO TO 5
  3   IF(AR.LT.0.) GO TO 25
      H=0.
      GO TO 5
 25   H=PI
      GO TO 5
 16   H=HPI
      GO TO 5
    4 H=ACOS(AR)
  5   HANG(I)=H
      TAUDA(I)=AMAX1(H/PI,0.0)
      IF(H.EQ.0.) GO TO 100
C     IF (I.EQ.1) WRITE (6,*) ' SS,CC,SIN(H),H,COSZ(I)'
      COSZ(I)=AMAX1((SS+CC*SIN(H)/H),0.0)
C     IF (I.EQ.1) WRITE (6,*) SS,CC,SIN(H),H,COSZ(I)
      GO TO 1
  100 COSZ(I)=0.0
    1 CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE RADZEN(FJD,DLT,SLAG,XLAT,DHR,NLNG,COSZRS,IDARK)
C ----------------------------------------------------------------------
C
CL            COMPUTE EFFECTIVE MEAN COSINE OF ZENITH ANGLE AND DAYLIGHT
CL            FRACTION FROM LATITUDE AND PARAMETERS GIVEN BY
CL            SUBROUTINE RADSLR
C
C             INPUT ARGUMENTS TO CIRCULAR FUNCTIONS ARE IN RADIANS
C
C-----------------------------------------------------------------------
C
CL            JULIAN DAY AND RELATED INFO FOR UPDATING BOUNDARY
CL            CONDITION DATA
C
      COMMON /COMDYR/
     +   CALDAY,  JYEAR,  JMNTH,   JDAY,  JHOUR,  JMINT,  ECCF
C
CJEN--------------------------------------------------------------------
      SAVE PI,TPI,PID24
C
      PARAMETER(N=144,JGMAX=88)
      REAL XLAT(0:JGMAX),COSZRS(N,0:JGMAX)
      LOGICAL IDARK(N,0:JGMAX)
      DATA PI/3.1415927/,TPI/6.2831853/,PID24/.13089969/
C
      CVPR=TPI/N
      GHA=FJD*TPI+SLAG
 
      DO 30 J=0,JGMAX
         SS=SIN(XLAT(J))*SIN(DLT)
         CC=COS(XLAT(J))*COS(DLT)
         DO 20 I = 1,N
            XLNG=CVPR*(I-1)
            HLOC=GHA+XLNG
            COSZRS(I,J)=SS+CC*COS(HLOC)
   20    CONTINUE
 30   CONTINUE
 
      DO 60 J=0,JGMAX
      DO 60 I=1,N
         IF (COSZRS(I,J).GE.0.) THEN
            IDARK(I,J) = .FALSE.
         ELSE
            IDARK(I,J) = .TRUE.
         END IF
 60   CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE CHEMINT(PROD1,PROD2,PROD,LOSS1,LOSS2,LOSS,
     ! THETA1,THETA2,THETA,TIME)
 
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (N=144,JGMAX=88,KTOP=19)
 
      REAL PROD1(0:JGMAX,KTOP), PROD2(0:JGMAX,KTOP), PROD(0:JGMAX,KTOP),
     !     LOSS1(0:JGMAX,KTOP), LOSS2(0:JGMAX,KTOP), LOSS(0:JGMAX,KTOP),
     !    THETA1(0:JGMAX,KTOP),THETA2(0:JGMAX,KTOP),THETA(0:JGMAX,KTOP)
 
      W2 =  MOD(TIME/86400.,15.)/15.
      WT2 = MOD(TIME/86400.,30.)/30.
 
      DO 10 K= 1,KTOP
         DO 5 J=0,JGMAX
            PROD(J,K) =  (1.-W2)*PROD1(J,K) + W2*PROD2 (J,K)
            LOSS(J,K) =  (1.-W2)*LOSS1(J,K) + W2*LOSS2 (J,K)
  5      CONTINUE
 10   CONTINUE
 
      DO 20 K= 1,KTOP
         DO 15 J=0,JGMAX
            THETA(J,K) =  (1.-WT2)*THETA1(J,K) + WT2*THETA2 (J,K)
 15      CONTINUE
 20   CONTINUE
 
      RETURN
      END
