      SUBROUTINE CLOUDS4(C MTX, TAC, SHJ, SHTJ, W, TF,
     1                   T, H, P SFC, PBLT, SHLT, TCV,
     2                   CLSTFAC, RADEQV, DZ, WCL, WCD, 
     3                   WCDW, WCDI, WCLW, WCLI, RADEQVW, 
     4                   RADEQVI, TACW, TACI, LTOP, CCLD, 
     5                   CVEC, CR1, CR2, HTHRESH, CRH,
     6                   GAMSAT, GAMFAC, TMOIST, RATIO, GAMMA,
     7                   H0, IRDM, CREF, EREF,
     8                   ILG, IL1, IL2, LEVP1, LEV, ILEV, 
     9                   MLEVP1, MLEV, MILEV) 
  
C     * MAY 09/91 - M.LAZARE. FINAL, FROZEN VERSION FOR GCM7. LIKE
C     *                       PREVIOUS INTERMEDIATE EXCEPT: 
C     *                       A. WATER/ICE CONTENT PARAMETERIZED IN 
C     *                          TERMS OF DEPLETION FACTORS (DFW,DFI) 
C     *                          IN TERMS OF HEIGHT, WITH UPPER BOUNDS
C     *                          (DZW,DZI). 
C     *                       B. FRACTAL SCALING REMOVED. 
C     *                       C. SCALING OF ICE CONTENT IN "CONVECTIVE
C     *                          REGIONS" REMOVED ("LTOP" NOW SET TO
C     *                          PBLT). 
C     *                       D. WCMIN NOW APPLIED (CORRECTLY) TO WATER 
C     *                          CONTENT AND NOT WATER PATH.
C     *                       E. THRESHOLD RELATIVE HUMIDTY NOW VARIABLE
C     *                          (HTHRESH) INSTEAD OF HM. IN PBLT,
C     *                          HTHRESH=H0=0.99. IN SHALLOW MIXING 
C     *                          REGION (DEFINED BY NEW PASSED "SHLT" 
C     *                          ARRAY), HTHRESH=H0=0.95. OTHERWISE, H0 
C     *                          VARIES BETWEEN H0HI(0.90) AND H0LOW
C     *                          (VALUE CHANGED TO 0.70) BASED ON STABILITY.
C     *                       F. "TFROW" NOW PASSED IN AND USED TO PROPERLY 
C     *                          DEFINE CLSTFAC.
  
C     * AUG 15/90 - M.LAZARE. INTERMEDIATE VERSION FOR GCM7. LIKE 
C     *                       PREVIOUS INITIAL EXCEPT:  
C     *                       A. H0=HM (PASSED) =0.95 AND H0LOW=0.75. 
C     *                       B. SLIGHT MODS TO STABILITY FACTOR. 
C     *                       C. WCD LIMITED TO VALUES GE. 0.01.
C     *                       D. ICE PHASE INTERPOLATION IN STATEMENT 
C     *                          FUNCTIONS NOW BETWEEN T1S->T2S, AS IN
C     *                          CONVEC7,SPWCON7,RELHUM3. 
C     *                       E. SIMILAR STATEMENT FUNCTIONS ADDED FOR
C     *                          CLOUD WATER CONTENT, CLOUD OPTICAL 
C     *                          DEPTH, CLOUD EMISSIVITY AND CLOUD WATER
C     *                          PATH.
C     *                       F. FRACTAL DIMENSION OF 0.7 USED ON TAC 
C     *                          INSTEAD OF WCL, ONLY IN CONVECTIVE 
C     *                          REGIONS (DEFINED BY LTOP). 
C     *                       G. FULLY VECTORIZED WITH NECESSARY WORK 
C     *                          ARRAYS ADDED.
C 
C     * AUG 28/89 - M.LAZARE. INITIAL VERSION FOR GCM7. LIKE
C     *                       PREVIOUS CLOUDS3 EXCEPT:  
C     *                       A. MAJOR CLEANUP TO REMOVE CODE CONCERNED 
C     *                          WITH DISTRIBUTING ZONALLY-AVERAGED 
C     *                          CLOUDS.
C     *                       B. SLIGHT MODIFICATIONS TO H0(K)
C     *                          CALCULATION. 
C     *                       C. NAM'S STABILITY MODIFICATION TO CLOUD
C     *                          AMOUNT INCLUDED. 
C     *                       D. PASSED MODEL OPTION PARAMETER "ICON" 
C     *                          USED TO DECIDE WHETHER TO PUT CLOUDS 
C     *                          IN LOWEST LAYER OR NOT.
C 
C     ***************** NOTE: THIS SUBROUTINE IS CALLED BY THE CLOUD ***
C     *                       DIAGNOSTIC PROGRAM CLDPRG2 AND ANY     ***
C     *                       CHANGES IN EITHER MUST BE MUTUTALLY    ***
C     *                       CONSISTENT!!!!!!!!!!!!! 
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), TF(ILG,MLEV)
      DIMENSION W(ILG,MILEV), T(ILG,MLEV), H(ILG,MLEV)
      DIMENSION P SFC(ILG), PBLT(ILG), SHLT(ILG), TCV(ILG)
  
