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



C     PROGRAM TRAN19
C
      IMPLICIT REAL (A-H,O-Z)
 
C     DRIVER PROGRAM TO SOLVE TRACER CONTINUITY EQN. IN A SIGMA-
C     PRESSURE COORDINATE SYSTEM. HORIZONTAL ADV IS DONE WITH VANLEER'S
C     SCHEME.  VERTICAL ADV IS DONE WITH PRATHER'S SCHEME.
 
C     LAST MAJOR UPDATE 6/24/91.
C
      PARAMETER(N=144,JGMAX=88,NCON=1,JGMAX1=JGMAX+1)
      PARAMETER (KSTRAT=11,KTROP=8,KTOP=KSTRAT+KTROP,NTAPES=61)
      PARAMETER (KTROP2=8)
      PARAMETER (NJG=N*JGMAX,NJG2=N*JGMAX1)
      PARAMETER (ILMM=N,IJMM=JGMAX,JDIAG=5,JDIAG2=3,KDIAG2=3,NVARS=10)
      PARAMETER (LSLT1=1,LSLT2=8,LMLT1=9,LMLT2=79,LHLT1=80,LHLT2=87)
C
C     N       # OF EAST-WEST GRID POINTS  0E TO (360-360/N)E
C     JGMAX+1 # OF NORTH-SOUTH GRID POINTS 0S TO 90N
C     JGMAX   # OF NORTH-SOUTH GRID EDGES 87S TO 87N
C     KTOP    # OF VERTICAL LEVELS
C     KSTRAT  # OF VERTICAL LEVELS IN STRATOSPHERE
C     KTROP   # OF VERTICAL LEVELS IN TROPOSPHERE
C
C     PHYSICAL SPACE ARRAYS ARE DEFINED
      REAL CONST(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON)
      REAL EWFLUX(N,JDIAG,KTOP,NCON), NSFLUX(N,JDIAG,KTOP,NCON),
     1 VRFLUX (N,JDIAG,KTOP,NCON), CHEMPROD(N,JDIAG,KTOP,NCON),
     1 PSFTERM(N,JDIAG,KTROP2,NCON),
     1 CHEMLOSS(N,JDIAG,KTOP,NCON)
 
      REAL CZ(N,0:JGMAX,KTOP,NCON), CZZ(N,0:JGMAX,KTOP,NCON)
      REAL CZ1(N,0:JGMAX,KTOP), CZZ1(N,0:JGMAX,KTOP)
 
      real test(n,0:jgmax)      

      REAL F1(N,0:JGMAX,KTOP),DPSF(N,0:JGMAX),COSZ(N,0:JGMAX)
      REAL HANG(0:JGMAX)
      REAL F3(N,0:JGMAX,KTOP),SCALE(NCON)
C*C*C*C*C*
C     REAL HNO3ZM(0:JGMAX,KTOP)
C
      LOGICAL IDARK(N,0:JGMAX),DIURNAL
      LOGICAL RSTAR,FIRST,WNDSAV,CHEM,DIFF,UPSTREAM,DIAG
      INTEGER NSTEPS(NTAPES),IST1(NTAPES),NOBS7(NTAPES),NOBS8(NTAPES)
      INTEGER LUNCHM(NCON),NDIAG(JDIAG),WINDSOURCE,GRID
      INTEGER INDEX(NJG2*KTOP),NOBS9(NTAPES)
      INTEGER ILUN(NTAPES)
C
      REAL PSU(N,0:JGMAX),PSV(N,JGMAX),PSNEW(N,0:JGMAX)
      REAL PS1(N,0:JGMAX), PS2(N,0:JGMAX),PI1(N,0:JGMAX),
     ! U1  (N,0:JGMAX,KTOP), U2  (N,0:JGMAX,KTOP),
     ! V1  (N,JGMAX,KTOP), V2  (N,JGMAX,KTOP),SGBOT(N,0:JGMAX),
     ! SGTROP (N,0:JGMAX,KTROP2),SGSTRAT(N,0:JGMAX,0:KSTRAT-1),
     ! KEL1(N,0:JGMAX,KTOP), KEL2(N,0:JGMAX,KTOP),
     !  DIVINT(N,0:JGMAX,2),PS1INV(N,0:JGMAX)
      REAL WTC(N,0:JGMAX,0:KTOP)
C
      REAL COSCEN(0:JGMAX), DLATR(0:JGMAX),
     ! TWOLAT(0:JGMAX), MCOR(N,0:JGMAX), PLON(N)
      REAL DLATR2(JGMAX),COSEDGE(JGMAX)
 
C     CHEMICAL PRODUCTION AND LOSS ARRAYS
      REAL PROD (0:JGMAX,KTOP,NCON), LOSS (0:JGMAX,KTOP,NCON)
      REAL PROD1(0:JGMAX,KTOP,NCON), LOSS1(0:JGMAX,KTOP,NCON)
      REAL PROD2(0:JGMAX,KTOP,NCON), LOSS2(0:JGMAX,KTOP,NCON)
 
      REAL CL(0:JGMAX,KTOP),CLRAT(0:JGMAX,KTOP)
      REAL THETA(0:JGMAX,KTOP),THETA1(0:JGMAX,KTOP)
      REAL THETA2(0:JGMAX,KTOP)
C
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),
     2 SIGMA(KTOP),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),
     3 SGINT(0:KTOP)
C
      COMMON / THREE/HTSCALE,HINV,PI
 
      COMMON/RDVAR/RSTAR,FIRST,ISAVE,ITAP1,IST1,NSTEPS,
     ! DT,NTSTRC,I1,B,CHEM,SCALE,WINDSOURCE,DIFF,DFCOEF,NPLOTS,DIURNAL
     !,GRID,iread0
 
      COMMON / RDVAR2/ DIAG
      COMMON / DIAGNOSTIC / DD(N,JDIAG2,KDIAG2,NVARS), NNLAT(JDIAG2),
     1 NNLEV(KDIAG2)
 
C     JANUARY AND FEBRUARY COMBINES NOBS
C     MAIN BODY
       DATA NOBS7/10,10,11,11,10,11,10,11,10,10,
     1            10,13,11,12, 1,11,12,11,11,11,
     1            11,12,11, 5, 4, 9,11,10, 9,10,
     1             6, 7,10,10,11,11, 9,12,13,13,
     1             8,15,15,16,16,16,16,16,16,16,
     1            16,18,15,14,13,14,13,13,14,13,
     1             9/
C     1989
      DATA NOBS8/ 10,10,11,13,14,13,12,14,13,12,
     1            13,12,15,13,13,13,13,12,13,15,
     1            16,12,13,13,13,14,13,14,14, 2,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0/
 
C     ERIC'S GCM
      DATA NOBS9/ 23,24,24, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0/
 
      DATA ILUN/  51,52,53,54,55,56,57,58,59,60,
     1            61,62,63,64,65,66,67,68,69,70,
     1            71,72,73,74,75,76,77,78,79,80,
     1            81,82,83,84,85,86,87,88,89,90,
     1            91,92,93,94,95,96,97,98,99,60,
     1            61,62,63,64,65,66,67,68,69,70,
     1            71/
 
      DATA NDIAG/ 44, 49, 59, 74, 84 /
 
      DATA NNLAT/ 44, 59, 74/
 
      DATA NNLEV/  9, 13, 15/
 
      DATA POLEDGE / -87. /
 
      TIN    = SECOND()
