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



      SUBROUTINE SETVER
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER(KSTRAT=11,KTROP=8,KTOP=KTROP+KSTRAT)
 
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),
     2 SIGMA(KTOP), PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),
     3 SGINT(0:KTOP)
 
      PTOP      = 0.305
      PINT      = 120.
      PTOP2     = LOG(PTOP)
      PINT2     = LOG(PINT)
      DIV       = (PTOP2 - PINT2) / KSTRAT
 
      DO 10 K=1,KSTRAT
         PRESS(K) =  EXP(PINT2+(K-1)*DIV+DIV*0.5)
 10   CONTINUE
 
      PRSLAY(0)      = PINT
      PRSLAY(KSTRAT) = PTOP
      DO 15 K=1,KSTRAT-1
         PRSLAY(K) = SQRT(PRESS(K)*PRESS(K+1))
 15   CONTINUE
 
C     DEFINE THE SIGMA LAYERS.
      DIV2 = 1. / KTROP
      DO 20 K=1,KTROP
         SIGMA(K) = (1.-DIV2/2.) - DIV2 * (K-1)
 20   CONTINUE
 
      DO 30 K=0,KTROP
         SGINT(K) = 1.0 - (K * DIV2)
 30   CONTINUE
      DO 35 K=KTROP+1,KTOP
         SGINT(K) = (PRSLAY(K-KTROP)-PINT) / (PINT-PTOP)
 35   CONTINUE
 
      DO 40 K= KTROP+1,KTOP
         SIGMA(K) = (PRESS(K-KTROP) - PINT) / (PINT-PTOP)
 40   CONTINUE
 
      DO 60 K=1,KSTRAT
         DELZ(K) = 7000. * LOG(PRESS(KSTRAT-1)/PRESS(KSTRAT))
 60   CONTINUE
 
      DSIG(0) = SIGMA(1) - 1.
      DO 80 K=1,KTOP-1
         DSIG(K) = SIGMA(K+1) -SIGMA(K   )
 80   CONTINUE
      DSIG(KTOP) =  -1.-SIGMA(KTOP)
C
      DO 100 K=1,KTOP
         DELTAS(K) = SGINT(K)-SGINT(K-1)
 100  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE PL3D(FIELD,PSF,IC,LUN)
      IMPLICIT REAL (A-H,O-Z)
C
C     DATA WILL BE TRANSFERRED FOR PLOTTING
C
      PARAMETER (N=144,JGMAX=88,JGMAX1=JGMAX+1)
      PARAMETER (KTROP=8,KSTRAT=11,KTOP=KTROP+KSTRAT)
      PARAMETER (JGMAX3=JGMAX1+2)
C
      REAL FIELD (N,0:JGMAX,KTOP), PSF(N,0:JGMAX)
      REAL ACYB(N,JGMAX3)
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
      IF (IC.EQ.1) THEN
         DO 20 IL=1,N
            ACYB(IL,1) = PSF(IL,0)
            ACYB(IL,2) = PSF(IL,0)
            ACYB(IL,JGMAX3-1)= PSF(IL,JGMAX)
            ACYB(IL,JGMAX3  )= PSF(IL,JGMAX)
 20      CONTINUE
 
         DO 40 J=3,JGMAX3-2
         J2 = J-2
         DO 40 IL=1,N
            ACYB(IL,J) = PSF(IL,J2)
 40      CONTINUE
 
         write (lun) acyb
      END IF
 
      DO 100 KK = 1,KTOP
         DO 60 IL=1,N
            ACYB(IL,1) = FIELD(IL,0,KK)
            ACYB(IL,2) = FIELD(IL,0,KK)
            ACYB(IL,JGMAX3-1)= FIELD(IL,JGMAX,KK)
            ACYB(IL,JGMAX3  )= FIELD(IL,JGMAX,KK)
 60      CONTINUE
 
         DO 80 J=3,JGMAX3-2
         J2 = J-2
         DO 80 IL=1,N
            ACYB(IL,J) = FIELD(IL,J2,KK)
 80      CONTINUE
 
         write (lun) acyb
 100  CONTINUE
C
      RETURN
      END
C
C
C
      SUBROUTINE PLXD(FIELD,LUN)
      IMPLICIT REAL (A-H,O-Z)
C
      PARAMETER (N=144,JGMAX=88,KTOP=19)
C
      REAL FIELD (N,JGMAX,KTOP)
C
      COMMON // ACYB(N,JGMAX)
C
      DO 100 KK = 1,KTOP
         DO 20 IJ=1,JGMAX
         DO 20 IL=1,N
            ACYB(IL,IJ) = FIELD(IL,IJ,KK)
 20      CONTINUE
         write (lun) acyb
 100  CONTINUE
C
      RETURN
      END
C
C
C
      SUBROUTINE PL2D(FIELD,LUN)
      IMPLICIT REAL (A-H,O-Z)
C
C     DATA WILL BE TRANSFERRED TO THE V07 FOR PLOTTING
C
      PARAMETER (N=144,JGMAX=88,JGMAX1=JGMAX+1)
C
      REAL FIELD (N,0:JGMAX), AIBM(N,JGMAX1)
 
      DO 20 J=1,JGMAX1
      DO 20 I=1,N
         AIBM(I,J) = FIELD(I,J)
 20   CONTINUE
C
      WRITE (LUN) AIBM
C
      RETURN
      END
C
C
C
      SUBROUTINE SURFPS(PS1,PSU,PSV)
      IMPLICIT REAL (A-H,O-Z)
C
C     PSU(1,J) IS AVG OF PS1(ILMM,J) AND PS1(1,J)
 
      PARAMETER (ILMM=144,IJMM=88)
 
      REAL PS1(ILMM,0:IJMM),PSU(ILMM,0:IJMM),PSV(ILMM,IJMM)
 
C     PSU IS DEFINED.
      DO 20 J=0,IJMM
      DO 20 I=1,ILMM
         IM=I-1
         IF (IM.LT.1) IM = ILMM
         PSU(I,J) = 0.5 * (PS1(IM,J) + PS1(I,J))
 20   CONTINUE
 
