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



C     PROGRAM TRAN19                                                            
C                                                                               
      IMPLICIT REAL (A-H,O-Z)                                                   
                                                                                
C     DRIVER PROGRAM TO SOLVE TRACER CONTINUITY EQN. IN A SIGMA-                
C     PRESSURE COORDINATE SYSTEM. ADVECTION IS PERFORMED WITH                   
C     PRATHER'S SCHEME.                                                         
                                                                                
C     FIRST VERSION 1/5/91 (MODIFIED VERSION OF MAIN2 F)                        
C     LATEST VERSION 3/26/91                                                    
C                                                                               
      PARAMETER(N=144,JGMAX=88,NCON=1,JGMAX1=JGMAX+1)                           
      PARAMETER (KSTRAT=11,KTROP=8,KTOP=KSTRAT+KTROP,NTAPES=61)                 
      PARAMETER (ILMM=N,IJMM=JGMAX)
C
C     PHYSICAL SPACE ARRAYS ARE DEFINED                                         
      REAL CONST(N,0:JGMAX,KTOP,NCON),F1(N,0:JGMAX,0:KTOP)                      
      REAL F3(N,0:JGMAX,KTOP),CHMASS(NCON),CHMPSF(NCON)                         
      REAL FD0(N,0:JGMAX,NCON),DSH(NCON),DNH(NCON),DTOT(NCON)                   
      REAL D0SH(NCON),D0NH(NCON),D0TOT(NCON)                                    
                                                                                
C     DIAGNOSTIC VARIABLES FOR CFC RUN                                          
      REAL NSFLUX(0:JGMAX,KTOP,NCON),VRFLUX(0:JGMAX,KTOP,NCON),                 
     1 CFCSOURC(0:JGMAX,NCON),CFCLOSS(0:JGMAX,NCON)                             
      REAL BOTFLUX(0:JGMAX,NCON)                                                
      REAL CX(N,0:JGMAX,KTOP,NCON),CXX(N,0:JGMAX,KTOP,NCON),                    
     1     CY(N,0:JGMAX,KTOP,NCON),CYY(N,0:JGMAX,KTOP,NCON),                    
     1     CZ(N,0:JGMAX,KTOP,NCON),CZZ(N,0:JGMAX,KTOP,NCON),                    
     1     CXY(N,0:JGMAX,KTOP,NCON),CXZ(N,0:JGMAX,KTOP,NCON),                   
     1     CYZ(N,0:JGMAX,KTOP,NCON),MASS(N,0:JGMAX,0:KTOP),                     
     1     MASS2(N,0:JGMAX,0:KTOP),DPSF(N,0:JGMAX)                              
                                                                                
      REAL SCALE(NCON)                                                          
C                                                                               
      LOGICAL RSTAR,FIRST,WNDSAV,CHEM,DIFF,UPSTREAM,DIAG,DIURNAL              
      INTEGER NSTEPS(NTAPES),IST1(NTAPES),NOBS7(NTAPES),NOBS8(NTAPES)           
      INTEGER NOBS9(NTAPES),LUNCHM(NCON),WINDSOURCE,GRID 
      INTEGER ILUN(NTAPES)                                                      
C                                                                               
      REAL PSU(N,0:JGMAX),PSV(N,JGMAX),PSNEW(N,0:JGMAX)                         
      REAL PS1(N,0:JGMAX), PS2(N,0:JGMAX),PI1(N,0:JGMAX),                       
     ! U1  (N,0:JGMAX,KTOP), U2  (N,0:JGMAX,KTOP),                              
     ! V1  (N,JGMAX,KTOP), V2  (N,JGMAX,KTOP),SGBOT(N,0:JGMAX),                 
     ! SGTROP (N,0:JGMAX,KTROP),SGSTRAT(N,0:JGMAX,0:KSTRAT-1),                  
     ! KEL1(N,0:JGMAX,KTOP), KEL2(N,0:JGMAX,KTOP),                              
     !  DIVINT(N,0:JGMAX,2),PS1INV(N,0:JGMAX)                                   
                                                                                
      REAL COSCEN(0:JGMAX), DLATR(0:JGMAX),                                     
     ! TWOLAT(0:JGMAX), MCOR(N,0:JGMAX), PLON(N)                                
      REAL DLATR2(JGMAX),COSEDGE(JGMAX)                                         
                                                                                
C     CHEMICAL PRODUCTION AND LOSS ARRAYS                                       
      REAL PROD(0:JGMAX,KTOP,NCON),LOSS(0:JGMAX,KTOP,NCON)                      
      REAL PROD1(0:JGMAX,KTOP,NCON),LOSS1(0:JGMAX,KTOP,NCON)                  
      REAL PROD2(0:JGMAX,KTOP,NCON),LOSS2(0:JGMAX,KTOP,NCON)
      real theta(0:jgmax,ktop), theta1(0:jgmax,ktop) 
      real theta2(0:jgmax,ktop)                   
      REAL CL(0:JGMAX,KTOP),CLRAT(0:JGMAX,KTOP)                                 
C     REAL PROD2(N,0:JGMAX,NCON), LOSS2(N,0:JGMAX,KTOP,NCON)                    
C                                                                               
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),                    
     2 SIGMA(KTOP),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),                    
     3 SGINT(0:KTOP)                                                            
C                                                                               
      COMMON / THREE/HTSCALE,HINV,PI                                            
                                                                                
      COMMON/RDVAR/RSTAR,FIRST,ISAVE,ITAP1,IST1,NSTEPS,                         
     ! DT,NTSTRC,I1,B,CHEM,SCALE,WINDSOURCE,DIFF,DFCOEF,NPLOTS,DIURNAL,
     ! GRID,IREAD0                    
                                                                                
      COMMON/RDVAR2/DIAG                                                 
                                                                                
C     JANUARY AND FEBRUARY COMBINES NOBS                                        
C     MAIN BODY                                                                 
       DATA NOBS7/10,10,11,11,10,11,10,11,10,10,                                
     1            10,13,11,12, 1,11,12,11,11,11,                                
     1            11,12,11, 5, 4, 9,11,10, 9,10,                                
     1             6, 7,10,10,11,11, 9,12,13,13,                                
     1             8,15,15,16,16,16,16,16,16,16,                                
     1            16,18,15,14,13,14,13,13,14,13,                                
     1             9/                                                           
C     1989                                                                      
      DATA NOBS8/ 10,10,11,13,14,13,12,14,13,12,                                
     1            13,12,15,13,13,13,13,12,13,15,                                
     1            16,12,13,13,13,14,13,14,14, 2,                                
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                
     1             0/                                                           
                                                                                
      DATA NOBS9/ 20, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                
     1             0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                
     1             0/                                                           
                                                                                
      DATA ILUN/  51,52,53,54,55,56,57,58,59,60,                                
     1            61,62,63,64,65,66,67,68,69,70,                                
     1            71,72,73,74,75,76,77,78,79,80,                                
     1            81,82,83,84,85,86,87,88,89,90,                                
     1            91,92,93,94,94,96,97,98,99,60,                                
     1            61,62,63,64,65,66,67,68,69,70,                                
     1            71/ 

      DATA POLEDGE/-87./
                                                          
      TIN    = SECOND()                                                         
C                                                                               
      WRITE (6,*) 'BEGIN EXECUTION'                                             
C     INITIALIZE TIMING PARAMETERS                                              
      TIMIN  = 0.                                                               
      TIMCHM = 0.                                                               
      TIMDYN = 0.                                                               
      TIMDIA = 0.                                                               
      TIMX   = 0.                                                               
      TIMY   = 0.                                                               
      TIMZ   = 0.                                                               
      TIMDFF  = 0.                                                              
C     DEFINE VARIABLES AND OUTPUT EXPERIMENT HEADER TO UNIT10                   
      CALL SET19                                                                
C                                                                               
      WRITE (6,*) ' RESTART = ',RSTAR,'FIRST=',FIRST,' ISAVE=',ISAVE            
      WRITE (6,*) ' ITAP1=',ITAP1, ' DT=',DT, ' NTSTRC=',NTSTRC                 
      WRITE (6,*) ' IST1(K)'                                                    
      WRITE (6,982) (IST1(IK),IK=1,NTAPES)                                      
      WRITE (6,*) ' NSTEPS(IK)'                                                 
      WRITE (6,982) (NSTEPS(IK),IK=1,NTAPES)                                    
 982  FORMAT (1X,6I6)                                                           
                                                                                
      PI        = 4.0   * ATAN(1.)                                              
      PIINV     = 1./ PI                                                        
      HTSCALE   = 7000.                                                         
      HINV      = 1. / HTSCALE                                                  
      A         =  6.371220E6                                                   
                                                                                
C     MOLECULAR WEIGHTS OF F-11 AND AIR. FORMULA FOR F-11 IS CCL3F              
      DMWF11 = 12.011 + 3 * 35.453 + 18.998                                     
      DMW    = 28.97                                                            
C     MASS MIXING RATIO = VOLUME MIXING RATIO * RAT                             
      RAT    = DMWF11 / DMW                                                     
*********************SPECIAL FIX********                                        
      ADPHI      = PI * A / 90.                                                 
****************************************                                        
                                                                                
      ICHEM     = -999                                                          
                                                                                