C
      WRITE (6,*) 'BEGIN EXECUTION'
C     INITIALIZE TIMING PARAMETERS
      TIMIN  = 0.
      TIMCHM = 0.
      TIMDYN = 0.
      TIMDIA = 0.
      TIMX   = 0.
      TIMY   = 0.
      TIMZ   = 0.
      TIMDFF  = 0.
C     DEFINE VARIABLES AND OUTPUT EXPERIMENT HEADER TO UNIT10
      CALL SET19
      WRITE (6,*) ' IS THIS A RUN WITH CHEMISTRY ',CHEM
      WRITE (6,*) ' IS DIURNAL CHEMISTRY USED ', DIURNAL
      WRITE (6,*) ' IS THIS A RESTART ', RSTAR
      WRITE (6,*) ' ITAP1=',ITAP1, ' DT=',DT, ' NTSTRC=',NTSTRC
      WRITE (6,*) ' IST1(K)'
      WRITE (6,982) (IST1(IK),IK=1,NTAPES)
      WRITE (6,*) ' NSTEPS(IK)'
      WRITE (6,982) (NSTEPS(IK),IK=1,NTAPES)
 982  FORMAT (1X,6I6)
 
      PI        = 4.0   * ATAN(1.)
      PIINV     = 1./ PI
      HTSCALE   = 7000.
      HINV      = 1. / HTSCALE
      A         =  6.371220E6
 
      ICHEM     = -999
 
C     LOGICAL UNIT NUMBERS ARE SPECIFIED.
      DO 22 IC=1,NCON
         LUNCHM(IC) = 48 + IC
 22   CONTINUE
      LUNCHM2 = 47
      LNOON     =  33
      LMID      =  34
      LUNDIAG1 = 35
      LUNDIAG2 = 359
      WNDSAV    =  .true.
      NDT       = 5
      TDT       =  DT / NTSTRC
C     TDT2      =  TDT * 0.5
****************SPECIAL
      TDT2 = TDT
      TDT2I     = 1. / TDT2
 
C     DETERMINE TIME INTO INTEGRATION
      TWIND1 = 0.
      IF (WINDSOURCE.EQ.1)  THEN
         IF (ITAP1.GT.I1) THEN
            DO 20 I=I1,ITAP1-1
               TWIND1 = TWIND1 + NOBS7(I) * DT
 20         CONTINUE
         END IF
      ELSE IF (WINDSOURCE.EQ.2) THEN
         IF (ITAP1.GT.I1) THEN
            DO 25 I=I1,ITAP1-1
               TWIND1 = TWIND1 + NOBS8(I) * DT
 25         CONTINUE
         END IF
      ELSE IF (WINDSOURCE.GE.3) THEN
         IF (ITAP1.GT.I1) THEN
            DO 26 I=I1,ITAP1-1
               TWIND1 = TWIND1 + NOBS9(I) * DT
 26         CONTINUE
         END IF
      END IF
 
      IF (RSTAR) TWIND1 = TWIND1 - DT
 
      IF (WINDSOURCE.EQ.1) THEN
C*C*C*   -1 IF 12Z DEC 1 IS 1ST OBS OF ITAP1 (1979).
C        -7 IF 12Z JAN 12 IS 7TH OBS OF ITAP1(1979)
C        -7 IF 12Z DEC  5 IS 7TH OBS OF ITAP1 (1979)
         TWIND1 = TWIND1 + (IST1(ITAP1)-1) * DT
      ELSE IF (WINDSOURCE.EQ.2) THEN
C        -3 SINCE 00Z DEC 28 IS 3RD OBS OF ITAP1 (1989).
         TWIND1 = TWIND1 + (IST1(ITAP1)-3) * DT
      ELSE IF (WINDSOURCE.GE.3) THEN
C        -1 SINCE 18Z JAN 28 IS 1ST OBS OF ITAP1 (1989).
         TWIND1 = TWIND1 + (IST1(ITAP1)-1) * DT
      END IF
      TWIND2 = TWIND1
 
      DLON      =  360. / N
      DTIM2     =  (DLON / 720.) * 8.64E4
      DTIM      =  (DLON / 360.) * 8.64E4
      PLON(1)   =  DTIM2
      DO 40 IL = 2,N
         PLON(IL) = 86400. - DTIM2 - (IL-2) * DTIM
 40   CONTINUE
 
C     VERTICAL DIFFERENCING ARRAYS ARE SET UP
      IF (GRID.EQ.1) THEN
         CALL SETVER
      ELSE
         CALL VERGCM
      END IF
 
      WRITE (6,*) ' PINT = ', PINT, 'PTOP = ', PTOP
      WRITE (6,*) ' SIGMA,DSIG, AND DELTAS'
      DO 60  K=1,KTOP
         WRITE (6,58) K, SIGMA(K), DSIG(K), DELTAS(K)
 58      FORMAT(1X,I4,3(1PE12.5))
 60   CONTINUE
 
C     GRID IS DEFINED
      D1 = -2. * POLEDGE / (JGMAX-1)
      D1O2 = D1 * 0.5
      DTR = PI /180.
      ADPHI = D1 * DTR * A
 
      DLATR(0)         = -89.99 * DTR
      COSCEN(0)        = COS(DLATR(0))
      DO 65 J=1,JGMAX-1
         DLATR(J) = (POLEDGE + (J-1)*D1 + D1O2) * DTR
         COSCEN(J) = COS(DLATR(J))
 65   CONTINUE
      DLATR(JGMAX)     =  89.99 * DTR
      COSCEN(JGMAX)    = COS(DLATR(JGMAX))
 
      DO 75 J=1,JGMAX
         DLATR2(J) = (POLEDGE + (J-1)* D1) * DTR
         COSEDGE(J) = COS(DLATR2(J))
 75   CONTINUE
C
C     MCOR A MIXING RATIO TO MASS CONVERSION TERM IS DETERMINED.
      TWOLAT(0) = 90. + POLEDGE
      DO 77 IJ=1,JGMAX-1
         TWOLAT(IJ) = D1
 77   CONTINUE
      TWOLAT(JGMAX) = 90. + POLEDGE
 
      DO  80 J=1,JGMAX-1
      DO  80 I=1,N
         MCOR(I,J) = COSCEN(J) * TWOLAT(J) * DTR
 80   CONTINUE
C
      DO 90 I=1,N
         MCOR(I,0)     = 0.5 * N * SIN((1.5*D1*DTR))**2
         MCOR(I,JGMAX) = 0.5 * N * SIN((1.5*D1*DTR))**2
 90   CONTINUE
 
      WRITE (6,*) ' TRANSPORT LATITUDES '
      DO 100 J=1,JGMAX
         WRITE (6,95) j,DLATR2(J)/dtr,DLATR(J)/dtr,COSEDGE(J),COSCEN(J)
 100  CONTINUE
 95   format(i4,f7.2,f7.2,e10.5,e10.5)
      WRITE (6,*)
      TIMIN = TIMIN + SECOND() - TIN
 
