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



      SUBROUTINE SAVER(CON1,FD0,TIME,IREAD)
      IMPLICIT REAL (A-H,O-Z)
C
C     NOTE: FD0 MADE 3-D ARRAY 3/6/91.
      PARAMETER (KTOP=19,NCON=1,N=144,JGMAX=88,JGMAX1=JGMAX+1)
      PARAMETER (NJG2=N*JGMAX1)
C
      REAL CON1(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON)
C
      LOGICAL DAY
      COMMON // ACYB2(N,0:JGMAX)
 
      REWIND 30
      DAY = .FALSE.
      IF (MOD(IREAD-1,40).EQ.0) DAY = .TRUE.
C
C     THE MIXING RATIO OF CONSTITUENT IS SAVED (LAYER BY LAYER) AND
C     CONSTITUENT  BY CONSTITUENT.
      DO 30 IC = 1, NCON
      DO 30 KK = 1, KTOP
         DO 20 IJ=0,JGMAX
         DO 20 IL=1,N
            ACYB2(IL,IJ) = CON1(IL,IJ,KK,IC)
  20     CONTINUE
         WRITE (30) ACYB2
         IF (DAY) WRITE (31) ACYB2
 30   CONTINUE
C
C     MIXING RATIO FOR LEVEL ZERO IS SAVED.
      WRITE (30) FD0
      IF (DAY) WRITE (31) FD0
 
C     THE TIME INTO THE INTEGRATION AND IREAD IS SAVED.
      WRITE (30) TIME
      IF (DAY) WRITE (31) TIME
      WRITE (30) IREAD
      IF (DAY) WRITE (31) IREAD
 
      RETURN
      END
C
C
C
      SUBROUTINE RESTAR (CON1,FD0,TIME,IREAD,ISAVE)
      IMPLICIT REAL (A-H,O-Z)
C
C     FD0 MADE 3-D ARRAY 3/6/91.
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,NJG2=N*(JGMAX+1))
 
      REAL CON1(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON)
C
      COMMON // ACYB2(N,0:JGMAX)
C
C     OBTAIN DESIRED RESTART DATA
      MRECS = NCON*KTOP  + 3
 
      DO 20 IOBS = 1,ISAVE
      DO 20 ISKIP= 1,MRECS
         READ(40) AAA
 20   CONTINUE
 
C     THE TRACER IS READ IN (LAYER BY LAYER) AND TRACER BY TRACER.
      DO 100 IC=1,NCON
      DO 100 KK = 1,KTOP
         READ (40) ACYB2
         DO 100 IJ=0,JGMAX
         DO 100 IL=1,N
            CON1(IL,IJ,KK,IC) = ACYB2(IL,IJ)
 100  CONTINUE
 
C     MIXING RATIO FOR LEVEL ZERO IS READ IN.
      READ (40) FD0
C     THE TIME INTO THE INTEGRATION AND IREAD IS READ IN
      READ (40) TIME
***************SPECIAL SECTION FOR READING (ETA TO MVS) RESTARTS
C      READ (40) TREAD
C      IREAD = TREAD+0.01
*****************END SPECIAL SECTION
      READ (40) IREAD
 
      RETURN
      END
C
C
C
      SUBROUTINE WDINT3 (U1,U2,V1,V2,
     ! PS1,PS2,KEL1,KEL2,TIME,TWIND1,TWIND2)
      IMPLICIT REAL (A-H,O-Z)
C
C     TIME INTERPOLATION IS PERFORMED ON THE INPUT FIELDS.
C     THE TIME INTERPOLATED FIELD IS RETURNED IN U1,V1,DIV1,PS1,+ KEL1
C
C
      PARAMETER (N=144,JGMAX=88,KTOP=19,NJG=N*JGMAX,NJG2=N*(JGMAX+1))
C
      REAL PS1(N,0:JGMAX), PS2(N,0:JGMAX),
     ! U1  (N,0:JGMAX,KTOP), U2  (N,0:JGMAX,KTOP),
     ! V1  (N,JGMAX,KTOP), V2  (N,JGMAX,KTOP),
     ! KEL1(N,0:JGMAX,KTOP), KEL2(N,0:JGMAX,KTOP)
 
      DELT = TWIND2 - TWIND1
      T2   = (TIME-TWIND1) / DELT
 
      DO 10 K= 1,KTOP
         DO 5 J=0,JGMAX
         DO 5 I=1,N
            U1  (I,J,K)= (1.-T2)* U1  (I,J,K)+T2 * U2 (I,J,K)
            KEL1(I,J,K)= (1.-T2)* KEL1(I,J,K)+T2* KEL2(I,J,K)
  5      CONTINUE
         DO 8 J=1,JGMAX
         DO 8 I=1,N
            V1  (I,J,K) = (1.-T2)*V1  (I,J,K)+T2 * V2  (I,J,K)
 8       CONTINUE
 10   CONTINUE
 
      DO 20 J=0,JGMAX
      DO 20 I=1,N
         PS1(I,J)  =      (1.-T2)*PS1 (I,J)+T2*PS2(I,J)
 20   CONTINUE
C
      RETURN
      END
C
C
C
      SUBROUTINE RDSCH(PSF,UWIND,VWIND,KEL,COSEDGE)
      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,NJG=N*(JGMAX+1))
      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), COSEDGE(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
 
      ZT = HT(KTOP)
C     ADJUST ZT SO VELOCITY FDS OF 30 AND 19 LEV MODS MATCH AT LOW LEVS
C     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. (set fac2 to 2 to speed up calculations)  
      fac2 = 2.  
      DO 60 K=1,KTOP
      DO 60 J=1,JGMAX
      DO  60 IL=1,N
         VWIND(IL,J,K) = -10.*HTSCALE * COSEDGE(J) * fac2 * 
     !                    (PI/ZT*COS(PI*HT2(K)/ZT) - (1./HTSCALE) *
     !                     SIN(PI*HT2(K) / ZT))
 60   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,NJG=N*(JGMAX+1))
 
      REAL CONST(N,0:JGMAX,KTOP,NCON), FD0(N,0:JGMAX)
 
C     MULTIPLY FOLLOWING BY 2 FOR "DOUBLE RESOLUTION RUN"
C     H0 = 7.5
C     K1 = 6
C     K2 = 9
C     TJ0 = 6.0 *2
C     J1 = 4    *2
C     J2 = 8    *2
      H0 = 9.0
      K1 = 7
      K2 = 11
      TJ0 = 10
      J1 = 8
      J2 = 12
 
      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 J=0,JGMAX
      DO 60 I=1,N
         FD0(I,J) = 0.
 60   CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE TINIT2(CONST,LUN1,LUN3,FD0)
      IMPLICIT REAL (A-H,O-Z)
