      SUBROUTINE CLOUDS3(C MTX, TAC, SHJ, SHTJ, 
     1                   CZON, Q, T, H, SIG HAF, SIG FUL, P SFC, TOPCV, 
     2                   EMICOEF, RADEQV, WCL, WCD, 
     3                   CVEC, HVEC, DMIN, DMAX,
     4                   HAVE, TAVE, DLGSIG,
     5                   ZRDM, CREF, EREF,
     6                   EST, HT, ALWC, XLENGTH,
     7                   H0, E, 
     8                   ILG, IL1, IL2, LEVP1, LEV, ILEV, COSL, 
     9                   MLEVP1, MLEV, MILEV) 

C     * NOV 25/94 - M.LAZARE. COSMETIC OPTIMIZATION FOR CLOUD OVERLAP,
C     *                       LIKE THAT DONE IN GCM10.
C     * JAN 30/89 - M.LAZARE. FINAL, FROZEN VERSION FOR GCM6. LIKE
C     *                       PREVIOUS EXCEPT:  
C     *                       A. IN-LINE FUNCTIONS FOR EST,HTVOCP.
C     *                       B. ICE-PHASE PARAMETERS PASSED IN COMMON
C     *                          BLOCK EPSICE FROM SPWCONH. 
C     *                       C. FRACTAL DIMENSION SET BACK TO 0.8. 
C     *                       D. ICE-PHASE BOUNDARY NOW AT TICE=-40 
C     *                          INSTEAD OF -30.
C     *                       E. UPPER ATMOSPHERE H0 INCREASED FROM 0.5 
C     *                          TO 0.85. 
C     * DEC 05/88 - M.LAZARE. NEW VERSION FOR GCM6. LIKE CLOUDS2 EXCEPT:  
C     *                       A. BETTS, PLATT ET AL USED FOR OPTICAL
C     *                          PROPERTIES (REQUIRES WORK ARRAYS TO
C     *                          VECTORIZE PROPERLY). 
C     *                       B. TOPO FIX REMOVED (DOESN'T PILE UP
C     *                          MOISTURE WITH HYBRID COORDINATE AND
C     *                          NEW MOISTURE VARIABLE).
C     *                       C. LOCAL CLOUDS SCHEME IMPLEMENTED (GJB). 
C     *                       D. LOCAL SIGMA ARRAYS NOW USED (AS
C     *                          REQUIRED BY NEW HYBRID FORMULATION). 
C     ***************** NOTE: THIS SUBROUTINE IS CALLED BY THE CLOUD ***
C     *                       DIAGNOSTIC PROGRAM CLDPROG AND ANY     ***
C     *                       CHANGES IN EITHER MUST BE MUTUTALLY    ***
C     *                       CONSISTENT!!!!!!!!!!!!! 
C     ******************************************************************
  
C     * J.P. BLANCHET. - MODIFIED NOV 13/87 TO REMOVE "BUG" IN DISMET=2 
C     *                  FOR "TOPO FIX".
C 
C     * OUTPUT ARRAYS...
  
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      DIMENSION C MTX(ILG,MLEVP1,MLEVP1), TAC(ILG,MLEV) 
  
C     * INPUT ARRAYS... 
  
      DIMENSION SHJ(ILG,MILEV), SHTJ(ILG,MLEV)
      DIMENSION Q(ILG,MLEV), T(ILG,MLEV), H(ILG,MLEV) 
      DIMENSION C ZON(MLEVP1,MLEVP1)
      DIMENSION SIG HAF(MILEV), SIG FUL(MLEV), P SFC(ILG), TOP CV(ILG)
  
C     * WORK ARRAYS...
  
      DIMENSION EMI COEF(ILG,MLEV), RAD EQV(ILG,MLEV) 
      DIMENSION WCL(ILG,MLEV), WCD(ILG,MLEV)
      DIMENSION C VEC(ILG), H VEC(ILG), D MIN(ILG), DMAX(ILG) 
      DIMENSION H AVE(MLEV), T AVE(MLEV), D LG SIG(MLEVP1)
      DIMENSION EST(ILG), HT(ILG), ALWC(ILG), XLENGTH(ILG)
      DIMENSION H0(ILG), E(ILG) 
      DIMENSION Z RDM(ILG), C REF(ILG), E REF(ILG)
  
      LOGICAL KLD GEN, KLD DIS, KLD LOC, L HIG
  
      COMMON /HTCP  / T1S,T2S,AI,BI,AW,BW,SLP 
      COMMON /PARAMS/ WW,TWW,RAYON,ASQ,GRAV,RGAS,RGOCP,RGOASQ,CPRES 
      COMMON /PARAMS/ RGASV,CPRESV
      COMMON /EPSICE/ AICE,BICE,TICE,QMIN 
      COMMON /EPS/ A, B, EPS 1, EPS 2 