C     LOGICAL UNIT NUMBERS ARE SPECIFIED.                                       
      DO 10 IC=1,NCON                                                           
         CHMASS(IC) = 0.                                                        
         LUNCHM(IC) = 48 + IC                                                   
 10   CONTINUE  
      LUNCHM2 = 47                                                           
      LNOON     =  33                                                           
      LMID      =  44                                                           
      WNDSAV    =  .FALSE.                                                      
      NDT       = 6                                                             
      NDT2      = 3                                                             
      TDT       =  DT / NTSTRC                                                  
      TDT2 = TDT                                                                
      TDT2I     = 1. / TDT2                                                     
      TDTPOL = TDT2 / NDT                                                       
      TDTPOL2= TDT2 / NDT2                                                      
                                                                                
C     DETERMINE TIME INTO INTEGRATION                                           
      TWIND1 = 0.                                                               
      IF (WINDSOURCE.EQ.1) THEN 
         IF (ITAP1.GT.I1) THEN                                                  
            DO 20 I=I1,ITAP1-1                                                  
               TWIND1 = TWIND1 + NOBS7(I) * DT                                  
 20         CONTINUE                                                            
         END IF                                                                 
      ELSE IF (WINDSOURCE.EQ.2) THEN                                           
         IF (ITAP1.GT.I1) THEN                                                  
            DO 30 I=I1,ITAP1-1                                                  
               TWIND1 = TWIND1 + NOBS8(I) * DT                                  
 30         CONTINUE                                                            
         END IF                                                                 
      ELSE IF (WINDSOURCE.GE.3) THEN                                           
         IF (ITAP1.GT.I1) THEN                                                  
            DO 26 I=I1,ITAP1-1                                                  
               TWIND1 = TWIND1 + NOBS9(I) * DT                                  
 26         CONTINUE                                                            
         END IF                                                                 
      END IF                                                                    
                                                                                
      IF (RSTAR) TWIND1 = TWIND1 - DT                                           
                                                                                
      IF (WINDSOURCE.EQ.1) THEN
C        -1 SINCE 12Z DEC 1 IS 1ST OBS OF ITAP1 (1979).                         
         TWIND1 = TWIND1 + (IST1(ITAP1)-1) * DT                                 
C        -3 SINCE 12Z JAN 14 IS 3RD OBS OF ITAP1 (1979).                        
C        TWIND1 = TWIND1 + (IST1(ITAP1)-3) * DT                                 
      ELSE IF (WINDSOURCE.EQ.2)	 THEN
C        -3 SINCE 00Z DEC 28 IS 3RD OBS OF ITAP1 (1989).                        
         TWIND1 = TWIND1 + (IST1(ITAP1)-3) * DT                                 
      ELSE IF (WINDSOURCE.GE.3)	 THEN
C        ADJUST AS NECESSARY
         TWIND1 = TWIND1 + (IST1(ITAP1)-1) * DT                                 
      END IF                                                                    
      TWIND2 = TWIND1                                                           
                                                                                
      DLON      =  360. / N                                                     
      DTIM2     =  (DLON / 720.) * 8.64E4                                       
      DTIM      =  (DLON / 360.) * 8.64E4                                       
      PLON(1)   =  DTIM2                                                        
      DO 40 IL = 2,N                                                            
         PLON(IL) = 86400. - DTIM2 - (IL-2) * DTIM                              
 40   CONTINUE                                                                  
                                                                                
C     VERTICAL DIFFERENCING ARRAYS ARE SET UP 
      IF (GRID.EQ.1) THEN                                  
         CALL SETVER 
      ELSE
         CALL VERGCM
      END IF                                                            
      WRITE (6,*) ' PINT = ', PINT, 'PTOP = ', PTOP                             
      WRITE (6,*) ' SIGMA,DSIG, AND DELTAS'                                     
      DO 50  K=1,KTOP                                                           
         WRITE (6,58) K, SIGMA(K), DSIG(K), DELTAS(K)                           
 58      FORMAT(1X,I4,3(1PE12.5))                                               
 50   CONTINUE                                                                  
                                                                                
C     GRID IS DEFINED (90. IS SPECIAL FIX).                                     
      D1 = -2.* POLEDGE/(JGMAX-1)
      D1O2 = D1 * 0.5                                               
      DTR = PI /180.                                                            
      ADPHI = D1 * DTR * A 
                                                    
      DLATR(0)         = -89.99 * DTR                                           
      COSCEN(0)        = COS(DLATR(0))                                          
      DO 60 J=1,JGMAX-1                                                         
         DLATR(J) = (POLEDGE+(J-1)*D1 + D1O2) * DTR 	
         COSCEN(J) = COS(DLATR(J))                                              
 60   CONTINUE                                                                  
      DLATR(JGMAX)     =  89.99 * DTR                                           
      COSCEN(JGMAX)    = COS(DLATR(JGMAX))                                      
                                                                                
      DO 70 J=1,JGMAX                                                           
         DLATR2(J) = (POLEDGE+(J-1)*D1) * DTR   
         COSEDGE(J) = COS(DLATR2(J))                                            
 70   CONTINUE                                                                  
      TWOLAT(0) = 90. + POLEDGE
      DO 80 IJ=1,JGMAX-1                                                        
         TWOLAT(IJ) = D1                                                        
 80   CONTINUE                                                                  
      TWOLAT(JGMAX) = 90. + POLEDGE
                                                                                
      DO  90 J=1,JGMAX-1                                                        
      DO  90 I=1,N                                                              
         MCOR(I,J) =  COSCEN(J) * TWOLAT(J) * DTR                         
 90   CONTINUE                                                                  
C                                                                               
      DPHIPOLE = 3. * PI / 180.                                                 
      DO 100 I=1,N                                                              
         MCOR(I,0)     = 0.5 * N * SIN((1.5*D1*DTR))**2                       
         MCOR(I,JGMAX) = 0.5 * N * SIN((1.5*D1*DTR))**2
 100  CONTINUE                                                                  
                                                                                
      WRITE (6,*) ' TRANSPORT LATITUDES '                                       
      DO 110 J=1,JGMAX                                                          
         WRITE (6,*) DLATR2(J), DLATR(J), COSEDGE(J), COSCEN(J)                 
 110  CONTINUE                                                                  
      WRITE (6,*)                                                               
      TIMIN = TIMIN + SECOND() - TIN                                            
                                                                                
C      TCH = SECOND()                                                           
C      IF (CHEM) THEN                                                           
*         COMMENT OUT AJRATE CALL UNLESS L4JAK.ZACONV.OUT CHEM USED             
C         CALL AJRATE(LUNCHM(1),LOSS(0,1,1),PROD(0,1,1))                        
C      END IF                                                                   
C      TIMCHM = TIMCHM + SECOND() - TCH                                         
                                                                                
      TDYN = SECOND()                                                           
C     INITIAL FIELDS ARE READ IN  AND TAPE IS POSITIONED                        
      IF (RSTAR) THEN                                                           
C        DIAGNOSTIC VARIABLES FOR CFC RUN READ IN.                              
         IF (DIAG) THEN                                                         
            CALL READDIAG(NSFLUX,VRFLUX,CFCSOURC,BOTFLUX,CFCLOSS,TIME,          
     1                    D0SH,D0NH,D0TOT)                                      
            WRITE (6,*) ' DIAGNOSTIC RESTART TIME ', TIME                       
         END IF                                                                 
C        ACCESS CONST FDS, TIME, AND TAPE POSITIONING PARAMETER                 
         CALL RESTARP2(CONST,CX,CXX,CY,CYY,CZ,CZZ,CXY,CXZ,CYZ,FD0,TIME,         
     1    IREAD,ISAVE)                                                          
                                                                                
C        ONE OBS IS READ BOTH IN THE 1ST RUN + IN THE SUBSEQUENT RESTART        
         IREAD = IREAD - 1                                                      
         WRITE (6,*) ' IREAD = ',IREAD                                          
           IF (WINDSOURCE.EQ.1) THEN                                            
              CALL RDOV19(IREAD,ILUN,LUN,ISK)                                   
              CALL RDSIM5(LUN,PS1,U1,V1,KEL1) 
           ELSE IF (WINDSOURCE.EQ.2) THEN
              CALL RDOV1989(IREAD,ILUN,LUN,ISK)                                 
              CALL RDSIM5(LUN,PS1,U1,V1,KEL1)  
           ELSE IF (WINDSOURCE.EQ.3) THEN
              CALL RDOVGCM(IREAD,ILUN,LUN,ISK)                                 
              CALL RDGCM(LUN,PS1,U1,V1,KEL1)  
           ELSE IF (WINDSOURCE.EQ.4) THEN
              CALL RDSCH(PS1,U1,V1,KEL1,COSEDGE)  
           ELSE IF (WINDSOURCE.EQ.5) THEN
              CALL WIND(U1,V1,COSCEN,B,DLATR,PS1,KEL1)  
           END IF                                                               
                                                                                
         IREAD = IREAD + 1                                                      
                                                                                
         TTEST = ABS(TIME - TWIND1)                                             
         IF (TTEST.GT.0.01) THEN                                                
            WRITE (6,*) ' TIMES ARE NOT EQUAL'                                  
            WRITE (6,*) ' TIME= ',TIME, ' TWIND1 = ',TWIND1                     
            AAA = 0.                                                            
            STOP                                                                
         END IF                                                                 
      ELSE                                                                      
C        ACCESS WIND, PRESSURE, AND TEMPERATURE FIELDS.                         
         IREAD = IREAD0                                                        
