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



      SUBROUTINE XADVP(C,CX,CXX,CY,CYY,CZ,CZZ,CXY,CXZ,CYZ,MASS,MASS2,           
     1 DT,U,L1,L2,COSCEN,PSCINV2,PSUEDGE2)                                      
                                                                                
C     EAST-WEST ADVECTION IS PERFORMED WITH PRATHER'S SCHEME.                   
C     FIRST VERSION 1/5/91   MODIFIED FOR N1 3/27/91.                           
                                                                                
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11,N1=n+0)               
                                                                                
      REAL CX(N1,0:JGMAX,KTOP),CY(N1,0:JGMAX,KTOP),CZ(N1,0:JGMAX,KTOP)          
      REAL CXX(N1,0:JGMAX,KTOP),CYY(N1,0:JGMAX,KTOP)                            
      REAL CXY(N1,0:JGMAX,KTOP),CXZ(N1,0:JGMAX,KTOP)                            
      REAL CZZ(N1,0:JGMAX,KTOP),CYZ(N1,0:JGMAX,KTOP)                            
      REAL C(N1,0:JGMAX,KTOP),MASS2(N1,0:JGMAX,0:KTOP)                          
      REAL U(N1,0:JGMAX,KTOP),COSCEN(0:JGMAX),PSCINV2(N1,0:JGMAX)               
      REAL PSUEDGE2(N1,0:JGMAX),MASS(N1,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)                   
                                                                                
      REAL DTOVDXL(0:JGMAX),PSCINV(N1,0:JGMAX,KTOP),                            
     1 FX(N1,KTOP),FY(N1,KTOP),FZ(N1,KTOP),                                     
     1 FXX(N1,KTOP),FYY(N1,KTOP),FZZ(N1,KTOP),                                  
     2 FXY(N1,KTOP),FXZ(N1,KTOP),FYZ(N1,KTOP),                                  
     2 F0(N1,KTOP),F0CON(N1,KTOP),PSUEDGE(N1,0:JGMAX,KTOP)                      
C                                                                               
C     SET CONSTANTS                                                             
C     WRITE (6,*) ' ENTER XADV'                                                 
      PI = 4.0 * ATAN(1.)                                                       
      A=6.37E6                                                                  
      DIST=2.*PI*A/N                                                            
      DO 20 J=L1,L2                                                             
         DTOVDXL(J) = DT / (DIST*COSCEN(J))                                     
 20   CONTINUE                                                                  
                                                                                
      DO 30 K=1,KTROP                                                           
      DO 30 J=L1,L2                                                             
      DO 30 I=1,N                                                               
         PSCINV(I,J,K) = PSCINV2(I,J)                                           
         PSUEDGE(I,J,K)= PSUEDGE2(I,J)                                          
 30   CONTINUE                                                                  
                                                                                
      DO 35 K=KTROP+1,KTOP                                                      
      DO 35 J=L1,L2                                                             
      DO 35 I=1,N                                                               
         PSCINV(I,J,K) = 1.                                                     
         PSUEDGE(I,J,K)= 1.                                                     
 35   CONTINUE                                                                  
                                                                                
C     ZERO TEMPORARY MOMENTS.                                                   
      DO 160  J=L1,L2                                                           
         DO 40 K=1,KTOP                                                         
         DO 40 I=1,N1                                                           
            FX (I,K)   = 0.                                                     
            FXX(I,K)   = 0.                                                     
            FY (I,K)   = 0.                                                     
            FYY(I,K)   = 0.                                                     
            FZ (I,K)   = 0.                                                     
            FZZ(I,K)   = 0.                                                     
            FXY(I,K)   = 0.                                                     
            FXZ(I,K)   = 0.                                                     
            FYZ(I,K)   = 0.                                                     
            F0CON(I,K) = 0.                                                     
            F0(I,K)   = 0.                                                      
 40      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 60 K=1,KTOP                                                         
         DO 60 I=1,N                                                            
            SLPMAX = C(I,J,K)                                                   
            S1MAX = 1.5* SLPMAX                                                 
            CX(I,J,K)=MIN(S1MAX,MAX(-S1MAX,CX(I,J,K)))                          
            CXX(I,J,K) = MIN((SLPMAX+SLPMAX-0.333*ABS(CX(I,J,K))),              
     1                 MAX(ABS(CX(I,J,K))-SLPMAX,CXX(I,J,K)))                   
            CXY(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CXY(I,J,K)))                    
            CXZ(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CXZ(I,J,K)))                    
 60      CONTINUE                                                               
C                                                                               
C        BEGIN TROPOSPHERIC LOOP                                                
         DO 100 I=1,N                                                           
            IP = I + 1                                                          
            IF (IP.GT.N) IP=1                                                   
            DO 80  K=1,KTOP                                                     
               IF (U(IP,J,K).GE.0.) THEN                                        
                 EPS=ABS(PSUEDGE(IP,J,K)*PSCINV(I,J,K)*                         
     1           U(IP,J,K)*DTOVDXL(J))                                          
C                IF (EPS.GE.1.0) WRITE (6,*) 'EPS GT 1 ',IP,J,K,EPS,            
C    1            U(IP,J,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*(CX(I,J,K)+                         
     1             (EPS1-EPS)*CXX(I,J,K)))                                      
                  FXS = EPSQ*(CX(I,J,K)+3.0*EPS1*CXX(I,J,K))                    
                  FXXS = EPS*EPSQ*CXX(I,J,K)                                    
                  FYS  = EPS*(CY(I,J,K)+EPS1*CXY(I,J,K))                        
                  FZS  = EPS*(CZ(I,J,K)+EPS1*CXZ(I,J,K))                        
                  FYYS = EPS *CYY(I,J,K)                                        
                  FZZS = EPS *CZZ(I,J,K)                                        
                  FXYS = EPSQ*CXY(I,J,K)                                        
                  FXZS = EPSQ*CXZ(I,J,K)                                        
                  FYZS = EPS *CYZ(I,J,K)                                        
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS * MASS2(I,J,K)                                     
                  F0CON(I,K) = F0CON(I,K) + TEMP                                
                  MASS2(I,J,K) = MASS2(I,J,K) -TEMP                             
                                                                                
C                 READJUST MOMENTS REMAINING IN BOX.                            
                  C(I,J,K) = C(I,J,K) - F0S                                     
                  CX(I,J,K)=EPS1Q*(CX(I,J,K)-3.0*EPS*CXX(I,J,K))                
                  CXX(I,J,K)=EPS1*EPS1Q*CXX(I,J,K)                              
                  CY (I,J,K) = CY(I,J,K)-FYS                                    
                  CZ (I,J,K) = CZ(I,J,K)-FZS                                    
                  CYY(I,J,K) = CYY(I,J,K)-FYYS                                  
                  CZZ(I,J,K) = CZZ(I,J,K)-FZZS                                  
                  CXY(I,J,K) = EPS1Q*CXY(I,J,K)                                 
                  CXZ(I,J,K) = EPS1Q*CXZ(I,J,K)                                 
                  CYZ(I,J,K) = CYZ(I,J,K)-FYZS                                  
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0(I,K) = F0(I,K) + F0S                                       
                  FX(I,K) = FX(I,K) + FXS                                       
                  FXX(I,K) = FXX(I,K) + FXXS                                    
                  FY(I,K) = FY(I,K) + FYS                                       
                  FYY(I,K) = FYY(I,K) + FYYS                                    
                  FZ(I,K) = FZ(I,K) + FZS                                       
                  FZZ(I,K) = FZZ(I,K) + FZZS                                    
                  FXY(I,K) = FXY(I,K) + FXYS                                    
                  FXZ(I,K) = FXZ(I,K) + FXZS                                    
                  FYZ(I,K) = FYZ(I,K) + FYZS                                    
                                                                                
               ELSE                                                             
                EPS=ABS(PSUEDGE(IP,J,K)*PSCINV(IP,J,K)                          
     1           *U(IP,J,K)*DTOVDXL(J))                                         
C                IF (EPS.GE.1.0) WRITE (6,*) 'EPS GT 1 ',IP,J,K,EPS,            
C    1            U(IP,J,K)                                                     
                  EPSQ = EPS*EPS                                                
                  EPS1 = 1.0 - EPS                                              
                  EPS1Q = EPS1 * EPS1                                           
C                 DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                   
                  F0S = EPS * (C(IP,J,K)-EPS1*(CX(IP,J,K)-                      
     1             (EPS1-EPS)*CXX(IP,J,K)))                                     
                  FXS=EPSQ*(CX(IP,J,K)-3.0*EPS1*CXX(IP,J,K))                    
                  FXXS=EPS*EPSQ*CXX(IP,J,K)                                     
                  FYS  = EPS*(CY(IP,J,K)-EPS1*CXY(IP,J,K))                      
                  FZS  = EPS*(CZ(IP,J,K)-EPS1*CXZ(IP,J,K))                      
                  FYYS = EPS *CYY(IP,J,K)                                       
                  FZZS = EPS *CZZ(IP,J,K)                                       
                  FXYS = EPSQ*CXY(IP,J,K)                                       
                  FXZS = EPSQ*CXZ(IP,J,K)                                       
                  FYZS = EPS *CYZ(IP,J,K)                                       
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS*MASS2(IP,J,K)                                      
                  F0CON(I,K) = F0CON(I,K) + TEMP                                
                  MASS2(IP,J,K) = MASS2(IP,J,K) - TEMP                          
                                                                                
C                 READJUST REMAINING MOMENTS.                                   
                  C(IP,J,K) = C(IP,J,K) - F0S                                   
                  CX(IP,J,K) = EPS1Q*(CX(IP,J,K)+3.0*EPS                        
     1             *CXX(IP,J,K))                                                
                  CXX(IP,J,K) = EPS1*EPS1Q*CXX(IP,J,K)                          
                  CY (IP,J,K) = CY(IP,J,K)-FYS                                  
                  CZ (IP,J,K) = CZ(IP,J,K)-FZS                                  
                  CYY(IP,J,K) = CYY(IP,J,K)-FYYS                                
                  CZZ(IP,J,K) = CZZ(IP,J,K)-FZZS                                
                  CXY(IP,J,K) = EPS1Q*CXY(IP,J,K)                               
                  CXZ(IP,J,K) = EPS1Q*CXZ(IP,J,K)                               
                  CYZ(IP,J,K) = CYZ(IP,J,K)-FYZS                                
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0  (I,K) = F0 (I,K) + F0S                                    
                  FX  (I,K) = FX (I,K) + FXS                                    
                  FXX (I,K) = FXX(I,K) + FXXS                                   
                  FY  (I,K) = FY (I,K) + FYS                                    
                  FYY (I,K) = FYY(I,K) + FYYS                                   
                  FZ  (I,K) = FZ (I,K) + FZS                                    
                  FZZ (I,K) = FZZ(I,K) + FZZS                                   
                  FXY (I,K) = FXY(I,K) + FXYS                                   
                  FXZ (I,K) = FXZ(I,K) + FXZS                                   
                  FYZ (I,K) = FYZ(I,K) + FYZS                                   
               END IF                                                           
 80         CONTINUE                                                            
 100     CONTINUE                                                               
                                                                                
C        RECONSTRUCT                                                            
         DO 140 I=1,N                                                           
            IP = I + 1                                                          
            IF (IP.GT.N) IP = 1                                                 
            DO 120 K=1,KTOP                                                     
               IF (U(IP,J,K).GE.0.) THEN                                        
                  EPS = F0CON(I,K) / (F0CON(I,K) + MASS2(IP,J,K))               
                  EPS1 = 1.-EPS                                                 
                  TM=EPS*C(IP,J,K)-EPS1*F0(I,K)                                 
                  C(IP,J,K) = C(IP,J,K) + F0(I,K)                               
                  CXX(IP,J,K)=EPS*EPS*FXX(I,K)+EPS1                            
     1             *EPS1*CXX(IP,J,K)+5.0*(EPS*EPS1                              
     1             *(CX(IP,J,K)-FX(I,K))-(EPS1-EPS)*TM)                         
                  CX(IP,J,K) = EPS*FX(I,K)+EPS1                                 
     1             *CX(IP,J,K)+ 3.0*TM                                          
                  CY(IP,J,K)  = CY(IP,J,K) + FY(I,K)                            
                  CYY(IP,J,K) = CYY(IP,J,K)+ FYY(I,K)                           
                  CZ(IP,J,K)  = CZ(IP,J,K) + FZ (I,K)                           
                  CZZ(IP,J,K) = CZZ(IP,J,K) + FZZ(I,K)                          
                  CXY(IP,J,K) = EPS*FXY(I,K)+EPS1*CXY(IP,J,K)+                  
     1             3.*(-EPS1*FY(I,K)+EPS*CY(IP,J,K))                            
                  CXZ(IP,J,K) = EPS*FXZ(I,K)+EPS1*CXZ(IP,J,K)+                  
     1             3.*(-EPS1*FZ(I,K)+EPS*CZ(IP,J,K))                            
                  CYZ(IP,J,K) = CYZ(IP,J,K) + FYZ(I,K)                          
                  MASS2(IP,J,K) = MASS2(IP,J,K) + F0CON(I,K)                    
                                                                                
               ELSE                                                             
                  EPS = F0CON(I,K) / (F0CON(I,K) + MASS2(I,J,K))                
                  EPS1 = 1. - EPS                                               
                  TM = -EPS*C(I,J,K) + EPS1*F0(I,K)                             
                  C(I,J,K) = C(I,J,K) + F0(I,K)                                 
                  CXX(I,J,K) = EPS*EPS*FXX(I,K)+                                
     1             EPS1*EPS1*CXX(I,J,K) + 5.*                                   
     1             (EPS*EPS1*(-CX(I,J,K)+FX(I,K))+                              
     1             (EPS1-EPS)*TM)                                               
                  CX(I,J,K) = EPS*FX(I,K)+EPS1*CX(I,J,K)+3.*TM                  
                  CY(I,J,K)  = CY(I,J,K) + FY(I,K)                              
                  CYY(I,J,K) = CYY(I,J,K)+ FYY(I,K)                             
                  CZ(I,J,K)  = CZ(I,J,K) + FZ (I,K)                             
                  CZZ(I,J,K) = CZZ(I,J,K) + FZZ(I,K)                            
                  CXY(I,J,K) = EPS*FXY(I,K)+EPS1*CXY(I,J,K)+                    
     1             3.*(EPS1*FY(I,K)-EPS*CY(I,J,K))                              
                  CXZ(I,J,K) = EPS*FXZ(I,K)+EPS1*CXZ(I,J,K)+                    
     1             3.*(EPS1*FZ(I,K)-EPS*CZ(I,J,K))                              
                  CYZ(I,J,K) = CYZ(I,J,K) + FYZ(I,K)                            
                  MASS2(I,J,K) = MASS2(I,J,K) + F0CON(I,K)                      
               END IF                                                           
 120        CONTINUE                                                            
 140     CONTINUE                                                               
 160  CONTINUE                                                                  
                                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE SGPRAP(MASS,MASS2,F1,SGBOT,SGTROP,SGSTRAT,DT,PS1)              
                                                                                