C 
C-----------------------------------------------------------------------
C     * STATEMENT FUNCTION DEFINITIONS. 
  
      TW(TTT)     = AW-BW*TTT 
      TI(TTT)     = AI-BI*TTT 
      HTVOCP(TTT) = MERGE( TW(TTT), 
     1              MERGE( TI(TTT), 
     2                     SLP*((TTT-T2S)*TW(TTT)+(T1S-TTT)*TI(TTT)), 
     3                     TTT.LE.T2S), 
     4                     TTT.GE.T1S)
  
      ESW(TTT)    = EXP(A-B/TTT)
      ESI(TTT)    = EXP(AICE-BICE/TTT)
      ESTEFF(TTT) = MERGE( ESW(TTT),
     1              MERGE( ESI(TTT),
     2                     SLP*((TTT-T2S)*ESW(TTT)+(T1S-TTT)*ESI(TTT)), 
     3                     TTT.LE.T2S), 
     4                     TTT.GE.T1S)
  
C........ PHYSICAL CONSTANTS: 
      DATA R O G /29.29E0/, P REF /101325.E0/, CP O R /3.4983E0/
  
C........ SCHEME PARAMETERS:  
      DATA WC MIN /0.0001E0/, CUT /0.01E0/
  
C........ PARAMETERS AFFECTING THE RESULTS: 
C 
C             FCT LAY O : VERTICAL FRACTION OF THE LAYER OCCUPIED BY CLD
C                         FOR UNRESOLVED LAYER ONLY 
C             SCALE : SCALE DEP. OF THE RH THRESHOLD (N TH ROOT OF COSL)
C                     ACCOUNT FOR REDUCED LONGITUDINAL SCALE LENGTH 
C                     NEAR THE POLE (IE. DC/DH -> INF AT THE POLES).
C             S LIMIT : DC/DH LIMIT (MINIMUM PHYSICAL VALUE = 1)
C                       LIMIT THE COVERAGE OF COLD CLOUDS (CIRRUS & POLES)
C             SHELL : SET LIMIT BETWEEN LOW AND HIGH R.H. REGIMES 
C                     ( 0 < SHELL < 4)
C             F     : FRACTION OF THE VARIANCE DUE TO CLOUDY PORTION. 
C             SCL OBS : SCALE FACTOR FOR OBSERVED LIQUID WATER CONTENT
C                        IN WATER CLOUDS. 
C 
      DATA FCT LAY O /0.5E0/, SCALE /0.33E0/, S LIMIT /1.25E0/,
     &  SHELL /2.E0/, 
     &  F /0.2E0/, SCL OBS /1.0E0/, DIS MET /2/, MATRIX /0/,
     &  WC FAC /10.E0/
  
C........ LOGICAL SWITCH FOR GENERATED / OBSERVED (FIXED) CLOUDINESS
C                 KLD GEN   =  .TRUE.  /  .FALSE. 
  
C........ LOGICAL SWITCH FOR DISTRIBUTED / ZONAL CLOUDINESS.
C                 KLD DIS   =   .TRUE.   / .FALSE.
  
C........ LOGICAL SWICH FOR LOCALY GENERATED CLOUDS (PARTICULAR CASE) 
C          KLD LOC = .TRUE. ==> KLD GEN = .TRUE. AND KLD DIS = .FALSE.
  
  
      DATA KLD GEN /.TRUE./ ,KLD DIS /.FALSE./, KLD LOC /.TRUE./
  
      IF (KLD LOC) THEN 
          KLD GEN = .TRUE.
          KLD DIS = .FALSE. 
      END IF
  
      IF (MATRIX .EQ. 1) THEN 
          KLD GEN = .FALSE. 
          KLD DIS = .FALSE. 
      END IF
  
C-----------------------------------------------------------------------
  
  
C........ INITIALIZE ARRAYS 
  
      LAY = LEV 
      POINTS = 1.E0 / FLOAT (IL2 - IL1 +1)
      I TEST = 1
