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



      SUBROUTINE ZADV(C2,CZ1,CZZ1,FD0,SGSTRAT,SGTROP,SGBOT,DT,F3,CZ,CZZ,
     1 WTC,DIAG,PS1)
 
C     VERTICAL ADVECTION IS PERFORMED WITH PRATHER'S SCHEME.
C     FIRST VERSION 10/5/90   COMPLETED? 10/25/90
 
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11)
      PARAMETER (KST1=KSTRAT-1,KTOPM1=KTOP-1)
      PARAMETER (JDIAG2=3,KDIAG2=3,NVARS=10,KTROP2=8)
 
      REAL C2(N,0:JGMAX,KTOP),CZ(N,0:JGMAX,KTOP),CZZ(N,0:JGMAX,KTOP)
      REAL SGSTRAT(N,0:JGMAX,0:KST1), SGTROP(N,0:JGMAX,KTROP2)
      REAL SGBOT(N,0:JGMAX),F3(N,0:JGMAX,KTOP),WTC(N,0:JGMAX,0:KTOP)
      REAL CZ1(N,0:JGMAX,KTOP),CZZ1(N,0:JGMAX,KTOP),FD0(N,0:JGMAX)
      REAL PS1(N,0:JGMAX)
      LOGICAL DIAG
 
      COMMON/VERT/DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),SIGMA(KTOP),
     1 PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)
 
      COMMON / DIAGNOSTIC / DD(N,JDIAG2,KDIAG2,NVARS), NNLAT(JDIAG2),
     1 NNLEV(KDIAG2)
 
      COMMON //F0(N,0:JGMAX,0:KTOPM1),FZ(N,0:JGMAX,0:KTOPM1),
     1 FZZ(N,0:JGMAX,0:KTOPM1),C(N,0:JGMAX,KTOP),FD0WTC(N,0:JGMAX),
     1 F0CON(N,0:JGMAX,0:KTOPM1),CON(N,0:JGMAX,KTOP)
C
C     ORIGINAL VALUES OF MEAN, 1ST, AND 2ND MOMENTS ARE SAVED.
C     THESE VALUES WILL BE USED WHEN THE HORIZONTAL ADVECTION
C     TERMS ARE CALCULATED. BACKGROUND MIXING RATIO IS SET TO 1.
      DO 10 K=1,KTOP
      DO 10 J=0,JGMAX
      DO 10 I=1,N
         C  (I,J,K)   = C2  (I,J,K) * WTC(I,J,K)
         CZ (I,J,K)   = CZ1 (I,J,K)
         CZZ(I,J,K)   = CZZ1(I,J,K)
         CON(I,J,K)   = WTC(I,J,K)
 10   CONTINUE
 
C     WRITE (6,*) ' C IN Z ADV '
C     WRITE (6,4093) (C(I,JGMAX,14),I=1,N)
 4093 FORMAT (10(1PE12.5))
 
      DO 15 J=0,JGMAX
      DO 15 I=1,N
         FD0WTC(I,J) = FD0(I,J) * WTC(I,J,0)
 15   CONTINUE
 
C     LIMIT 1ST AND 2ND MOMENTS TO ASSURE THAT UPDATED MEAN FD IS
C     GREATER THAN OR EQUAL TO ZERO AFTER VERTICAL TRANSPORT.
C     SEE EQN. (39) OF PRATHER: NUMERICAL ADVECTION
      DO 20 K=1,KTOP
      DO 20 J=0,JGMAX
      DO 20 I=1,N
         SLPMAX = 0.
         IF (C(I,J,K).GT.0.) SLPMAX = C(I,J,K)
         S1MAX = 1.5* SLPMAX
         CZ(I,J,K)=MIN(S1MAX,MAX(-S1MAX,CZ(I,J,K)))
         CZZ(I,J,K) = MIN((SLPMAX+SLPMAX-0.333*ABS(CZ(I,J,K))),
     1                 MAX(ABS(CZ(I,J,K))-SLPMAX,CZZ(I,J,K)))
 20   CONTINUE
C
C     ZERO THE TEMPORARY MOMENTS.
      DO 40 K=0,KTOPM1
      DO 40 J=0,JGMAX
      DO 40 I=1,N
         F0(I,J,K) = 0.
         FZ(I,J,K) = 0.
         FZZ(I,J,K) = 0.
         F0CON(I,J,K) = 0.
 40   CONTINUE
 
