      SUBROUTINE CLDPRPC(Q,T,P,Z,S,MU,EU,DU,MD,ED,SD,QD,MC,QU,SU,ZF,QST,   
     1                   HMN,HSAT,ALPHA,SHAT,QL, 
     2                   EUO,EDO,DZO,DUO,MDO,EPS0O,                    
     3                   JB,LEL,JT,MX,J0,JD,WRK,ITRPHS,
     4                   ILEV,ILG,IL1G,IL2G,C0FAC,RD,GRAV,CP,MSG)       
C
C     * FEB 02/2017 - M.LAZARE.    LOOP 378 MODIFIED TO SUPPORT BOUNDS
C     *                            ARRAY CHECKING EXPLICITLY.
C     * FEB 16/2009 - M.LAZARE.    NEW VERSION FOR GCM15H:
C     *                            - ADD ACCURACY RESTRICTION ON EPS0. 
C     * MAR 27/2007 - M.LAZARE.    PREVIOUS VERSION CLDPRPB FOR GCM15G:
C     *                            - WEIGHT INCREASED FROM 0.25 TO 0.75.
C     *                            - SELECTIVE PROMOTION OF VARIABLES
C     *                              IN REAL*8 MODE TO WORK IN 32-BIT.  
C     *                            - C0FAC NOW PASSED IN FROM CONVECTION
C     *                              DRIVER, TO ENSURE CONSISTENT USEAGE
C     *                              IN CLDPRP AND CONTRA.
C     *                            - WEIGHT CHANGED FROM 0.50 TO 0.25
C     *                              FOR TUNING PURPOSES. 
C     * JUN 15/2006 - M.LAZARE.    PREVIOUS VERSION CLDPRPA FOR GCM15F:
C     *                            - NOW USES INTERNAL WORK ARRAYS
C     *                              INSTEAD OF PASSED WORKSPACE.    
C     *                            - "ZERO","ONE","FMAX" DATA CONSTANTS
C     *                              ADDED AND USED IN CALL TO 
C     *                              INTRINSICS. 
C     * MAY 08/2006 - M.LAZARE/    PREVIOUS VERSION CLDPRPA FOR GCM15E:
C     *               K.VONSALZEN. - CALCULATION OF SATURATION VAPOUR
C     *                              PRESSURE MODIFIED TO WORK OVER
C     *                              COMPLETE VERTICAL DOMAIN, BY
C     *                              SPECIFYING LIMIT AS ES->P.
C     *                            - C0 INCREASED FROM 4.E-3 TO 6.E-3
C     *                              FOR REDUCED CLOUD FORCING IN TROPICS.
C     * DEC 12/2005 - M.LAZARE/    PREVIOUS VERSION CLDPRP9 FOR GCM15D:
C     *               K.VONSALZEN. 

      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
C     * I/O FIELDS:
C
      REAL  , DIMENSION(ILG,ILEV) :: Q,T,P,Z,S,MU,EU,DU,MD,ED,SD,QD,
     1                               MC,QU,SU,ZF,QST,HMN,HSAT,ALPHA,
     2                               SHAT,QL,EUO,DZO,DUO,MDO,EDO
      REAL  , DIMENSION(ILG)      :: EPS0O,WRK
      INTEGER,DIMENSION(ILG)      :: JB,LEL,JT,MX,J0,JD
C                                                                             
C     * INTERNAL WORK FIELDS:                                                         
C
C
C     *** THE FOLLOWING MUST BE PROMOTED FOR ACCURACY!!! *********
C
      REAL*8, DIMENSION(ILG,ILEV) :: IPRM,I1,I2,IHAT,I3,IDAG,I4,F,EPS
      REAL*8, DIMENSION(ILG)      :: HMIN,EXPDIF,EXPNUM,FTEMP,HMAX,ZUEF
      REAL*8  FMAX,ZERO
C     ************************************************************
C
      REAL, DIMENSION(ILG,ILEV) :: GAMMA,DZ,HU,HD,
     1                             QSTHAT,HSTHAT,GAMHAT,CU,QDS 
      REAL, DIMENSION(ILG)      :: EPS0,RMUE,ZDEF,FACT,TU,RL,EST,
     2                             TOTPCP,TOTEVP,ALFA,ZFD,ZFB,
     3                             DENOM,HUB,IRESET,JLCL  