C     THE VERTICAL VELOCITY IS DETERMINED RESIDUALLY FOR USE WITH               
C     PRATHER'S SCHEME.                                                         
C     INITIAL VERSION OCTOBER 25, 1990                                          
C     3-D PRATHER VERSION 1/10/91                                               
                                                                                
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11,N1=n+0)               
                                                                                
      REAL SGBOT(N1,0:JGMAX), SGTROP(N1,0:JGMAX,KTROP),                         
     1 SGSTRAT(N1,0:JGMAX,0:KSTRAT-1),PS1(N1,0:JGMAX),                          
     3 MASS(N1,0:JGMAX,0:KTOP),MASS2(N1,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)                   
                                                                                
      REAL P2(N1,0:JGMAX),P0(N1,0:JGMAX),P1(N1,0:JGMAX),                        
     1 F1(N1,0:JGMAX,KTOP)                                                      
                                                                                
      NJG2 = N*(JGMAX+1)                                                        
                                                                                
      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                                                                       
      KM1 = K-1                                                                 
      DO 50 J=0,JGMAX                                                           
      DO 50 I=1,N                                                               
         IF (F1(I,J,K).GT.0.) THEN                                              
            SGBOT(I,J) = F1(I,J,K) * DELTAS(K)*DT1/MASS2(I,J,KM1)               
         ELSE                                                                   
            SGBOT(I,J) = F1(I,J,K) * DELTAS(K)*DT1/MASS2(I,J,K  )               
         END IF                                                                 
 50   CONTINUE                                                                  
                                                                                
      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      (MASS2(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       MASS2(I,J,KM1)                                                     
         ELSE                                                                   
            SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(K) * DT1/                      
     1       MASS2(I,J,K)                                                       
         END IF                                                                 
 60   CONTINUE                                                                  
                                                                                
      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         (MASS2(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          MASS2(I,J,KM1)                                                  
            ELSE                                                                
               SGTROP(I,J,KM1) = F1(I,J,K) * DELTAS(K) * DT1 /                  
     1          MASS2(I,J,K)                                                    
            END IF                                                              
 80      CONTINUE                                                               
 100  CONTINUE                                                                  
                                                                                
      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       (MASS2(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       MASS2(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       MASS2(I,J,K)                                                       
            SGTROP(I,J,KM1) = SGSTRAT(I,J,KSM1) * CONV                          
         END IF                                                                 
 120  CONTINUE                                                                  
                                                                                
      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       (MASS2(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          MASS2(I,J,KM1)                                                  
            ELSE                                                                
               SGSTRAT(I,J,KSM1) = F1(I,J,K) * DELTAS(K) * DT1/                 
     1          MASS2(I,J,K)                                                    
            END IF                                                              
 140     CONTINUE                                                               
 160  CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE YADVP(C,CX,CXX,CY,CYY,CZ,CZZ,CXY,CXZ,CYZ,MASS,MASS2,           
     1 DT,V,COSCEN,COSEDGE,ADPHI,PSCINV,PSVEDGE)                                
                                                                                
C     NORTH-SOUTH ADVECTION IS PERFORMED WITH PRATHER'S SCHEME.                 
C     FIRST VERSION 1/7/91   COMPLETED?                                         
                                                                                
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11,N1=n+0)               
                                                                                
      REAL CX(N1,0:JGMAX,KTOP),CY(N1,0:JGMAX,KTOP),CZ(N1,0:JGMAX,KTOP),         
     1 CXX(N1,0:JGMAX,KTOP),CYY(N1,0:JGMAX,KTOP),CZZ(N1,0:JGMAX,KTOP),          
     1 CXY(N1,0:JGMAX,KTOP),CXZ(N1,0:JGMAX,KTOP),CYZ(N1,0:JGMAX,KTOP)           
      REAL C(N1,0:JGMAX,KTOP),MASS2(N1,0:JGMAX,0:KTOP)                          
      REAL V(N1,JGMAX,KTOP),COSCEN(0:JGMAX),PSCINV(N1,0:JGMAX)                  
      REAL PSVEDGE(N1,JGMAX),COSEDGE(JGMAX),MASS(N1,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)                   
                                                                                
      REAL                                                                      
     1 FX(N1,0:JGMAX),FY(N1,0:JGMAX),FZ(N1,0:JGMAX),                            
     1 FYY(N1,0:JGMAX),FXX(N1,0:JGMAX),FZZ(N1,0:JGMAX),                         
     2 FXY(N1,0:JGMAX),FYZ(N1,0:JGMAX),FXZ(N1,0:JGMAX),                         
     2 F0(N1,0:JGMAX),F0CON(N1,0:JGMAX)                                         
C                                                                               
C     SET CONSTANTS                                                             
      A = 6.37E6                                                                
C     DTOVDXJ = DT / (A*SIN(ADPHI/A))                                           
      DTOVDXJ=DT/ADPHI                                                          
      PI = 4.0 * ATAN(1.)                                                       
      DPHIPOLE = 1.5 * ADPHI/A                                                  
      PC = 2. * DT / (A*N*SIN(DPHIPOLE))                                        
      DPOLE = SIN(DPHIPOLE) * N * 0.5                                           
                                                                                
C     ZERO TEMPORARY MOMENTS.                                                   
      DO 160  K=1,KTROP                                                         
         DO 40 J=0,JGMAX                                                        
         DO 40 I=1,N1                                                           
            FY(I,J)   = 0.                                                      
            FYY(I,J)   = 0.                                                     
            FX(I,J)   = 0.                                                      
            FXX(I,J)   = 0.                                                     
            FZ (I,J)   = 0.                                                     
            FZZ(I,J)   = 0.                                                     
            FXY(I,J)   = 0.                                                     
            FYZ(I,J)   = 0.                                                     
            FXZ(I,J)   = 0.                                                     
            F0CON(I,J) = 0.                                                     
            F0(I,J)   = 0.                                                      
 40      CONTINUE                                                               
                                                                                
C        LIMIT 1ST AND 2ND MOMENTS TO ASSURE THAT UPDATED MEAN FD IS            
C        GREATER THAN OR EQUAL TO ZERO AFTER TRANSPORT.                         
C        SEE EQN. (39) OF PRATHER: NUMERICAL ADVECTION                          
         DO 60 J=0,JGMAX                                                        
         DO 60 I=1,N                                                            
            SLPMAX = C(I,J,K)                                                   
            S1MAX = 1.5* SLPMAX                                                 
            CY(I,J,K)=MIN(S1MAX,MAX(-S1MAX,CY(I,J,K)))                          
            CYY(I,J,K) = MIN((SLPMAX+SLPMAX-0.333*ABS(CY(I,J,K))),              
     1                 MAX(ABS(CY(I,J,K))-SLPMAX,CYY(I,J,K)))                   
            CXY(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CXY(I,J,K)))                    
            CYZ(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CYZ(I,J,K)))                    
 60      CONTINUE                                                               
C                                                                               
C        ENSURE THAT HIGHER ORDER MOMENTS ARE ZERO IN POLAR CAPS.               
         DO 70 I=1,N                                                            
            CX(I,0    ,K) = 0.                                                  
            CX(I,JGMAX,K) = 0.                                                  
            CY(I,0    ,K) = 0.                                                  
            CY(I,JGMAX,K) = 0.                                                  
            CZ(I,0    ,K) = 0.                                                  
            CZ(I,JGMAX,K) = 0.                                                  
            CXX(I,0    ,K) = 0.                                                 
            CXX(I,JGMAX,K) = 0.                                                 
            CYY(I,0    ,K) = 0.                                                 
            CYY(I,JGMAX,K) = 0.                                                 
            CZZ(I,0    ,K) = 0.                                                 
            CZZ(I,JGMAX,K) = 0.                                                 
            CXY(I,0    ,K) = 0.                                                 
            CXY(I,JGMAX,K) = 0.                                                 
            CXZ(I,0    ,K) = 0.                                                 
            CXZ(I,JGMAX,K) = 0.                                                 
            CYZ(I,0    ,K) = 0.                                                 
            CYZ(I,JGMAX,K) = 0.                                                 
 70      CONTINUE                                                               
                                                                                
C        SOUTH POLE.                                                            
         J=0                                                                    
         JP = J + 1                                                             
         F0SA = 0.                                                              
         TEMPA = 0.                                                             
         F0SB = 0.                                                              
         TEMPB = 0.                                                             
         DO 75  I=1,N                                                           
            IF (V(I,JP,K).GE.0.) THEN                                           
C              CALCULATE FLUX OUT OF POLAR CAP.                                 
                F0S= PSVEDGE(I,JP)*PSCINV(I,J)*V(I,JP,K)*DT                     
     1           *C(I,J,K)    /(DPOLE*A)                                        
                TEMP=PSVEDGE(I,JP)*PSCINV(I,J)*V(I,JP,K)*DT                     
     1           *MASS2(I,J,K)/(DPOLE*A)                                        
                F0SA = F0SA + F0S                                               
                TEMPA = TEMPA + TEMP                                            
               F0(I,J)=F0(I,J)+F0S                                              
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
                                                                                
            ELSE                                                                
               EPS=ABS(PSVEDGE(I,JP)*PSCINV(I,JP)*V(I,JP,K)*DTOVDXJ             
     1                       *COSEDGE(JP)/COSCEN(JP))                           
               EPSQ = EPS*EPS                                                   
               EPS1 = 1.0 - EPS                                                 
               EPS1Q = EPS1 * EPS1                                              
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                      
               F0S = EPS * (C(I,JP,K)-EPS1*(CY(I,JP,K)-                         
     1          (EPS1-EPS)*CYY(I,JP,K)))                                        
               F0SB = F0SB + F0S                                                
               FYS=EPSQ*(CY(I,JP,K)-3.0*EPS1*CYY(I,JP,K))                       
               FYYS=EPS*EPSQ*CYY(I,JP,K)                                        
               FXS  = EPS*(CX(I,JP,K)-EPS1*CXY(I,JP,K))                         
               FZS  = EPS*(CZ(I,JP,K)-EPS1*CYZ(I,JP,K))                         
               FXXS = EPS *CXX(I,JP,K)                                          
               FZZS = EPS *CZZ(I,JP,K)                                          
               FXYS = EPSQ*CXY(I,JP,K)                                          
               FYZS = EPSQ*CYZ(I,JP,K)                                          
               FXZS = EPS *CXZ(I,JP,K)                                          
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS*MASS2(I,JP,K)                                         
               TEMPB = TEMPB + TEMP                                             
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
               MASS2(I,JP,K) = MASS2(I,JP,K) - TEMP                             
                                                                                
C              READJUST REMAINING MOMENTS.                                      
               C(I,JP,K) = C(I,JP,K) - F0S                                      
               CY(I,JP,K) = EPS1Q*(CY(I,JP,K)+3.0*EPS                           
     1          *CYY(I,JP,K))                                                   
               CYY(I,JP,K) = EPS1*EPS1Q*CYY(I,JP,K)                             
               CX(I,JP,K) = CX(I,JP,K)-FXS                                      
               CZ (I,JP,K) = CZ(I,JP,K)-FZS                                     
               CXX(I,JP,K) = CXX(I,JP,K)-FXXS                                   
               CZZ(I,JP,K) = CZZ(I,JP,K)-FZZS                                   
               CXY(I,JP,K) = EPS1Q*CXY(I,JP,K)                                  
               CYZ(I,JP,K) = EPS1Q*CYZ(I,JP,K)                                  
               CXZ(I,JP,K) = CXZ(I,JP,K)-FXZS                                   
                                                                                
C              SET TEMPORARY MOMENTS                                            
               F0  (I,J) = F0 (I,J) + F0S                                       
               FY(I,J) = FY(I,J) + FYS                                          
               FYY(I,J) = FYY(I,J) + FYYS                                       
               FX(I,J) = FX(I,J) + FXS                                          
               FXX(I,J) = FXX(I,J) + FXXS                                       
               FZ  (I,J) = FZ (I,J) + FZS                                       
               FZZ (I,J) = FZZ(I,J) + FZZS                                      
               FXY(I,J) = FXY(I,J) + FXYS                                       
               FYZ(I,J) = FYZ(I,J) + FYZS                                       
               FXZ(I,J) = FXZ(I,J) + FXZS                                       
            END IF                                                              
 75      CONTINUE                                                               
C        ADJUST MASS IN POLAR CAP.                                              
                                                                                
         DO 79 I=1,N                                                            
            C(I,J,K) = C(I,J,K) - F0SA                                          
            MASS2(I,J,K) = MASS2(I,J,K) - TEMPA                                 
 79      CONTINUE                                                               
                                                                                
         DO 100 J=1,JGMAX-2                                                     
            JP = J + 1                                                          
            DO 80  I=1,N                                                        
               IF (V(I,JP,K).GE.0.) THEN                                        
                  EPS=ABS(PSVEDGE(I,JP)*PSCINV(I,J)*V(I,JP,K)*DTOVDXJ*          
     1             COSEDGE(JP)/COSCEN(J))                                       
                  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*(CY(I,J,K)+                         
     1             (EPS1-EPS)*CYY(I,J,K)))                                      
                  FYS = EPSQ*(CY(I,J,K)+3.0*EPS1*CYY(I,J,K))                    
                  FYYS = EPS*EPSQ*CYY(I,J,K)                                    
                  FXS  = EPS*(CX(I,J,K)+EPS1*CXY(I,J,K))                        
                  FZS  = EPS*(CZ(I,J,K)+EPS1*CYZ(I,J,K))                        
                  FXXS = EPS *CXX(I,J,K)                                        
                  FZZS = EPS *CZZ(I,J,K)                                        
                  FXYS = EPSQ*CXY(I,J,K)                                        
                  FYZS = EPSQ*CYZ(I,J,K)                                        
                  FXZS = EPS *CXZ(I,J,K)                                        
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS * MASS2(I,J,K)                                     
                  F0CON(I,J) = F0CON(I,J) + TEMP                                
                  MASS2(I,J,K) = MASS2(I,J,K) -TEMP                             
                                                                                
C                 READJUST MOMENTS REMAINING IN BOX.                            
                  C(I,J,K) = C(I,J,K) - F0S                                     
                  CY(I,J,K)=EPS1Q*(CY(I,J,K)-3.0*EPS*CYY(I,J,K))                
                  CYY(I,J,K)=EPS1*EPS1Q*CYY(I,J,K)                              
                  CX(I,J,K) = CX(I,J,K)-FXS                                     
                  CZ (I,J,K) = CZ(I,J,K)-FZS                                    
                  CXX(I,J,K) = CXX(I,J,K)-FXXS                                  
                  CZZ(I,J,K) = CZZ(I,J,K)-FZZS                                  
                  CXY(I,J,K) = EPS1Q*CXY(I,J,K)                                 
                  CYZ(I,J,K) = EPS1Q*CYZ(I,J,K)                                 
                  CXZ(I,J,K) = CXZ(I,J,K)-FXZS                                  
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0(I,J) = F0(I,J) + F0S                                       
                  FY(I,J) = FY(I,J) + FYS                                       
                  FYY(I,J) = FYY(I,J) + FYYS                                    
                  FX(I,J) = FX(I,J) + FXS                                       
                  FXX(I,J) = FXX(I,J) + FXXS                                    
                  FZ(I,J) = FZ(I,J) + FZS                                       
                  FZZ(I,J) = FZZ(I,J) + FZZS                                    
                  FXY(I,J) = FXY(I,J) + FXYS                                    
                  FYZ(I,J) = FYZ(I,J) + FYZS                                    
                  FXZ(I,J) = FXZ(I,J) + FXZS                                    
                                                                                
               ELSE                                                             
                  EPS=ABS(PSVEDGE(I,JP)*PSCINV(I,JP)*V(I,JP,K)*DTOVDXJ*         
     1             COSEDGE(JP)/COSCEN(JP))                                      
                  EPSQ = EPS*EPS                                                
                  EPS1 = 1.0 - EPS                                              
                  EPS1Q = EPS1 * EPS1                                           
C                 DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                   
                  F0S = EPS * (C(I,JP,K)-EPS1*(CY(I,JP,K)-                      
     1             (EPS1-EPS)*CYY(I,JP,K)))                                     
                  FYS=EPSQ*(CY(I,JP,K)-3.0*EPS1*CYY(I,JP,K))                    
                  FYYS=EPS*EPSQ*CYY(I,JP,K)                                     
                  FXS  = EPS*(CX(I,JP,K)-EPS1*CXY(I,JP,K))                      
                  FZS  = EPS*(CZ(I,JP,K)-EPS1*CYZ(I,JP,K))                      
                  FXXS = EPS *CXX(I,JP,K)                                       
                  FZZS = EPS *CZZ(I,JP,K)                                       
                  FXYS = EPSQ*CXY(I,JP,K)                                       
                  FYZS = EPSQ*CYZ(I,JP,K)                                       
                  FXZS = EPS *CXZ(I,JP,K)                                       
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS*MASS2(I,JP,K)                                      
                  F0CON(I,J) = F0CON(I,J) + TEMP                                
                  MASS2(I,JP,K) = MASS2(I,JP,K) - TEMP                          
                                                                                
C                 READJUST REMAINING MOMENTS.                                   
                  C(I,JP,K) = C(I,JP,K) - F0S                                   
                  CY(I,JP,K) = EPS1Q*(CY(I,JP,K)+3.0*EPS                        
     1             *CYY(I,JP,K))                                                
                  CYY(I,JP,K) = EPS1*EPS1Q*CYY(I,JP,K)                          
                  CX(I,JP,K) = CX(I,JP,K)-FXS                                   
                  CZ (I,JP,K) = CZ(I,JP,K)-FZS                                  
                  CXX(I,JP,K) = CXX(I,JP,K)-FXXS                                
                  CZZ(I,JP,K) = CZZ(I,JP,K)-FZZS                                
                  CXY(I,JP,K) = EPS1Q*CXY(I,JP,K)                               
                  CYZ(I,JP,K) = EPS1Q*CYZ(I,JP,K)                               
                  CXZ(I,JP,K) = CXZ(I,JP,K)-FXZS                                
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0  (I,J) = F0 (I,J) + F0S                                    
                  FY(I,J) = FY(I,J) + FYS                                       
                  FYY(I,J) = FYY(I,J) + FYYS                                    
                  FX(I,J) = FX(I,J) + FXS                                       
                  FXX(I,J) = FXX(I,J) + FXXS                                    
                  FZ  (I,J) = FZ (I,J) + FZS                                    
                  FZZ (I,J) = FZZ(I,J) + FZZS                                   
                  FXY(I,J) = FXY(I,J) + FXYS                                    
                  FYZ(I,J) = FYZ(I,J) + FYZS                                    
                  FXZ(I,J) = FXZ(I,J) + FXZS                                    
               END IF                                                           
 80         CONTINUE                                                            
 100     CONTINUE                                                               
                                                                                
C        NORTH POLE                                                             
         J=JGMAX-1                                                              
         JP = J + 1                                                             
C        ADJUST MASS IN POLAR CAP.                                              
         F0SC = 0.                                                              
         F0SD = 0.                                                              
         TEMPC = 0.                                                             
         TEMPD = 0.                                                             
         DO 105  I=1,N                                                          
            IF (V(I,JP,K).GE.0.) THEN                                           
                EPS=ABS(PSVEDGE(I,JP)*PSCINV(I,J)*DTOVDXJ*V(I,JP,K)             
     1          *COSEDGE(JP)/COSCEN(J))                                         
               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*(CY(I,J,K)+                            
     1          (EPS1-EPS)*CYY(I,J,K)))                                         
               F0SC = F0SC + F0S                                                
               FYS = EPSQ*(CY(I,J,K)+3.0*EPS1*CYY(I,J,K))                       
               FYYS = EPS*EPSQ*CYY(I,J,K)                                       
               FXS  = EPS*(CX(I,J,K)+EPS1*CXY(I,J,K))                           
               FZS  = EPS*(CZ(I,J,K)+EPS1*CYZ(I,J,K))                           
               FXXS = EPS *CXX(I,J,K)                                           
               FZZS = EPS *CZZ(I,J,K)                                           
               FXYS = EPSQ*CXY(I,J,K)                                           
               FYZS = EPSQ*CYZ(I,J,K)                                           
               FXZS = EPS *CXZ(I,J,K)                                           
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS * MASS2(I,J,K)                                        
               TEMPC = TEMPC + TEMP                                             
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
               MASS2(I,J,K) = MASS2(I,J,K) -TEMP                                
                                                                                
C              READJUST MOMENTS REMAINING IN BOX.                               
               C(I,J,K) = C(I,J,K) - F0S                                        
               CY(I,J,K)=EPS1Q*(CY(I,J,K)-3.0*EPS*CYY(I,J,K))                   
               CYY(I,J,K)=EPS1*EPS1Q*CYY(I,J,K)                                 
               CX(I,J,K) = CX(I,J,K)-FXS                                        
               CZ (I,J,K) = CZ(I,J,K)-FZS                                       
               CXX(I,J,K) = CXX(I,J,K)-FXXS                                     
               CZZ(I,J,K) = CZZ(I,J,K)-FZZS                                     
               CXY(I,J,K) = EPS1Q*CXY(I,J,K)                                    
               CYZ(I,J,K) = EPS1Q*CYZ(I,J,K)                                    
               CXZ(I,J,K) = CXZ(I,J,K)-FXZS                                     
                                                                                
C              SET TEMPORARY MOMENTS                                            
               F0(I,J) = F0(I,J) + F0S                                          
               FY(I,J) = FY(I,J) + FYS                                          
               FYY(I,J) = FYY(I,J) + FYYS                                       
               FX(I,J) = FX(I,J) + FXS                                          
               FXX(I,J) = FXX(I,J) + FXXS                                       
               FZ(I,J) = FZ(I,J) + FZS                                          
               FZZ(I,J) = FZZ(I,J) + FZZS                                       
               FXY(I,J) = FXY(I,J) + FXYS                                       
               FYZ(I,J) = FYZ(I,J) + FYZS                                       
               FXZ(I,J) = FXZ(I,J) + FXZS                                       
                                                                                
            ELSE                                                                
               F0S=  ABS(PSVEDGE(I,JP)*PSCINV(I,JP)*DT*V(I,JP,K)*               
     1          C(I,JP,K)/(A*DPOLE))                                            
               TEMP= ABS(PSVEDGE(I,JP)*PSCINV(I,JP)*DT*V(I,JP,K)*               
     1          MASS2(I,JP,K)/(A*DPOLE))                                        
               F0SD = F0SD + F0S                                                
               TEMPD = TEMPD + TEMP                                             
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
C              SET TEMPORARY MOMENTS                                            
               F0(I,J)=F0(I,J)+ F0S                                             
            END IF                                                              
 105     CONTINUE                                                               
                                                                                
         DO 179 I=1,N                                                           
            C(I,JP,K) = C(I,JP,K) - F0SD                                        
            MASS2(I,JP,K) = MASS2(I,JP,K) - TEMPD                               
 179     CONTINUE                                                               
                                                                                
C        RECONSTRUCT TROPOSPHERE.                                               
         J=0                                                                    
         JP = J + 1                                                             
         DO 478 I=1,N                                                           
            C(I,J,K) = C(I,J,K) + F0SB                                          
            MASS2(I,J,K) = MASS2(I,J,K) + TEMPB                                 
 478     CONTINUE                                                               
                                                                                
         DO 110 I=1,N                                                           
            IF (V(I,JP,K).GE.0.) THEN                                           
               EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,JP,K))                  
               EPS1 = 1.-EPS                                                    
               TM=EPS*C(I,JP,K)-EPS1*F0(I,J)                                    
               C(I,JP,K) = C(I,JP,K) + F0(I,J)                                  
               CYY(I,JP,K)=EPS*EPS*FYY(I,J)+EPS1                               
     1          *EPS1*CYY(I,JP,K)+5.0*(EPS*EPS1                                 
     1          *(CY(I,JP,K)-FY(I,J))-(EPS1-EPS)*TM)                            
               CY(I,JP,K) = EPS*FY(I,J)+EPS1                                    
     1          *CY(I,JP,K)+ 3.0*TM                                             
               CX(I,JP,K)  = CX(I,JP,K) + FX(I,J)                               
               CXX(I,JP,K) = CXX(I,JP,K)+ FXX(I,J)                              
               CZ(I,JP,K)  = CZ(I,JP,K) + FZ (I,J)                              
               CZZ(I,JP,K) = CZZ(I,JP,K) + FZZ(I,J)                             
               CXY(I,JP,K) = EPS*FXY(I,J)+EPS1*CXY(I,JP,K)+                     
     1          3.*(-EPS1*FX(I,J)+EPS*CX(I,JP,K))                               
               CYZ(I,JP,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,JP,K)+                     
     1          3.*(-EPS1*FZ(I,J)+EPS*CZ(I,JP,K))                               
               CXZ(I,JP,K) = CXZ(I,JP,K) + FXZ(I,J)                             
               MASS2(I,JP,K) = MASS2(I,JP,K) + F0CON(I,J)                       
                                                                                
            END IF                                                              
 110     CONTINUE                                                               
                                                                                
         DO 140 J=1,JGMAX-2                                                     
            JP = J + 1                                                          
            DO 120 I=1,N                                                        
               IF (V(I,JP,K).GE.0.) THEN                                        
                  EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,JP,K))               
                  EPS1 = 1.-EPS                                                 
                  TM=EPS*C(I,JP,K)-EPS1*F0(I,J)                                 
                  C(I,JP,K) = C(I,JP,K) + F0(I,J)                               
                  CYY(I,JP,K)=EPS*EPS*FYY(I,J)+EPS1                            
     1             *EPS1*CYY(I,JP,K)+5.0*(EPS*EPS1                              
     1             *(CY(I,JP,K)-FY(I,J))-(EPS1-EPS)*TM)                         
                  CY(I,JP,K) = EPS*FY(I,J)+EPS1                                 
     1             *CY(I,JP,K)+ 3.0*TM                                          
                  CX(I,JP,K)  = CX(I,JP,K) + FX(I,J)                            
                  CXX(I,JP,K) = CXX(I,JP,K)+ FXX(I,J)                           
                  CZ(I,JP,K)  = CZ(I,JP,K) + FZ (I,J)                           
                  CZZ(I,JP,K) = CZZ(I,JP,K) + FZZ(I,J)                          
                  CXY(I,JP,K) = EPS*FXY(I,J)+EPS1*CXY(I,JP,K)+                  
     1             3.*(-EPS1*FX(I,J)+EPS*CX(I,JP,K))                            
                  CYZ(I,JP,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,JP,K)+                  
     1             3.*(-EPS1*FZ(I,J)+EPS*CZ(I,JP,K))                            
                  CXZ(I,JP,K) = CXZ(I,JP,K) + FXZ(I,J)                          
                  MASS2(I,JP,K) = MASS2(I,JP,K) + F0CON(I,J)                    
                                                                                
               ELSE                                                             
                  EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,J,K))                
                  EPS1 = 1. - EPS                                               
                  TM = -EPS*C(I,J,K) + EPS1*F0(I,J)                             
                  C(I,J,K) = C(I,J,K) + F0(I,J)                                 
                  CYY(I,J,K) = EPS*EPS*FYY(I,J)+                                
     1             EPS1*EPS1*CYY(I,J,K) + 5.*                                   
     1             (EPS*EPS1*(-CY(I,J,K)+FY(I,J))+                              
     1             (EPS1-EPS)*TM)                                               
                  CY(I,J,K) = EPS*FY(I,J)+EPS1*CY(I,J,K)+3.*TM                  
                  CX(I,J,K)  = CX(I,J,K) + FX(I,J)                              
                  CXX(I,J,K) = CXX(I,J,K)+ FXX(I,J)                             
                  CZ(I,J,K)  = CZ(I,J,K) + FZ (I,J)                             
                  CZZ(I,J,K) = CZZ(I,J,K) + FZZ(I,J)                            
                  CXY(I,J,K) = EPS*FXY(I,J)+EPS1*CXY(I,J,K)+                    
     1             3.*(EPS1*FX(I,J)-EPS*CX(I,J,K))                              
                  CYZ(I,J,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,J,K)+                    
     1             3.*(EPS1*FZ(I,J)-EPS*CZ(I,J,K))                              
                  CXZ(I,J,K) = CXZ(I,J,K) + FXZ(I,J)                            
                  MASS2(I,J,K) = MASS2(I,J,K) + F0CON(I,J)                      
               END IF                                                           
 120        CONTINUE                                                            
 140     CONTINUE                                                               
                                                                                
         J=JGMAX-1                                                              
         JP = J + 1                                                             
         DO 183 I=1,N                                                           
            C(I,JP,K) = C(I,JP,K) + F0SC                                        
            MASS2(I,JP,K) = MASS2(I,JP,K) + TEMPC                               
 183     CONTINUE                                                               
                                                                                
         DO 145 I=1,N                                                           
            IF (V(I,JP,K).LT.0.) THEN                                           
               EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,J,K))                   
               EPS1 = 1. - EPS                                                  
               TM = -EPS*C(I,J,K) + EPS1*F0(I,J)                                
               C(I,J,K) = C(I,J,K) + F0(I,J)                                    
               CYY(I,J,K) = EPS*EPS*FYY(I,J)+                                   
     1          EPS1*EPS1*CYY(I,J,K) + 5.*                                      
     1          (EPS*EPS1*(-CY(I,J,K)+FY(I,J))+                                 
     1          (EPS1-EPS)*TM)                                                  
               CY(I,J,K) = EPS*FY(I,J)+EPS1*CY(I,J,K)+3.*TM                     
               CX(I,J,K)  = CX(I,J,K) + FX(I,J)                                 
               CXX(I,J,K) = CXX(I,J,K)+ FXX(I,J)                                
               CZ(I,J,K)  = CZ(I,J,K) + FZ (I,J)                                
               CZZ(I,J,K) = CZZ(I,J,K) + FZZ(I,J)                               
               CXY(I,J,K) = EPS*FXY(I,J)+EPS1*CXY(I,J,K)+                       
     1          3.*(EPS1*FX(I,J)-EPS*CX(I,J,K))                                 
               CYZ(I,J,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,J,K)+                       
     1          3.*(EPS1*FZ(I,J)-EPS*CZ(I,J,K))                                 
               CXZ(I,J,K) = CXZ(I,J,K) + FXZ(I,J)                               
               MASS2(I,J,K) = MASS2(I,J,K) + F0CON(I,J)                         
            END IF                                                              
 145     CONTINUE                                                               
 160  CONTINUE                                                                  
                                                                                