C     WC FAC = MAX (10. * COS L, 1.)
  
      DO 70 I = 1, LEV P 1
            DO 60 J = 1, LEV P 1
                  DO 50 K = IL1, IL2
                        C MTX(K,J,I) = 0.E0 
   50             CONTINUE
                  IF (KLD GEN) C ZON(I,J) = 0.E0
   60       CONTINUE
   70 CONTINUE
  
      CV AVE = 0.E0 
      DO 78 K = IL1, IL2
          CV AVE = CV AVE + MAX(FLOAT(ILEV) -TOP CV(K), 0.E0) 
          WCL(K,1) = 0.E0 
          WCD(K,1) = 0.E0 
          TAC(K,1) = 0.E0 
          EMI COEF(K,1) = 0.07E0
   78 CONTINUE
      CV AVE = CV AVE * POINTS
  
C........ FIND MEAN RELATIVE HUMIDITY (H AVE), MEAN TEMPERATURE (T AVE) 
C         AND LOCAL INFRARED EMISSIVITY COEFFICIENT (M2/G). 
  
      DO 85 J = 2, LEV
          T AVE(J) = 0.E0 
          H AVE(J) = 0.E0 
          DO 76 K = IL1, IL2
              T AVE(J) = T AVE(J) + T(K,J)
              H AVE(J) = H AVE(J) + H(K,J)
              EMI COEF(K,J) = MERGE(0.07E0,
     1                        MERGE(0.15E0, 2.667E-3*T(K,J) - 0.54338E0,
     2                              T(K,J) .GE. 260.E0),
     3                              T(K,J) .LE. 230.E0) 
   76     CONTINUE
          T AVE(J) = T AVE(J) * POINTS
          H AVE(J) = H AVE(J) * POINTS
   85 CONTINUE
  
C........ COMPUTES LIQUID WATER CONTENT (WC). 
  
      DO 90 J = 2, LEV
          J M 1 = J - 1 
          DO 88 K = IL1, IL2
C . . . . EVALUATE L.W.C. FROM THEORY BY BETTS AND HARSVARDAN (1987). 
C         GAM SAT IS FROM CONVEC ROUTINE = T * (1 - GAM S/GAM D)
              EST(K) = ESTEFF(T(K,J)) 
              HT(K)  = HTVOCP(T(K,J)) 
C . . . . MID LAYER PRESSURE IN [PASCAL]. 
              P = P SFC(K) * SHJ(K,JM1) 
              QST = EPS1*EST(K)/(0.01E0*P-EPS2*EST(K))
              GA = HT(K) * QST / (RGOCP * T(K,J)) 
              GB = GA * EPS1 * HT(K) / T(K,J) 
              GAMSAT = T(K,J) * (GB - GA) / (1.E0 + GB) 
C . . . . GAMMA W = (D THETA/DP) ON THETA ES CONST (BETTS EQ(4))
              GAM W = (PREF/P)**RGOCP * RGOCP/P * GAMSAT
C . . . . AIR DENSITY [KG/M3].
              RHO = P / (RGAS * T(K,J)) 
C . . . . LWC = (CP*1000[G/KG]/L)*(T/THETA)*GAMW*RHO*DP*MIXING. 
C             UNITS [G/M3]. (DP = 3000 PA AND MIXING = 0.56433) 
              CLWC = 1000.E0 / HT(K)
              WC OBS = CLWC * (P/PREF)**RGOCP * GAMW * RHO * 1693.E0
C . . . . FIT OF MEAN WC FROM OBSERVATIONS (FEIGELSON, 1977)
C             WC OBS = (4.35E-5 * T(K,J) - 0.0174) * T(K,J) + 1.734 
C             WC OBS = MAX (WC OBS, 0.01) 
C . . . . FIT OF MINIMUM LIQUID WATER CONTENT (HEYMSFIELD, 1977)
C             WC LOW = 8.64E-18 * EXP (0.1462 * T(K,J)) 
C . . . . SET LIQUID WATER CONTENT IN STRATIFORM CLOUDS.
C             WC STA = MIN (WC OBS, WC LOW) 
              WC STA = WC OBS 
C . . . . SET LIQUID WATER CONTENT IN CONVECTIVE REGIONS. 
C             WC CON = MIN (WC OBS, WC LOW * WC FAC)
C*****        WC CON = WC OBS * WC FAC
C . . . . SELECT STRATIFORM OR CONVECTIVE VALUES OF LWC 
              WC = WC STA 
