      SUBROUTINE COND11(ISULF,XROW,TH,QV,RH,ZCLF,ZCDN,DSHJ,SHJ,SHTJ,
     1                  ZRFL,ZSFL,PRESSG,SO4CROW,QLWC,
     2                  QIWC,TF,DQLDT,ZFEVAP,ZMRATEP,ZFSNOW,ZMLWC,
     3                  CLRFR,ZFRAIN,CVAR,CVDU,CVSG,PBLT,QC,RHC,
     4                  CLRFS,ZFSUBL,QCWVAR,
     5                  ALMC,ALMX,ZTMST,NTRAC,ILG,IL1,IL2,ILEV,LEV,LEVS,
     6                  IPAM,ISAVLS,ILWC,IIWC,ISO4,ISSA,ISSC,ICOSP,
     7                  RMIXROL, SMIXROL, RREFROL, SREFROL,
     8                  AGGROL,AUTROL,CNDROL,DEPROL,
     9                  EVPROL,FRHROL,FRKROL,FRSROL,
     A                  MLTIROL,MLTSROL,RACLROL,RAINROL,
     B                  SACIROL,SACLROL,SNOWROL,SUBROL,SEDIROL,
     C                  QTNC,HMNNC)
C
C     * FEB 25/2015 - M.LAZARE/      New version for gcm18:
C     *               K.VONSALZEN:   - Comment-out lower bound on total
C     *                                water (don't think we need it
C     *                                any more due to implemented
C     *                                improvements and bugfixes elsewhere).
C     *                              - Modifications to autoconversion
C     *                                and ice crystal fall velocity for
C     *                                both bulk and PLA aerosols.
C     *                              - Removed unused ZCR.
C     *                              - Removed unused GCROW.
C     *                              - ZRFL and ZSFL now output in
C     *                                proper MKS units, so no conversion
C     *                                done at end (loop 721 removed).
C     * NOV 14/2013 - M.LAZARE.      Previous version cond10 for gcm17:
C     *                              - UICE lowered from 5000. to 4000.
C     *                              - Autoconversion efficiency for
C     *                                PLA lowered from 1.00 to 0.65.
C     * JUL 31/2013 - M.LAZARE.      Cosmetic changes (since QTN.ge.0. in
C     *                              our tests) to:
C     *                              - Print out warning if QTN<0 before
C     *                                calling statcld5, then limit
C     *                                it to be non-negative.
C     *                              - Pass in LVL and ICALL to statcld5
C     *                                to aid in future debugging.
C     *                              - Correct calculation of QTNC,HMNNC
C     *                                for diagnostic "ISAVLS" calculations.
C     * JUN 26/2013 - K.VONSALZEN/   New version for gcm17:
C     *               M.LAZARE/      - Implement support for new PAM
C     *               J.COLE.          option (IPAM passed in).
C     *                              - No statcld above 10 mbs in order
C     *                                to prevent spurious crashes.
C     *                              - Tuning of UICE from 4000. to 5000.
C     *                              - Revised call to new STATCLD5.
C     *                              - Calculate and save extra cloud
C     *                                microphysics field, SEDIROL
C     *                                (ice crystal sedimentation).
C     *                              - Revised second indirect effect
C     *                                implemented for PAM.
C     *                              - Pass {RAIN,SNOW} to physics
C     *                                separately, instead of total.
C     *                              - MCICA integer switch removed
C     *                                (since always is true).
C     * APR 29/2012 - K.VONSALZEN/   PREVIOUS VERSION COND9 FOR GCM16:
C     *                M.LAZARE/     - ZRNFT,ZCLRR,ZCLFR,ZSNFT,ZCLRS,
C     *                J.COLE.         ZCLFS INITIALIZED IN LOOP 202 
C     *                                INSTEAD OF 2D LOOP 201.
C     *                              - TUNING CHANGES FOR FCAUT
C     *                                (2.5 INSTEAD OF 1.3) AND UICE
C     *                                (4000 INSTEAD OF 6000).
C     *                              - CALLS NEW STATCLD4 INSTEAD OF
C     *                                STATCLD3.
C     *                              - REVISED CALCULATION OF RH 
C     *                                USING SATURATION VAPOUR PRESSURE
C     *                                INSTEAD OF SATURATION SPECIFIC
C     *                                HUMIDITY.
C     *                              - PROVIDE THE GRIDBOX MEAN SNOW
C     *                                AND RAIN MIXING RATIO TO COSP.
C     *                              - STATISTICAL CLOUD SCHEME ONLY
C     *                                CALLED FOR L>LEV1 (DEFINED IN RAD
C     *                                AND PASSED IN ITOPLW COMMON BLOCK)
C     *                                AND EST<1.01*P.
C     * APR 21/2010 - K.VONSALZEN/   PREVIOUS VERSION COND8 FOR GCM15I:
C     *               M.LAZARE/      - ADD DIAGNOSTIC ARRAYS {RMIXROL,
C     *               J.COLE.          SMIXROL,RREFROL,SREFROL} FOR
C     *                                EXTRA DIAGNOSTIC OUTPUT ARRAYS
C     *                                UNDER CONTROL OF "ICOSP". 
C     *                              - SLIGHT TWEAK TO INDIRECT EFFECT
C     *                                FOR BETTER AEROSOL FORCING (PCDNC).
C     * FEB 17/2009 - K.VONSALZEN/   PREVIOUS VERSION COND7 FOR GCM15H:
C     *               M.LAZARE.      - BUGFIX TO DEFINE ZMRATEP,ZMLWC
C     *                                IN PROPER PORTION OF CODE.
C     *                              - "REMOVAL" OF SECOND INDIRECT
C     *                                EFFECT BY SETTING CONSTANT
C     *                                CLOUD DROPLET DENSITY IN 
C     *                                AUTOCONVERSION.
C     *                              - SET ACCRETION FACTOR TO UNITY
C     *                                AND TUNING OF AUTOCONVERSION
C     *                                FACTOR TO GIVE GOOD BUDGETS (FOR
C     *                                MCICA). 
C     *                              - REVISED CLOUD DROPLET DENSITY
C     *                                (PCDNC) CALCULATION, CONSISTENT
C     *                                WITH SHALLOW CONVECTION.
C     * DEC 20/2007 - K.VONSALZEN/   PREVIOUS VERSION COND6 FOR GCM15G:
C     *               M.LAZARE.      - SWITCH "MCICA" PASSED IN AND USED
C     *                                TO DEFINE FACTOR IN FRONT OF
C     *                                AUTOCONVERSION AND AUTOCORRELATION
C     *                                (LARGER FOR MCICA).
C     *                              - MANY BUGFIXES, SUCH AS IN EVAPORATION
C     *                                OF RAIN AND SUBLIMATION OF SNOW, TO
C     *                                PROPERLY CONSERVE.
C     *                              - USE THE ROECKEL WATER VAPOUR 
C     *                                PROBABILITY CURVE.
C     *                              - CALLS NEW STATCLD3.
C     *                              - USES CLEAR-SKY AND ALL-SKY MIXING
C     *                                LENGTHS (ALMC AND ALMX, RESPECTIVELY)
C     *                                TO DEFINE A CONSISTENT TOTAL
C     *                                MIXING LENGTH (ALMIX) PASSED TO
C     *                                STATCLD3.
C     *                              - PASSES IN ADELT=2.*DELT FROM PHYSICS
C     *                                AS ZTMST WHICH IS USED THROUGHOUT
C     *                                ROUTINE IN PLACE OF HARD-CODED 2.*DELT.
C     *                              - NOW DEFINES {QLWC,QIWC} AT THE END
C     *                                OF THE ROUTINE INSTEAD OF AT THE
C     *                                BEGINNING, FOR SUBSEQUENT USE IN THE
C     *                                PHYSICS.
C     * MAR 27/2007 - M.LAZARE/      PREVIOUS FINAL VERSION COND5 FOR GCM15F.
C     * JAN 15/2007 - K.VONSALZEN.   - UICE INCREASED FROM 2000 TO 4000.
C     *                              - MODIFIED CALL TO STATCLD2, IN
C     *                                CONJUNCTION WITH CHANGES TO ADD
C     *                                "QCWVAR".  
C     *                              - INITIALIZE ZSAUT,ZSACI,FRSTEMP
C     *                                FOR "XTRALS" OPTION.
C     * NOV 28/2006 - M.LAZARE.      NEW ENHANCEMENTS BASED ON GCM15E.2:
C     *                              - CORRECT CALCULATION OF CLEAR-SKY
C     *                                RELATIVE HUMIDITY.
C     *                              - CORRECT LARGE-SCALE BUDGETS
C     *                                (%DF XTRALS"). 
C     * JUN 19/2006 - M.LAZARE.      NEW VERSION FOR GCM15F:
C     *                              - USE VARIABLE INSTEAD OF 
C     *                                CONSTANT IN INTRINSICS SUCH AS "MAX",
C     *                                SO THAT CAN COMPILE IN 32-BIT MODE
C     *                                WITH REAL*8.  
C     *                              - UNUSED "EPSICE" COMMON BLOCK 
C     *                                REMOVED.
C     *                              - WORK ARRAYS NOW LOCAL.
C     *                              - CALLS NEW VERSION STATCLD2. 
C     *                              - DEFINE ALMIX,DHLDZ,DRWDZ,RRL,CPM
C     *                                OUTSIDE OF CONDITIONAL TEST ON
C     *                                ES<P, TO AVOID NAN'S ENTERING
C     *                                STATCLD.
C     * MAY 06/2006 - K.VONSALZEN/   PREVIOUS VERSION COND4 FOR GCM15E:
C     *               M.LAZARE.    - USE SUBROUTINE STATCLD TO PERFORM
C     *                              UNRAVELLING AND STATISTICAL CLOUD
C     *                              CALCULATIONS.
C     *                            - DEFINE ESAT OVER COMPLETE VERTICAL
C     *                              DOMAIN (SPECIAL HANDLING FOR ES>P)
C     *                              SO THAT CAN REMOVE IF CONDITIONS
C     *                              ON CHTOP PASSED FROM PHYSICS (NO
C     *                              LONGER NEEDED) AND DON'T REQUIRE
C     *                              "MTOP".
C     *                            - CLEANUP TO REMOVE ALL OLD UNNECESSARY
C     *                              ARRAYS.
C     *                            - INCLUDE CALCULATION OF CLEAR-SKY
C     *                              RELATIVE HUMIDITY (RHC) AND SPECIFIC
C     *                              HUMIDITY (QC) WHICH EVOLVE FROM
C     *                              CLEAR-SKY SATURATION SPECIFIC 
C     *                              HUMIDITY PASSED OUT FROM STATISTICAL
C     *                              CLOUD SCHEME (SSH).
C
C**** *COND4* - COMPUTES LARGE-SCALE WATER PHASE CHANGES AND CLOUD COVER. 
C                                                                        
C     SUBJECT.                                                           
C     --------                                                           
C                                                                        
C          THIS ROUTINE COMPUTES THE PHYSICAL TENDENCIES OF THE FOUR
C     PROGNOSTIC VARIABLES T, Q, AND XL, XI  DUE TO WATER PHASE  
C     CHANGES (CONDENSATION, EVAPORATION OF FALLING PRECIPITATION IN     
C     UNSATURATED LAYERS AND MELTING OR FREEZING OF THE FALLING WATER)   
C     AND PRECIPITAION FORMATION (AUTOCONVERSION AND ACCRETION IN THE
C     WARM PHASE AND AGGREGATION AND ACCRETION IN THE COLD PHASE).
C                                                                        
C     RAIN, SNOWFALL, SURFACE FLUXES, CLOUD WATER AND CLOUD COVER        
C     LATER TO BE USED FOR SOIL PROCESSES AND RADIATION ARE STORED.      
C                                                                        
C**   INTERFACE.                                                         
C     ----------                                                         
C                                                                        
C     *CALL* *COND*                                                      
C                                                                        
C     INPUT ARGUMENTS.                                                   
C     ----- ----------                                                   
C  - 2D                                                                  
C  QV     : SPECIFIC HUMIDITY 
C  TH     : TEMPERATURE 
C  XROW(X,X,ILWC)    : CLOUD WATER (KG/KG)
C  XROW(X,X,IIWC)    : CLOUD ICE   (KG/KG)                                    
C  SHJ    : SIGMA LEVELS
C
C  PRESSG  : SURFACE PRESSURE
C  PCP   : STRATIFORM PRECIPITAION
C                                                                        
C                                                                        
C     REFERENCES.                                                        
C     ----------                                                         
C                                                                        
C     LOHMANN AND ROECKNER, 1995: MPI-REPORT NO. 179
C     LEVKOV ET AL., 1992: BEITR. PHYS. ATMOS. 35-58. (ICE PHASE)
C     BEHENG, 1994: ATMOS. RES. 193-206.              (WARM PHASE)
C     ECHAM3/4-DESCRIPTION, 1992: MPI-REPORT NO. 93   (CLOUD COVER,
C                   CONDENSATION,EVAPORATION)
C                                                                        
C     AUTHOR.                                                            
C     -------                                                            
C     U.LOHMANN     MPI-HAMBURG  1995     
C                                                                        
C     MODIFICATIONS.                                                     
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
C     * ARRAYS SHARED WITH PHYSICS.
C
      REAL   XROW(ILG,LEV,NTRAC)
      REAL   TH  (ILG,ILEV), QV  (ILG,ILEV), RH  (ILG,ILEV),
     1       ZCLF(ILG,ILEV), ZCDN(ILG,ILEV), TF(ILG,ILEV),
     2       DSHJ(ILG,ILEV), SHJ (ILG,ILEV), SHTJ(ILG,LEV)
      REAL   ZRFL   (ILG), ZSFL  (ILG), PRESSG(ILG)
      REAL   SO4CROW(ILG,ILEV),QLWC(ILG,ILEV),QIWC(ILG,ILEV)
C
      REAL    DQLDT  (ILG,ILEV), ZFEVAP(ILG,ILEV),
     1        ZMRATEP(ILG,ILEV), ZFSNOW(ILG,ILEV), 
     2        ZMLWC  (ILG,ILEV), QC    (ILG,ILEV),
     3        RHC    (ILG,ILEV), CVAR  (ILG,ILEV),
     4        ALMX  (ILG,ILEV), ALMC   (ILG,ILEV),
     4        CVDU   (ILG,ILEV), CVSG  (ILG,ILEV),
     5        QCWVAR (ILG,ILEV), CLRFR (ILG,ILEV),
     6        ZFRAIN (ILG,ILEV), CLRFS (ILG,ILEV),
     7        ZFSUBL (ILG,ILEV)
      REAL    PBLT   (ILG)

C
C     * DIAGNOSTIC OUTPUT ARRAYS, UNDER CONTROL OF "ISAVLS".
C
      REAL   AGGROL (ILG,ILEV),  AUTROL (ILG,ILEV),  CNDROL (ILG,ILEV)
      REAL   DEPROL (ILG,ILEV),  EVPROL (ILG,ILEV),  FRHROL (ILG,ILEV)
      REAL   FRKROL (ILG,ILEV),  FRSROL (ILG,ILEV),  MLTIROL(ILG,ILEV)
      REAL   MLTSROL(ILG,ILEV),  RACLROL(ILG,ILEV),  RAINROL(ILG,ILEV)
      REAL   SACIROL(ILG,ILEV),  SACLROL(ILG,ILEV),  SNOWROL(ILG,ILEV)
      REAL   SUBROL (ILG,ILEV),  SEDIROL(ILG,ILEV),  HMNNC  (ILG,ILEV), 
     1       QTNC   (ILG,ILEV)
C
C     * DIAGNOSTIC OUTPUT ARRAYS, UNDER CONTROL OF "ICOSP".
C
      REAL, DIMENSION(ILG,ILEV), INTENT(OUT) ::
     1 RMIXROL, ! Rain mixing ratio (kg/kg)
     2 SMIXROL, ! Snow mixing ratio (kg/kg)
     3 RREFROL, ! Effective raindrop size (meters) 
     4 SREFROL  ! Effective snow particle size (meters) 
C
C     * INTERNAL WORK ARRAYS.
C
      REAL  ,  DIMENSION(ILG,ILEV)  :: ZCRAUT,ZFRAC,ZXIB,Z,P,HMNN,QTN
      REAL  ,  DIMENSION(ILG) :: ZDP,ZEVP,ZFRL,ZIMLT,ZSMLT,ZSACL
      REAL  ,  DIMENSION(ILG) :: ZRPR,ZSPR,ZSUB,ZXLB
      REAL  ,  DIMENSION(ILG) :: ZRFLN,ZSFLN,ZQNL,ZTFALL,ZFALLO,ZFALLI
      REAL  ,  DIMENSION(ILG) :: ZSACI,ZSAUT,QCW,SSH,RRL,CPM,ALMIX
      REAL  ,  DIMENSION(ILG) :: DHLDZ,DRWDZ,PMBS,FLAGEST,RAT,VTR,ZRNFT
      REAL  ,  DIMENSION(ILG) :: ZSNFT,ZCLRR,ZCLRS,ZCLFR,ZCLFS,ZCLFRM
      REAL  ,  DIMENSION(ILG) :: ZCLFSM,ZRCLR,ZSCLR
C
      LOGICAL LO,LO1,LO2
      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 /PARAM1/ PI,     RVORD, TFREZ, HS,   HV,   DAYLNT
      COMMON /EPS   / A,B,EPS1,EPS2    
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
C     * PARAMETERS USED IN NEW SATURATION VAPOUR PRESSURE FORMULATION.
C
      COMMON /ESTWI/ RW1,RW2,RW3,RI1,RI2,RI3
                                                                      
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.-UUU)*TI(TTT)               
                                                                     
