      SUBROUTINE CLOUDS18(CCLD, CLD, RH, CLW, CIC, 
     1                    WCDW, WCDI, RADEQVW, RADEQVI, CLDCDD, ETA,
     3                    SHJ,  DZ,   QLWC, QIWC, T, H,   
     4                    ZCLF, ZCDN, SCLF, SCDN, SLWC, 
     5                    PSFC, TCV, BCYLOAD,
     6                    IL1, IL2, ILG, LAY, LEV, IPAM) 
C
C     * APR 26/2018 - J.Cole Removed calculation of CLWT, CICT, CLDT
C     * JUN 25/2013 - K.Vonsalzen/New version for gcm17:
C     *               J.Li.       - Implement PAM option.
C     *                           - Implement semi-direct
C     *                             effect updates. ETA calculated
C     *                             and passed out, RHO promoted to
C     *                             an internal array.
C     * FEB 12/2009 - JASON COLE.  PREVIOUS VERSION CLOUDS17 FOR
C     *                            GCM15H/GCM15I/GCM16:
C     *                            - HARD-CODED VALUE OF 75.46 FOR
C     *                              LIQUID EQUIVILENT RADIUS IS
C     *                              NOW PROPERLY DECOMPOSED INTO
C     *                              PIFAC=62.035 AND THE "TUNABLE"
C     *                              FACTOR BETA (INCREASED TO 1.3
C     *                             TO GIVE BETTER AGREEMENT WITH
C     *                             OBSERVATIONS OF LIQUID EQUIVILENT
C     *                             RADIUS). 
C     * OCT 15/2007 - M.LAZARE.    PREVIOUS VERSION CLOUDS16 FOR GCM15G:
C     *                            - DZ NOW INPUT ARRAY FROM PHYSICS
C     *                              AND NOT CALCULATED INSIDE, SO
C     *                              SHTJ NOT REQUIRED AND THUS REMOVED.
C     *                            - WCDW AND WCDI ARE NOW PHYSICS
C     *                              WORK ARRAYS AND ARE THUS PASSED
C     *                              OUT. WE USE INSTEAD WCLW AND
C     *                              WCLI WORK ARRAYS TO CALCULATE
C     *                              TOTAL CLOUD AT THE END.
C     *                            - CODE RELATED TO CLOUD OPTICAL
C     *                              PROPERTIES REMOVED.
C     *                            - COSMETIC RE-ORGANIZATION OF INPUT/
C     *                              OUTPUT FIELDS.
C     *                            - REMOVAL OF UNUSED PARAMETERS.
C     * JAN 13/2007 - M.LAZARE.    PREVIOUS VERSION CLOUDS15 FOR GCM15F.
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
C     * OUTPUT ARRAYS... 
C
      REAL   CCLD   (ILG,LAY), CLD    (ILG,LAY), RH     (ILG,LAY),
     1       CLW    (ILG,LAY), CIC    (ILG,LAY),
     2       WCDW   (ILG,LAY), WCDI   (ILG,LAY),
     3       RADEQVW(ILG,LAY), RADEQVI(ILG,LAY), CLDCDD (ILG,LAY),
     4       ETA    (ILG,LAY)


C     * INPUT ARRAYS... 

      REAL   SHJ (ILG,LAY), DZ  (ILG,LAY)
      REAL   QLWC(ILG,LAY), QIWC(ILG,LAY)
      REAL   T   (ILG,LAY), H   (ILG,LAY)   
      REAL   ZCLF(ILG,LAY), ZCDN(ILG,LAY), BCYLOAD(ILG,LAY)
      REAL   SCLF(ILG,LAY), SCDN(ILG,LAY), SLWC(ILG,LAY)
                                   
      REAL   PSFC(ILG), TCV(ILG)
C          
C     * INTERNAL WORK ARRAYS... 
C
      REAL, DIMENSION(ILG,LAY) :: WCL, WCD, WCLW, WCLI, RHO
C
      COMMON /PARAMS/ WW, TWW, RAYON, ASQ, GRAV, RGAS, RGOCP, 
     1                RGOASQ, CPRES, RGASV, CPRESV
      COMMON /PARAM4/ ACTFRC
C
C     * THERE MUST BE NO CLOUD WITHIN FIRST "LEV1" LAYERS FOR RADIATION
C     * TO WORK PROPERLY. THIS IS DEFINED IN THE "TOPLW" SUBROUTINE
C     * CALLED AT THE BEGINNING OF THE MODEL.
C
      COMMON /ITOPLW/ LEV1
C
      DATA WCMIN /0.00001/, CUT /0.001/
      DATA ZCDNMIN /1./, ZCDNMAX /500./
      DATA THIRD /0.3333333333333/
      DATA CLDMIN /1.E-9/
      DATA RIEFFMAX /50./
      DATA RWEFFMIN,RWEFFMAX /2., 30./
      DATA ZERO,ONE /0., 1./
C
C     * COEFFICIENTS FOR LIQUID EQUIVILENT RADIUS.
C     * NOTE THAT THESE VALUES ARE DEFINED IN CLD_PARTICLE_SIZE_CLD_GEN2.
C     * MAKE SURE THEY ARE CONSISTENT!
C
      DATA PIFAC / 62.035 /
      DATA BETA  / 1.3 / 