C
      COMMON /EPS   / A,B,EPS1,EPS2                                          
      COMMON /HTCP  / TFREEZ,T2S,AI,BI,AW,BW,SLP
C
C     * PARAMETERS USED IN NEW SATURATION VAPOUR PRESSURE FORMULATION.
C
      COMMON /ESTWI/ RW1,RW2,RW3,RI1,RI2,RI3
C
C     * STATEMENT FUNCTION TO CALCULATE SATURATION VAPOUR PRESSURE
C     * OVER WATER.
C
      ESW(TTT)    = EXP(RW1+RW2/TTT)*TTT**RW3
C
      DATA RRL     / 2.501E6/
      DATA ZERO    / 0.     /
      DATA ONE     / 1.     /
      DATA FMAX    / 0.0002 /
C----------------------------------------------------------------------
      ZEPS=1.E-33                      
      EPSLIM=0.001                                   
      DO 10 K=MSG+1,ILEV-1                                          
      DO 10 IL=IL1G,IL2G                                                      
        DZ(IL,K)=ZF(IL,K)-ZF(IL,K+1)                                           
   10 CONTINUE                                                                 
C
      DO 50 K=MSG+1,ILEV                                                      
      DO 50 IL=IL1G,IL2G                                                      
        MU(IL,K)=0.                                                           
        F(IL,K)=0.                                                            
        EPS(IL,K)=0.                                                          
        EU(IL,K)=0.                                                           
        DU(IL,K)=0.                                                           
        QL(IL,K)=0.                                                           
        CU(IL,K)=0.                                                           
        QDS(IL,K)=Q(IL,K)                                                     
        MD(IL,K)=0.                                                           
        ED(IL,K)=0.                                                           
        SD(IL,K)=S(IL,K)                                                      
        QD(IL,K)=Q(IL,K)                                                      
        MC(IL,K)=0.                                                           
        HD(IL,K)=0.
        I1(IL,K)=0.
        I2(IL,K)=0.
        I3(IL,K)=0.
        I4(IL,K)=0.
        IHAT(IL,K)=0.
        IDAG(IL,K)=0.
        IPRM(IL,K)=0.
        QU(IL,K)=Q(IL,K)                                                      
        SU(IL,K)=S(IL,K)                                                      
        ETMP=ESW(T(IL,K))
        ESTREF=P(IL,K)*(1.-EPSLIM)/(1.-EPSLIM*EPS2)
        IF ( ETMP.LT.ESTREF ) THEN
          EST(IL)=ETMP
        ELSE
          EST(IL)=ESTREF
        ENDIF
        IF(K.GE.LEL(IL))                  THEN
          QST(IL,K)=EPS1*EST(IL)/(P(IL,K)-EST(IL))                             
        ELSE
          QST(IL,K) = Q(IL,K)
        ENDIF
        GAMMA(IL,K)=QST(IL,K)*(1.+QST(IL,K)/EPS1)*EPS1*RRL                     
     1              /(RD*T(IL,K)**2)*RRL/CP                                    
        HMN (IL,K)=CP*T(IL,K)+GRAV*Z(IL,K)+RRL*Q  (IL,K)                      
        HSAT(IL,K)=CP*T(IL,K)+GRAV*Z(IL,K)+RRL*QST(IL,K)                       
        HU(IL,K)=HMN(IL,K)                                                    
   50 CONTINUE