C     PSV IS DEFINED.
      DO 40 J=1,IJMM
         JM1 = J-1
         DO 30 IL=1,ILMM
            PSV(IL,J) = 0.5 * (PS1(IL,JM1)+PS1(IL,J))
 30      CONTINUE
 40   CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE RDOV1989(ISKIP,ILUN,LUN,ISK)
      IMPLICIT REAL (A-H,O-Z)
C
C     CREATED FROM SKELETON OF RDOV19 ON SEPTEMBER 12, 1988
C     REVISED 07 JUNE 1990.
C
      PARAMETER (NTAPES=61,KTOP=19)
C
      INTEGER ILUN(NTAPES)
C
C     BEGIN
C
C     TAPE 1  8   OBS (1-8)     TAPE18  12 OBS (213-224)
C     TAPE 2  10  OBS (9-18)    TAPE19  13 OBS (225-237)
C     TAPE 3  11  OBS (19-29)   TAPE20  15 OBS (238-252)
C     TAPE 4  13  OBS (30-42)   TAPE21  16 OBS (253-268)
C     TAPE 5  14  OBS (43-56)   TAPE22  12 OBS (269-280)
C     TAPE 6  13  OBS (57-69)   TAPE23  13 OBS (281-293)
C     TAPE 7  12  OBS (70-81)   TAPE24  13 OBS (294-306)
C     TAPE 8  14  OBS (82-95)   TAPE25  13 OBS (307-319)
C     TAPE 9  13  OBS (96-108)  TAPE26  14 OBS (320-333)
C     TAPE 10 12  OBS (109-120) TAPE27  13 OBS (334-346)
C     TAPE 11 13  OBS (121-133) TAPE28  14 OBS (347-360)
C     TAPE 12 12  OBS (134-145) TAPE29  14 OBS (361-374)
C     TAPE 13 15  OBS (146-160) TAPE30   2 OBS (375-376)
C     TAPE 14 13  OBS (161-173)
C     TAPE 15 13  OBS (174-186)
C     TAPE 16 13  OBS (187-199)
C     TAPE 17 13  OBS (200-212)
C
C     NVAR VARIABLES ON KLEV LEVELS, A HEADER AND PSFC EXIST EACH OBS
C     (NOTICE NVAR=3 AS SPECIFIC HUMIDITY IS NOT USED IN VVV DECK)
      NVAR = 3
C
      NRECS = 2  + NVAR * KTOP
C
      IF (ISKIP.LT.8) THEN
         LUN  = ILUN (1)
C        2 IS ADDED SINCE EXP BEGINS 00Z DEC 28, BUT DATA BEG 12Z DEC 27
         ISK  = ISKIP+2
      ELSE IF (ISKIP.LT.18) THEN
         LUN  = ILUN (2)
         ISK  = ISKIP - 8
      ELSE IF (ISKIP.LT.29) THEN
         LUN  = ILUN (3)
         ISK  = ISKIP - 18
      ELSE IF (ISKIP.LT.42) THEN
         LUN  = ILUN (4)
         ISK  = ISKIP - 29
      ELSE IF (ISKIP.LT.56) THEN
         LUN  = ILUN (5)
         ISK  = ISKIP - 42
      ELSE IF (ISKIP.LT.69) THEN
         LUN = ILUN(6)
         ISK = ISKIP - 56
      ELSE IF (ISKIP.LT.81) THEN
         LUN  = ILUN (7)
         ISK  = ISKIP - 69
      ELSE IF (ISKIP.LT.95) THEN
         LUN  = ILUN (8)
         ISK  = ISKIP - 81
      ELSE IF (ISKIP.LT.108) THEN
         LUN  = ILUN (9)
         ISK  = ISKIP - 95
      ELSE IF (ISKIP.LT.120) THEN
         LUN  = ILUN (10)
         ISK  = ISKIP - 108
      ELSE IF (ISKIP.LT.133) THEN
         LUN = ILUN(11)
         ISK = ISKIP - 120
      ELSE IF (ISKIP.LT.145) THEN
         LUN = ILUN(12)
         ISK = ISKIP - 133
      ELSE IF (ISKIP.LT.160) THEN
         LUN = ILUN(13)
         ISK = ISKIP - 145
      ELSE IF (ISKIP.LT.173) THEN
         LUN  = ILUN (14)
         ISK  = ISKIP - 160
      ELSE IF (ISKIP.LT.186) THEN
         LUN  = ILUN (15)
         ISK  = ISKIP - 173
      ELSE IF (ISKIP.LT.199) THEN
         LUN  = ILUN (16)
         ISK  = ISKIP - 186
      ELSE IF (ISKIP.LT.212) THEN
         LUN  = ILUN (17)
         ISK  = ISKIP - 199
      ELSE IF (ISKIP.LT.224) THEN
         LUN  = ILUN (18)
         ISK  = ISKIP - 212
      ELSE IF (ISKIP.LT.237) THEN
         LUN  = ILUN (19)
         ISK  = ISKIP - 224
      ELSE IF (ISKIP.LT.252) THEN
         LUN  = ILUN (20)
         ISK  = ISKIP - 237
      ELSE IF (ISKIP.LT.268) THEN
         LUN  = ILUN (21)
         ISK  = ISKIP - 252
      ELSE IF (ISKIP.LT.280) THEN
         LUN  = ILUN (22)
         ISK  = ISKIP - 268
      ELSE IF (ISKIP.LT.293) THEN
         LUN  = ILUN (23)
         ISK  = ISKIP - 280
      ELSE IF (ISKIP.LT.306) THEN
         LUN  = ILUN (24)
         ISK  = ISKIP - 293
      ELSE IF (ISKIP.LT.319) THEN
         LUN  = ILUN (25)
         ISK  = ISKIP - 306
      ELSE IF (ISKIP.LT.333) THEN
         LUN  = ILUN (26)
         ISK  = ISKIP - 319
      ELSE IF (ISKIP.LT.346) THEN
         LUN  = ILUN (27)
         ISK  = ISKIP - 333
      ELSE IF (ISKIP.LT.360) THEN
         LUN = ILUN(28)
         ISK = ISKIP - 346
      ELSE IF (ISKIP.LT.374) THEN
         LUN = ILUN(29)
         ISK = ISKIP - 360
      ELSE IF (ISKIP.LT.376) THEN
         LUN = ILUN(30)
         ISK = ISKIP - 374
      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 NEGS(FIELD,IC,TIME)
      IMPLICIT REAL (A-H,O-Z)
 