C     * WORK ARRAYS...
  
      DIMENSION CLSTFAC(ILG,MLEV), RAD EQV(ILG,MLEV), DZ(ILG,MLEV)
      DIMENSION WCL(ILG,MLEV), WCD(ILG,MLEV)
      DIMENSION WCDW(ILG,MLEV), WCDI(ILG,MLEV)
      DIMENSION WCLW(ILG,MLEV), WCLI(ILG,MLEV)
      DIMENSION RADEQVW(ILG,MLEV), RADEQVI(ILG,MLEV)
      DIMENSION TACW(ILG,MLEV), TACI(ILG,MLEV)
  
      DIMENSION C VEC(ILG), CR1(ILG), CR2(ILG), HTHRESH(ILG), CRH(ILG)
      DIMENSION GAMSAT(ILG), GAMFAC(ILG), TMOIST(ILG), RATIO(ILG) 
      DIMENSION H0(ILG), CCLD(ILG), I RDM(ILG), C REF(ILG), E REF(ILG)
      DIMENSION GAMMA(ILG)
  
      INTEGER LTOP(ILG) 
  
      COMMON /ADJPCP/ HC,HF,HM,AA,DEPTH,LHEAT,MOIADJ,MOIFLX 
      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,EPS1,EPS2 
C 
C-----------------------------------------------------------------------
C     * STATEMENT FUNCTION DEFINITIONS. 
C     * WATER->ICE ONSET ASSUMED FROM 0->-20 DEGREES CELCIUS, IN LINE 
C     * WITH SEVERAL INDEPENDANT OBSERVATIONS.
  
C     * A) COMPUTES THE EFFECTIVE CLOUD EMISSIVITY. 
  
      EMMW(UUU)                 = 1.E0 - EXP(-0.14E0 * UUU) 
      EMMI(VVV)                 = 1.E0 - EXP(-0.75E0 * VVV) 
      EMI(TTT,UUU,VVV)          = MERGE( EMMW(UUU), 
     1                            MERGE( EMMI(VVV), 
     2         SLP*((TTT-T2S)*EMMW(UUU)+(T1S-TTT)*EMMI(VVV)), 
     3                      TTT.LE.T2S),
     4                      TTT.GE.T1S) 
  
C     * B) COMPUTES THE RATIO OF LATENT HEAT OF VAPORIZATION OF 
C     *    WATER OR ICE TO THE SPECIFIC HEAT OF AIR AT CONSTANT 
C     *    PRESSURE CP. 
  
      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) 
  
C     * C) COMPUTES THE SATURATION VAPOUR PRESSURE OVER WATER OR ICE. 
  
      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     * D) COMPUTES THE EFFECTIVE CLOUD WATER CONTENT.
  
      WCDEFF(TTT,UUU,VVV)         = MERGE( UUU, 
     1                              MERGE( VVV, 
     2                   SLP*((TTT-T2S)*UUU+(T1S-TTT)*VVV), 
     3                   TTT.LE.T2S), 
     4                   TTT.GE.T1S)
  
C     * E) COMPUTES THE EFFECTIVE CLOUD OPTICAL DEPTH.
  
      TACEFF(TTT,UUU,VVV)         = MERGE( UUU, 
     1                              MERGE( VVV, 
     2          SLP*((TTT-T2S)*UUU+(T1S-TTT)*VVV),
     3          TTT.LE.T2S),
     4          TTT.GE.T1S) 
  
