Received: from SGCCP.GSFC.NASA.GOV by spica.npac.syr.edu (4.1/I-1.98K)
	id AA02931; Tue, 2 Feb 93 08:55:29 EST
Date:    Tue, 2 Feb 1993 8:56:14 EST
From: ALLEN@SGCCP.GSFC.NASA.GOV
Message-Id: <930202085614.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     TRACER CONTINUITY EQN IS SOLVED WITH VANLEER'S SCHEME IN A
C            SIGMA-PRESSURE COORDINATE SYSTEM.
C
C     DEVECTORIZED VERSION (30 JUNE 1990)
C     CRAY VERSION (1 AUG 1990)  DPSF READ IN (5 MAR 91)
C     LATEST UPDATE 21 JAN 91 
C
      PARAMETER(N=144,JGMAX=88,NCON=1,JGMAX1=JGMAX+1)
      PARAMETER (KSTRAT=11,KTROP=8,KTOP=KSTRAT+KTROP,NTAPES=61)
      PARAMETER (ILMM=N,IJMM=JGMAX)
C     FOR COMPATABILITY ADD 
      PARAMETER(JDIAG=5,JDIAG2=3,KDIAG2=3,NVARS=10) 
C
C     PHYSICAL SPACE ARRAYS ARE DEFINED
      REAL CONST(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON)
 
      REAL F1(N,0:JGMAX,KTOP)
      REAL F3(N,0:JGMAX,KTOP),SCALE(NCON)
C
      LOGICAL RSTAR,FIRST,WNDSAV,CHEM,DIFF,UPSTREAM
      INTEGER NSTEPS(NTAPES),IST1(NTAPES),NOBS7(NTAPES),NOBS8(NTAPES)
      INTEGER LUNCHM(NCON),GRID,WINDSOURCE 
      INTEGER ILUN(NTAPES),NOBS9(NTAPES) 
C
      REAL PSU(N,0:JGMAX),PSV(N,JGMAX),PSNEW(N,0:JGMAX)
      REAL DPSF(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,KTROP),SGSTRAT(N,0:JGMAX,0:KSTRAT-1),
     ! KEL1(N,0:JGMAX,KTOP), KEL2(N,0:JGMAX,KTOP),PS1INV(N,0:JGMAX)
C
      REAL COSCEN(0:JGMAX), DLATR(0:JGMAX),
     ! TWOLAT(0:JGMAX), MCOR(N,0:JGMAX), PLON(N)
      REAL DLATR2(JGMAX),COSEDGE(JGMAX)

C     FOR COMPATABILITY WITH MAINPRA AND MAIN2 ADD
      LOGICAL IDARK(N,0:JGMAX), DIURNAL,DIAG 
      INTEGER NDIAG(JDIAG)
      REAL CHEMPROD(N,JDIAG,KTOP,NCON),CHEMLOSS(N,JDIAG,KTOP,NCON)        
 
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 THETA1(0:JGMAX,KTOP), THETA2(0:JGMAX,KTOP)
      REAL THETA (0:JGMAX,KTOP) 

      REAL CL(0:JGMAX,KTOP),CLRAT(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/
      DATA NOBS9/ 20, 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, 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/ 
C
      DATA POLEDGE / -87./

      DATA NDIAG/44,49,59,74,84/

      DATA NNLAT/ 44,59,74/

      DATA NNLEV/ 9,13,15/ 
C
      TIN    = SECOND()
C
      WRITE (6,*) 'BEGIN EXECUTION'
C     INITIALIZE TIMING PARAMETERS
      TIMIN  = 0.
      TIMCHM = 0.
      TIMDYN = 0.
      TIMDIA = 0.
      TIMX   = 0.
      TIMXH  = 0.
      TIMY   = 0.
      TIMZ   = 0.
      TIMDFF  = 0.
C     DEFINE VARIABLES AND OUTPUT EXPERIMENT HEADER TO UNIT10
      CALL SET19
C
      WRITE (6,*) ' RESTART = ',RSTAR,'FIRST=',FIRST,' ISAVE=',ISAVE
      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
*********************SPECIAL FIX********
      ADPHI      = PI * A / 90.
****************************************
 
      ICHEM     = -999
 
C     LOGICAL UNIT NUMBERS ARE SPECIFIED.
      DO 22 IC=1,NCON
         LUNCHM(IC) = 48 + IC
 22   CONTINUE
      LUNCHM2 = 47 
      LNOON     =  33
      LMID      =  44
      WNDSAV    =  .FALSE.
 
      NDT       = 5
      TDT       =  DT / NTSTRC
      TDT2      =  TDT * 0.5
      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     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
      DPHIPOLE = 3. * PI / 180.
      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,*) DLATR2(J), DLATR(J), COSEDGE(J), COSCEN(J)
 100  CONTINUE
      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 RESTAR(CONST,FD0,TIME,IREAD,ISAVE)
         WRITE (6,*) ' IREAD = ',IREAD
 