C        IREAD = 176  (IREAD FOR 12Z JAN 14 1979 START)                      
         IF (WINDSOURCE.EQ.1) THEN   
            CALL RDOV19(IREAD,ILUN,LUN,ISK)
            CALL RDSIM5(LUN,PS1,U1,V1,KEL1)
         ELSE IF (WINDSOURCE.EQ.2) THEN
            CALL RDOV1989(IREAD,ILUN,LUN,ISK)                                
            CALL RDSIM5(LUN,PS1,U1,V1,KEL1)
         ELSE IF (WINDSOURCE.EQ.3) THEN
            CALL RDOVGCM(IREAD,ILUN,LUN,ISK)                                 
            CALL RDGCM(LUN,PS1,U1,V1,KEL1)  
         ELSE IF (WINDSOURCE.EQ.4) THEN
            CALL RDSCH(PS1,U1,V1,KEL1,COSEDGE)  
            CALL BLOB(CONST,FD0(1,0,1)) 
         ELSE IF (WINDSOURCE.EQ.5) THEN
            CALL WIND(U1,V1,COSCEN,B,DLATR,PS1,KEL1)  
            CALL PATTERN(COSCEN,CONST,FD0(1,0,1)) 
         END IF
         IST1(ITAP1) = IST1(ITAP1) + 1
         IREAD = IREAD + 1 
                                                                                
C        DEFINE INITIAL CONDITION                                               
C        SET 1ST, 2ND, AND CROSS MOMENTS TO ZERO.                               
         DO 120 IC=1,NCON                                                       
         DO 120 K=1,KTOP                                                        
         DO 120 J=0,JGMAX                                                       
         DO 120 I=1,N                                                           
            CX   (I,J,K,IC) = 0.                                                
            CXX  (I,J,K,IC) = 0.                                                
            CY   (I,J,K,IC) = 0.                                                
            CYY  (I,J,K,IC) = 0.                                                
            CZ   (I,J,K,IC) = 0.                                                
            CZZ  (I,J,K,IC) = 0.                                                
            CXY  (I,J,K,IC) = 0.                                                
            CXZ  (I,J,K,IC) = 0.                                                
            CYZ  (I,J,K,IC) = 0.                                                
 120     CONTINUE   
C                                                            
C        CONSTITUENT IC'S ARE DEFINED.                                       
         LUNHNO = 45                                                         
         IC = 1                                                              
C        CALL TINIT(CONST(1,0,1,IC),LUNHNO,SCALE(IC),FD0(1,0,IC),IC)         
C        SPECIAL HCL IC IS DEFINED.                                          
C        DO 130 IC=2,NCON                                                    
C           DO 125 K=1,KTOP                                                  
C           DO 125 J=0,JGMAX                                                 
C           DO 125 I=1,N                                                     
C              FD0(I,J,IC) = FD0(I,J,1)                                      
C              CONST(I,J,K,IC) = CONST(I,J,K,1)                              
C125        CONTINUE                                                         
C130     CONTINUE                                                            
C        SPECIAL TRACER IC IS DEFINED.                                       
C        CALL INITTEST(CONST,FD0)                                            
C        1979 O3 AND HNO3 IC'S ARE DEFINED.                                  
C        CALL TINIT2(CONST,LUNHNO,16,FD0(1,0,1))                             
C        F-11 IC IS DEFINED FOR CFC RUN.                                     
******************************                    
         CALL INITO3(CONST(1,0,1,IC),LUNHNO,SCALE(IC),FD0(1,0,IC))
C        CALL INITCFC(CONST,FD0)                                 
C        DO 2133 IC=1,NCON                                                   
C        DO 2133 K=1,KTOP                                                    
C        DO 2133 J=0,JGMAX                                                   
C        DO 2133 I=1,N                                                       
C           FD0(I,J,IC) = 10.                                                
C           CONST(I,J,K,IC) = 10.                                            
C2133    CONTINUE                     
C                                      
         DO 133 IC=1,NCON                                                    
         DO 133 K=1,KTOP                                                     
         DO 133 J=0,JGMAX                                                    
            NSFLUX(J,K,IC) = 0.                                              
            VRFLUX(J,K,IC) = 0.                                              
 133     CONTINUE                                                            
         DO 234 IC=1,NCON                                                    
         DO 234 J=0,JGMAX                                                    
            CFCSOURC(J,IC) = 0.                                              
            BOTFLUX (J,IC) = 0.                                              
            CFCLOSS(J,IC)  = 0.                                              
 234     CONTINUE                                                            
                                                                                
         TIME   =  0.0                                                          
         IF (TIME.NE.TWIND1) THEN                                               
            WRITE (6,*) ' TIMES ARE NOT EQUAL'                                  
            STOP                                                                
         END IF                                                                 
      END IF                                                                    
                                                                                
C     TERM PROPORTIONAL TO BACKGROUND MASS IS CALC FOR EACH PT.                 
      CALL BKGNMASS(MCOR,PS1,MASS)                                              
C     THE MASS OF CONSTITUENTS ARE DETERMINED.                                  
      WRITE (6,314)  TWIND1/ 86400.                                             
 314  FORMAT (1X,' MASS OF CONSTS AT TIME ',F9.3)                               
      IF (RSTAR) THEN                                                           
         CALL MASS4(CONST,MASS,6,CHMASS,DSH,DNH,DTOT)                           
      ELSE                                                                      
         CALL MASS4(CONST,MASS,6,CHMASS,D0SH,D0NH,D0TOT)                        
         D0SH = D0SH * RAT                                                      
         D0NH = D0NH * RAT                                                      
         D0TOT= D0TOT* RAT                                                      
      END IF                                                                    
      TIMDYN = TIMDYN + SECOND() - TDYN                                         
C     BEGIN TIME STEPPING                                                       
      ICOUNT = -1                                                               
      DO 380 ITAPE = ITAP1,NTAPES                                               
      LUN = ILUN(ITAPE)                                                         
      DO 370 ISTDY = IST1(ITAPE),NSTEPS(ITAPE)                                  
C        READ IN WIND FIELD AT NEXT AVAIL TIME STEP                             
                                                                                
         TDYN = SECOND()                                                        
         IF (WINDSOURCE.EQ.1) THEN   
            CALL RDSIM5(LUN,PS2,U2,V2,KEL2)
         ELSE IF (WINDSOURCE.EQ.2) THEN
            CALL RDSIM5(LUN,PS2,U2,V2,KEL2)
         ELSE IF (WINDSOURCE.EQ.3) THEN
            CALL RDGCM(LUN,PS2,U2,V2,KEL2)  
         ELSE IF (WINDSOURCE.EQ.4) THEN
            CALL RDSCH(PS2,U2,V2,KEL2,COSEDGE)  
         ELSE IF (WINDSOURCE.EQ.5) THEN
            CALL WIND(U2,V2,COSCEN,B,DLATR,PS2,KEL2)  
         END IF
         TWIND2 = TWIND2 + DT                                                   
         IREAD = IREAD + 1                                                      
         TIMDYN = TIMDYN + SECOND() - TDYN                                      
C                                                                               
*************READ IN CHEMISTRY FIELDS****************                           
C        FOR ALL BUT THE N2O EXPERIMENT USE FOLLOWING CHEM CODE                 
          TCHM = SECOND()                                                       
          IF (CHEM) THEN                                                        
             DAY0 = -4.                                                         
             IF (WINDSOURCE.EQ.1) DAY0 = 0. 
C            CHEM SKIP PARAMETER FOR HCL,N2O, AND NOY EXPS.                     
             ICHEM2 = INT((DAY0 + TIME/86400.)/15.)                             
C            CHEM SKIP PARAMETER FOR O3                                         
C            WARNING INPUT O3 FILE DOES NOT CONTAIN HEADER!!!!                  
C            ICHEM2 = INT((DAY0+TIME/86400.)/30.)                               
C                                                                               
C            NOTICE NO CHEM READ IN FOR CFC RUN.                                
             IF (ICHEM2.NE.ICHEM) THEN
                CALL O3IN(PROD1,PROD2,LOSS1,LOSS2,THETA1,THETA2,
     1           ICHEM2,LUNCHM,LUNCHM2,DIURNAL)                                
C               CALL CHEMIN(PROD,LOSS,LUNCHM,ICHEM2)                            
C               CALL HCLIN(PROD,LOSS,CL,CLRAT,LUNCHM,ICHEM2)                    
C               CALL N2OIN (PROD(0,1,1),LOSS(0,1,1),LUNCHM(1),ICHEM2)           
C               CALL NOYIN (PROD(0,1,2),LOSS(0,1,2),LUNCHM(2),ICHEM2)           
                ICHEM = ICHEM2                                                  
             END IF                                                             
          END IF                                                                
         TIMCHM = TIMCHM + SECOND() - TCHM                                      
         DO 330 ISTCO = 1,NTSTRC  
C        SPECIAL SECTION FOR O3 CALCULATION
         TCHM=SECOND()
         IF (CHEM) CALL CHEMINT(PROD1,PROD2,PROD,LOSS1,LOSS2,LOSS,
     1    THETA1,THETA2,THETA,TIME) 
          TIMCHM = TIMCHM + SECOND() - TCHM 
                                                                                
         TDYN = SECOND()                                                        
         TIME2 = TIME + TDT                                                     