C **********************************************************************
C *                                                                    *
C *   INITIAL CONDITION FOR O3 TRANSPORT. 3-DIMENSIONAL (1979 IC)      *
C **********************************************************************
C ----------------------------------------------------------------------
C     PART1 : DEALING WITH O3
C     THE O3 FIELD IS READ IN FROM O3
C
C     THE INPUT FILE CONTAINS  O3 MIXING RATIOS
C     (MOLECULES OF O3*10**6 / MOLECULES OF AIR) WITH LEVEL 1 AT
C     91.5 MB (16.74 KM), DELZ = 3.8 KM, AND LEVEL 11 AT 54.8 KM.
C     DATA IS ON R15 GAUSSIAN LATS.
C     THE DATA RUNS FROM 0.0 TO 360 BY 5 DEGREES (73 LON PTS).
C
C     THE FOLLOWING STEPS ARE EXECUTED.
C
C     1. THE DATA IS READ IN FROM O3XXX01.DATA
C     3. THE DATA IS INTERPOLATED LONGITUDINALLY ONTO N LONGITUDES.
C     4. DATA IS INTERP ONTO R30 GAUSS LATS.
C     5. THE TROPOSPHERIC HNO3 FD IS READ IN.
C
C
      PARAMETER (IMP=72,JMP=40,KTROP=8,KSTRAT=11,KTOP=KTROP+KSTRAT)
      PARAMETER (N=144,JGMAX=88,NCON=1,JGG=90)
      PARAMETER (NJG2=N*(JGMAX+1),NJGM1=N*(JGMAX-1))
C
      REAL CONST(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX)
C
      REAL DLAT(JMP), DGLAT(0:JGG),DLON1(IMP+1), DLON2(N)
      LOGICAL CYC
C
      COMMON //  C(IMP+1,JMP,KSTRAT),WORKA(N,JMP,KSTRAT), FD(IMP+1),
     ! HNO3TROP(KTROP),FD1(JMP,KTROP),FD2(0:JGG  ,KTROP),
     ! WORKB(N,0:JGG,KSTRAT),XIN(IMP+1),XOUT(N),YIN(JMP),YOUT(JGG+1)
C
      CHARACTER*80 A1
 
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),
     2 SIGMA(KTOP),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),
     3 SGINT(0:KTOP)
 
      COMMON/THREE/HTSCALE,HINV,PI
 
      DATA DLAT/-86.6,-82.2,-77.8,-73.3,-68.0,-64.4,-60.0,-55.5,-51.1,
     1 -46.7,-42.2,-37.8,-33.3,-28.9,-24.4,-20.0,-15.6,-11.1,-6.67,
     2 -2.22,2.22,6.67,11.1,15.6,20.0,24.4,28.9,33.3,37.8,42.2,46.7,
     3 51.1,55.5,60.0,64.4,68.0,73.3,77.8,82.2,86.6/
 
C     STEP1. READ IN O3 DATA.
C
C     THE VALUES ARE ON THE GAUSSIAN LATITUDES OF AN R15 RUN.
C     11 VERTICAL LEVELS. THE BOTTOM LAYER IS AT 91.5 MB
C     DELZ IS 3.8 KM
C
      READ (LUN1,19) A1
 19   FORMAT(A80)
      WRITE (6,19) A1
      DO 20 K=1,KSTRAT
      DO 20 J=1,JMP
         READ (LUN1,610) (FD(I),I=1,IMP+1)
         DO 10 IL=1,IMP
            C(IL,J,K) =  FD(IL) * 1.E-6
 10      CONTINUE
         C(IMP+1,J,K) = C(1,J,K)
   20 CONTINUE
  610 FORMAT (9F8.3)
C
C     STEP 3: E-W INTERPOLATION IS PERFORMED.
      D1 = 360. / IMP
      D2 = 360. / N
      DO 30 IL=1,IMP+1
         DLON1(IL) = 0. + (IL-1) * D1
 30   CONTINUE
      DO 40 IL=1,N
         DLON2(IL ) = 0. + (IL-1)*D2
 40   CONTINUE
 
      CALL INTERX3D(C,IMP+1,JMP,KSTRAT,WORKA,N,DLON1,DLON2,XIN,XOUT)
C
C     STEP 4: N-S INTERPOLATION ONTO OUTPUT LATS
      D2 = 180. / JGG
      DO 50 IJ=0,JGG
         DGLAT(IJ) = -90.+IJ*D2
 50   CONTINUE
      CALL INTERY3D(WORKA,N,JMP,KSTRAT,WORKB,JGG+1,DLAT,DGLAT,YIN,YOUT)
 
      DO 60 K=KTROP+1,KTOP
         DO 55 IL=1,N
            CONST(IL,0    ,K,1 ) = WORKB(IL,0  ,K-KTROP)
            CONST(IL,JGMAX,K,1) =  WORKB(IL,JGG,K-KTROP)
 55      CONTINUE
 
         DO 58 IJ=1,JGMAX-1
         IJ1 = IJ + 1
         DO 58 IL=1,N
            CONST(IL,IJ,K,1) =     WORKB(IL,IJ1,  K-KTROP)
 58      CONTINUE
 60   CONTINUE
 
C     STEP 5: TROPOSPHERIC O3 IS READ IN.
      READ (LUN3,19) A1
      DO 63 J=1,JMP
         READ(LUN3,44) JJJ, (HNO3TROP(K),K=1,KTROP)
         DO 62 K=1,KTROP
            FD1(J,K) = HNO3TROP(K)
  62     CONTINUE
  63  CONTINUE
  44  FORMAT(I4,8E9.2)
 
      CALL INTERX2D(FD1,JMP,KTROP,FD2,JGG+1,DLAT,DGLAT,YIN,YOUT)
      DO 80 K=1,KTROP
         DO 65 IL=1,N
            CONST(IL,0    ,K,1) = FD2(0  ,K)
            CONST(IL,JGMAX,K,1) = FD2(JGG,K)
 65      CONTINUE
         DO 70 J=1,JGMAX-1
         JP1 = J+1
         DO 70 IL=1,N
            CONST(IL,J,K,1) = FD2(JP1,K)
 70      CONTINUE
 80   CONTINUE
 
      DO 90 IJ=0,JGMAX
      DO 90 IL=1,N
         FD0(IL,IJ)  = CONST(IL,IJ,1,1)
 90   CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE PSFCHN(DIV,DIVINT,PSF,DPSF)
      IMPLICIT REAL(A-H,O-Z)
 
C     CREATION DATE  12/19/89
 
C     THE LOCAL DERIVATIVE OF SURFACE PRESSURE WRT TIME IS CALCULATED.
C
C     SEE EQNS 3.3-3.5 OF PAGE VI-8 OF GLAS GCM DOCUMENTATION
C
C     TROPOSPHERE   SIGMA = (P-PINT) / (PSFC-PINT)
C     STRATOSPHERIC SIGMA = (P-PINT) / (PINT-PTOP)
C
C     SELECTED INPUT FIELDS
C     SIG                SIGMA LEVELS
C     PSF                THE SURFACE PRESSURE
C     PINT               TOP OF TROPOSPHERE
C     PTOP               TOP OF ATMOSPHERE
C     DIVINT(1)          TROPOSPHERIC DIVERGENCE AT INTERFACE.
C     DIVINT(2)          STRATOSPHERIC DIVERGENCE AT INTERFACE.
C     OUTPUT FIELD
C     DPSF               LOCAL RATE OF CHANGE OF PSF WRT TIME.
C
      PARAMETER (ILMM=144,IJMM=88,IJMMP1=IJMM+1)
      PARAMETER (KTROP=8,KSTRAT=11,ILIJ=ILMM*IJMMP1,ILIJM=ILMM*IJMM)
      PARAMETER (IKMM=KSTRAT+KTROP,IKMMP1=IKMM+1,IKMMM1=IKMM-1)
 
      REAL PSF(ILMM,0:IJMM),
     1 DPSF(ILMM,0:IJMM),DIV(ILMM,0:IJMM,IKMM)
      REAL DIVINT(ILMM,0:IJMM,2)