C     * F) COMPUTES THE EFFECTIVE CLOUD WATER PATH. 
  
      WCLEFF(TTT,UUU,VVV)         = MERGE( UUU, 
     1                              MERGE( VVV, 
     2          SLP*((TTT-T2S)*UUU+(T1S-TTT)*VVV),
     3          TTT.LE.T2S),
     4          TTT.GE.T1S) 
  
C     * G) COMPUTES THE CRITICAL SATURATION RELATIVE HUMIDITY.
  
      CR(HM,HHH,AA)     =     MERGE( AA*(HHH-HM)**3,
     &  AA*(2.E0-HM-HHH)**3+HHH-1.E0,
     &  HHH.LE.1.E0) 
      CRIRLH(HM,HHH,AA) = HHH-MERGE( CR(HM,HHH,AA), HHH-1.E0,
     &  HM+HHH.LE.2)
  
C........ PHYSICAL CONSTANTS: 
      DATA P REF /101325.E0/
  
C........ SCHEME PARAMETERS:  
      DATA WC MIN /0.0001E0/, CUT /0.01E0/, H0LOW /0.70E0/, H0HI/0.90E0/
      DATA DFW /1.00E0/, DFI /1.00E0/, DFWUNS /0.03E0/, DFIUNS/0.03E0/
C      * WHEN ATMOSPHERIC Q FIELD IMPROVED (INCREASED), USE THE 
C      * FOLLOWING INSTEAD: 
C 
C     DATA DZW /300./, DZI /150./ 
      DATA DZW /300.E0/, DZI /300.E0/ 
C-----------------------------------------------------------------------
  
  
C........ INITIALIZE ARRAYS 
  
      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
   60       CONTINUE
   70 CONTINUE
  
      DO 80 K = IL1, IL2
            TAC(K,1) = 0.E0 
            LTOP(K)  = NINT(PBLT(K))
   80 CONTINUE
  
C     * COMPUTES LIQUID WATER CONTENT (WCD).
C     * EVALUATE L.W.C. FROM THEORY BY BETTS AND HARSVARDAN (1987). 
C     * LWC = (CP*1000[G/KG]/L)*(T/THETA)*GAMW*RHO*DPMIX (UNITS [G/M3]).
C     * GAMMA W = (D THETA/DP) ON THETA ES CONST (BETTS EQ(4)). 
C     * GAM SAT IS FROM CONVEC ROUTINE = T * (1 - GAM S/GAM D). 
C 
C     * WATER CONTENT CALCULATION IS BASED ON "DEPLETION" FACTORS FOR 
C     * WATER AND ICE PHASES, EACH HAVING STABLE (CLSTFAC APPROACHING 
C     * INFINITY) AND UNSTABLE (CLSTFAC APPROACHING ZERO) LIMITS. THE 
C     * UNSTABLE LIMITS ARE LOWER DUE TO ENTRAINMENT OF DRIER AIR IN
C     * THE MIXING PROCESS. 
C 
C     * CLOUDS ARE "EXPANDED" IN THE HORIZONTAL UNDER STABLE CONDITIONS 
C     * BUT REDUCED IN VERTICAL EXTENT ACCORDINGLY, IN ORDER TO CONSERVE
C     * MASS. THIS IS ACCOMPLISHED USING THE CLOUD STABILITY PARAMETER
C     * "CLSTFAC".
  
      DO 90 J = LEV, 2, -1
        DO 88 K = IL1, IL2
          RATIO(K)=LOG(SHTJ(K,J)/SHTJ(K,J-1))
          GAMMA(K)=T(K,J)*(1.E0-LOG(TF(K,J)/TF(K,J-1))/
     &      (RATIO(K)*RGOCP))
          GAMMA(K)=MAX(GAMMA(K),0.E0) 
          EST = ESTEFF(T(K,J))
          HT  = HTVOCP(T(K,J))
          P = P SFC(K) * SHJ(K,J-1) 
          QST = EPS1*EST/(0.01E0*P-EPS2*EST)
          GA = HT * QST / (RGOCP * T(K,J))
          GB = GA * EPS1 * HT / T(K,J)
          GAMSAT(K) = T(K,J) * (GB - GA) / (1.E0 + GB)
          GAM W = (PREF/P)**RGOCP * RGOCP/P * GAMSAT(K) 
          RHO = P / (RGAS * T(K,J)) 
          CLWC = 1000.E0 / HT 
          CLSTFAC(K,J)=(MAX(GAMMA(K)-GAMSAT(K),0.E0)/GAMSAT(K))**2
          DZ(K,J) = (RGAS/GRAV) * T(K,J) * LOG(SHTJ(K,J)/SHTJ(K,J-1))
          DPMIXW = RHO * GRAV * MIN(
     &      DZ(K,J)*(DFWUNS+DFW*CLSTFAC(K,J))/(1.E0+CLSTFAC(K,J)), 
     &      DZW)
          DPMIXI = RHO * GRAV * MIN(
     &      DZ(K,J)*(DFIUNS+DFI*CLSTFAC(K,J))/(1.E0+CLSTFAC(K,J)), 
     &      DZI)
          WCDW(K,J) = CLWC * (P/PREF)**RGOCP * GAMW * RHO * DPMIXW
          WCDI(K,J) = CLWC * (P/PREF)**RGOCP * GAMW * RHO * DPMIXI
          WCDW(K,J)=MAX(WCDW(K,J),WCMIN)
          WCDI(K,J)=MAX(WCDI(K,J),WCMIN)
          WCD(K,J) = WCDEFF(T(K,J),WCDW(K,J),WCDI(K,J)) 
   88   CONTINUE
   90 CONTINUE
