      SUBROUTINE CLOUD6(C,TAUC,  SIGMA,LEVS,LEV,LEVP,XLAT,
     1                 ILAT,IFINE,NLAT,NLATP1,NTYPE,RNDOFF,AMT,BASE,TOP,
     2                 FINEGR,P,Z,TO,IZ,OCCUPY,LBL,W,TC,BNEB, 
     3                 ALAT,OVRLAP,LWC,LWP,LB,LT,LXTND) 
C 
C     * MAY 24/88 - M.LAZARE (BASED ON CLOUDC WITH GENERAL DIMENSIONS). 
C     * MAY 20/83 - R.LAPRISE,J.P.BLANCHET. 
C 
C     * DETERMINE THE MEAN ZONAL AMOUNT OF CLOUD IN  EACH LAYER FOR A GI
C     * HORIZONTAL AND VERTICAL GRID. 
C     *     C(ILAT,1LEV,2LEV)  = NEBULOSITY BETWEEN LAYER 1  AND 2 ,
C     *    EQUIVALENT TO THE MATRIX CCC IN IR CALCULATION WHEN EMISSIVIT
C     *     TAUC(ILAT,LEVS) = RESULTING CLOUD OPTICAL THICKNESS FOR EXTI
C     *     TOTNEB(ILAT) = RESULTING TOTAL NEBULOSITY.
C     *     SIGMA(LEVS) = VECTOR OF SIGMA COORDINATES (GT. 0.). 
C     *     XLAT(ILAT) = VECTOR OF LATITUDE COORDINATES (BETWEEN -90 AND
C     *     RNDOFF = 0, NO MODIFICATION TO ORIGINAL DATA, OTHERWISE CLOU
C     *              ADJUSTED SO, 
C     *                1-OCCUPENCY OF A TYPE IS CONSIDERED IF 1/2 OF THE
C     *                  IS OCCUPIED OR MORE. 
C     *                2-SMALL CLOUDS IN 2LAYERS ARE PLACED IN THE MOST 
C     *                  REPRESENTATIVE LAYER (MOST COVERED)
C     * NOTE... THE PROGRAM USES A FINE GRID (DZ=0.1KM) TO RESOLVE THE
C     *         THE VERTICAL POSITION OF CLOUDS.  IN EACH LAYER THE 
C     *         RANDOM OVERLAP APPROXIMATION IS MADE FOR CLOUDS THAT ARE
C     *         VERTICALLY SEPARATED, OTHERWISE THEIR NEBULOSITY ARE
C     *         ADDED.  THE TOTAL NEBULOSITY OF THE SKY MAY DIFFER FROM 
C     *         THE SUM OF THE NEBULOSITY OF EACH LAYERS DUE TO POSSIBLE
C     *         OVERLAP, THEREFOR THE TOTAL NEBULOSITY IS EVALUATED SEPA
C     *         PROBALITY OF MULTIPLE OVERLAP (MORE THAN 2 CLOUDS) IS 
C     *         IGNORED. THE OPTICAL THICKNESS IS CALCULATED FROM A 
C     *         PRAMETRZATION BY STEPHENS (1978) JAS 35 PP.2123-2124. 
C 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL TAUC(ILAT,LEVS),SIGMA(LEVS),XLAT(ILAT) 
      REAL AMT(NLAT,NTYPE),BASE(NLAT,NTYPE),TOP(NLAT,NTYPE) 
      REAL P(LEVS),Z(LEVS),ALAT(NLATP1),W(LEVS),TC(NLAT,LEVS) 
      REAL OVRLAP(NTYPE,NTYPE),LWC(NTYPE),LWP(NTYPE),TO(LEV)
C 
      INTEGER IZ(LEV),LB(NTYPE),LT(NTYPE),LXTND(NTYPE)
      INTEGER FINEGR(IFINE,NTYPE),RNDOFF,LBL(ILAT),OCCUPY(LEVS,NTYPE) 
C 
C     * LEVEL 2,C 
      REAL C(LEVP,LEVP,ILAT)
C     * LEVEL 2,B 
      REAL BNEB(NLAT,LEVS,LEVS) 
C 
C    **** STEPHENS EQ.10A PP.2124.
C 
      T(X)=10.E0**(0.2633E0+1.7095E0*LOG(LOG10(X))) 
C 
C    **** LMS REGRESSION   WT(T) (FEIGELSON,1977) 
C 
      WT(X)=0.0000435E0*X*X-0.0174E0*X+1.734E0
C-----------------------------------------------------------------------
C 
C    **** HEIGHT SCALE FROM STANDARD US ATMOSPHERE (1962) 
C 
      IZ(1)=110 
      DO 100 I=1,LEVS 
      P(I)=1013.E0*SIGMA(I) 
      Z(I)=-LOG  (P(I)/1018.0E0)/0.1256E0
      IF(P(I).LT.540.5E0) Z(I)=-LOG  (P(I)/1109.2E0)/0.1425E0
      IF(P(I).LT.265.0E0) Z(I)=-LOG  (P(I)/1272.0E0)/0.1565E0
      IF(P(I).LT.121.1E0) Z(I)=-LOG  (P(I)/1214.7E0)/0.1541E0
      IF(SIGMA(I).EQ.1.0E0) Z(I)=0.E0 
      IZ(I+1)=Z(I)*10.E0+0.5E0
      IF(IZ(I+1).GT.110) IZ(I+1)=110
      Z(I)=IZ(I+1)/10.E0
      TO(I)=288.1E0-6.45E0*Z(I) 
      IF(Z(I).GT.11.E0) TO(I)=216.6E0 
  100 CONTINUE
      TO(LEV)=288.1E0 
C 
C    **** DEFINITION OF THE LATITUDE GRID 
C 
      ALAT(1)=-85.E0
      DO 115 I=2,18 
  115 ALAT(I)=ALAT(I-1)+10.E0 
      ALAT(19)=90.0E0 
C 
C    **** CALCULATION OF THE NEBULOSITY OF EACH LAYER.
C 
      DO 180 J=1,18 
C 
C    **** DEFINE A FINER GRID (DZ=0.1KM) BETWEEN 0 (SFC) AND 11KM.
C 
      DO 135 K=1,NTYPE
      DO 125 I=1,110
  125 FINEGR(I,K)=0 
      LBASE=BASE(J,K)*10.E0+1.5E0 
      LTOP=TOP(J,K)*10.E0+0.5E0 
      IF(LTOP.LT.LBASE) LTOP=LBASE
      LB(K)=LBASE 
      LT(K)=LTOP
      LXTND(K)=LTOP-LBASE 
      DO 130 I=LBASE,LTOP 
  130 FINEGR(I,K)=1 
  135 CONTINUE
C 
C    **** CALCULATION OF FRACTIONNAL OVERLAP. 
C         CORRECTION FOR DATA FROM SASAMORI (I.E. NO OVERLAP) . 
C 
      DO 145 K=1,NTYPE
      DO 140 KI=K,NTYPE 
      OVRLAP(K,KI)=0.E0 
      IF(LB(K).GT.LT(KI).OR.LT(K).LT.LB(KI))
     1               OVRLAP(K,KI)=AMT(J,K)*AMT(J,KI)
      OVRLAP(KI,K)=OVRLAP(K,KI) 
  140 CONTINUE
  145 CONTINUE
C 
C    **** DETERMINE THE TYPE OF CLOUD PRESENT IN THE LAYER. 
C 
      DO 157 N=1,LEVS 
C 
      DO 147 K=1,6
  147 OCCUPY(N,K)=0 
C 
      IZ1=IZ(N) 
      IZ2=IZ(N+1)+1 
      IDIF=IZ1-IZ2
      IF(IDIF.LE.0) GO TO 157 
      DO 156 K=1,NTYPE
C 
C        TESTS TO ROUND OFF CLOUD AMOUNT IN THE LAYER.
C 
      IF(RNDOFF.EQ.0) GO TO 152 
      ISUM=0
      DO 151 I=IZ2,IZ1
  151 ISUM=ISUM+FINEGR(I,K) 
      X=IDIF/2.E0 
      Y=LXTND(K)/2.E0 
      SUM=ISUM
      IF(SUM.LT.X.AND.SUM.LT.Y) GO TO 155 
C 
  152 CONTINUE
      DO 153 I=IZ2,IZ1
      IF(FINEGR(I,K).NE.0)  GO TO 154 
  153 CONTINUE
      GO TO 155 
  154 OCCUPY(N,K)=1 
  155 CONTINUE
  156 CONTINUE
  157 CONTINUE
C 
C    **** NEBULOSITY OF THE LAYER.
C 
      DO 165 N=1,LEVS 
      DO 164 M=N,LEVS 
      BNEB(J,N,M)=0.E0
      DO 163 K=1,NTYPE
      DO 158 N1=N,M 
      IF(OCCUPY(N1,K).NE.0) GO TO 159 
  158 CONTINUE
      GO TO 163 
  159 BNEB(J,N,M)=BNEB(J,N,M)+AMT(J,K)
      DO 162 KI=K,NTYPE 
      DO 160 N1=N,M 
      IF(OCCUPY(N1,KI).NE.0) GO TO 161
  160 CONTINUE
      GO TO 162 
  161 BNEB(J,N,M)=BNEB(J,N,M)-OVRLAP(K,KI)
  162 CONTINUE
  163 CONTINUE
      BNEB(J,M,N)=BNEB(J,N,M) 
  164 CONTINUE
  165 CONTINUE
C 
C    **** CALCUL OF THE OPTICAL THICKNESS FOR EXTINCTION BY CLOUDS. 
C 
      DO 178 N=1,LEVS 
      W(N)=0.E0 
      TC(J,N)=0.E0
      IZ1=IZ(N) 
      IZ2=IZ(N+1)+1 
      DO 172 K=1,NTYPE
      LWP(K)=0.E0 
      DO 171 I=IZ2,IZ1
  171 LWP(K)=LWP(K)+FINEGR(I,K)*LWC(K)*100.E0 
      W(N)=W(N)+LWP(K)*AMT(J,K) 
  172 CONTINUE
      IF(BNEB(J,N,N).NE.0.E0) GO TO 174 
      W(N)=0.E0 
      GO TO 175 
  174 W(N)=W(N)/BNEB(J,N,N) 
  175 CONTINUE
C 
C    **** NOTE... WITH THE FOLOWING MODIFICATION, TAUC WILL CONTAIN 
C         LIQUID WATER PATH RATHER THAN CLOUD OPTICAL THICKNESS.
C         CONSEQUENTLY, A TEMPERATURE DEPENDANCE CAN BE INCLUDED IN 
C         PHYSICS.  FOR PTICAL THICKNESS, USE TC(J,N)=T(W(N)) 
C 
      WTAV=(WT(TO(N))+WT(TO(N+1)))/2.E0 
      IF(WTAV.LT.0.02E0) WTAV=0.02E0
      IF(W(N).GT.2.E0) TC(J,N)=W(N) 
  178 CONTINUE
C 
  180 CONTINUE
C 
C    **** LOCATION OF INTERPOLATION POINTS ALONG THE LATITUDE AXIS. 
C 
      DO 210 J=1,ILAT 
      DO 190 I=1,19 
      IBUF=I
      IF(XLAT(J).LE.ALAT(I)) GO TO 193
  190 CONTINUE
  191 WRITE(6,6001) 
      GO TO 999 
  193 LBL(J)=IBUF-1 
      IF(XLAT(J).LE.-85.E0) LBL(J)=LBL(J)+1 
      IF(XLAT(J).GE.+85.E0) LBL(J)=LBL(J)-1 
      IF(LBL(J).GE.1.AND.LBL(J).LE.17) GO TO 200
      WRITE(6,6001) 
      GO TO 999 
  200 CONTINUE
  210 CONTINUE
C 
C    **** LATITUDE INTERPOLATION. 
C 
      DO 250 J=1,ILAT 
      L=LBL(J)
C 
      DO 230 N=1,LEVS 
      IN=N+1
      C(N+1,N+1,J)=0.E0 
C 
      DO 220 II  =IN,LEV
      I=II-1
      B1=BNEB(L,N,I)
      B2=BNEB(L+1,N,I)
      CALL LININT(ALAT(L),B1,ALAT(L+1),B2,XLAT(J),RESULT) 
      RESULT=MAX(0.E0,RESULT) 
      C( N+1,II+1,J) =RESULT
      C(II+1, N+1,J) =RESULT
  220 CONTINUE
      CALL LININT(ALAT(L),TC(L,N),ALAT(L+1),TC(L+1,N),XLAT(J),TAUC(J,N))
      IF(TAUC(J,N).LT.0.E0) TAUC(J,N)=0.0E0 
  230 CONTINUE
      C(LEVP,LEVP,J)=0.E0 
C 
C     * FILL IN MOON LEVEL COLUMN AND ROW.
C 
      C(1,1,J)=0.E0 
      DO 240 L=2,LEVP 
      C(1,L,J)=C(2,L,J) 
  240 C(L,1,J)=C(L,2,J) 
  250 CONTINUE
C 
  999 RETURN
C-----------------------------------------------------------------------
 6001 FORMAT(////,5X,'**** INTERPOLATION ATTEMPTED OUTSIDE THE RANGE OF 
     1THE LATITUDE GRID (-90 TO 90 DEG OF LAT)',////) 
      END 