C      TCH = SECOND()
C      IF (CHEM) THEN
*         COMMENT OUT AJRATE CALL UNLESS L4JAK.ZACONV.OUT CHEM USED
C         CALL AJRATE(LUNCHM(1),LOSS(0,1,1),PROD(0,1,1))
C      END IF
C      TIMCHM = TIMCHM + SECOND() - TCH
 
      TDYN = SECOND()
C     INITIAL FIELDS ARE READ IN  AND TAPE IS POSITIONED
      IF (RSTAR) THEN
C        ACCESS CONST FDS, TIME, AND TAPE POSITIONING PARAMETER
         CALL RESTARP(CONST,CZ,CZZ,FD0,TIME,IREAD,ISAVE)
C        CALL RESTARPB(CONST,CZ,CZZ,FD0,TIME,IREAD,ISAVE,HNO3ZM)
C        ONE OBS IS READ BOTH IN THE 1ST RUN + IN THE SUBSEQUENT RESTART
         IREAD = IREAD - 1
         WRITE (6,*) ' IREAD = ',IREAD
         IF (WINDSOURCE.EQ.1) THEN
            CALL RDOV19(IREAD,ILUN,LUN,ISK)
            CALL RDSIM5(LUN,PS1,U1,V1,KEL1)
         ELSE IF (WINDSOURCE.EQ.2) THEN
            CALL RDOV1989(IREAD,ILUN,LUN,ISK)
            CALL RDSIM5(LUN,PS1,U1,V1,KEL1)
         ELSE IF (WINDSOURCE.EQ.3) THEN
            CALL RDOVGCM(IREAD,ILUN,LUN,ISK)
            CALL RDGCM(LUN,PS1,U1,V1,KEL1)
         ELSE IF (WINDSOURCE.EQ.4) THEN
            CALL RDSCH(PS1,U1,V1,KEL1,COSEDGE)
         ELSE IF (WINDSOURCE.EQ.5) THEN
            CALL WIND(U1,V1,COSCEN,B,DLATR,PS1,KEL1)
         END IF
 
         IREAD = IREAD + 1
 
         TTEST = ABS(TIME - TWIND1)
         IF (TTEST.GT.0.01) THEN
            WRITE (6,*) ' TIMES ARE NOT EQUAL'
            WRITE (6,*) ' TIME= ',TIME, ' TWIND1 = ',TWIND1
            AAA = 0.
            STOP
         END IF
      ELSE
         iread = iread0 
          IF (WINDSOURCE.EQ.1) THEN
             CALL RDOV19(IREAD,ILUN,LUN,ISK)
             CALL RDSIM5(LUN,PS1,U1,V1,KEL1)
          ELSE IF (WINDSOURCE.EQ.2) THEN
             CALL RDOV1989(IREAD,ILUN,LUN,ISK)
             CALL RDSIM5(LUN,PS1,U1,V1,KEL1)
          ELSE IF (WINDSOURCE.EQ.3) THEN
             CALL RDOVGCM(IREAD,ILUN,LUN,ISK)
             CALL RDGCM(LUN,PS1,U1,V1,KEL1)
          ELSE IF (WINDSOURCE.EQ.4) THEN
             CALL RDSCH(PS1,U1,V1,KEL1,COSEDGE)
             CALL BLOB(CONST,FD0(1,0,1))
          ELSE IF (WINDSOURCE.EQ.5) THEN
             CALL WIND(U1,V1,COSCEN,B,DLATR,PS1,KEL1)
             CALL PATTERN(COSCEN,CONST,FD0(1,0,1))
          END IF
C          DO 7004 J=0,JGMAX
C          DO 7004 I=1,N
C             PS1(I,J) = 1000.
C 7004     CONTINUE
 
         IST1(ITAP1) = IST1(ITAP1) + 1
         IREAD =  IREAD + 1
C        INITIALIZE CONSTITUENTS.
         LUNHNO = 45
         IC = 1
C        CALL TINIT(CONST(1,0,1,IC),LUNHNO,SCALE(IC),FD0(1,0,IC),IC)
         CALL INITO3(CONST(1,0,1,IC),LUNHNO,SCALE(IC),FD0(1,0,IC))
C        CALL HNO3IC(LUNHNO,FD0(1,0,IC),CONST(1,0,1,IC),HNO3ZM)
C         CALL ERIC(FD0(1,0,IC),CONST(1,0,1,IC),DLATR)
         DO 4 K=1,KTOP
         DO 4 J=0,JGMAX
         DO 4 I=1,N
ccc            CONST(I,J,K,1) = 10. 
ccc            FD0  (I,J,  1) = 10.  
            CZ   (I,J,K,1) = 0.
            CZZ  (I,J,K,1) = 0.
   4    CONTINUE
****************HCL INITIAL CONDITION************
C        DO 5 IC=2,NCON
C        DO 5 K=1,KTOP
C        DO 5 J=0,JGMAX
C        DO 5 I=1,N
C           CONST(I,J,K,IC) = CONST(I,J,K,1)
C           CZ   (I,J,K,IC) = CZ   (I,J,K,1)
C           CZZ  (I,J,K,IC) = CZZ  (I,J,K,1)
C           FD0(I,J,IC) = FD0(I,J,1)
C  5    CONTINUE
**********END HCL INITIAL CONDITION**************
C        INITIAL CONDITIONS FOR OLD 1979 O3 EXPERIMENT (ALSO 79 HNO3)
C         CALL TINIT2(CONST,LUNHNO,16,FD0(1,0,1))
 
         TIME   =  0.0
         IF (TIME.NE.TWIND1) THEN
            WRITE (6,*) ' TIMES ARE NOT EQUAL'
            WRITE (6,*) ' TIME= ',TIME, ' TWIND1 = ',TWIND1
            STOP
         END IF
 
         IF (DIAG) THEN
            DO 130 IC=1,NCON
            DO 130 K=1,KTROP
            DO 130 J=1,JDIAG
            DO 130 I=1,N
               EWFLUX  (I,J,K,IC) = 0.
               NSFLUX  (I,J,K,IC) = 0.
               VRFLUX  (I,J,K,IC) = 0.
               CHEMLOSS(I,J,K,IC) = 0.
               CHEMPROD(I,J,K,IC) = 0.
               PSFTERM (I,J,K,IC) = 0.
 130       CONTINUE
            DO 132 IC=1,NCON
            DO 132 K=KTROP+1,KTOP
            DO 132 J=1,JDIAG
            DO 132 I=1,N
               EWFLUX  (I,J,K,IC) = 0.
               NSFLUX  (I,J,K,IC) = 0.
               VRFLUX  (I,J,K,IC) = 0.
               CHEMLOSS(I,J,K,IC) = 0.
               CHEMPROD(I,J,K,IC) = 0.
 132       CONTINUE
         END IF
 
      END IF
 
      TIMDYN = TIMDYN + SECOND() - TDYN