C     BEGIN Z ADVECTION FOR LEVEL 1.
      K=0
      KPLUS = K + 1
      DO 60 J=0,JGMAX
      DO 60 I=1,N
         IF (SGBOT(I,J).LT.0.) THEN
            EPS = ABS(SGBOT(I,J)*DT/DELTAS(1))
            EPSQ = EPS * EPS
            EPS1 = 1.0 - EPS
            EPS1Q = EPS1 * EPS1
C           DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
            F0S = EPS* FD0WTC(I,J)
            FZS = 0.
            FZZS = 0.
C           DETERMINE BACKGROUND MIXING RATIO CHANGE
            TEMP = EPS*CON(I,J,1)
            F0CON(I,J,K) = F0CON(I,J,K) + TEMP
C           MOMENTS FOR LEVEL ZERO ARE NOT ADJUSTED.
C           SET TEMPORARY MOMENTS
            F0(I,J,K) = F0(I,J,K) + F0S
            FZ(I,J,K) = FZ(I,J,K) + FZS
            FZZ(I,J,K) = FZZ(I,J,K) + FZZS
         ELSE
            EPS = ABS(SGBOT(I,J)*DT/DELTAS(1))
            EPSQ = EPS*EPS
            EPS1 = 1.0 - EPS
            EPS1Q = EPS1 * EPS1
C           DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
            F0S = EPS * (C(I,J,KPLUS)-EPS1*(CZ(I,J,KPLUS)-
     1       (EPS1-EPS)*CZZ(I,J,KPLUS)))
            FZS=EPSQ*(CZ(I,J,KPLUS)-3.0*EPS1*CZZ(I,J,KPLUS))
            FZZS=EPS*EPSQ*CZZ(I,J,KPLUS)
C           DETERMINE BACKGROUND MIXING RATIO CHANGE
            TEMP = EPS*CON(I,J,KPLUS)
            F0CON(I,J,K) = F0CON(I,J,K) + TEMP
            CON(I,J,KPLUS) = CON(I,J,KPLUS) - TEMP
C           READJUST REMAINING MOMENTS.
            C(I,J,KPLUS) = C(I,J,KPLUS) - F0S
            CZ(I,J,KPLUS) = EPS1Q*(CZ(I,J,KPLUS)+3.0*EPS
     1       *CZZ(I,J,KPLUS))
            CZZ(I,J,KPLUS) = EPS1*EPS1Q*CZZ(I,J,KPLUS)
C           SET TEMPORARY MOMENTS
            F0(I,J,K) = F0(I,J,K) + F0S
            FZ(I,J,K) = FZ(I,J,K) + FZS
            FZZ(I,J,K) = FZZ(I,J,K) + FZZS
         END IF
 60   CONTINUE
C
C     BEGIN TROPOSPHERIC LOOP
      DO 100 K=1,KTROP-1
         KPLUS = K + 1
         DO 80  J=0,JGMAX
         DO 80  I=1,N
            IF (SGTROP(I,J,K).LT.0.) THEN
               EPS = ABS(SGTROP(I,J,K)*DT/DELTAS(K))
               EPSQ = EPS * EPS
               EPS1 = 1.0 - EPS
               EPS1Q = EPS1 * EPS1
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
               F0S = EPS* (C(I,J,K)+EPS1*(CZ(I,J,K)+
     1          (EPS1-EPS)*CZZ(I,J,K)))
               FZS = EPSQ*(CZ(I,J,K)+3.0*EPS1*CZZ(I,J,K))
               FZZS = EPS*EPSQ*CZZ(I,J,K)
C              DETERMINE BACKGROUND MIXING RATIO CHANGE
               TEMP = EPS*CON(I,J,K)
               F0CON(I,J,K) = F0CON(I,J,K) + TEMP
               CON(I,J,K) = CON(I,J,K) -TEMP
C              READJUST MOMENTS REMAINING IN BOX.
               C(I,J,K) = C(I,J,K) - F0S
               CZ(I,J,K)=EPS1Q*(CZ(I,J,K)-3.0*EPS*CZZ(I,J,K))
               CZZ(I,J,K)=EPS1*EPS1Q*CZZ(I,J,K)
C              SET TEMPORARY MOMENTS
               F0(I,J,K) = F0(I,J,K) + F0S
               FZ(I,J,K) = FZ(I,J,K) + FZS
               FZZ(I,J,K) = FZZ(I,J,K) + FZZS
            ELSE
               EPS = ABS(SGTROP(I,J,K)*DT/DELTAS(KPLUS))
               EPSQ = EPS*EPS
               EPS1 = 1.0 - EPS
               EPS1Q = EPS1 * EPS1
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
               F0S = EPS * (C(I,J,KPLUS)-EPS1*(CZ(I,J,KPLUS)-
     1          (EPS1-EPS)*CZZ(I,J,KPLUS)))
               FZS=EPSQ*(CZ(I,J,KPLUS)-3.0*EPS1*CZZ(I,J,KPLUS))
               FZZS=EPS*EPSQ*CZZ(I,J,KPLUS)