C     THE NUMBER OF GRID PTS WITH NEGATIVE MIXING RATIOS IS DETERMINED
C     AND BINNED BY LATITUDE AND LEVEL.
 
      PARAMETER (N=144,JGMAX=88,KTOP=19)
 
      REAL FIELD(N,0:JGMAX,KTOP)
 
      INTEGER  NEG1(0:JGMAX), NEG2(KTOP)
 
      NEG3 = 0
      DO 20 K=1,KTOP
         NEG2(K) = 0
 20   CONTINUE
      DO 40 J=0,JGMAX
         NEG1(J) = 0
 40   CONTINUE
 
      T2 = TIME / 86400.
 
      DO 60 K=1,KTOP
      DO 60 J=0,JGMAX
      DO 60 I=1,N
         IF (FIELD(I,J,K).LT.0.) THEN
            NEG3    = NEG3    + 1
            NEG2(K) = NEG2(K) + 1
            NEG1(J) = NEG1(J) + 1
         END IF
  60  CONTINUE
 
      IF (NEG3.NE.0) THEN
         WRITE (6,29) T2
  29     FORMAT(1X,' MODEL TIME IN DAYS ',F6.2)
         WRITE (6,30) IC,NEG3
  30     FORMAT (1X,'CONSTITUENT ',I3,' .LT. ZERO ',I6,' TIMES ')
         WRITE (6,50) (NEG2(K),K=1,KTOP)
  50     FORMAT (1X,19I6)
         WRITE (6,70) (NEG1(J),J=0,JGMAX)
  70     FORMAT(1X,10I6)
      END IF
      RETURN
      END
C
C
C
      SUBROUTINE INTERZ3D(FD1,N,JGMAX,KIN,FD2,KOUT,XIN,XOUT,
     1 YIN,YOUT)
C
      IMPLICIT REAL (A-H,O-Z)
 
      REAL FD1 (N,JGMAX,KIN),FD2(N,JGMAX,KOUT)
      REAL XIN(KIN), XOUT(KOUT),YIN(KIN),YOUT(KOUT)
C
      DO  80 J=1,JGMAX
         DO  60 I=1,N
            DO  20 K=1,KIN
               YIN(K) = FD1(I,J,K)
  20        CONTINUE
            CALL INTERP(XIN,YIN,KIN,XOUT,YOUT,KOUT)
            DO  40 K=1,KOUT
               FD2(I,J,K) = YOUT(K)
  40        CONTINUE
  60     CONTINUE
  80  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE INTERY3D(FD1,N,JMP,KOUT,FD2,JGMAX,XIN,XOUT,YIN,
     1 YOUT)
C
      IMPLICIT REAL (A-H,O-Z)
 
      REAL FD1 (N,JMP,KOUT),FD2(N,JGMAX,KOUT)
      REAL XIN(JMP), XOUT(JGMAX),YIN(JMP),YOUT(JGMAX)
C
      DO  80 K=1,KOUT
         DO  60 I=1,N
            DO  20 J=1,JMP
               YIN(J) = FD1(I,J,K)
  20        CONTINUE
            CALL INTERP(XIN,YIN,JMP,XOUT,YOUT,JGMAX)
            DO  40 J=1,JGMAX
               FD2(I,J,K) = YOUT(J)
  40        CONTINUE
  60     CONTINUE
  80  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE INTERX3D(FD1,IMP,JGMAX,KOUT,FD2,N,XIN,XOUT,YIN,YOUT)
C
      IMPLICIT REAL (A-H,O-Z)
 
      REAL FD1 (IMP,JGMAX,KOUT),FD2(N,JGMAX,KOUT)
      REAL XIN(IMP), XOUT(N),YIN(IMP),YOUT(N)
C
      DO  80 K=1,KOUT
         DO  60 J=1,JGMAX
            DO  20 I=1,IMP
               YIN(I) = FD1(I,J,K)
  20        CONTINUE
            CALL INTERP(XIN,YIN,IMP,XOUT,YOUT,N)
            DO  40 I=1,N
               FD2(I,J,K) = YOUT(I)
  40        CONTINUE
  60     CONTINUE
  80  CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE INTERY2D(FD1,N,JMP,FD2,JGMAX,XIN,XOUT,YIN,YOUT)
C
      IMPLICIT REAL (A-H,O-Z)
 
      REAL FD1 (N,JMP),FD2(N,JGMAX)
      REAL XIN(JMP), XOUT(JGMAX),YIN(JMP),YOUT(JGMAX)
C
         DO  60 I=1,N
            DO  20 J=1,JMP
               YIN(J) = FD1(I,J)
  20        CONTINUE
            CALL INTERP(XIN,YIN,JMP,XOUT,YOUT,JGMAX)
            DO  40 J=1,JGMAX
               FD2(I,J) = YOUT(J)
  40        CONTINUE
  60     CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE INTERX2D(FD1,IMP,JGMAX,FD2,N,XIN,XOUT,
     1 YIN,YOUT)
C
      IMPLICIT REAL (A-H,O-Z)
 
      REAL FD1 (IMP,JGMAX),FD2(N,JGMAX)
      REAL XIN(IMP), XOUT(N),YIN(IMP),YOUT(N)
C
         DO  60 J=1,JGMAX
            DO  20 I=1,IMP
               YIN(I) = FD1(I,J)
  20        CONTINUE
            CALL INTERP(XIN,YIN,IMP,XOUT,YOUT,N)
            DO  40 I=1,N
               FD2(I,J) = YOUT(I)
  40        CONTINUE
  60     CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE INTER2(XIN,YIN,N,XOUT,YOUT,M)
C
      IMPLICIT REAL (A-H,O-Z)
 