C        INTERPOLATE TO FIND DYNAMIC FIELDS AT CURRENT TIME.                    
         CALL WDINT4 (U1,U2,V1,V2,                                              
     !     PS1,PS2,KEL1,KEL2,TIME,TWIND1,TWIND2,PSNEW,TIME2)                    
         DO 135 J=0,JGMAX                                                       
         DO 135 I=1,N                                                           
            DPSF(I,J) = PSNEW(I,J) - PS1(I,J)                                   
 135     CONTINUE                                                               
                                                                                
         TWIND1 = TIME                                                          
         TIME   = TIME + TDT                                                    
         CALL SURFPS(PS1,PSU,PSV)                                               
                                                                                
C        DETERMINE "PI" ALONG THE BOX EDGES. AND 1/"PI" IN BOX CENT.            
         DO 140 IJ=1,IJMM                                                       
         DO 140 IL=1,ILMM                                                       
            PSV(IL,IJ) = PSV(IL,IJ) - PINT                                      
 140     CONTINUE                                                               
         DO 150 IJ=0,IJMM                                                       
         DO 150 IL=1,ILMM                                                       
            PSU(IL,IJ) = PSU(IL,IJ) - PINT                                      
            PS1INV(IL,IJ) = 1. / (PS1(IL,IJ)-PINT)                              
 150     CONTINUE                                                               
                                                                                
         TIMDYN = TIMDYN +SECOND() - TDYN                                       
                                                                                
         TDIA = SECOND()                                                        
                                                                                
C        THE CONSTITUENT FIELDS ARE SAVED.                                      
         MPLOT = MOD(ISTCO-1,NTSTRC)                                            
         IF ((MPLOT.EQ.0).AND.(FIRST)) THEN                                     
            ICOUNT = ICOUNT + 1                                                 
            LUNOUT = 10 + ICOUNT/NPLOTS                                         
            ILON = 0                                                            
            CALL HEADER (TWIND1,TDT,LUNOUT,ILON)                                
            DO 160 IC=1,NCON                                                    
               CALL PL3D (CONST(1,0,1,IC),PS1,IC,LUNOUT)                        
 160        CONTINUE                                                            
         ELSE IF (ISTCO.EQ.1) THEN                                              
            FIRST=.TRUE.                                                        
         END IF                                                                 
C                                                                               
C        THE LONG NEAREST LOCAL NOON + MIDNIGHT IS PLOTTED EACH HR.             
**       IF (MOD(TIME,3600.).EQ.0.0) THEN                                       
**          DO 170 IC=1,NCON                                                    
**             CALL PLNON2(CONST(1,0,1,IC),PS1,IC,TWIND1,TDT,PLON,              
**   !          LNOON,LMID)                                                     
**170       CONTINUE                                                            
**       END IF                                                                 
         TIMDIA = TIMDIA + SECOND() - TDIA                                      
                                                                                
         DO 260 IC=1,NCON                                                       
            TDYN = SECOND()                                                     
C           SET BACKGROUND MASS PRIME TO BACKGROUND MASS.                       
            DO 175 K=0,KTOP                                                     
            DO 175 J=0,JGMAX                                                    
            DO 175 I=1,N                                                        
               MASS2(I,J,K) = MASS(I,J,K)                                       
 175        CONTINUE                                                            
                                                                                
C           CONVERT FROM MIXING RATIO TO MASS.                                  
            DO 180 K=1,KTOP                                                     
            DO 180 J=0,JGMAX                                                    
            DO 180 I=1,N                                                        
               CONST(I,J,K,IC) = CONST(I,J,K,IC) * MASS(I,J,K)                  
 180        CONTINUE                                                            
                                                                                
C            WRITE (6,*) ' MASS BEFORE X ADVECTION '                           
C            CALL MASS5(CONST(1,0,1,IC),6)                                      
                                                                                
            TX   = SECOND()                                                     
C           SOUTHERN LATITUDE X ADVECTION IS CALCULATED.                        
            L1 = 1                                                              
            L2 = 2                                                              
 2004       FORMAT(10E12.5)                                                     
            DO 190 JCOUNT=1,NDT                                                 
               CALL XADVP(CONST(1,0,1,IC),CX(1,0,1,IC),CXX(1,0,1,IC),           
     1          CY (1,0,1,IC),CYY(1,0,1,IC),CZ (1,0,1,IC),CZZ(1,0,1,IC),        
     1          CXY(1,0,1,IC),CXZ(1,0,1,IC),CYZ(1,0,1,IC),MASS,                 
     1          MASS2,TDTPOL,U1,L1,L2,COSCEN,PS1INV,PSU)                        
 190        CONTINUE                                                            
            L1 = 3                                                              
            L2 = 9                                                              
            DO 195 JCOUNT=1,NDT2                                                
               CALL XADVP(CONST(1,0,1,IC),CX(1,0,1,IC),CXX(1,0,1,IC),           
     1          CY (1,0,1,IC),CYY(1,0,1,IC),CZ (1,0,1,IC),CZZ(1,0,1,IC),        
     1          CXY(1,0,1,IC),CXZ(1,0,1,IC),CYZ(1,0,1,IC),MASS,                 
     1          MASS2,TDTPOL2,U1,L1,L2,COSCEN,PS1INV,PSU)                       
 195        CONTINUE                                                            
C           MIDLATITUDE X ADVECTION IS CALCULATED.                              
            L1 = 10                                                             
            L2 = 75                                                             
            CALL XADVP(CONST(1,0,1,IC),CX(1,0,1,IC),CXX(1,0,1,IC),              
     1       CY (1,0,1,IC),CYY(1,0,1,IC),CZ (1,0,1,IC),CZZ(1,0,1,IC),           
     1       CXY(1,0,1,IC),CXZ(1,0,1,IC),CYZ(1,0,1,IC),MASS,                    
     1       MASS2,TDT2,U1,L1,L2,COSCEN,PS1INV,PSU)                             
                                                                                
C           NORTHERN LATITUDE X ADVECTION IS CALCULATED.                        
            L1 = 76                                                             
            L2 = 85                                                             
            DO 200 JCOUNT=1,NDT2                                                
               CALL XADVP(CONST(1,0,1,IC),CX(1,0,1,IC),CXX(1,0,1,IC),           
     1          CY (1,0,1,IC),CYY(1,0,1,IC),CZ (1,0,1,IC),CZZ(1,0,1,IC),        
     1          CXY(1,0,1,IC),CXZ(1,0,1,IC),CYZ(1,0,1,IC),MASS,                 
     1          MASS2,TDTPOL2,U1,L1,L2,COSCEN,PS1INV,PSU)                       
 200        CONTINUE                                                            
            L1 = 86                                                             
            L2 = 87                                                             
            DO 205 JCOUNT=1,NDT                                                 
               CALL XADVP(CONST(1,0,1,IC),CX(1,0,1,IC),CXX(1,0,1,IC),           
     1          CY (1,0,1,IC),CYY(1,0,1,IC),CZ (1,0,1,IC),CZZ(1,0,1,IC),        
     1          CXY(1,0,1,IC),CXZ(1,0,1,IC),CYZ(1,0,1,IC),MASS,                 
     1          MASS2,TDTPOL,U1,L1,L2,COSCEN,PS1INV,PSU)                        
 205        CONTINUE                                                            
            TIMX = TIMX + SECOND() - TX                                         
                                                                                
C           WRITE (6,*) ' MASS AFTER X ADVECTION '                              
C           CALL MASS5(CONST(1,0,1,IC),6)                                       
            TY = SECOND()                                                       
C           2. CALCULATE Y ADVECTION.                                           
                                                                                
C           MASS BEFORE Y ADVECTION IS NEEDED FOR CFC RUN.                      
            IF (DIAG) THEN                                                      
               DO 208 K=1,KTOP                                                  
               DO 208 J=0,JGMAX                                                 
               DO 208 I=1,N                                                     
                  F3(I,J,K) = CONST(I,J,K,IC)                                   
 208           CONTINUE                                                         
            END IF                                                              
C                                                                               
            CALL YADVP(CONST(1,0,1,IC),CX(1,0,1,IC),CXX(1,0,1,IC),              
     1       CY (1,0,1,IC),CYY(1,0,1,IC),CZ (1,0,1,IC),CZZ(1,0,1,IC),           
     1       CXY(1,0,1,IC),CXZ(1,0,1,IC),CYZ(1,0,1,IC),MASS,MASS2,              
     1       TDT2,V1,COSCEN,COSEDGE,ADPHI,PS1INV,PSV)                           
                                                                                
C           NORTH-SOUTH MASS FLUX IS CALCULATED FOR CFC RUN.                    
            IF (DIAG) THEN                                                      
               CALL DIAGNOSTIC(F3,CONST(1,0,1,IC),1,NSFLUX(0,1,IC),             
     1          VRFLUX(0,1,IC))                                                 
            END IF                                                              
                                                                                
C           WRITE (6,*) ' MASS AFTER Y ADVECTION '                              
C           CALL MASS5(CONST(1,0,1,IC),6)                                       
            TIMY = TIMY + SECOND() - TY                                         
            TIMDYN = TIMDYN + SECOND() - TDYN                                   
                                                                                
            IF (IC.EQ.1) THEN                                                   
               TDYN = SECOND()                                                  
               TZ   = SECOND()                                                  