C              DETERMINE BACKGROUND MIXING RATIO CHANGE
               TEMP = EPS*CON(I,J,KPLUS)
               F0CON(I,J,K) = F0CON(I,J,K) + TEMP
               CON(I,J,KPLUS) = CON(I,J,KPLUS) - TEMP
C              READJUST REMAINING MOMENTS.
               C(I,J,KPLUS) = C(I,J,KPLUS) - F0S
               CZ(I,J,KPLUS) = EPS1Q*(CZ(I,J,KPLUS)+3.0*EPS
     1          *CZZ(I,J,KPLUS))
               CZZ(I,J,KPLUS) = EPS1*EPS1Q*CZZ(I,J,KPLUS)
C              SET TEMPORARY MOMENTS
               F0(I,J,K) = F0(I,J,K) + F0S
               FZ(I,J,K) = FZ(I,J,K) + FZS
               FZZ(I,J,K) = FZZ(I,J,K) + FZZS
            END IF
 80      CONTINUE
 100  CONTINUE
 
C     INTERFACE SECTION.
      IF ((KTROP.GT.0).AND.(KSTRAT.GT.0)) THEN
      K=KTROP
      K2 = K-KTROP
      KPLUS = K + 1
      DO 105  J=0,JGMAX
      DO 105  I=1,N
         IF (SGTROP(I,J,K).LT.0.) THEN
            EPS = ABS(SGTROP(I,J,K)*DT/DELTAS(K))
            EPSQ = EPS * EPS
            EPS1 = 1.0 - EPS
            EPS1Q = EPS1 * EPS1
C           DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
            F0S = EPS* (C(I,J,K)+EPS1*(CZ(I,J,K)+
     1       (EPS1-EPS)*CZZ(I,J,K)))
            FZS = EPSQ*(CZ(I,J,K)+3.0*EPS1*CZZ(I,J,K))
            FZZS = EPS*EPSQ*CZZ(I,J,K)
C           DETERMINE BACKGROUND MIXING RATIO CHANGE
            TEMP = EPS*CON(I,J,K)
            F0CON(I,J,K) = F0CON(I,J,K) + TEMP
            CON(I,J,K) = CON(I,J,K) -TEMP
C           READJUST MOMENTS REMAINING IN BOX.
            C(I,J,K) = C(I,J,K) - F0S
            CZ(I,J,K)=EPS1Q*(CZ(I,J,K)-3.0*EPS*CZZ(I,J,K))
            CZZ(I,J,K)=EPS1*EPS1Q*CZZ(I,J,K)
C           SET TEMPORARY MOMENTS
            F0(I,J,K) = F0(I,J,K) + F0S
            FZ(I,J,K) = FZ(I,J,K) + FZS
            FZZ(I,J,K) = FZZ(I,J,K) + FZZS
         ELSE
            EPS = ABS(SGSTRAT(I,J,K2)*DT/DELTAS(KPLUS))
            EPSQ = EPS*EPS
            EPS1 = 1.0 - EPS
            EPS1Q = EPS1 * EPS1
C           DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
            F0S = EPS * (C(I,J,KPLUS)-EPS1*(CZ(I,J,KPLUS)-
     1       (EPS1-EPS)*CZZ(I,J,KPLUS)))
C           IF (F0S.LT.0.) WRITE (6,*) ' WARNING ',I,J,K
            FZS=EPSQ*(CZ(I,J,KPLUS)-3.0*EPS1*CZZ(I,J,KPLUS))
            FZZS=EPS*EPSQ*CZZ(I,J,KPLUS)
C           DETERMINE BACKGROUND MIXING RATIO CHANGE
            TEMP = EPS*CON(I,J,KPLUS)
            F0CON(I,J,K) = F0CON(I,J,K) + TEMP
            CON(I,J,KPLUS) = CON(I,J,KPLUS) - TEMP
C           READJUST REMAINING MOMENTS.
            C(I,J,KPLUS) = C(I,J,KPLUS) - F0S
            CZ(I,J,KPLUS) = EPS1Q*(CZ(I,J,KPLUS)+3.0*EPS
     1       *CZZ(I,J,KPLUS))
            CZZ(I,J,KPLUS) = EPS1*EPS1Q*CZZ(I,J,KPLUS)