C     * C) COMPUTES THE SATURATION VAPOUR PRESSURE OVER WATER OR ICE.
C
      ESW(TTT)    = EXP(RW1+RW2/TTT)*TTT**RW3
      ESI(TTT)    = EXP(RI1+RI2/TTT)*TTT**RI3
c     ESTEFF(TTT,UUU) = (1.-UUU)*ESW(TTT) + UUU*ESI(TTT)   
C
      DATA ZERO,ONE /0., 1./
      DATA ZRIMIN,ZRFLNL_MAX/1.E-5, 0.028/
C=======================================================================
C
C*    SECURITY PARAMETERS.
C     --------------------
C
      ZSEC=1.E-9
      ZSECFRL=1.E-7
      ZEPCLC=1.E-2
C
C*    COMPUTATIONAL CONSTANTS.
C     ------------- ----------
C
C         1. INPUT
 100  CONTINUE
C
C     * XU AND RANDALL CONSTANTS (COMMENTED-OUT ONES ARE THEIR
C     * ORIGINAL VALUES!).
C
      ZCGAMMA=-0.49
      ZCXU=100.
      ZCP=0.25
C
      API=3.1415926535897
      RHOH2O=1000.
C          RHOSNO ---> BULK DENSITY OF SNOW
      RHOSNO=100.
      CP=1004.5
      ZMELTP2=T1S+2.
      PREF=101325.
      ZRT=0.9
      NEX=4
      ZRS=0.975
      ZCLOIA=1.0E+02
      ZRHSC=0.4
      ZSATSC=0.8
      ZTHOMI=238.16
      ZN0S=3.E6
      ZRHOI=500.