C     BEGIN STRATOSPHERIC LOOP                                                  
                                                                                
      DO 260 K=KTROP+1,KTOP                                                     
         DO 180 J=0,JGMAX                                                       
         DO 180 I=1,N1                                                          
            FY(I,J)   = 0.                                                      
            FYY(I,J)   = 0.                                                     
            FX(I,J)   = 0.                                                      
            FXX(I,J)   = 0.                                                     
            FZ (I,J)   = 0.                                                     
            FZZ(I,J)   = 0.                                                     
            FXY(I,J)   = 0.                                                     
            FYZ(I,J)   = 0.                                                     
            FXZ(I,J)   = 0.                                                     
            F0CON(I,J) = 0.                                                     
            F0  (I,J)  = 0.                                                     
 180      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 190 J=0,JGMAX                                                       
         DO 190 I=1,N                                                           
            SLPMAX = C(I,J,K)                                                   
            S1MAX = 1.5* SLPMAX                                                 
            CY(I,J,K)=MIN(S1MAX,MAX(-S1MAX,CY(I,J,K)))                          
            CYY(I,J,K) = MIN((SLPMAX+SLPMAX-0.333*ABS(CY(I,J,K))),              
     1                 MAX(ABS(CY(I,J,K))-SLPMAX,CYY(I,J,K)))                   
            CXY(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CXY(I,J,K)))                    
            CYZ(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CYZ(I,J,K)))                    
 190     CONTINUE                                                               
                                                                                