C        ONE OBS IS READ BOTH IN THE 1ST RUN + IN THE SUBSEQUENT RESTART
         IREAD = IREAD - 1
         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
C***     WILL SKIP OVER DATA UNTIL 00Z DEC 28 ***
C        THE WIND, PRESSURE, AND TEMPERATURE FIELDS ARE READ IN
         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
         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)
c          DO 9009 K=1,KTOP
c          DO 9009 J=0,JGMAX
c          DO 9009 I=1,N
c             CONST(I,J,K,IC) = 10.
c             FD0(I,J,IC) = 10.
c 9009     CONTINUE
         CALL INITO3(CONST(1,0,1,IC),LUNHNO,SCALE(IC),FD0(1,0,IC))  
*****************************************************
C2299    FORMAT(10E12.5)
C        LUNHNO = 45
C        IC = 2
C         CALL TINIT(CONST(1,0,1,IC),LUNHNO,SCALE(IC),FD0(1,0,IC),IC)
C         CALL BLOB(CONST(1,0,1,1))
C         NEED 3-D FD0.
C           FD0(1) = 10.
C          FD0(3) = FD0(1)
C        INITIAL CONDITIONS FOR 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'
            STOP
         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
         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
             DAY0 = -4.
             IF (WINDSOURCE.EQ.1) DAY0 = 0.
C            CHEM SKIP PARAMETER FOR HCL EXPS.
             ICHEM2 = INT((DAY0 + TIME/86400.)/15.)
C            CHEM SKIP PARAMETER FOR O3,NOY, AND N2O EXPS.
C            WARNING INPUT O3 FILE DOES NOT CONTAIN HEADER!!!!
C            ICHEM2 = INT((DAY0+TIME/86400.)/15.)
             IF (ICHEM2.NE.ICHEM) THEN
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)
                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) 
                ICHEM = ICHEM2
             END IF
          END IF
         TIMCHM = TIMCHM + SECOND() - TCHM
         DO 460 ISTCO = 1,NTSTRC  
 
C        SPECIAL SECTION FOR O3 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) 
         TIMCHM = TIMCHM + SECOND() - TCHM 
C    
         TDYN = SECOND()
         TIME2 = TIME + TDT
C        INTERPOLATE TO FIND DYNAMIC FIELDS AT CURRENT TIME.
         CALL WDINT4 (U1,U2,V1,V2,
     !     PS1,PS2,KEL1,KEL2,TIME,TWIND1,TWIND2,PSNEW,TIME2)
         DO 101 J=0,JGMAX
         DO 101 I=1,N
            DPSF(I,J) = PSNEW(I,J) - PS1(I,J)
 101     CONTINUE
 
         TWIND1 = TIME
         TIME   = TIME + TDT
         CALL SURFPS(PS1,PSU,PSV)
 
         DO 102 IJ=1,IJMM
         DO 102 IL=1,ILMM
            PSV(IL,IJ) = PSV(IL,IJ) - PINT
 102     CONTINUE
 
         DO 105 IJ=0,IJMM
         DO 105 IL=1,ILMM
            PSU(IL,IJ) = PSU(IL,IJ) - PINT
            PS1INV(IL,IJ) = 1. / (PS1(IL,IJ)-PINT)
 105     CONTINUE
C        VERTICAL WIND IS CALCULATED BY DETERMINING THE X AND Y
         TX = SECOND()