C              CALCULATE VERTICAL VELOCITY RESIDUALLY.                          
***************************                                                     
C              CALCULATE NEG OF BACKGROUND MASS CHANGE DUE TO HOR ADV.          
               DO 2100 K=1,KTOP                                                 
               DO 2100 J=0,JGMAX                                                
               DO 2100 I=1,N                                                    
                  F3(I,J,K) = MASS(I,J,K) - MASS2(I,J,K)                        
 2100          CONTINUE                                                         
C               WRITE (6,*) ' RESIDUAL BEFORE PSF ADDED ON '                   
C               WRITE (6,988) (F3(1,JJ,1),JJ=0,JGMAX)                          
C 988           FORMAT (10E12.5)                                               
                                                                                
C              CALCULATE NEG OF BACKGROUND MASS CHNGE DUE TO PSF CHNGE.         
               DO 2105 K=1,KTROP                                                
               DO 2105 J=0,JGMAX                                                
               DO 2105 I=1,N                                                    
                  F3(I,J,K) = F3(I,J,K) + DPSF(I,J)*MASS2(I,J,K)*               
     1              PS1INV(I,J)                                                 
 2105          CONTINUE                                                         
C              WRITE (6,*) ' RESIDUAL AFTER PSF ADDED ON '                      
C              WRITE (6,988) (F3(1,JJ,1),JJ=0,JGMAX)                            
                                                                                
               CALL SGPRAP(MASS,MASS2,F3,SGBOT,SGTROP,SGSTRAT,TDT2,PS1)         
               TIMZ = TIMZ + SECOND() - TZ                                      
               TIMDYN = TIMDYN + SECOND() - TDYN                                
                                                                                
               TDIA = SECOND()                                                  
C              U,V,SG, AND TEMPERATURE FIELDS ARE SAVED.                        
               IF ((WNDSAV).AND.(ISTCO.EQ.1)) THEN                              
                  ILON = 0                                                      
                  LUNWND = 90 + (ICOUNT+1)/8                                    
                  CALL HEADER (TWIND1,TDT,LUNWND,ILON)                          
                  DO 210 IJ=0,IJMM                                              
                  DO 210 IL=1,ILMM                                              
                     F3(IL,IJ,1) = SGBOT(IL,IJ) * (PS1(IL,IJ)-PINT)             
 210              CONTINUE                                                      
                  DO 220  K=1,KTROP                                             
                     KP = K + 1                                                 
                     DO 220 IJ=0,IJMM                                           
                     DO 220 IL=1,ILMM                                           
                        F3(IL,IJ,KP)=SGTROP(IL,IJ,K)*(PS1(IL,IJ)-PINT)          
 220              CONTINUE                                                      
                  DO 230   K=1,KSTRAT-1                                         
                     KP = K + KTROP+1                                           
                     DO 230 IJ=0,IJMM                                           
                     DO 230 IL=1,ILMM                                           
                        F3(IL,IJ,KP)=SGSTRAT(IL,IJ,K) * (PINT-PTOP)             
 230              CONTINUE                                                      
                  CALL PL3D(U1,PS1,1,LUNWND)                                    
                  CALL PLXD(V1,LUNWND)                                          
                  CALL PL3D(F3,PS1,2,LUNWND)                                    
               END IF                                                           
               TIMDIA = TIMDIA + SECOND() - TDIA                                
            END IF                                                              
                                                                                
C           ADJUST FOR SURFACE PRESSURE CHANGE                                  
             DO 240 K=1,KTROP                                                   
             DO 240 IJ=0,IJMM                                                   
             DO 240 IL=1,ILMM                                                   
                CONST(IL,IJ,K,IC)= CONST(IL,IJ,K,IC)-DPSF(IL,IJ)                
     1        * CONST(IL,IJ,K,IC) * PS1INV(IL,IJ)                               
  240        CONTINUE                                                           
                                                                                
C           CALCULATE CHANGE DUE TO VERTICAL ADV.                               
            TDYN = SECOND()                                                     
            TZ = SECOND()                                                       
C           MASS BEFORE VERTICAL ADVECTION IS NEEDED FOR CFC RUN.               
            IF (DIAG) THEN                                                      
               DO 235 K=1,KTOP                                                  
               DO 235 J=0,JGMAX                                                 
               DO 235 I=1,N                                                     
                  F3(I,J,K) = CONST(I,J,K,IC)                                   
 235           CONTINUE                                                         
            END IF                                                              
            CALL ZADVP(CONST(1,0,1,IC),CX(1,0,1,IC),CXX(1,0,1,IC),              
     1       CY (1,0,1,IC),CYY(1,0,1,IC),CZ (1,0,1,IC),CZZ(1,0,1,IC),           
     1       CXY(1,0,1,IC),CXZ(1,0,1,IC),CYZ(1,0,1,IC),MASS,                    
     1       MASS2,FD0(1,0,IC),TDT2,SGSTRAT,SGTROP,SGBOT,CHMASS(IC),            
     1       BOTFLUX(0,IC))                                                     
                                                                                
C            WRITE (6,*) ' BOTTOM FLUX TERMS ',CHMASS(IC)                        
C           WRITE (6,*) (BOTFLUX(J,IC),J=0,JGMAX)                               
                                                                                
C           VERTICAL FLUX IS CALCULATED FOR CFC RUN.                            
            IF (DIAG) THEN                                                      
               CALL DIAGNOSTIC(F3,CONST(1,0,1,IC),2,NSFLUX(0,1,IC),             
     1          VRFLUX(0,1,IC))                                                 
            END IF                                                              
            TIMZ = TIMZ + SECOND() - TZ                                         
****************************                                                    
C           WRITE (6,*) ' CHMASS(IC) ', IC,CHMASS(IC)                           
C           WRITE (6,*) ' MASS AFTER Z ADVECTION ',IC                           
C           CALL MASS5(CONST(1,0,1,IC),6)                                       
                                                                                
C           TERM PROPORTIONAL TO BACKGROUND MASS IS CALC FOR EACH PT.           
            CALL BKGNMASS(MCOR,PSNEW,MASS)                                      
C           CONVERT FROM MASS TO MIXING RATIO.                                  
            DO 250 K=1,KTOP                                                     
            DO 250 J=0,JGMAX                                                    
            DO 250 I=1,N                                                        
               CONST(I,J,K,IC) = CONST(I,J,K,IC) / MASS(I,J,K)                  
 250        CONTINUE                                                            
                                                                                
            TIMDYN = TIMDYN + SECOND() - TDYN                                   
 260     CONTINUE                                                               
                                                                                
C        CHEMISTRY IS PERFORMED.                                                
         TCH = SECOND()                                                         
         IF (CHEM) THEN                                                         
            DO 270 IC=1,NCON                                                    
               CALL NEGS(CONST(1,0,1,IC),IC,TIME)                               
 270        CONTINUE                                                            
            DO 280 IC=1,NCON                                                    
            DO 280 K=1,KTOP                                                     
            DO 280 J=0,JGMAX                                                    
            DO 280 I=1,N                                                        
               IF (CONST(I,J,K,IC).LT.0.0)                                      
     1          CONST(I,J,K,IC)=0.0                                             
 280        CONTINUE                                                            
                                                                                
C           1.SPECIAL F-11 CHEMISTRY FOR CFC RUN.                               
C               CALL SOURCCFC(COSCEN,PS1,PROD2,MASS,CFCSOURC,TDT)              
C               CALL LOSSCFC(MASS,LOSS2,COSCEN,CONST,CFCLOSS,TDT)              
C               CALL CHEM3DDO(CONST,PROD2,LOSS2,TDT,MASS,CHMASS)               
C           CALL HCLCHEM(CONST,KEL1,PROD,LOSS,CL,CLRAT,TDT)                     
C           CALL CHEMDO(CONST,PROD,LOSS,TDT)                        
            CALL CHEMO3(CONST,PROD,LOSS,TDT,CHEMPROD,CHEMLOSS,NDIAG,
     1       DIAG,KEL1,THETA,DIURNAL,IDARK,TIME)       
         END IF                                                                 
         TIMCHM = TIMCHM + SECOND() - TCH                                       
                                                                                
         TIMD = SECOND()                                                        
         IF (DIFF) THEN                                                         
            DO 320 IC=1,NCON                                                    
               CALL DIFFX(CONST(1,0,1,IC),COSCEN,TDT2,F1,PS1INV,DFCOEF)         
               DO 290 K=1,KTOP                                                  
               DO 290 J=0,JGMAX                                                 
               DO 290 I=1,N                                                     
                  CONST(I,J,K,IC) = CONST(I,J,K,IC) +                           
     1              F1(I,J,K)                                                   
 290           CONTINUE                                                         
               CALL DIFFY(CONST(1,0,1,IC),COSCEN,COSEDGE,TDT,ADPHI,F1,          
     1            PS1INV,DFCOEF)                                                
               DO 300 K=1,KTOP                                                  
               DO 300 J=0,JGMAX                                                 
               DO 300 I=1,N                                                     
                  CONST(I,J,K,IC) = CONST(I,J,K,IC) +                           
     1              F1(I,J,K)                                                   
 300           CONTINUE                                                         
               CALL DIFFX(CONST(1,0,1,IC),COSCEN,TDT2,F1,PS1INV,DFCOEF)         
               DO 310 K=1,KTOP                                                  
               DO 310 J=0,JGMAX                                                 
               DO 310 I=1,N                                                     
                  CONST(I,J,K,IC) = CONST(I,J,K,IC) +                           
     1              F1(I,J,K)                                                   
 310           CONTINUE                                                         
 320        CONTINUE                                                            
         END IF                                                                 
         TIMDFF = TIMDFF + SECOND() - TIMD                                      
 330  CONTINUE                                                                  
                                                                                