C        ASSUME WELL MIXED POLAR CAPS                                           
         DO 192 I=1,N                                                           
            CX(I,0    ,K) = 0.                                                  
            CX(I,JGMAX,K) = 0.                                                  
            CY(I,0    ,K) = 0.                                                  
            CY(I,JGMAX,K) = 0.                                                  
            CZ(I,0    ,K) = 0.                                                  
            CZ(I,JGMAX,K) = 0.                                                  
            CXX(I,0    ,K) = 0.                                                 
            CXX(I,JGMAX,K) = 0.                                                 
            CYY(I,0    ,K) = 0.                                                 
            CYY(I,JGMAX,K) = 0.                                                 
            CZZ(I,0    ,K) = 0.                                                 
            CZZ(I,JGMAX,K) = 0.                                                 
            CXY(I,0    ,K) = 0.                                                 
            CXY(I,JGMAX,K) = 0.                                                 
            CXZ(I,0    ,K) = 0.                                                 
            CXZ(I,JGMAX,K) = 0.                                                 
            CYZ(I,0    ,K) = 0.                                                 
            CYZ(I,JGMAX,K) = 0.                                                 
 192     CONTINUE                                                               
                                                                                
C        SOUTH POLE.                                                            
         J=0                                                                    
         JP = J + 1                                                             
         F0SA = 0.                                                              
         F0SB = 0.                                                              
         TEMPA = 0.                                                             
         TEMPB = 0.                                                             
         DO 193  I=1,N                                                          
            IF (V(I,JP,K).GE.0.) THEN                                           
C              CALCULATE FLUX OUT OF POLAR CAP.                                 
               F0S= V(I,JP,K)*DT                                                
     1          *C(I,J,K)    /(DPOLE*A)                                         
               TEMP=V(I,JP,K)*DT                                                
     1          *MASS2(I,J,K)/(DPOLE*A)                                         
               F0SA = F0SA + F0S                                                
               TEMPA = TEMPA + TEMP                                             
               F0(I,J)=F0(I,J)+F0S                                              
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
                                                                                
            ELSE                                                                
               EPS=ABS(V(I,JP,K)*DTOVDXJ*COSEDGE(JP)/COSCEN(JP))                
               EPSQ = EPS*EPS                                                   
               EPS1 = 1.0 - EPS                                                 
               EPS1Q = EPS1 * EPS1                                              
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                      
               F0S = EPS * (C(I,JP,K)-EPS1*(CY(I,JP,K)-                         
     1          (EPS1-EPS)*CYY(I,JP,K)))                                        
               F0SB = F0SB + F0S                                                
               FYS=EPSQ*(CY(I,JP,K)-3.0*EPS1*CYY(I,JP,K))                       
               FYYS=EPS*EPSQ*CYY(I,JP,K)                                        
               FXS  = EPS*(CX(I,JP,K)-EPS1*CXY(I,JP,K))                         
               FZS  = EPS*(CZ(I,JP,K)-EPS1*CYZ(I,JP,K))                         
               FXXS = EPS *CXX(I,JP,K)                                          
               FZZS = EPS *CZZ(I,JP,K)                                          
               FXYS = EPSQ*CXY(I,JP,K)                                          
               FYZS = EPSQ*CYZ(I,JP,K)                                          
               FXZS = EPS *CXZ(I,JP,K)                                          
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS*MASS2(I,JP,K)                                         
               TEMPB = TEMPB + TEMP                                             
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
               MASS2(I,JP,K) = MASS2(I,JP,K) - TEMP                             
                                                                                
C              READJUST REMAINING MOMENTS.                                      
               C(I,JP,K) = C(I,JP,K) - F0S                                      
               CY(I,JP,K) = EPS1Q*(CY(I,JP,K)+3.0*EPS                           
     1          *CYY(I,JP,K))                                                   
               CYY(I,JP,K) = EPS1*EPS1Q*CYY(I,JP,K)                             
               CX(I,JP,K) = CX(I,JP,K)-FXS                                      
               CZ (I,JP,K) = CZ(I,JP,K)-FZS                                     
               CXX(I,JP,K) = CXX(I,JP,K)-FXXS                                   
               CZZ(I,JP,K) = CZZ(I,JP,K)-FZZS                                   
               CXY(I,JP,K) = EPS1Q*CXY(I,JP,K)                                  
               CYZ(I,JP,K) = EPS1Q*CYZ(I,JP,K)                                  
               CXZ(I,JP,K) = CXZ(I,JP,K)-FXZS                                   
                                                                                
C              SET TEMPORARY MOMENTS                                            
               F0  (I,J) = F0 (I,J) + F0S                                       
               FY(I,J) = FY(I,J) + FYS                                          
               FYY(I,J) = FYY(I,J) + FYYS                                       
               FX(I,J) = FX(I,J) + FXS                                          
               FXX(I,J) = FXX(I,J) + FXXS                                       
               FZ  (I,J) = FZ (I,J) + FZS                                       
               FZZ (I,J) = FZZ(I,J) + FZZS                                      
               FXY(I,J) = FXY(I,J) + FXYS                                       
               FYZ(I,J) = FYZ(I,J) + FYZS                                       
               FXZ(I,J) = FXZ(I,J) + FXZS                                       
            END IF                                                              
 193     CONTINUE                                                               
                                                                                
*********                                                                       
         DO 279 I=1,N                                                           
            C(I,J,K) = C(I,J,K) - F0SA                                          
            MASS2(I,J,K) = MASS2(I,J,K) - TEMPA                                 
 279     CONTINUE                                                               
                                                                                
         DO 200 J=1,JGMAX-2                                                     
            JP = J + 1                                                          
C           TO VECTORIZE 195 IF IS BROKEN INTO 2 PARTS.  I AM NOT               
C           SURE WHY THIS IS NECESSARY.                                         
            DO 195  I=1,N                                                       
               IF (V(I,JP,K).GE.0.) THEN                                        
                  EPS = ABS(V(I,JP,K)*DTOVDXJ*COSEDGE(JP)/COSCEN(J))            
                  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*(CY(I,J,K)+                         
     1             (EPS1-EPS)*CYY(I,J,K)))                                      
                  FYS = EPSQ*(CY(I,J,K)+3.0*EPS1*CYY(I,J,K))                    
                  FYYS = EPS*EPSQ*CYY(I,J,K)                                    
                  FXS  = EPS*(CX(I,J,K)+EPS1*CXY(I,J,K))                        
                  FZS  = EPS*(CZ(I,J,K)+EPS1*CYZ(I,J,K))                        
                  FXXS = EPS *CXX(I,J,K)                                        
                  FZZS = EPS *CZZ(I,J,K)                                        
                  FXYS = EPSQ*CXY(I,J,K)                                        
                  FYZS = EPSQ*CYZ(I,J,K)                                        
                  FXZS = EPS *CXZ(I,J,K)                                        
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS * MASS2(I,J,K)                                     
                  F0CON(I,J) = F0CON(I,J) + TEMP                                
                  MASS2(I,J,K) = MASS2(I,J,K) -TEMP                             
                                                                                