C                                                                       
      DO 75 IL=IL1G,IL2G                                                      
        HSTHAT(IL,MSG+1)=HSAT(IL,MSG+1)                                       
        QSTHAT(IL,MSG+1)=QST(IL,MSG+1)                                        
        GAMHAT(IL,MSG+1)=GAMMA(IL,MSG+1)                                      
        TOTPCP(IL)  =0.                                                       
        TOTEVP(IL)  =0.                                                       
        DZ(IL,ILEV) =ZF(IL,ILEV)
        RL(IL)=RRL
        JT(IL)=LEL(IL)                                                        
        EPS0(IL)=0.
        ALFA(IL)=0.
        EXPDIF(IL)=0.
        HUB(IL)=HU(IL,ILEV)
        J0(IL)=ILEV
        JD(IL)=ILEV
        JLCL(IL)=0.
        HMIN(IL)=1.E6
        HMAX(IL)=HMN(IL,MX(IL))
        ZFB(IL)=ZF(IL,JB(IL)) 
        IRESET(IL)=0.                             
   75 CONTINUE
C
      DO 100 K=MSG+2,ILEV                                                      
      DO 100 IL=IL1G,IL2G                                                      
       QSTHAT(IL,K)=QST(IL,K)
       GAMHAT(IL,K)=GAMMA(IL,K)
       IF(K.GE.LEL(IL))                                             THEN
        IF(QST(IL,K-1).NE.QST(IL,K))                            THEN          
          QSTHAT(IL,K)=LOG(QST(IL,K-1)/QST(IL,K))*QST(IL,K-1)*QST(IL,K)      
     1                 /(QST(IL,K-1)-QST(IL,K))                                
        ENDIF
C                                                                              
        IF(GAMMA(IL,K-1).NE.GAMMA(IL,K))                        THEN          
          GAMHAT(IL,K)=LOG(GAMMA(IL,K-1)/GAMMA(IL,K))                          
     1            *GAMMA(IL,K-1)*GAMMA(IL,K)/(GAMMA(IL,K-1)-GAMMA(IL,K))       
        ENDIF     
       ENDIF         
       HSTHAT(IL,K)=CP*SHAT(IL,K)+RRL*QSTHAT(IL,K)
  100 CONTINUE
C                                                                             
CCC   **********************************************************              
CCC   FIND THE LEVEL OF MINIMUM HSAT, WHERE DETRAINMENT STARTS.
CCC   INITIALIZE CERTAIN ARRAYS INSIDE THE CLOUD PLUME.                    
CCC   **********************************************************              
C                                                                             
      DO 150 J=MSG+1,ILEV                                                     
      DO 150 IL=IL1G,IL2G                                                     
        IF(HSAT(IL,J).LE.HMIN(IL).AND.(J.GE.JT(IL).AND.J.LE.JB(IL)))THEN      
           HMIN(IL)=HSAT(IL,J)                                                 
           J0(IL)  =J                                                         
        ENDIF
        IF(J.GE.JT(IL) .AND. J.LE.JB(IL))                           THEN      
           F(IL,J)= 0.                                                       
          HU(IL,J)= HMAX(IL)                                                
        ENDIF                                                                 
  150 CONTINUE
C
C     * DEFINE BOUNDS FOR DETRAINMENT LEVEL.
C     * RE-INITIALIZE "HMIN" FOR ENSUING CALCULATIONS.
C                                                   
      DO 175 IL=IL1G,IL2G  
        J0(IL)=MIN(J0(IL),JB(IL)-2)
        J0(IL)=MAX(J0(IL),JT(IL)+2) 
        HMIN(IL)=1.E6                                                         
  175 CONTINUE      
C                                                                             
CCCC  *********************************************************               
CCCC  COMPUTE TAYLOR SERIES FOR APPROXIMATE EPS(Z) BELOW                      
CCCC  *********************************************************               
C                                                                             
      DO 275 J=ILEV-1,MSG+1,-1                                                
      DO 275 IL=IL1G,IL2G                                                     
        IF(J.LT.JB(IL) .AND. J.GE.JT(IL))                           THEN      
          I1  (IL,J) = I1(IL,J+1)+(HMAX(IL)-HMN(IL,J))*DZ(IL,J)             
          IHAT(IL,J) = 0.5*(I1(IL,J+1)+I1(IL,J))                              
          I2  (IL,J) = I2(IL,J+1)+IHAT(IL,J)*DZ(IL,J)                         
          IDAG(IL,J) = 0.5*(I2(IL,J+1)+I2(IL,J))                              
          I3  (IL,J) = I3(IL,J+1)+IDAG(IL,J)*DZ(IL,J)                         
          IPRM(IL,J) = 0.5*(I3(IL,J+1)+I3(IL,J))                              
          I4  (IL,J) = I4(IL,J+1)+IPRM(IL,J)*DZ(IL,J)                         
        ENDIF                                                                 
  275 CONTINUE       