C
C
C     * DEFINE FACTORS FOR AUTOCONVERSION, ICE CRYSTAL FALL VELOCITY,
C     * ACCRETION AND SECOND INDIRECT EFFECT, BASED ON WHETHER PAM 
C     * APPROACH IS BEING USED OR NOT.
C
      IF ( IPAM .NE. 1 ) THEN
        I2NDIE=1
        COAUT=3.
        FACAUT=1.3E+09/1.08E+10
        FACACC=15.
        UICEFAC=4200.
      ELSE
        I2NDIE=1
        COAUT=3.
        FACAUT=2.E+09/1.08E+10
        FACACC=10.
        UICEFAC=4200.
      ENDIF

      COACC=1.15
      COTMP1=COAUT-2.
      COTMP2=COACC-1.
      ZCSAUT=50.
      ZCDNMIN=1.E6
      ZQTMST=1./ZTMST
      ZCONS2=1./(ZTMST*GRAV)
      ZEAC=0.7
      WCMIN=0.00001
      EPSLIM=0.001
C
C
C          SCE    ---> SCHMIDT NUMBER
      SCE= 0.64
C          ZNU     ---> KINEMATIC VISCOSITY OF AIR
      ZNU= 0.1346E-4
C          THERCO ---> THERMAL CONDUCTIVITY OF AIR (AT 0 CELSIUS)
      THERCO= 2.43E-2
C          PSI    ---> WATER VAPOR DIFFUSIVITY 
      PSI= 0.211E-4
C          NZEROS ---> INTERCEPT PARAMETER FOR SNOW DISTRIBUTION
      ZNZEROS= 3.0E6
C          A19    ---> C-PARAMETER IN SNOW FALL VELOCITY
      A19= 4.8361
      FACT = 0.4495*(SCE**0.33)*((A19/ZNU)**0.5)
      YFR=0.50
      EXPM1=EXP(-1.)
      FACTE=1./EPS1-1.
C          RNZEROS ---> INTERCEPT PARAMETER FOR RAIN DISTRIBUTION
      RNZEROS= 8.0E6 ! m-4
C
C*         2.  QUANTITIES NEEDED FOR
C*             ---------- ------ ---
C*                CONDENSATION AND PRECIPITATION CALCULATIONS.
C*                ------------ --- ------------- ------------
C
C*         2.0     INITIALIZE FIELDS.
C
      DO 201 JK=1,ILEV
      DO 201 IL=IL1,IL2
        RH     (IL,JK)=0.
        ZCLF   (IL,JK)=0.
        ZFEVAP (IL,JK)=0.
        ZMRATEP(IL,JK)=0.
        ZFSNOW (IL,JK)=0.
        ZMLWC  (IL,JK)=0.
        ZCRAUT (IL,JK)=1.
        QC     (IL,JK)=0.
        RHC    (IL,JK)=0.
        ZFSUBL (IL,JK)=0.
        ZFRAIN (IL,JK)=0.
        CLRFR  (IL,JK)=0.
        CLRFS  (IL,JK)=0.
        IF(ISAVLS.NE.0) THEN
          RAINROL(IL,JK)=0.
          SNOWROL(IL,JK)=0.
        ENDIF
  201 CONTINUE
C
      DO 202 IL=IL1,IL2
        ZRFL(IL)=0.
        ZSFL(IL)=0.
        ZSMLT(IL)=0.
        ZRFLN(IL)=0.
        ZSFLN(IL)=0.
        ZFALLI(IL)=0.
        ZFALLO(IL)=0.
        ZRNFT(IL)=0.
        ZCLRR(IL)=1.
        ZCLFR(IL)=0.
        ZSNFT(IL)=0.
        ZCLRS(IL)=1.
        ZCLFS(IL)=0.
 202  CONTINUE
C
C*         2.1    CALCULATE WATER/ICE FRACTION.
C
      MSG=ILEV-LEVS
      DO 210 JK=ILEV,1,-1
      DO 210 IL=IL1,IL2
        P(IL,JK)   = SHJ(IL,JK) * PRESSG(IL)
        IF( XROW(IL,JK+1,IIWC).GT.ZSECFRL ) THEN
          ZFRAC(IL,JK)=1./(1.+XROW(IL,JK+1,ILWC)/XROW(IL,JK+1,IIWC))
        ELSE
          IF ( XROW(IL,JK+1,ILWC).GT.ZSECFRL ) THEN
            ZFRAC(IL,JK)=0.
          ELSE
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 = MERGE( 1.,    
     1                  0.0059+0.9941*EXP(-0.003102*(T1S-TH(IL,JK))**2),
     2                  TH(IL,JK).GE.T1S )  
            ZFRAC(IL,JK)=1.-FRACW
          ENDIF
        ENDIF 
  210 CONTINUE
C                                                                       
C     * CALCULATE LOCAL HEIGHT (M) AND VARIANCE.
C                             
      ROG=RGAS/GRAV                                                     
      DO 30 IL=IL1,IL2                                                  
        Z(IL,ILEV)     = ROG*TF(IL,ILEV)*LOG(PRESSG(IL)/P(IL,ILEV))  
   30 CONTINUE                                                          
      ILEVM=ILEV-1 
C
      DO 40 JK=ILEVM,1,-1
      DO 40 IL=IL1,IL2                                 
        Z(IL,JK)=Z (IL,JK+1)+ROG*TF(IL,JK)*LOG(P(IL,JK+1)/P(IL,JK))   
   40 CONTINUE
C
      DSR=1.
      CSIGMA=0.2
      DO 50 JK=1,ILEV
        JKP = JK+1
        DO 45 IL=IL1,IL2                                 
C
C         * TOTAL WATER MIXING RATIO AND LIQUID WATER STATIC ENERGY,
C         * AND DRY STATIC ENERGY.
C 
          QCWX=(XROW(IL,JKP,ILWC)+XROW(IL,JKP,IIWC))*DSR
          QR=QV(IL,JK)*DSR
          QTN(IL,JK)=QR+QCWX
          IF ( QTN(IL,JK).LE.0. ) CALL WRN('COND10',-1)
c         QTN(IL,JK)=MAX(QTN(IL,JK),0.) 
          DST=CPRES*TH(IL,JK)+GRAV*Z(IL,JK)
          RRLX=(1.-ZFRAC(IL,JK))*HV+ZFRAC(IL,JK)*HS
          HMNN(IL,JK)=DST-RRLX*QCWX
   45   CONTINUE
   50 CONTINUE
C
C     * CLOUD FRACTION, EVAPORATION, AND CONDENSATION FROM STATISTICAL
C     * CLOUD SCHEME.
C
      ICVSG=1
      ISUBG=1
      DO 240 JK=1,ILEV
       JKP = JK+1
       DO 200 IL=IL1,IL2
        PMBS(IL)=0.01*P(IL,JK)
        RRL(IL)=(1.-ZFRAC(IL,JK))*HV+ZFRAC(IL,JK)*HS
        CPM(IL)=CPRES
C
        ALMIX(IL)=MAX(ALMC(IL,JK),ALMX(IL,JK),10.)
C
        IF(JK.GT.1)                              THEN
          DZ       =Z   (IL,JK-1)-Z  (IL,JK)
          DHLDZ(IL)=(HMNN(IL,JK-1)-HMNN(IL,JK))/DZ
          DRWDZ(IL)=(QTN (IL,JK-1)-QTN (IL,JK))/DZ
        ELSE
          DHLDZ(IL)=0.
          DRWDZ(IL)=0.
        ENDIF  
C
        TX=TH(IL,JK)
        EST=(1.-ZFRAC(IL,JK))*ESW(TX)+ZFRAC(IL,JK)*ESI(TX)
        IF(JK.GT.LEV1 .AND.
     1     EST.LT.PMBS(IL) .AND. PMBS(IL).GE.10.)    THEN
          FLAGEST(IL)=1.
        ELSE
          FLAGEST(IL)=0.
        ENDIF
  200  CONTINUE
C
       ICALL=1
       CALL STATCLD5(QCW,ZCLF(1,JK),CVAR(1,JK),ZCRAUT(1,JK),
     1               QCWVAR(1,JK),SSH,CVSG(1,JK),QTN(1,JK),
     2               HMNN(1,JK),ZFRAC(1,JK),CPM,
     3               PMBS,Z(1,JK),RRL,CVDU(1,JK),FLAGEST,
     4               ALMIX,DHLDZ,DRWDZ,I2NDIE,
     5               GRAV,ZTMST,ILEV,ILG,IL1,IL2,ICVSG,ISUBG,
     6               JK,ICALL                               )