C                 READJUST MOMENTS REMAINING IN BOX.                            
                  C(I,J,K) = C(I,J,K) - F0S                                     
                  CY(I,J,K)=EPS1Q*(CY(I,J,K)-3.0*EPS*CYY(I,J,K))                
                  CYY(I,J,K)=EPS1*EPS1Q*CYY(I,J,K)                              
                  CX(I,J,K) = CX(I,J,K)-FXS                                     
                  CZ (I,J,K) = CZ(I,J,K)-FZS                                    
                  CXX(I,J,K) = CXX(I,J,K)-FXXS                                  
                  CZZ(I,J,K) = CZZ(I,J,K)-FZZS                                  
                  CXY(I,J,K) = EPS1Q*CXY(I,J,K)                                 
                  CYZ(I,J,K) = EPS1Q*CYZ(I,J,K)                                 
                  CXZ(I,J,K) = CXZ(I,J,K)-FXZS                                  
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0(I,J) = F0(I,J) + F0S                                       
                  FY(I,J) = FY(I,J) + FYS                                       
                  FYY(I,J) = FYY(I,J) + FYYS                                    
                  FX(I,J) = FX(I,J) + FXS                                       
                  FXX(I,J) = FXX(I,J) + FXXS                                    
                  FZ(I,J) = FZ(I,J) + FZS                                       
                  FZZ(I,J) = FZZ(I,J) + FZZS                                    
                  FXY(I,J) = FXY(I,J) + FXYS                                    
                  FYZ(I,J) = FYZ(I,J) + FYZS                                    
                  FXZ(I,J) = FXZ(I,J) + FXZS                                    
               END IF                                                           
 195        CONTINUE                                                            
                                                                                
            DO 198 I=1,N                                                        
               IF (V(I,JP,K).LT.0.) THEN                                        
                  EPS = ABS(V(I,JP,K)*DTOVDXJ*COSEDGE(JP)/COSCEN(JP))           
                  EPSQ = EPS*EPS                                                
                  EPS1 = 1.0 - EPS                                              
                  EPS1Q = EPS1 * EPS1                                           
C                 DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                   
                  F0S = EPS * (C(I,JP,K)-EPS1*(CY(I,JP,K)-                      
     1             (EPS1-EPS)*CYY(I,JP,K)))                                     
                  FYS=EPSQ*(CY(I,JP,K)-3.0*EPS1*CYY(I,JP,K))                    
                  FYYS=EPS*EPSQ*CYY(I,JP,K)                                     
                  FXS  = EPS*(CX(I,JP,K)-EPS1*CXY(I,JP,K))                      
                  FZS  = EPS*(CZ(I,JP,K)-EPS1*CYZ(I,JP,K))                      
                  FXXS = EPS *CXX(I,JP,K)                                       
                  FZZS = EPS *CZZ(I,JP,K)                                       
                  FXYS = EPSQ*CXY(I,JP,K)                                       
                  FYZS = EPSQ*CYZ(I,JP,K)                                       
C                 FXZS = EPS *CXZ(I,JP,K)                                       
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS*MASS2(I,JP,K)                                      
                  F0CON(I,J) = F0CON(I,J) + TEMP                                
                  MASS2(I,JP,K) = MASS2(I,JP,K) - TEMP                          
                                                                                
C                 READJUST REMAINING MOMENTS.                                   
                  C(I,JP,K) = C(I,JP,K) - F0S                                   
                  CY(I,JP,K) = EPS1Q*(CY(I,JP,K)+3.0*EPS                        
     1             *CYY(I,JP,K))                                                
                  CYY(I,JP,K) = EPS1*EPS1Q*CYY(I,JP,K)                          
                  CX(I,JP,K) = CX(I,JP,K)-FXS                                   
                  CZ (I,JP,K) = CZ(I,JP,K)-FZS                                  
                  CXX(I,JP,K) = CXX(I,JP,K)-FXXS                                
                  CZZ(I,JP,K) = CZZ(I,JP,K)-FZZS                                
                  CXY(I,JP,K) = EPS1Q*CXY(I,JP,K)                               
                  CYZ(I,JP,K) = EPS1Q*CYZ(I,JP,K)                               
                  CXZ(I,JP,K) = CXZ(I,JP,K)-FXZS                                
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0  (I,J) = F0 (I,J) + F0S                                    
                  FY(I,J) = FY(I,J) + FYS                                       
                  FYY(I,J) = FYY(I,J) + FYYS                                    
                  FX(I,J) = FX(I,J) + FXS                                       
                  FXX(I,J) = FXX(I,J) + FXXS                                    
                  FZ  (I,J) = FZ (I,J) + FZS                                    
                  FZZ (I,J) = FZZ(I,J) + FZZS                                   
                  FXY(I,J) = FXY(I,J) + FXYS                                    
                  FYZ(I,J) = FYZ(I,J) + FYZS                                    
                  FXZ(I,J) = FXZ(I,J) + FXZS                                    
               END IF                                                           
 198        CONTINUE                                                            
 200     CONTINUE                                                               
                                                                                
C        NORTH POLE.                                                            
         J=JGMAX-1                                                              
         JP = J + 1                                                             
C        SUBTRACT STUFF LEAVING THE POLAR CAP.                                  
         F0SC = 0.                                                              
         F0SD = 0.                                                              
         TEMPC = 0.                                                             
         TEMPD = 0.                                                             
         DO 202  I=1,N                                                          
            IF (V(I,JP,K).GE.0.) THEN                                           
               EPS=ABS(DTOVDXJ*V(I,JP,K)*COSEDGE(JP)/COSCEN(J))                 
               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*(CY(I,J,K)+                            
     1          (EPS1-EPS)*CYY(I,J,K)))                                         
               F0SC = F0SC + F0S                                                
               FYS = EPSQ*(CY(I,J,K)+3.0*EPS1*CYY(I,J,K))                       
               FYYS = EPS*EPSQ*CYY(I,J,K)                                       
               FXS  = EPS*(CX(I,J,K)+EPS1*CXY(I,J,K))                           
               FZS  = EPS*(CZ(I,J,K)+EPS1*CYZ(I,J,K))                           
               FXXS = EPS *CXX(I,J,K)                                           
               FZZS = EPS *CZZ(I,J,K)                                           
               FXYS = EPSQ*CXY(I,J,K)                                           
               FYZS = EPSQ*CYZ(I,J,K)                                           
               FXZS = EPS *CXZ(I,J,K)                                           
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS * MASS2(I,J,K)                                        
               TEMPC = TEMPC + TEMP                                             
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
               MASS2(I,J,K) = MASS2(I,J,K) -TEMP                                
                                                                                
C              READJUST MOMENTS REMAINING IN BOX.                               
               C(I,J,K) = C(I,J,K) - F0S                                        
               CY(I,J,K)=EPS1Q*(CY(I,J,K)-3.0*EPS*CYY(I,J,K))                   
               CYY(I,J,K)=EPS1*EPS1Q*CYY(I,J,K)                                 
               CX(I,J,K) = CX(I,J,K)-FXS                                        
               CZ (I,J,K) = CZ(I,J,K)-FZS                                       
               CXX(I,J,K) = CXX(I,J,K)-FXXS                                     
               CZZ(I,J,K) = CZZ(I,J,K)-FZZS                                     
               CXY(I,J,K) = EPS1Q*CXY(I,J,K)                                    
               CYZ(I,J,K) = EPS1Q*CYZ(I,J,K)                                    
               CXZ(I,J,K) = CXZ(I,J,K)-FXZS                                     
                                                                                
C              SET TEMPORARY MOMENTS                                            
               F0(I,J) = F0(I,J) + F0S                                          
               FY(I,J) = FY(I,J) + FYS                                          
               FYY(I,J) = FYY(I,J) + FYYS                                       
               FX(I,J) = FX(I,J) + FXS                                          
               FXX(I,J) = FXX(I,J) + FXXS                                       
               FZ(I,J) = FZ(I,J) + FZS                                          
               FZZ(I,J) = FZZ(I,J) + FZZS                                       
               FXY(I,J) = FXY(I,J) + FXYS                                       
               FYZ(I,J) = FYZ(I,J) + FYZS                                       
               FXZ(I,J) = FXZ(I,J) + FXZS                                       
                                                                                
            ELSE                                                                
               F0S=  ABS(DT*V(I,JP,K)*                                          
     1          C(I,JP,K)/(A*DPOLE))                                            
               TEMP=  ABS(DT*V(I,JP,K)*                                         
     1          MASS2(I,JP,K)/(A*DPOLE))                                        
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               F0SD = F0SD + F0S                                                
               TEMPD = TEMPD + TEMP                                             
               F0CON(I,J) = F0CON(I,J) + TEMP                                   
C              SET TEMPORARY MOMENTS                                            
               F0(I,J)=F0(I,J)+ F0S                                             
             END IF                                                             
 202     CONTINUE                                                               
                                                                                
         DO 379 I=1,N                                                           
            C(I,JP,K) = C(I,JP,K) - F0SD                                        
            MASS2(I,JP,K) = MASS2(I,JP,K) - TEMPD                               
 379     CONTINUE                                                               
                                                                                
C        RECONSTRUCT STRATOSPHERE.                                              
         J=0                                                                    
         JP = J + 1                                                             
         DO 579 I=1,N                                                           
            C(I,J,K) = C(I,J,K) + F0SB                                          
            MASS2(I,J,K) = MASS2(I,J,K) + TEMPB                                 
 579     CONTINUE                                                               
         DO 210 I=1,N                                                           
            IF (V(I,JP,K).GE.0.) THEN                                           
               EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,JP,K))                  
               EPS1 = 1.-EPS                                                    
               TM=EPS*C(I,JP,K)-EPS1*F0(I,J)                                    
               C(I,JP,K) = C(I,JP,K) + F0(I,J)                                  
               CYY(I,JP,K)=EPS*EPS*FYY(I,J)+EPS1                               
     1          *EPS1*CYY(I,JP,K)+5.0*(EPS*EPS1                                 
     1          *(CY(I,JP,K)-FY(I,J))-(EPS1-EPS)*TM)                            
               CY(I,JP,K) = EPS*FY(I,J)+EPS1                                    
     1          *CY(I,JP,K)+ 3.0*TM                                             
               CX(I,JP,K)  = CX(I,JP,K) + FX(I,J)                               
               CXX(I,JP,K) = CXX(I,JP,K)+ FXX(I,J)                              
               CZ(I,JP,K)  = CZ(I,JP,K) + FZ (I,J)                              
               CZZ(I,JP,K) = CZZ(I,JP,K) + FZZ(I,J)                             
               CXY(I,JP,K) = EPS*FXY(I,J)+EPS1*CXY(I,JP,K)+                     
     1          3.*(-EPS1*FX(I,J)+EPS*CX(I,JP,K))                               
               CYZ(I,JP,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,JP,K)+                     
     1          3.*(-EPS1*FZ(I,J)+EPS*CZ(I,JP,K))                               
               CXZ(I,JP,K) = CXZ(I,JP,K) + FXZ(I,J)                             
               MASS2(I,JP,K) = MASS2(I,JP,K) + F0CON(I,J)                       
            END IF                                                              
 210     CONTINUE                                                               
                                                                                
         DO 240 J=1,JGMAX-2                                                     
            JP = J + 1                                                          
            DO 220 I=1,N                                                        
               IF (V(I,JP,K).GE.0.) THEN                                        
                  EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,JP,K))               
                  EPS1 = 1.-EPS                                                 
                  TM=EPS*C(I,JP,K)-EPS1*F0(I,J)                                 
                  C(I,JP,K) = C(I,JP,K) + F0(I,J)                               
                  CYY(I,JP,K)=EPS*EPS*FYY(I,J)+EPS1
     1             *EPS1*CYY(I,JP,K)+5.0*(EPS*EPS1                              
     1             *(CY(I,JP,K)-FY(I,J))-(EPS1-EPS)*TM)                         
                  CY(I,JP,K) = EPS*FY(I,J)+EPS1                                 
     1             *CY(I,JP,K)+ 3.0*TM                                          
                  CX(I,JP,K)  = CX(I,JP,K) + FX(I,J)                            
                  CXX(I,JP,K) = CXX(I,JP,K)+ FXX(I,J)                           
                  CZ(I,JP,K)  = CZ(I,JP,K) + FZ (I,J)                           
                  CZZ(I,JP,K) = CZZ(I,JP,K) + FZZ(I,J)                          
                  CXY(I,JP,K) = EPS*FXY(I,J)+EPS1*CXY(I,JP,K)+                  
     1             3.*(-EPS1*FX(I,J)+EPS*CX(I,JP,K))                            
                  CYZ(I,JP,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,JP,K)+                  
     1             3.*(-EPS1*FZ(I,J)+EPS*CZ(I,JP,K))                            
                  CXZ(I,JP,K) = CXZ(I,JP,K) + FXZ(I,J)                          
                  MASS2(I,JP,K) = MASS2(I,JP,K) + F0CON(I,J)                    
                                                                                
               ELSE                                                             
                  EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,J,K))                
                  EPS1 = 1. - EPS                                               
                  TM = -EPS*C(I,J,K) + EPS1*F0(I,J)                             
                  C(I,J,K) = C(I,J,K) + F0(I,J)                                 
                  CYY(I,J,K) = EPS*EPS*FYY(I,J)+                                
     1             EPS1*EPS1*CYY(I,J,K) + 5.*                                   
     1             (EPS*EPS1*(-CY(I,J,K)+FY(I,J))+                              
     1             (EPS1-EPS)*TM)                                               
                  CY(I,J,K) = EPS*FY(I,J)+EPS1*CY(I,J,K)+3.*TM                  
                  CX(I,J,K)  = CX(I,J,K) + FX(I,J)                              
                  CXX(I,J,K) = CXX(I,J,K)+ FXX(I,J)                             
                  CZ(I,J,K)  = CZ(I,J,K) + FZ (I,J)                             
                  CZZ(I,J,K) = CZZ(I,J,K) + FZZ(I,J)                            
                  CXY(I,J,K) = EPS*FXY(I,J)+EPS1*CXY(I,J,K)+                    
     1             3.*(EPS1*FX(I,J)-EPS*CX(I,J,K))                              
                  CYZ(I,J,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,J,K)+                    
     1             3.*(EPS1*FZ(I,J)-EPS*CZ(I,J,K))                              
                  CXZ(I,J,K) = CXZ(I,J,K) + FXZ(I,J)                            
                  MASS2(I,J,K) = MASS2(I,J,K) + F0CON(I,J)                      
               END IF                                                           
 220        CONTINUE                                                            
 299     FORMAT(2I4,(9E12.5))                                                   
 240     CONTINUE                                                               
                                                                                
         J=JGMAX-1                                                              
         JP = J + 1                                                             
         DO 322 I=1,N                                                           
            C(I,JP,K) = C(I,JP,K) + F0SC                                        
            MASS2(I,JP,K) = MASS2(I,JP,K) + TEMPC                               
 322     CONTINUE                                                               
         DO 245 I=1,N                                                           
            IF (V(I,JP,K).LT.0.) THEN                                           
               EPS = F0CON(I,J) / (F0CON(I,J) + MASS2(I,J,K))                   
               EPS1 = 1. - EPS                                                  
               TM = -EPS*C(I,J,K) + EPS1*F0(I,J)                                
               C(I,J,K) = C(I,J,K) + F0(I,J)                                    
               CYY(I,J,K) = EPS*EPS*FYY(I,J)+                                   
     1          EPS1*EPS1*CYY(I,J,K) + 5.*                                      
     1          (EPS*EPS1*(-CY(I,J,K)+FY(I,J))+                                 
     1          (EPS1-EPS)*TM)                                                  
               CY(I,J,K) = EPS*FY(I,J)+EPS1*CY(I,J,K)+3.*TM                     
               CX(I,J,K)  = CX(I,J,K) + FX(I,J)                                 
               CXX(I,J,K) = CXX(I,J,K)+ FXX(I,J)                                
               CZ(I,J,K)  = CZ(I,J,K) + FZ (I,J)                                
               CZZ(I,J,K) = CZZ(I,J,K) + FZZ(I,J)                               
               CXY(I,J,K) = EPS*FXY(I,J)+EPS1*CXY(I,J,K)+                       
     1          3.*(EPS1*FX(I,J)-EPS*CX(I,J,K))                                 
               CYZ(I,J,K) = EPS*FYZ(I,J)+EPS1*CYZ(I,J,K)+                       
     1          3.*(EPS1*FZ(I,J)-EPS*CZ(I,J,K))                                 
               CXZ(I,J,K) = CXZ(I,J,K) + FXZ(I,J)                               
               MASS2(I,J,K) = MASS2(I,J,K) + F0CON(I,J)                         
            END IF                                                              
 245     CONTINUE                                                               
 260  CONTINUE                                                                  
                                                                                