C
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      COMMON/THREE/ HTSCALE, HINV,PI
 
      COMMON //  DIV2(ILMM,0:IJMM,0:IKMM),WK1(ILMM,0:IJMM)
C
C     DIVERGENCE AT SFC IS SET TO ZERO.
      DO 5 IJ=0,IJMM
      DO 5 IL=1,ILMM
         DIV2(IL,IJ,0) = 0.
 5    CONTINUE
C     DIVERGENCE IS INTERPOLATED ONTO MODEL LAYERS.
      DO 10 IK=1,KTROP-1
         IKP = IK +1
         DO 7 IJ=0,IJMM
         DO 7 IL=1,ILMM
            DIV2(IL,IJ,IK) = 0.5 * (DIV(IL,IJ,IK)+DIV(IL,IJ,IKP))
 7       CONTINUE
 10   CONTINUE
      DO 12 IJ=0,IJMM
      DO 12 IL=1,ILMM
         DIV2(IL,IJ,KTROP) = DIVINT(IL,IJ,2)
 12   CONTINUE
      DO 15 IK=KTROP+1,IKMMM1
         IKP = IK +1
         DO 13 IJ=0,IJMM
         DO 13 IL=1,ILMM
           DIV2(IL,IJ,IK)=0.5*(DIV(IL,IJ,IK)+DIV(IL,IJ,IKP))
 13      CONTINUE
 15   CONTINUE
C     DIVERGENCE AT TOP OF MODEL IS SET TO ZERO.
      DO 20 IJ=0,IJMM
      DO 20 IL=1,ILMM
         DIV2(IL,IJ,IKMM) = 0.
 20   CONTINUE
      PI1INV         = 1.0 / (PINT-PTOP)
      DO 223 K=KTROP,IKMMM1
         DO 200 IJ=0,IJMM
         DO 200 IL=1,ILMM
            DIV2(IL,IJ,K)=DIV2(IL,IJ,K)*PRSLAY(K-KTROP)*PI1INV*HINV
 200     CONTINUE
 223  CONTINUE
 
C     INTEGRATE DOWN TO TROPOSPHERE-STRATOSPHERE INTERFACE USING
C     MODIFIED TRAPEZOIDAL RULE OR SIMPSON'S RULE.
      IPTS=KSTRAT+1
      K3 = IKMM-IPTS+1
      CALL FINT(IPTS,DELZ(1)*(IPTS-1),DIV2(1,0,K3),WK1)
      DO 250 IJ=0,IJMM
      DO 250 IL=1,ILMM
         WK1(IL,IJ) = -WK1(IL,IJ)
 250  CONTINUE
C     LOCAL DER OF PSF WRT TIME IS DETERMINED.
      DO 260 IJ=0,IJMM
      DO 260 IL=1,ILMM
         DIV2(IL,IJ,KTROP) = DIVINT(IL,IJ,1)
 260  CONTINUE
      IPTS = 9
      K2 = KTROP-IPTS+1
      CALL FINT(IPTS,DELTAS(1)*(IPTS-1),DIV2(1,0,K2),DPSF(1,0))
      DO 280 IJ=0,IJMM
      DO 280 IL=1,ILMM
         DPSF(IL,IJ) = DPSF(IL,IJ)+WK1(IL,IJ)
 280  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE SIGADV(D,DT,SGTROP,SGSTRAT,SGBOT,DNEW,FD0)
      IMPLICIT REAL(A-H,O-Z)
 
C     6 JANUARY 1990 MODIFIED TO QUICKEN
C
C     VERTICAL FLUX TERM IN CONTINUITY EQN. IS CALCULATED.
 
C     D     MIXING RATIO OF CONSTITUENT * PI WHERE PI =
C           (PINT-PTOP) IN STRATOSPHERE AND (PSFC-PINT) IN TROPOSPHERE.
C     DT    TIME STEP
C     SG    SIGMA DOT IS DEFINED AT INTERFACES BETWEEN LEVELS AND IS
C           ZERO AT THE TOP AND BOTTOM INTERFACE. I.E. SG(0)=SG(IKMM)=0.
C     DNEW  FLUX TERM.
C
      PARAMETER (ILMM=144,IJMM=88,IJMMP1=IJMM+1)
      PARAMETER (KTROP=8,KSTRAT=11,IKMM=KTROP+KSTRAT)
      PARAMETER (IKMMM1=IKMM-1,ILIJ=ILMM*IJMMP1,KST1=KSTRAT-1)
 
      REAL SGTROP(ILMM,0:IJMM,KTROP),D(ILMM,0:IJMM,IKMM)
      REAL DNEW(ILMM,0:IJMM,IKMM),SGSTRAT(ILMM,0:IJMM,0:KST1)
      REAL SGBOT(ILMM,0:IJMM),FD0(ILMM,0:IJMM)
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     2 SIGMA(IKMM), PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),
     3 SGINT(0:IKMM)
 
      COMMON // TOP(ILMM,0:IJMM,IKMM), BOT(ILMM,0:IJMM,IKMM),
     1 FLUX1(ILMM,0:IJMM,0:KTROP),FLUX2(ILMM,0:IJMM,0:KST1),
     1 DELTA(ILMM,0:IJMM,IKMM),TDTDSG(IKMM)
C
      DO 20 K=1,IKMM
         TDTDSG(K) = DT / DELTAS(K)
         DO 10 IJ=0,IJMM
         DO 10 IL=1,ILMM
            DNEW(IL,IJ,K) = 0.
 10      CONTINUE
 20   CONTINUE
C
      DO 40 IK=2,IKMMM1
         IKM = IK-1
         IKP = IK+1
         DO 30 IJ=0,IJMM
         DO 30 IL=1,ILMM
            TOP(IL,IJ,IK) = (D(IL,IJ,IK) - D(IL,IJ,IKP)) *
     1                 (D(IL,IJ,IKM) - D(IL,IJ,IK))
            BOT(IL,IJ,IK) = D(IL,IJ,IKM) - D(IL,IJ,IKP)
 30      CONTINUE
 40   CONTINUE
 
      DO 50 IJ=0,IJMM
      DO 50 IL=1,ILMM
         TOP(IL,IJ,1)    = 0.
         TOP(IL,IJ,IKMM) = 0.
         BOT(IL,IJ,1)    = 1.
         BOT(IL,IJ,IKMM) = 1.
 50   CONTINUE
