      SUBROUTINE CLOUD6O(CMTX, TAC, 
     1                   SHJ, SHTJ, T, H, PSFC, 
     2                   EMICOEF, RADEQV, WCL, WCD, CVEC,
     3                   HVEC, ZRDM, CREF, EREF, EST, 
     4                   HT, ALWC, XLENGTH, H0, E, 
     5                   ILG, IL1, IL2, LEVP1, LEV, ILEV, 
     6                   MLEVP1, MLEV, MILEV) 

C     * JAN 23/95 - M.LAZARE. OPTIMIZED VERSION OF CLOUDS3 WITH OLD
C     *                       FIXED/DISTRIBUTED CLOUDS REMOVED.
C     * NOV 25/94 - M.LAZARE. PREVIOUS VERSION CLOUDS3 FOR MODEL VERSION
C     *                       GCM6.
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     ******************************************************************
  
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

C     * OUTPUT ARRAYS...
  
      DIMENSION CMTX(ILG,MLEVP1,MLEVP1), TAC(ILG,MLEV) 
  
C     * INPUT ARRAYS... 
  
      DIMENSION SHJ(ILG,MILEV), SHTJ(ILG,MLEV)
      DIMENSION T(ILG,MLEV), H(ILG,MLEV) 
      DIMENSION PSFC(ILG)
  
C     * WORK ARRAYS...
  
      DIMENSION EMICOEF(ILG,MLEV), RADEQV(ILG,MLEV) 
      DIMENSION WCL(ILG,MLEV), WCD(ILG,MLEV)
      DIMENSION CVEC(ILG), HVEC(ILG)
      DIMENSION ZRDM(ILG), CREF(ILG), EREF(ILG)
      DIMENSION EST(ILG), HT(ILG), ALWC(ILG), XLENGTH(ILG)
      DIMENSION H0(ILG), E(ILG) 
  
      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. 
  
      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 ROG /29.29/, PREF /101325./, CPOR /3.4983/
  
C........ SCHEME PARAMETERS:  
      DATA WCMIN /0.0001/, CUT /0.01/, FCTLAYO /0.5/
  
C........ PARAMETERS AFFECTING THE RESULTS: 
C 
C             FCT LAY O : VERTICAL FRACTION OF THE LAYER OCCUPIED BY CLD
C                         FOR UNRESOLVED LAYER ONLY 
C-----------------------------------------------------------------------
C........ INITIALIZE ARRAYS 
  
      DO 70 I = 1, LEVP1
            DO 60 J = 1, LEVP1
                  DO 50 K = IL1, IL2
                        C MTX(K,J,I) = 0. 
   50             CONTINUE
   60       CONTINUE
   70 CONTINUE
  
      DO 72 K = IL1, IL2
          WCL(K,1) = 0. 
          WCD(K,1) = 0. 
          TAC(K,1) = 0. 
          EMI COEF(K,1) = 0.07
 72    CONTINUE
  
C........ FIND LOCAL INFRARED EMISSIVITY COEFFICIENT (M2/G). 
  
      DO 85 J = 2, LEV
          DO 76 K = IL1, IL2
              EMI COEF(K,J) = MERGE(0.07,
     1                        MERGE(0.15, 2.667E-3*T(K,J) - 0.54338,
     2                              T(K,J) .GE. 260.),
     3                              T(K,J) .LE. 230.) 
   76     CONTINUE
   85 CONTINUE
  
C........ COMPUTES LIQUID WATER CONTENT (WC). 
  
      DO 90 J = 2, LEV
          J M 1 = J - 1 
C 
          DO 88 K = IL1, IL2
C
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)
C
              EST(K) = ESTEFF(T(K,J)) 
              HT(K)  = HTVOCP(T(K,J)) 
C
C . . . . MID LAYER PRESSURE IN [PASCAL]. 
C
              P = P SFC(K) * SHJ(K,JM1) 
              QST = EPS1*EST(K)/(0.01*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. + GB) 
C
C . . . . GAMMA W = (D THETA/DP) ON THETA ES CONST (BETTS EQ(4))
C
              GAM W = (PREF/P)**RGOCP * RGOCP/P * GAMSAT
C
C . . . . AIR DENSITY [KG/M3].
C
              RHO = P / (RGAS * T(K,J)) 
C
C . . . . LWC = (CP*1000[G/KG]/L)*(T/THETA)*GAMW*RHO*DP*MIXING. 
C         UNITS [G/M3]. (DP = 3000 PA AND MIXING = 0.56433) 
C
              CLWC = 1000. / HT(K)
              WCD(K,J) = CLWC * (P/PREF)**RGOCP * GAMW * RHO * 1693.
