      SUBROUTINE STATCLD5(QCW,ZCLF,SIGMA,ZCRAUT,QCWVAR,SSH,CVSG,
     1                    QT,HMN,ZFRAC,CPM,P,Z,RRL,CONS,CRH,
     2                    ALMIX,DHLDZ,DRWDZ,I2NDIE,
     3                    GRAV,ZTMST,ILEV,ILG,IL1,IL2,ICVSG,ISUBG,
     4                    LVL,ICALL                              )
C
C     * AUG 16/2013 - K.VONSALZEN. Cosmetic change to abort condition
C     *                            where now is QCW.GT.QT instead of
C     *                            QCW.GE.QT. This does not change
C     *                            the answer but prevents infrequent
C     *                            crashes with small downward
C     *                            moisture fluxes causing very small
C     *                            negative total water.
C     * JUL 31/2013 - M.LAZARE.    Cosmetic change to pass in LVL
C     *                            and ICALL to aid in future debugging.
C     * JUN 25/2013 - K.VONSALZEN. New version for gcm17:
C     *                            - Second indirect effect parameter
C     *                              I2NDIE passed in and used.
C     *                            - Function ZCR replaces look-up
C     *                              table ZCR to calculate autoconversion
C     *                              efficiency.
C     * APR 29/2012 - K.VONSALZEN. PREVIOUS VERSION STATCLD4 FOR GCM16:
C     *                            ELIMINATE CLOUDS AND CONDENSATE IF
C     *                            CLOUD FRACTION OR CONDENSATE MIXING
C     *                            RATIO TOO LOW TO ENSURE CONSISTENT
C     *                            RESULTS FOR CLOUD FRACTION AND CLOUD
C     *                            CONDENSATE.
C     * MAR 25/2009 - L.SOLHEIM.   REVISED COSMETIC CHANGE:
C     *                            - USE REAL*4 OR REAL*8 ARGUMENT TO
C     *                              "ERF", DEPENDING ON THE VALUE OF
C     *                              "MACHINE", SO WILL WORK SEAMLESSLY ON
C     *                              PGI COMPILER AS WELL
C     * DEC 19/2007 - K.VONSALZEN/ PREVIOUS VERSION STATCLD3 FOR GCM15G/H/I:
C     *               M.LAZARE.    - ELIMINATE UNPHYSICAL RESULTS
C     *                              BY NOT PERMITTING THE INFERRED
C     *                              WATER VAPOUR MIXING RATIO (QCW)
C     *                              TO EXCEED A THRESHOLD (DEFINED
C     *                              BY QCWMAX).
C     *                            - PASSES IN ADELT FROM PHYSICS
C     *                              AS ZTMST INSTEAD OF ZTMST=2.*DELT.
C     * JAN 13/2007 - K.VONSALZEN. CALCULATE NORMALIZED CLOUD WATER
C     *                            VARIANCE TO BE USED IN CLOUD
C     *                            INHOMOGENEITY (QCWVAR).
C     * JUN 19/2006 - M.LAZARE.  NEW VERSION FOR GCM15F:
C     *                          - COSMETIC: 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     * MAY 06/2006 - M.LAZARE/  PREVIOUS VERSION STATCLD FOR GCM15E.
C     *               K.VONSALZEN:

C
C     * THIS SUBROUTINE PERFORMS THE BULK OF THE STATISTICAL CLOUD
C     * SCHEME BASED ON CHABOUREAU AND BECHTHOLD (C&B). IT TAKES AS
C     * INPUT LIQUID STATIC ENERGY AND TOTAL WATER, AND OUTPUTS
C     * CONDENSED WATER/ICE, CLOUD FRACTION, THE VARIANCE OF MELLORS
C     * VARIABLE ("S"), THE EFFICIENCY FACTOR FOR AUTOCONVERSION
C     * FROM CLOUD WATER TO RAIN WATER, AND THE SATURATION SPECIFIC
C     * HUMIDITY.
C
C********************************************************************
C     * DICTIONARY OF VARIABLES:
C
C     * OUTPUT:
C     * ------
C
C     * QCW:     TOTAL CONDENSED WATER/ICE (BOTH PHASES) IN KG/KG.
C     * ZCLF:    CLOUD FRACTION (DIMENSIONLESS)
C     * SIGMA:   VARIANCE OF MELLOR'S VARIABLE S IN KG/KG
C     * ZCRAUT:  EFFICIENCY OF AUTOCONVERSION FROM CLOUD WATER TO
C     *          RAIN WATER.
C     * QCWVAR:  NORMALIZED CLOUD WATER VARIANCE (<W'**2>/<W>**2)
C     *          (DIMENSIONLESS).
C     * SSH:     SATURATION SPECIFIC HUMIDITY IN KG/KG
C     * CVSG:    CONTRIBUTION TO TOTAL VARIANCE OF S FROM CONVECTIVE
C     *          PROCESSES IN KG/KG (ONLY IF SWITCH ICVSG=1,
C     *          OTHERWISE IT IS INPUT!).
C
C     * INPUT:
C     * -----
C
C     * QT:     TOTAL WATER (BOTH PHASES) IN KG/KG.
C     * HMN:    TOTAL LIQUID STATIC ENERGY IN JOULES/KG.
C     * ZFRAC:  FRACTION OF ICE PHASE.
C     * CPM:    SPECIFIC HEAT IN JOULES/(KG-DEGK).
C     * P:      MID-LAYER PRESSURE IN MBS.
C     * Z:      MID-LAYER HEIGHT ABOVE GROUND IN METRES.
C     * PBLT:   LAYER INDEX OF TOP OF PBL.
C     * RRL:    WATER/ICE WEIGHTED LATENT HEAT IN JOULES/KG.
C     * CONS:   FIELD CONTRIBUTING TO VARIANCE FROM CONVECTION.
C     * CRH:    FLAG (0/1) TO CONTROL CALCULATION.
C     * ALMIX:  MIXING LENGTH FOR C&B PARAMETERIZATION IN METRES.
C     * DHLDZ:  TOTAL WATER VERTICAL GRADIENT IN (M-1).
C     * DRWDZ:  LIQUID STATIC ENERGY GRADIENT IN JOULES/(KG-M). 
C********************************************************************
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      REAL  , DIMENSION(ILG), INTENT(OUT) :: QCW,ZCLF,SIGMA,ZCRAUT,
     1                                       QCWVAR,SSH,CVSG
      REAL  , DIMENSION(ILG), INTENT(IN)  :: HMN,QT,P,Z,CPM,ZFRAC,
     1                                       CONS,CRH,RRL,ALMIX,DHLDZ,
     2                                       DRWDZ
C
      REAL*4 :: R4VAL
      REAL*8 :: R8VAL
C
      COMMON /EPS  / AP,BP,EPS1,EPS2
      COMMON /ESTWI/ RW1,RW2,RW3,RI1,RI2,RI3
C
C     * COMMON BLOCK TO HOLD WORD SIZE.
C
      INTEGER MACHINE,INTSIZE
      COMMON /MACHTYP/ MACHINE,INTSIZE
C
      ESW(TTT)    = EXP(RW1+RW2/TTT)*TTT**RW3
      ESI(TTT)    = EXP(RI1+RI2/TTT)*TTT**RI3
C
      DATA ZERO,ONE /0., 1./
C==============================================================================================
C     * PARAMETERS.
C
      ZEPCLC=1.E-2
      ZSEC=1.E-9
      CSIGMA=0.2
      YFR=0.50
      SQRT2=SQRT(2.)
      PI=3.141592653589793
      SQRT2PI=SQRT(2.*PI)
      Q1MIN1=-3.
      Q1MIN2=TAN((ZEPCLC-0.5)/0.36)/1.55
      Q1MIN=MAX(Q1MIN1,Q1MIN2)
      Q1MINM=Q1MIN/2.
      Q1MAX=2.
      Q1MAXM=Q1MAX/2.
      EXPM1=EXP(-1.)
      QDELT=2./30.
      TAUSIG=21600.
      SIGFAC=EXP(-ZTMST/TAUSIG)
C
C     * INITIALIZATION.
C
      ILBAD=0
      DO IL=IL1,IL2
        QCW   (IL)=0.
        ZCLF  (IL)=0.
        ZCRAUT(IL)=1.
        SIGMA (IL)=0.
        SSH   (IL)=0.
      ENDDO
C
      DO IL=IL1,IL2
        IF ( CRH(IL).NE.0. ) THEN
C
C         * TEMPERATURE UNDER CLEAR-SKY CONDITIONS AS FIRST GUESS.
C
          TX=(HMN(IL)-GRAV*Z(IL))/CPM(IL)
          TL=TX
C
C         * WATER VAPOUR SATURATION MIXING RATIO UNDER CLEAR-SKY CONDITIONS.
C
          EST=(1.-ZFRAC(IL))*ESW(TX)+ZFRAC(IL)*ESI(TX)
          RSS=EPS1*EST/( P(IL)-EPS2*EST )
          DRSSDT=EPS1*P(IL)
     1             *((1.-ZFRAC(IL))*ESW(TX)*(RW3-RW2/TX)
     2                  +ZFRAC(IL) *ESI(TX)*(RI3-RI2/TX))
     3              /(TX*(P(IL)-EPS2*EST)**2)
          RSSTL=RSS
          DRSSDTL=DRSSDT
C
C         * ACCOUNT FOR CLOUD WATER/ICE IF SATURATION. USE
C         * SERIES REVERSION AND ITERATE.
C
          IF ( RSS.LT.QT(IL) ) THEN
C
C           * FIRST ITERATION.
C
            HMNPR=CPM(IL)*TX+GRAV*Z(IL)-HMN(IL)
            DHMNPR=CPM(IL)
            RC=QT(IL)-RSS
            HMNPR=HMNPR-RRL(IL)*RC
            DHMNPR=DHMNPR+RRL(IL)*DRSSDT
            TX=TX-HMNPR/DHMNPR
            EST=(1.-ZFRAC(IL))*ESW(TX)+ZFRAC(IL)*ESI(TX)
            RSS=EPS1*EST/( P(IL)-EPS2*EST )
            DRSSDT=EPS1*P(IL)
     1             *((1.-ZFRAC(IL))*ESW(TX)*(RW3-RW2/TX)
     2                  +ZFRAC(IL) *ESI(TX)*(RI3-RI2/TX))
     3              /(TX*(P(IL)-EPS2*EST)**2)
C
C           * SECOND ITERATION.
C
            HMNPR=CPM(IL)*TX+GRAV*Z(IL)-HMN(IL)
            DHMNPR=CPM(IL)
            RC=QT(IL)-RSS
            HMNPR=HMNPR-RRL(IL)*RC
            DHMNPR=DHMNPR+RRL(IL)*DRSSDT
            TX=TX-HMNPR/DHMNPR
            EST=(1.-ZFRAC(IL))*ESW(TX)+ZFRAC(IL)*ESI(TX)
            RSS=EPS1*EST/( P(IL)-EPS2*EST )
            DRSSDT=EPS1*P(IL)
     1             *((1.-ZFRAC(IL))*ESW(TX)*(RW3-RW2/TX)
     2                  +ZFRAC(IL) *ESI(TX)*(RI3-RI2/TX))
     3              /(TX*(P(IL)-EPS2*EST)**2)
C
C           * THIRD ITERATION.
C
            HMNPR=CPM(IL)*TX+GRAV*Z(IL)-HMN(IL)
            DHMNPR=CPM(IL)
            RC=QT(IL)-RSS
            HMNPR=HMNPR-RRL(IL)*RC
            DHMNPR=DHMNPR+RRL(IL)*DRSSDT
            TX=TX-HMNPR/DHMNPR
            EST=(1.-ZFRAC(IL))*ESW(TX)+ZFRAC(IL)*ESI(TX)
            RSS=EPS1*EST/( P(IL)-EPS2*EST )
            DRSSDT=EPS1*P(IL)
     1             *((1.-ZFRAC(IL))*ESW(TX)*(RW3-RW2/TX)
     2                  +ZFRAC(IL) *ESI(TX)*(RI3-RI2/TX))
     3              /(TX*(P(IL)-EPS2*EST)**2)
C
C           * TEMPERATURE-DERIVATIVE OF RSS FOR T=TL. THE DERIVATIVE
C           * IS USED IN A TAYLOR SERIES APPROXIMATION TO DETERMINE THE
C           * WATER VAPOUR SATURATION MIXING RATIO AT TEMPERATURE TX
C           * (I.E. AMBIENT TEMPERATURE) ABOUT TL (I.E. LIQUID WATER
C           * TEMPERATURE).
C
            IF( TX.NE.TL .AND. RSS > RSSTL ) THEN
              DRSSDT=(RSS-RSSTL)/(TX-TL)
            ELSE
              DRSSDT=DRSSDTL
              RSS=RSSTL
            ENDIF
          ENDIF
C
C         * SAVE FINAL ITTERATION OF SATURATION SPECIFIC HUMIDITY FOR
C         * USE ELSEWHERE IN CALCULATION RELATIVE HUMIDITIES.
C 
          SSH(IL)=RSS
          IF (ISUBG.EQ.1)                                           THEN
C
C           * A AND B PARAMETERS FOR STATISTICAL CLOUD SCHEME.
C
            APA=1./(1.+RRL(IL)*DRSSDT/CPM(IL))
            APB=APA*DRSSDT
C
C           * DETERMINE VARIANCE BASED ON LOCAL APPROACH (CHABOUREAU AND
C           * BECHTHOLD).
C
            SIGMACB=CSIGMA*ALMIX(IL)*
     1              ABS(APA*DRWDZ(IL)-APB*DHLDZ(IL)/CPM(IL))
C
            IF(ICVSG.EQ.1)                                          THEN         
C
C             * CALCULATE CONVECTIVE CONTRIBUTION TO VARIANCE.
C             * THIS MUST BE OUTPUT BECAUSE OF "AGEING".
C             * NOTE ** THIS IS ONLY DONE ONCE NEAR THE BEGINNING OF THE
C             *         PHYSICS, DUE TO THE "AGEING" CALCULATION (IE
C             *         ONE DOESN'T WANT MULTIPLE AGEING PER TIMESTEP FOR
C             *         EACH PHYSICS ROUTINE WHICH CALLS THIS SUBROUTINE.
C             *         THIS IS CONTROLLED BY THE SWITCH "ICVSG").
C
              IF(CONS(IL).GT.0.) THEN
                CVSG(IL)=3.E-03*CONS(IL)/APA
cold            CVSG(IL)=1.E-05*APA
              ELSE
                CVSG(IL)=CVSG(IL)*SIGFAC
              ENDIF
            ENDIF
C
C           * ADD CONVECTIVE CONTRIBUTION TO VARIANCE.
C
            SIGMASUM=SIGMACB+CVSG(IL)
C
C           * LIMIT VARIANCE TO AVOID UNPHYSICAL CASES WITH LARGE VARIANCES
C           * AN LOW MEAN VALUES FOR THE TOTAL WATER PROBABILITY DISTRIBUTION.
C
            SIGMAMAX=YFR*APA*QT(IL)/2.
            SIGMA(IL)=MIN(SIGMASUM,SIGMAMAX)
          ENDIF
C
C         * CLOUD FRACTION AND CLOUD WATER CONTENT BASED ON CHABOUREAU
C         * AND BECHTHOLD'S STATISTICAL CLOUD SCHEME AND AUTOCONVERSION
C         * SCALING FACTOR BASED ON MELLORS'S STATISTICAL CLOUD SCHEME
C         * (FOR THE 2.47TH MOMENT OF THE TOTAL WATER PROBABILITY
C         * DISTRIBUTION).
C
          QCWT=(QT(IL)-RSSTL)
          IF ( ISUBG .EQ. 1 .AND. SIGMA(IL) .GT. ZSEC ) THEN
            QCWT=QCWT*APA
            Q1=QCWT/SIGMA(IL)
            Q1M=Q1/2.
            INDX=NINT((Q1/2.+2.)/QDELT)+1
            IF ( Q1.LT.Q1MIN ) THEN
              ZCLF(IL)=0.
              QCW (IL)=0.
              ZCRAUT(IL)=ZCR(Q1MINM,I2NDIE)
              QCWVAR(IL)=0.
            ELSE IF ( Q1.GE.Q1MIN .AND. Q1.LT.0. ) THEN
              ZCLF(IL)=0.5+0.36*ATAN(1.55*Q1)
              QCW (IL)=SIGMA(IL)*EXP(1.2*Q1-1.)
              ZCRAUT(IL)=ZCR(Q1M,I2NDIE)
C
              IF (MACHINE.EQ.2) THEN
                R4VAL=Q1M/SQRT2
                AR=.5*(1.+ERF(R4VAL))
              ELSE
                R8VAL=Q1M/SQRT2
                AR=.5*(1.+ERF(R8VAL))
              ENDIF
C
              AEXP=EXP(-.5*Q1M**2)
              QL=AR*Q1M+AEXP/SQRT2PI
              QLVAR=AR*(1.+(Q1M-QL)**2)+AEXP*(Q1M/SQRT2PI-QL)
              QCWVAR(IL)=QLVAR/QL**2
            ELSE IF ( Q1.GE.0. .AND. Q1.LE.Q1MAX ) THEN
              ZCLF(IL)=0.5+0.36*ATAN(1.55*Q1)
              QCW (IL)=SIGMA(IL)*(EXPM1+0.66*Q1+0.086*Q1**2)
              ZCRAUT(IL)=ZCR(Q1M,I2NDIE)
C
              IF (MACHINE.EQ.2) THEN
                R4VAL=Q1M/SQRT2
                AR=.5*(1.+ERF(R4VAL))
              ELSE
                R8VAL=Q1M/SQRT2
                AR=.5*(1.+ERF(R8VAL))
              ENDIF
C
              AEXP=EXP(-.5*Q1M**2)
              QL=AR*Q1M+AEXP/SQRT2PI
              QLVAR=AR*(1.+(Q1M-QL)**2)+AEXP*(Q1M/SQRT2PI-QL)
              QCWVAR(IL)=QLVAR/QL**2
            ELSE
              ZCLF(IL)=1.
              QCW (IL)=QT(IL)-RSS
              ZCRAUT(IL)=1.
              QCWVAR(IL)=0.
            ENDIF
C
C           * IMPOSE UPPER BOUND ON CLOUD WATER. THIS IS NECESSARY SINCE
C           * THE EMPIRICAL FITS ON WHICH THE STATISTICAL CLOUD SCHEME IS
C           * IS BASED MAY PRODUCE UNPHYSICAL RESULTS UNDER CERTAIN
C           * CIRCUMSTANCES (I.E. EXCESSIVELY LARGE CLOUD WATER CONTENTS).
C           * THE APPROACH HERE IS NOT TO ALLOW THE INFERRED WATER VAPOUR
C           * MIXING RATIO IN THE CLEAR-SKY TO BECOME SMALLER THAN A
C           * GIVEN THRESHOLD (ZSEC).
C
            QCWMAX=MAX(QT(IL)-((1.-ZCLF(IL))*ZSEC+ZCLF(IL)*RSS),0.)
            IF ( Q1.GE.Q1MIN .AND. Q1.LE.Q1MAX
     1                                   .AND. QCW(IL).GE.QCWMAX ) THEN
              QCW(IL)=QCWMAX
            ENDIF
          ELSE
            IF ( QCWT .LT. 0. ) THEN
              ZCLF(IL)=0.
              QCW (IL)=0.
              ZCRAUT(IL)=ZCR(Q1MINM,I2NDIE)
              QCWVAR(IL)=0.
            ELSE
              ZCLF(IL)=1.
              QCW (IL)=QT(IL)-RSS
              ZCRAUT(IL)=1.
              QCWVAR(IL)=0.
            ENDIF
          ENDIF
          QCW(IL)=MAX(MIN(QCW(IL),QT(IL)),ZERO)
          ZCLF(IL)=MAX(MIN(ZCLF(IL),ONE),ZERO)
C
C         * ELIMINATE CLOUDSAND CONDENSATE IF CLOUD FRACTION OR
C         * CONDENSATE MIXING RATIO TOO LOW TO ENSURE CONSISTENT
C         * RESULTS FOR CLOUD FRACTION AND CLOUD CONDENSATE.
C
          IF ( ZCLF(IL).LT.0.01 .OR. QCW(IL).LT.1.E-12 ) THEN
            ZCLF(IL)=0.
            QCW (IL)=0.
            ZCRAUT(IL)=ZCR(Q1MINM,I2NDIE)
            QCWVAR(IL)=0.
          ENDIF
C
C         * TRACK UNPHYSICAL VALUES IN PREPARATION FOR ABORTING
C         * (DONE OUTSIDE LOOP IN ORDER TO KEEP OPTIMIZATION).
C
          IF ( ISUBG.EQ.1 .AND. QCW(IL).GT.QT(IL) .AND. ILBAD.EQ.0 )
     1                                                           THEN
            ILBAD=IL
            SIGCBBAD=SIGMACB
            SIGMAXBAD=SIGMAMAX
            QCWTBAD=QCWT 
          ENDIF
        ENDIF
      ENDDO
C
C     * PRINT OUT PROBLEM POINTS AND ABORT, IF THEY EXIST.
C
      IF(ILBAD.NE.0)                                            THEN
        IL=ILBAD
        PRINT*,'ICALL,LVL,IL,QCW,QT,RSS=',ICALL,LVL,ILBAD,QCW(IL),
     1          QT(IL),SSH(IL)
        PRINT*,'SIGMA,SIGMACB,CVSG,SIGMAMAX,ZCLF,ZCRAUT,QCWT=',
     1          SIGMA(IL),SIGCBBAD,CVSG(IL),SIGMAXBAD,ZCLF(IL),
     2          ZCRAUT(IL),QCWTBAD
        CALL XIT('STATCLD5',-1)
      ENDIF
C
      RETURN
      END