C
      DO 60 K=1,IKMM
         DO 55 IJ=0,IJMM
         DO 55 IL=1,ILMM
            IF ((TOP(IL,IJ,K).GT.0.0).AND.(BOT(IL,IJ,K).NE.0.0)) THEN
               DELTA(IL,IJ,K)=2.*TOP(IL,IJ,K)/BOT(IL,IJ,K)
            ELSE
               DELTA(IL,IJ,K) = 0.0
            END IF
 55      CONTINUE
 60   CONTINUE
 
C     VERTICAL FLUX AT INTERIOR LEVELS OF TROPOSPHERE
      DO 80 IK=1,KTROP
         IM = IK - 1
         IP = IK + 1
         DO 75 IJ=0,IJMM
         DO 75 IL=1,ILMM
            IF (SGTROP(IL,IJ,IK) .GE. 0.0)  THEN
               FLUX1(IL,IJ,IK) = D(IL,IJ,IP) +
     1          (.5* (1.+SGTROP(IL,IJ,IK) * TDTDSG(IP)))
     1                * DELTA(IL,IJ,IP)
            ELSE
               FLUX1(IL,IJ,IK)=D(IL,IJ,IK)-(.5*(1.-SGTROP(IL,IJ,IK)
     1          *TDTDSG(IK))) * DELTA(IL,IJ,IK)
            END IF
 75      CONTINUE
 80   CONTINUE
C
      DO 85 IK=2,KTROP
         IM = IK - 1
         IP = IK + 1
         DO 82 IJ=0,IJMM
         DO 82 IL=1,ILMM
            DNEW(IL,IJ,IK) =
     1         TDTDSG(IK) * (SGTROP(IL,IJ,IM)*FLUX1(IL,IJ,IM) -
     1         SGTROP(IL,IJ,IK)*FLUX1(IL,IJ,IK))
 82      CONTINUE
  85  CONTINUE
C
C     VERTICAL FLUX FOR LEVEL 1 IS DETERMINED.
      IK = 1
      IM = IK - 1
      IP = IK + 1
      DO 88 IJ=0,IJMM
      DO 88 IL=1,ILMM
         IF (SGBOT(IL,IJ).GE.0.0) THEN
            FLUX1(IL,IJ,IM)=D(IL,IJ,IK)+(.5*(1.+SGBOT(IL,IJ)
     1       * TDTDSG(IK)))* DELTA(IL,IJ,IK)
         ELSE
            FLUX1(IL,IJ,IM)= FD0(IL,IJ)
         END IF
         DNEW(IL,IJ,IK)=TDTDSG(IK)*(SGBOT(IL,IJ)*FLUX1(IL,IJ,IM)
     1    -SGTROP(IL,IJ,IK)*FLUX1(IL,IJ,IK))
 88   CONTINUE
 
C     VERTICAL FLUX IN STRATOSPHERE.
      DO 100 IK=KTROP,IKMMM1
         IM = IK - 1
         IP = IK + 1
         IMM = IK-1-KTROP
         IPP = IK+1-KTROP
         IKK = IK  -KTROP
         DO 90 IJ=0,IJMM
         DO 90 IL=1,ILMM
            IF (SGSTRAT(IL,IJ,IKK) .GE. 0.0) THEN
               FLUX2(IL,IJ,IKK) = D(IL,IJ,IP) +
     1         (.5* (1.+SGSTRAT(IL,IJ,IKK) * TDTDSG(IP)))
     1                * DELTA(IL,IJ,IP)
            ELSE
               FLUX2(IL,IJ,IKK)=D(IL,IJ,IK)-
     1         (.5*(1.-SGSTRAT(IL,IJ,IKK)*TDTDSG(IK)))*DELTA(IL,IJ,IK)
            END IF
 90       CONTINUE
 100   CONTINUE
C
      DO 105 IK=KTROP+1,IKMMM1
         IM = IK - 1
         IP = IK + 1
         IMM = IK-1-KTROP
         IPP = IK+1-KTROP
         IKK = IK  -KTROP
         DO 102 IJ=0,IJMM
         DO 102 IL=1,ILMM
            DNEW(IL,IJ,IK) =
     1       TDTDSG(IK) * (SGSTRAT(IL,IJ,IMM)*FLUX2(IL,IJ,IMM) -
     1       SGSTRAT(IL,IJ,IKK)*FLUX2(IL,IJ,IKK))
 102     CONTINUE
 105  CONTINUE
 
C     VERTICAL VELOCITY AT TOP OF MODEL IS ZERO SO
      IK = IKMM
      IMM = IK-1-KTROP
      DO 120 IJ=0,IJMM
      DO 120 IL=1,ILMM
         DNEW(IL,IJ,IK) =
     1    TDTDSG(IK)*SGSTRAT(IL,IJ,IMM)*FLUX2(IL,IJ,IMM)
 120  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE SGDOT3(SGTROP,SGSTRAT,SGBOT,DNEW,PS1,PS1INV,COSLAT,TDT)
      IMPLICIT REAL(A-H,O-Z)
 
C     19 DECEMBER 1989
C     VERTICAL VELOCITY FIELD IS DETERMINED AS A RESIDUAL BY FORCING
C     A CONSTANT FIELD TO REMAIN CONSTANT AFTER TRANSPORT.
 
C     SGTROP       VERTICAL VELOCITY IN TROPOSPHERE.
C     SGSTRAT      VERTICAL VELOCITY IN STRATOSPHERE.
C     SGBOT        VERTICAL VELOCITY AT SURFACE
C     DNEW         -(RESIDUAL AFT HORIZONTAL FLUXES CALCULATED.)
C
      PARAMETER (ILMM=144,IJMM=88,KSTRAT=11,IJMMP1=IJMM+1,KTROP=8)
      PARAMETER (IKMM=KTROP+KSTRAT)
      PARAMETER (IKMMM1=IKMM-1,ILIJ=ILMM*IJMMP1,KST1=KSTRAT-1)
      PARAMETER (IJMMM1=IJMM-1)
 
      REAL SGTROP(ILMM,0:IJMM,KTROP),PS1INV(ILMM,0:IJMM)
      REAL DNEW(ILMM,0:IJMM,IKMM),SGSTRAT(ILMM,0:IJMM,0:KST1)
      REAL SGBOT(ILMM,0:IJMM),COSLAT(0:IJMM),PS1(ILMM,0:IJMM)
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     2 SIGMA(IKMM), PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),
     3 SGINT(0:IKMM)
 
      COMMON // TDTDSG(IKMM),SGAVG2(ILMM,0:IJMM)
C
      PI = 4.0 * ATAN (1.)
      DILMM = 1. / ILMM
 
      DO 20 K=1,IKMM
         TDTDSG(K) = DELTAS(K) / TDT
 20   CONTINUE
C
      IK = IKMM
      IMM = IKMM-KTROP-1
      DO 30 IJ=0,IJMM
      DO 30 IL=1,ILMM
         SGSTRAT(IL,IJ,IMM) = DNEW(IL,IJ,IK) * TDTDSG(IK)
 30   CONTINUE