C     THE MASS OF CONSTITUENTS ARE DETERMINED.                                  
c      TDIA = SECOND()                                                        
c      WRITE (6,314)  TIME/ 86400.                                            
c      CALL MASS4(CONST,MASS,6,CHMASS,DSH,DNH,DTOT)                           
c      TIMDIA = TIMDIA + SECOND() - TDIA                                      
      TDYN = SECOND()                                                           
      CALL SAVERP2(CONST,CX,CXX,CY,CYY,CZ,CZZ,CXY,CXZ,CYZ,FD0,TIME,             
     1 IREAD)                                                                   
C     DIAGNOSTIC FIELDS ARE SAVED FOR CFC RUN.                                  
      IF (DIAG)                                                                 
     1  CALL SAVEDIAG(NSFLUX,VRFLUX,CFCSOURC,BOTFLUX,CFCLOSS,TIME,              
     1 D0SH,D0NH,D0TOT,DSH,DNH,DTOT)                                            
C                                                                               
C     THE WIND,TEMPERATURE, + PRESSURE ARRAYS ARE REARRANGED                    
      DO 340 K=1,KTOP                                                           
      DO 340 J=1,JGMAX                                                          
      DO 340 I=1,N                                                              
         V1(I,J,K) = V2(I,J,K)                                                  
 340  CONTINUE                                                                  
                                                                                
      DO 350 K=1,KTOP                                                           
      DO 350 J=0,JGMAX                                                          
      DO 350 I=1,N                                                              
         U1  (I,J,K) = U2  (I,J,K)                                              
         KEL1(I,J,K) = KEL2(I,J,K)                                              
 350  CONTINUE                                                                  
                                                                                
      DO 360 J=0,JGMAX                                                          
      DO 360 I=1,N                                                              
         PS1(I,J) = PS2(I,J)                                                    
 360  CONTINUE                                                                  
      TWIND1 = TWIND2                                                           
      TIMDYN = TIMDYN + SECOND() - TDYN                                         
 370  CONTINUE                                                                  
 380  CONTINUE                                                                  
C                                                                               
      TDIA = SECOND()                                                           
C     MASS OF CONSTITUENT IS DETERMINED.                                        
      WRITE (6,314)  TIME/86400.                                                
      CALL MASS4(CONST,MASS,6,CHMASS,DSH,DNH,DTOT)                              
*     DO 390 IC = 1, NCON                                                       
C        THE CONSTITUENT FIELD AT LOCAL NOON IS PLOTTED OUT                     
*        CALL PLNON2(CONST(1,0,1,IC),PS1,IC,TIME,TDT,PLON,LNOON,LMID,           
*    !    FILFOR,DIFF)                                                          
* 390 CONTINUE                                                                  
C                                                                               
      ICOUNT = ICOUNT + 1                                                       
      LUNOUT = 10 + ICOUNT/NPLOTS                                               
      ILON = 0                                                                  
      CALL HEADER(TIME,TDT,LUNOUT,ILON)                                         
      DO 400 IC= 1,NCON                                                         
         CALL PL3D (CONST(1,0,1,IC),PS1,IC,LUNOUT)                              
 400  CONTINUE                                                                  
      TIMDIA = TIMDIA + SECOND() - TDIA                                         
C                                                                               
      DO 410 IC=1,NCON                                                          
      DO 410 K=1,KTOP                                                           
         DMAX = -1.E15                                                          
         DMIN =  1.E15                                                          
         DO 405 J=0,JGMAX                                                       
         DO 405 I=1,N                                                           
            IF (CONST(I,J,K,IC).GT.DMAX) THEN                                   
               DMAX  = CONST(I,J,K,IC)                                          
               JMAX  = J                                                        
               IMAX  = I                                                        
            ELSE IF (CONST(I,J,K,IC).LT.DMIN) THEN                              
               DMIN  = CONST(I,J,K,IC)                                          
               JMIN  = J                                                        
               IMIN  = I                                                        
            END IF                                                              
 405     CONTINUE                                                               
         WRITE (6,618) IMAX,JMAX,K,DMAX                                         
         WRITE (6,619) IMIN,JMIN,K,DMIN                                         
 618     FORMAT(' LONG = ',I4,' LAT = ',I4,' HT = ',I4,                         
     1   ' MAX = ',1PE15.7)                                                     
 619     FORMAT(' LONG = ',I4,' LAT = ',I4,' HT = ',I4,                         
     1   ' MIN = ',1PE15.7)                                                     
 410  CONTINUE                                                                  
                                                                                
      TIMIN = TIMIN + SECOND() - TIN                                            
      TIMCHM = TIMCHM*100. / TIMIN                                              
      TIMDYN = TIMDYN*100. / TIMIN                                              
      TIMDIA = TIMDIA*100. / TIMIN                                              
      TIMX =   TIMX  *100. / TIMIN                                              
      TIMY =   TIMY  *100. / TIMIN                                              
      TIMZ =   TIMZ  *100. / TIMIN                                              
      TIMDFF = TIMDFF *100./ TIMIN                                              
      WRITE (6,*) ' TOTAL TIME OF JOB ',TIMIN                                   
      WRITE (6,211) TIMCHM                                                      
      WRITE (6,212) TIMDYN                                                      
      WRITE (6,214) TIMDIA                                                      
      WRITE (6,215) TIMX                                                        
      WRITE (6,216) TIMY                                                        
      WRITE (6,217) TIMZ                                                        
      WRITE (6,219) TIMDFF                                                      
 211  FORMAT(1X,' CHEMISTRY TIME',F7.2)                                         
 212  FORMAT(1X,' DYNAMICAL TIME',F7.2)                                         
 214  FORMAT(1X,' DIAGNOSTIC TIME',F7.2)                                        
 215  FORMAT(1X,' X ADV     TIME',F7.2)                                         
 216  FORMAT(1X,' Y ADV     TIME',F7.2)                                         
 217  FORMAT(1X,' Z ADV     TIME',F7.2)                                         
 219  FORMAT(1X,' DIFFUSION TIME',F7.2)                                         
                                                                                
      STOP                                                                      
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE BKGNMASS(MCOR,PS1,MASS)                                        
                                                                                
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11)                      
                                                                                
      REAL MCOR(N,0:JGMAX),MASS(N,0:JGMAX,0:KTOP)                               
      REAL PS1(N,0:JGMAX)                                                       
                                                                                
      COMMON / VERT/DELTAS(KTOP), DSIG(0:KTOP),DELZ(KSTRAT),SIGMA(KTOP),        
     1 PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)                   
                                                                                
      G = 9.8                                                                   
      DO 10 K=1,KTROP                                                           
      DO 10 J=0,JGMAX                                                           
      DO 10 I=1,N                                                               
         MASS(I,J,K) =100.*(PS1(I,J)-PINT)*(SGINT(K-1)-SGINT(K))/G              
 10   CONTINUE                                                                  
                                                                                
      DO 15 J=0,JGMAX                                                           
      DO 15 I=1,N                                                               
         MASS(I,J,0) = MASS(I,J,1)                                              
 15   CONTINUE                                                                  
                                                                                
      DO 30 K=KTROP+1,KTOP                                                      
         K2 = K - KTROP                                                         
         DO 20 J=0,JGMAX                                                        
         DO 20 I=1,N                                                            
             MASS(I,J,K) =  100. * (PRSLAY(K2-1) - PRSLAY(K2)) / G              
 20      CONTINUE                                                               
 30   CONTINUE                                                                  
                                                                                
      DO 50 K=0,KTOP                                                            
      DO 50 J=0,JGMAX                                                           
      DO 50 I=1,N                                                               
         MASS(I,J,K) = MCOR(I,J) * MASS(I,J,K)                                  
 50   CONTINUE                                                                  
                                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE SAVERP2(CON1,CX,CXX,CY,CYY,CZ,CZZ,CXY,CXZ,CYZ,                 
     1  FD0,TIME,IREAD)                                                         
      IMPLICIT REAL (A-H,O-Z)                                                   
C                                                                               
      PARAMETER (KTOP=19,NCON=1,N=144,JGMAX=88,JGMAX1=JGMAX+1)                  
C                                                                               
      REAL CON1(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON),FD00(1:NCON)           
      REAL CX  (N,0:JGMAX,KTOP,NCON),CXX(N,0:JGMAX,KTOP,NCON)                   
      REAL CY  (N,0:JGMAX,KTOP,NCON),CYY(N,0:JGMAX,KTOP,NCON)                   
      REAL CZ  (N,0:JGMAX,KTOP,NCON),CZZ(N,0:JGMAX,KTOP,NCON)                   
      REAL CXY (N,0:JGMAX,KTOP,NCON),CXZ(N,0:JGMAX,KTOP,NCON)                   
      REAL CYZ (N,0:JGMAX,KTOP,NCON)                                            
C                                                                               
      LOGICAL DAY                                                               
                                                                                
      REWIND 30                                                                 
      DAY = .FALSE.                                                             
      IF (MOD(IREAD-1,40).EQ.0) DAY = .TRUE.                                    