C           SET TEMPORARY MOMENTS
            F0(I,J,K) = F0(I,J,K) + F0S
            FZ(I,J,K) = FZ(I,J,K) + FZS
            FZZ(I,J,K) = FZZ(I,J,K) + FZZS
         END IF
 105  CONTINUE
      END IF
 
C     BEGIN STRATOSPHERIC LOOP
      DO 140 K=KTROP+1,KTOP-1
         KPLUS = K + 1
         K2 = K-KTROP
         DO 120 J=0,JGMAX
         DO 120 I=1,N
            IF (SGSTRAT(I,J,K2).LT.0.) THEN
               EPS = ABS(SGSTRAT(I,J,K2)*DT/DELTAS(K))
               EPSQ = EPS * EPS
               EPS1 = 1.0 - EPS
               EPS1Q = EPS1 * EPS1
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
               F0S = EPS* (C(I,J,K)+EPS1*(CZ(I,J,K)+
     1          (EPS1-EPS)*CZZ(I,J,K)))
               FZS = EPSQ*(CZ(I,J,K)+3.0*EPS1*CZZ(I,J,K))
               FZZS = EPS*EPSQ*CZZ(I,J,K)
C              DETERMINE BACKGROUND MIXING RATIO CHANGE
               TEMP = EPS*CON(I,J,K)
               F0CON(I,J,K) = F0CON(I,J,K) + TEMP
               CON(I,J,K) = CON(I,J,K) -TEMP
C              READJUST MOMENTS REMAINING IN BOX.
               C(I,J,K) = C(I,J,K) - F0S
               CZ(I,J,K)=EPS1Q*(CZ(I,J,K)-3.0*EPS*CZZ(I,J,K))
               CZZ(I,J,K)=EPS1*EPS1Q*CZZ(I,J,K)
C              SET TEMPORARY MOMENTS
               F0(I,J,K) = F0(I,J,K) + F0S
               FZ(I,J,K) = FZ(I,J,K) + FZS
               FZZ(I,J,K) = FZZ(I,J,K) + FZZS
    9          FORMAT(I3,6(1PE12.5),/,8(1PE12.5))
            ELSE
               EPS = ABS(SGSTRAT(I,J,K2)*DT/DELTAS(KPLUS))
               EPSQ = EPS*EPS
               EPS1 = 1.0 - EPS
               EPS1Q = EPS1 * EPS1
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.
               F0S = EPS * (C(I,J,KPLUS)-EPS1*(CZ(I,J,KPLUS)-
     1          (EPS1-EPS)*CZZ(I,J,KPLUS)))
               FZS=EPSQ*(CZ(I,J,KPLUS)-3.0*EPS1*CZZ(I,J,KPLUS))
               FZZS=EPS*EPSQ*CZZ(I,J,KPLUS)
C              DETERMINE BACKGROUND MIXING RATIO CHANGE
               TEMP = EPS*CON(I,J,KPLUS)
               F0CON(I,J,K) = F0CON(I,J,K) + TEMP
               CON(I,J,KPLUS) = CON(I,J,KPLUS) - TEMP
C              READJUST REMAINING MOMENTS.
               C(I,J,KPLUS) = C(I,J,KPLUS) - F0S
               CZ(I,J,KPLUS) = EPS1Q*(CZ(I,J,KPLUS)+3.0*EPS
     1          *CZZ(I,J,KPLUS))
               CZZ(I,J,KPLUS) = EPS1*EPS1Q*CZZ(I,J,KPLUS)
C              SET TEMPORARY MOMENTS
               F0(I,J,K) = F0(I,J,K) + F0S
               FZ(I,J,K) = FZ(I,J,K) + FZS
               FZZ(I,J,K) = FZZ(I,J,K) + FZZS
            END IF
 120     CONTINUE
 140  CONTINUE
 