C 
C     * CLOUD AMOUNT LOOP.
C 
      DO 120 J = 2, LEV 
        DO 100 K = IL1, IL2 
C 
C          * GENERATE LOCAL CLOUDS. 
C          * ASSUME THAT THE PLANETARY BOUNDARY LAYER (THE TOP LEVEL
C          * OF WHICH IS DEFINED IN VRTDFS_) HAS ENOUGH TURBULENCE
C          * TO NEGATE THE ONSET RELATIVE HUMIDITY DEPENDENCE ON
C          * STABILITY. 
C 
           IF(J-1.LE.NINT(PBLT(K)))                           THEN
              HTHRESH(K)=HM 
           ELSE 
              HTHRESH(K)=0.99E0 
           ENDIF
           IF(H(K,J).GT.HTHRESH(K))                          THEN 
              CRH(K)   = CRIRLH(HTHRESH(K),H(K,J),AA) 
           ELSE 
              CRH(K)   = MIN(H(K,J),1.E0) 
           ENDIF
C 
           IF(J-1.LE.NINT(SHLT(K)))                                 THEN
             H0(K)=(H0HI+H0LOW*CLSTFAC(K,J))/(1.E0+CLSTFAC(K,J))
           ELSE IF(J-1.GT.NINT(SHLT(K)) .AND. J-1.LE.NINT(PBLT(K))) THEN
             H0(K)=HM 
           ELSE 
             H0(K)=0.99E0 
           ENDIF
           CCLD(K)=(MAX(CRH(K)-H0(K),0.E0))/(1.E0-H0(K))
C 
C          * EXPAND AREAL CLOUD COVERAGE UNDER STABLE CONDITIONS, BUT 
C          * REDUCE VERTICAL EXTENT ACCORDINGLY, IN ORDER TO CONSERVE 
C          * MASS.
C 
           CVEC(K)=CCLD(K)*(1.E0+CLSTFAC(K,J))/
     &       (1.E0+CCLD(K)*CLSTFAC(K,J))
           IF(CVEC(K).GT.CUT) DZ(K,J)=DZ(K,J)*CCLD(K)/CVEC(K) 
C 
C          * CALCULATE LIQUID WATER PATH (G/M2) AND EQUIVALENT RADIUS 
C          * (MICRONS), AS WELL AS OPTICAL DEPTH, FOR BOTH PHASES.
C 
           WCLW(K,J) = WCDW(K,J) * DZ(K,J)
           WCLI(K,J) = WCDI(K,J) * DZ(K,J)
C 
           RADEQVW(K,J) = 11.E0*WCDW(K,J) + 4.E0
           RADEQVI(K,J) = MAX(MIN(74.6E0+(T(K,J)-273.16E0),
     &       74.6E0) , 5.E0)