C
      DO 350 J=MSG+1,ILEV                                                     
      DO 350 IL=IL1G,IL2G                                                     
        IF( (J.GE.J0(IL) .AND. J.LE.JB(IL)) .AND.                             
     1      HMN(IL,J).LE.HMIN(IL)                )                  THEN      
          HMIN(IL)=HMN(IL,J)
          EXPDIF(IL)=HMAX(IL)-HMIN(IL)                                      
        ENDIF                                                                 
  350 CONTINUE       
C                                                                             
CCCC  *********************************************************               
CCCC  COMPUTE APPROXIMATE EPS(Z) USING ABOVE TAYLOR SERIES                    
CCCC  *********************************************************               
C                    
      DO 375 J=MSG+2,ILEV                                                     
      DO 375 IL=IL1G,IL2G
        IF(J.LT.JT(IL) .OR. J.GT.JB(IL))                            THEN      
          I1(IL,J)=0.                                                         
          EXPNUM(IL)=0.                                                       
        ELSE                                                                  
          EXPNUM(IL)=HMAX(IL)- ( HSAT(IL,J-1)*(ZF(IL,J)-Z(IL,J))            
     1          +HSAT(IL,J)*(Z(IL,J-1)-ZF(IL,J)) )/(Z(IL,J-1)-Z(IL,J))        
        ENDIF                                                                 
        IF( (EXPDIF(IL).GT.100. .AND. EXPNUM(IL).GT.0.) .AND.                 
     1   I1(IL,J).GT.EXPNUM(IL)*DZ(IL,J))                        THEN         
          FTEMP(IL)=EXPNUM(IL)/I1(IL,J)                                       
          F(IL,J)=FTEMP(IL)                                                   
     1            +I2(IL,J)/I1(IL,J)*FTEMP(IL)**2                             
     2            +(2.*I2(IL,J)**2-I1(IL,J)*I3(IL,J))/I1(IL,J)**2             
     3            *FTEMP(IL)**3                                               
     4            +(-5.*I1(IL,J)*I2(IL,J)*I3(IL,J)+5.*I2(IL,J)**3             
     5            +I1(IL,J)**2*I4(IL,J))/I1(IL,J)**3*FTEMP(IL)**4             
          F(IL,J)=MAX(F(IL,J),ZERO)                                             
          F(IL,J)=MIN(F(IL,J),FMAX)                  
        ENDIF                                                                 
  375 CONTINUE      
C
        DO 378 IL=IL1G,IL2G                                                   
          IF(F(IL,J0(IL)).LT.1.E-6 .AND. J0(IL).LT.JB(IL) .AND.               
     1       F(IL,MIN(J0(IL)+1,ILEV)) .GT. F(IL,J0(IL)) )          THEN
             J0(IL)=J0(IL)+1         
          ENDIF
  378   CONTINUE     
C
       DO 380 J=MSG+2,ILEV                                                    
       DO 380 IL=IL1G,IL2G                                                    
         IF(J.GE.JT(IL) .AND. J.LE.J0(IL))                         THEN       
            F(IL,J)=MAX(F(IL,J),F(IL,J-1))                                    
         ENDIF                                                                
  380 CONTINUE       
C                                                                           
      DO 425 J=ILEV,MSG+1,-1                                                  
      DO 425 IL=IL1G,IL2G
        IF(J.GE.JT(IL) .AND. J.LT.J0(IL))                          THEN
          EPS(IL,J)=F(IL,J)               
        ELSE IF(J.GE.J0(IL) .AND. J.LE.JB(IL))                     THEN      
          EPS(IL,J)=F(IL,J0(IL))                                             
        ENDIF
  425 CONTINUE                                                                