C     BEGIN TEDIOUS FLUX CALCULATION.  FOR CONSISTENCY THE MIXING
C     RATIO FLUX IS ESTIMATED.
C     WRITE (6,*) ' DIAG = ', DIAG,KDIAG2,JDIAG2
      IF (DIAG) THEN
         DO 158 KK=1,KDIAG2
            K2 = NNLEV(KK)
            IF (K2.EQ.0) THEN
               DO 145 JJ=1,JDIAG2
                  DO 145 I=1,N
                  IF (SGBOT(I,NNLAT(JJ)).LT.0.) THEN
                     DD(I,JJ,KK,10)=F0(I,NNLAT(JJ),K2)
                  ELSE
                     DD(I,JJ,KK,10)=-F0(I,NNLAT(JJ),K2)
                  END IF
 145           CONTINUE
            ELSE IF (K2.LE.KTROP) THEN
               DO 148 JJ=1,JDIAG2
               DO 148 I=1,N
                  IF (SGTROP(I,NNLAT(JJ),K2).LT.0.) THEN
                     DD(I,JJ,KK,10)=F0(I,NNLAT(JJ),K2)
                  ELSE
                     DD(I,JJ,KK,10)=-F0(I,NNLAT(JJ),K2)
                  END IF
 148           CONTINUE
            ELSE IF (K2.GT.KTROP) THEN
               DO 152 JJ=1,JDIAG2
               DO 152 I=1,N
                  IF (SGSTRAT(I,NNLAT(JJ),K2).LT.0.) THEN
                     DD(I,JJ,KK,10)=F0(I,NNLAT(JJ),K2)
                  ELSE
                     DD(I,JJ,KK,10)=-F0(I,NNLAT(JJ),K2)
                  END IF
 152           CONTINUE
            END IF
 158     CONTINUE
      END IF
 
C     WRITE (6,*)
C     RECONSTRUCT LOWEST LAYER.
      K=0
C     WRITE (6,*) ' BEGIN RECONSTRUCT FOR LEVEL ',K
      KPLUS = K + 1
      DO 160 J=0,JGMAX
      DO 160 I=1,N
         IF (SGBOT(I,J).LT.0.) THEN
            EPS = F0CON(I,J,K) / (F0CON(I,J,K) + CON(I,J,KPLUS))
            EPS1 = 1.-EPS
            TM=EPS*C(I,J,KPLUS)-EPS1*F0(I,J,K)
            C(I,J,KPLUS) = C(I,J,KPLUS) + F0(I,J,K)
            CZZ(I,J,KPLUS) = EPS*EPS*FZZ(I,J,K) +
     1       EPS1*EPS1*CZZ(I,J,KPLUS)+5.0*(EPS*EPS1
     1       *(CZ(I,J,KPLUS)-FZ(I,J,K))-(EPS1-EPS)*TM)
            CZ(I,J,KPLUS) = EPS*FZ(I,J,K)+EPS1*CZ(I,J,KPLUS)
     1       + 3.0*TM
            CON(I,J,KPLUS) = CON(I,J,KPLUS) + F0CON(I,J,K)
         END IF
 160  CONTINUE
 
 
C     RECONSTRUCT TROPOSPHERIC BOXES
      DO 500 K=1,KTROP-1
         KPLUS = K + 1
         DO 480 J=0,JGMAX
         DO 480 I=1,N
            IF (SGTROP(I,J,K).LT.0.) THEN
               EPS = F0CON(I,J,K) / (F0CON(I,J,K) + CON(I,J,KPLUS))
               EPS1 = 1.-EPS
               TM=EPS*C(I,J,KPLUS)-EPS1*F0(I,J,K)
               C(I,J,KPLUS) = C(I,J,KPLUS) + F0(I,J,K)
               CZZ(I,J,KPLUS)=EPS*EPS*FZZ(I,J,K)+EPS1*
     1          *EPS1*CZZ(I,J,KPLUS)+5.0*(EPS*EPS1
     1          *(CZ(I,J,KPLUS)-FZ(I,J,K))-(EPS1-EPS)*TM)
               CZ(I,J,KPLUS) = EPS*FZ(I,J,K)+EPS1
     1          *CZ(I,J,KPLUS)+ 3.0*TM
               CON(I,J,KPLUS) = CON(I,J,KPLUS) + F0CON(I,J,K)
            ELSE
               EPS = F0CON(I,J,K) / (F0CON(I,J,K) + CON(I,J,K))
               EPS1 = 1. - EPS
               TM = -EPS*C(I,J,K) + EPS1*F0(I,J,K)
               C(I,J,K) = C(I,J,K) + F0(I,J,K)
               CZZ(I,J,K) = EPS*EPS*FZZ(I,J,K)+
     1          EPS1*EPS1*CZZ(I,J,K) + 5.*
     1          (EPS*EPS1*(-CZ(I,J,K)+FZ(I,J,K))+
     1          (EPS1-EPS)*TM)
               CON(I,J,K) = CON(I,J,K) + F0CON(I,J,K)
               CZ(I,J,K) = EPS*FZ(I,J,K)+EPS1*CZ(I,J,K)+3.*TM
            END IF
 480     CONTINUE
 500  CONTINUE
 