C    GIVEN A SET OF INPUT COORDINATES, XIN AND YIN, AND A SET OF
C    DESIRED ABSCISSA VALUES, XOUT, WHICH ARE NOT, IN GENERAL,
C    EQUAL TO XIN, AT WHICH YOU WANT TO KNOW THE CORRESPONDING
C    ORDINATES, YOUT, THIS SUBROUTINE WILL DO A THIRD ORDER POLYNOMIAL
C    FIT WHILE XOUT LIES WITHIN THE RANGE OF XIN, AND A LINEAR
C    EXTRAPOLATION IF XOUT LIES OUTSIDE THE RANGE OF XIN.
C
C    XIN AND YIN HAVE DIMENSION N, XOUT AND YOUT HAVE DIMENSION M.
C
C    GENEROUSLY DONATED BY MARK SCHOEBERL
C
      DIMENSION XIN(N),YIN(N),XOUT(M),YOUT(M),Q(4),R(4),S(3),T(2),E(3)
C
      DO 90 J=1,M
      IF(XOUT(J) .GE. XIN(N)) GO TO 70
      IF(XOUT(J) .LE. XIN(1)) GO TO 80
      IFLAG=0
      DO 20 K=1,N
      IF(IFLAG .EQ. 1) GO TO 20
      IF(XOUT(J) .LE. XIN(K)) IFLAG=1
      KSAVE=K
   20 CONTINUE
      LQ=KSAVE-2
      IF(LQ .LE. 1) LQ=1
      IR=N-3
      IF(LQ .GE. IR) LQ=N-3
      DO 30 L=1,4
      LN=L+LQ-1
      Q(L)=XIN(LN)
      R(L)=YIN(LN)
   30 CONTINUE
      DO 40 L=1,3
      S(L)=(R(L+1)-R(L))/(Q(L+1)-Q(L))
   40 CONTINUE
      DO 50 L=1,2
      T(L)=(S(L+1)-S(L))/(Q(L+2)-Q(L))
   50 CONTINUE
      D=(T(2)-T(1))/(Q(4)-Q(1))
      DO 60 L=1,3
      E(L)=XOUT(J)-Q(L)
   60 CONTINUE
      YOUT(J)=R(1)+S(1)*E(1)+T(1)*E(1)*E(2)+D*E(1)*E(2)*E(3)
      GO TO 90
   70 YOUT(J)=(YIN(N)-YIN(N-1))/(XIN(N)-XIN(N-1))*
     +        (XOUT(J)-XIN(N-1))+YIN(N-1)
      GO TO 90
   80 YOUT(J)=(YIN(2)-YIN(1))/(XIN(2)-XIN(1))*
     +        (XOUT(J)-XIN(1))+YIN(1)
   90 CONTINUE
      RETURN
 
      END
C
C
C
      SUBROUTINE INTERP2(XIN,YIN,N,XOUT,YOUT,M)
C
      IMPLICIT REAL (A-H,O-Z)
 
C     XIN  INPUT LEVELS       XOUT   OUTPUT LEVELS
C     YIN  INPUT VALUES       YOUT   OUTPUT VALUES
C     N    # OF INPUT LEVS    M      # OF OUTPUT LEVS
C     LINEAR INTERPOLATION IS PERFORMED FOR
C                       XIN(1) .LT. XOUT(J). LT. XIN(N)
C
      REAL XIN(N),YIN(N),XOUT(M),YOUT(M)
C
      DO 120 J=1,M
         IF (XOUT(J).LE.XIN(1)) THEN
C           YOUT(J)=(YIN(2)-YIN(1))/(XIN(2)-XIN(1))*
C    +       (XOUT(J)-XIN(1))+YIN(1)
            YOUT(J) = YIN(1)
 
         ELSE IF (XOUT(J).GE.XIN(N)) THEN
C           YOUT(J)=(YIN(N)-YIN(N-1))/(XIN(N)-XIN(N-1))*
C    +       (XOUT(J)-XIN(N-1))+YIN(N-1)
            YOUT(J) = YIN(N)
 
         ELSE
            DO 80 JJ=2,N
               JJJ = JJ
               IF (XOUT(J).LE.XIN(JJ)) GOTO 100
  80        CONTINUE
 100        YOUT(J) = (YIN(JJJ)-YIN(JJJ-1))/(XIN(JJJ)-XIN(JJJ-1))*
     !       (XOUT(J)-XIN(JJJ-1)) + YIN(JJJ-1)
         END IF
  120 CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE INTERP(XIN,YIN,N,XOUT,YOUT,M)
C
      IMPLICIT REAL (A-H,O-Z)
 
C     XIN  INPUT LEVELS       XOUT   OUTPUT LEVELS
C     YIN  INPUT VALUES       YOUT   OUTPUT VALUES
C     N    # OF INPUT LEVS    M      # OF OUTPUT LEVS
C     LINEAR INTERPOLATION IS PERFORMED FOR
C                       XIN(1) .LT. XOUT(J). LT. XIN(N)
C
      REAL XIN(N),YIN(N),XOUT(M),YOUT(M)
      INTEGER LOC(M)
C
      DO 20 J=1,M
         LOC(J) = ISRCHFGT(N,XIN,1,XOUT(J))
 20   CONTINUE
 
      DO 30 J=1,M
         IF (LOC(J).EQ.1) THEN
            YOUT(J) = YIN(1)
         ELSE IF (LOC(J).EQ.(N+1)) THEN
            YOUT(J) = YIN(N)
         ELSE
            YOUT(J)=(YIN(LOC(J))-YIN(LOC(J)-1))/
     !       (XIN(LOC(J))-XIN(LOC(J)-1))*
     !       (XOUT(J)-XIN(LOC(J)-1)) + YIN(LOC(J)-1)
         END IF
  30  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE RDSIM5(LUN1,PS,U,V,KEL)
      IMPLICIT REAL (A-H,O-Z)