C     WELL MIXED POLAR BOXES.                                                   
      DO 275 K=1,KTOP                                                           
         DO 270 I=1,N                                                           
            CX(I,0    ,K) = 0.                                                  
            CX(I,JGMAX,K) = 0.                                                  
            CY(I,0    ,K) = 0.                                                  
            CY(I,JGMAX,K) = 0.                                                  
            CZ(I,0    ,K) = 0.                                                  
            CZ(I,JGMAX,K) = 0.                                                  
            CXX(I,0    ,K) = 0.                                                 
            CXX(I,JGMAX,K) = 0.                                                 
            CYY(I,0    ,K) = 0.                                                 
            CYY(I,JGMAX,K) = 0.                                                 
            CZZ(I,0    ,K) = 0.                                                 
            CZZ(I,JGMAX,K) = 0.                                                 
            CXY(I,0    ,K) = 0.                                                 
            CXY(I,JGMAX,K) = 0.                                                 
            CXZ(I,0    ,K) = 0.                                                 
            CXZ(I,JGMAX,K) = 0.                                                 
            CYZ(I,0    ,K) = 0.                                                 
            CYZ(I,JGMAX,K) = 0.                                                 
 270     CONTINUE                                                               
 275  CONTINUE                                                                  
                                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE ZADVP(C,CX,CXX,CY,CYY,CZ,CZZ,CXY,CXZ,CYZ,MASS,MASS2,           
     1 FD0,DT,SGSTRAT,SGTROP,SGBOT,CHMASS,BOTFLUX)                              
                                                                                
C     VERTICAL ADVECTION IS PERFORMED WITH PRATHER'S SCHEME.                    
C     FIRST VERSION 1/8/91   COMPLETED?                                         
                                                                                
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11,KTM1=KTOP-1)          
      PARAMETER (N1=n+0)                                                        
                                                                                
      REAL CX(N1,0:JGMAX,KTOP),CY(N1,0:JGMAX,KTOP),CZ(N1,0:JGMAX,KTOP),         
     1 CXX(N1,0:JGMAX,KTOP),CYY(N1,0:JGMAX,KTOP),CZZ(N1,0:JGMAX,KTOP),          
     1 CXY(N1,0:JGMAX,KTOP),CXZ(N1,0:JGMAX,KTOP),CYZ(N1,0:JGMAX,KTOP)           
      REAL C(N1,0:JGMAX,KTOP),MASS2(N1,0:JGMAX,0:KTOP)                          
      REAL SGSTRAT(N1,0:JGMAX,0:KSTRAT-1),SGTROP(N1,0:JGMAX,KTROP)              
      REAL SGBOT(N1,0:JGMAX),MASS(N1,0:JGMAX,0:KTOP),FD0(N1,0:JGMAX)            
      REAL BOTFLUX(0:JGMAX)                                                     
                                                                                
      COMMON/VERT/DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),SIGMA(KTOP),           
     1 PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)                   
                                                                                
      REAL                                                                      
     1  FX(0:JGMAX,0:KTM1),FY(0:JGMAX,0:KTM1),FZ(0:JGMAX,0:KTM1),               
     1  FXX(0:JGMAX,0:KTM1),FYY(0:JGMAX,0:KTM1),FZZ(0:JGMAX,0:KTM1),            
     2  FXY(0:JGMAX,0:KTM1),FXZ(0:JGMAX,0:KTM1),FYZ(0:JGMAX,0:KTM1),            
     2  F0(0:JGMAX,0:KTM1),F0CON(0:JGMAX,0:KTM1),FD0WTC(N,0:JGMAX)              
C                                                                               
      COMMON // WEIGHT(0:JGMAX)                                                 
                                                                                
      DO 10 J=1,JGMAX-1                                                         
         WEIGHT(J) = 1.                                                         
 10   CONTINUE                                                                  
      WEIGHT(0)     = 1. / N                                                    
      WEIGHT(JGMAX) = 1. / N                                                    
                                                                                
      DO 260  I=1,N                                                             
C        ZERO TEMPORARY MOMENTS.                                                
         DO 20 K=0,KTM1                                                         
         DO 20 J=0,JGMAX                                                        
            FZ(J,K)   = 0.                                                      
            FZZ(J,K)   = 0.                                                     
            FY(J,K)   = 0.                                                      
            FYY(J,K)   = 0.                                                     
            FX(J,K)   = 0.                                                      
            FXX(J,K)   = 0.                                                     
            FYZ(J,K)   = 0.                                                     
            FXZ(J,K)   = 0.                                                     
            FXY(J,K)   = 0.                                                     
            F0CON(J,K) = 0.                                                     
            F0(J,K)   = 0.                                                      
 20      CONTINUE                                                               
                                                                                
C        CONVERT MIXING RATIO TO MASS IN LAYER ZERO.                            
         DO 40 J=0,JGMAX                                                        
            FD0WTC(I,J) = FD0(I,J) * MASS(I,J,0)                                
 40      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 60 K=1,KTOP                                                         
         DO 60 J=0,JGMAX                                                        
            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)))                   
            CYZ(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CYZ(I,J,K)))                    
            CXZ(I,J,K) = MIN(SLPMAX,MAX(-SLPMAX,CXZ(I,J,K)))                    
 60      CONTINUE                                                               
C                                                                               
C        LOWEST LEVEL.                                                          
         K=0                                                                    
         KP = K + 1                                                             
         DO 70  J=0,JGMAX                                                       
            IF (SGBOT(I,J).LT.0.) THEN                                          
               EPS=ABS(SGBOT(I,J)*DT/DELTAS(1))                                 
               IF (EPS.GE.1.0) WRITE (6,*) 'K=0,J=',J,SGBOT(I,J),EPS            
               EPSQ = EPS * EPS                                                 
               EPS1 = 1.0 - EPS                                                 
               EPS1Q = EPS1 * EPS1                                              
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                      
               F0S = EPS* FD0WTC(I,J)                                           
               CHMASS = CHMASS + F0S * WEIGHT(J)                                
               BOTFLUX(J) = BOTFLUX(J) + F0S * WEIGHT(J)                        
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS * MASS2(I,J,K)                                        
               F0CON(J,K) = F0CON(J,K) + TEMP                                   
C              SET TEMPORARY MOMENTS                                            
               F0(J,K) = F0(J,K) + F0S                                          
                                                                                
            ELSE                                                                
               EPS=ABS(SGBOT(I,J)*DT/DELTAS(1))                                 
               IF (EPS.GE.1.) WRITE (6,*) ' EPS GE 1 ',EPS,I,J,K                
               EPSQ = EPS*EPS                                                   
               EPS1 = 1.0 - EPS                                                 
               EPS1Q = EPS1 * EPS1                                              
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                      
               F0S = EPS * (C(I,J,KP)-EPS1*(CZ(I,J,KP)-                         
     1          (EPS1-EPS)*CZZ(I,J,KP)))                                        
               CHMASS = CHMASS - F0S * WEIGHT(J)                                
               BOTFLUX(J) = BOTFLUX(J) - F0S * WEIGHT(J)                        
               FZS=EPSQ*(CZ(I,J,KP)-3.0*EPS1*CZZ(I,J,KP))                       
               FZZS=EPS*EPSQ*CZZ(I,J,KP)                                        
               FYS  = EPS*(CY(I,J,KP)-EPS1*CYZ(I,J,KP))                         
               FXS  = EPS*(CX(I,J,KP)-EPS1*CXZ(I,J,KP))                         
               FYYS = EPS *CYY(I,J,KP)                                          
               FXXS = EPS *CXX(I,J,KP)                                          
               FYZS = EPSQ*CYZ(I,J,KP)                                          
               FXZS = EPSQ*CXZ(I,J,KP)                                          
               FXYS = EPS *CXY(I,J,KP)                                          
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS*MASS2(I,J,KP)                                         
               F0CON(J,K) = F0CON(J,K) + TEMP                                   
               MASS2(I,J,KP) = MASS2(I,J,KP) - TEMP                             
                                                                                
C              READJUST REMAINING MOMENTS.                                      
               C(I,J,KP) = C(I,J,KP) - F0S                                      
               CZ(I,J,KP) = EPS1Q*(CZ(I,J,KP)+3.0*EPS                           
     1          *CZZ(I,J,KP))                                                   
               CZZ(I,J,KP) = EPS1*EPS1Q*CZZ(I,J,KP)                             
               CY(I,J,KP) = CY(I,J,KP)-FYS                                      
               CX(I,J,KP) = CX(I,J,KP)-FXS                                      
               CYY(I,J,KP) = CYY(I,J,KP)-FYYS                                   
               CXX(I,J,KP) = CXX(I,J,KP)-FXXS                                   
               CYZ(I,J,KP) = EPS1Q*CYZ(I,J,KP)                                  
               CXZ(I,J,KP) = EPS1Q*CXZ(I,J,KP)                                  
               CXY(I,J,KP) = CXY(I,J,KP)-FXYS                                   
                                                                                
C              SET TEMPORARY MOMENTS                                            
               F0(J,K) = F0(J,K) + F0S                                          
               FZ(J,K) = FZ(J,K) + FZS                                          
               FZZ(J,K) = FZZ(J,K) + FZZS                                       
               FY(J,K) = FY (J,K) + FYS                                         
               FYY(J,K) = FYY(J,K) + FYYS                                       
               FX(J,K) = FX(J,K) + FXS                                          
               FXX(J,K) = FXX(J,K) + FXXS                                       
               FYZ(J,K) = FYZ(J,K) + FYZS                                       
               FXZ(J,K) = FXZ(J,K) + FXZS                                       
               FXY(J,K) = FXY(J,K) + FXYS                                       
            END IF                                                              
 70      CONTINUE                                                               
                                                                                
C        BEGIN TROPOSPHERIC LOOP                                                
         DO 100 K=1,KTROP-1                                                     
C           WRITE (6,*) ' K = ',K                                               
            KP = K + 1                                                          
            DO 80  J=0,JGMAX                                                    
               IF (SGTROP(I,J,K).LT.0.) THEN                                    
                  EPS=ABS(SGTROP(I,J,K)*DT/DELTAS(K))                           