C
       DO 220 IL=IL1,IL2
        QCWI=ZFRAC(IL,JK)*QCW(IL)/DSR
        QCWL=(1.-ZFRAC(IL,JK))*QCW(IL)/DSR
C
C       * CONDENSATION AND EVAPORATION.
C
        ZCND=QCWL-XROW(IL,JKP,ILWC)
        ZDEP=QCWI-XROW(IL,JKP,IIWC)
        ZCND=MAX(ZCND,-XROW(IL,JKP,ILWC))
        ZDEP=MAX(ZDEP,-XROW(IL,JKP,IIWC))
        ASUM=ZCND+ZDEP
        ASCL=1.
        IF ( ASUM.GT.0. ) THEN
          ASCL=MIN(MAX((QV(IL,JK)-ZSEC)/ASUM,0.),1.)
        ENDIF
        ZLVDCP=HTVOCP(TH(IL,JK),1.)
        ZLSDCP=HTVOCP(TH(IL,JK),0.)
        ZCND=ASCL*ZCND
        ZDEP=ASCL*ZDEP
C
C       * UPDATE TEMPERATURE, VAPOUR, AND CONDENSATE.
C
        TH(IL,JK)=TH(IL,JK)+ZCND*ZLVDCP+ZDEP*ZLSDCP
        QV(IL,JK)=QV(IL,JK)-ZCND-ZDEP
        XROW(IL,JKP,ILWC)=XROW(IL,JKP,ILWC)+ZCND
        XROW(IL,JKP,IIWC)=XROW(IL,JKP,IIWC)+ZDEP
        IF(ISAVLS.NE.0) THEN
          CNDROL (IL,JK)=ZCND*ZQTMST
          DEPROL (IL,JK)=ZDEP*ZQTMST
        ENDIF
C
C       * RELATIVE HUMIDITY FOR ALL- AND CLEAR-SKY. THE RELATIVE
C       * HUMIDITY IS CALCULATED BASED ON THE VAPOUR PRESSURE IF THE
C       * SATURATION VAPOUR PRESSURE EXCEEDS THE TOTAL AIR PRESSURE,
C       * I.E. FOR CONDITIONS THAT DO NOT PERMIT FORMATION OF CLOUDS.
C       * OTHERWISE, THE RELATIVE HUMIDITY IS CALCULATED BASED ON
C       * SPECIFIC HUMIDITY (SEE CHAPTER 4 IN EMANUEL'S TEXTBOOK ON
C       * CONVECTION) IN ORDER TO PERMIT CALCULATIONS FOR CLEAR-SKY.
C
        IF (FLAGEST(IL).LE.0.)                         THEN
          QC(IL,JK)=QV(IL,JK)
          TX=TH(IL,JK)
          EST=(1.-ZFRAC(IL,JK))*ESW(TX)+ZFRAC(IL,JK)*ESI(TX)
          EV=PMBS(IL)*QV(IL,JK)/((1.-QV(IL,JK))*EPS1+QV(IL,JK))
          RHC(IL,JK)=MIN(MAX(EV/EST,0.),1.)
          RH (IL,JK)=RHC(IL,JK)
        ELSE
          IF(ZCLF(IL,JK).LT..99) THEN
            QC (IL,JK)=(QV(IL,JK)-ZCLF(IL,JK)*SSH(IL))/(1.-ZCLF(IL,JK))
            RHC(IL,JK)=MIN(MAX((QC(IL,JK)*(1.+SSH(IL)*FACTE))
     1                      /(SSH(IL)*(1.+QC(IL,JK)*FACTE)),0.),1.)
            RH (IL,JK)=ZCLF(IL,JK)+(1.-ZCLF(IL,JK))*RHC(IL,JK)
          ELSE
            QC (IL,JK)=QV(IL,JK)
            RH (IL,JK)=1.
            RHC(IL,JK)=RH(IL,JK)
          ENDIF
        ENDIF
C
C       * DIAGNOSE LIQUID WATER STATIC ENERGY AND TOTAL WATER
C       * IN THE CLEAR-SKY ENVIRONMENT OF THE CLOUD.
C
        QTNC(IL,JK)=QC(IL,JK)*DSR
        DST=CPRES*TH(IL,JK)+GRAV*Z(IL,JK)
        HMNNC(IL,JK)=DST
  220  CONTINUE
  240 CONTINUE 
C
      IF ( IPAM .NE. 1 ) THEN
C
C       * MASS SCALING FACTOR. ASSUME THAT ALL SULPHATE AEROSOL
C       * IS IN THE FORM OF SO4.
C     
        SO4MASS=96.058/32.064
        SSAMASS=22.990/58.443
        SSCMASS=22.990/58.443
C
C       * CALCULATE CLOUD DROPLET NUMBER CONCENTRATION (IN DROPLETS/M**3)
C       * USING SEMI-EMPIRICAL PARAMETERIZATION.
C
        DO 300 L=1,ILEV
        DO 300 IL=IL1,IL2
          ZRHO0=P(IL,L)/(RGAS*TH(IL,L))
          CSO4=XROW(IL,L+1,ISO4)*SO4MASS*ZRHO0*1.E+09
          PCDNC=100.*CSO4**0.2
          PCDNC=PCDNC*0.6
          ZCDN(IL,L)=MAX(1.E+06*PCDNC,ZCDNMIN)
  300   CONTINUE
      ENDIF
C
C     * SEDIMENTATION OF ICE CRYSTALS.
C
      DO JK=LEV1,ILEV 
        JKP=JK+1
        DO IL=IL1,IL2
C
C         * DENSITY OF WET AIR AND LAYER DEPTH.
C
          ZRHO0=P(IL,JK)/(RGAS*TH(IL,JK))
          DZ=LOG(SHTJ(IL,JK+1)/SHTJ(IL,JK))*RGAS*TH(IL,JK)/GRAV
C
C         * IN-CLOUD ICE CONTENT FOR EQUIVILENT RADIUS
C         * CALCULATION (KEPT AS SCALAR TO AVOID UNNECESSARY
C         * MEMORY MOVEMENT!)
C
          ZXIBIC=XROW(IL,JKP,IIWC)
          IF(ZCLF(IL,JK).GT.ZEPCLC.AND.XROW(IL,JKP,IIWC).GT.ZSEC)
     1       ZXIBIC=XROW(IL,JKP,IIWC)/ZCLF(IL,JK)
C
C         * ICE PARTICLE SIZE.
C
          WCFAC=MAX(ZXIBIC*ZRHO0*1000.,WCMIN)
          ZRIEFF0=83.8*WCFAC**0.216
          ZRIH0=-2261.236+SQRT(5113188.044+2809.*ZRIEFF0**3)
          ZRI=MAX(ZRIMIN,1.E-6*ZRIH0**(1./3.))
C
C         * TERMINAL FALL VELOCITY.
C
          UICE =UICEFAC*2.*ZRI*(1.3/ZRHO0)**0.35*ZRHO0
C
C         * ICE TENDENCY FROM FALLING OF ICE.
C
          ZFALLO(IL)=ZFALLI(IL)
          ZFALLI(IL)=MIN(UICE,.99*DZ/ZTMST)*XROW(IL,JKP,IIWC)*ZRHO0
          IF(JK.EQ.ILEV) ZFALLI(IL)=0.  !ensure zero flux at surface
          ZTFALL(IL)=GRAV*(ZFALLO(IL)-ZFALLI(IL))
     1              /(DSHJ(IL,JK)*PRESSG(IL))
C
C         * UPDATE ICE MIXING RATIO.
C
          XROW(IL,JKP,IIWC)=XROW(IL,JKP,IIWC)+ZTFALL(IL)*ZTMST
          IF(ISAVLS.NE.0) THEN
             SEDIROL(IL,JK) = ZTFALL(IL)
          END IF
        ENDDO
      ENDDO
C
      DO 700 JK=LEV1,ILEV
      JKP=JK+1
      ZSR=0.
      ZSS=0.
      DO 410 IL=IL1,IL2
      ZLVDCP=HTVOCP(TH(IL,JK),ONE)
      ZLSDCP=HTVOCP(TH(IL,JK),ZERO)
      ZLFDCP=ZLSDCP-ZLVDCP
      ZFRL(IL)=0.
      ZRPR(IL)=0.
      ZSPR(IL)=0.
      ZIMLT(IL)=0.
      ZSACL(IL)=0.
      ZSAUT(IL)=0.
      ZSACI(IL)=0.
      RAT(IL)=1.
C
C     * RAIN AREA FRACTION FROM MAXIMUM/RANDOM CLOUD OVERLAP ASSUMPTION.
C     * ONLY CLOUDS THAT ARE SUBJECT TO WARM RAIN PROCESSES ARE 
C     * CONSIDERED.
C
       ZCLFRM(IL)=ZCLFR(IL)
       IF ( ZRFLN(IL) .GT. 1.E-20 ) THEN
         ZCLFR(IL)=ZCLF(IL,JK)
       ELSE
         ZCLFR(IL)=0.
       ENDIF
       IF ( ZCLFRM(IL).LT.1. ) THEN
         ZCLRR(IL)=ZCLRR(IL)*(1.-MAX(ZCLFR(IL),ZCLFRM(IL)))
     1                      /(1.-ZCLFRM(IL))
       ENDIF
       ZRNFT(IL)=1.-ZCLRR(IL)
