      SUBROUTINE CLOUDS7(C MTX, TAC, SHJ, SHTJ, TF, T, H,               
     1                   CG, P SFC, PBLT, TCV, GC,
     2                   CLSTFAC, RADEQV, DZ, WCL, WCD,                 
     3                   WCDW, WCDI, WCLW, WCLI, RADEQVW,               
     4                   RADEQVI, FRACW, DUMMY, LTOP, CCLD,             
     5                   CVEC, CR1, CR2, HTHRESH, CRH,                  
     6                   GAMSAT, H0LIQ, TMOIST, RATIO, GAMMA,           
     7                   H0, ZRDM, CREF, EREF,                          
     8                   ILG, IL1, IL2, LEVP1, LEV, ILEV,               
     9                   MLEVP1, MLEV, MILEV)
C
C     * SEP 06/95 - H.BARKER, SLIGHT MODIFICATION OF BACKSCATTERING 
C     *             M.LAZARE. ASSYMMETRY COEFFICIENTS.
C
C     * AUG 16/95 - M.LAZARE. PREVIOUS VERSION CLOUDS6 FOR GCM10.
C
C     * JAN 12/94 - M.LAZARE. PREVIOUS VERSION CLOUDS5 FOR GCM8/GCM9.
C                                                                                 
C     ***************** NOTE: THIS SUBROUTINE IS CALLED BY THE CLOUD ***          
C     *                       DIAGNOSTIC PROGRAM CLDPRG4 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), CG(ILG,2,MLEV)
                                                                        
C     * INPUT ARRAYS...                                                           
                                                                        
      DIMENSION SHJ(ILG,MILEV), SHTJ(ILG,MLEV), TF(ILG,MLEV)            
      DIMENSION T(ILG,MLEV), H(ILG,MLEV)                                
      DIMENSION P SFC(ILG), PBLT(ILG), TCV(ILG), GC(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 FRACW(ILG,MLEV), DUMMY(ILG,MLEV)                        
                                                                        
      DIMENSION CVEC(ILG), CR1(ILG), CR2(ILG), HTHRESH(ILG), CRH(ILG)   
      DIMENSION GAMSAT(ILG), H0LIQ(ILG), TMOIST(ILG), RATIO(ILG)        
      DIMENSION H0(ILG), CCLD(ILG), Z RDM(ILG), CREF(ILG), EREF(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     * STATEMENT FUNCTION DEFINITIONS.                                           
                                                                        
C     * A) 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     * 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,UUU) = UUU*TW(TTT) + (1.E0-UUU)*TI(TTT)                
                                                                        
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,UUU) = UUU*ESW(TTT) + (1.E0-UUU)*ESI(TTT)              
                                                                        
C     * D) COMPUTES THE EFFECTIVE WATER/ICE CLOUD PROPERTY.                       
                                                                        
      WIEFF(TTT,UUU,VVV) = TTT*UUU + (1.E0-TTT)*VVV                     
                                                                        
C........ PHYSICAL CONSTANTS:                                                     
      DATA P REF /101325.E0/                                            
                                                                        
C........ SCHEME PARAMETERS:                                                      
      DATA WCMIN /0.00001E0/, CUT /0.001E0/, H0ICE /0.70E0/             
      DATA DZW /100.E0/, DZI /100.E0/, FCTLAYO /0.5E0/
      DATA RHOW /1000.E0/, RHOI /780.E0/      
      DATA CDDL /500.E0/, CDDW /100.E0/, THIRD /0.3333333333333E0/
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               
C
      DO 80 K = IL1, IL2                                                
            TAC(K,1) = 0.E0                                             
   80 CONTINUE                                                          
C                             
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                           
      DO 90 J = 2, LEV                                                  
        DO 85 K = IL1, IL2                                              