C                                                                               
C     SAVE MIXING RATIO AND MOMENTS                                             
      WRITE (30) CON1                                                           
      WRITE (30) CX                                                             
      WRITE (30) CXX                                                            
      WRITE (30) CY                                                             
      WRITE (30) CYY                                                            
      WRITE (30) CZ                                                             
      WRITE (30) CZZ                                                            
      WRITE (30) CXY                                                            
      WRITE (30) CXZ                                                            
      WRITE (30) CYZ                                                            
C     SAVE LEVEL ZERO MIXING RATIO.                                             
***********FIX FOR 1-D FD0 ARRAYS**************                                 
      DO 30 IC=1,NCON                                                           
         FD00(IC) = 0.                                                          
         DO 20 J=0,JGMAX                                                        
         DO 20 I = 1,N                                                          
            FD00(IC) = FD00(IC) + FD0(I,J,IC)                                   
 20      CONTINUE                                                               
         FD00(IC) = FD00(IC) / (N*(JGMAX+1))                                    
 30   CONTINUE                                                                  
         WRITE (30) FD00                                                        
***********************************************                                 
                                                                                
C     WRITE (30) FD0                                                            
C     SAVE TIME INTO INTEGRATION AND IREAD.                                     
      WRITE (30) TIME                                                           
      WRITE (30) IREAD                                                          
                                                                                
C     SAVE MIXING RATIO AND MOMENTS                                             
      IF (DAY) THEN                                                             
         WRITE (31) CON1                                                        
         WRITE (31) CX                                                          
         WRITE (31) CXX                                                         
         WRITE (31) CY                                                          
         WRITE (31) CYY                                                         
         WRITE (31) CZ                                                          
         WRITE (31) CZZ                                                         
         WRITE (31) CXY                                                         
         WRITE (31) CXZ                                                         
         WRITE (31) CYZ                                                         
C        SAVE LEVEL ZERO MIXING RATIO.                                          
****************************************************FIX                         
         WRITE (31) FD00                                                        
C        WRITE (31) FD0                                                         
C        SAVE TIME INTO INTEGRATION AND IREAD.                                  
         WRITE (31) TIME                                                        
         WRITE (31) IREAD                                                       
      END IF                                                                    
                                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE RESTARP2(CON1,CX,CXX,CY,CYY,CZ,CZZ,CXY,CXZ,CYZ,                
     1  FD0,TIME,IREAD,ISAVE)                                                   
                                                                                
      IMPLICIT REAL (A-H,O-Z)                                                   
C                                                                               
C     THE INTEGRATION IS RESTARTED BY ACCESSING THE VALUES AT T AND             
C     T-TDT OF THE CONST. THE TIME AND IREAD ARE ALSO READ IN.                  
C                                                                               
      PARAMETER (KTOP=19,NCON=1,N=144,JGMAX=88)
                                                                                
      REAL CON1(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON),FD00(1:NCON)           
      REAL CX  (N,0:JGMAX,KTOP,NCON),CXX(N,0:JGMAX,KTOP,NCON)                   
      REAL CY  (N,0:JGMAX,KTOP,NCON),CYY(N,0:JGMAX,KTOP,NCON)                   
      REAL CZ  (N,0:JGMAX,KTOP,NCON),CZZ(N,0:JGMAX,KTOP,NCON)                   
      REAL CXY (N,0:JGMAX,KTOP,NCON),CXZ(N,0:JGMAX,KTOP,NCON)                   
      REAL CYZ (N,0:JGMAX,KTOP,NCON)                                            
C                                                                               
C     MEAN + 9 MOMENTS +FD0 + TIME + IREAD = 13                                 
      MRECS = 13                                                                
                                                                                
      DO 20 IOBS = 1,ISAVE                                                      
      DO 20 ISKIP= 1,MRECS                                                      
         READ(40) AAA                                                           
 20   CONTINUE                                                                  
                                                                                
C     READ IN MIXING RATIO AND MOMENTS                                          
      READ (40) CON1                                                            
      READ (40) CX                                                              
      READ (40) CXX                                                             
      READ (40) CY                                                              
      READ (40) CYY                                                             
      READ (40) CZ                                                              
      READ (40) CZZ                                                             
      READ (40) CXY                                                             
      READ (40) CXZ                                                             
      READ (40) CYZ                                                             
C     READ IN LEVEL ZERO MIXING RATIO.                                          
C     READ (40) FD0                                                             
*************************************FIX                                        
      READ (40) FD00                                                            
      DO 40 IC=1,NCON                                                           
      DO 40 J=0,JGMAX                                                           
      DO 40 I=1,N                                                               
         FD0(I,J,IC) = FD00(IC)                                                 
 40   CONTINUE                                                                  
****************************************                                        
C     READ IN TIME INTO INTEGRATION AND IREAD.                                  
      READ (40) TIME                                                            
      READ (40) IREAD                                                           
                                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE RDSCH(PSF,UWIND,VWIND,KEL,COSLAT2)                             
      IMPLICIT REAL (A-H,O-Z)                                                   
C                                                                               
C     THE V FIELD IS SPECIFIED.  THE U FIELD IS SET TO ZERO.                    
                                                                                
      PARAMETER (N=144,JGMAX=88)
      PARAMETER (KSTRAT=11,KTROP=8,KTOP=KTROP+KSTRAT)                           
                                                                                
      REAL VWIND(N,  JGMAX,KTOP)                                                
      REAL UWIND(N,0:JGMAX,KTOP)                                                
      REAL HT(KTOP), HT2(KTOP), KEL(N,0:JGMAX,KTOP)                             
      REAL PSF(N,0:JGMAX), COSLAT2(JGMAX)                                       
                                                                                
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),                    
     1 SIG(KTOP),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)         
                                                                                
      COMMON/THREE/HTSCALE,HINV,PI                                              
                                                                                
      DO 10 IJ=0,JGMAX                                                          
      DO 10 IL=1,N                                                              
         PSF(IL,IJ) = 1000.                                                     
 10   CONTINUE                                                                  
                                                                                
      DO 20 K=1,KTROP                                                           
         HT(K) = -HTSCALE * LOG(.001 * (SIG(K)*(PSF(1,1)-PINT)+PINT))           
 20   CONTINUE                                                                  
C                                                                               
      DO 25 K=KTROP+1,KTOP                                                      
         HT(K) = -HTSCALE * LOG(.001* PRESS(K-KTROP))                           
 25   CONTINUE                                                                  
                                                                                
      HT2(1) = HT(1) / 2.                                                       
      DO 30 K=2,KTOP                                                            
         HT2(K) = (HT(K) + HT(K-1)) / 2.                                        
 30   CONTINUE                                                                  
                                                                                
C     ZT = HT(KTOP)                                                             
C     ADJUST ZT SO VELOCITY FDS OF 30 AND 19 LEV MODS MATCH AT LOW LEVS         
      ZT = 56666.4                                                              
      A  = 6.371E6                                                              
C                                                                               
C     SET THE U FIELD TO ZERO.                                                  
      DO 40 K=1,KTOP                                                            
         DO 35 IJ=0,JGMAX                                                       
         DO 35 IL=1,N                                                           
            UWIND(IL,IJ,K)  =   0.                                              
            KEL  (IL,IJ,K)  = 300.                                              
 35      CONTINUE                                                               
 40   CONTINUE                                                                  
C                                                                               
C     SPECIFY V  FIELD.                                                         
      DO 60 K=1,KTOP                                                            
      DO 60 J=1,JGMAX                                                           
      DO  60 IL=1,N                                                             
         VWIND(IL,J,K) = -10.*HTSCALE * COSLAT2(J) *                            
     !                    (PI/ZT*COS(PI*HT2(K)/ZT) - (1./HTSCALE) *             
     !                     SIN(PI*HT2(K) / ZT))                                 
 60   CONTINUE                                                                  
                                                                                
C     DO 62 K=1,KTOP                                                            
C        J=1                                                                    
C        DO 61 IL=1,N                                                           
C           VWIND(IL,J,K) = 0.                                                  
C61      CONTINUE                                                               
C62   CONTINUE                                                                  
C                                                                               
C     DO 64 K=1,KTOP                                                            
C        J=JGMAX                                                                
C        DO 63 IL=1,N                                                           
C           VWIND(IL,J,K) = 0.                                                  
C63      CONTINUE                                                               
C64   CONTINUE                                                                  
                                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE BLOB(CONST,FD0)                                                
      IMPLICIT REAL (A-H,O-Z)                                                   
C                                                                               
                                                                                
      PARAMETER (N=144,JGMAX=88,KTOP=19,NCON=1)                 
                                                                                
      REAL CONST(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON)                       
                                                                                
C     MULTIPLY FOLLOWING BY 2 FOR "DOUBLE RESOLUTION RUN"                       
      H0 = 7.5                                                                  
      K1 = 6                                                                    
      K2 = 9                                                                    
      TJ0 = 6.0 *2                                                              
      J1 = 4    *2                                                              
      J2 = 8    *2                                                              
                                                                                
      DO 20 IC=1,NCON                                                           
      DO 20 K=1,KTOP                                                            
      DO 20 IJ=0,JGMAX                                                          
      DO 20 IL=1,N                                                              
         CONST(IL,IJ,K,IC) = 0.                                                 
 20   CONTINUE                                                                  
                                                                                
      DO 40 IC=1,NCON                                                           
      DO 40 K = K1,K2                                                           
      DO 40 J = J1,J2                                                           
         RAD = SQRT(((J-TJ0)**2 + (2*(K-H0))**2))                               
         DD  = EXP(-RAD * RAD/16.)                                              
         DO 30 IL=1,N                                                           
            CONST(IL,J,K,IC) = 40. * DD                                         
 30      CONTINUE                                                               
 40   CONTINUE                                                                  
                                                                                
      DO 60 IC=1,NCON                                                           
      DO 60 J=0,JGMAX                                                           
      DO 60 I=1,N                                                               
         FD0(I,J,IC) = CONST(I,J,1,IC)                                          
 60   CONTINUE                                                                  
                                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE MASS4(CON,MASS,LUN,CHMASS,DSH,DNH,DTOT)                        
      IMPLICIT REAL (A-H,O-Z)                                                   