C     RECONSTRUCT IN TROPOSPHERE-STRATOSPHERE INTERFACE REGION.
      IF ((KTROP.GT.0).AND.(KSTRAT.GT.0)) THEN
      K=KTROP
      KPLUS = K + 1
      K2 = K-KTROP
      DO 210 J=0,JGMAX
      DO 210 I=1,N
         IF (SGTROP(I,J,K).LT.0.) THEN
            EPS = F0CON(I,J,K) / (F0CON(I,J,K) + CON(I,J,KPLUS))
            EPS1 = 1.-EPS
            TM=EPS*C(I,J,KPLUS)-EPS1*F0(I,J,K)
            C(I,J,KPLUS) = C(I,J,KPLUS) + F0(I,J,K)
            CZZ(I,J,KPLUS)=EPS*EPS*FZZ(I,J,K)+EPS1*
     1       *EPS1*CZZ(I,J,KPLUS)+5.0*(EPS*EPS1
     1       *(CZ(I,J,KPLUS)-FZ(I,J,K))-(EPS1-EPS)*TM)
            CZ(I,J,KPLUS) = EPS*FZ(I,J,K)+EPS1
     1       *CZ(I,J,KPLUS)+ 3.0*TM
            CON(I,J,KPLUS) = CON(I,J,KPLUS) + F0CON(I,J,K)
         ELSE
            EPS = F0CON(I,J,K) / (F0CON(I,J,K) + CON(I,J,K))
            EPS1 = 1. - EPS
            TM = -EPS*C(I,J,K) + EPS1*F0(I,J,K)
            C(I,J,K) = C(I,J,K) + F0(I,J,K)
            CZZ(I,J,K) = EPS*EPS*FZZ(I,J,K)+
     1       EPS1*EPS1*CZZ(I,J,K) + 5.*
     1       (EPS*EPS1*(-CZ(I,J,K)+FZ(I,J,K))+
     1       (EPS1-EPS)*TM)
            CZ(I,J,K) = EPS*FZ(I,J,K)+EPS1*CZ(I,J,K)+3.*TM
            CON(I,J,K) = CON(I,J,K) + F0CON(I,J,K)
         END IF
 210  CONTINUE
      END IF
 
C     RECONSTRUCT STRATOSPHERIC BOXES
      DO 240 K=KTROP+1,KTOPM1
         K2 = K-KTROP
         KPLUS = K + 1
         DO 220 J=0,JGMAX
         DO 220 I=1,N
            IF (SGSTRAT(I,J,K2).LT.0.) THEN
               EPS = F0CON(I,J,K) / (F0CON(I,J,K) + CON(I,J,KPLUS))
               EPS1 = 1.-EPS
               TM=EPS*C(I,J,KPLUS)-EPS1*F0(I,J,K)
               C(I,J,KPLUS) = C(I,J,KPLUS) + F0(I,J,K)
               CZZ(I,J,KPLUS)=EPS*EPS*FZZ(I,J,K)+EPS1*
     1          *EPS1*CZZ(I,J,KPLUS)+5.0*(EPS*EPS1
     1          *(CZ(I,J,KPLUS)-FZ(I,J,K))-(EPS1-EPS)*TM)
               CZ(I,J,KPLUS) = EPS*FZ(I,J,K)+EPS1
     1          *CZ(I,J,KPLUS)+ 3.0*TM
               CON(I,J,KPLUS) = CON(I,J,KPLUS) + F0CON(I,J,K)
 701          FORMAT(I4,9(1PE12.5))
            ELSE
               EPS = F0CON(I,J,K) / (F0CON(I,J,K) + CON(I,J,K))
               EPS1 = 1. - EPS
               TM = -EPS*C(I,J,K) + EPS1*F0(I,J,K)
               C(I,J,K) = C(I,J,K) + F0(I,J,K)
               CZZ(I,J,K) = EPS*EPS*FZZ(I,J,K)+
     1          EPS1*EPS1*CZZ(I,J,K) + 5.*
     1          (EPS*EPS1*(-CZ(I,J,K)+FZ(I,J,K))+
     1          (EPS1-EPS)*TM)
               CZ(I,J,K) = EPS*FZ(I,J,K)+EPS1*CZ(I,J,K)+3.*TM
               CON(I,J,K) = CON(I,J,K) + F0CON(I,J,K)
            END IF
 220     CONTINUE
 240  CONTINUE
 
      DO 260 K=1,KTOP
      DO 260 J=0,JGMAX
      DO 260 I=1,N
         F3(I,J,K) = C(I,J,K)/WTC(I,J,K)  - C2(I,J,K)
 260  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE SGPRA(F1,SGBOT,SGTROP,SGSTRAT,DT,PS1,WTC)
 