C                                                                             
CCC   ****************************************************************        
CCC   SPECIFY THE UPDRAFT MASS FLUX MU, ENTRAINMENT EU, DETRAINMENT DU        
CCC   AND MOIST STATIC ENERGY HU.                                             
CCC   HERE AND BELOW MU, EU,DU, MD AND ED ARE ALL NORMALIZED BY MB            
CCC   ****************************************************************        
C                                                                             
      DO 475 IL=IL1G,IL2G
        EPS0(IL)=F(IL,J0(IL))    
        IF(EPS0(IL).LT.1.E-6) EPS0(IL)=0.                                
        IF(EPS0(IL).GT.0.)                                          THEN  
          MU(IL,JB(IL))=1.
          HUB(IL)      =HU(IL,JB(IL))
          EU(IL,JB(IL))=EPS0(IL)/2.                                        
        ENDIF                                                              
  475 CONTINUE
C                                                                          
      DO 480 J=ILEV-1,MSG+1,-1
      DO 480 IL=IL1G,IL2G                                                      
        IF( EPS0(IL).GT.0. .AND. (J.GE.JT(IL).AND.J.LT.JB(IL)) )    THEN       
          ZUEF(IL)=ZF(IL,J)-ZFB(IL)                                         
          RMUE(IL)=(1./EPS0(IL))*(EXP(EPS(IL,J+1)*ZUEF(IL))-1.)/ZUEF(IL)       
          MU(IL,J)=(1./EPS0(IL))*(EXP(EPS(IL,J)*ZUEF(IL))-1.)/ZUEF(IL)         
          EU(IL,J)=(RMUE(IL)-MU(IL,J+1))/DZ(IL,J)                              
          DU(IL,J)=(RMUE(IL)-MU(IL,J))/DZ(IL,J)                                
        ENDIF                                                                  
  480 CONTINUE       
C                                    
      DO 500 J=ILEV-1,MSG+1,-1                                             
      DO 500 IL=IL1G,IL2G                                                     
        IF( (J.GE.LEL(IL) .AND. J.LT.JB(IL)) .AND.
     1    IRESET(IL).EQ.0. .AND. EPS0(IL).GT.0.)                    THEN       
          IF(MU(IL,J).LT.0.01 .OR. HU(IL,J+1).GT.HUB(IL))        THEN          
            HU(IL,J)=HUB(IL)                                                
            MU(IL,J)=0.             
            JT(IL)=J+1
            IRESET(IL)=1.                                               
C           PRINT *,'JT DUE TO TOO MUCH DETRAINMENT  ===',JT(IL),          
C    1              'IL= ',IL        
          ELSE                                                                
            HU(IL,J)=MU(IL,J+1)/MU(IL,J)*HU(IL,J+1)+DZ(IL,J)/MU(IL,J)*        
     1               (EU(IL,J)*HMN(IL,J)-DU(IL,J)*HSAT(IL,J))
            IF (HU(IL,J).LE.HSTHAT(IL,J)   .AND.                               
     1          HU(IL,J+1).GT.HSTHAT(IL,J+1))         THEN
              JT(IL)=J
              IRESET(IL)=2.                                                    
C             PRINT *,'JT DUE TO OVERSHOOTING ===',JT(IL),' IL= ',IL          
            ENDIF
          ENDIF
        ENDIF                                                                 
  500 CONTINUE       
C
      DO 600 J=ILEV-1,MSG+1,-1
      DO 600 IL=IL1G,IL2G                                                     
        IF( J.GE.LEL(IL) .AND. J.LE.JT(IL)                                    
     1                     .AND. EPS0(IL).GT.0. )                   THEN      
          MU(IL,J)=0.                                                         
          EU(IL,J)=0.                                                         
          DU(IL,J)=0.                                                         
          HU(IL,J)=HUB(IL)                                              
        ENDIF                                                                 
        IF(J.EQ.JT(IL) .AND. EPS0(IL).GT.0.)                        THEN      
          DU(IL,J)=MU(IL,J+1)/DZ(IL,J)                                        
        ENDIF                                                                 
  600 CONTINUE       