C **********************************************************************
C *   PS,U,V, AND KEL ARE READ IN                                      *
C **********************************************************************
C ----------------------------------------------------------------------
C     LAST UPDATE 6/21/90
C
      PARAMETER (IMP=72,JMP=45,JGG=90)
      PARAMETER (N=144,JGMAX=88,KTOP=19,KSTRAT=11)
      PARAMETER (IMJM1=IMP*(JMP+1),IMJM=IMP*JMP)
C
      REAL U(N,0:JGMAX,KTOP), V(N,JGMAX,KTOP),
     1 PS(N,0:JGMAX),KEL(N,0:JGMAX,KTOP)
 
      LOGICAL CYC
C
      REAL HDR(3), FD1(IMP,0:JMP), FIELD1(IMP,JMP)
 
      COMMON // HDR2(3), DLATRA0(JMP),DLATRA1(JGG),DLON0B(IMP+1),
     1 DLATRB0(0:JMP),DLATRB1(0:JGG),DLON0(IMP),DLON1(N),
     2 FD2(IMP,0:JMP),FD3(N,0:JMP),FD4(N,0:JGG),XIN(IMP+1),XOUT(N),
     3 FIELD2(IMP,JMP),FIELD3(N,JMP),FIELD4(N,JGG),F(IMP+1,0:JMP),
     4 YIN(0:JMP),YOUT(0:JGG),YIN1(JMP),YOUT1(JGG),F2(IMP+1,JMP),
     5 DLON0J(IMP+1),DLON0K(N)
 
      WRITE (6,*) ' LUN1 IN RDSIM5 = ', LUN1
 
C     READ IN HDR
c     READ (LUN1) HDR
c     CALL Q9IC32(HDR,HDR2,3,IERR)
       read (lun1) hdr2
      WRITE (6,11) HDR2(1), HDR2(2), HDR2(3)
 11   FORMAT(1X,3F12.2)
 
      D0 = 180. / JMP
      D1 = 180. / JGG
      D2 = 360. / IMP
      D3 = 360. / N
      DO 20 J=1,JMP
         DLATRA0(J) = -88.0 + (J-1)*D0
 20   CONTINUE
      DO 30 J=1,JGG
         DLATRA1(J) = -89.0 + (J-1)*D1
 30   CONTINUE
      DO 40 J=0,JMP
         DLATRB0(J) = -90.0 + J*D0
 40   CONTINUE
      DO 50 J=0,JGG
         DLATRB1(J) = -90.0 + J*D1
 50   CONTINUE
      DO 60 I=1,IMP
         DLON0(I) =   0.0 + (I-1)*D2
         DLON0J(I) = -2.5 + (I-1)*D2
 60   CONTINUE
      DO 65 IL=1,IMP
         DLON0B(IL) = DLON0(IL)
 65   CONTINUE
      DLON0B(IMP+1) = 360.
      DLON0J(IMP+1) = 360. + DLON0J(1)
      DO 70 I=1,N
         DLON1(I) =   0.0 + (I-1)*D3
         DLON0K(I) = -1.25 + (I-1)*D3
 70   CONTINUE
******************************
c     READ (LUN1) FD1
       read (lun1) fd2
      DO 72 J=0,JMP
         DO 71 IL=1,IMP
            F(IL,J) = FD2(IL,J)
 71      CONTINUE
         F(IMP+1,J    ) = FD2(1,J    )
 72   CONTINUE
C     WRITE (6,*) ' F BEFORE X INTERPOLATION'
C     WRITE (6,2230) (F(I,10),I=1,IMP)
 2230 FORMAT(10E12.5)
      CALL INTERX2D(F  ,IMP+1,JMP+1,FD3,N   ,DLON0B ,DLON1  ,XIN,XOUT)
C     WRITE (6,*) 'FD3 AFTER X INTERPOLATION '
C     WRITE (6,2230) (FD3(I,10),I=1,N)
C     WRITE (6,*) 'FD3 BEFORE Y INTERPOLATION '
C     WRITE (6,2230) (FD3(1,J),J=0,JMP)
      CALL INTERY2D(FD3,N   ,JMP+1,FD4,JGG+1,DLATRB0,DLATRB1,YIN,YOUT)
C     WRITE (6,*) 'FD4 AFTER Y INTERPOLATION '
C     WRITE (6,2230) (FD4(1,J),J=0,JGG)
*********************************
      SUMSP = 0.
      SUMNP = 0.
      DO 80 IL=1,N
         SUMSP = SUMSP + FD4(IL,0)
         SUMNP = SUMNP + FD4(IL,JGG)
 80   CONTINUE
      SUMSP = SUMSP / N
      SUMNP = SUMNP / N
 
      DO 90 IL=1,N
         PS(IL,   0)  =  SUMSP
         PS(IL,JGMAX) =  SUMNP
 90   CONTINUE
 
      DO 100 IJ=1,JGMAX-1
         IJ1 = IJ + 1
         DO 95 IL=1,N
            PS(IL,IJ) =   FD4(IL,IJ1)
 95      CONTINUE
 100  CONTINUE
 
      DO 160 K=1,KTOP
******************************
c        READ (LUN1) FD1
c        CALL Q9IC32(FD1,FD2,IMJM1,IERR)
         read (lun1) fd2
         DO 115 J=0,JMP
            DO 112 IL=1,IMP
               F(IL    ,J) = FD2(IL,J)
 112        CONTINUE
            F(IMP+1,J    ) = FD2(1,J    )
 115     CONTINUE
         CALL INTERX2D(F,IMP+1,JMP+1,FD3,N ,DLON0J ,DLON0K ,XIN,XOUT)
         CALL INTERY2D(FD3,N,JMP+1,FD4,JGG+1,DLATRB0,DLATRB1,YIN,YOUT)
*********************************
 
         SUMSP = 0.
         SUMNP = 0.
         DO 120 IL=1,N
            SUMSP = SUMSP + FD4(IL,0)   / N
            SUMNP = SUMNP + FD4(IL,JGG) / N
 120     CONTINUE
 
         DO 130 IL=1,N
            U(IL,0,K)     = SUMSP
            U(IL,JGMAX,K) = SUMNP
 130     CONTINUE
 
         DO 140 IJ=1,JGMAX-1
         IJ1 = IJ + 1
         DO 140 IL=1,N
            U(IL,IJ,K) =   FD4(IL,IJ1)
 140     CONTINUE
 160  CONTINUE