C                                                                               
      PARAMETER (N=144,JGMAX=88,KTROP=8,KSTRAT=11,KTOP=KTROP+KSTRAT,            
     ! NCON=1,JG2=JGMAX/2)                                                      
C                                                                               
      REAL CON(N,0:JGMAX,KTOP,NCON),DMAS(KTOP,NCON),DMTOTL(NCON)                
      REAL MASS(N,0:JGMAX,0:KTOP),CHMASS(NCON)                                  
      REAL DSH(NCON), DNH(NCON), DTOT(NCON)                                     
C                                                                               
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),                    
     2 SIG(KTOP), P(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)            
                                                                                
C     MAIN BODY                                                                 
C     MIXING RATIO IS CONVERTED TO  MASS AND SUMMED.                            
C                                                                               
      DO 160 IC = 1,NCON                                                        
         DTOT(IC) = 0.0                                                         
         DSH(IC)  = 0.0                                                         
         DNH(IC)  = 0.0                                                         
         DO 140 K=1,KTOP                                                        
            SUM = 0.                                                            
            SUMS=0.                                                             
            SUMN=0.                                                             
            DO 130 J=0,JG2-1                                                    
            FAC = 1.                                                            
            IF (J.EQ.0) FAC=1./N                                                
            DO 130 I=1,N                                                        
               SUM = SUM + MASS(I,J,K) * CON(I,J,K,IC) * FAC                    
               SUMS= SUMS+ MASS(I,J,K) * CON(I,J,K,IC) * FAC                    
 130        CONTINUE                                                            
            DO 135 J=JG2+1,JGMAX                                                
            FAC = 1.                                                            
            IF (J.EQ.JGMAX) FAC=1./N                                            
            DO 135 I=1,N                                                        
               SUM = SUM + MASS(I,J,K) * CON(I,J,K,IC) * FAC                    
               SUMN= SUMN+ MASS(I,J,K) * CON(I,J,K,IC) * FAC                    
 135        CONTINUE                                                            
            J = JG2                                                             
            FAC = 1.                                                            
            DO 138 I=1,N                                                        
               SUM = SUM + MASS(I,J,K) * CON(I,J,K,IC) * FAC                    
               SUMN= SUMN+ MASS(I,J,K) * CON(I,J,K,IC) * FAC * 0.5              
               SUMS= SUMS+ MASS(I,J,K) * CON(I,J,K,IC) * FAC * 0.5              
 138        CONTINUE                                                            
            DMAS(K,IC) = SUM                                                    
            DTOT(IC) = DTOT(IC) + SUM                                           
            DSH(IC) = DSH(IC) + SUMS                                            
            DNH(IC) = DNH(IC) + SUMN                                            
 140     CONTINUE                                                               
C                                                                               
         WRITE (LUN,11) IC,DTOT(IC),CHMASS(IC)                                  
         WRITE (LUN,22) (DMAS(K,IC),K=1,KTOP)                                   
         WRITE (6,*)                                                            
 160  CONTINUE                                                                  
C                                                                               
C      IF (NCON.EQ.12) THEN                                                     
C         DMTOTL(11) = 0.                                                       
C         DMTOTL(12) = 0.                                                       
C         DO 180 K=1,KTOP                                                       
C            DMAS(K,11) = DMAS(K,2)+DMAS(K,3)+2.*DMAS(K,4)+DMAS(K,6)            
C     !                + DMAS(K,10)                                             
C            DMAS(K,12) = DMAS(K,5)+ DMAS(K,6)+DMAS(K,7)+DMAS(K,8)              
C            DMTOTL(11) = DMTOTL(11) + DMAS(K,11)                               
C            DMTOTL(12) = DMTOTL(12) + DMAS(K,12)                               
C 180     CONTINUE                                                              
C         WRITE (LUN,12) DMTOTL(11)                                             
C         WRITE (LUN,22) (DMAS(K,11),K=1,KTOP)                                  
C         WRITE (LUN,13) DMTOTL(12)                                             
C         WRITE (LUN,22) (DMAS(K,12),K=1,KTOP)                                  
C      END IF                                                                   
                                                                                
      RETURN                                                                    
 11   FORMAT (1X,' TOTAL SUM OF CONST ',I3,' IN ATM',E20.10,' MASS ',           
     1 'CHANGE DUE TO NET SFC FLUX ',E20.10)                                    
C12   FORMAT(1X,' SUM OF NITROGEN SPECIES IN ATM ',E15.5)                       
C13   FORMAT(1X,' SUM OF CHLORINE SPECIES IN ATM ',E15.5)                       
 22   FORMAT (1X,' SUM BY LEVEL ',(7E15.7))                                     
      END                                                                       
C                                                                               
C                                                                               
C                                                                               
      SUBROUTINE MASS5(CON,LUN)                                                 
      IMPLICIT REAL (A-H,O-Z)                                                   
C                                                                               
      PARAMETER (N=144,JGMAX=88,KTROP=8,KSTRAT=11,KTOP=KTROP+KSTRAT)            
      PARAMETER (JG2 = JGMAX/2)                                                 
C                                                                               
      REAL CON(N,0:JGMAX,KTOP),DMAS(KTOP),DMTOTL,WEIGHT(0:JGMAX)                
C                                                                               
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),                    
     2 SIG(KTOP), P(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)            
                                                                                
      DO 10 J=1,JGMAX-1                                                         
         WEIGHT(J) = 1.                                                         
 10   CONTINUE                                                                  
      WEIGHT(0) = 1./N                                                          
      WEIGHT(JGMAX) = 1. / N                                                    
                                                                                
C     MAIN BODY                                                                 
C                                                                               
         DMTOTL = 0.0                                                           
         SUMSH = 0.                                                             
         SUMNH = 0.                                                             
         DO 140 K=1,KTOP                                                        
            SUM = 0.                                                            
            DO 130 J=0,JG2-1                                                    
            DO 130 I=1,N                                                        
               SUM   = SUM   +  CON(I,J,K) * WEIGHT(J)                          
               SUMSH = SUMSH +  CON(I,J,K) * WEIGHT(J)                          
 130        CONTINUE                                                            
                                                                                
            DO 133 J=JG2+1,JGMAX                                                
            DO 133 I=1,N                                                        
               SUM = SUM +  CON(I,J,K) * WEIGHT(J)                              
               SUMNH = SUMNH +  CON(I,J,K) * WEIGHT(J)                          
 133        CONTINUE                                                            
                                                                                
            J=JG2                                                               
            DO 135 I=1,N                                                        
               SUMNH = SUMNH + CON(I,J,K) * WEIGHT(J) * 0.5                     
               SUMSH = SUMSH + CON(I,J,K) * WEIGHT(J) * 0.5                     
               SUM   = SUM   + CON(I,J,K) * WEIGHT(J)                           
 135        CONTINUE                                                            
            DMAS(K) = SUM                                                       
            DMTOTL = DMTOTL + SUM                                               
 140     CONTINUE                                                               
C               
         ic = 999                                                             
         WRITE (LUN,11) IC,DMTOTL,SUMSH,SUMNH                                   
         WRITE (LUN,22) (DMAS(K),K=1,KTOP)                                      
         WRITE (6,*)                                                            
C                                                                               
C      IF (NCON.EQ.12) THEN                                                     
C         DMTOTL(11) = 0.                                                       
C         DMTOTL(12) = 0.                                                       
C         DO 180 K=1,KTOP                                                       
C            DMAS(K,11) = DMAS(K,2)+DMAS(K,3)+2.*DMAS(K,4)+DMAS(K,6)            
C     !                + DMAS(K,10)                                             
C            DMAS(K,12) = DMAS(K,5)+ DMAS(K,6)+DMAS(K,7)+DMAS(K,8)              
C            DMTOTL(11) = DMTOTL(11) + DMAS(K,11)                               
C            DMTOTL(12) = DMTOTL(12) + DMAS(K,12)                               
C 180     CONTINUE                                                              
C         WRITE (LUN,12) DMTOTL(11)                                             
C         WRITE (LUN,22) (DMAS(K,11),K=1,KTOP)                                  
C         WRITE (LUN,13) DMTOTL(12)                                             
C         WRITE (LUN,22) (DMAS(K,12),K=1,KTOP)                                  
C      END IF                                                                   
                                                                                
      RETURN                                                                    
 11   FORMAT (1X,' TOT SM OF CONST ',I3,' IN ATM',E20.10,' SH SM ',             
     1 E20.10,' NH SM ',E20.10)                                                 
C12   FORMAT(1X,' SUM OF NITROGEN SPECIES IN ATM ',E15.5)                       
C13   FORMAT(1X,' SUM OF CHLORINE SPECIES IN ATM ',E15.5)                       
 22   FORMAT (1X,' SUM BY LEVEL ',(7E15.7))                                     
      END                                                                       