C                                                                              
CCC   ******************************************************************      
CCC   SPECIFY DOWNDRAFT PROPERTIES (NO DOWNDRAFTS IF JD.GE.JB).               
C     SCALE DOWN DOWNWARD MASS FLUX PROFILE SO THAT NET FLUX                  
C     (UP-DOWN) AT CLOUD BASE IN NOT NEGATIVE.                                 
CCC   ************************************************************            
C                    
      DO 625 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0.)                                          THEN
          JD(IL)=MAX(J0(IL),JT(IL)+1)
          ZFD(IL)=ZF(IL,JD(IL)) 
          HD(IL,JD(IL))=HMN(IL,JD(IL)-1)
          FACT(IL)=2.*EPS0(IL)*(ZFD(IL)-ZFB(IL))  
          IF(FACT(IL).GT.0.)                                      THEN      
            ALFA(IL)=FACT(IL)/(EXP(FACT(IL))-1.)            
          ELSE
            ALFA(IL)=0.
          ENDIF                                                            
        ENDIF
  625 CONTINUE       
C
      DO 630 J=MSG+1,ILEV                                           
      DO 630 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0.)                                          THEN   
          ZDEF(IL)=ZFD(IL)-ZF(IL,J)
          IF((J.EQ.JD(IL)) .AND. FACT(IL).GT.0.)                  THEN
            MD(IL,J)=-1.
          ENDIF
          IF((J.EQ.JB(IL)) .AND. FACT(IL).GT.0.)                  THEN
            MD(IL,J)=-1.
          ENDIF
          IF((J.LT.JB(IL)) .AND. ZDEF(IL).GT.0.)                  THEN
            MD(IL,J)=-1.
          ENDIF
        ENDIF
  630 CONTINUE       
C
      DO 650 J=MSG+2,ILEV    
      DO 650 IL=IL1G,IL2G                                                     
        IF((J.GT.JD(IL) .AND. J.LE.JB(IL)) .AND. EPS0(IL).GT.0.)    THEN      
          ED(IL,J-1)=(MD(IL,J-1)-MD(IL,J))/DZ(IL,J-1)                         
          HD(IL,J)=MD(IL,J-1)/MD(IL,J)*HD(IL,J-1)                             
     1             -DZ(IL,J-1)/MD(IL,J)*ED(IL,J-1)*HMN(IL,J-1)                
        ENDIF                                                                 
  650 CONTINUE       
C                                                                             
CCC   *********************************************************               
CCC   CALCULATE UPDRAFT AND DOWNDRAFT EFFECTS ON Q AND S.                     
CCC   *********************************************************               
C                    
      DO 670 K=MSG+1,ILEV             
      DO 670 IL=IL1G,IL2G                                                     
        IF((K.GE.JD(IL) .AND. K.LE.JB(IL)) .AND. EPS0(IL).GT.0.               
     1                                     .AND. JD(IL).LT.JB(IL))  THEN    
C         SD(IL,K) = SHAT(IL,K)                                               
C    1             +              (HD(IL,K)-HSTHAT(IL,K))/                    
C    2               (CP    *(1.+GAMHAT(IL,K)))                               
          QDS(IL,K) = QSTHAT(IL,K)                                            
     1             + GAMHAT(IL,K)*(HD(IL,K)-HSTHAT(IL,K))/                    
     2               (RRL*(1.+GAMHAT(IL,K)))                                   
        ENDIF                                                                 
  670 CONTINUE
C
      DO 690 J=ILEV-1,MSG+2,-1
       DO 680 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0. .AND. JLCL(IL).EQ.0.)                     THEN       
          IF(J.GT.JT(IL) .AND. J.LT.JB(IL))     THEN          
            SU(IL,J)=MU(IL,J+1)/MU(IL,J)*SU(IL,J+1)                           
     1              +DZ(IL,J)/MU(IL,J)*(EU(IL,J)-DU(IL,J))*S(IL,J)            
            QU(IL,J)=MU(IL,J+1)/MU(IL,J)*QU(IL,J+1)                           
     1        +DZ(IL,J)/MU(IL,J)*(EU(IL,J)*Q(IL,J)-DU(IL,J)*QST(IL,J))
            TU(IL)=SU(IL,J)-GRAV/CP*ZF(IL,J)                                   
            ETMP=ESW(TU(IL))
            ESTREF=P(IL,J)*(1.-EPSLIM)/(1.-EPSLIM*EPS2)
            IF ( ETMP.LT.ESTREF ) THEN
              ESTU=ETMP
            ELSE
              ESTU=ESTREF
            ENDIF
            QSTU=EPS1*ESTU/((P(IL,J)+P(IL,J-1))/2.-ESTU)                      
            IF(QU(IL,J) .GE. QSTU)                     THEN            
              JLCL(IL)=REAL(J)                                                 
            ENDIF
          ENDIF
        ENDIF                                                            
 680   CONTINUE                                                                
 690  CONTINUE      