C
C     PROCESS V WIND FIELD
      DO 180 K=1,KTOP
******************************
c        READ(LUN1) FIELD1
c        CALL Q9IC32(FIELD1,FIELD2,IMJM,IERR)
         read (lun1) field2
         DO 165 J=1,JMP
            DO 162 IL=1,IMP
               F2(IL    ,J) = FIELD2(IL,J)
 162        CONTINUE
            F2(IMP+1,J    ) = FIELD2(1,J    )
 165     CONTINUE
         CALL INTERX2D(F2,IMP+1,JMP,FIELD3,N,DLON0B ,DLON1  ,XIN,XOUT)
         CALL INTERY2D(FIELD3,N,JMP,FIELD4,JGG,DLATRA0,DLATRA1,YIN1,
     1    YOUT1)
*********************************
         DO 170 IJ=1,JGMAX
         IJ1 = IJ + 1
         DO 170 IL=1,N
            V(IL,IJ,K) =  FIELD4(IL,IJ1)
 170     CONTINUE
 180   CONTINUE
C
      DO 240 K=1,KTOP
******************************
c        READ (LUN1) FD1
c        CALL Q9IC32(FD1,FD2,IMJM1,IERR)
            read (lun1) fd2
         DO 190 J=0,JMP
            DO 185 IL=1,IMP
               F(IL    ,J) = FD2(IL,J)
 185        CONTINUE
            F(IMP+1,J    ) = FD2(1,J    )
 190     CONTINUE
         CALL INTERX2D(F,IMP+1,JMP+1,FD3,N ,DLON0B ,DLON1  ,XIN,XOUT)
         CALL INTERY2D(FD3,N,JMP+1,FD4,JGG+1,DLATRB0,DLATRB1,YIN,YOUT)
*********************************
         SUMSP = 0.
         SUMNP = 0.
         DO 200 IL=1,N
            SUMSP = SUMSP + FD4(IL,0)   / N
            SUMNP = SUMNP + FD4(IL,JGG) / N
 200     CONTINUE
 
         DO 210 IL=1,N
            KEL(IL,0    ,K) = SUMSP
            KEL(IL,JGMAX,K) = SUMNP
 210     CONTINUE
 
         DO 220 IJ=1,JGMAX-1
         IJ1 = IJ + 1
         DO 220 IL=1,N
            KEL(IL,IJ,K) =   FD4(IL,IJ1)
 220     CONTINUE
 240  CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE TINIT(CONST,LUN1,SCALE,FD0,IC)
      IMPLICIT REAL (A-H,O-Z)
C **********************************************************************
C *   INITIAL CONDITION FOR NOZ TRANSPORT. 3-DIMENSIONAL.  6/89        *
C     INTERPOLATED ONTO LEER GRID  FOR TEST RUN 11/89                  *
C **********************************************************************
C ----------------------------------------------------------------------
C
      PARAMETER (N=144,JGG=90,IMP=96,JG2=78,KTROP=8,JGMAX=88)
      PARAMETER (KSTRAT=11,KTOP=KTROP+KSTRAT)
      PARAMETER (IMJM=IMP*JG2)
C
      REAL CONST(N,0:JGMAX,KTOP),DGLAT3(JG2)
C
      REAL FIELD(IMP,JG2),FD0(N,0:JGMAX)
 
      COMMON //  FD(IMP,JG2),DLON1(IMP),DLON2(N),
     1 DLAT(0:JGG  ),FD3(N,JG2,KTOP),FD4(N,0:JGG,KTOP),DLON10(IMP+1),
     2 F(IMP+1,JG2,KTOP),XIN(IMP+1),XOUT(N),YIN(JG2),YOUT(0:JGG)
C
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),
     2 SIG(KTOP),P(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:KTOP)
 
      COMMON/THREE/HTSCALE,HINV,PI
 
      DATA DGLAT3/ -88.245, -85.971, -83.684, -81.394, -79.102, -76.810,
     !    -74.518, -72.226, -69.933, -67.640, -65.348, -63.055, -60.762,
     !    -58.469, -56.177, -53.884, -51.591, -49.298, -47.005, -44.712,
     !    -42.419, -40.126, -37.833, -35.541, -33.248, -30.955, -28.662,
     !    -26.369, -24.076, -21.783, -19.490, -17.197, -14.904, -12.611,
     !    -10.318,  -8.025,  -5.732,  -3.439,  -1.146,   1.146,   3.439,
     !      5.732,   8.025,  10.318,  12.611,  14.904,  17.197,  19.490,
     !     21.783,  24.076,  26.369,  28.662,  30.955,  33.248,  35.541,
     !     37.833,  40.126,  42.419,  44.712,  47.005,  49.298,  51.591,
     !     53.884,  56.177,  58.469,  60.762,  63.055,  65.348,  67.640,
     !     69.933,  72.226,  74.518,  76.810,  79.102,  81.394,  83.684,
     !     85.971,  88.245/
 
      REWIND (LUN1)
C     INTERPOLATE.
      D1 = 360./IMP
      D2 = 360./N
      DO 5 IL=1,IMP
         DLON1(IL) = 0.0 + (IL-1)*D1
 5    CONTINUE
      DO 7 IL=1,IMP
         DLON10(IL) = DLON1(IL)
 7    CONTINUE
      DLON10(IMP+1) = DLON1(1) + 360.
      DO 10 IL=1,N
         DLON2(IL) = 0.0 + (IL-1)*D2
 10   CONTINUE
      DY = 180. / JGG
      DO 15 IJ=0,JGG
         DLAT(IJ) = -90.+ IJ*DY
 15   CONTINUE
 
C     READ IN PSF
      READ (LUN1) AAA
C     READ IN INITIAL CONDITION
      DO 20 K=1,KTOP
 
          read (lun1) fd
         DO 18 J=1,JG2
            DO 16 IL=1,IMP
               F(IL,J  ,K) = FD(IL,J)
 16         CONTINUE
            F(IMP+1,J,K)      = FD(1,J)
 18      CONTINUE
 20   CONTINUE
 