C
C      * FRACTION OF GRID CELL FOR WHICH CLEAR-SKY IS AFFECTED BY RAIN.
C
       ZRCLR(IL)=MIN(MAX(ZRNFT(IL)-ZCLFR(IL),0.),1.)
C
C     * SNOW AREA FRACTION FROM MAXIMUM/RANDOM CLOUD OVERLAP ASSUMPTION.
C     * ONLY CLOUDS THAT ARE AFFECTED BY SNOW ARE CONSIDERED.
C
       ZCLFSM(IL)=ZCLFS(IL)
       IF ( ZSFLN(IL) .GT. 1.E-20 ) THEN
         ZCLFS(IL)=ZCLF(IL,JK)
       ELSE
         ZCLFS(IL)=0.
       ENDIF
       IF ( ZCLFSM(IL).LT.1. ) THEN
         ZCLRS(IL)=ZCLRS(IL)*(1.-MAX(ZCLFS(IL),ZCLFSM(IL)))
     1                      /(1.-ZCLFSM(IL))
       ENDIF
       ZSNFT(IL)=1.-ZCLRS(IL)
C
C      * FRACTION OF GRID CELL FOR WHICH CLEAR-SKY IS AFFECTED BY SNOW.
C
       ZSCLR(IL)=MIN(MAX(ZSNFT(IL)-ZCLFS(IL),0.),1.)
C
C*       3.0    CALCULATION OF CLOUD WATER AND CLOUD ICE IN THE
C*              CLOUDY AND CLOUD FREE PART OF THE GRID BOX FOR T-1
C
      ZXLB(IL)=XROW(IL,JK+1,ILWC)
      ZXIB(IL,JK)=XROW(IL,JK+1,IIWC)

      ZXLMA=ZXLB(IL)
      ZXIMA=ZXIB(IL,JK)
      ZQNL(IL)=ZCDN(IL,JK)

      IF(ZCLF(IL,JK).GT.ZEPCLC.AND.XROW(IL,JK+1,ILWC).GT.ZSEC) THEN
       ZXLB(IL)=XROW(IL,JK+1,ILWC)/ZCLF(IL,JK)
       ZXLMA=0.
      ENDIF
      IF(ZCLF(IL,JK).GT.ZEPCLC.AND.XROW(IL,JK+1,IIWC).GT.ZSEC) THEN
       ZXIB(IL,JK)=XROW(IL,JK+1,IIWC)/ZCLF(IL,JK)
       ZXIMA=0.
      ENDIF
C                                                                       
C*       3.1    MELTING OF INCOMING SNOW AND
C
      ZDP(IL)=DSHJ(IL,JK)*PRESSG(IL)
      ZCONS=ZCONS2*ZDP(IL)/ZLFDCP
      ZSNMLT=MIN(ZSFL(IL),ZCONS*MAX(ZERO,(TH(IL,JK)-ZMELTP2)))
      ZSMLT(IL)=ZSMLT(IL)+ZSNMLT
      ZRFLN(IL)=ZRFL(IL)+ZSNMLT
      ZSFLN(IL)=MAX(ZSFL(IL)-ZSNMLT,ZERO)
      ZDT=-ZSNMLT/ZCONS
      TH(IL,JK)=TH(IL,JK)+ZDT
      IF(ISAVLS.NE.0) THEN
        MLTSROL(IL,JK)=ZSNMLT*GRAV/ZDP(IL)
      ENDIF
C
C*       3.2    INSTANTENOUS MELTING OF CLOUD ICE AND
C
      IF (TH(IL,JK).GE.T1S) THEN
       ZIMLT(IL)=ZXIB(IL,JK)*ZCLF(IL,JK)
       TH(IL,JK)=TH(IL,JK)-ZLFDCP*ZIMLT(IL)
       ZXLB(IL)=ZXLB(IL)+ZXIB(IL,JK)
       ZXIB(IL,JK)=0.
       IF(ISAVLS.NE.0) THEN
         MLTIROL(IL,JK)=ZIMLT(IL)*ZQTMST
       ENDIF
      ENDIF 
C
C*       3.3    INSTANTENOUS FREEZING OF CLOUD WATER
C
      IF (TH(IL,JK).LE.ZTHOMI) THEN
       ZIMLT(IL)=-ZXLB(IL)*ZCLF(IL,JK)
       TH(IL,JK)=TH(IL,JK)-ZLFDCP*ZIMLT(IL)
       ZXIB(IL,JK)=ZXIB(IL,JK)+ZXLB(IL)
       ZXLB(IL)=0.
       IF(ISAVLS.NE.0) THEN
         FRHROL (IL,JK)=-ZIMLT(IL)*ZQTMST
       ENDIF
      ENDIF
      ZXLB(IL)=ZXLB(IL)+(1.-ZFRAC(IL,JK))*DQLDT(IL,JK)
      ZXIB(IL,JK)=ZXIB(IL,JK)+ZFRAC(IL,JK)*DQLDT(IL,JK)
C
C*       3.6    FREEZING OF CLOUD WATER BETWEEN 238 AND 273 K
C
      ZRHO0=P(IL,JK)/(RGAS*TH(IL,JK))
      LO=ZXLB(IL)*ZCLF(IL,JK).GE.ZSEC
     1   .AND.TH(IL,JK).LT.T1S.AND.TH(IL,JK).GT.ZTHOMI
      IF (LO) THEN
       ZFRL(IL)=100.*(EXP(0.66*(T1S-TH(IL,JK)))-1.)*
     1  ZRHO0*0.001/ZQNL(IL)*ZTMST