C . . . . LOCATE UPPER TROPOSPHERE ABOVE CONVECTIVE REGIONS.
C*****        L HIG = (J-1 .GE. NINT(TOP CV(K)))
C*****        IF(L HIG) WC = WC CON 
              WCD(K,J) = WC 
C . . . . LAYER THICKNESS IN METERS.
              DZ = R O G * T(K,J) * LOG(SHTJ(K,J)/SHTJ(K,JM1)) 
C . . . . CORRECTION FOR EXCESS MOISTURE OVER HIGH TOPOGRAPHY (SIGMA CO)
C             FCT LAY = P SFC(K) / P REF
              FCT LAY = 1.E0
C . . . . INTEGRATED LIQUID WATER PATH
              WCL(K,J) = MAX (WC * DZ * FCT LAY, WC MIN)
C . . . . CORRECTION FOR ALBEDO PARADOX USING EMPIRICAL FRACTAL DIMENSION 
C         OF 1.6D (**1.6/2) FROM GABRIEL,LOVEJOY ET AL (1987) 
              WCL(K,J) = WCL(K,J) ** 0.8E0
C. . . . EQUIVALENT LIQUID DROPLETS RADIUS (UM) (FOUQUART, 1985)
C        FOR CIRRUS (T<TICE) THE ICE CRYSTAL LENGTH IS FROM HEYMSFIELD
C        (1977) AND THE EQUIV RADIUS FROM PLATT AND HARSHVARDHAN (1988).
              IF(T(K,J).GE.TICE) THEN 
                  RAD EQV(K,J) = 11.E0 * WCD(K,J) + 4.E0
              ELSE
                  ALWC(K)=LOG10(WCD(K,J))
                  XLENGTH(K) = 1.E-3 * (0.698E0 + 0.366E0 * ALWC(K) + 
     &              0.122E0 * ALWC(K)**2 + 0.0136E0 * ALWC(K)**3)
                  RAD EQV(K,J) = 0.2E0 * 0.0285E6 *
     &              XLENGTH(K) ** 0.786E0 
              ENDIF 
   88     CONTINUE
   90 CONTINUE
  
C........ EVALUATE THE ZONAL MEAN CLOUDINESS (GENERATES CLOUDS).
  
      IF (KLD GEN .AND. .NOT. KLD LOC) THEN 
  
          DO 92 J = 2, LEV
              J M 1 = J - 1 
              J P 1 = J + 1 
  
C . . . . CACULATE THE VARIANCE OF HUMIDITY IN LATITUDE BELT
  
              C MIN 3 = 0.E0
              H VAR   = 0.E0
              DO 91 K = IL1, IL2
                   H VAR = H VAR + (H(K,J) - H AVE(J)) ** 2 
                  IF (H(K,J) .GT. 0.99E0) C MIN 3 = C MIN 3 + H(K,J)
   91         CONTINUE
              C MIN 3 = C MIN 3 * POINTS
              H VAR = H VAR * POINTS
  
C . . . . LINEAR C VS (H,T) FORMULATION 
  
            IF (H AVE(J) .LT. 1.E0) THEN
              C MIN 1 = F * H VAR / ((1.E0 - H AVE(J))**2 + F * H VAR)
            ELSE
              C MIN 1 = 1.0E0 
            END IF
              C MIN 2 = SHELL * H AVE(J) * (H AVE(J) * (1.E0-H AVE(J))) 
              C MIN = MIN (C MIN 1, C MIN 2)
              C MIN = MAX (C MIN, C MIN 3)
  
              ETAO = 1.6E0 / SQRT(COS L)
C             ETA = MAX(EXP (0.0477866 * T AVE(J) - 11.005759), ETAO) 
              ETA = MAX(0.77E0 * COS L * EXP(3.3E0 *
     &          SIG HAF(J-1)), ETAO) 
              C ZO = MAX (H AVE(J) / (1.E0 + ETA * (1.E0 - H AVE(J))),
     1                       C MIN) 
  
C . . . . KEEP GENERATED CLOUDS.
  
              C ZON(J,JP1) = MERGE(C ZO, 0.E0, C ZO .GE. CUT)
              C ZON(JP1,J) = C ZON(J,JP1) 
  
   92     CONTINUE
  
C . . . . DO NOT ALLOW FOR CLOUD IN THE LOWEST MODEL LAYER
C         (DUE TO EXCESSIVE MOISTURE IN THE SOME GCM VERSIONS)
  
          C ZON(LEVP1,LEV) = 0.E0 
          C ZON(LEV,LEVP1) = 0.E0 
  
      END IF
  