C     INTERPOLATE ONTO TRANSPORT GRID.
C     WRITE (6,*) ' F BEFORE X INT'
*************************
C     WRITE (6,2299) (F(II,60,17),II=1,IMP+1)
C2299 FORMAT(10E12.5)
      CALL INTERX3D(F,IMP+1,JG2,KTOP,FD3,N,DLON10,DLON2,XIN,XOUT)
C     WRITE (6,*) ' F AFT X INT'
C     WRITE (6,2299) (FD3(II,60,17),II=1,N)
C     WRITE (6,*) ' F BEFORE Y INTERPOLATION'
C     WRITE (6,2230) (FD3(1,IJ,17),IJ=1,JG2)
 2230 FORMAT(10E12.5)
      CALL INTERY3D(FD3,N,JG2,KTOP,FD4,JGG+1,DGLAT3,DLAT,YIN,YOUT)
C     WRITE (6,*) ' F AFT Y INT'
C     WRITE (6,2299) (FD4(1,IJ,17),IJ=0,JGG)
 
*************************
      DO 25 K=1,KTOP
         DO 22 IL=1,N
            CONST(IL,0    ,K) = FD4(IL,0  ,K)
            CONST(IL,JGMAX,K) = FD4(IL,JGG,K)
 22      CONTINUE
         DO 24 IJ=1,JGMAX-1
         IJ1 = IJ + 1
         DO 24 IL=1,N
            CONST(IL,IJ,K) = FD4(IL,IJ1,K)
 24      CONTINUE
 25   CONTINUE
 
*******SPECIAL SECTION NOY RUN ONLY***************
      IF (IC.EQ.2) THEN
         DO 30 K=1,KTOP
            DO 27 IJ=0,JGMAX
            DO 27 IL=1,N
               CONST(IL,IJ,K) = 22.4 - 0.074*CONST(IL,IJ,K)
 27         CONTINUE
            DO 28 IJ=0,JGMAX
            DO 28 IL=1,N
               IF (CONST(IL,IJ,K).LT.0.0) CONST(IL,IJ,K) = 0.
               IF (CONST(IL,IJ,K).GT.20.0) CONST(IL,IJ,K) = 20.
 28         CONTINUE
 30      CONTINUE
      END IF
**************************************************
 
C     NOTE: CONVERSION FROM PPB TO MIXING RATIO
C     AVERAGE TO DETERMINE MIXING RATIO IN POLAR BOXES.
      DO 40 K=1,KTOP
         DO 35 J=0,JGMAX
         DO 35 I=1,N
            CONST(I,J,    K) = CONST(I,J,K) * SCALE
 35      CONTINUE
         SUMSP = 0.
         SUMNP = 0.
         DO 38 IL=1,N
            SUMSP = SUMSP + CONST(IL,0    ,K)
            SUMNP = SUMNP + CONST(IL,JGMAX,K)
 38      CONTINUE
         DO 39 IL=1,N
            CONST(IL,0    ,K) = SUMSP / N
            CONST(IL,JGMAX,K) = SUMNP / N
 39      CONTINUE
 40   CONTINUE
 
      DO 50 IJ=0,JGMAX
      DO 50 IL=1,N
         FD0(IL,IJ) =  CONST(IL,IJ,1)
 50   CONTINUE
 
C     WRITE (6,*) ' CONST AT END OF SUB'
C     WRITE (6,2299) (CONST(II,60,17),II=1,N)
      RETURN
      END
C
C
C
      SUBROUTINE PATTERN (COSCEN,FIELD,FD0)
      IMPLICIT REAL (A-H,O-Z)
 
      PARAMETER (N=144,JGMAX=88,KTOP=19)
 
      REAL COSCEN(0:JGMAX),FD0(N,0:JGMAX)
      REAL FIELD(N,0:JGMAX,KTOP),FD(N,0:JGMAX),DLONG(N)
 
      PI = 4.0 * ATAN(1.)
      D2 = 2. * PI / N
      DO 10 IL=1,N
         DLONG(IL) = (IL-1)* D2
 10   CONTINUE
      N2 = N / 2
      RCAP = 14. *  PI / 128
      WRITE (6,*) ' ABOUT TO BEGIN 15 LOOP'
      DO 15 IJ=0,JGMAX
      DO 15 IL=1,N
         FD(IL,IJ) = 0.
 15   CONTINUE
      DO 20 J=0,JGMAX
      DO 20 I=N2+1,N
         R = ACOS(COS(DLONG(I)-1.5*PI)*COSCEN(J))
         IF (R.LT.RCAP) FD(I,J) = 0.5*(1.+COS(PI*R/RCAP))
 20   CONTINUE
      WRITE (6,*) ' AFTER 20 LOOP '
 
      DO 40 K=5,KTOP
      DO 40 J=0,JGMAX
      DO 40 I=1,N
         FIELD(I,J,K) = FD(I,J)
C        FIELD(I,J,K) = 10.
 40   CONTINUE
 
      DO 60 IJ=0,JGMAX
      DO 60 IL=1,N
         FD0(IL,IJ) =  FIELD(IL,IJ,1)
 60   CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE MASS3(CON,PS1,MCOR,LUN)
 
      PARAMETER (N=144,JGMAX=88,KTOP=19,KTROP=8,KSTRAT=11,NCON=1)
 
      REAL MCOR(N,0:JGMAX),MASS(N,0:JGMAX,KTOP)
      REAL CON(N,0:JGMAX,KTOP,NCON)
      REAL PS1(N,0:JGMAX), DMAS(KTOP,NCON), DMTOTL(NCON)
 
      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
