      SUBROUTINE XADV(DENJK,COSLAT,DT,U1,F1,PS1INV,PSU,L1,L2,UPSTREAM,
     1 DIAG)
      IMPLICIT REAL(A-H,O-Z)
 
C     FLUX CALCULATIONS VIA VANLEER'S SCHEME. FOR A DISCUSSION OF
C     THIS SCHEME SEE HAWLEY,SMARR, AND WILSON'S JUNE 1984 ASTROPHYSICALL
C     JL. SUPPLEMENT SERIES PAPER.
 
      PARAMETER (ILMM=144,IJMM=88,ILMMM1=ILMM-1)
      PARAMETER (IJMMP1=IJMM+1,IJMMM1=IJMM-1,KTROP=8,KSTRAT=11)
      PARAMETER (IKMM=KTROP+KSTRAT)
      PARAMETER (JDIAG2=3,KDIAG2=3,NVARS=10)
 
C     EAST-WEST ADVECTION IS PERFORMED.
 
      REAL DENJK(ILMM,0:IJMM,IKMM), COSLAT(0:IJMM)
      REAL U1(ILMM,0:IJMM,IKMM),PSU(ILMM,0:IJMM)
      REAL F1(ILMM,0:IJMM,IKMM),PS1INV(ILMM,0:IJMM)
      LOGICAL UPSTREAM,DIAG
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      COMMON / DIAGNOSTIC / DD(ILMM,JDIAG2,KDIAG2,NVARS),NNLAT(JDIAG2),
     1 NNLEV(KDIAG2)
 
      REAL DTOVDXL(0:IJMM), BOT(ILMM),
     ! DNEW(ILMM),D(ILMM), DELTA(ILMM), TOP(ILMM),FLUX(ILMM)