C........ LOOPS ON REFERENCE LAYER (SKIP THE MOON LAYER (J=1)). 
  
      DO 120 J = 2, LAY 
            J P 1 = J + 1 
            I     = J P 1 
            I M 1 = J 
  
C........ INITIALIZES CLOUD AND HUMIDITY VECTORS OR GENERATES LOCAL 
C         CLOUDS. (DO NOT ALLOW CLOUDS FOR MID-LAYER LOCATED ABOVE SIGMA
C         = 0.1 FOR "WARM" STRATOSPHERIC TEMPERATURES, AS OBSERVED).
  
         CLD Z = C ZON(J,I) 
         IF (KLD LOC) CLD Z = 1.E0
         IF (SIG HAF(J-1) .LT. 0.1E0.AND.T AVE(J) .GT. 200.E0)
     &     CLD Z = 0.E0 
  
         IF (CLD Z .GT. 0.E0) THEN
  
C . . . . GENERATE LOCAL CLOUDS INSTEAD OF ZONAL & DISTRIBUTED CLOUDS.
C         NOTE: "DC/DH" DETERMINES THE THRESHOLD FOR CLOUD IN GRID BOX. 
  
            IF (KLD LOC) THEN 
              DO 99 K = IL1, IL2
C 
C               * C = (H - H0)/(1-H0), H0 = A *SIG + B
C               * EXCEPT H0=CONSTANT ABOVE MID-MASS POINT ABOVE GROUND. 
C 
                H0(K)=0.90E0-0.125E0*(0.915E0-SHJ(K,J-1)) 
                IF(SHJ(K,J-1) .LE. 0.5E0) H0(K)=0.85E0
C 
                IF( (H(K,IM1) - H0(K)).GT.0.0E0) THEN 
                  CVEC(K) = (H(K,IM1) - H0(K))/(1.E0-H0(K)) 
                ELSE
                  CVEC(K) = 0.E0
                ENDIF 
C 
                IF(J.EQ.LAY) CVEC(K) = 0.0E0
   99         CONTINUE
  
            ELSE
C . . . . PREPARES VECTORS TO DISTRIBUTES ZONAL CLOUDS. 
C        NOTE: TOPO IS A TEMPORARY FIX TO REMOVE THE SYSTEMATIC 
C              BIAS DUE TO SIGMA COORDINATE AND LEADING TO
C              EXCESSIVE RELATIVE HUMIDITY ABOVE HIGH TERRAINS. 
  
              H BAR = 0.E0
              DO 100 K = IL1, IL2 
C               TOPO = CVMGT(1., (PSFC(K)/PREF)**0.5, PSFC(K).GT.90000.)
                TOPO = 1.E0 
                H VEC(K) = H(K,IM1) * TOPO
                C VEC(K) = CLD Z
                H BAR = H BAR + H VEC(K)
  100         CONTINUE
              H BAR = H BAR * POINTS
            END IF
  
C........ DISTRIBUTE CLOUDS IN THE LONGITUDE DIRECTION. 
  
            IF (KLD DIS) THEN 
              METHOD = DIS MET
              IF (METHOD .EQ. 0) THEN 
                     METHOD = 1 
                     IF (J .GE. ILEV) METHOD = 2
C                    FRACTN = SIG HAF(J-1)
              END IF
  
              CALL DISTRIB (C VEC, H VEC, CLD Z, H BAR, ILG, IL1, IL2,
     1                      I TEST, C REF, METHOD, FRACTN)
            END IF
  
C........ STORE CLOUDINESS IN UPPER TRIANGULAR CLOUD MATRIX.
  
            DO 110 K = IL1, IL2 
                  C MTX(K,J,I) = MERGE(CVEC(K),0.E0,CVEC(K).GE.CUT) 
  110       CONTINUE
  
         END IF 
  
  120 CONTINUE
  
C........ INTEGRATE LIQUID WATER PATH 
  
      ZF = 0.E0 
      DO 200 K = IL1, IL2 
          FCT LAY = MERGE(FCT LAY O, ZF, C MTX(K,LEV,LEVP1) .GT. CUT)
          WCL(K,LEV) = WCL(K,LEV) * FCT LAY 
          TAC(K,LEV) = 1.5E0 * WCL(K,LEV) / RAD EQV(K,LEV)