C        ADVECTION TERMS AND THE PSF CHANGE TERM ON A CONSTANT FIELD
C        1. CALCULATE X ADV TERM ON CONSTANT FIELD AND PUT IN F1.
         DO 110 K=1,KTOP
            DO 107 IJ=0,IJMM
            DO 107 IL=1,ILMM
               F3(IL,IJ,K) = 1.0
 107        CONTINUE
 110     CONTINUE
C        CHANGE DUE TO X ADV AT SOUTHERN LATS IS RETURNED IN F1
         L1 = 1
         L2 = 8
 2004    FORMAT(10E12.5)
         UPSTREAM = .FALSE.
         CALL XADVHLAT(F3,COSCEN,TDT2,U1,F1,PS1INV,PSU,L1,L2,NDT,
     1    UPSTREAM,.FALSE.)
         TIMXH = TIMXH + SECOND() - TX
 
C        CHANGE DUE TO X ADV FOR MIDDLE LATS IS RETURNED IN F1.
         L1 = 9
         L2 = 79
         UPSTREAM = .FALSE.
         CALL XADV(F3,COSCEN,TDT2,U1,F1,PS1INV,PSU,L1,L2,UPSTREAM,
     1    .FALSE.)
C        CHANGE DUE TO X ADV AT NORTHERN LATS RET IN F1.
         TXH = SECOND()
         L1 = 80
         L2 = 87
         UPSTREAM = .FALSE.
         CALL XADVHLAT(F3,COSCEN,TDT2,U1,F1,PS1INV,PSU,L1,L2,NDT,
     1    UPSTREAM,.FALSE.)
         TIMX = TIMX + SECOND() - TX
         TIMXH = TIMXH + SECOND() - TXH
 2002    FORMAT(10(1PE12.5))
 
         TY = SECOND()
C        2. CALCULATE Y ADVECTION TERM AND PLACE IN F3.
         CALL YADV2(COSCEN,COSEDGE,V1,ADPHI,F3,PS1INV,PSV,TDT2)
         TIMY = TIMY + SECOND() - TY
 
C        SUM X AND Y ADVECTION TERMS
         DO 120 K=1,KTOP
            DO 115 IJ=0,IJMM
            DO 115 IL=1,ILMM
               F1(IL,IJ,K) = -(F1(IL,IJ,K) + F3(IL,IJ,K))
 115        CONTINUE
 120     CONTINUE
 
         DO 125 K=1,KTROP
         DO 125 J=0,JGMAX
         DO 125 I=1,N
            F1(I,J,K) = F1(I,J,K) + DPSF(I,J) / PS1(I,J)
 125     CONTINUE
 
         TZ = SECOND()
C        CALCULATE VERTICAL VELOCITY AS A RESIDUAL.
         CALL SGDOT3(SGTROP,SGSTRAT,SGBOT,F1,PS1,PS1INV,COSCEN,TDT2)
         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
             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
             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
             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)
            CALL MASS3(CONST,PS1,MCOR,6)
         END IF
 
C        THE CONSTITUENT FIELDS ARE SAVED.
         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 260 IC=1,NCON
               CALL PL3D (CONST(1,0,1,IC),PS1,IC,LUNOUT)
 260        CONTINUE
         ELSE IF (ISTCO.EQ.1) THEN
            FIRST=.TRUE.
         END IF
C
C        THE LONG NEAREST LOCAL NOON + MIDNIGHT IS PLOTTED EACH HR.
**         IF (MOD(TIME,3600.).EQ.0.0) THEN
**            DO 280 IC=1,NCON
**               CALL PLNON2(CONST(1,0,1,IC),PS1,IC,TWIND1,TDT,PLON,
**   !          LNOON,LMID)
**280         CONTINUE
**         END IF
         TIMDIA = TIMDIA + SECOND() - TDIA
 
         TDYN = SECOND()
         DO 350 ILOOP = 1,2
         DO 340 IC=1,NCON