*******  BEGIN W CONSTRAINT SECTION 1
C      DCON = (PI/IJMM) * 0.5
C      SGAVG= 0.
C      DO 40 IJ=1,IJMMM1
C         SUM = 0.
C         DO 35 IL=1,ILMM
C            SUM = SUM + SGSTRAT(IL,IJ,IMM)
C35       CONTINUE
C         SUM = SUM * DILMM
C         SGAVG=SGAVG+ SUM * COSLAT(IJ)
C 40   CONTINUE
C      DO 45 IJ=0,IJMM
C      DO 45 IL=1,ILMM
C         SGSTRAT(IL,IJ,IMM) = SGSTRAT(IL,IJ,IMM)-DCON*SGAVG
C45    CONTINUE
******     END W CONSTRAINT SECTION 1
      DO 80 IK=IKMMM1,KTROP+1,-1
         IKK = IK - KTROP
         IMM = IK - KTROP-1
         DO 50 IJ=0,IJMM
         DO 50 IL=1,ILMM
         SGSTRAT(IL,IJ,IMM)=SGSTRAT(IL,IJ,IKK)+DNEW(IL,IJ,IK)*
     1       TDTDSG(IK)
 50      CONTINUE
*******  BEGIN  W CONSTRAINT SECTION 2
C         SGAVG= 0.
C         DO 60 IJ=1,IJMMM1
C            SUM = 0.
C            DO 55 IL=1,ILMM
C              SUM = SUM + SGSTRAT(IL,IJ,IMM)
C55          CONTINUE
C            SGAVG=SGAVG+SUM*DILMM*COSLAT(IJ)
C  60     CONTINUE
C         DO 65 IJ=0,IJMM
C         DO 65 IL=1,ILMM
C            SGSTRAT(IL,IJ,IMM)=SGSTRAT(IL,IJ,IMM)-DCON*SGAVG
C65       CONTINUE
*****    END  W CONSTRAINT SECTION 2
  80  CONTINUE
C
C     DYNAMIC INTERFACE CONDITION
      DO 85 IJ=0,IJMM
      DO 85 IL=1,ILMM
         SGTROP(IL,IJ,KTROP)=SGSTRAT(IL,IJ,0) * (PINT-PTOP)
     1                       *PS1INV(IL,IJ)
 85   CONTINUE
      DO 100 IK=KTROP,2,-1
         IM = IK - 1
         DO 88 IJ=0,IJMM
         DO 88 IL=1,ILMM
            SGTROP(IL,IJ,IM)=SGTROP(IL,IJ,IK)+DNEW(IL,IJ,IK)*TDTDSG(IK)
 88      CONTINUE
*******  BEGIN W CONSTRAINT SECTION 3
C         SGAVG= 0.
C         DO 90 IJ=1,IJMMM1
C            SUM = 0.
C            DO 89 IL=1,ILMM
C               SUM=SUM+SGTROP(IL,IJ,IM)*DILMM*(PS1(IL,IJ)-PINT)
C89          CONTINUE
C            SGAVG=SGAVG+SUM*COSLAT(IJ)
C  90     CONTINUE
C         DO 95 IJ=0,IJMM
C         DO 95 IL=1,ILMM
C            SGAVG2(IL,IJ) = DCON*SGAVG*PS1INV(IL,IJ)
C            SGTROP(IL,IJ,IM)=SGTROP(IL,IJ,IM) - SGAVG2(IL,IJ)
C95       CONTINUE
*****    END  W CONSTRAINT SECTION 3
 100  CONTINUE
 
      IK=1
      DO 105 IJ=0,IJMM
      DO 105 IL=1,ILMM
         SGBOT(IL,IJ)=SGTROP(IL,IJ,IK)+DNEW(IL,IJ,IK)*TDTDSG(IK)
 105  CONTINUE
*******  BEGIN W CONSTRAINT SECTION 4
C      SGAVG= 0.
C      DO 110 IJ=1,IJMMM1
C         SUM = 0.
C         DO 108 IL=1,ILMM
C            SUM = SUM + SGBOT(IL,IJ)*DILMM*(PS1(IL,IJ)-PINT)
C108      CONTINUE
C         SGAVG=SGAVG+SUM*COSLAT(IJ)
C 110  CONTINUE
C      DO 120 IJ=0,IJMM
C      DO 120 IL=1,ILMM
C         SGAVG2(IL,IJ) = DCON*SGAVG*PS1INV(IL,IJ)
C         SGBOT(IL,IJ) = SGBOT(IL,IJ) - SGAVG2(IL,IJ)
C120   CONTINUE
****  END  W CONSTRAINT SECTION 4
 
      RETURN
      END
C
C
C
      SUBROUTINE DIVERT(U1,V1,COSLAT,COSLAT2,PS1,PSU,PSV,DIVINT,DPHI)
      IMPLICIT REAL(A-H,O-Z)
C
C     THE DIVERGENCE IS DETERMINED AT THE TROPOSPHERE-STRATOSPHERE
C     INTERFACE.  DIVINT(I,J,1) IS VALUE WHEN APPROACHING THE INTERFACE
C     FROM BELOW. DIVINT(I,J,2) IS VALUE WHEN APPROACHING THE INTERFACE
C     FROM ABOVE.
 
      PARAMETER(ILMM=144,IL2=ILMM/2,IL2P=IL2+1,ILMMM1=ILMM-1)
      PARAMETER(KTROP=8,KSTRAT=11,IKMM=KTROP+KSTRAT)
      PARAMETER(IJMM=88,IJMMP1=IJMM+1,IJMMM1=IJMM-1)
      PARAMETER(ILIJ=ILMM*IJMMP1,ILIJM=ILMM*IJMM)
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      REAL RDNEW(ILMM)
      REAL U1(ILMM,0:IJMM,IKMM),V1(ILMM,IJMM,IKMM),
     1 PSU(ILMM,0:IJMM), PSV(ILMM,IJMM),DIVINT(ILMM,0:IJMM,2)
      REAL COSLAT(0:IJMM),COSLAT2(IJMM),PS1(ILMM,0:IJMM)
 
      COMMON // UINT(ILMM,0:IJMM),VINT(ILMM,0:IJMM),WK1(ILMM,0:IJMM),
     1 WK2(ILMM,0:IJMM),DTOVDXL(0:IJMM),
     ! DNEW(ILMM), DNEW2(ILMM,0:IJMM)
C
      A = 6.37E6
      PI=4.0*ATAN(1.)
C     SPECIFY RADIUS OF POLAR BOX
      DPHI2 = PI  / 60.
      POLCON = 2./(A*ILMM*SIN(DPHI2))
 