C     BEGIN TIME STEPPING
      ICOUNT = -1
      DO 520 ITAPE = ITAP1,NTAPES
      LUN = ILUN(ITAPE)
      DO 500 ISTDY = IST1(ITAPE),NSTEPS(ITAPE)
C        READ IN WIND FIELD AT NEXT AVAIL TIME STEP
 
         TDYN = SECOND()
         IF ((WINDSOURCE.EQ.1).OR.(WINDSOURCE.EQ.2)) THEN
            CALL RDSIM5(LUN,PS2,U2,V2,KEL2)
         ELSE IF (WINDSOURCE.EQ.3) THEN
            CALL RDGCM (LUN,PS2,U2,V2,KEL2)
         ELSE IF (WINDSOURCE.EQ.4) THEN
            CALL RDSCH(PS2,U2,V2,KEL2,COSEDGE)
         ELSE IF (WINDSOURCE.EQ.5) THEN
            CALL WIND(U2,V2,COSCEN,B,DLATR,PS2,KEL2)
         END IF
C          DO 4002 J=0,JGMAX
C          DO 4002 I=1,N
C             PS2(I,J) = 1000.
C 4002     CONTINUE
 
         TWIND2 = TWIND2 + DT
         IREAD = IREAD + 1
         TIMDYN = TIMDYN + SECOND() - TDYN
C
*************READ IN CHEMISTRY FIELDS****************
C        FOR ALL BUT THE N2O EXPERIMENT USE FOLLOWING CHEM CODE
          TCHM = SECOND()
          IF (CHEM) THEN
C            DAY0 = -4 FOR DEC 28 STARTUP.  DAY0 = 0 FOR JAN 1 STARTUP.
C            AND DAY0 = 4 FOR JAN 5 STARTUP.
             DAY0 = -4. 
            IF (WINDSOURCE.EQ.1) DAY0 = 0.
**************************
c             DAY0 = 4.
C            CHEM SKIP PARAMETER FOR O3,HCL,N2O, AND NOY EXPS.
c             day0 = -4. + 30 
             ICHEM2 = INT((DAY0 + TIME/86400.)/15.)
C            WARNING OXPLF.DAT INPUT O3 FILE DOES NOT CONTAIN HEADER!!!!
C            MOST RECENT (4/10/91) O3 FILE DOES CONTAIN HEADER
C             ICHEM2 = INT((DAY0+TIME/86400.)/30.)
 
             IF (ICHEM2.NE.ICHEM) THEN
                CALL O3IN(PROD1(0,1,1),PROD2(0,1,1),LOSS1(0,1,1),
     1       LOSS2(0,1,1),THETA1,THETA2,ICHEM2,LUNCHM,LUNCHM2,DIURNAL)
C               CALL HNO3IN(DLATR,PROD,PI)
C               CALL CHEMIN(PROD,LOSS,LUNCHM,ICHEM2,LUNCHM2,THETA)
C               CALL HCLIN(PROD,LOSS,CL,CLRAT,LUNCHM,ICHEM2)
C               CALL N2OIN (PROD(0,1,1),LOSS(0,1,1),LUNCHM(1),ICHEM2)
C               CALL NOYIN (PROD(0,1,2),LOSS(0,1,2),LUNCHM(2),ICHEM2)
                ICHEM = ICHEM2
             END IF
          END IF
         TIMCHM = TIMCHM + SECOND() - TCHM
         LUNDIAG2 = LUNDIAG2 + 1
 
         DO 460 ISTCO = 1,NTSTRC
c         DO 460 ISTCO = 1,1
 
         TCHM = SECOND()
         IF (DIURNAL) THEN
C           DATE EXPERIMENT BEGINS IS HARDWIRED IN.
            MCDATE = 781201 + (43200.+TIME+0.01) / 86400.
            MMM   = TIME + 0.01
            MCSEC =  MOD(43200+MMM,86400)
            JYEAR = 1900 + MCDATE/10000
            JMNTH = MOD(MCDATE/100,100)
            JDAY  = MOD(MCDATE,100)
            JHOUR = MCSEC/3600
            DHR = 24.
            JMINT = MOD(MCSEC/60,60)
C           NCAR CCM SUBS ARE CALLED TO FIND WHERE IT IS DAYLIGHT.
            CALL JULDAY(JYEAR,JMNTH,JDAY,JHOUR,JMINT,JD,FJD)
            CALL RADSLR(JD,FJD,R,DLT,ALP,SLAG,DLATR)
            CALL RADZEN(FJD,DLT,SLAG,DLATR,DHR,N,COSZ,IDARK)
            IF (ISTCO.EQ.1) THEN
               WRITE (6,*)
               DO 135 J=0,JGMAX
                  WRITE (6,136) (IDARK(I,J),I=1,N,2)
 135           CONTINUE
 136           FORMAT(72L1)
            END IF
         END IF
 
C        SPECIAL SECTION FOR OZONE CALCULATION.
         TCHM=SECOND()
         IF (CHEM) CALL CHEMINT(PROD1(0,1,1),PROD2(0,1,1),PROD(0,1,1),
     1  LOSS1(0,1,1),LOSS2(0,1,1),LOSS(0,1,1),THETA1,THETA2,THETA,TIME)
c         if (istco.eq.1) write (96) prod,loss,theta
         TIMCHM = TIMCHM + SECOND() - TCHM
 
         TDYN = SECOND()
         TIME2 = TIME + TDT
C        INTERPOLATE TO FIND DYNAMIC FIELDS AT CURRENT TIME (N).
C        NOTE: PSNEW IS THE SURFACE PRESSURE AT TIME N+1! BY USING
C        PSNEW WE AVOID INTEGRATING TO DETERMINE THE CHANGE IN PSF.
         CALL WDINT4 (U1,U2,V1,V2,
     !     PS1,PS2,KEL1,KEL2,TIME,TWIND1,TWIND2,PSNEW,TIME2)
         DO 140 J=0,JGMAX
         DO 140 I=1,N
            DPSF(I,J) = PSNEW(I,J) - PS1(I,J)
 140     CONTINUE
 
         TWIND1 = TIME
         TIME   = TIME + TDT
C        INTERPOLATE IN SPACE TO DETERMINE SFC PRESSURE AT EDGES
         CALL SURFPS(PS1,PSU,PSV)
 
C        DETERMINE "PI" ALONG THE BOX EDGES. AND 1/"PI" IN BOX CENTERS
C        "PI" IS NEEDED TO ENSURE CONSERVATION OF MASS IN TROPOSPHERE
C        WHERE THE HEIGHT OF BOXES VARIES WITH X AND Y.
         DO 145 IJ=1,IJMM
         DO 145 IL=1,ILMM
            PSV(IL,IJ) = PSV(IL,IJ) - PINT
 145     CONTINUE
         DO 150 IJ=0,IJMM
         DO 150 IL=1,ILMM
            PSU(IL,IJ) = PSU(IL,IJ) - PINT
            PS1INV(IL,IJ) = 1. / (PS1(IL,IJ)-PINT)
 150     CONTINUE

C        CALCULATE VERTICAL WIND RESIDUALLY.
         TX = SECOND()