C     THE VERTICAL VELOCITY IS DETERMINED RESIDUALLY FOR USE WITH
C     PRATHER'S SCHEME.
C     OCTOBER 25, 1990
 
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11)
      PARAMETER (KTROP2=8)
 
C     F1  INPUT MIXING RATIO CHANGE NEEDED TO SATISFY CONTINUITY
C         OUTPUT MASS FLUX NEEDED TO SATISFY CONTINUITY
C
      REAL SGBOT(N,0:JGMAX), SGTROP(N,0:JGMAX,KTROP2),
     1 SGSTRAT(N,0:JGMAX,0:KSTRAT-1),PS1(N,0:JGMAX),
     1 F1(N,0:JGMAX,KTOP),  WTC(N,0:JGMAX,0:KTOP)
 
      COMMON / VERT/DELTAS(KTOP), DSIG(0:KTOP),DELZ(KSTRAT),SIGMA(KTOP),
     1 PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)
 
      COMMON // P2(N,0:JGMAX), P0(N,0:JGMAX), P1(N,0:JGMAX)
 
C     REQUIRED MIXING RATIO CHANGE IS CONVERTED TO TERM PROPORTIONAL TO
C     REQUIRED MASS CHANGE BY MULTIPLYING BY TERM PROPORTIONAL TO
C     DELZ * PRESSURE/ TEMPERATURE
      DO 10 K=1,KTROP
      DO 10 J=0,JGMAX
      DO 10 I=1,N
C        P2(I,J)  = PINT + SGINT(K  )*(PS1(I,J)-PINT)
C        P0(I,J)  = PINT + SGINT(K-1)*(PS1(I,J)-PINT)
C        P1(I,J)  = PINT + SIGMA(K)  *(PS1(I,J)-PINT)
         WTC(I,J,K) = (PS1(I,J)-PINT)*(SGINT(K-1)-SGINT(K))
         F1(I,J,K) = F1(I,J,K)*WTC(I,J,K)
 10   CONTINUE
 
C     DETERMINE DELZ IN THE STRATOSPHERE AND THEN VOL OF EACH BOX.
C     NOTE: -.543 = LN( P(TOP OF VOL)/P(BOTTOM OF VOL) ) IN STRAT
      DO 20 K=KTROP+1,KTOP
         K2 = K - KTROP
         DO 15 J=0,JGMAX
         DO 15 I=1,N
C           WTC(I,J,K) = 0.543 * PRESS (K2)
            WTC(I,J,K) = PRSLAY(K2-1)-PRSLAY(K2)
            F1(I,J,K) = F1(I,J,K) * WTC(I,J,K)
 15      CONTINUE
 20   CONTINUE
 
      DO 30 J=0,JGMAX
      DO 30 I=1,N
         WTC(I,J,0) = WTC(I,J,1)
 30   CONTINUE
 
      DT1 = 1. / DT
C     DETERMINE FLUXES REQUIRED TO KEEP CONSTANT FIELD CONSTANT.
      DO 40 K=KTOP,2,-1
         KM1 = K-1
         DO 35 J=0,JGMAX
         DO 35 I=1,N
            F1(I,J,KM1) = F1(I,J,KM1) + F1(I,J,K)
 35      CONTINUE
 40   CONTINUE
 