C******************************************************************
C      * NOTE THAT THE FOLLOWING EXPRESSION IS ESSENTIALLY
C      * BOUNDING ZFRL SO THAT IT APPROACHES ZXLB AS IT GETS
C      * BIG ENOUGH (IE DOESN'T REMOVE MORE THAN ZXLB).
C      * WE PROBABLY SHOULD REPLACE THIS BY THE EXPONENTIAL
C      * EXPRESSION USED ELSEWHERE (TO TEST?)
C
C*******************************************************************
       ZFRL(IL)=ZFRL(IL)*ZXLB(IL)*ZXLB(IL)/(1.+ZFRL(IL)*ZXLB(IL))
       IF(ISAVLS.NE.0) THEN
         FRKTEMP=ZFRL(IL)*ZQTMST
       ENDIF
       ZRADL=(0.75*ZXLB(IL)*ZRHO0/(API*RHOH2O*ZQNL(IL)))**(1./3.)
       ZF1=4.*API*ZRADL*ZQNL(IL)*2.E5*(270.16-TH(IL,JK))/ZRHO0
       ZF1=MAX(ZERO,ZF1)
       ZFRL(IL)=ZFRL(IL)+ZTMST*1.4E-20*ZF1
       ZFRL(IL)=MAX(ZERO,MIN(ZFRL(IL),ZXLB(IL)))
       ZXLB(IL)=ZXLB(IL)-ZFRL(IL)
       ZXIB(IL,JK)=ZXIB(IL,JK)+ZFRL(IL)
       ZFRL(IL)=ZCLF(IL,JK)*ZFRL(IL)
       ZXLB(IL)=MAX(ZXLB(IL),ZERO)
       IF(ISAVLS.NE.0) THEN
         FRSTEMP=1.4E-20*ZF1
         XSUM=ZCLF(IL,JK)*(FRKTEMP+FRSTEMP)*ZTMST
         IF ( XSUM > 0. ) THEN
           FRKROL(IL,JK)=FRKTEMP*(ZFRL(IL)/XSUM)*ZCLF(IL,JK)
           FRSROL(IL,JK)=FRSTEMP*(ZFRL(IL)/XSUM)*ZCLF(IL,JK)
         ENDIF
       ENDIF
      ENDIF
C
C     * NEW TEMPERATURE AND MOISTURE.
C
      ZLVDCP=HTVOCP(TH(IL,JK),ONE)
      ZLSDCP=HTVOCP(TH(IL,JK),ZERO)
      ZLFDCP=ZLSDCP-ZLVDCP
      TH(IL,JK)=TH(IL,JK)+ZLFDCP*ZFRL(IL)
  410 CONTINUE
C
C*       4.  CLOUD PHYSICS AND PRECIPITATION FLUXES.
C            ----- ------- --- ------------- ------
C
      DO 425 IL=IL1,IL2
C
C      * 4.1 WARM CLOUD PHYSICS (AFTER KHAIROUTDINOV AND KOGAN, 2000).
C            
       ZRHO0=P(IL,JK)/(RGAS*TH(IL,JK))
       ZQRHO0=1.3/ZRHO0
C
C      * RAIN WATER MIXING RATIO FROM RAINFALL RATE (ROTSTAYN, 1997)
C
       ZXRP1=0.
       VTR(IL)=0.
       AKR=141.4
       IF ( ZRNFT(IL).GT.ZEPCLC ) THEN
         ZRFLNL=ZRFLN(IL)/ZRNFT(IL)
         ZRFLNL=MIN(ZRFLNL,ZRFLNL_MAX)
         IF ( ZRFLNL.GT.1.E-20 ) THEN
           ALAMBR=714.*(1.2/ZRHO0)**(1./9.)*ZRFLNL**(-0.22)
           VTR(IL)=1.94*AKR*SQRT(1.2/(ZRHO0*ALAMBR))
           ZXRP1=ZRFLNL/(ZRHO0*VTR(IL))
         ENDIF
       ENDIF

       IF (ZXLB(IL)*ZCLF(IL,JK).GE.ZSEC)            THEN
         ZSEC2=ZSEC
         CDDFIX=50.
         IF ( I2NDIE .EQ. 1 ) THEN
           RVOL=1.E+06*(3.*ZXLB(IL)/(4.*API*RHOH2O*ZQNL(IL)))**(1./3.)
           BETA6=((RVOL+3.)/RVOL)**(1./3.)
           R6=BETA6*RVOL
           R6C=7.5/(ZXLB(IL)**(1./6.)*SQRT(R6))
           FACAUTN=1.08*FACAUT
           EP=FACAUTN*1.E+10*(BETA6**6)
           IF ( R6 >= R6C ) THEN
             FAAUT=ZCRAUT(IL,JK)*EP/ZQNL(IL)
           ELSE
             FAAUT=0.
           ENDIF
         ELSE
           FAAUT=FACAUT*ZCRAUT(IL,JK)*1350.
           FAAUT=FAAUT*(CDDFIX)**(-1.79)
         ENDIF
C
C        *** INITIAL GUESS.
C
         ZQEST1=(ZXLB(IL)**(1.-COAUT)
     1         -ZTMST*FAAUT*(1.-COAUT))**(1./(1.-COAUT))
         ZQEST=MIN(MAX(ZQEST1,ZSEC2),ZXLB(IL))
         FAACC=0.
         IF ( ZXRP1.GT.ZSEC2 ) THEN
           FAACC=FACACC*67.*ZXRP1**COACC
           ZQEST2=(ZXLB(IL)**(1.-COACC)
     1           -ZTMST*FAACC*(1.-COACC))**(1./(1.-COACC))
           ZQEST=MAX(MIN(ZQEST,ZQEST2),ZSEC2)
           APK=-LOG(ZQEST/ZXLB(IL))/ZTMST
           IF ( ZQEST.LT.ZXLB(IL) .AND. APK.GT.0. ) THEN
             ALPHA=((1.-ZQEST/ZXLB(IL))*ZXLB(IL)/(APK*ZTMST)-ZQEST)
     1            /(ZXLB(IL)-ZQEST)
             ALPHA=MIN(MAX(ALPHA,ZERO),ONE)
C
C            *** FIRST ITERATION.
C
             AQHAT =ALPHA*ZXLB(IL)+(1.-ALPHA)*ZQEST
             AQHAT1=AQHAT**COTMP1
             AQHAT2=AQHAT**COTMP2
             ATERM1=FAAUT*AQHAT1*ZXLB(IL)+FAACC*AQHAT2
             ATERM2=FAAUT*AQHAT1*ZQEST   +FAACC*AQHAT2
             AEXPF =EXP(-FAACC*AQHAT2*ZTMST)
             AFUN=ATERM1*ZQEST-AEXPF*ATERM2*ZXLB(IL)
             DQHDQ=1.-ALPHA
             ATERM3=FAAUT*COTMP1*AQHAT1*DQHDQ/AQHAT
             ATERM4=FAACC*COTMP2*AQHAT2*DQHDQ/AQHAT
             ADFUN=ATERM1+ZQEST*(ATERM3*ZXLB(IL)+ATERM4)
     1            +AEXPF*ZXLB(IL)*(ZTMST*ATERM4*ATERM2
     2                          -ATERM3*ZQEST+FAAUT*AQHAT1+ATERM4)
             ZQEST=MAX(ZQEST-AFUN/ADFUN,ZSEC2)
C
C            *** SECOND ITERATION.
C
             AQHAT =ALPHA*ZXLB(IL)+(1.-ALPHA)*ZQEST
             AQHAT1=AQHAT**COTMP1
             AQHAT2=AQHAT**COTMP2
             ATERM1=FAAUT*AQHAT1*ZXLB(IL)+FAACC*AQHAT2
             ATERM2=FAAUT*AQHAT1*ZQEST   +FAACC*AQHAT2
             AEXPF =EXP(-FAACC*AQHAT2*ZTMST)
             AFUN=ATERM1*ZQEST-AEXPF*ATERM2*ZXLB(IL)
             DQHDQ=1.-ALPHA
             ATERM3=FAAUT*COTMP1*AQHAT1*DQHDQ/AQHAT
             ATERM4=FAACC*COTMP2*AQHAT2*DQHDQ/AQHAT
             ADFUN=ATERM1+ZQEST*(ATERM3*ZXLB(IL)+ATERM4)
     1            +AEXPF*ZXLB(IL)*(ZTMST*ATERM4*ATERM2
     2                          -ATERM3*ZQEST+FAAUT*AQHAT1+ATERM4)
             ZQEST=MAX(ZQEST-AFUN/ADFUN,ZSEC2)
C
C            *** THIRD ITERATION.
C
             AQHAT =ALPHA*ZXLB(IL)+(1.-ALPHA)*ZQEST
             AQHAT1=AQHAT**COTMP1
             AQHAT2=AQHAT**COTMP2
             ATERM1=FAAUT*AQHAT1*ZXLB(IL)+FAACC*AQHAT2
             ATERM2=FAAUT*AQHAT1*ZQEST   +FAACC*AQHAT2
             AEXPF =EXP(-FAACC*AQHAT2*ZTMST)
             AFUN=ATERM1*ZQEST-AEXPF*ATERM2*ZXLB(IL)
             DQHDQ=1.-ALPHA
             ATERM3=FAAUT*COTMP1*AQHAT1*DQHDQ/AQHAT
             ATERM4=FAACC*COTMP2*AQHAT2*DQHDQ/AQHAT
             ADFUN=ATERM1+ZQEST*(ATERM3*ZXLB(IL)+ATERM4)
     1            +AEXPF*ZXLB(IL)*(ZTMST*ATERM4*ATERM2
     2                          -ATERM3*ZQEST+FAAUT*AQHAT1+ATERM4)
             ZQEST=MAX(ZQEST-AFUN/ADFUN,ZSEC2)
           ENDIF
         ENDIF
         ZRPR(IL)=MIN(ZCLF(IL,JK)*ZXLB(IL),
     1                MAX(-ZCLF(IL,JK)*(ZQEST-ZXLB(IL)),ZERO))
         IF(ISAVLS.NE.0) THEN
           APK=-LOG(ZQEST/ZXLB(IL))/ZTMST
           IF ( APK > 0. ) THEN
             FA=FAAUT*(ZXLB(IL)**COAUT)*(1.-EXP(-APK*COAUT*ZTMST))
     1         /(APK*COAUT*ZTMST)
             FB=FAACC*(ZXLB(IL)**COACC)*(1.-EXP(-APK*COACC*ZTMST))
     1         /(APK*COACC*ZTMST)
             IF ( (FA+FB) /= 0. ) RAT(IL)=FA/(FA+FB)
           ENDIF
         ENDIF
       ELSE
         ZRPR(IL)=ZCLF(IL,JK)*ZXLB(IL)
       ENDIF
       IF(ZCLF(IL,JK).GT.ZEPCLC) THEN
         ZMRATEP(IL,JK)=ZRPR(IL)/ZCLF(IL,JK)
         ZMLWC(IL,JK)=ZXLB(IL)
       ENDIF
C
C      * CORRECT, IF SUM OF PROCESSES SO FAR RESULTS IN NEGATIVE CLOUD
C      * LIQUID WATER CONTENT.
C
       ZXLTE=-ZRPR(IL)+ZIMLT(IL)-ZFRL(IL)
       IF ( XROW(IL,JK+1,ILWC)+ZXLTE .LT. 0. .AND. ZRPR(IL).GT. 0.
     1                                                           ) THEN
          ACOR=(ZIMLT(IL)-ZFRL(IL)+XROW(IL,JK+1,ILWC))/ZRPR(IL)
          ZRPR(IL)=MIN(MAX(ACOR,ZERO),ONE)*ZRPR(IL)
       ENDIF
       IF(ISAVLS.NE.0) THEN
         RACLROL(IL,JK)=ZRPR(IL)*ZQTMST*(1.-RAT(IL))
         AUTROL (IL,JK)=ZRPR(IL)*ZQTMST*RAT(IL)
       ENDIF
C
C        * 4.2  COLD CLOUD PHYSICS
C 
C        AGGREGATION OF CLOUD ICE TO SNOW.
C 
       IF (ZXIB(IL,JK)*ZCLF(IL,JK).GT.ZSEC)          THEN
         WCFAC=MAX(ZXIB(IL,JK)*ZRHO0*1000.,WCMIN)
         ZRIEFF0=83.8*WCFAC**0.216
         ZRIH0=-2261.236+SQRT(5113188.044+2809.*ZRIEFF0**3)
         ZRI=MAX(ZRIMIN,1.E-6*ZRIH0**(1./3.))
         IF (ZRI.LE.0.9999E-4)             THEN
           ZC1=17.5*ZXIB(IL,JK)*ZRHO0/ZRHOI*ZQRHO0**0.33
           ZDT2=-6./ZC1*LOG10(ZRI*1.E4)
           ZSAUT(IL)=ZCSAUT*ZTMST/ZDT2
         ELSE
C
C          * FOR SUITABLY LARGE ICE CRYSTALS, THE PROCESS
C          * IS EFFECTIVE ENOUGH TO CONVERT EVERYTHING.
C          * SO, WE ESSENTIALLY SET ZSAUT TO A LARGE
C          * ENOUGH NUMBER SO THAT THE BOUNDED EXPRESSION
C          * BELOW REMOVES ALL OF THE ICE.
C 
           ZSAUT(IL)=ZTMST
         ENDIF
         ZSAUT(IL)=ZXIB(IL,JK)*(1.-EXP(-ZSAUT(IL)))
C
C        * ACCRETION OF SNOW BY ICE CRYSTALS (BASED ON ROTSTAYN, 1997)
C
         IF (ZSFLN(IL).GE.ZSEC) THEN
           ZCOLEFF=EXP(0.025*(TH(IL,JK)-T1S))  
           ZLAMS=1.6E+03*10.**(0.023*(T1S-TH(IL,JK)))
           ZSACI(IL)=-ZCOLEFF*ZLAMS*ZSFLN(IL)*ZTMST
     1                                 /(2.*RHOSNO*ZCLF(IL,JK))
           ZSACI(IL)=ZXIB(IL,JK)*(1.-EXP(ZSACI(IL)))
         ELSE
           ZSACI(IL)=0.
         ENDIF
C
         ZSPR(IL)=ZCLF(IL,JK)*MAX(ZERO,MIN(ZSACI(IL)+ZSAUT(IL),
     1            ZXIB(IL,JK)))
       ELSE
         ZSPR(IL)=ZCLF(IL,JK)*ZXIB(IL,JK)
       ENDIF
C
C      * CORRECT, IF SUM OF PROCESSES SO FAR RESULTS IN NEGATIVE CLOUD
C      * ICE WATER CONTENT.
C
       ZXITE=-ZSPR(IL)-ZIMLT(IL)+ZFRL(IL)
       IF ( XROW(IL,JK+1,IIWC)+ZXITE .LT. 0. .AND. ZSPR(IL).GT. 0.
     1                                                           ) THEN
          ACOR=(-ZIMLT(IL)+ZFRL(IL)+XROW(IL,JK+1,IIWC))/ZSPR(IL)
          ZSPR(IL)=MIN(MAX(ACOR,ZERO),ONE)*ZSPR(IL)
       ENDIF
C
       IF(ISAVLS.NE.0) THEN
C
C        * WEIGHT TERMS TO GET TRUE PICTURE.
C
         XSUM=ZCLF(IL,JK)*(ZSACI(IL)+ZSAUT(IL))
         IF(XSUM  > 0.)                                    THEN
           AGGROL (IL,JK)=ZCLF(IL,JK)*ZSAUT(IL)*(ZSPR(IL)/XSUM)*ZQTMST
           SACIROL(IL,JK)=ZCLF(IL,JK)*ZSACI(IL)*(ZSPR(IL)/XSUM)*ZQTMST
         ENDIF
       ENDIF
C
C      * ACCRETION ON OF SNOW BY CLOUD DROPLETS (ROTSTAYN, 1997)
C
       IF (ZXLB(IL)*ZCLF(IL,JK).GT.ZSEC.AND.ZSFLN(IL).GE.ZSEC) THEN
         ZLAMS=1.6E+03*10.**(0.023*(T1S-TH(IL,JK)))
         ZSACL(IL)=-ZEAC*ZLAMS*ZSFLN(IL)*ZTMST
     1                                 /(2.*RHOSNO*ZCLF(IL,JK))
         ZSACL(IL)=ZXLB(IL)*(1.-EXP(ZSACL(IL)))
         ZSACL(IL)=MAX(ZERO,MIN(ZCLF(IL,JK)*ZSACL(IL),
     1             (ZCLF(IL,JK)*ZXLB(IL)-ZRPR(IL))))
C
C        * CORRECT, IF SUM OF PROCESSES SO FAR RESULTS IN NEGATIVE CLOUD
C        * LIQUID WATER CONTENT.
C
         ZXLTE=-ZRPR(IL)+ZIMLT(IL)-ZFRL(IL)-ZSACL (IL)
         IF ( XROW(IL,JK+1,ILWC)+ZXLTE .LT. 0. .AND. ZSACL(IL).GT. 0.
     1                                                           ) THEN
            ACOR=(-ZRPR(IL)+ZIMLT(IL)-ZFRL(IL)+XROW(IL,JK+1,ILWC))
     1          /ZSACL(IL)
            ZSACL(IL)=MIN(MAX(ACOR,ZERO),ONE)*ZSACL(IL)
         ENDIF
         IF(ISAVLS.NE.0) THEN
           SACLROL(IL,JK)=ZSACL(IL)*ZQTMST
         ENDIF
       ENDIF
 425  CONTINUE
C
      DO 500 IL=IL1,IL2
        ZEVP(IL)=0.
        ZSUB(IL)=0.
C
C*      5.1   NEW PRECIPITATION FLUXES AFTER CONSIDERATION OF
C*            COALESCENCE AND SEDIMENTATION PROCESSES.
C
        ZZDRR=MAX(ZERO,ZCONS2*ZDP(IL)*ZRPR(IL))
        ZZDRS=MAX(ZERO,ZCONS2*ZDP(IL)*(ZSPR(IL)+ZSACL(IL)))
        ZRFLN(IL)=ZRFLN(IL)+ZZDRR
        ZSFLN(IL)=ZSFLN(IL)+ZZDRS
  500 CONTINUE
C
C*    5.2     IF THERE IS NO PRECIPITATION AVOID EVAPORATION.
C
      DO 510 IL=IL1,IL2
        ZSR=ZSR+ZRFLN(IL)
 510  CONTINUE
      IF (ZSR.GT.0.)                                             THEN
C
C*     5.3     EVAPORATION OF RAIN
C*             ----------- -- ----
C
       DO 530 IL=IL1,IL2
        IF ( ZRFLN(IL).GT.0. .AND. ZRCLR(IL).GT.0. ) THEN
          ZRFLNL=ZRFLN(IL)/ZRNFT(IL)
          ZRFLNL=MIN(ZRFLNL,ZRFLNL_MAX)
          ZRHO0=P(IL,JK)/(RGAS*TH(IL,JK))
          RHODZ=ZDP(IL)/GRAV
C
          ETMP=ESW(TH(IL,JK))
          ESTREF=0.01*P(IL,JK)*(1.-EPSLIM)/(1.-EPSLIM*EPS2)
          IF ( ETMP.LT.ESTREF ) THEN
            ESAT=ETMP
          ELSE
            ESAT=ESTREF
          ENDIF
C
          QSW=EPS1*ESAT/(0.01*P(IL,JK)-EPS2*ESAT)
          RHCW=MIN(MAX((QC(IL,JK)*(1.+QSW*FACTE))
     1                /(QSW*(1.+QC(IL,JK)*FACTE)),0.),1.)
          SUSATW=RHCW-1.
          ZXSEC=1.-1.E-12
          ZKA=0.024
          ZDV=2.21/P(IL,JK)
          ZLV=ZLVDCP*CP
          ZAST=ZLV*(ZLV/(RGASV*TH(IL,JK))-1.)/(ZKA*TH(IL,JK))
          ZBST=RGASV*TH(IL,JK)/(ZDV*ESAT)
          ZZEPR1=870.*SUSATW*ZRFLNL**0.61/(SQRT(ZRHO0)*(ZAST+ZBST))
          ZZEPR=ZZEPR1*RHODZ
          ZPLIM=-ZXSEC*ZRFLN(IL)/ZRNFT(IL)
          ITOTEVP=0
          IF ( ZZEPR .LT. ZPLIM ) THEN
            ZZEPR=ZPLIM
            ITOTEVP=1
          ENDIF
          ZEVP(IL)=-ZZEPR*ZTMST*ZRCLR(IL)/RHODZ
          ZEVP(IL)=MIN(ZEVP(IL),
     1           MAX(ZXSEC*((1.-RHCW)*QC(IL,JK))*ZRCLR(IL),ZERO))
          ZEVP(IL)=MAX(ZEVP(IL),ZERO)
          ZZEPR=-ZEVP(IL)*RHODZ/ZTMST
          ZFRAIN(IL,JK)=ZFRAIN(IL,JK)+ZRFLN(IL)
          CLRFR(IL,JK)=ZRCLR(IL)
          IF ( ITOTEVP .EQ. 1 ) THEN
            ZFEVAP(IL,JK)=MIN(MAX(0.,-ZZEPR/ZRFLN(IL)),1.)
          ENDIF
          ZRFLN(IL)=ZRFLN(IL)+ZZEPR
          IF(ISAVLS.NE.0) THEN
             EVPROL (IL,JK)=ZEVP(IL)*ZQTMST
             RAINROL(IL,JK)=ZRFLN(IL)-RAINROL(IL,JK-1)
          ENDIF
        ENDIF
  530  CONTINUE
      ENDIF
C
      DO 535 IL=IL1,IL2
        ZSS=ZSS+ZSFLN(IL)
 535  CONTINUE
      IF (ZSS.GT.0.)                                             THEN
C
C*        5.4     SUBLIMATION OF SNOW
C*                ----------- -- ----
C
       DO 540 IL=IL1,IL2
         IF(TH(IL,JK).LE.ZMELTP2 .AND. ZSFLN(IL).GT.ZSEC
     1                                     .AND. ZSCLR(IL).GT.0.)  THEN
           ZRHO0=P(IL,JK)/(RGAS*TH(IL,JK))
           RHODZ=ZDP(IL)/GRAV
           IF ( ZSNFT(IL).GT.ZEPCLC ) THEN
             ZSFLNL=ZSFLN(IL)/ZSNFT(IL)
             ZXSP1=(1./ZRHO0)*(ZSFLNL/3.23)**(1./1.17)
           ELSE
             ZXSP1=0.
           ENDIF
           ETMP=ESI(TH(IL,JK))
           ESTREF=0.01*P(IL,JK)*(1.-EPSLIM)/(1.-EPSLIM*EPS2)
           IF ( ETMP.LT.ESTREF ) THEN
             ESAT=ETMP
           ELSE
             ESAT=ESTREF
           ENDIF
           QSI=EPS1*ESAT/(0.01*P(IL,JK)-EPS2*ESAT)

           RHCI=MIN(MAX((QC(IL,JK)*(1.+QSI*FACTE))
     1                 /(QSI*(1.+QC(IL,JK)*FACTE)),0.),1.)
           SUSATI=RHCI-1.
           CLAMBS=(ZRHO0*ZXSP1/(API*RHOSNO*ZN0S))**0.25
           ZLSDCP=HTVOCP(TH(IL,JK),ZERO)
           CFAC4C     = 0.78 * CLAMBS**2.
     1                + FACT * (1.3/ZRHO0)**.25 * CLAMBS**2.625
           AAAA       = 2.*API*SUSATI
           BBBB       = ZLSDCP**2. / (THERCO*RGASV*(TH(IL,JK)**2.))
           CCCC       = 1. / (ZRHO0 * QSI * PSI)
           COEFF      = ZNZEROS * AAAA/ (ZRHO0 * (BBBB + CCCC))
           ZXSEC=1.-1.E-12
           ZZEPS=COEFF*CFAC4C*RHODZ
           ZPLIM=-ZXSEC*ZSFLN(IL)/ZSNFT(IL)
           ITOTSUB=0
           IF ( ZZEPS .LT. ZPLIM ) THEN
             ZZEPS=ZPLIM
             ITOTSUB=1
           ENDIF
           ZSUB(IL)=-ZZEPS*ZTMST*ZSCLR(IL)/RHODZ
           ZSUB(IL)=MIN(ZSUB(IL),
     1           MAX(ZXSEC*((1.-RHCI)*QC(IL,JK))*ZSCLR(IL),ZERO))
           ZSUB(IL)=MAX(ZSUB(IL),ZERO)
           ZZEPS=-ZSUB(IL)*RHODZ/ZTMST
           ZFSNOW(IL,JK)=ZFSNOW(IL,JK)+ZSFLN(IL)
           CLRFS(IL,JK)=ZSCLR(IL)
           IF ( ITOTSUB .EQ. 1 ) THEN
             ZFSUBL(IL,JK)=MIN(MAX(0.,-ZZEPS/ZSFLN(IL)),1.)
           ENDIF
           ZSFLN(IL)=ZSFLN(IL)+ZZEPS
           IF(ISAVLS.NE.0) THEN
             SUBROL (IL,JK)=ZSUB(IL)*ZQTMST
             SNOWROL(IL,JK)=ZSFLN(IL)-SNOWROL(IL,JK-1)
           ENDIF
         ENDIF
  540  CONTINUE
      ENDIF
C
C*      6.    INCREMENTATION OF T,Q AND X TENDENCIES AND SURFACE FLUXES
C             -------------- -- - - --- - ---------- --- ------- ------
C
      DO 600 IL=IL1,IL2
C
C*      6.1     MODIFICATION OF THE T,Q AND X TENDENCIES,
C*                 CLOUD COVER FOR DIAGNOSTICS AND RADIATION.
C
        ZLVDCP=HTVOCP(TH(IL,JK),ONE)
        ZLSDCP=HTVOCP(TH(IL,JK),ZERO)
        ZQTE=ZEVP(IL)+ZSUB(IL)
        ZTTE=-ZLVDCP*ZEVP(IL)-ZLSDCP*ZSUB(IL)
     1       +(ZLSDCP-ZLVDCP)*ZSACL(IL)
C
        ZXLTE=-ZRPR(IL)+ZIMLT(IL)-ZFRL(IL)-ZSACL (IL)
        ZXITE=-ZSPR(IL)-ZIMLT(IL)+ZFRL(IL)
C
        ZXITE=ZXITE + ZFRAC(IL,JK)     * DQLDT(IL,JK)
        ZXLTE=ZXLTE + (1.-ZFRAC(IL,JK))* DQLDT(IL,JK)
C
        TH(IL,JK)=TH(IL,JK)+ZTTE
        QV(IL,JK)=QV(IL,JK)+ZQTE
        XROW(IL,JK+1,ILWC)=XROW(IL,JK+1,ILWC)+ZXLTE
        XROW(IL,JK+1,IIWC)=XROW(IL,JK+1,IIWC)+ZXITE
C
C       * CORRECTION: AVOID SMALL OR NEGATIVE CLOUD WATER AND ICE.
C
        ZXLOLD=XROW(IL,JK+1,ILWC)
        LO=(XROW(IL,JK+1,ILWC).LT.ZSEC)
        IF(LO) XROW(IL,JK+1,ILWC)=0.
        ZDXLCOR=(XROW(IL,JK+1,ILWC)-ZXLOLD)
C 
        ZXIOLD=XROW(IL,JK+1,IIWC)
        LO=(XROW(IL,JK+1,IIWC).LT.ZSEC)
        IF(LO) XROW(IL,JK+1,IIWC)=0.
        ZDXICOR=(XROW(IL,JK+1,IIWC)-ZXIOLD)
C
        QV(IL,JK)=QV(IL,JK)-ZDXLCOR-ZDXICOR
        TH(IL,JK)=TH(IL,JK)+ZDXLCOR*ZLVDCP+ZDXICOR*ZLSDCP
        IF ( ZCLF(IL,JK).LT.ZSEC ) ALMC(IL,JK)=0.
  600 CONTINUE
C
C*         6.4     SWAP OF THE FLUXES AND END OF THE VERTICAL LOOP.
C
      DO 641 IL=IL1,IL2
        ZRFL(IL)=ZRFLN(IL)
        ZSFL(IL)=ZSFLN(IL)
  641 CONTINUE
C
C     * DEFINE CLOUD WATER AND ICE CONTENTS TO BE SUBSEQUENTLY USED
C     * IN CLOUDS.
C
      DO 650 IL=IL1,IL2
        IF(ZCLF(IL,JK).GT.0.)                        THEN
          QLWC(IL,JK)=XROW(IL,JK+1,ILWC)
          QIWC(IL,JK)=XROW(IL,JK+1,IIWC)
        ELSE
          QLWC(IL,JK)=0.
          QIWC(IL,JK)=0.
        ENDIF
  650 CONTINUE
C
C     * EXTRA COSP FIELDS.
C
      IF (ICOSP .EQ. 1) THEN
        DO IL = IL1,IL2
C
C         * Rain parameters.
C
          ZRHO0=P(IL,JK)/(RGAS*TH(IL,JK))
          ZXRP1=0.
          VTR(IL)=0.
          AKR=141.4
          RMIXROL(IL,JK)=0.0
          RREFROL(IL,JK) = 0.0
          IF ( ZRNFT(IL).GT.ZEPCLC ) THEN
            ZRFLNL=ZRFLN(IL)/ZRNFT(IL)
            ZRFLNL=MIN(ZRFLNL,ZRFLNL_MAX)
            IF ( ZRFLNL.GT.1.E-20) THEN
              ALAMBR  = 714.*(1.2/ZRHO0)**(1./9.)*ZRFLNL**(-0.22)
              VTR(IL) = 1.94*AKR*SQRT(1.2/(ZRHO0*ALAMBR))
              ZXRP1   = ZRFLNL/(ZRHO0*VTR(IL))
              IF (ZXRP1 .GE. ZSEC) THEN ! Do I need this check?
                RMIXROL(IL,JK) = ZXRP1*ZRNFT(IL)
                RREFROL(IL,JK) = (3.0/2.0)/ALAMBR
              ENDIF
            ENDIF
          ENDIF
C            
C         * Snow parameters.
C
          SMIXROL(IL,JK)=0.0
          SREFROL(IL,JK)=0.0
          IF ( ZSNFT(IL).GT.ZEPCLC ) THEN
            ZSFLNL=ZSFLN(IL)/ZSNFT(IL)
            ZXSP1=(1./ZRHO0)*(ZSFLNL/3.23)**(1./1.17)
            IF (ZXSP1 .GE. ZSEC) THEN ! Do I need this check?
              SMIXROL(IL,JK) = ZXSP1*ZSNFT(IL)
C
C             * Use Rotstayn, 1997 approach to set slope and
C             * intercept parameter (Eq. 27 and 28).
C
              ZLAMS = 1.6E+03*10.**(0.023*(T1S-TH(IL,JK)))
              ZN0F  = (ZRHO0*ZXSP1/API*RHOSNO)*ZLAMS**4
C
C             * Compute snow particle size, noting there is 
C             * quite abit of uncertainity what a mean diameter
C             * means for a snow particle.
C
              SREFROL(IL,JK) = (3.0/2.0)/ZLAMS
            ENDIF
          ELSE
            ZXSP1=0.
          ENDIF
        END DO ! IL
      ENDIF

  700 CONTINUE
C                                                                    
      RETURN
      END