C 
           TACW(K,J) = 1.5E0*WCLW(K,J)/RADEQVW(K,J) 
           TACI(K,J) = 1.5E0*WCLI(K,J)/RADEQVI(K,J) 
C 
C          * STORE CLOUDINESS IN UPPER TRIANGULAR CLOUD MATRIX. 
C 
           C MTX(K,J,J+1) = MERGE(CVEC(K),0.E0,CVEC(K).GE.CUT)
  
  100   CONTINUE
  120 CONTINUE
  
C........ DETERMINE THE OVERLAP (FULL OR RANDOM OVERLAP). 
  
      DO 190 I = 1, LEV 
          I P 1 = I + 1 
          I P 2 = I + 2 
  
          DO 150 K = IL1, IL2 
              I RDM(K) = I P 1
  150     CONTINUE
  
          DO 180 J = IP2, LEVP1 
              J M 1 = J - 1 
  
C . . . . INDIRECT ADDRESSING FIRST.
  
              DO 160 K = IL1, IL2 
                  C REF(K) = C MTX(K,I,IRDM(K)) 
  160         CONTINUE
  
              DO 170 K = IL1, IL2 
                  C MAX = MAX (C MTX(K,JM1,J), C MTX(K,I,JM1))
                  C RDM = 1.E0 - (1.E0 - C MTX(K,JM1,J)) *
     &              (1.E0 - C REF(K)) 
                  C RDM = MAX (C RDM, C MTX(K,I,JM1)) 
                  IF(CMTX(K,JM1,J) .LE. CUT) I RDM(K) = J
                  C MTX(K,I,J) = MERGE(C MAX, C RDM, I RDM(K).EQ.IP1) 
  170         CONTINUE
  180     CONTINUE
  190 CONTINUE
C 
C     * CALCULATE RESULTING CLOUD EMISSIVITY AND OPTICAL DEPTH. 
C 
C     * THEN RE-CALCULATE THE EFFECTIVE LIQUID WATER PATH AND EQUIVILENT
C     * RADIUS FOR DIAGNOSTIC PURPOSES. 
C     * FINAL, RESULTING TOTAL LIQUID WATER PATH ONLY FOR LIQUID PHASE. 
C 
      DO 230 J = 2, LEV 
         J P 1 = J + 1
         DO 200 K = IL1, IL2
            E = EMI(T(K,J),WCLW(K,J),TACI(K,J)) 
            C MTX(K,JP1,J) = C MTX(K,J,JP1) * E 
            TAC(K,J) = TACEFF(T(K,J),TACW(K,J),TACI(K,J)) 
            WCL(K,J) = WCLEFF(T(K,J),WCLW(K,J),WCLI(K,J)) 
            RAD EQV(K,J) = 1.5E0 * WCL(K,J) / TAC(K,J)
            IF(CMTX(K,J,JP1).LE.CUT)   TAC(K,J) = 0.E0
            WCL(K,J) = WCLEFF(T(K,J),WCLW(K,J),0.E0)
  200    CONTINUE 
  230 CONTINUE
  
C........ STORE (EMISSIVITY * CLOUDINESS) IN LOWER TRIANGULAR MATRIX. 
  
       DO 260 I = 1, LEV
           I P 1 = I + 1
           I P 2 = I + 2
           DO 235 K = IL1, IL2
               I RDM(K) = I P 1 
  235      CONTINUE 
           DO 250 J = IP2, LEVP1
              J M 1 = J - 1 
              DO 238 K = IL1, IL2 
                  E REF(K) = C MTX(K,IRDM(K),I) 
  238         CONTINUE
              DO 240 K = IL1, IL2 
                  E MAX = MAX (C MTX(K,J,JM1), C MTX(K,JM1,I))
                  E RDM = 1.E0 - (1.E0 - C MTX(K,J,JM1)) *
     &              (1.E0 - E REF(K)) 
                  E RDM = MAX (E RDM, C MTX(K,JM1,I)) 
                  IF(CMTX(K,JM1,J) .LE. CUT) I RDM(K) = J
                  C MTX(K,J,I) = MERGE(E MAX, E RDM, I RDM(K).EQ.IP1) 
  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 
      RETURN
  
C-----------------------------------------------------------------------
  
      END 