C     STEP1: U AND V ARE DETERMINED AT THE INTERFACE.
      DO 20 IJ=0,IJMM
      DO 20 IL=1,ILMM
         WK1(IL,IJ) = LOG(SIG(KTROP)*PSU(IL,IJ)+PINT)
         WK2(IL,IJ) = LOG(PRESS(1)) - WK1(IL,IJ)
 20   CONTINUE
      AA = LOG(PRESS(1)) - LOG(PINT)
 
      DO 40 IJ=0,IJMM
      DO 40 IL=1,ILMM
         UINT(IL,IJ) = LOG(PINT) - WK1(IL,IJ)
         UINT(IL,IJ) = 1. / (WK2(IL,IJ)) *(U1(IL,IJ,KTROP)* AA
     1               + U1(IL,IJ,KTROP+1)*UINT(IL,IJ))
 40   CONTINUE
 
      DO 60 IJ=1,IJMM
      DO 60 IL=1,ILMM
         WK1(IL,IJ) = LOG(SIG(KTROP)*PSV(IL,IJ)+PINT)
         WK2(IL,IJ) = LOG(PRESS(1)) - WK1(IL,IJ)
         VINT(IL,IJ) = LOG(PINT) - WK1(IL,IJ)
          VINT(IL,IJ) = 1./(WK2(IL,IJ))*(V1(IL,IJ,KTROP)*AA
     1               + V1(IL,IJ,KTROP+1)*VINT(IL,IJ))
 60   CONTINUE
 
C     STEP2: EAST-WEST PART OF DIVERGENCE IS DET FOR TROP SIDE OF INT
      PI=4.0 * ATAN(1.)
      A = 6.37E6
      DIST = 2. * PI * A / ILMM
      DO 80 IJ=0,IJMM
         DTOVDXL(IJ) = 1.0 / (DIST* COSLAT(IJ))
 80   CONTINUE
      IK = 1
      DO 100 IL=1,ILMM
         DIVINT(IL,0,IK) = 0.
 100  CONTINUE
      DO 140 IJ=1,IJMMM1
         DO 120 ILL = 1,ILMM
            IM = ILL - 1
            IP = ILL+ 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
            TERM3 = 1.0
            TERM4 = 1.0
            DNEW(ILL) = DTOVDXL(IJ) *
     1                       (UINT(IP,IJ)*TERM3*PSU(IP,IJ)-
     1                        UINT(ILL,IJ)*TERM4*PSU(ILL,IJ))
         DIVINT(ILL,IJ,IK) = DNEW(ILL)
 120     CONTINUE
 140  CONTINUE
      DO 160 IL=1,ILMM
         DIVINT(IL,IJMM,IK) = 0.
 160  CONTINUE
 
C     STEP3: EAST-WEST PART OF DIVERGENCE IS DET FOR STRAT SIDE OF INT
      IK = 2
      DO 180 IL=1,ILMM
         DIVINT(IL,0,IK) = 0.
 180  CONTINUE
      DO  220 IJ=1,IJMMM1
         DO  200 IIIL = 1,ILMM
            IM = IIIL - 1
            IP = IIIL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
            TERM3 = 1.0
            TERM4 = 1.0
            DNEW(IIIL) =  DTOVDXL(IJ) * (PINT-PTOP) *
     1       (UINT(IP,IJ)*TERM3 - UINT(IIIL,IJ)*TERM4)
         DIVINT(IIIL,IJ,IK) = DNEW(IIIL)
 200     CONTINUE
 220  CONTINUE
      DO 240 IL=1,ILMM
         DIVINT(IL,IJMM,IK) = 0.
 240  CONTINUE
 
C     STEP4: NORTH-SOUTH PT OF DIVERGENCE IS DET FOR TROP SIDE OF INT
C     INITIALIZE CONSTANTS
      DTOVDXJ = 1. / DPHI
      DD = 1. / ILMM
      IK = 1
C     FLUX CALCULATIONS.
      DO 280 IJ = 1,IJMMM1
         IM = IJ - 1
         IP = IJ + 1
         TERM1 = 1.0
         TERM2 = 1.0
         DO 260 IL=1,ILMM
            DNEW2(IL,IJ) = (DTOVDXJ/COSLAT(IJ))*
     1      (VINT(IL,IP)*TERM1*PSV(IL,IP)*COSLAT2(IP)-
     1      VINT(IL,IJ)*TERM2*PSV(IL,IJ)*COSLAT2(IJ))
 260     CONTINUE
 280  CONTINUE
 
C     FLUX INTO POLAR BOXES IS (FLUX OUT OF NEIGHBORING BOX INTO
C     POLAR BOX) - (FLUX INTO NEIGHBORING BOX FROM POLAR BOX)
C     WELL MIXED POLAR BOX APPROXIMATION IS USED.
      IJ = 1
      SUMSP = 0.
      SUMNP = 0.
      DO 300 IL=1,ILMM
         SUMSP = SUMSP + VINT(IL,1)
         SUMNP = SUMNP + VINT(IL,IJMM)
 300  CONTINUE
      DO 320 IL=1,ILMM
         DNEW2(IL,0) = -POLCON*SUMSP
         DNEW2(IL,IJMM) = POLCON*SUMNP
 320  CONTINUE
 
      DO 340 IJ=0,IJMM
      DO 340 IL=1,ILMM
         DIVINT(IL,IJ,IK) = DIVINT(IL,IJ,IK) + DNEW2(IL,IJ)
 340  CONTINUE
 
C     STEP5: NORTH-SOUTH PT OF DIVERGENCE IS DET FOR STRAT SIDE OF INT
      IK = 2
      DO 380 IJ=1,IJMMM1
         IM = IJ - 1
         IP = IJ + 1
         TERM1 = 1.0
         TERM2 = 1.0
         DO 360 IL=1,ILMM
            DNEW2(IL,IJ) =
     1      (DTOVDXJ/COSLAT(IJ))*(PINT-PTOP)*(VINT(IL,IP)*TERM1*
     1        COSLAT2(IP)-VINT(IL,IJ)*TERM2*COSLAT2(IJ))
 360     CONTINUE
 380  CONTINUE
 
C     FLUX INTO POLAR BOXES IS (FLUX OUT OF NEIGHBORING BOX INTO
C     POLAR BOX) - (FLUX INTO NEIGHBORING BOX FROM POLAR BOX)
C     WELL MIXED POLAR BOX APPROXIMATION IS USED.
      IJ = 1
      SUMSP = 0.
      SUMNP = 0.
      DO 400 IL=1,ILMM
         SUMSP = SUMSP + VINT(IL,1)
         SUMNP = SUMNP + VINT(IL,IJMM)
 400  CONTINUE
      DO 420 IL=1,ILMM
         DNEW2(IL,0) =    -POLCON*SUMSP
         DNEW2(IL,IJMM) =  POLCON*SUMNP
 420  CONTINUE
 
C     RELOAD ARRAY
      DO 440 IJ=0,IJMM
      DO 440 IL=1,ILMM
         DIVINT(IL,IJ,IK) = DIVINT(IL,IJ,IK) + DNEW2(IL,IJ)
 440  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE FINT(NPTS,DZ,F,FINT3)
      IMPLICIT REAL(A-H,O-Z)
 
C     NPTS  NUMBER OF PTS TO INTEGRATE OVER.
C     DZ    RANGE OF INTEGRATION LIMITS
C     F     VALUES OF INTEGRATED FUNCTION
C     FINT3 VALUE OF INTEGRAL.
 