C
      DO 700 J=MSG+1,ILEV      
      DO 700 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0.)                                          THEN      
          IF(J.EQ.JB(IL))                                     THEN        
            QU(IL,J)=Q(IL,MX(IL))                                             
            SU(IL,J)=(HU(IL,J)-RRL*QU(IL,J))/CP                               
          ENDIF                                                              
          IF(J.GT.JT(IL) .AND. J.LE.NINT(JLCL(IL)))           THEN          
            SU(IL,J) = SHAT(IL,J)                                             
     1               +              (HU(IL,J)-HSTHAT(IL,J))/                   
     2               (CP    *(1.+GAMHAT(IL,J)))                               
            QU(IL,J) = QSTHAT(IL,J)                                            
     1               + GAMHAT(IL,J)*(HU(IL,J)-HSTHAT(IL,J))/                   
     2               (RRL*(1.+GAMHAT(IL,J)))
          ENDIF
        ENDIF                                                                 
  700 CONTINUE       
C
      DO 750 J=ILEV-1,MSG+1,-1                                    
      DO 750 IL=IL1G,IL2G                                                     
        IF(EPS0(IL) .GT. 0.)                                      THEN        
          IF(J.GT.JT(IL) .AND. J.LT.JB(IL))       THEN        
            CU(IL,J)=((MU(IL,J)*SU(IL,J)-MU(IL,J+1)*SU(IL,J+1))/DZ(IL,J)
     1              -(EU(IL,J)-DU(IL,J))*S(IL,J))/(RRL/CP)
          ENDIF
        ENDIF  
  750 CONTINUE       
C                                                                              
      BETA=1.                                                                
      DO 800 J=ILEV,MSG+2,-1                                                 
      DO 800 IL=IL1G,IL2G
        IF( J.GE.JT(IL) .AND. J.LT.JB(IL) .AND. EPS0(IL).GT.0.               
     1                  .AND. MU(IL,J).GE.0.0 )                     THEN     
          FAC=MAX(T(IL,JB(IL))-T(IL,J),0.)/
     1        MAX(T(IL,JB(IL))-TFREEZ, ONE)
          C0=C0FAC*MIN(FAC**2,ONE)
          DENOM(IL)=(MU(IL,J)+DZ(IL,J)*BETA*DU(IL,J)                         
     1          +DZ(IL,J)*C0*(BETA*MU(IL,J)+(1.-BETA)*MU(IL,J+1))*BETA)      
          IF(DENOM(IL).GT.0.)                                 THEN          
             QL(IL,J)=1./DENOM(IL)                                            
     1        *( MU(IL,J+1)*QL(IL,J+1)-DZ(IL,J)*DU(IL,J)*(1.-BETA)            
     2        *QL(IL,J+1)+DZ(IL,J)*CU(IL,J)-DZ(IL,J)*C0*(BETA*MU(IL,J)        
     3        +(1.-BETA)*MU(IL,J+1))*(1.-BETA)*QL(IL,J+1))                    
          ELSE                                                                
             QL(IL,J)=0.                                                      
          ENDIF                                                               
          TOTPCP(IL)=TOTPCP(IL)+DZ(IL,J)*(CU(IL,J)                            
     1               -DU(IL,J)*(BETA*QL(IL,J)+(1.-BETA)*QL(IL,J+1)))          
        ENDIF                                                                 
  800 CONTINUE       