C           CHANGE DUE TO XADV IS DETERMINED AND PUT IN F1.
            TX = SECOND()
            L1 = 1
            L2 = 8
            NLATS=L2-L1+1
            UPSTREAM = .FALSE.
            CALL XADVHLAT(CONST(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       L1,L2,NDT,UPSTREAM,.FALSE.)
            TIMXH = TIMXH + SECOND() - TX
            L1 = 9
            L2 = 79
            NLATS=L2-L1+1
            UPSTREAM = .FALSE.
           CALL XADV(CONST(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1       L1,L2,UPSTREAM,.FALSE.)
            TXH = SECOND()
            L1 = 80
            L2 = 87
            NLATS=L2-L1+1
            UPSTREAM = .FALSE.
            CALL XADVHLAT(CONST(1,0,1,IC),COSCEN,TDT2,U1,F1,PS1INV,PSU,
     1                                 L1,L2,NDT,UPSTREAM,.FALSE.)
            TIMX = TIMX + SECOND() - TX
            TIMXH = TIMXH + SECOND() - TXH
 
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,.FALSE.)
 
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 SIGADV(CONST(1,0,1,IC),TDT2,SGTROP,SGSTRAT,SGBOT,F3,
     1      FD0(1,0,IC))
            TIMZ = TIMZ + SECOND() - TZ
 
C          ADD IN ADVECTION TERMS AND ADJUST FOR SURFACE PRESSURE CHANGE
           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
 
 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,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) 
         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)
c        WRITE (6,*) ' DIFFUSION X '
c        WRITE (6,2004) (F1(5,J,8),J=0,JGMAX)
C        WRITE (6,2004) (F1(5,J,9),J=0,JGMAX)
               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)
c        WRITE (6,*) ' DIFFUSION Y '
c        WRITE (6,2004) (F1(5,J,8),J=0,JGMAX)
C        WRITE (6,2004) (F1(5,J,9),J=0,JGMAX)
               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)
c        WRITE (6,*) ' DIFFUSION X '
c        WRITE (6,2004) (F1(5,J,8),J=0,JGMAX)
C        WRITE (6,2004) (F1(5,J,9),J=0,JGMAX)
               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
 460  CONTINUE
 
      TDYN = SECOND()
      CALL SAVER (CONST,FD0,TIME,IREAD)
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)
      WRITE (6,*) ' FINAL MIXING RATIOS AT TIME ', TIME/86400.
      DO 530 K=1,KTOP
         WRITE (6,*) ' K = ', K
         WRITE (6,531) (CONST(1,J,K,1),J=0,JGMAX,2)
 530  CONTINUE
 531  FORMAT(12E10.3)
 
*     DO 540 IC = 1, NCON
C        THE CONSTITUENT FIELD AT LOCAL NOON IS PLOTTED OUT
*        CALL PLNON2(CONST(1,0,1,IC),PS1,IC,TIME,TDT,PLON,LNOON,LMID,
*    !    FILFOR,DIFF)
* 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
      WRITE (6,618)
 618  FORMAT('  K',T12,'MAXIMUM',T32,'MINIMUM')
      do 620 k=1,ktop
      amax=0.
      amin = 1000.
      do 600 j=0,jgmax
      do 600 i=1,n
         amax = max(const(i,j,k,1),amax)
         AMIN = MIN(CONST(I,J,K,1),AMIN)
 600  continue
      WRITE (6,619) K,AMAX,AMIN
 619  FORMAT(I3,2E20.10)
 620  continue
      TIMIN = TIMIN + SECOND() - TIN
      TIMCHM = TIMCHM*100. / TIMIN
      TIMDYN = TIMDYN*100. / TIMIN
      TIMDIA = TIMDIA*100. / TIMIN
      TIMX =   TIMX  *100. / TIMIN
      TIMXH = TIMXH  * 100. / TIMIN
      TIMY =   TIMY  *100. / TIMIN
      TIMZ =   TIMZ  *100. / TIMIN
      TIMDPSF= TIMDPSF*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,209) TIMXH
      WRITE (6,216) TIMY
      WRITE (6,217) TIMZ
      WRITE (6,218) TIMDPSF
      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)
 209  FORMAT(1X,' X ADV HIGH LATITUDES TIME',F7.2)
 216  FORMAT(1X,' Y ADV     TIME',F7.2)
 217  FORMAT(1X,' Z ADV     TIME',F7.2)
 218  FORMAT(1X,' DPSF      TIME',F7.2)
 219  FORMAT(1X,' DIFFUSION TIME',F7.2)
 
      STOP
      END