C                 IF (EPS.GT.1.) WRITE (6,*) K,J,EPS,SGTROP(I,J,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)                                    
                  FYS  = EPS*(CY(I,J,K)+EPS1*CYZ(I,J,K))                        
                  FXS  = EPS*(CX(I,J,K)+EPS1*CXZ(I,J,K))                        
                  FYYS = EPS *CYY(I,J,K)                                        
                  FXXS = EPS *CXX(I,J,K)                                        
                  FYZS = EPSQ*CYZ(I,J,K)                                        
                  FXZS = EPSQ*CXZ(I,J,K)                                        
                  FXYS = EPS *CXY(I,J,K)                                        
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS * MASS2(I,J,K)                                     
                  F0CON(J,K) = F0CON(J,K) + TEMP                                
                  MASS2(I,J,K) = MASS2(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)                              
                  CY(I,J,K) = CY(I,J,K)-FYS                                     
                  CX(I,J,K) = CX(I,J,K)-FXS                                     
                  CYY(I,J,K) = CYY(I,J,K)-FYYS                                  
                  CXX(I,J,K) = CXX(I,J,K)-FXXS                                  
                  CYZ(I,J,K) = EPS1Q*CYZ(I,J,K)                                 
                  CXZ(I,J,K) = EPS1Q*CXZ(I,J,K)                                 
                  CXY(I,J,K) = CXY(I,J,K)-FXYS                                  
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0(J,K) = F0(J,K) + F0S                                       
                  FZ(J,K) = FZ(J,K) + FZS                                       
                  FZZ(J,K) = FZZ(J,K) + FZZS                                    
                  FY(J,K) = FY(J,K) + FYS                                       
                  FYY(J,K) = FYY(J,K) + FYYS                                    
                  FX(J,K) = FX(J,K) + FXS                                       
                  FXX(J,K) = FXX(J,K) + FXXS                                    
                  FYZ(J,K) = FYZ(J,K) + FYZS                                    
                  FXZ(J,K) = FXZ(J,K) + FXZS                                    
                  FXY(J,K) = FXY(J,K) + FXYS                                    
                                                                                
               ELSE                                                             
                  EPS=ABS(SGTROP(I,J,K)*DT/DELTAS(KP))                          
C                 IF (EPS.GT.1.) WRITE (6,*) K,J,EPS,SGTROP(I,J,K)              
                  EPSQ = EPS*EPS                                                
                  EPS1 = 1.0 - EPS                                              
                  EPS1Q = EPS1 * EPS1                                           
C                 DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                   
                  F0S = EPS * (C(I,J,KP)-EPS1*(CZ(I,J,KP)-                      
     1             (EPS1-EPS)*CZZ(I,J,KP)))                                     
                  FZS=EPSQ*(CZ(I,J,KP)-3.0*EPS1*CZZ(I,J,KP))                    
                  FZZS=EPS*EPSQ*CZZ(I,J,KP)                                     
                  FYS  = EPS*(CY(I,J,KP)-EPS1*CYZ(I,J,KP))                      
                  FXS  = EPS*(CX(I,J,KP)-EPS1*CXZ(I,J,KP))                      
                  FYYS = EPS *CYY(I,J,KP)                                       
                  FXXS = EPS *CXX(I,J,KP)                                       
                  FYZS = EPSQ*CYZ(I,J,KP)                                       
                  FXZS = EPSQ*CXZ(I,J,KP)                                       
                  FXYS = EPS *CXY(I,J,KP)                                       
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS*MASS2(I,J,KP)                                      
                  F0CON(J,K) = F0CON(J,K) + TEMP                                
                  MASS2(I,J,KP) = MASS2(I,J,KP) - TEMP                          
                                                                                
C                 READJUST REMAINING MOMENTS.                                   
                  C(I,J,KP) = C(I,J,KP) - F0S                                   
                  CZ(I,J,KP) = EPS1Q*(CZ(I,J,KP)+3.0*EPS                        
     1             *CZZ(I,J,KP))                                                
                  CZZ(I,J,KP) = EPS1*EPS1Q*CZZ(I,J,KP)                          
                  CY(I,J,KP) = CY(I,J,KP)-FYS                                   
                  CX(I,J,KP) = CX(I,J,KP)-FXS                                   
                  CYY(I,J,KP) = CYY(I,J,KP)-FYYS                                
                  CXX(I,J,KP) = CXX(I,J,KP)-FXXS                                
                  CYZ(I,J,KP) = EPS1Q*CYZ(I,J,KP)                               
                  CXZ(I,J,KP) = EPS1Q*CXZ(I,J,KP)                               
                  CXY(I,J,KP) = CXY(I,J,KP)-FXYS                                
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0(J,K) = F0(J,K) + F0S                                       
                  FZ(J,K) = FZ(J,K) + FZS                                       
                  FZZ(J,K) = FZZ(J,K) + FZZS                                    
                  FY(J,K) = FY (J,K) + FYS                                      
                  FYY(J,K) = FYY(J,K) + FYYS                                    
                  FX(J,K) = FX(J,K) + FXS                                       
                  FXX(J,K) = FXX(J,K) + FXXS                                    
                  FYZ(J,K) = FYZ(J,K) + FYZS                                    
                  FXZ(J,K) = FXZ(J,K) + FXZS                                    
                  FXY(J,K) = FXY(J,K) + FXYS                                    
               END IF                                                           
 80         CONTINUE                                                            
 100     CONTINUE                                                               
                                                                                
C        INTERFACE SECTION.                                                     
         K=KTROP                                                                
         K2 = K-KTROP                                                           
         KP = K + 1                                                             
         DO 120 J=0,JGMAX                                                       
            IF (SGTROP(I,J,K).LT.0.) THEN                                       
               EPS=ABS(SGTROP(I,J,K)*DT/DELTAS(K))                              
C                 IF (EPS.GT.1.) WRITE (6,*) K,J,EPS,SGTROP(I,J,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)                                       
               FYS  = EPS*(CY(I,J,K)+EPS1*CYZ(I,J,K))                           
               FXS  = EPS*(CX(I,J,K)+EPS1*CXZ(I,J,K))                           
               FYYS = EPS *CYY(I,J,K)                                           
               FXXS = EPS *CXX(I,J,K)                                           
               FYZS = EPSQ*CYZ(I,J,K)                                           
               FXZS = EPSQ*CXZ(I,J,K)                                           
               FXYS = EPS *CXY(I,J,K)                                           
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS * MASS2(I,J,K)                                        
               F0CON(J,K) = F0CON(J,K) + TEMP                                   
               MASS2(I,J,K) = MASS2(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)                                 
               CY(I,J,K) = CY(I,J,K)-FYS                                        
               CX(I,J,K) = CX(I,J,K)-FXS                                        
               CYY(I,J,K) = CYY(I,J,K)-FYYS                                     
               CXX(I,J,K) = CXX(I,J,K)-FXXS                                     
               CYZ(I,J,K) = EPS1Q*CYZ(I,J,K)                                    
               CXZ(I,J,K) = EPS1Q*CXZ(I,J,K)                                    
               CXY(I,J,K) = CXY(I,J,K)-FXYS                                     
                                                                                
C              SET TEMPORARY MOMENTS                                            
               F0(J,K) = F0(J,K) + F0S                                          
               FZ(J,K) = FZ(J,K) + FZS                                          
               FZZ(J,K) = FZZ(J,K) + FZZS                                       
               FY(J,K) = FY(J,K) + FYS                                          
               FYY(J,K) = FYY(J,K) + FYYS                                       
               FX(J,K) = FX(J,K) + FXS                                          
               FXX(J,K) = FXX(J,K) + FXXS                                       
               FYZ(J,K) = FYZ(J,K) + FYZS                                       
               FXZ(J,K) = FXZ(J,K) + FXZS                                       
               FXY(J,K) = FXY(J,K) + FXYS                                       
                                                                                
            ELSE                                                                
               EPS=ABS(SGSTRAT(I,J,K2)*DT/DELTAS(KP))                           
C                 IF (EPS.GT.1.) WRITE (6,*) K,J,EPS,SGTROP(I,J,K)              
               EPSQ = EPS*EPS                                                   
               EPS1 = 1.0 - EPS                                                 
               EPS1Q = EPS1 * EPS1                                              
C              DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                      
               F0S = EPS * (C(I,J,KP)-EPS1*(CZ(I,J,KP)-                         
     1          (EPS1-EPS)*CZZ(I,J,KP)))                                        
               FZS=EPSQ*(CZ(I,J,KP)-3.0*EPS1*CZZ(I,J,KP))                       
               FZZS=EPS*EPSQ*CZZ(I,J,KP)                                        
               FYS  = EPS*(CY(I,J,KP)-EPS1*CYZ(I,J,KP))                         
               FXS  = EPS*(CX(I,J,KP)-EPS1*CXZ(I,J,KP))                         
               FYYS = EPS *CYY(I,J,KP)                                          
               FXXS = EPS *CXX(I,J,KP)                                          
               FYZS = EPSQ*CYZ(I,J,KP)                                          
               FXZS = EPSQ*CXZ(I,J,KP)                                          
               FXYS = EPS *CXY(I,J,KP)                                          
                                                                                
C              DETERMINE BACKGROUND MIXING RATIO CHANGE                         
               TEMP = EPS*MASS2(I,J,KP)                                         
               F0CON(J,K) = F0CON(J,K) + TEMP                                   
               MASS2(I,J,KP) = MASS2(I,J,KP) - TEMP                             
                                                                                
C              READJUST REMAINING MOMENTS.                                      
               C(I,J,KP) = C(I,J,KP) - F0S                                      
               CZ(I,J,KP) = EPS1Q*(CZ(I,J,KP)+3.0*EPS                           
     1          *CZZ(I,J,KP))                                                   
               CZZ(I,J,KP) = EPS1*EPS1Q*CZZ(I,J,KP)                             
               CY(I,J,KP) = CY(I,J,KP)-FYS                                      
               CX(I,J,KP) = CX(I,J,KP)-FXS                                      
               CYY(I,J,KP) = CYY(I,J,KP)-FYYS                                   
               CXX(I,J,KP) = CXX(I,J,KP)-FXXS                                   
               CYZ(I,J,KP) = EPS1Q*CYZ(I,J,KP)                                  
               CXZ(I,J,KP) = EPS1Q*CXZ(I,J,KP)                                  
               CXY(I,J,KP) = CXY(I,J,KP)-FXYS                                   
                                                                                
C              SET TEMPORARY MOMENTS                                            
               F0(J,K) = F0(J,K) + F0S                                          
               FZ(J,K) = FZ(J,K) + FZS                                          
               FZZ(J,K) = FZZ(J,K) + FZZS                                       
               FY(J,K) = FY (J,K) + FYS                                         
               FYY(J,K) = FYY(J,K) + FYYS                                       
               FX(J,K) = FX(J,K) + FXS                                          
               FXX(J,K) = FXX(J,K) + FXXS                                       
               FYZ(J,K) = FYZ(J,K) + FYZS                                       
               FXZ(J,K) = FXZ(J,K) + FXZS                                       
               FXY(J,K) = FXY(J,K) + FXYS                                       
            END IF                                                              
 120     CONTINUE                                                               
                                                                                
C        BEGIN STRATOSPHERIC LOOP                                               
         DO 160 K=KTROP+1,KTOP-1                                                
C           WRITE (6,*) ' K = ', K                                              
            KP = K + 1                                                          
            K2 = K-KTROP                                                        
            DO 140 J=0,JGMAX                                                    
               IF (SGSTRAT(I,J,K2).LT.0.) THEN                                  
                  EPS=ABS(SGSTRAT(I,J,K2)*DT/DELTAS(K))                         
C                 IF (EPS.GT.1.) WRITE (6,*) K,J,EPS,SGTROP(I,J,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)                                    
                  FYS  = EPS*(CY(I,J,K)+EPS1*CYZ(I,J,K))                        
                  FXS  = EPS*(CX(I,J,K)+EPS1*CXZ(I,J,K))                        
                  FYYS = EPS *CYY(I,J,K)                                        
                  FXXS = EPS *CXX(I,J,K)                                        
                  FYZS = EPSQ*CYZ(I,J,K)                                        
                  FXZS = EPSQ*CXZ(I,J,K)                                        
                  FXYS = EPS *CXY(I,J,K)                                        
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS * MASS2(I,J,K)                                     
                  F0CON(J,K) = F0CON(J,K) + TEMP                                
                  MASS2(I,J,K) = MASS2(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)                              
                  CY(I,J,K) = CY(I,J,K)-FYS                                     
                  CX(I,J,K) = CX(I,J,K)-FXS                                     
                  CYY(I,J,K) = CYY(I,J,K)-FYYS                                  
                  CXX(I,J,K) = CXX(I,J,K)-FXXS                                  
                  CYZ(I,J,K) = EPS1Q*CYZ(I,J,K)                                 
                  CXZ(I,J,K) = EPS1Q*CXZ(I,J,K)                                 
                  CXY(I,J,K) = CXY(I,J,K)-FXYS                                  
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0(J,K) = F0(J,K) + F0S                                       
                  FZ(J,K) = FZ(J,K) + FZS                                       
                  FZZ(J,K) = FZZ(J,K) + FZZS                                    
                  FY(J,K) = FY(J,K) + FYS                                       
                  FYY(J,K) = FYY(J,K) + FYYS                                    
                  FX(J,K) = FX(J,K) + FXS                                       
                  FXX(J,K) = FXX(J,K) + FXXS                                    
                  FYZ(J,K) = FYZ(J,K) + FYZS                                    
                  FXZ(J,K) = FXZ(J,K) + FXZS                                    
                  FXY(J,K) = FXY(J,K) + FXYS                                    
                                                                                
               ELSE                                                             
                  EPS=ABS(SGSTRAT(I,J,K2)*DT/DELTAS(KP))                        