C        1. CALCULATE X ADV TERM ON CONSTANT FIELD (F3) AND PUT IN F1.
         DO 160 K=1,KTOP
            DO 155 IJ=0,IJMM
            DO 155 IL=1,ILMM
               F3(IL,IJ,K) = 1.0
 155        CONTINUE
 160     CONTINUE
 
C        CHANGE DUE TO X ADV AT SOUTHERNMOST 8 LATS IS PUT IN F1.
 2004    FORMAT(10E12.5)
         UPSTREAM = .FALSE.
         CALL XADVHLAT(F3,COSCEN,TDT2,U1,F1,PS1INV,PSU,LSLT1,LSLT2,NDT,
     1    UPSTREAM,.FALSE.)
 
C        CHANGE DUE TO X ADV AT MIDDLE LATS IS PUT IN F1
         UPSTREAM = .FALSE.
         CALL XADV(F3,COSCEN,TDT2,U1,F1,PS1INV,PSU,LMLT1,LMLT2,UPSTREAM,
     1    .FALSE.)
 
C        CHANGE DUE TO X ADV AT NORTHERNMMOST 8 LATS IS PUT IN F1.
         UPSTREAM = .FALSE.
         CALL XADVHLAT(F3,COSCEN,TDT2,U1,F1,PS1INV,PSU,LHLT1,LHLT2,NDT,
     1    UPSTREAM,.FALSE.)
         TIMX = TIMX + SECOND() - TX
 2002    FORMAT(10(1PE12.5))
 
         TY = SECOND()
C        2. CALCULATE Y ADVECTION ON CONSTANT FIELD AND PLACE IN F3.
C        SUBROUTINE YADV REDUCES TO YADV2 FOR A CONSTANT FIELD.
         CALL YADV2(COSCEN,COSEDGE,V1,ADPHI,F3,PS1INV,PSV,TDT2,POLEDGE)
         TIMY = TIMY + SECOND() - TY
 
C        SUM X AND Y ADVECTION TERMS.
         DO 180 K=1,KTOP
            DO 175 IJ=0,IJMM
            DO 175 IL=1,ILMM
               F1(IL,IJ,K) = -(F1(IL,IJ,K) + F3(IL,IJ,K))
 175        CONTINUE
 180     CONTINUE
C        ADD IN SURFACE PRESSURE CHANGE TERM.
         DO 185 K=1,KTROP
         DO 185 J=0,JGMAX
         DO 185 I=1,N
            F1(I,J,K) = F1(I,J,K) + DPSF(I,J) /  PS1(I,J)
 185     CONTINUE
 
         TZ = SECOND()
C        DETERMINE VERTICAL VELOCITY RESIDUALLY.
         CALL SGPRA(F1,SGBOT,SGTROP,SGSTRAT,TDT2,PS1,WTC)
         TIMZ = TIMZ + SECOND() - TZ
         TIMDYN = TIMDYN + SECOND() - TDYN
 
         TDIA = SECOND()
C        U,V,SG, AND TEMPERATURE FIELDS ARE SAVED.
         IF ((WNDSAV).AND.(ISTCO.EQ.1)) THEN
             ILON = 0
             LUNWND = 90 + (ICOUNT+1)/8
****       SPECIAL DIAGNOSTIC CHANGE
             CALL HEADER (TWIND1,TDT,LUNWND,ILON)
             DO 190 IJ=0,IJMM
             DO 190 IL=1,ILMM
                F1(IL,IJ,1) = SGBOT(IL,IJ) * (PS1(IL,IJ)-PINT)
 190         CONTINUE
             DO 220  K=1,KTROP
                KP = K + 1
                DO 220 IJ=0,IJMM
                DO 220 IL=1,ILMM
                   F1(IL,IJ,KP)=SGTROP(IL,IJ,K)*(PS1(IL,IJ)-PINT)
 220         CONTINUE
 
             IF (KTROP.GT.0) THEN
                DO 240   K=1,KSTRAT-1
                   KP = K + KTROP+1
                   DO 240 IJ=0,IJMM
                   DO 240 IL=1,ILMM
                      F1(IL,IJ,KP)=SGSTRAT(IL,IJ,K) * (PINT-PTOP)
 240            CONTINUE
            ELSE
                DO 241   K=1,KSTRAT-1
                   KP = K + KTROP+1
                   DO 241 IJ=0,IJMM
                   DO 241 IL=1,ILMM
                      F1(IL,IJ,KP)=SGSTRAT(IL,IJ,K) * (PS1(IL,IJ)-PINT)
 241            CONTINUE
            END IF
 
*******************SPECIAL DIAGNOSTIC CHANGE******** 
             CALL PL3D(U1,PS1,1,LUNWND)
             CALL PLXD(V1,LUNWND)
             CALL PL3D(F1,PS1,2,LUNWND)
****************************************************
         END IF
 
C        THE MASS OF CONSTITUENTS ARE DETERMINED.
         IF (ISTCO.EQ.1) THEN
            WRITE (6,314)  TWIND1/ 86400.
 314        FORMAT (1X,' MASS OF CONSTS AT TIME ',F9.3)
C           WRITE (6,*) ' ABOUT TO CALL MASS3'
            CALL MASS3(CONST,PS1,MCOR,6)
         END IF
 
C        THE CONSTITUENT FIELDS ARE SAVED.
          MPLOT = MOD(ISTCO-1,72)
C         MPLOT = MOD(ISTCO-1,24)
         IF ((MPLOT.EQ.0).AND.(FIRST)) THEN
            ICOUNT = ICOUNT + 1
            LUNOUT = 10 + ICOUNT/NPLOTS
            ILON = 0
            CALL HEADER (TWIND1,TDT,LUNOUT,ILON)
            DO 2601 IC=1,NCON
               CALL PL3D (CONST(1,0,1,IC),PS1,IC,LUNOUT)
 2601        CONTINUE
         ELSE IF (ISTCO.EQ.1) THEN
            FIRST=.TRUE.
            LUNOUT = 10
         END IF
C
C        THE LONG NEAREST LOCAL NOON + MIDNIGHT IS PLOTTED EACH HR.
C           IF (MOD(TIME+0.01,3600.).LT.0.1) THEN
C              DO 280 IC=1,NCON
C                 CALL PLNON2(CONST(1,0,1,IC),PS1,IC,TWIND1,TDT,PLON,
C     !          LNOON,LMID)
C  280         CONTINUE
C           END IF
         TIMDIA = TIMDIA + SECOND() - TDIA
 
         TDYN = SECOND()
         DO 350 ILOOP = 1,1
         DO 340 IC=1,NCON
 
C           DIAGNOSTIC TERMS ARE SAVED.
            IF (DIAG) THEN
               DO 282 IVAR=1,NVARS
               DO 282 KD=1,KDIAG2
               DO 282 JD=1,JDIAG2
               DO 282 I=1,N
                  DD(I,JD,KD,IVAR) = 0.
 282           CONTINUE
 
C              SAVE MIXING RATIO AT TIME N.
               DO 283 KD=1,KDIAG2
               DO 283 JD=1,JDIAG2
               DO 283 I=1,N
                  DD(I,JD,KD,5) = CONST(I,NNLAT(JD),NNLEV(KD),IC)
 283           CONTINUE
 