C
C . . . . LAYER THICKNESS IN METERS.
C
              DZ = ROG * T(K,J) * LOG(SHTJ(K,J)/SHTJ(K,JM1)) 
C
C . . . . INTEGRATED LIQUID WATER PATH
C
              WCL(K,J) = MAX (WCD(K,J)*DZ, WCMIN)
C
C . . . . CORRECTION FOR ALBEDO PARADOX USING EMPIRICAL FRACTAL DIMENSION 
C         OF 1.6D (**1.6/2) FROM GABRIEL,LOVEJOY ET AL (1987) 
C
              WCL(K,J) = WCL(K,J) ** 0.8
C
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).
C
              IF(T(K,J).GE.TICE) THEN 
                  RAD EQV(K,J) = 11. * WCD(K,J) + 4.
              ELSE
                  ALWC(K)=LOG10(WCD(K,J))
                  XLENGTH(K) = 1.E-3 * (0.698 + 0.366 * ALWC(K) + 
     1                         0.122 * ALWC(K)**2 + 0.0136 * ALWC(K)**3)
                  RAD EQV(K,J) = 0.2 * 0.0285E6 * XLENGTH(K) ** 0.786 
              ENDIF 
C
C . . . . GENERATE LOCAL CLOUDS.
C         C = (H - H0)/(1-H0), H0 = A *SIG + B
C         EXCEPT H0=CONSTANT ABOVE MID-MASS POINT ABOVE GROUND. 
C         DO NOT ALLOW CLOUDS FOR MID-LAYER LOCATED ABOVE SIGMA=0.1
C         FOR "WARM" STRATOSPHERIC TEMPERATURES, AS OBSERVED.
C  
          IF (SHJ(K,J-1).GE.0.1 .OR. T(K,J).LE.200.)    THEN
C 
                  H0(K)=0.90-0.125*(0.915-SHJ(K,J-1)) 
                  IF(SHJ(K,J-1) .LE. 0.5) H0(K)=0.85
C 
                  IF( (H(K,J) - H0(K)).GT.0.0) THEN 
                    CVEC(K) = (H(K,J) - H0(K))/(1.-H0(K)) 
                  ELSE
                    CVEC(K) = 0.
                  ENDIF 
C 
                  IF(J.EQ.LEV) CVEC(K) = 0.0
C
C........ STORE CLOUDINESS IN UPPER TRIANGULAR CLOUD MATRIX.
C  
                  C MTX(K,J,J+1) = MERGE(CVEC(K),0.,CVEC(K).GE.CUT) 
              ENDIF
   88     CONTINUE
   90 CONTINUE
  
C........ INTEGRATE LIQUID WATER PATH 
  
      ZF = 0. 
      DO 200 K = IL1, IL2 
          FCTLAY = MERGE(FCTLAYO, ZF, C MTX(K,LEV,LEVP1) .GT. CUT)
          WCL(K,LEV) = WCL(K,LEV) * FCTLAY 
          TAC(K,LEV) = 1.5 * WCL(K,LEV) / RADEQV(K,LEV)
C
C . . . . PARAMETRIZATION BY PLATT AND HARSVARDHAN (1988) FOR CIRRUS: 
C
          IF (T(K,LEV) .LT. TICE) THEN
              E(K) = 1. - EXP(-0.75 * TAC(K,LEV)) 
          ELSE
              E(K) = 1. - EXP(-EMICOEF(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 = CMTX(K,JP1,JP2)
              CTOTAL = MIN (CMTX(K,J,JP1), C, CMTX(K,JP2,JP3)) 
              CSINGL = C - CTOTAL 
              CD = MAX (C, .0001) 
              FCTLAY = MAX((CTOTAL + FCTLAYO * CSINGL)/ CD, ZF)
              WCL(K,JP1) = WCL(K,JP1) * FCTLAY 
              TAC(K,JP1) = 1.5 * WCL(K,JP1) / RADEQV(K,JP1)
              IF (T(K,JP1) .LT. TICE) THEN
                  E(K) = 1. - EXP(-0.75 * TAC(K,JP1)) 
              ELSE
                  E(K) = 1. - 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.
               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.
                 IF(ZRDM(K).EQ.0.)                                  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                   1. - (1. - C MTX(K,JM1,J)) * (1. - C REF(K)) ) 
                    C MTX(K,J,I) = MAX ( C MTX(K,JM1,I),              
     1                   1. - (1. - C MTX(K,J,JM1)) * (1. - 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.999) 
  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 
      RETURN
      END 