C                 IF (EPS.GT.1.) WRITE (6,*) K,J,EPS,SGTROP(I,J,K)              
                  EPSQ = EPS*EPS                                                
                  EPS1 = 1.0 - EPS                                              
                  EPS1Q = EPS1 * EPS1                                           
C                 DETERMINE MOMENTS FOR MASS LEAVING THE BOX.                   
                  F0S = EPS * (C(I,J,KP)-EPS1*(CZ(I,J,KP)-                      
     1             (EPS1-EPS)*CZZ(I,J,KP)))                                     
                  FZS=EPSQ*(CZ(I,J,KP)-3.0*EPS1*CZZ(I,J,KP))                    
                  FZZS=EPS*EPSQ*CZZ(I,J,KP)                                     
                  FYS  = EPS*(CY(I,J,KP)-EPS1*CYZ(I,J,KP))                      
                  FXS  = EPS*(CX(I,J,KP)-EPS1*CXZ(I,J,KP))                      
                  FYYS = EPS *CYY(I,J,KP)                                       
                  FXXS = EPS *CXX(I,J,KP)                                       
                  FYZS = EPSQ*CYZ(I,J,KP)                                       
                  FXZS = EPSQ*CXZ(I,J,KP)                                       
                  FXYS = EPS *CXY(I,J,KP)                                       
                                                                                
C                 DETERMINE BACKGROUND MIXING RATIO CHANGE                      
                  TEMP = EPS*MASS2(I,J,KP)                                      
                  F0CON(J,K) = F0CON(J,K) + TEMP                                
                  MASS2(I,J,KP) = MASS2(I,J,KP) - TEMP                          
                                                                                
C                 READJUST REMAINING MOMENTS.                                   
                  C(I,J,KP) = C(I,J,KP) - F0S                                   
                  CZ(I,J,KP) = EPS1Q*(CZ(I,J,KP)+3.0*EPS                        
     1             *CZZ(I,J,KP))                                                
                  CZZ(I,J,KP) = EPS1*EPS1Q*CZZ(I,J,KP)                          
                  CY(I,J,KP) = CY(I,J,KP)-FYS                                   
                  CX(I,J,KP) = CX(I,J,KP)-FXS                                   
                  CYY(I,J,KP) = CYY(I,J,KP)-FYYS                                
                  CXX(I,J,KP) = CXX(I,J,KP)-FXXS                                
                  CYZ(I,J,KP) = EPS1Q*CYZ(I,J,KP)                               
                  CXZ(I,J,KP) = EPS1Q*CXZ(I,J,KP)                               
                  CXY(I,J,KP) = CXY(I,J,KP)-FXYS                                
                                                                                
C                 SET TEMPORARY MOMENTS                                         
                  F0(J,K) = F0(J,K) + F0S                                       
                  FZ(J,K) = FZ(J,K) + FZS                                       
                  FZZ(J,K) = FZZ(J,K) + FZZS                                    
                  FY(J,K) = FY (J,K) + FYS                                      
                  FYY(J,K) = FYY(J,K) + FYYS                                    
                  FX(J,K) = FX(J,K) + FXS                                       
                  FXX(J,K) = FXX(J,K) + FXXS                                    
                  FYZ(J,K) = FYZ(J,K) + FYZS                                    
                  FXZ(J,K) = FXZ(J,K) + FXZS                                    
                  FXY(J,K) = FXY(J,K) + FXYS                                    
               END IF                                                           
 140        CONTINUE                                                            
 160     CONTINUE                                                               
                                                                                
C        RECONSTRUCT BOTTOM LEVEL.                                              
         K=0                                                                    
         KP = K + 1                                                             
         DO 170 J=0,JGMAX                                                       
            IF (SGBOT(I,J).LT.0.) THEN                                          
               EPS = F0CON(J,K) / (F0CON(J,K) + MASS2(I,J,KP))                  
               EPS1 = 1.-EPS                                                    
               TM=EPS*C(I,J,KP)-EPS1*F0(J,K)                                    
               C(I,J,KP) = C(I,J,KP) + F0(J,K)                                  
               CZZ(I,J,KP)=EPS*EPS*FZZ(J,K)+EPS1                               
     1          *EPS1*CZZ(I,J,KP)+5.0*(EPS*EPS1                                 
     1          *(CZ(I,J,KP)-FZ(J,K))-(EPS1-EPS)*TM)                            
               CZ(I,J,KP) = EPS*FZ(J,K)+EPS1                                    
     1          *CZ(I,J,KP)+ 3.0*TM                                             
               CY(I,J,KP)  = CY(I,J,KP) + FY(J,K)                               
               CYY(I,J,KP) = CYY(I,J,KP)+ FYY(J,K)                              
               CX(I,J,KP)  = CX(I,J,KP) + FZ (J,K)                              
               CXX(I,J,KP) = CXX(I,J,KP) + FXX(J,K)                             
               CYZ(I,J,KP) = EPS*FYZ(J,K)+EPS1*CYZ(I,J,KP)+                     
     1          3.*(-EPS1*FY(J,K)+EPS*CY(I,J,KP))                               
               CXZ(I,J,KP) = EPS*FXZ(J,K)+EPS1*CXZ(I,J,KP)+                     
     1          3.*(-EPS1*FX(J,K)+EPS*CX(I,J,KP))                               
               CXY(I,J,KP) = CXY(I,J,KP) + FXY(J,K)                             
               MASS2(I,J,KP) = MASS2(I,J,KP) + F0CON(J,K)                       
                                                                                
            END IF                                                              
 170     CONTINUE                                                               
                                                                                
C        RECONSTRUCT TROPOSPHERE.                                               
         DO 200 K=1,KTROP                                                       
            KP = K + 1                                                          
            DO 180 J=0,JGMAX                                                    
               IF (SGTROP(I,J,K).LT.0.) THEN                                    
                  EPS = F0CON(J,K) / (F0CON(J,K) + MASS2(I,J,KP))               
                  EPS1 = 1.-EPS                                                 
                  TM=EPS*C(I,J,KP)-EPS1*F0(J,K)                                 
                  C(I,J,KP) = C(I,J,KP) + F0(J,K)                               
                  CZZ(I,J,KP)=EPS*EPS*FZZ(J,K)+EPS1                            
     1             *EPS1*CZZ(I,J,KP)+5.0*(EPS*EPS1                              
     1             *(CZ(I,J,KP)-FZ(J,K))-(EPS1-EPS)*TM)                         
                  CZ(I,J,KP) = EPS*FZ(J,K)+EPS1                                 
     1             *CZ(I,J,KP)+ 3.0*TM                                          
                  CY(I,J,KP)  = CY(I,J,KP) + FY(J,K)                            
                  CYY(I,J,KP) = CYY(I,J,KP)+ FYY(J,K)                           
                  CX(I,J,KP)  = CX(I,J,KP) + FX (J,K)                           
                  CXX(I,J,KP) = CXX(I,J,KP) + FXX(J,K)                          
                  CYZ(I,J,KP) = EPS*FYZ(J,K)+EPS1*CYZ(I,J,KP)+                  
     1             3.*(-EPS1*FY(J,K)+EPS*CY(I,J,KP))                            
                  CXZ(I,J,KP) = EPS*FXZ(J,K)+EPS1*CXZ(I,J,KP)+                  
     1             3.*(-EPS1*FX(J,K)+EPS*CX(I,J,KP))                            
                  CXY(I,J,KP) = CXY(I,J,KP) + FXY(J,K)                          
                  MASS2(I,J,KP) = MASS2(I,J,KP) + F0CON(J,K)                    
                                                                                
               ELSE                                                             
                  EPS = F0CON(J,K) / (F0CON(J,K) + MASS2(I,J,K))                
                  EPS1 = 1. - EPS                                               
                  TM = -EPS*C(I,J,K) + EPS1*F0(J,K)                             
                  C(I,J,K) = C(I,J,K) + F0(J,K)                                 
                  CZZ(I,J,K) = EPS*EPS*FZZ(J,K)+                                
     1             EPS1*EPS1*CZZ(I,J,K) + 5.*                                   
     1             (EPS*EPS1*(-CZ(I,J,K)+FZ(J,K))+                              
     1             (EPS1-EPS)*TM)                                               
                  CZ(I,J,K) = EPS*FZ(J,K)+EPS1*CZ(I,J,K)+3.*TM                  
                  CY(I,J,K)  = CY(I,J,K) + FY(J,K)                              
                  CYY(I,J,K) = CYY(I,J,K)+ FYY(J,K)                             
                  CX(I,J,K)  = CX(I,J,K) + FX (J,K)                             
                  CXX(I,J,K) = CXX(I,J,K) + FXX(J,K)                            
                  CYZ(I,J,K) = EPS*FYZ(J,K)+EPS1*CYZ(I,J,K)+                    
     1             3.*(EPS1*FY(J,K)-EPS*CY(I,J,K))                              
                  CXZ(I,J,K) = EPS*FXZ(J,K)+EPS1*CXZ(I,J,K)+                    
     1             3.*(EPS1*FX(J,K)-EPS*CX(I,J,K))                              
                  CXY(I,J,K) = CXY(I,J,K) + FXY(J,K)                            
                  MASS2(I,J,K) = MASS2(I,J,K) + F0CON(J,K)                      
               END IF                                                           
 180        CONTINUE                                                            
 200     CONTINUE                                                               
                                                                                
C        RECONSTRUCT STRATOSPHERE.                                              
         DO 240 K=KTROP+1,KTM1                                                  
            KP = K + 1                                                          
            K2 = K-KTROP                                                        
            DO 220 J=0,JGMAX                                                    
               IF (SGSTRAT(I,J,K2).LT.0.) THEN                                  
                  EPS = F0CON(J,K) / (F0CON(J,K) + MASS2(I,J,KP))               
                  EPS1 = 1.-EPS                                                 
                  TM=EPS*C(I,J,KP)-EPS1*F0(J,K)                                 
                  C(I,J,KP) = C(I,J,KP) + F0(J,K)                               
                  CZZ(I,J,KP)=EPS*EPS*FZZ(J,K)+EPS1                            
     1             *EPS1*CZZ(I,J,KP)+5.0*(EPS*EPS1                              
     1             *(CZ(I,J,KP)-FZ(J,K))-(EPS1-EPS)*TM)                         
                  CZ(I,J,KP) = EPS*FZ(J,K)+EPS1                                 
     1             *CZ(I,J,KP)+ 3.0*TM                                          
                  CY(I,J,KP)  = CY(I,J,KP) + FY(J,K)                            
                  CYY(I,J,KP) = CYY(I,J,KP)+ FYY(J,K)                           
                  CX(I,J,KP)  = CX(I,J,KP) + FX (J,K)                           
                  CXX(I,J,KP) = CXX(I,J,KP) + FXX(J,K)                          
                  CYZ(I,J,KP) = EPS*FYZ(J,K)+EPS1*CYZ(I,J,KP)+                  
     1             3.*(-EPS1*FY(J,K)+EPS*CY(I,J,KP))                            
                  CXZ(I,J,KP) = EPS*FXZ(J,K)+EPS1*CXZ(I,J,KP)+                  
     1             3.*(-EPS1*FX(J,K)+EPS*CX(I,J,KP))                            
                  CXY(I,J,KP) = CXY(I,J,KP) + FXY(J,K)                          
                  MASS2(I,J,KP) = MASS2(I,J,KP) + F0CON(J,K)                    
                                                                                
               ELSE                                                             
                  EPS = F0CON(J,K) / (F0CON(J,K) + MASS2(I,J,K))                
                  EPS1 = 1. - EPS                                               
                  TM = -EPS*C(I,J,K) + EPS1*F0(J,K)                             
                  C(I,J,K) = C(I,J,K) + F0(J,K)                                 
                  CZZ(I,J,K) = EPS*EPS*FZZ(J,K)+                                
     1             EPS1*EPS1*CZZ(I,J,K) + 5.*                                   
     1             (EPS*EPS1*(-CZ(I,J,K)+FZ(J,K))+                              
     1             (EPS1-EPS)*TM)                                               
                  CZ(I,J,K) = EPS*FZ(J,K)+EPS1*CZ(I,J,K)+3.*TM                  
                  CY(I,J,K)  = CY(I,J,K) + FY(J,K)                              
                  CYY(I,J,K) = CYY(I,J,K)+ FYY(J,K)                             
                  CX(I,J,K)  = CX(I,J,K) + FX (J,K)                             
                  CXX(I,J,K) = CXX(I,J,K) + FXX(J,K)                            
                  CYZ(I,J,K) = EPS*FYZ(J,K)+EPS1*CYZ(I,J,K)+                    
     1             3.*(EPS1*FY(J,K)-EPS*CY(I,J,K))                              
                  CXZ(I,J,K) = EPS*FXZ(J,K)+EPS1*CXZ(I,J,K)+                    
     1             3.*(EPS1*FX(J,K)-EPS*CX(I,J,K))                              
                  CXY(I,J,K) = CXY(I,J,K) + FXY(J,K)                            
                  MASS2(I,J,K) = MASS2(I,J,K) + F0CON(J,K)                      
               END IF                                                           
 220        CONTINUE                                                            
 240     CONTINUE                                                               
 260  CONTINUE                                                                  
                                                                                
      RETURN                                                                    
      END                                                                       