C=======================================================================
C     * INITIALIZE ARRAYS AND CALCULATE LAYER THICKNESS.
C
      DO 80 J = 1, LAY
      DO 80 K = IL1, IL2
          WCDW(K,J)=0.
          WCDI(K,J)=0.
          CCLD(K,J)=0.
          CLD (K,J)=0.
          ETA (K,J)=0.
   80 CONTINUE  

C----------------------------------------------------------------------C
C     * CLOUD AMOUNT AND WATER/ICE CONTENTS.
C----------------------------------------------------------------------C
      DO 95 J = 1, LAY
        JP1 = J + 1
        DO 85 K = IL1, IL2
          P = PSFC(K) * SHJ(K,J)     
          RHO(K,J) = P / (RGAS * T(K,J))
C
          IF ( IPAM == 1 ) THEN
            TCLF=MIN(MAX(ZCLF(K,J),ZERO),ONE)
          ELSE
            TCLF=MIN(MAX(ZCLF(K,J),SCLF(K,J),ZERO),ONE)
          ENDIF
          IF (TCLF.GE.CUT) THEN
            IF ( IPAM == 1 ) THEN
              WCDW(K,J)=QLWC(K,J)*RHO(K,J)*1000./TCLF
            ELSE
              WCDW(K,J)=(QLWC(K,J)+SLWC(K,J))*RHO(K,J)*1000./TCLF
            ENDIF
            WCDI(K,J) = QIWC(K,J) * RHO(K,J) * 1000. / TCLF
          ELSE
            WCDW(K,J)=0.
            WCDI(K,J)=0.
          ENDIF
C
C         * GENERATE LOCAL CLOUDS.
C
          CVOL=TCLF
          IF(J.GT.LEV1) CCLD(K,J)=MIN(CVOL,ONE) 

C
C         * ASSIGN CCLD AND H TO CLD AND RH
C
          CLD(K,J)  =  CCLD(K,J)
          RH (K,J)  =  H(K,J)

   85   CONTINUE
   95 CONTINUE     

C----------------------------------------------------------------------C
C     * CLOUD OPTICAL PROPERTIES.
C----------------------------------------------------------------------C
      DO 190 J = 1, LAY
        JP1 = J + 1
        DO 150 K = IL1, IL2    
C
C         * DEFINE WATER/ICE PATHS, OPTICAL DEPTH AND EMMISSIVITY.
C
          WCD(K,J)  = WCDW(K,J) + WCDI(K,J)
          CLW(K,J)  = QLWC(K,J)
          IF ( IPAM /= 1 ) CLW(K,J) = CLW(K,J) + SLWC(K,J)
          CIC(K,J)  = QIWC(K,J)
          WCLW(K,J) = WCDW(K,J) * DZ(K,J)
          WCLI(K,J) = WCDI(K,J) * DZ(K,J) 
          WCL(K,J)  = WCLW(K,J)
C          
C         * CALCULATE EQUIVALENT RADIUS (MICRONS) FOR BOTH PHASES,
C         * ASSUMING GAMMA SIZE DISTRIBUTION.
C         * NOTE THAT CLOUD DROPLET DENSITY NOW INCLUDES EFFECT OF
C         * SHALLOW CLOUDS!
C
          IF ( IPAM == 1 ) THEN
            PCDNC=(ZCLF(K,J)*ZCDN(K,J))/MAX(CCLD(K,J),CLDMIN)
          ELSE
            PCDNC=(ZCLF(K,J)*ZCDN(K,J)+SCLF(K,J)*SCDN(K,J))
     1         /MAX(CCLD(K,J),CLDMIN)
          ENDIF
          CDD=MAX(MIN(ZCDNMAX,PCDNC*1.E-6),ZCDNMIN)
          CLDCDD(K,J) = CDD
C
          IF ( IPAM == 1 ) THEN
C
C           * DISPERSION EFFECT (PENG AND LOHMANN, 2003).
C
            BETAS = (1.22+0.00084*CDD) 
          ELSE
            BETAS = BETA 
          ENDIF
          RADEQVW(K,J) = BETAS*PIFAC*(WCDW(K,J)/CDD)**THIRD 
          RADEQVW(K,J) = MAX(MIN(RADEQVW(K,J),RWEFFMAX),RWEFFMIN)
C
          WCFAC=MAX(WCDI(K,J),WCMIN)
          ZRIEFF=83.8*WCFAC**0.216
          RADEQVI(K,J) = MIN(ZRIEFF,RIEFFMAX)
C
C----------------------------------------------------------------------C
C     ETA (VOLUME RATIO OF BC IN CLOUD DROPLET) 0.9 = BC DENSITY       C
C     FACTOR 1000, SINCE WCDW IN UNITS G/M^3, WCDW IS THE GRID MEAN    C
C     RESULTS, WCDW/CCLD IS THE CLOUD LOCATED AREA MEAN                C
C----------------------------------------------------------------------C
C
          IF (WCDW(K,J) .GT. 1.E-05)                                THEN
            ETA(K,J) =  1000. * ACTFRC * RHO(K,J) * CCLD(K,J) * 
     1                  BCYLOAD(K,J) / (0.9 * WCDW(K,J))
          ELSE
            ETA(K,J) =  0.0
          ENDIF
  150   CONTINUE
  190 CONTINUE   
C
      RETURN
      END