C                                                                                 
C         * COMPUTE THE FRACTIONAL PROBABILITY OF WATER PHASE                     
C         *  EXISTING AS A FUNCTION OF TEMPERATURE (FROM ROCKEL,                  
C      '  *  RASCHKE AND WEYRES, BEITR. PHYS. ATMOSPH., 1991.)                    
C                                                                                 
          FRACW(K,J) = MERGE( 1.E0,                                     
     &      0.0059E0+0.9941E0*EXP(-0.003102E0*(T1S-T(K,J))**2), 
     &      T(K,J).GE.T1S )                               
                                                                        
          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),FRACW(K,J))                               
          HT  = HTVOCP(T(K,J),FRACW(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 * DZW                                     
          DPMIXI = RHO * GRAV * 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) = WIEFF(FRACW(K,J),WCDW(K,J),WCDI(K,J))              
C                                                                                 
C         * GENERATE LOCAL CLOUDS.                                                
C                                                                                 
          H0(K)=WIEFF(FRACW(K,J),HM,H0ICE)                              
          CCLD(K)=(MAX(H(K,J)-H0(K),0.E0))/(1.E0-H0(K))                 
C                                                                                 
C         * CALCULATE EQUIVALENT RADIUS (MICRONS) FOR BOTH PHASES.
C         * WATER PHASE PORTION ASSUMES GAMMA SIZE DISTRIBUTION
C         * PARAMETER=3 AND ESTIMATES FOR CLOUD DROP NUMBER DENSITY
C         * IN (CM**-3) FOR LAND AND OCEAN/ICE GIVEN BY CDDL AND
C         * CDDW, RESPECTIVELY.
C
          IF(GC(K).GE.-0.5E0)     THEN
            CDD=CDDW
          ELSE
            CDD=CDDL
          ENDIF                                                         
          RADEQVW(K,J) = 75.46E0*(WCDW(K,J)/CDD)**THIRD
          RADEQVI(K,J) = 10.E0
C
C         * DEFINE CLOUD ASSYMETRY FACTOR FOR BOTH SOLAR SPECTRAL
C         * INTERVALS. WATER PHASE FOLLOWS SLINGO (1989) FOR PURE WATER
C         * SPHERES BUT SLIGHTLY REDUCED. VALUE OF ICE CLOUDS IS A
C         * SPECIFIED CONSTANT AND USUAL ICE/WATER MIX ASSUMED.
C
          CGW1 = 0.829E0 + 0.0025E0*RADEQVW(K,J)
          CGW2 = 0.782E0 + 0.0050E0*RADEQVW(K,J)
          CGI  = 0.750E0
          CG(K,1,J) = WIEFF(FRACW(K,J),CGW1,CGI)
          CG(K,2,J) = WIEFF(FRACW(K,J),CGW2,CGI)
C                                                                                 
C         * STORE CLOUDINESS IN UPPER TRIANGULAR CLOUD MATRIX.                    
C                                                                                 
          C MTX(K,J,J+1) = MERGE(CCLD(K),0.E0,CCLD(K).GE.CUT)           
                                                                        
   85   CONTINUE                                                        
   90 CONTINUE                                                          
C                                                                                 
C     * CALCULATE RESULTING PATH, CLOUD EMMISSIVITY AND OPTICAL DEPTH.            
C                                                                                 
C     * THEN RE-CALCULATE THE EFFECTIVE LIQUID WATER PATH AND EQUIVILENT          
C     * RADIUS FOR DIAGNOSTIC PURPOSES.                                           
C                   
      DO 230 J = 2, LEV                                                 
         J M 1 = J - 1                                                  
         J P 1 = J + 1                                                  
         J P 2 = J + 2                                                  
         DO 200 K = IL1, IL2                                            
            C = C MTX(K,J,JP1)                                          
            IF(J.EQ.LEV)                                        THEN    
              FCTLAY = FCTLAYO                                          
            ELSE                                                        
              IF(CMTX(K,JM1,J).EQ.0.E0.AND.CMTX(K,JP1,JP2).EQ.0.E0) THEN
                FCTLAY=FCTLAYO                                          
              ELSE IF(CMTX(K,JM1,J).NE.0.E0.AND.CMTX(K,JP1,JP2).NE.0.E0)
     &            THEN
                FCTLAY=1.E0                                             
              ELSE                                                      
                FCTLAY=1.5E0*FCTLAYO                                    
              ENDIF                                                     
            ENDIF                                                       
            WCLW(K,J) = WCDW(K,J) * DZ(K,J) * FCTLAY
            WCLI(K,J) = WCDI(K,J) * DZ(K,J) * FCTLAY
            TACW = 1.5E0*WCLW(K,J)/(RADEQVW(K,J)*0.001E0*RHOW)
            TACI = 1.5E0*WCLI(K,J)/(RADEQVI(K,J)*0.001E0*RHOI) 
            TAC(K,J) = WIEFF(FRACW(K,J),TACW,TACI)
            IF(CMTX(K,J,JP1).LE.CUT)   TAC(K,J) = 0.E0                  
            RAD EQV(K,J) = WIEFF(FRACW(K,J),RADEQVW(K,J),RADEQVI(K,J))  
            EW = 1.E0 - EXP(-0.75E0 * TACW)           
            EI = TACI*(1.E0-0.5E0*TACI*LOG((2.E0+TACI)/TACI))
            E = WIEFF(FRACW(K,J),EW,EI)
            C MTX(K,JP1,J) = C MTX(K,J,JP1) * E                         
            WCL(K,J) = WIEFF(FRACW(K,J),WCLW(K,J),0.E0)                 
  200    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-----------------------------------------------------------------------          
      RETURN
      END