C     ORIGINAL STOLEN FROM ERIC NASH
 
      PARAMETER (N=144,JGMAX=88,NJG=N*(JGMAX+1))
      PARAMETER (KTROP=8,KSTRAT=11,KTOP=KTROP+KSTRAT)
 
      REAL F(N,0:JGMAX,NPTS), FINT3(N,0:JGMAX)
 
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),
     2 SIG(KTOP),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)
 
      REAL SUM(N,0:JGMAX), SUM1(N,0:JGMAX),
     ! SUM2(N,0:JGMAX), SUM4(N,0:JGMAX), RESULT(N,0:JGMAX)
 
      DO 20 IJ=0,JGMAX
      DO 20 IIIL=1,N
         FINT3(IIIL,IJ) =0.
 20   CONTINUE
      ODNM1=DZ/(NPTS-1.)
      OD3=1./3.
      IF  (MOD(NPTS,2).EQ.0)  THEN
          IF  (NPTS.EQ.2)  THEN
              DO 40 IJ=0,JGMAX
              DO 40 IIIL=1,N
                 FINT3(IIIL,IJ) =.5*DZ*(F(IIIL,IJ,1)+F(IIIL,IJ,2))
 40           CONTINUE
              RETURN
          ENDIF
          DO 60 IJ=0,JGMAX
          DO 60 IL=1,N
             RESULT(IL,IJ) =.375*ODNM1*
     1       (F(IL,IJ,1)+F(IL,IJ,4)+3.*(F(IL,IJ,2)+F(IL,IJ,3)))
             FINT3(IL,IJ)=RESULT(IL,IJ)
 60       CONTINUE
          IF  (NPTS.EQ.4)  RETURN
          DO 80 IJ=0,JGMAX
          DO 80 IL=1,N
             SUM4(IL,IJ) =0.
 80       CONTINUE
          DO  120 I=5,NPTS-1,2
              DO 100 IJ=0,JGMAX
              DO 100 IL=1,N
                 SUM4(IL,IJ) = SUM4(IL,IJ) +F(IL,IJ,I)
 100          CONTINUE
120       CONTINUE
          DO 140 IJ=0,JGMAX
          DO 140 IL=1,N
             SUM2(IL,IJ) =0.
 140      CONTINUE
          DO  160 I=6,NPTS-2,2
          DO 160 IJ=0,JGMAX
          DO 160 IL=1,N
              SUM2(IL,IJ) =SUM2(IL,IJ) +F(IL,IJ,I)
160       CONTINUE
          DO 180 IJ=0,JGMAX
          DO 180 IL=1,N
             FINT3(IL,IJ) =OD3*ODNM1*
     !                (F(IL,IJ,4)+F(IL,IJ,NPTS)+2.*SUM2(IL,IJ)+
     !                 4.*SUM4(IL,IJ)) + RESULT(IL,IJ)
 180      CONTINUE
          RETURN
          ENDIF
      DO 200 IJ=0,JGMAX
      DO 200 IL=1,N
         SUM(IL,IJ)  =0.
         SUM1(IL,IJ) =0.
 200  CONTINUE
         DO  220  I=2,NPTS-1
         DO 220 IJ=0,JGMAX
         DO 220 IL=1,N
            SUM(IL,IJ) =SUM(IL,IJ) +F(IL,IJ,I)
220      CONTINUE
      DO   225  I=3,NPTS-2,2
      DO 225 IJ=0,JGMAX
      DO 225 IL=1,N
          SUM1(IL,IJ) =SUM1(IL,IJ) +F(IL,IJ,I)
225   CONTINUE
 
      DO 235 IJ=0,JGMAX
      DO 235 IL=1,N
         FINT3(IL,IJ) =OD3*ODNM1*
     !    (4.*SUM(IL,IJ)+F(IL,IJ,1)+F(IL,IJ,NPTS)-2.*SUM1(IL,IJ))
 235  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE DIFFY(DENJK,COSLAT,COSLAT2,DT,ADPHI,F1,PS1INV,DFCOEF)
 
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (ILMM=144,IJMM=88,IJMMM1=IJMM-1)
      PARAMETER (IJMMP1=IJMM+1,ILIJ=ILMM*IJMMP1)
      PARAMETER (KTROP=8,KSTRAT=11)
      PARAMETER (IKMM=KTROP+KSTRAT,ILIJM=ILMM*IJMMM1)
 
      LOGICAL DIFSTRAT
      REAL RDNEW(ILMM)
      REAL DENJK(ILMM,0:IJMM,IKMM), COSLAT(0:IJMM)
      REAL PS1INV(ILMM,0:IJMM)
      REAL F1(ILMM,0:IJMM,IKMM),COSLAT2(IJMM)
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      COMMON //  FLUX(ILMM,IJMM), PS11(ILMM,0:IJMM)
 
C     INITIALIZE CONSTANTS
      DIFSTRAT = .FALSE.
      DTOVDXJ  = DT / ADPHI
      A= 6.37E6
      PI = 4.0 * ATAN(1.)
C     RADIUS OF POLAR BOX IS SET TO 3 DEGREES.
      DPHI2 = 3.* PI  / 180.
      ADPHI2 = A * DPHI2
 
C     DETERMINE PSF-PINT BY INVERTING 1. / (PSF-PINT)  (BRILLIANT)
 
      DO 20 IJ=0,IJMM
      DO 20 IL=1,ILMM
         PS11(IL,IJ) = 1. / PS1INV(IL,IJ)
 20   CONTINUE
 
      POLCON = 2.*DT/(A*ILMM*SIN(DPHI2))
      DD = 1. / ILMM
      DO 100 IK=1,KTROP
         DO 30 IL=1,ILMM
            FLUX(IL,1)= -DFCOEF*(DENJK(IL,1,IK)*PS11(IL,1)
     1       -DENJK(IL,0,IK)*PS11(1,0))/ (ADPHI*0.5+ADPHI2)
 30      CONTINUE
         DO 40 IJ=2,IJMMM1
         DO 40 IL=1,ILMM
            IM = IJ - 1
            FLUX(IL,IJ)= -DFCOEF*(DENJK(IL,IJ,IK)*PS11(IL,IJ)
     1       -DENJK(IL,IM,IK)*PS11(IL,IM))/ADPHI
 40      CONTINUE
         DO 50 IL=1,ILMM
            FLUX(IL,IJMM)=-DFCOEF*(DENJK(IL,IJMM,IK)
     1       *PS11(IL,IJMM)-DENJK(IL,IJMMM1,IK)*PS11(IL,IJMMM1))/
     2       (ADPHI*0.5+ADPHI2)
 50      CONTINUE
C
        DO 65 IJ=1,IJMMM1
        DO 65 IL=1,ILMM
           IM = IJ - 1
           IP = IJ + 1
           F1(IL,IJ,IK) =-(DTOVDXJ*PS1INV(IL,IJ)/COSLAT(IJ))*
     1     (FLUX(IL,IP)*COSLAT2(IP)-FLUX(IL,IJ)*COSLAT2(IJ))
 65     CONTINUE
 