C              CALCULATE MAGNITUDE OF PSF CHANGE TERM FOR DIAG PURPOSES.
               DO 284 K=1,KTROP
               DO 284 J=0,JGMAX
               DO 284 I=1,N
                  F3(I,J,K) = -CONST(I,J,K,IC)*DPSF(I,J)/PS1(I,J)
 284           CONTINUE
               CALL DIAG2(F3,CONST(1,0,1,IC),4,EWFLUX(1,1,1,IC),
     1         NSFLUX(1,1,1,IC),VRFLUX(1,1,1,IC),
     1         PSFTERM(1,1,1,IC),NDIAG)
               DO 287 K=1,KDIAG2
                  IF (NNLEV(K).LE.KTROP) THEN
                     DO 285 J=1,JDIAG2
                     DO 285 I=1,N
                        DD(I,J,K,4) = F3(I,NNLAT(J),NNLEV(K))
 285                 CONTINUE
                  END IF
 287           CONTINUE
            END IF
 
C           CHANGE DUE TO XADV IS DETERMINED AND PUT IN F1.
            TX = SECOND()
            UPSTREAM = .FALSE.
            CALL XADVHLAT(CONST(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       LSLT1,LSLT2,NDT,UPSTREAM,DIAG)
            UPSTREAM = .FALSE.
           CALL XADV(CONST(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       LMLT1,LMLT2,UPSTREAM,DIAG)
            UPSTREAM = .FALSE.
            CALL XADVHLAT(CONST(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       LHLT1,LHLT2,NDT,UPSTREAM,DIAG)
            TIMX = TIMX + SECOND() - TX
 
C           CHANGE DUE TO YADV IS DETERMINED AND PUT IN F3.
            TY = SECOND()
            UPSTREAM = .FALSE.
            CALL YADV(CONST(1,0,1,IC),COSCEN,COSEDGE,TDT2,V1,ADPHI,
     1      F3,PS1INV,PSV,UPSTREAM,DIAG,POLEDGE)
 
            IF (DIAG) THEN
               CALL DIAG2(F1,CONST(1,0,1,IC),1,EWFLUX(1,1,1,IC),
     1          NSFLUX(1,1,1,IC),VRFLUX(1,1,1,IC),PSFTERM(1,1,1,IC),
     1          NDIAG)
 
               CALL DIAG2(F3,CONST(1,0,1,IC),2,EWFLUX(1,1,1,IC),
     1          NSFLUX(1,1,1,IC),VRFLUX(1,1,1,IC),PSFTERM(1,1,1,IC),
     1          NDIAG)
 
               DO 288 K=1,KDIAG2
               DO 288 J=1,JDIAG2
               DO 288 I=1,N
                  DD(I,J,K,1) = F1(I,NNLAT(J),NNLEV(K))
 288           CONTINUE
               DO 289 K=1,KDIAG2
               DO 289 J=1,JDIAG2
               DO 289 I=1,N
                  DD(I,J,K,2) = F3(I,NNLAT(J),NNLEV(K))
 289           CONTINUE
            END IF
 
C           WRITE (6,*) ' X ADV TERM '
C           WRITE (6,2002) (F1(II,JGMAX,14),II=1,N)
C           WRITE (6,*) ' Y ADV TERM '
C           WRITE (6,2002) (F3(II,JGMAX,14),II=1,N)
 
C           ADD X AND Y ADVECTION TERMS TOGETHER.
            DO 290 K=1,KTOP
            DO 290 J=0,JGMAX
            DO 290 I=1,N
               F1(I,J,K) = F1(I,J,K) + F3(I,J,K)
 290        CONTINUE
 
            TIMY = TIMY + SECOND() - TY
 
C           CHANGE DUE TO VERTICAL ADV IS PUT INTO F3.
            TZ = SECOND()
            CALL ZADV(CONST(1,0,1,IC),CZ(1,0,1,IC),CZZ(1,0,1,IC),
     1       FD0(1,0,IC),SGSTRAT,SGTROP,SGBOT,TDT2,F3,CZ1,CZZ1,WTC,
     1       DIAG,PS1)
C           WRITE (6,*) ' SGSTRAT (5) '
C           WRITE (6,2002) (SGSTRAT(II,JGMAX,5),II=1,N)
C           WRITE (6,*) ' SGSTRAT (6) '
C           WRITE (6,2002) (SGSTRAT(II,JGMAX,6),II=1,N)
C           WRITE (6,*) ' Z ADV TERM '
C           WRITE (6,2002) (F3(II,JGMAX,14),II=1,N)
            TIMZ = TIMZ + SECOND() - TZ
 
            IF (DIAG) THEN
               CALL DIAG2(F3,CONST(1,0,1,IC),3,EWFLUX(1,1,1,IC),
     1          NSFLUX(1,1,1,IC),VRFLUX(1,1,1,IC),PSFTERM(1,1,1,IC),
     1          NDIAG)
               DO 295 K=1,KDIAG2
               DO 295 J=1,JDIAG2
               DO 295 I=1,N
                  DD(I,J,K,3) = F3(I,NNLAT(J),NNLEV(K))
 295           CONTINUE
            END IF
 
C          ADVECTION TERMS ARE ADDED IN. ADJ FOR PSF CHANGE IS MADE
           DO 300 K=1,KTROP
           DO 300 IJ=0,IJMM
           DO 300 IL=1,ILMM
              CONST(IL,IJ,K,IC)= CONST(IL,IJ,K,IC)+ F1(IL,IJ,K) +
     1         F3(IL,IJ,K) - CONST(IL,IJ,K,IC)*DPSF(IL,IJ) / PS1(IL,IJ)
 300       CONTINUE
 
            DO 320 K=KTROP+1,KTOP
            DO 320 IJ=0,IJMM
            DO 320 IL=1,ILMM
               CONST(IL,IJ,K,IC)=CONST(IL,IJ,K,IC) +
     1            F1(IL,IJ,K)+ F3(IL,IJ,K)
 320        CONTINUE
C           WRITE (6,*) ' CONST TERM AFTER ADV'
C           WRITE (6,2002) (CONST(II,JGMAX,14,IC),II=1,N)
 
C           CONVERT MOMENTS TO MIXING RATIO UNITS.
            DO 322 K=1,KTOP
            DO 322 J=0,JGMAX
            DO 322 I=1,N
               CZ  (I,J,K,IC) = CZ  (I,J,K,IC) / WTC(I,J,K)
               CZZ (I,J,K,IC) = CZZ (I,J,K,IC) / WTC(I,J,K)
               CZ1 (I,J,K   ) = CZ1 (I,J,K   ) / WTC(I,J,K)
               CZZ1(I,J,K   ) = CZZ1(I,J,K   ) / WTC(I,J,K)
 322        CONTINUE
C           WRITE (6,*) ' CZ '
C           WRITE (6,2002) (CZ(II,JGMAX,14,IC),II=1,N)
C           WRITE (6,*) ' CZZ TERM '
C           WRITE (6,2002) (CZZ(II,JGMAX,14,IC),II=1,N)
C           WRITE (6,*) ' CZ1'
C           WRITE (6,2002) (CZ1(II,JGMAX,14),II=1,N)
C           WRITE (6,*) ' CZZ1 TERM '
C           WRITE (6,2002) (CZZ1(II,JGMAX,14),II=1,N)
 
C           EAST-WEST ADVECTION OF 1ST VERTICAL MOMENT IS CALCULATED.
            TX = SECOND()
            UPSTREAM = .TRUE.
            CALL XADVHLAT(CZ(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       LSLT1,LSLT2,NDT,UPSTREAM,.FALSE.)
            CALL XADV(CZ(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       LMLT1,LMLT2,UPSTREAM,.FALSE.)
            CALL XADVHLAT(CZ(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1                         LHLT1,LHLT2,NDT,UPSTREAM,.FALSE.)
            TIMX = TIMX + SECOND() - TX
 
C           NORTH-SOUTH ADVECTION OF FIRST VERTICAL MOMENT.
            TY = SECOND()
            CALL YADV(CZ(1,0,1,IC),COSCEN,COSEDGE,TDT2,V1,ADPHI,
     1      F3,PS1INV,PSV,UPSTREAM,.FALSE.,POLEDGE)
 
C           WRITE (6,*) ' F1 FOR 1ST MOMENT ADV '
C           WRITE (6,2002) (F1(II,JGMAX,14),II=1,N)
C           WRITE (6,*) ' F3 FOR 1ST MOMENT ADV '
C           WRITE (6,2002) (F3(II,JGMAX,14),II=1,N)
C          ADD X AND Y ADV OF 1ST MOM TO 1ST MOM AFTER ADDING IN VERT
C          ADV TERM. IN THE TROPOSPHERE ALSO ADD IN SFC PRESSURE TERM.
           DO 325 K=1,KTROP
           DO 325 IJ=0,IJMM
           DO 325 IL=1,ILMM
              CZ(IL,IJ,K,IC)= CZ1(IL,IJ,K)+ F1(IL,IJ,K) +
     1         F3(IL,IJ,K) - CZ1(IL,IJ,K)*DPSF(IL,IJ) / PS1(IL,IJ)
 325       CONTINUE
           DO 328 K=KTROP+1,KTOP
           DO 328 IJ=0,IJMM
           DO 328 IL=1,ILMM
              CZ(IL,IJ,K,IC)=CZ1(IL,IJ,K) +
     1         F1(IL,IJ,K)+ F3(IL,IJ,K)
 328       CONTINUE
C           WRITE (6,*) ' CZ AFTER ADV '
C           WRITE (6,2002) (CZ(II,JGMAX,14,IC),II=1,N)
 
C           EAST-WEST ADVECTION OF SECOND VERTICAL MOMENT IS CALCULATED.
            TX = SECOND()
            UPSTREAM = .TRUE.
            CALL XADVHLAT(CZZ(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       LSLT1,LSLT2,NDT,UPSTREAM,.FALSE.)
            CALL XADV(CZZ(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       LMLT1,LMLT2,UPSTREAM,.FALSE.)
            CALL XADVHLAT(CZZ(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1                            LHLT1,LHLT2,NDT,UPSTREAM,.FALSE.)
            TIMX = TIMX + SECOND() - TX
 
C           CALCULATE NORTH-SOUTH ADVECTION OF 2ND VERTICAL MOMENT.
            TY = SECOND()
            CALL YADV(CZZ(1,0,1,IC),COSCEN,COSEDGE,TDT2,V1,ADPHI,
     1      F3,PS1INV,PSV,UPSTREAM,.FALSE.,POLEDGE)
 
C           WRITE (6,*) ' F1 FOR 2ND MOMENT ADV '
C           WRITE (6,2002) (F1(II,JGMAX,14),II=1,N)
C           WRITE (6,*) ' F3 FOR 2ND MOMENT ADV '
C           WRITE (6,2002) (F3(II,JGMAX,14),II=1,N)
 
C           ADD X AND Y ADV OF 1ST MOM TO 1ST MOM AFTER ADDING IN VERT
C           ADV TERM. IN THE TROPOSPHERE ALSO ADD IN SFC PRESSURE TERM.
            DO 330 K=1,KTROP
            DO 330 IJ=0,IJMM
            DO 330 IL=1,ILMM
              CZZ(IL,IJ,K,IC)= CZZ1(IL,IJ,K)+ F1(IL,IJ,K) +
     1         F3(IL,IJ,K) - CZZ1(IL,IJ,K)*DPSF(IL,IJ) / PS1(IL,IJ)
 330        CONTINUE
            DO 335 K=KTROP+1,KTOP
            DO 335 IJ=0,IJMM
            DO 335 IL=1,ILMM
               CZZ(IL,IJ,K,IC)=CZZ1(IL,IJ,K) +
     1          F1(IL,IJ,K)+ F3(IL,IJ,K)
 335        CONTINUE
C           WRITE (6,*) ' CZZ AFTER ADV '
C           WRITE (6,2002) (CZZ(II,JGMAX,14,IC),II=1,N)
C           CONVERT MOMENTS BACK TO MASS.
            DO 337 K=1,KTOP
            DO 337 J=0,JGMAX
            DO 337 I=1,N
               CZ  (I,J,K,IC) = CZ  (I,J,K,IC) * WTC(I,J,K)
               CZZ (I,J,K,IC) = CZZ (I,J,K,IC) * WTC(I,J,K)
 337        CONTINUE
 340     CONTINUE
 350     CONTINUE
         TIMDYN = TIMDYN + SECOND() - TDYN
 
C        CHEMISTRY IS PERFORMED.
         TCH = SECOND()
         IF (CHEM) THEN
            DO 360 IC=1,NCON
               CALL NEGS(CONST(1,0,1,IC),IC,TIME)
 360        CONTINUE
            DO 380 IC=1,NCON
            DO 380 K=1,KTOP
            DO 380 J=0,JGMAX
            DO 380 I=1,N
               IF (CONST(I,J,K,IC).LT.0.0)
     1          CONST(I,J,K,IC)=0.0
 380        CONTINUE
C           CALL HCLCHEM(CONST,KEL1,PROD,LOSS,CL,CLRAT,TDT)
C           CALL CHEMDO(CONST,PROD,LOSS,TDT)
            CALL CHEMO3(CONST,PROD,LOSS,TDT,CHEMPROD,CHEMLOSS,NDIAG,
     1        DIAG,KEL1,THETA,DIURNAL,IDARK,TIME)
C           CALL HNO3CHEM(CONST,HNO3ZM,PROD)
         END IF
         TIMCHM = TIMCHM + SECOND() - TCH
 
         TIMD = SECOND()
         IF (DIFF) THEN
            DO 450 IC=1,NCON
               CALL DIFFX(CONST(1,0,1,IC),COSCEN,TDT2,F1,PS1INV,DFCOEF)
               DO 390 K=1,KTOP
               DO 390 J=0,JGMAX
               DO 390 I=1,N
                  CONST(I,J,K,IC) = CONST(I,J,K,IC) +
     1              F1(I,J,K)
 390           CONTINUE
               CALL DIFFY(CONST(1,0,1,IC),COSCEN,COSEDGE,TDT,ADPHI,F1,
     1            PS1INV,DFCOEF)
               DO 420 K=1,KTOP
               DO 420 J=0,JGMAX
               DO 420 I=1,N
                  CONST(I,J,K,IC) = CONST(I,J,K,IC) +
     1              F1(I,J,K)
 420           CONTINUE
               CALL DIFFX(CONST(1,0,1,IC),COSCEN,TDT2,F1,PS1INV,DFCOEF)
               DO 440 K=1,KTOP
               DO 440 J=0,JGMAX
               DO 440 I=1,N
                  CONST(I,J,K,IC) = CONST(I,J,K,IC) +
     1              F1(I,J,K)
 440           CONTINUE
 450        CONTINUE
         END IF
         TIMDFF = TIMDFF + SECOND() - TIMD
 
C        SAVE DIAGNOSTIC FIELDS.
         IF (DIAG) THEN
            WRITE (LUNDIAG2/10) TIME
            WRITE (LUNDIAG2/10) DD
         END IF
 460  CONTINUE
 
      TDYN = SECOND()
C*C*C*
C     CALL SAVERPB(CONST,CZ,CZZ,FD0,TIME,IREAD,HNO3ZM)
      CALL SAVERP(CONST,CZ,CZZ,FD0,TIME,IREAD)
      IF (DIAG) CALL SVDIAG(EWFLUX,NSFLUX,VRFLUX,CHEMPROD,CHEMLOSS,
     1 PSFTERM,TIME,LUNDIAG1)
C     THE WIND,TEMPERATURE, + PRESSURE ARRAYS ARE REARRANGED
      DO 470 K=1,KTOP
      DO 470 J=1,JGMAX
      DO 470 I=1,N
         V1(I,J,K) = V2(I,J,K)
 470  CONTINUE
 
      DO 480 K=1,KTOP
      DO 480 J=0,JGMAX
      DO 480 I=1,N
         U1  (I,J,K) = U2  (I,J,K)
         KEL1(I,J,K) = KEL2(I,J,K)
 480  CONTINUE
 
      DO 485 J=0,JGMAX
      DO 485 I=1,N
         PS1(I,J) = PS2(I,J)
 485  CONTINUE
      TWIND1 = TWIND2
      TIMDYN = TIMDYN + SECOND() - TDYN
 500  CONTINUE
 520  CONTINUE
C
      TDIA = SECOND()
C     MASS OF CONSTITUENT IS DETERMINED.
      WRITE (6,314)  TIME/86400.
      CALL MASS3(CONST,PS1,MCOR,6)
 
      DO 530 K=1,KTOP
         WRITE (6,*) ' K = ', K
         WRITE (6,531) (CONST(1,J,K,1),J=0,JGMAX,4)
 530  CONTINUE
 531  FORMAT(12E10.3)
C      DO 540 IC = 1, NCON
C        THE CONSTITUENT FIELD AT LOCAL NOON IS PLOTTED OUT
C         CALL PLNON2(CONST(1,0,1,IC),PS1,IC,TIME,TDT,PLON,LNOON,LMID,
C     !    FILFOR,DIFF)
C  540 CONTINUE
C
      ICOUNT = ICOUNT + 1
      LUNOUT = 10 + ICOUNT/NPLOTS
      ILON = 0
      CALL HEADER(TIME,TDT,LUNOUT,ILON)
      DO 560 IC= 1,NCON
         CALL PL3D (CONST(1,0,1,IC),PS1,IC,LUNOUT)
 560  CONTINUE
      TIMDIA = TIMDIA + SECOND() - TDIA
C
C     DO 620 IC=1,NCON
C     DO 620 K=1,KTOP
C        IMAX=ISMAX(NJG2,CONST(1,0,K,IC),1)
C        IMIN=ISMIN(NJG2,CONST(1,0,K,IC),1)
C        ICNT= (IC-1)*NJG2*KTOP + (K-1)*NJG2
C        I22= (IMAX-1) / N
C        I33= (IMIN-1) / N
C        J2 = IMAX - I22 * N
C        J3 = IMIN - I33 * N
C        WRITE (6,618) J2,I22,K,CONST(ICNT+IMAX)
C        WRITE (6,619) J3,I33,K,CONST(ICNT+IMIN)
C618     FORMAT(' LONG = ',I4,' LAT = ',I4,' HT = ',I4,
C    1   ' MAX = ',1PE15.7)
C619     FORMAT(' LONG = ',I4,' LAT = ',I4,' HT = ',I4,
C    1   ' MIN = ',1PE15.7)
C620  CONTINUE
 
      DO 620 IC=1,NCON
      DO 620 K=1,KTOP
         DMAX = -1.E15
         DMIN =  1.E15
         DO 605 J=0,JGMAX
         DO 605 I=1,N
            IF (CONST(I,J,K,IC).GT.DMAX) THEN
               DMAX  = CONST(I,J,K,IC)
               JMAX  = J
               IMAX  = I
            ELSE IF (CONST(I,J,K,IC).LT.DMIN) THEN
               DMIN  = CONST(I,J,K,IC)
               JMIN  = J
               IMIN  = I
            END IF
 605     CONTINUE
         WRITE (6,618) IMAX,JMAX,K,DMAX
         WRITE (6,619) IMIN,JMIN,K,DMIN
 618     FORMAT(' LONG = ',I4,' LAT = ',I4,' HT = ',I4,
     1   ' MAX = ',1PE15.7)
 619     FORMAT(' LONG = ',I4,' LAT = ',I4,' HT = ',I4,
     1   ' MIN = ',1PE15.7)
 620  CONTINUE
 
      TIMIN = TIMIN + SECOND() - TIN
      TIMCHM = TIMCHM*100. / TIMIN
      TIMDYN = TIMDYN*100. / TIMIN
      TIMDIA = TIMDIA*100. / TIMIN
      TIMX =   TIMX  *100. / TIMIN
      TIMY =   TIMY  *100. / TIMIN
      TIMZ =   TIMZ  *100. / TIMIN
      TIMDFF = TIMDFF *100./ TIMIN
      WRITE (6,*) ' TOTAL TIME OF JOB ',TIMIN
      WRITE (6,211) TIMCHM
      WRITE (6,212) TIMDYN
      WRITE (6,214) TIMDIA
      WRITE (6,215) TIMX
      WRITE (6,216) TIMY
      WRITE (6,217) TIMZ
      WRITE (6,219) TIMDFF
 211  FORMAT(1X,' CHEMISTRY TIME',F7.2)
 212  FORMAT(1X,' DYNAMICAL TIME',F7.2)
 214  FORMAT(1X,' DIAGNOSTIC TIME',F7.2)
 215  FORMAT(1X,' X ADV     TIME',F7.2)
 216  FORMAT(1X,' Y ADV     TIME',F7.2)
 217  FORMAT(1X,' Z ADV     TIME',F7.2)
 219  FORMAT(1X,' DIFFUSION TIME',F7.2)
 
      STOP
      END