C                    
      DDRAT=0.
      DO 850 J=MSG+1,ILEV
      DO 850 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0. .AND. J.EQ.JD(IL))                        THEN
          QD(IL,J)=QDS(IL,J) - DDRAT*(QDS(IL,J) -Q(IL,J-1))
          SD(IL,J)=(HD(IL,J)-RRL*QD(IL,J) )/CP
        ENDIF
  850 CONTINUE       
C
      DO 900 J=MSG+1,ILEV-1                                                  
      DO 900 IL=IL1G,IL2G                                                      
        IF((J.GE.JD(IL) .AND. J.LT.JB(IL)) .AND. EPS0(IL).GT.0.               
     1                                     .AND. JD(IL).LT.JB(IL))  THEN      
          QD(IL,J+1)=QDS(IL,J+1) - DDRAT*(QDS(IL,JD(IL))-QD(IL,JD(IL)))
          SD(IL,J+1)=(HD(IL,J+1)-RRL*QD(IL,J+1))/CP                            
          TOTEVP(IL)=TOTEVP(IL)-DZ(IL,J)*ED(IL,J)*Q(IL,J)                     
        ENDIF                                                                
  900 CONTINUE       
C
C     * SAVE CERTAIN FIELDS TO BE USED IN SUBROUTINE CONTRA3 FOR
C     * SULFUR CYCLE, RATHER THAN HAVING TO RECOMPUTE THEM.
C
      DO 923 J=MSG+1,ILEV                                                 
      DO 923 IL=IL1G,IL2G                                                 
         DZO(IL,J) = DZ(IL,J)                                             
         EUO(IL,J) = EU(IL,J)                
         EDO(IL,J) = ED(IL,J)                
         DUO(IL,J) = DU(IL,J)
         MDO(IL,J) = MD(IL,J)
  923 CONTINUE                                                            
C
      DO 925 IL=IL1G,IL2G     
        EPS0O(IL) = EPS0(IL)                                               
        IF(EPS0(IL).GT.0.)                                          THEN
          TOTEVP(IL)=TOTEVP(IL)+MD(IL,JD(IL))*Q(IL,JD(IL)-1)                
     1                         -MD(IL,JB(IL))*QD(IL,JB(IL))                 
        ELSE
C
C         * RESET INDICES TO ILEV TO THAT MB WILL BE ZERO FOR REST OF CODE.
C       
          J0(IL)=ILEV
          JD(IL)=ILEV
          JT(IL)=ILEV
        ENDIF
  925 CONTINUE       
C
      WEIGHT=0.75   
      DO 950 J=MSG+2,ILEV
      DO 950 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0. .AND. TOTEVP(IL).GT.0.)  THEN                   
          RATIO = TOTPCP(IL)/TOTEVP(IL)  
          MD(IL,J)=MD(IL,J)*WEIGHT*RATIO/(1.+WEIGHT*RATIO)                   
          ED(IL,J)=ED(IL,J)*WEIGHT*RATIO/(1.+WEIGHT*RATIO)
        ELSE
          MD(IL,J) = 0.
          ED(IL,J) = 0.                    
        ENDIF            
  950 CONTINUE
C
C     * SAVE FRACTION OF RAIN WHICH EVAPORATES IN THE DOWNDRAFTS.
C
      DO 960 IL=IL1G,IL2G                                               
         WRK(IL)=0.
         IF (EPS0(IL).GT.0. .AND. TOTEVP(IL).NE.0.) THEN
            WRK(IL)=MIN(ONE, WEIGHT*TOTEVP(IL)                    
     1                   /(TOTEVP(IL)+WEIGHT*TOTPCP(IL)))
         ENDIF
 960  CONTINUE
C
      DO 975 J=MSG+1,ILEV                                                     
      DO 975 IL=IL1G,IL2G                                                     
        IF(EPS0(IL).GT.0. .AND. (J.GE.JT(IL) .AND. J.LE.JB(IL)))    THEN
          MC(IL,J)=MU(IL,J)+MD(IL,J)                                          
        ENDIF                                                                 
  975 CONTINUE       
C
      RETURN                                                            
      END       