C     REQUIRED MIXING RATIO CHANGE IS CONVERTED TO TERM PROPORTIONAL TO
C     REQUIRED MASS CHANGE BY MULTIPLYING BY DX * DY * DP
      DO 10 K=1,KTROP
      DO 10 J=0,JGMAX
      DO 10 I=1,N
         MASS(I,J,K) = (PS1(I,J)-PINT) *(SGINT(K-1)-SGINT(K))
 10   CONTINUE
 
      DO 30 K=KTROP+1,KTOP
         K2 = K - KTROP
         DO 20 J=0,JGMAX
         DO 20 I=1,N
             MASS(I,J,K) =  PRSLAY(K2-1) - PRSLAY(K2)
 20      CONTINUE
 30   CONTINUE
 
      DO 50 K=1,KTOP
      DO 50 J=0,JGMAX
      DO 50 I=1,N
         MASS(I,J,K) = 100. * MCOR(I,J) * MASS(I,J,K)/ G
 50   CONTINUE
 
      DO 160 IC = 1,NCON
         DMTOTL(IC) = 0.0
         DO 140 K=1,KTOP
            SUM = 0.
            DO 130 J=0,JGMAX
            FAC = 1.
            IF ((J.EQ.0.).OR.(J.EQ.JGMAX)) FAC=1./N
            DO 130 I=1,N
               SUM = SUM + MASS(I,J,K) * CON(I,J,K,IC) * FAC
 130        CONTINUE
            DMAS(K,IC) = SUM
            DMTOTL(IC) = DMTOTL(IC) + SUM
 140     CONTINUE
C
         WRITE (LUN,11) IC,DMTOTL(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',E15.5)
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 ',(8E13.5))
      END
C
C
C
      SUBROUTINE WDINT4 (U1,U2,V1,V2,
     ! PS1,PS2,KEL1,KEL2,TIME,TWIND1,TWIND2,PSNEW,TIME2)
      IMPLICIT REAL (A-H,O-Z)
C
      PARAMETER (KSTRAT=11)
      PARAMETER (N=144,JGMAX=88,KTOP=19)
 
      COMMON / VERT/ DELTAS(KTOP),DSIG(0:KTOP),DELZ(KSTRAT),
     2 SIGMA(KTOP),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),
     3 SGINT(0:KTOP)
 
C     TIME INTERPOLATION IS PERFORMED ON THE INPUT FIELDS.
C     THE TIME INTERPOLATED FIELD IS RETURNED IN U1,V1,DIV1,PS1,+ KEL1
C
C
      REAL PS1(N,0:JGMAX), PS2(N,0:JGMAX),PSNEW(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
      TNEW = (TIME2-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)
         PSNEW(I,J) =  (1.-TNEW)*PS1(I,J) + TNEW*PS2(I,J)
 20   CONTINUE
C
      RETURN
      END
C
C
C
      SUBROUTINE INITCFC(CONST,FD0)
 
C     CONST    INITIAL CONDITION FOR F-11 IN PPV.
 
      PARAMETER (N=144,JGMAX=88,IKMM=19,KSTRAT=11,KTROP=8,KTOP=19)
 
      REAL CONST(N,0:JGMAX,KTOP),FD0(N,0:JGMAX)
 
      COMMON / VERT/ DELTAS(IKMM),DSIG(0:IKMM),DELZ(KSTRAT),
     1 SIG(IKMM),PRESS(KSTRAT),PTOP,PINT,PRSLAY(0:KSTRAT),SGINT(0:IKMM)
 
      J2 = 33
      DO 20 K=1,KTROP
      DO 20 J=0,J2
      DO 20 I=1,N
         CONST(I,J,K) = 220 * 1.E-12
 20   CONTINUE
 
      J1 = 34
      J2 = 54
      DO 30 K=1,KTROP
      DO 30 J=J1,J2
      FAC = (J-J1)* 15 / (J2-J1)
      DO 30 I=1,N
         CONST(I,J,K) = (220 + FAC) * 1.E-12
 30   CONTINUE
 
      J1 = 55
      DO 40 K=1,KTROP
      DO 40 J=J1,JGMAX
      DO 40 I=1,N
         CONST(I,J,K) = 235 * 1.E-12
 40   CONTINUE
 
C     STRATOSPHERIC VALUES FALL OFF LINEARLY WITH PRESSURE.
      DO 60 K=KTROP+1,KTOP
      DO 60 J=0,JGMAX
      DO 60 I=1,N
         CONST(I,J,K) = CONST(I,J,KTROP) * PRESS(K-KTROP) / 100.
 60   CONTINUE
 
      DO 80 J=0,JGMAX
      DO 80 I=1,N
         FD0(I,J) = CONST(I,J,1)
 80   CONTINUE
 
      RETURN
      END
C
C
C
      SUBROUTINE INITO3(CONST,LUN1,SCALE,FD0)
 
C **********************************************************************
C *   INITIAL CONDITION FOR 1979 O3 CALCULATION.   4/12/91.            *
C **********************************************************************
C ----------------------------------------------------------------------
      PARAMETER (N=144,JGMAX=88,KTOP=19)
C
      REAL CONST(N,0:JGMAX,KTOP),FD0(N,0:JGMAX)
 
      COMMON // AAA(N,0:JGMAX,KTOP)
 
      REWIND (LUN1)
C     READ IN INITIAL CONDITION
      READ (LUN1)  AAA
 
      DO 50 K=1,KTOP
      DO 50 IJ=0,JGMAX
      DO 50 IL=1,N
         CONST(IL,IJ,K) =  AAA(IL,IJ,K) * SCALE
 50   CONTINUE
 
      DN = 1. / N
      IJ = 0
      DO 70 K=1,KTOP
         SUM = 0.
         DO 60 IL=1,N
            SUM = SUM + CONST(IL,IJ,K) * DN
 60      CONTINUE
 
         DO 65 IL=1,N
            CONST(IL,IJ,K) = SUM
 65      CONTINUE
 70   CONTINUE
 
      IJ = JGMAX
      DO 85 K=1,KTOP
         SUM = 0.
         DO 80 IL=1,N
            SUM = SUM + CONST(IL,IJ,K) * DN
 80      CONTINUE
         DO 82 IL=1,N
            CONST(IL,IJ,K) = SUM
 82      CONTINUE
 85   CONTINUE
 
      DO 90 IJ=0,JGMAX
      DO 90 IL=1,N
         FD0(IL,IJ) = CONST(IL,IJ,1)
 90   CONTINUE
 
      RETURN
      END