C
C     NOTE: U(1) IS EAST-WEST COMP OF WIND AT LON (-360./(2.*ILMM)
C     DENJK IS THE CHEMICAL NUMBER DENSITY FOR EACH SPECIES
C     JSPD IS NUMBER OF CHEMICAL SPECIES
C
      PI=4.0 * ATAN(1.)
      A = 6.37E6
      DIST = 2. * PI * A / ILMM
C
      DO 20 IJ=0,IJMM
         DTOVDXL(IJ) = DT  / (DIST* COSLAT(IJ))
 20   CONTINUE
 
      DO 200 IK=1,KTROP
         DO 40 IL=1,ILMM
            F1(IL,0   ,IK) = 0.
            F1(IL,IJMM,IK) = 0.
 40      CONTINUE
 
      DO 180 IJ=L1,L2
         DO 60 IL=1,ILMM
            D(IL) = DENJK(IL,IJ,IK)
 60      CONTINUE
         IF (.NOT.UPSTREAM) THEN
            DO 80 IIIL=1,ILMM
               IM = IIIL - 1
               IP = IIIL + 1
               IF (IM .LT. 1) IM = ILMM
               IF (IP.GT. ILMM) IP= 1
               TOP(IIIL) = (D(IIIL)-D(IM)) * (D(IP)- D(IIIL))
               BOT(IIIL) = D(IP) - D(IM)
               IF ((TOP(IIIL).GT.0.).AND.(BOT(IIIL).NE.0.0)) THEN
                  DELTA(IIIL) = 2. * TOP(IIIL) / BOT(IIIL)
               ELSE
                  DELTA(IIIL) = 0.
               END IF
 80         CONTINUE
         ELSE
            DO 100 IIIL=1,ILMM
               TOP(IIIL) = 0.
               BOT(IIIL) = 1.
               DELTA(IIIL) = 0.
 100        CONTINUE
         END IF
         DO 120 IL = 1,ILMM
            IM = IL - 1
            IP = IL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
c           IF (U1(IP,IJ,IK).GE.0.0) THEN
c              FLUX(IP)  =  D(IL) + (.5 *
c    1         (1.- U1(IP,IJ,IK) * DTOVDXL(IJ))) * DELTA(IL)
c           ELSE
c              FLUX(IP)  =  D(IP)
c    1         - (.5 * (1.+ U1(IP,IJ,IK)*DTOVDXL(IJ)))*DELTA(IP)
c           END IF
            flux(ip) = cvmgp(d(il)+(.5*(1.-u1(ip,ij,ik)*dtovdxl(ij)))*
     1       delta(il),d(ip)-(.5*(1.+u1(ip,ij,ik)*dtovdxl(ij)))*
     1       delta(ip),u1(ip,ij,ik))
 120     CONTINUE
 
C        CALCULATE EAST-WEST FLUX.
         IF (DIAG) THEN
            DO 140 KK=1,KDIAG2
            DO 140 JJ=1,JDIAG2
               IF ((IJ.EQ.NNLAT(JJ)).AND.(IK.EQ.NNLEV(KK))) THEN
                  DO 130 IL=1,ILMM
                     DD(IL,JJ,KK,8) = U1(IL,IJ,IK)*FLUX(IL)*PSU(IL,IJ)
 130              CONTINUE
               END IF
 140       CONTINUE
         END IF
 
         DO 160 IL=1,ILMM
            IM = IL - 1
            IP = IL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
            DNEW(IL) = -DTOVDXL(IJ) * PS1INV(IL,IJ)*
     1                       (U1(IP,IJ,IK)*FLUX(IP) *PSU(IP,IJ)-
     1                        U1(IL,IJ,IK)*FLUX(IL) *PSU(IL,IJ))
            F1(IL,IJ,IK) = DNEW(IL)
 160     CONTINUE
 180  CONTINUE
 200  CONTINUE
 
      DO 380 IK=KTROP+1,IKMM
      DO 220 IL=1,ILMM
         F1(IL,0,IK) = 0.
         F1(IL,IJMM,IK) = 0.
 220  CONTINUE
      DO 360 IJ=L1,L2
         DO 240 IL=1,ILMM
            D(IL) = DENJK(IL,IJ,IK)
 240     CONTINUE
         IF (.NOT.UPSTREAM) THEN
            DO 260 IIIL=1,ILMM
               IM = IIIL - 1
               IP = IIIL + 1
               IF (IM .LT. 1) IM = ILMM
               IF (IP.GT. ILMM) IP= 1
               TOP(IIIL) = (D(IIIL)-D(IM)) * (D(IP)- D(IIIL))
               BOT(IIIL) = D(IP) - D(IM)
               IF ((TOP(IIIL).GT.0.).AND.(BOT(IIIL).NE.0.0)) THEN
                  DELTA(IIIL) = 2. * TOP(IIIL) / BOT(IIIL)
               ELSE
                  DELTA(IIIL) = 0.
               END IF
 260        CONTINUE
         ELSE
            DO 280 IIIL=1,ILMM
               TOP(IIIL) = 0.
               BOT(IIIL) = 1.
               DELTA(IIIL) = 0.
 280        CONTINUE
         END IF
 
         DO 300 IIIL = 1,ILMM
            IM = IIIL - 1
            IP = IIIL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
c           IF (U1(IP,IJ,IK).GE.0.0) THEN
c              FLUX(IP)  =  D(IIIL) + (.5 *
c    1         (1.- U1(IP,IJ,IK) * DTOVDXL(IJ))) * DELTA(IIIL)
c           ELSE
c              FLUX(IP)  =  D(IP)
c    1         - (.5 * (1.+ U1(IP,IJ,IK)*DTOVDXL(IJ)))*DELTA(IP)
c           END IF
            flux(ip) = cvmgp(d(iiil)+(.5*(1.-u1(ip,ij,ik)*dtovdxl(ij)))*
     1       delta(iiil),d(ip)-(.5*(1.+u1(ip,ij,ik)*dtovdxl(ij)))*
     2         delta(ip),u1(ip,ij,ik))
 300     CONTINUE
 
C        CALCULATE EAST-WEST FLUX.
         IF (DIAG) THEN
            DO 335 KK=1,KDIAG2
            DO 335 JJ=1,JDIAG2
               IF ((IJ.EQ.NNLAT(JJ)).AND.(IK.EQ.NNLEV(KK))) THEN
                  DO 330 IL=1,ILMM
                     DD(IL,JJ,KK,8) = U1(IL,IJ,IK)*FLUX(IL)
 330              CONTINUE
               END IF
 335       CONTINUE
         END IF
 
         DO 340  IIIL=1,ILMM
            IM = IIIL - 1
            IP = IIIL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
 
            DNEW(IIIL) = -DTOVDXL(IJ)
     1       * (U1(IP,IJ,IK)*FLUX(IP)- U1(IIIL,IJ,IK)*FLUX(IIIL))
            F1(IIIL,IJ,IK) = DNEW(IIIL)
 340     CONTINUE
 360  CONTINUE
 380  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE XADVHLAT(DENJK,COSLAT,DT,U1,F1,PS1INV,PSU,L1,L2,NDT,
     1 UPSTREAM,DIAG)
      IMPLICIT REAL(A-H,O-Z)
 
C     FLUX CALCULATIONS VIA VANLEER'S SCHEME. FOR A DISCUSSION OF
C     THIS SCHEME SEE HAWLEY,SMARR, AND WILSON'S JUNE 1984 ASTROPHYSICAL
C     JL. SUPPLEMENT SERIES PAPER.
 
C     EAST-WEST ADVECTION IS PERFORMED AT HIGH LATITUDES WHERE
C     INSTABILITY MAY BE A PROBLEM. CHANGE TERM (F1) OVER TIME DT IS
C     DETERMINED BY SUBDIVIDING DT INTO NDT INTERVALS
C
C     WARNING!! AS CURRENTLY CODE ONLY PART OF F1 IS DETERMINED DURING
C     SUBROUTINE CALL.
 
      PARAMETER (ILMM=144,IJMM=88,ILMMM1=ILMM-1)
      PARAMETER (IJMMP1=IJMM+1,IJMMM1=IJMM-1,KTROP=8,KSTRAT=11)
      PARAMETER (ILIJ=ILMM*IJMMP1,IKMM=KTROP+KSTRAT)
      PARAMETER (JDIAG2=3,KDIAG2=3,NVARS=10)
 
C     EAST-WEST ADVECTION IS PERFORMED.
 
      REAL  DENJK(ILMM,0:IJMM,IKMM), COSLAT(0:IJMM)
      REAL U1(ILMM,0:IJMM,IKMM),PSU(ILMM,0:IJMM)
      REAL F1(ILMM,0:IJMM,IKMM),PS1INV(ILMM,0:IJMM)
      LOGICAL UPSTREAM,DIAG
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      COMMON / DIAGNOSTIC / DD(ILMM,JDIAG2,KDIAG2,NVARS), NNLAT(JDIAG2),
     1 NNLEV(KDIAG2)
 
      COMMON //  DTOVDXL(0:IJMM), BOT(ILMM),
     ! DNEW(ILMM),D(ILMM), DELTA(ILMM), TOP(ILMM),FLUX(ILMM)
C
C     NOTE: U(1) IS EAST-WEST COMP OF WIND AT LON (-360./(2.*ILMM)
C     DENJK IS THE CHEMICAL NUMBER DENSITY FOR EACH SPECIES
C     JSPD IS NUMBER OF CHEMICAL SPECIES
C
      PI=4.0 * ATAN(1.)
      A = 6.37E6
      DIST = 2. * PI * A / FLOAT(ILMM)
C
      DT2 = DT / NDT
      DO 10 IJ=0,IJMM
         DTOVDXL(IJ) = DT2 / (DIST* COSLAT(IJ))
 10   CONTINUE
 
C     COPY HIGH LATITUDE DENSITY ONLY TO OUTPUT ARRAY.
      ILIJ2 = (L2-L1+1)*ILMM
      DO 20 IK=1,IKMM
      DO 20 IJ=L1,L2
      DO 20 IL=1,ILMM
         F1(IL,IJ,IK) = DENJK(IL,IJ,IK)
 20   CONTINUE
 
      DO 120 ISTEP=1,NDT
      DO 110 IK=1,KTROP
      DO 100 IJ=L1,L2
         DO 82 IL=1,ILMM
            D(IL) = F1(IL,IJ,IK)
 82      CONTINUE
 
         IF (.NOT.UPSTREAM) THEN
            DO 85 IIIL=1,ILMM
               IM = IIIL - 1
               IP = IIIL + 1
               IF (IM .LT. 1) IM = ILMM
               IF (IP.GT. ILMM) IP= 1
               TOP(IIIL) = (D(IIIL)-D(IM)) * (D(IP)- D(IIIL))
               BOT(IIIL) = D(IP) - D(IM)
               IF ((TOP(IIIL).GT.0.).AND.(BOT(IIIL).NE.0.0)) THEN
                  DELTA(IIIL) = 2. * TOP(IIIL) / BOT(IIIL)
               ELSE
                  DELTA(IIIL) = 0.
               END IF
 85         CONTINUE
         ELSE
C           SET TOP AND BOT FOR UPSTREAM DIFFERENCING.
            DO 88 IIIL=1,ILMM
              TOP(IIIL) = 0.
              BOT(IIIL) = 1.
              DELTA(IIIL) = 0.
 88         CONTINUE
         END IF
 
         DO 90 IL = 1,ILMM
            IM = IL - 1
            IP = IL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
c           IF (U1(IP,IJ,IK).GE.0.0) THEN
c              FLUX(IP)  =  D(IL) + (.5 *
c    1         (1.- U1(IP,IJ,IK) * DTOVDXL(IJ))) * DELTA(IL)
c           ELSE
c              FLUX(IP)  =  D(IP)
c    1         - (.5 * (1.+ U1(IP,IJ,IK)*DTOVDXL(IJ)))*DELTA(IP)
c           END IF
            flux(ip) = cvmgp(d(il)+(.5*(1.-u1(ip,ij,ik)*dtovdxl(ij)))*
     1        delta(il),d(ip)-(.5*(1.+u1(ip,ij,ik)*dtovdxl(ij)))
     1        *delta(ip),u1(ip,ij,ik))
  90     CONTINUE
 
C        CALCULATE EAST-WEST FLUX.
         IF ((ISTEP.EQ.NDT).AND.(DIAG)) THEN
            DO 93  KK=1,KDIAG2
            DO 93  JJ=1,JDIAG2
               IF ((IJ.EQ.NNLAT(JJ)).AND.(IK.EQ.NNLEV(KK))) THEN
                  DO 92  IL=1,ILMM
                     DD(IL,JJ,KK,8) = U1(IL,IJ,IK)*FLUX(IL)*PSU(IL,IJ)
 92               CONTINUE
               END IF
 93        CONTINUE
         END IF
 
         DO 95 IL=1,ILMM
            IM = IL - 1
            IP = IL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
            DNEW(IL) = -DTOVDXL(IJ) * PS1INV(IL,IJ)*
     1                       (U1(IP,IJ,IK)*FLUX(IP) *PSU(IP,IJ)-
     1                        U1(IL,IJ,IK)*FLUX(IL) *PSU(IL,IJ))
 95      CONTINUE
 
C        CHANGE TERM IS RETURNED IF ISTEP EQUALS NDT
         IF (ISTEP.NE.NDT) THEN
            DO 99 IL=1,ILMM
               F1(IL,IJ,IK) = DNEW(IL) + F1(IL,IJ,IK)
 99         CONTINUE
         ELSE
            DO 98 IL=1,ILMM
               F1(IL,IJ,IK) = DNEW(IL) +(F1(IL,IJ,IK)-
     1          DENJK(IL,IJ,IK))
 98         CONTINUE
         END IF
 100  CONTINUE
 110  CONTINUE
 120  CONTINUE
 
      DO 240 ISTEP=1,NDT
      DO 210 IK=KTROP+1,IKMM
      DO 200 IJ=L1,L2
         DO 182 IL=1,ILMM
            D(IL) = F1(IL,IJ,IK)
 182     CONTINUE
         IF (.NOT.UPSTREAM) THEN
            DO 185 IIIL=1,ILMM
               IM = IIIL - 1
               IP = IIIL + 1
               IF (IM .LT. 1) IM = ILMM
               IF (IP.GT. ILMM) IP= 1
               TOP(IIIL) = (D(IIIL)-D(IM)) * (D(IP)- D(IIIL))
               BOT(IIIL) = D(IP) - D(IM)
               IF ((TOP(IIIL).GT.0.).AND.(BOT(IIIL).NE.0.0)) THEN
                  DELTA(IIIL) = 2. * TOP(IIIL) / BOT(IIIL)
               ELSE
                  DELTA(IIIL) = 0.
               END IF
 185        CONTINUE
         ELSE
            DO 188 IIIL=1,ILMM
               TOP(IIIL) = 0.
               BOT(IIIL) = 1.
               DELTA(IIIL) = 0.
 188        CONTINUE
         END IF
 
         DO 190 IIIL = 1,ILMM
            IM = IIIL - 1
            IP = IIIL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
c           IF (U1(IP,IJ,IK).GE.0.0) THEN
c              FLUX(IP)  =  D(IIIL) + (.5 *
c    1         (1.- U1(IP,IJ,IK) * DTOVDXL(IJ))) * DELTA(IIIL)
c           ELSE
c              FLUX(IP)  =  D(IP)
c    1         - (.5 * (1.+ U1(IP,IJ,IK)*DTOVDXL(IJ)))*DELTA(IP)
c           END IF
            flux(ip) = cvmgp(d(iiil)+(.5*(1.-u1(ip,ij,ik)*dtovdxl(ij)))
     1            *delta(iiil),d(ip)-(.5*(1.+u1(ip,ij,ik)*dtovdxl(ij)))
     1            *delta(ip),u1(ip,ij,ik))
 190     CONTINUE
 
C        CALCULATE EAST-WEST FLUX.
         IF ((ISTEP.EQ.NDT).AND.(DIAG)) THEN
            DO 193  KK=1,KDIAG2
            DO 193  JJ=1,JDIAG2
               IF ((IJ.EQ.NNLAT(JJ)).AND.(IK.EQ.NNLEV(KK))) THEN
                  DO 192  IL=1,ILMM
                     DD(IL,JJ,KK,8) = U1(IL,IJ,IK)*FLUX(IL)
 192              CONTINUE
               END IF
 193        CONTINUE
         END IF
 
         DO 195  IIIL=1,ILMM
            IM = IIIL - 1
            IP = IIIL + 1
            IF (IM .LT. 1) IM = ILMM
            IF (IP.GT. ILMM) IP= 1
            DNEW(IIIL) = -DTOVDXL(IJ)
     1       * (U1(IP,IJ,IK)*FLUX(IP)- U1(IIIL,IJ,IK)*FLUX(IIIL))
 195     CONTINUE
 
 
C        CHANGE TERM IS RETURNED IF ISTEP EQUALS NDT
         IF (ISTEP.NE.NDT) THEN
            DO 199 IL=1,ILMM
               F1(IL,IJ,IK) = DNEW(IL) + F1(IL,IJ,IK)
 199        CONTINUE
         ELSE
            DO 198 IL=1,ILMM
               F1(IL,IJ,IK) = DNEW(IL) +(F1(IL,IJ,IK)-
     1         DENJK(IL,IJ,IK))
 198        CONTINUE
         END IF
 200  CONTINUE
 210  CONTINUE
 240  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE YADV(DENJK,COSLAT,COSLAT2,DT,VFLD,DPHI,F1,PS1INV,PSV,
     1 UPSTREAM,DIAG,POLEDGE)
      IMPLICIT REAL(A-H,O-Z)
 
      PARAMETER (ILMM=144,IJMM=88,IJMMM1=IJMM-1)
      PARAMETER (IL=ILMM,IJMMP1=IJMM+1,ILIJ=ILMM*IJMMP1)
      PARAMETER (IL2=IL/2,IL2P=IL2+1,KTROP=8,KSTRAT=11)
      PARAMETER (IKMM=KTROP+KSTRAT)
      PARAMETER (JDIAG2=3,KDIAG2=3,NVARS=10)
 
C     NORTH-SOUTH ADVECTION
      REAL RDNEW(ILMM)
      REAL DENJK(ILMM,0:IJMM,IKMM), COSLAT(0:IJMM)
      REAL VFLD(ILMM,IJMM,IKMM),PS1INV(ILMM,0:IJMM)
      REAL F1(ILMM,0:IJMM,IKMM),COSLAT2(IJMM)
      REAL PSV(ILMM,IJMM)
      LOGICAL UPSTREAM,DIAG
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
      COMMON / DIAGNOSTIC / DD(ILMM,JDIAG2,KDIAG2,NVARS), NNLAT(JDIAG2),
     1 NNLEV(KDIAG2)
 
      COMMON //  DELTA(ILMM,0:IJMM),FLUX(ILMM,IJMM),
     1    TOP(ILMM,0:IJMM),BOT(ILMM,0:IJMM),
     2    DNEW(ILMM,0:IJMM),D(ILMM,0:IJMM)
 
***** FORMULATION ASSUMES  DENJK(I,0,IK,JS) + DENJK(I,IJMM,IK,JS) ARE
***   GRID BOXES CENTERED AT THE SOUTH AND NORTH POLES. NORTH-SOUTH
***   VELOCITIES DEFINED AT BOX EDGES.
 
C     INITIALIZE CONSTANTS
      DTOVDXJ  = DT / DPHI
      A= 6.37E6
      PI = 4.0 * ATAN(1.)
C     RADIUS OF POLAR BOX IS SET.
      DPHI2 = (90. + POLEDGE) * PI  / 180.
 
      POLCON = 2.*DT/(A*ILMM*SIN(DPHI2))
      DDD = 1. / ILMM
      DO 300 IK=1,KTROP
         DO 20 IJ=0,IJMM
         DO 20 IIIL=1,ILMM
            D(IIIL,IJ) = DENJK(IIIL,IJ,IK)
 20      CONTINUE
 
C        SLOPE OF DENSITY FIELD IS DEFINED AT INTERIOR POINTS.
         IF (.NOT.UPSTREAM) THEN
            DO 60 IJ= 1,IJMMM1
               IM = IJ-1
               IP = IJ+1
               DO 40 IIIL=1,ILMM
                  TOP(IIIL,IJ) = (D(IIIL,IJ) - D(IIIL,IM)) *
     1                           (D(IIIL,IP) - D(IIIL,IJ))
                  BOT(IIIL,IJ) = D(IIIL,IP) - D(IIIL,IM)
 40           CONTINUE
 60         CONTINUE
 
            DO 80 IIIL=1,ILMM
               TOP(IIIL,0)    = 0.
               BOT(IIIL,0)    = 1.
               TOP(IIIL,IJMM) = 0.
               BOT(IIIL,IJMM) = 1.
 80         CONTINUE
C
            DO 100 IJ=0,IJMM
            DO 100 IIIL=1,ILMM
               IF((TOP(IIIL,IJ).GT.0.0).AND.(BOT(IIIL,IJ).NE.0.0)) THEN
                  DELTA(IIIL,IJ) = 2. * TOP(IIIL,IJ) / BOT(IIIL,IJ)
               ELSE
                  DELTA(IIIL,IJ) = 0.0
               END IF
 100        CONTINUE
         ELSE
C           SET TOP AND BOT FOR UPSTREAM DIFFERENCING.
            DO 120 IJ=0,IJMM
            DO 120 IIIL=1,ILMM
               TOP(IIIL,IJ) = 0.
               BOT(IIIL,IJ) = 1.
               DELTA(IIIL,IJ) = 0.
 120        CONTINUE
         END IF
 
C        FLUX CALCULATIONS.
         DO 160 IJ = 0,IJMMM1
            IM = IJ - 1
            IP = IJ + 1
            DO 140 IIIL=1,ILMM
               IF (VFLD(IIIL,IP,IK).GE.0.0) THEN
                  FLUX(IIIL,IP) = D(IIIL,IJ)+(0.5*(1.-VFLD(IIIL,IP,IK)*
     1             DTOVDXJ)) * DELTA(IIIL,IJ)
               ELSE
                  FLUX(IIIL,IP) = D(IIIL,IP)
     1             -(.5*(1.+VFLD(IIIL,IP,IK)*DTOVDXJ))*DELTA(IIIL,IP)
               END IF
 140        CONTINUE
 160     CONTINUE
C
C        NORTH-SOUTH FLUX IS CALCULATED.
         IF (DIAG) THEN
            DO 180 KK=1,KDIAG2
               IF (IK.EQ.NNLEV(KK)) THEN
                  DO 170 JJ=1,JDIAG2
                  JN = NNLAT(JJ)
                  DO 170 IIIL=1,ILMM
                     DD(IIIL,JJ,KK,9)=VFLD(IIIL,JN,IK) *
     1                FLUX(IIIL,JN) * PSV(IIIL,JN)
 170              CONTINUE
               END IF
 180        CONTINUE
         END IF
 
         DO 220 IJ=1,IJMMM1
            IM = IJ - 1
            IP = IJ + 1
            DO 200 IIIL=1,ILMM
               DNEW(IIIL,IJ) =-(DTOVDXJ/COSLAT(IJ))*PS1INV(IIIL,IJ)*
     1        (VFLD(IIIL,IP,IK)*FLUX(IIIL,IP)*PSV(IIIL,IP)*COSLAT2(IP)-
     1         VFLD(IIIL,IJ,IK)*FLUX(IIIL,IJ)*PSV(IIIL,IJ)*COSLAT2(IJ))
 200        CONTINUE
 220      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.
         SUMSP = 0.
         SUMNP = 0.
         DO 240 IIIL=1,ILMM
            SUMSP = SUMSP +  VFLD(IIIL,1,IK)   *FLUX(IIIL,1)
            SUMNP = SUMNP +  VFLD(IIIL,IJMM,IK)*FLUX(IIIL,IJMM)
 240     CONTINUE
         DO 260 IIIL=1,ILMM
            DNEW(IIIL,0) =   -POLCON*SUMSP
            DNEW(IIIL,IJMM) = POLCON*SUMNP
 260     CONTINUE
 
C        RELOAD ARRAY
         DO 280 IJ=0,IJMM
         DO 280 IIIL=1,ILMM
            F1(IIIL,IJ,IK) = DNEW(IIIL,IJ)
 280     CONTINUE
 300  CONTINUE
 
C     BEGIN CALCULATIONS FOR STRATOSPHERE.
      DO 600 IK=KTROP+1,IKMM
         DO 320 IJ = 0,IJMM
         DO 320 IIIL=1,ILMM
            D(IIIL,IJ) = DENJK(IIIL,IJ,IK)
 320     CONTINUE
C        SLOPE OF DENSITY FIELD IS DEFINED AT INTERIOR POINTS.
         IF (.NOT.UPSTREAM) THEN
            DO 360 IJ= 1,IJMMM1
               IM = IJ-1
               IP = IJ+1
               DO 340 IIIL=1,ILMM
                  TOP(IIIL,IJ) = (D(IIIL,IJ) - D(IIIL,IM)) *
     1             (D(IIIL,IP) - D(IIIL,IJ))
                  BOT(IIIL,IJ) = D(IIIL,IP) - D(IIIL,IM)
 340           CONTINUE
 360        CONTINUE
 
            DO 380 IIIL=1,ILMM
               TOP(IIIL,0)  = 0.
               BOT(IIIL,0)  = 1.
               TOP(IIIL,IJMM) = 0.
               BOT(IIIL,IJMM) = 1.
 380        CONTINUE
 
            DO 400 IJ=0,IJMM
            DO 400 IIIL=1,ILMM
              IF ((TOP(IIIL,IJ).GT.0.0).AND.(BOT(IIIL,IJ).NE.0.0))THEN
                 DELTA(IIIL,IJ) = 2. * TOP(IIIL,IJ) / BOT(IIIL,IJ)
              ELSE
                 DELTA(IIIL,IJ) = 0.0
              END IF
 400        CONTINUE
         ELSE
C           SET TOP AND BOT FOR UPSTREAM DIFFERENCING
            DO 420 IJ=0,IJMM
            DO 420 IIIL=1,ILMM
               TOP(IIIL,IJ) = 0.
               BOT(IIIL,IJ) = 1.
               DELTA(IIIL,IJ) = 0.
 420        CONTINUE
         END IF
C       FLUX CALCULATIONS.
        DO 460 IJ = 0,IJMMM1
           IM = IJ - 1
           IP = IJ + 1
           DO 440 IIIL=1,ILMM
              IF (VFLD(IIIL,IP,IK).GE.0.0)  THEN
                 FLUX(IIIL,IP) = D(IIIL,IJ)+ (0.5*(1.-VFLD(IIIL,IP,IK)*
     1            DTOVDXJ)) * DELTA(IIIL,IJ)
              ELSE
                 FLUX(IIIL,IP) = D(IIIL,IP)
     1          -(.5*(1.+VFLD(IIIL,IP,IK)*DTOVDXJ)) * DELTA(IIIL,IP)
              END IF
 440       CONTINUE
 460    CONTINUE
C
C        NORTH-SOUTH FLUX IS CALCULATED.
         IF (DIAG) THEN
            DO 490 KK=1,KDIAG2
               IF (IK.EQ.NNLEV(KK)) THEN
                  DO 480 JJ=1,JDIAG2
                  JM=NNLAT(JJ)
                  DO 480 IIIL=1,ILMM
                     DD(IIIL,JJ,KK,9)=VFLD(IIIL,JM,IK) *
     1                FLUX(IIIL,JM)
 480              CONTINUE
               END IF
 490        CONTINUE
         END IF
 
        DO 520 IJ=1,IJMMM1
           IM = IJ - 1
           IP = IJ + 1
           DO 500 IIIL=1,ILMM
              DNEW(IIIL,IJ) =
     1                    -(DTOVDXJ/COSLAT(IJ))*
     1                  (VFLD(IIIL,IP,IK)*FLUX(IIIL,IP)*COSLAT2(IP)
     1                  -VFLD(IIIL,IJ,IK)*FLUX(IIIL,IJ)*COSLAT2(IJ))
 500       CONTINUE
 520    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 540 IIIL=1,ILMM
            SUMSP = SUMSP +  VFLD(IIIL,IJ,IK)   *FLUX(IIIL,IJ)
            SUMNP = SUMNP +  VFLD(IIIL,IJMM,IK)*FLUX(IIIL,IJMM)
 540     CONTINUE
         DO 560 IIIL=1,ILMM
            DNEW(IIIL,0)   =  -POLCON*SUMSP
            DNEW(IIIL,IJMM) =  POLCON*SUMNP
 560     CONTINUE
 
C        RELOAD ARRAY
         DO 580 IJ=0,IJMM
         DO 580 IIIL=1,ILMM
            F1(IIIL,IJ,IK) = DNEW(IIIL,IJ)
 580     CONTINUE
 600  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE YADV2(COSLAT,COSLAT2,VFLD,DPHI,F1,PS1INV,PSV,TDT,
     1 POLEDGE)
      IMPLICIT REAL(A-H,O-Z)
 
      PARAMETER (ILMM=144,IJMM=88,IJMMM1=IJMM-1)
      PARAMETER (IL=ILMM,IJMMP1=IJMM+1,ILIJ=ILMM*IJMMP1)
      PARAMETER (IL2=IL/2,IL2P=IL2+1,KTROP=8,KSTRAT=11)
      PARAMETER (IKMM=KTROP+KSTRAT)
 
C     NORTH-SOUTH ADVECTION
      REAL RDNEW(ILMM)
      REAL COSLAT(0:IJMM)
      REAL VFLD(ILMM,IJMM,IKMM),PS1INV(ILMM,0:IJMM)
      REAL F1(ILMM,0:IJMM,IKMM),COSLAT2(IJMM)
      REAL PSV(ILMM,IJMM)
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      COMMON //  DNEW(ILMM,0:IJMM)
 
C     INITIALIZE CONSTANTS
      DTOVDXJ = TDT / DPHI
      A = 6.37E6
      PI=4.0*ATAN(1.)
      DPHI2 = (90. + POLEDGE) * PI / 180.
      POLCON = 2.* TDT/(A*ILMM*SIN(DPHI2))
 
      DDD = 1. / ILMM
      DO 100 IK=1,KTROP
         DO 60 IJ = 1,IJMMM1
            IM = IJ - 1
            IP = IJ + 1
            TERM1 = 1.0
            TERM2 = 1.0
            DO 50 IIIL=1,IL
               DNEW(IIIL,IJ)=-(DTOVDXJ/COSLAT(IJ))*PS1INV(IIIL,IJ)*
     1         (VFLD(IIIL,IP,IK)*TERM1*PSV(IIIL,IP)*COSLAT2(IP)-
     1          VFLD(IIIL,IJ,IK)*TERM2*PSV(IIIL,IJ)*COSLAT2(IJ))
 50         CONTINUE
 60     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
        SUM = 0.
        DO 70 IIIL=1,IL
           SUM = SUM + VFLD(IIIL,IJ,IK)
 70     CONTINUE
        DO 80 IIIL=1,IL
           DNEW(IIIL,0) = -POLCON*SUM
 80     CONTINUE
        IJ = IJMMM1
        IP = IJMM
        SUM = 0.
        DO 85 IIIL=1,IL
           SUM = SUM + VFLD(IIIL,IP,IK)
 85     CONTINUE
        DO 90 IIIL=1,IL
           DNEW(IIIL,IJMM) = POLCON*SUM
 90     CONTINUE
 
C       RELOAD ARRAY
        DO 95 IJ=0,IJMM
        DO 95 IIIL=1,ILMM
           F1(IIIL,IJ,IK) = DNEW(IIIL,IJ)
 95     CONTINUE
 100  CONTINUE
 
C     BEGIN CALCULATIONS FOR STRATOSPHERE.
      DO 200 IK=KTROP+1,IKMM
 
C       FLUX CALCULATIONS.
        DO 160 IJ = 1,IJMMM1
           IM = IJ - 1
           IP = IJ + 1
           TERM1 = 1.0
           TERM2 = 1.0
           DO 140 IIIL=1,IL
              DNEW(IIIL,IJ) =   -(DTOVDXJ/COSLAT(IJ)) *
     1         (VFLD(IIIL,IP,IK)*TERM1*COSLAT2(IP)
     1         -VFLD(IIIL,IJ,IK)*TERM2*COSLAT2(IJ))
 140       CONTINUE
 160    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
        SUM = 0.
        DO 170 IIIL=1,IL
           SUM = SUM + VFLD(IIIL,IJ,IK)
 170    CONTINUE
        DO 180 IIIL=1,IL
           DNEW(IIIL,0) = -POLCON*SUM
 180    CONTINUE
        IJ = IJMMM1
        IP = IJMM
        SUM = 0.
        DO 185 IIIL=1,IL
           SUM = SUM + VFLD(IIIL,IP,IK)
 185    CONTINUE
        DO 190 IIIL=1,IL
           DNEW(IIIL,IJMM) = POLCON*SUM
 190    CONTINUE
 
C       RELOAD ARRAY
        DO 195 IJ=0,IJMM
        DO 195 IIIL=1,IL
           F1(IIIL,IJ,IK) = DNEW(IIIL,IJ)
 195    CONTINUE
 200  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE RDOV19(ISKIP,ILUN,LUN,ISK)
      IMPLICIT REAL(A-H,O-Z)
C
C     CREATED FROM SKELETON OF RDOV19 ON SEPTEMBER 12, 1988
C     LAST REVISION: 9/27/90
C
      PARAMETER (NTAPES=61,KTOP=19)
C
      INTEGER ILUN(NTAPES)
C
C     TAPE INVENTORY
C     TAPE 01 10 OBS (1-10)         TAPE 02 10 OBS (11-20)
C     TAPE 03 11 OBS (21-31)        TAPE 04 11 OBS (32-42)
C     TAPE 05 10 OBS (43-52)        TAPE 06 11 OBS (53-63)
C     TAPE 07 10 OBS (64-73)        TAPE 08 11 OBS (74-84)
C     TAPE 09 10 OBS (85-94)        TAPE 10 10 OBS (95-104)
C     TAPE 11 10 OBS (105-114)      TAPE 12 13 OBS (115-127)
C     TAPE 13 11  OBS (128-138)     TAPE 14 12 OBS (139-150)
C     TAPE 15  2  OBS (151)         TAPE 16 11 OBS (152-162)
C     TAPE 17 12  OBS (163-174)     TAPE 18 11 OBS (175-185)
C     TAPE 19 11  OBS (186-196)     TAPE 20 11 OBS (197-207)
C     TAPE 21 11  OBS (208-218)     TAPE 22 12 OBS (219-230)
C     TAPE 23 11  OBS (231-241)     TAPE 24  5 OBS (242-246)
C     TAPE 25  4  OBS (247-250)     TAPE 26  9 OBS (251-259)
C     TAPE 27 11  OBS (260-270)     TAPE 28 10 OBS (271-280)
C     TAPE 29  9  OBS (281-289)     TAPE 30 10 OBS (290-299)
C     TAPE 31  6  OBS (300-305)     TAPE 32  7 OBS (306-312)
C     TAPE 33 10  OBS (313-322)     TAPE 34 10 OBS (323-332)
C     TAPE 35 11  OBS (333-343)     TAPE 36 11 OBS (344-354)
C     TAPE 37  9  OBS (355-363)     TAPE 38 12 OBS (364-375)
C     TAPE 39 13  OBS (376-388)     TAPE 40 13 OBS (389-401)
C     TAPE 41  8  OBS (402-409)     TAPE 42 15 OBS (410-424)
C     TAPE 43 15  OBS (425-439)     TAPE 44 16 OBS (440-455)
C     TAPE 45 16  OBS (456-471)     TAPE 46 16 OBS (472-487)
C     TAPE 47 16  OBS (488-503)     TAPE 48 16 OBS (504-519)
C     TAPE 49 16  OBS (520-535)     TAPE 50 16 OBS (536-551)
C     TAPE 51 16  OBS (552-567)     TAPE 52 18 OBS (568-585)
C     TAPE 53 15  OBS (586-600)     TAPE 54 14 OBS (601-614)
C     TAPE 55 13  OBS (615-627)     TAPE 56 14 OBS (628-641)
C     TAPE 57 13  OBS (642-654)     TAPE 58 13 OBS (655-667)
C     TAPE 59 14  OBS (668-681)     TAPE 60 13 OBS (682-694)
C     TAPE 61  9  OBS (695-703)
C
C
C     NVAR VARIABLES ON KLEV LEVELS, A HEADER AND PSFC EXIST EACH OBS
C     (NOTICE NVAR=4 AS SPECIFIC HUMIDITY IS NOT USED.
      NVAR = 3
C
      NRECS = 2  + NVAR * KTOP
C
      IF (ISKIP.LT.10) THEN
         LUN  = ILUN (1)
         ISK  = ISKIP
      ELSE IF (ISKIP.LT.20) THEN
         LUN  = ILUN (2)
         ISK  = ISKIP - 10
      ELSE IF (ISKIP.LT.31) THEN
         LUN  = ILUN (3)
         ISK  = ISKIP - 20
      ELSE IF (ISKIP.LT.42) THEN
         LUN  = ILUN (4)
         ISK  = ISKIP - 31
      ELSE IF (ISKIP.LT.52) THEN
         LUN  = ILUN (5)
         ISK  = ISKIP - 42
      ELSE IF (ISKIP.LT.63) THEN
         LUN = ILUN(6)
         ISK = ISKIP - 52
      ELSE IF (ISKIP.LT.73) THEN
         LUN  = ILUN (7)
         ISK  = ISKIP - 63
      ELSE IF (ISKIP.LT.84) THEN
         LUN  = ILUN (8)
         ISK  = ISKIP - 73
      ELSE IF (ISKIP.LT.94) THEN
         LUN  = ILUN (9)
         ISK  = ISKIP - 84
      ELSE IF (ISKIP.LT.104) THEN
         LUN = ILUN(10)
         ISK = ISKIP - 94
      ELSE IF (ISKIP.LT.114) THEN
         LUN = ILUN(11)
         ISK = ISKIP - 104
      ELSE IF (ISKIP.LT.127) THEN
         LUN = ILUN(12)
         ISK = ISKIP - 114
      ELSE IF (ISKIP.LT.138) THEN
         LUN  = ILUN (13)
         ISK  = ISKIP - 127
      ELSE IF (ISKIP.LT.150) THEN
         LUN  = ILUN (14)
         ISK  = ISKIP - 138
      ELSE IF (ISKIP.LT.151) THEN
         LUN  = ILUN (15)
         ISK  = ISKIP - 150
      ELSE IF (ISKIP.LT.162) THEN
         LUN  = ILUN (16)
         ISK  = ISKIP - 151
      ELSE IF (ISKIP.LT.174) THEN
         LUN  = ILUN (17)
         ISK  = ISKIP - 162
      ELSE IF (ISKIP.LT.185) THEN
         LUN = ILUN(18)
         ISK = ISKIP - 174
      ELSE IF (ISKIP.LT.196) THEN
         LUN  = ILUN (19)
         ISK  = ISKIP - 185
      ELSE IF (ISKIP.LT.207) THEN
         LUN  = ILUN (20)
         ISK  = ISKIP - 196
      ELSE IF (ISKIP.LT.218) THEN
         LUN  = ILUN (21)
         ISK  = ISKIP - 207
      ELSE IF (ISKIP.LT.230) THEN
         LUN  = ILUN (22)
         ISK  = ISKIP - 218
      ELSE IF (ISKIP.LT.241) THEN
         LUN = ILUN(23)
         ISK = ISKIP - 230
      ELSE IF (ISKIP.LT.246) THEN
         LUN = ILUN(24)
         ISK = ISKIP - 241
      ELSE IF (ISKIP.LT.250) THEN
         LUN = ILUN(25)
         ISK = ISKIP - 246
      ELSE IF (ISKIP.LT.259) THEN
         LUN = ILUN(26)
         ISK = ISKIP - 250
      ELSE IF (ISKIP.LT.270) THEN
         LUN  = ILUN (27)
         ISK  = ISKIP - 259
      ELSE IF (ISKIP.LT.280) THEN
         LUN  = ILUN (28)
         ISK  = ISKIP - 270
      ELSE IF (ISKIP.LT.289) THEN
         LUN  = ILUN (29)
         ISK  = ISKIP - 280
      ELSE IF (ISKIP.LT.299) THEN
         LUN  = ILUN (30)
         ISK  = ISKIP - 289
      ELSE IF (ISKIP.LT.305) THEN
         LUN = ILUN(31)
         ISK = ISKIP - 299
      ELSE IF (ISKIP.LT.312) THEN
         LUN = ILUN(32)
         ISK = ISKIP - 305
      ELSE IF (ISKIP.LT.322) THEN
         LUN = ILUN(33)
         ISK = ISKIP - 312
      ELSE IF (ISKIP.LT.332) THEN
         LUN  = ILUN (34)
         ISK  = ISKIP - 322
      ELSE IF (ISKIP.LT.343) THEN
         LUN  = ILUN (35)
         ISK  = ISKIP - 332
      ELSE IF (ISKIP.LT.354) THEN
         LUN  = ILUN (36)
         ISK  = ISKIP - 343
      ELSE IF (ISKIP.LT.363) THEN
         LUN  = ILUN (37)
         ISK  = ISKIP - 354
      ELSE IF (ISKIP.LT.375) THEN
         LUN = ILUN(38)
         ISK = ISKIP - 363
      ELSE IF (ISKIP.LT.388) THEN
         LUN = ILUN(39)
         ISK = ISKIP - 375
      ELSE IF (ISKIP.LT.401) THEN
         LUN = ILUN(40)
         ISK = ISKIP - 388
      ELSE IF (ISKIP.LT.409) THEN
         LUN  = ILUN (41)
         ISK  = ISKIP - 401
      ELSE IF (ISKIP.LT.424) THEN
         LUN  = ILUN (42)
         ISK  = ISKIP - 409
      ELSE IF (ISKIP.LT.439) THEN
         LUN  = ILUN (43)
         ISK  = ISKIP - 424
      ELSE IF (ISKIP.LT.455) THEN
         LUN  = ILUN (44)
         ISK  = ISKIP - 439
      ELSE IF (ISKIP.LT.471) THEN
         LUN  = ILUN (45)
         ISK  = ISKIP - 455
      ELSE IF (ISKIP.LT.487) THEN
         LUN  = ILUN (46)
         ISK  = ISKIP - 471
      ELSE IF (ISKIP.LT.503) THEN
         LUN = ILUN(47)
         ISK = ISKIP-487
      ELSE IF (ISKIP.LT.519) THEN
         LUN = ILUN(48)
         ISK = ISKIP-503
      ELSE IF (ISKIP.LT.535) THEN
         LUN = ILUN(49)
         ISK = ISKIP-519
      ELSE IF (ISKIP.LT.551) THEN
         LUN = ILUN(50)
         ISK = ISKIP-535
      ELSE IF (ISKIP.LT.567) THEN
         LUN = ILUN(51)
         ISK = ISKIP-551
      ELSE IF (ISKIP.LT.585) THEN
         LUN  = ILUN (52)
         ISK  = ISKIP - 567
      ELSE IF (ISKIP.LT.600) THEN
         LUN  = ILUN (53)
         ISK  = ISKIP - 585
      ELSE IF (ISKIP.LT.614) THEN
         LUN  = ILUN (54)
         ISK  = ISKIP - 600
      ELSE IF (ISKIP.LT.627) THEN
         LUN  = ILUN (55)
         ISK  = ISKIP - 614
      ELSE IF (ISKIP.LT.641) THEN
         LUN  = ILUN (56)
         ISK  = ISKIP - 627
      ELSE IF (ISKIP.LT.654) THEN
         LUN = ILUN(57)
         ISK = ISKIP-641
      ELSE IF (ISKIP.LT.667) THEN
         LUN = ILUN(58)
         ISK = ISKIP-654
      ELSE IF (ISKIP.LT.681) THEN
         LUN = ILUN(59)
         ISK = ISKIP-667
      ELSE IF (ISKIP.LT.694) THEN
         LUN = ILUN(60)
         ISK = ISKIP-681
      ELSE IF (ISKIP.LT.703) THEN
         LUN = ILUN(61)
         ISK = ISKIP-694
      ELSE
         STOP 'ISK LG'
      END IF
C
      DO 20 I = 1, ISK
      DO 20 IREC = 1,NRECS
         READ(LUN) AAA
 20   CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE SAVERP(CON1,CZ,CZZ,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)
      REAL CZ  (N,0:JGMAX,KTOP,NCON),CZZ(N,0:JGMAX,KTOP,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     THE FIRST MOMENT OF CONSTITUENT IS SAVED (LAYER BY LAYER) AND
C     CONSTITUENT  BY CONSTITUENT.
      DO 50 IC = 1, NCON
      DO 50 KK = 1, KTOP
         DO 40 IJ=0,JGMAX
         DO 40 IL=1,N
            ACYB2(IL,IJ) = CZ(IL,IJ,KK,IC)
  40     CONTINUE
         WRITE (30) ACYB2
         IF (DAY) WRITE (31) ACYB2
 50   CONTINUE
 
C     THE SECOND MOMENT OF CONSTITUENT IS SAVED (LAY BY LAY) AND
C     CONSTITUENT  BY CONSTITUENT.
      DO  70 IC = 1, NCON
      DO  70 KK = 1, KTOP
         DO  60 IJ=0,JGMAX
         DO  60 IL=1,N
            ACYB2(IL,IJ) = CZZ(IL,IJ,KK,IC)
  60     CONTINUE
         WRITE (30) ACYB2
         IF (DAY) WRITE (31) ACYB2
 70   CONTINUE
 
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 SAVERPB(CON1,CZ,CZZ,FD0,TIME,IREAD,HNO3ZM)
      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)
      REAL CZ  (N,0:JGMAX,KTOP,NCON),CZZ(N,0:JGMAX,KTOP,NCON)
      REAL HNO3ZM(0:JGMAX,KTOP)
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     THE FIRST MOMENT OF CONSTITUENT IS SAVED (LAYER BY LAYER) AND
C     CONSTITUENT  BY CONSTITUENT.
      DO 50 IC = 1, NCON
      DO 50 KK = 1, KTOP
         DO 40 IJ=0,JGMAX
         DO 40 IL=1,N
            ACYB2(IL,IJ) = CZ(IL,IJ,KK,IC)
  40     CONTINUE
         WRITE (30) ACYB2
         IF (DAY) WRITE (31) ACYB2
 50   CONTINUE
 
C     THE SECOND MOMENT OF CONSTITUENT IS SAVED (LAY BY LAY) AND
C     CONSTITUENT  BY CONSTITUENT.
      DO  70 IC = 1, NCON
      DO  70 KK = 1, KTOP
         DO  60 IJ=0,JGMAX
         DO  60 IL=1,N
            ACYB2(IL,IJ) = CZZ(IL,IJ,KK,IC)
  60     CONTINUE
         WRITE (30) ACYB2
         IF (DAY) WRITE (31) ACYB2
 70   CONTINUE
 
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
 
      WRITE (30) HNO3ZM
      IF (DAY) WRITE (31) HNO3ZM
 
      RETURN
      END
C
C
C
      SUBROUTINE RESTARP(CON1,CZ,CZZ,FD0,TIME,IREAD,ISAVE)
      IMPLICIT REAL (A-H,O-Z)
C
C     NOTE: 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))
C     PARAMETER (N1=N+1)
      REAL CON1(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON)
      REAL CZ(N,0:JGMAX,KTOP,NCON), CZZ(N,0:JGMAX,KTOP,NCON)
C
      COMMON // ACYB2(N,0:JGMAX)
C
C     OBTAIN DESIRED RESTART DATA
      MRECS = 3*NCON*KTOP  + 3
      LUNIN = 44
 
      DO 20 IOBS = 1,ISAVE
      DO 20 ISKIP= 1,MRECS
         READ(LUNIN) 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 (LUNIN) ACYB2
         DO 100 IJ=0,JGMAX
         DO 100 IL=1,N
            CON1(IL,IJ,KK,IC) = ACYB2(IL,IJ)
 100  CONTINUE
 
C     THE 1ST MOM IS RD IN (LAYER BY LAYER) AND TRACER BY TRACER.
      DO 120 IC=1,NCON
      DO 120 KK = 1,KTOP
         READ (LUNIN) ACYB2
         DO 120 IJ=0,JGMAX
         DO 120 IL=1,N
            CZ(IL,IJ,KK,IC) = ACYB2(IL,IJ)
 120  CONTINUE
 
C     THE 2ND MOM IS RD IN (LAYER BY LAYER) AND TRACER BY TRACER.
      DO 140 IC=1,NCON
      DO 140 KK = 1,KTOP
         READ (LUNIN) ACYB2
         DO 140 IJ=0,JGMAX
         DO 140 IL=1,N
            CZZ(IL,IJ,KK,IC) = ACYB2(IL,IJ)
 140  CONTINUE
 
C     MIXING RATIO FOR LEVEL ZERO IS READ IN.
      READ (LUNIN) FD0
 
C     THE TIME INTO THE INTEGRATION AND IREAD IS READ IN
      READ (LUNIN) TIME
***************SPECIAL SECTION FOR READING (ETA TO MVS) RESTARTS
C      READ (LUNIN) TREAD
C      IREAD = TREAD+0.01
*****************END SPECIAL SECTION
      READ (LUNIN) IREAD
 
      RETURN
      END
C
C
C
      SUBROUTINE RESTARPB(CON1,CZ,CZZ,FD0,TIME,IREAD,ISAVE,HNO3ZM)
      IMPLICIT REAL (A-H,O-Z)
C
C     NOTE: 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))
C     PARAMETER (N1=N+1)
      REAL CON1(N,0:JGMAX,KTOP,NCON),FD0(N,0:JGMAX,NCON)
      REAL CZ(N,0:JGMAX,KTOP,NCON), CZZ(N,0:JGMAX,KTOP,NCON)
      REAL HNO3ZM(0:JGMAX,KTOP)