C . . . . PARAMETRIZATION BY PLATT AND HARSVARDHAN (1988) FOR CIRRUS: 
          IF (T(K,LEV) .LT. TICE) THEN
              E(K) = 1.E0 - EXP(-0.75E0 * TAC(K,LEV)) 
          ELSE
              E(K) = 1.E0 - EXP(-EMI COEF(K,LEV) * WCL(K,LEV))
          ENDIF 
          C MTX(K,LEVP1,LEV) = C MTX(K,LEV,LEVP1) * E(K)
  200 CONTINUE
  
      DO 230 J = 1, LEV - 2 
          J P 1 = J + 1 
          J P 2 = J + 2 
          J P 3 = J + 3 
          DO 210 K = IL1, IL2 
              C = C MTX(K,JP1,JP2)
              C TOTAL = MIN (CMTX(K,J,JP1), C, CMTX(K,JP2,JP3)) 
              C SINGL = C - C TOTAL 
              CD = MAX (C, 0.0001E0) 
              FCT LAY = MAX((C TOTAL + FCT LAY O * C SINGL)/ CD, ZF)
              WCL(K,JP1) = WCL(K,JP1) * FCT LAY 
              TAC(K,JP1) = 1.5E0 * WCL(K,JP1) / RAD EQV(K,JP1)
              IF (T(K,JP1) .LT. TICE) THEN
                  E(K) = 1.E0 - EXP(-0.75E0 * TAC(K,JP1)) 
              ELSE
                  E(K) = 1.E0 - EXP(-EMI COEF(K,JP1) * WCL(K,JP1))
              ENDIF 
              C MTX(K,JP2,JP1) = C MTX(K,JP1,JP2) * E(K)
  210     CONTINUE
  230 CONTINUE
C
C........ DETERMINE THE OVERLAP (FULL OR RANDOM OVERLAP).                         
C........ STORE (EMISSIVITY * CLOUDINESS) IN LOWER TRIANGULAR MATRIX.             
C
       DO 260 I = 1, LEV                                                
           I P 1 = I + 1                                                
           I P 2 = I + 2                                                
           DO 235 K = IL1, IL2                                          
               ZRDM(K)  = 0.E0
               C REF(K) = C MTX(K,I,IP1)
               E REF(K) = C MTX(K,IP1,I)
  235      CONTINUE                                                     
           DO 250 J = IP2, LEVP1                                        
              J M 1 = J - 1                                             
              DO 240 K = IL1, IL2
                 IF(CMTX(K,JM1,J) .LE. CUT) ZRDM(K) = 1.E0
                 IF(ZRDM(K).EQ.0.E0)          THEN
                    C MTX(K,I,J) = MAX (C MTX(K,JM1,J), C MTX(K,I,JM1)) 
                    C MTX(K,J,I) = MAX (C MTX(K,J,JM1), C MTX(K,JM1,I))
                 ELSE 
                    C MTX(K,I,J) = MAX ( C MTX(K,I,JM1),
     &                1.E0 - (1.E0 - C MTX(K,JM1,J)) *
     &                (1.E0 - C REF(K)) ) 
                    C MTX(K,J,I) = MAX ( C MTX(K,JM1,I),              
     &                1.E0 - (1.E0 - C MTX(K,J,JM1)) *
     &                (1.E0 - E REF(K)) )
                    IF(CMTX(K,JM1,J) .LE. CUT)            THEN
                       C REF(K) = C MTX(K,I,J)                          
                       E REF(K) = C MTX(K,J,I)
                    ENDIF
                 ENDIF
  240         CONTINUE                                                  
  250     CONTINUE                                                      
  260 CONTINUE
C 
C........ LIMIT RESULTING CLOUD MATRIX VALUES TO AVOID NUMERICAL PROBLEMS.
C 
      DO 320 J = 1, LEVP1 
      DO 320 I = 1, LEVP1 
          DO 310 K = IL1, IL2 
              C MTX(K,J,I) = MIN(C MTX(K,J,I), 0.999E0) 
  310     CONTINUE
  320 CONTINUE
C 
C     * SCALE LIQUID WATER BY CLOUD AMOUNT. 
C 
      DO 340 J=2,LEV
          DO 330 K=IL1,IL2
              WCD(K,J) = WCD(K,J)*C MTX(K,J,J+1)
  330     CONTINUE
  340 CONTINUE
C 
C///////////////////////////////////////////////////////////////////////
      RETURN
  
C-----------------------------------------------------------------------
  
      END 