C     DETERMINE VELOCITIES AT SURFACE.
      K=1
      DO 50 J=0,JGMAX
      DO 50 I=1,N
         SGBOT(I,J) = F1(I,J,K) * DELTAS(1) * DT1 /
     1    WTC(I,J,K)
 50   CONTINUE
 
      IF (KTROP.GT.0) THEN
         K=2
         KM1 = K-1
         DO 60 J=0,JGMAX
         DO 60 I=1,N
            IF ((F1(I,J,K).GT.0.).AND.(F1(I,J,KM1).LT.0.)) THEN
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(KM1)*DT1/
     1         (WTC(I,J,KM1) * (1.-ABS(SGBOT(I,J)*DT/DELTAS(KM1))))
            ELSE IF (F1(I,J,K).GT.0.) THEN
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(KM1) * DT1 /
     1          WTC(I,J,KM1)
            ELSE
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(K) * DT1/
     1          WTC(I,J,K)
            END IF
 60      CONTINUE
      ELSE
         K=2
         KM1 = K-1
         DO 70 J=0,JGMAX
         DO 70 I=1,N
            IF ((F1(I,J,K).GT.0.).AND.(F1(I,J,KM1).LT.0.)) THEN
               SGSTRAT(I,J,KM1) = F1(I,J,K) * DELTAS(KM1)*DT1/
     1         (WTC(I,J,KM1) * (1.-ABS(SGBOT(I,J)*DT/DELTAS(KM1))))
            ELSE IF (F1(I,J,K).GT.0.) THEN
               SGSTRAT(I,J,KM1) = F1(I,J,K) * DELTAS(KM1) * DT1 /
     1          WTC(I,J,KM1)
            ELSE
               SGSTRAT(I,J,KM1) = F1(I,J,K) * DELTAS(K) * DT1/
     1          WTC(I,J,K)
            END IF
 70      CONTINUE
      END IF
 
      DO 100 K=3,KTROP
         KM1 = K-1
         KM2 = K-2
         DO 80 J=0,JGMAX
         DO 80 I=1,N
            IF ((F1(I,J,K).GT.0.).AND.(F1(I,J,KM1).LT.0.)) THEN
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(KM1)*DT1/
     1         (WTC(I,J,KM1)*(1.-ABS(SGTROP(I,J,KM2)*DT/DELTAS(KM1))))
            ELSE IF (F1(I,J,K).GT.0.) THEN
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(KM1) * DT1 /
     1          WTC(I,J,KM1)
            ELSE
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(K) * DT1 /
     1          WTC(I,J,K)
            END IF
 80      CONTINUE
 100  CONTINUE
 
      IF ((KTROP.GT.0).AND.(KSTRAT.GT.0)) THEN
         K=KTROP+1
         KS = K-KTROP
         KSM1 = KS - 1
         KM1 = K-1
         KM2 = K-2
         DO 120 J=0,JGMAX
         DO 120 I=1,N
            IF ((F1(I,J,K).GT.0.).AND.(F1(I,J,KM1).LT.0.)) THEN
               CONV = (PS1(I,J)-PINT) / (PINT-PTOP)
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(KM1)*DT1/
     1          (WTC(I,J,KM1)*(1.-ABS(SGTROP(I,J,KM2)*DT/DELTAS(KM1))))
               SGSTRAT(I,J,KSM1) = SGTROP(I,J,KM1) * CONV
            ELSE IF (F1(I,J,K).GT.0.) THEN
               CONV = (PS1(I,J)-PINT) / (PINT-PTOP)
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(KM1) * DT1/
     1          WTC(I,J,KM1)
               SGSTRAT(I,J,KSM1) = SGTROP(I,J,KM1) * CONV
            ELSE
               CONV = (PINT-PTOP) / (PS1(I,J)-PINT)
               SGSTRAT(I,J,KSM1) = F1(I,J,K) * DELTAS(K) * DT1/
     1          WTC(I,J,K)
               SGTROP(I,J,KM1) = SGSTRAT(I,J,KSM1) * CONV
            END IF
 120     CONTINUE
      ELSE
         DO 130 J=0,JGMAX
         DO 130 I=1,N
            SGSTRAT(I,J,0) = SGBOT(I,J)
 130     CONTINUE
      END IF
 
      DO 160 K=KTROP+2,KTOP
         KM1 = K-1
         KS = K-KTROP
         KSM1 = KS -1
         KSM2 = KS-2
         DO 140 J=0,JGMAX
         DO 140 I=1,N
            IF ((F1(I,J,K).GT.0.).AND.(F1(I,J,KM1).LT.0.)) THEN
               SGSTRAT(I,J,KSM1)=F1(I,J,K)*DELTAS(KM1)*DT1/
     1         (WTC(I,J,KM1)*(1.-ABS(SGSTRAT(I,J,KSM2)*DT/DELTAS(KM1))))
            ELSE IF (F1(I,J,K).GT.0.) THEN
               SGSTRAT(I,J,KSM1) = F1(I,J,K) * DELTAS(KM1) * DT1/
     1          WTC(I,J,KM1)
            ELSE
               SGSTRAT(I,J,KSM1) = F1(I,J,K) * DELTAS(K) * DT1/
     1          WTC(I,J,K)
            END IF
 140     CONTINUE
 160  CONTINUE
      RETURN
      END