C
      COMMON // ACYB2(N,0:JGMAX)
C
C     OBTAIN DESIRED RESTART DATA
      MRECS = 3*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     THE 1ST MOM IS RD IN (LAYER BY LAYER) AND TRACER BY TRACER.
      DO 120 IC=1,NCON
      DO 120 KK = 1,KTOP
         READ (40) ACYB2
         DO 120 IJ=0,JGMAX
         DO 120 IL=1,N
            CZ(IL,IJ,KK,IC) = ACYB2(IL,IJ)
 120  CONTINUE
 
C     THE 2ND MOM IS RD IN (LAYER BY LAYER) AND TRACER BY TRACER.
      DO 140 IC=1,NCON
      DO 140 KK = 1,KTOP
         READ (40) ACYB2
         DO 140 IJ=0,JGMAX
         DO 140 IL=1,N
            CZZ(IL,IJ,KK,IC) = ACYB2(IL,IJ)
 140  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
      READ (40) HNO3ZM
 
      RETURN
      END
C
C
C
      SUBROUTINE WIND(U,V,COSLAT,B,DG,PS,KEL)
      IMPLICIT REAL (A-H,O-Z)
 
C     VELOCITY FIELDS ARE DEFINED ON P. 37 OF WILLIAMSON + RASCH'S 1988
C     PAPER "TWO-DIMENSIONAL SEMI-LAGRANGIAN TRANSPORT WITH SPACE
C     PRESERVING INTERPOLATION".
 
      PARAMETER (N=144,JGMAX=88,KTOP=19,NJG=N*JGMAX,NJG2=N*(JGMAX+1))
 
      REAL COSLAT(0:JGMAX)
      REAL U(N,0:JGMAX,KTOP), V(N,JGMAX,KTOP),
     1 DL(N), DG(0:JGMAX),U2(N,0:JGMAX),V2(N,JGMAX), DL2(N)
      REAL PS(N,0:JGMAX),KEL(N,0:JGMAX,KTOP)
 
      DO 5 IJ=0,JGMAX
      DO 5 IL=1,N
         PS(IL,IJ) = 1000.
 5    CONTINUE
 
      DO 10 K=1,KTOP
      DO 10 IJ=0,JGMAX
      DO 10 IL=1,N
         KEL(IL,IJ,K) = 273.
 10   CONTINUE
 
      A = 6.37E6
      PI = 4.0 * ATAN(1.)
      D2 = (360. / N)
      F0 =  -D2/2.
 
      DTR = PI / 180.
      DO 12 IL=1,N
         DL(IL) = (IL-1)*D2 * DTR
         DL2(IL) = (F0 + (IL-1)*D2) * DTR
 12   CONTINUE
C     U1 IS DEFINED AS PI/128 RADIANS PER TIME STEP
      U1 = 2. * PI * A  / (256. * 900.)
C     SPECIAL TEST
C      U1 = U1 * .66666666666666666666667
C
      J = 0
      DO 15 I=1,N
        U2(I,J) =U1*(COS(B)*COSLAT(J)+SIN(B)*SIN(DG(J))*COS(DL2(I)))
 15   CONTINUE
 
      DO 20 J=1,JGMAX
      DO 20 I=1,N
        U2(I,J) =U1*(COS(B)*COSLAT(J)+SIN(B)*SIN(DG(J))*COS(DL2(I)))
        V2(I,J) =-U1 * SIN(B) * SIN(DL(I))
  20  CONTINUE
 
      DO 40 K=1,KTOP
      DO 40 IJ=0,JGMAX
      DO 40 IL=1,N
         U(IL,IJ,K) = U2(IL,IJ)
 40   CONTINUE
 
      DO 50 K=1,KTOP
      DO 50 IJ=1,JGMAX
      DO 50 IL=1,N
         V(IL,IJ,K) = V2(IL,IJ)
 50   CONTINUE
      RETURN
      END