C       FLUX INTO POLAR BOXES IS (FLUX OUT OF NEIGHBORING BOX INTO
C       POLAR BOX) - (FLUX INTO NEIGHBORING BOX FROM POLAR BOX)
C       WELL MIXED POLAR BOX APPROXIMATION IS USED.
        DO 70 IL=1,ILMM
           F1(IL,0,IK) = FLUX(IL,1)
           F1(IL,IJMM,IK) = FLUX(IL,IJMM)
 70     CONTINUE
 
        SUMSP = 0.
        SUMNP = 0.
        DO 80 IL=1,ILMM
           SUMSP = SUMSP-POLCON*F1(IL,0,IK)
           SUMNP = SUMNP+POLCON*F1(IL,IJMM,IK)
 80     CONTINUE
 
        DO 90 IL=1,ILMM
           F1(IL,0,IK)    = SUMSP
           F1(IL,IJMM,IK) = SUMNP
 90     CONTINUE
 100  CONTINUE
 
C     BEGIN CALCULATIONS FOR STRATOSPHERE.
      IF (DIFSTRAT) THEN
         DO 200 IK=KTROP+1,IKMM
            DO 110 IL=1,ILMM
               FLUX(IL,1)= -DFCOEF*(DENJK(IL,1,IK)-DENJK(IL,0,IK))/
     1           (ADPHI*0.5+ADPHI2)
 110        CONTINUE
            DO 120 IJ=2,IJMMM1
            DO 120 IL=1,ILMM
               IM = IJ - 1
               FLUX(IL,IJ)=-DFCOEF*(DENJK(IL,IJ,IK)-DENJK(IL,IM,IK))
     1                        /ADPHI
 120        CONTINUE
            DO 130 IL=1,ILMM
             FLUX(IL,IJMM)=-DFCOEF*
     1       (DENJK(IL,IJMM,IK)-DENJK(IL,IJMMM1,IK))/(ADPHI*0.5+ADPHI2)
 130        CONTINUE
 
            DO 165 IJ=1,IJMMM1
            DO 165 IL=1,ILMM
               IM = IJ - 1
               IP = IJ + 1
               F1(IL,IJ,IK) =-(DTOVDXJ/COSLAT(IJ))*
     1         (FLUX(IL,IP)*COSLAT2(IP)-FLUX(IL,IJ)*COSLAT2(IJ))
 165        CONTINUE
 
C           FLUX INTO POLAR BOXES IS (FLUX OUT OF NEIGHBORING BOX INTO
C           POLAR BOX) - (FLUX INTO NEIGHBORING BOX FROM POLAR BOX)
C           WELL MIXED POLAR BOX APPROXIMATION IS USED.
 
            DO 175 IL=1,ILMM
               F1(IL,0,IK) =    FLUX(IL,1)
               F1(IL,IJMM,IK) = FLUX(IL,IJMM)
 175        CONTINUE
 
            SUMSP = 0.
            SUMNP = 0.
            DO 180 IL=1,ILMM
               SUMSP = SUMSP-POLCON*F1(IL,0,IK)
               SUMNP = SUMNP+POLCON*F1(IL,IJMM,IK)
 180        CONTINUE
 
            DO 190 IL=1,ILMM
               F1(IL,0,IK)    = SUMSP
               F1(IL,IJMM,IK) = SUMNP
 190        CONTINUE
 200     CONTINUE
      ELSE
            DO 195 IK=KTROP+1,IKMM
            DO 195 IJ=0,IJMM
            DO 195 IL=1,ILMM
               F1(IL,IJ,IK) = 0.
 195        CONTINUE
      END IF
 
      RETURN
      END
C
C
C
      SUBROUTINE DIFFX(DENJK,COSLAT,DT,F1,PS1INV,DFCOEF)
 
C     DFCOEF*(2ND DERIVATIVE OF MIXING RATIO WRT X IS CALCULATED)
C
C     DENJK    = MIXING RATIO
C     DFCOEF   = DIFFUSION COEFFICIENT
C     COSLAT=  = COSINE OF LATITUDE IN CENTER OF BOX
C     DT       = TIME STEP
C     F1       = CHANGE DUE TO DIFFUSION
C     PS1INV   = 1./(PSF-PINT)  WHERE PSF = PSF IN CENTER OF BOX
C
C     FIRST DRAFT 5/31/90
C
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (ILMM=144,IJMM=88,ILMMM1=ILMM-1,NLATS2=87)
      PARAMETER (IJMMP1=IJMM+1,IJMMM1=IJMM-1,KTROP=8,KSTRAT=11)
      PARAMETER (IKMM=KTROP+KSTRAT,L1=1,L2=NLATS2)
 
      LOGICAL DIFSTRAT
      REAL DENJK(ILMM,0:IJMM,IKMM), COSLAT(0:IJMM)
      REAL F1(ILMM,0:IJMM,IKMM),PS1INV(ILMM,0:IJMM)
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      COMMON // PSI3(ILMM,0:IJMM),DTOVDXL(0:IJMM)
C
      PI=4.0 * ATAN(1.)
      A = 6.37E6
      DIST = 2. * PI * A / FLOAT(ILMM)
C
      DO 60 IJ=0,IJMM
      DO 60 IL=1,ILMM
         PSI3(IL,IJ)=1./(PS1INV(IL,IJ))
 60   CONTINUE
 
C     NOTE DENOMINATOR IS SQUARED.
      DO 100 IJ=0,IJMM
         DTOVDXL(IJ) =DT/(DIST*COSLAT(IJ))**2
 100  CONTINUE
 
C     BEGIN CALCULATION OF TROPOSPHERIC DIFFUSION.
      DO 110 IK=1,KTROP
      DO 110 IJ=1,IJMMM1
         DO 105 IL=1,ILMM
            IM = IL - 1
            IP = IL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
            F1(IL,IJ,IK)=DFCOEF*DTOVDXL(IJ)* PS1INV(IL,IJ)*
     1      (PSI3(IP,IJ)*DENJK(IP,IJ,IK) -2.*
     1       PSI3(IL,IJ)*DENJK(IL,IJ,IK)    +
     1       PSI3(IM,IJ)*DENJK(IM,IJ,IK))
 105     CONTINUE
 110  CONTINUE
 
C     POSSIBLE CALCULATION OF STRATOSPHERIC DIFFUSION BEGINS HERE.
      DIFSTRAT = .FALSE.
      IF (DIFSTRAT) THEN
         DO 145 IK=KTROP+1,IKMM
         DO 145 IJ=1,IJMMM1
            DO 140 IIIL=1,ILMM
               IM = IIIL - 1
               IP = IIIL + 1
               IF (IM .LT. 1) IM = ILMM
               IF (IP.GT. ILMM) IP= 1
               F1(IL,IJ,IK)=DFCOEF*DTOVDXL(IJ)*
     1         (DENJK(IP,IJ,IK)-2.*DENJK(IL,IJ,IK)+DENJK(IM,IJ,IK))
 140        CONTINUE
 145     CONTINUE
      ELSE
         DO 150 IK=KTROP+1,IKMM
         DO 150 IJ=1,IJMMM1
         DO 150 IL=1,ILMM
            F1(IL,IJ,IK) = 0.
 150     CONTINUE
      END IF
 
      DO 220 IK=1,IKMM
      DO 220 IL=1,ILMM
         F1(IL,0   ,IK) = 0.
         F1(IL,IJMM,IK) = 0.
 220  CONTINUE
 
      RETURN
      END
