!
! $Header$
!
      SUBROUTINE conccm (dtime,paprs,pplay,t,q,conv_q,
     s                   d_t, d_q, rain, snow, kbascm, ktopcm)
c
      IMPLICIT none
c======================================================================
c Auteur(s): Z.X. Li (LMD/CNRS) date: le 14 mars 1996
c Objet: Schema simple (avec flux de masse) pour la convection 
c        (schema standard du modele NCAR CCM2)
c======================================================================
#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "YOETHF.h"
c
c Entree:
      REAL dtime              ! pas d'integration
      REAL paprs(klon,klev+1) ! pression inter-couche (Pa)
      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
      REAL t(klon,klev)       ! temperature (K)
      REAL q(klon,klev)       ! humidite specifique (g/g)
      REAL conv_q(klon,klev)  ! taux de convergence humidite (g/g/s)
c Sortie:
      REAL d_t(klon,klev)     ! incrementation temperature
      REAL d_q(klon,klev)     ! incrementation vapeur
      REAL rain(klon)         ! pluie (mm/s)
      REAL snow(klon)         ! neige (mm/s)
      INTEGER kbascm(klon)    ! niveau du bas de convection
      INTEGER ktopcm(klon)    ! niveau du haut de convection
c
      REAL pt(klon,klev)
      REAL pq(klon,klev)
      REAL pres(klon,klev)
      REAL dp(klon,klev)
      REAL zgeom(klon,klev)
      REAL cmfprs(klon)
      REAL cmfprt(klon)
      INTEGER ntop(klon)
      INTEGER nbas(klon)
      INTEGER i, k
      REAL zlvdcp, zlsdcp, zdelta, zz, za, zb
c
      LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
      PARAMETER (usekuo=.TRUE.)
c
      REAL d_t_bis(klon,klev)
      REAL d_q_bis(klon,klev)
      REAL rain_bis(klon)
      REAL snow_bis(klon)
      INTEGER ibas_bis(klon)
      INTEGER itop_bis(klon)
      REAL d_ql_bis(klon,klev)
      REAL rneb_bis(klon,klev)
c
c initialiser les variables de sortie (pour securite)
      DO i = 1, klon
         rain(i) = 0.0
         snow(i) = 0.0
         kbascm(i) = 0
         ktopcm(i) = 0
      ENDDO
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
         d_q(i,k) = 0.0
      ENDDO
      ENDDO
c
c preparer les variables d'entree (attention: l'ordre des niveaux 
c verticaux augmente du haut vers le bas)
      DO k = 1, klev
      DO i = 1, klon
         pt(i,k) = t(i,klev-k+1)
         pq(i,k) = q(i,klev-k+1)
         pres(i,k) = pplay(i,klev-k+1)
         dp(i,k) = paprs(i,klev+1-k)-paprs(i,klev+1-k+1)
      ENDDO
      ENDDO
      DO i = 1, klon
         zgeom(i,klev) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
     .                      * (paprs(i,1)-pplay(i,1))
      ENDDO
      DO k = 2, klev
      DO i = 1, klon
         zgeom(i,klev+1-k) = zgeom(i,klev+1-k+1)
     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
     .                   * (pplay(i,k-1)-pplay(i,k))
      ENDDO
      ENDDO
c
      CALL cmfmca(dtime, pres, dp, zgeom, pt, pq,
     $                  cmfprt, cmfprs, ntop, nbas)
C
      DO k = 1, klev
      DO i = 1, klon
         d_q(i,klev+1-k) = pq(i,k) - q(i,klev+1-k) 
         d_t(i,klev+1-k) = pt(i,k) - t(i,klev+1-k)
      ENDDO
      ENDDO
c
      DO i = 1, klon
         rain(i) = cmfprt(i) * rhoh2o
         snow(i) = cmfprs(i) * rhoh2o
         kbascm(i) = klev+1 - nbas(i)
         ktopcm(i) = klev+1 - ntop(i)
      ENDDO
c
      IF (usekuo) THEN
      CALL conkuo(dtime, paprs, pplay, t, q, conv_q,
     s            d_t_bis, d_q_bis, d_ql_bis, rneb_bis,
     s            rain_bis, snow_bis, ibas_bis, itop_bis)
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = d_t(i,k) + d_t_bis(i,k)
         d_q(i,k) = d_q(i,k) + d_q_bis(i,k)
      ENDDO
      ENDDO
      DO i = 1, klon
         rain(i) = rain(i) + rain_bis(i)
         snow(i) = snow(i) + snow_bis(i)
         kbascm(i) = MIN(kbascm(i),ibas_bis(i))
         ktopcm(i) = MAX(ktopcm(i),itop_bis(i))
      ENDDO
      DO k = 1, klev ! eau liquide convective est
      DO i = 1, klon ! dispersee dans l'air
         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q(i,k))
         zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q(i,k))
         zdelta = MAX(0.,SIGN(1.,RTT-t(i,k)))
         zz = d_ql_bis(i,k) ! re-evap. de l'eau liquide
         zb = MAX(0.0,zz)
         za = - MAX(0.0,zz) * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
         d_t(i,k) = d_t(i,k) + za
         d_q(i,k) = d_q(i,k) + zb
      ENDDO
      ENDDO
      ENDIF
c
      RETURN
      END
      SUBROUTINE cmfmca(deltat, p, dp, gz,
     $                  tb, shb,
     $                  cmfprt, cmfprs, cnt, cnb)
      IMPLICIT none
C-----------------------------------------------------------------------
C Moist convective mass flux procedure:
C If stratification is unstable to nonentraining parcel ascent,
C complete an adjustment making use of a simple cloud model
C 
C Code generalized to allow specification of parcel ("updraft")
C properties, as well as convective transport of an arbitrary
C number of passive constituents (see cmrb array).
C----------------------------Code History-------------------------------
C Original version:  J. J. Hack, March 22, 1990
C Standardized:      J. Rosinski, June 1992
C Reviewed:          J. Hack, G. Taylor, August 1992
c Adaptation au LMD: Z.X. Li, mars 1996 (reference: Hack 1994,
c                    J. Geophys. Res. vol 99, D3, 5551-5568). J'ai
c                    introduit les constantes et les fonctions thermo-
c                    dynamiques du Centre Europeen. J'ai elimine le
c                    re-indicage du code en esperant que cela pourra
c                    simplifier la lecture et la comprehension.
C-----------------------------------------------------------------------
#include "dimensions.h"
#include "dimphy.h"
      INTEGER pcnst ! nombre de traceurs passifs
      PARAMETER (pcnst=1)
C------------------------------Arguments--------------------------------
C Input arguments
C
      REAL deltat                 ! time step (seconds)
      REAL p(klon,klev)           ! pressure
      REAL dp(klon,klev)          ! delta-p
      REAL gz(klon,klev)          ! geopotential (a partir du sol)
c
      REAL thtap(klon)            ! PBL perturbation theta
      REAL shp(klon)              ! PBL perturbation specific humidity 
      REAL pblh(klon)             ! PBL height (provided by PBL routine)
      REAL cmrp(klon,pcnst)       ! constituent perturbations in PBL
c
c Updated arguments:
c
      REAL tb(klon,klev)         ! temperature (t bar)
      REAL shb(klon,klev)        ! specific humidity (sh bar)
      REAL cmrb(klon,klev,pcnst) ! constituent mixing ratios (cmr bar)
C
C Output arguments
C
      REAL cmfdt(klon,klev)    ! dT/dt due to moist convection
      REAL cmfdq(klon,klev)    ! dq/dt due to moist convection
      REAL cmfmc(klon,klev )   ! moist convection cloud mass flux
      REAL cmfdqr(klon,klev)   ! dq/dt due to convective rainout 
      REAL cmfsl(klon,klev )   ! convective lw static energy flux
      REAL cmflq(klon,klev )   ! convective total water flux
      REAL cmfprt(klon)        ! convective precipitation rate
      REAL cmfprs(klon)        ! convective snowfall rate
      REAL qc(klon,klev)       ! dq/dt due to rainout terms
      INTEGER cnt(klon)        ! top level of convective activity   
      INTEGER cnb(klon)        ! bottom level of convective activity
C------------------------------Parameters-------------------------------
      REAL c0         ! rain water autoconversion coefficient
      PARAMETER (c0=1.0e-4)
      REAL dzmin       ! minimum convective depth for precipitation
      PARAMETER (dzmin=0.0)
      REAL betamn      ! minimum overshoot parameter
      PARAMETER (betamn=0.10)
      REAL cmftau      ! characteristic adjustment time scale
      PARAMETER (cmftau=3600.)
      INTEGER limcnv   ! top interface level limit for convection
      PARAMETER (limcnv=1)
      REAL tpmax       ! maximum acceptable t perturbation (degrees C)
      PARAMETER (tpmax=1.50)
      REAL shpmax      ! maximum acceptable q perturbation (g/g)
      PARAMETER (shpmax=1.50e-3)
      REAL tiny        ! arbitrary small num used in transport estimates
      PARAMETER (tiny=1.0e-36)
      REAL eps         ! convergence criteria (machine dependent)
      PARAMETER (eps=1.0e-13)
      REAL tmelt       ! freezing point of water(req'd for rain vs snow)
      PARAMETER (tmelt=273.15)
      REAL ssfac ! supersaturation bound (detrained air)
      PARAMETER (ssfac=1.001)
C
C---------------------------Local workspace-----------------------------
      REAL gam(klon,klev)     ! L/cp (d(qsat)/dT)
      REAL sb(klon,klev)      ! dry static energy (s bar)
      REAL hb(klon,klev)      ! moist static energy (h bar)
      REAL shbs(klon,klev)    ! sat. specific humidity (sh bar star)
      REAL hbs(klon,klev)     ! sat. moist static energy (h bar star)
      REAL shbh(klon,klev+1)  ! specific humidity on interfaces
      REAL sbh(klon,klev+1)   ! s bar on interfaces
      REAL hbh(klon,klev+1)   ! h bar on interfaces
      REAL cmrh(klon,klev+1)  ! interface constituent mixing ratio 
      REAL prec(klon)         ! instantaneous total precipitation
      REAL dzcld(klon)        ! depth of convective layer (m)
      REAL beta(klon)         ! overshoot parameter (fraction)
      REAL betamx             ! local maximum on overshoot
      REAL eta(klon)          ! convective mass flux (kg/m^2 s)
      REAL etagdt             ! eta*grav*deltat
      REAL cldwtr(klon)       ! cloud water (mass)
      REAL rnwtr(klon)        ! rain water  (mass)
      REAL sc  (klon)         ! dry static energy   ("in-cloud")
      REAL shc (klon)         ! specific humidity   ("in-cloud")
      REAL hc  (klon)         ! moist static energy ("in-cloud")
      REAL cmrc(klon)         ! constituent mix rat ("in-cloud")
      REAL dq1(klon)          ! shb  convective change (lower lvl)
      REAL dq2(klon)          ! shb  convective change (mid level)
      REAL dq3(klon)          ! shb  convective change (upper lvl)
      REAL ds1(klon)          ! sb   convective change (lower lvl)
      REAL ds2(klon)          ! sb   convective change (mid level)
      REAL ds3(klon)          ! sb   convective change (upper lvl)
      REAL dcmr1(klon)        ! cmrb convective change (lower lvl)
      REAL dcmr2(klon)        ! cmrb convective change (mid level)
      REAL dcmr3(klon)        ! cmrb convective change (upper lvl)
      REAL flotab(klon)       ! hc - hbs (mesure d'instabilite)
      LOGICAL ldcum(klon)     ! .true. si la convection existe
      LOGICAL etagt0          ! true if eta > 0.0
      REAL dt                 ! current 2 delta-t (model time step)
      REAL cats     ! modified characteristic adj. time
      REAL rdt      ! 1./dt
      REAL qprime   ! modified specific humidity pert.
      REAL tprime   ! modified thermal perturbation
      REAL pblhgt   ! bounded pbl height (max[pblh,1m])
      REAL fac1     ! intermediate scratch variable
      REAL shprme   ! intermediate specific humidity pert.
      REAL qsattp   ! saturation mixing ratio for 
C                   !  thermally perturbed PBL parcels 
      REAL dz       ! local layer depth
      REAL b1       ! bouyancy measure in detrainment lvl
      REAL b2       ! bouyancy measure in condensation lvl
      REAL g     ! bounded vertical gradient of hb
      REAL tmass ! total mass available for convective exchange
      REAL denom ! intermediate scratch variable
      REAL qtest1! used in negative q test (middle lvl) 
      REAL qtest2! used in negative q test (lower lvl) 
      REAL fslkp ! flux lw static energy (bot interface)
      REAL fslkm ! flux lw static energy (top interface)
      REAL fqlkp ! flux total water (bottom interface)
      REAL fqlkm ! flux total water (top interface)
      REAL botflx! bottom constituent mixing ratio flux
      REAL topflx! top constituent mixing ratio flux
      REAL efac1 ! ratio cmrb to convectively induced change (bot lvl)
      REAL efac2 ! ratio cmrb to convectively induced change (mid lvl)
      REAL efac3 ! ratio cmrb to convectively induced change (top lvl)
c
      INTEGER i,k  ! indices horizontal et vertical
      INTEGER km1  ! k-1 (index offset)
      INTEGER kp1  ! k+1 (index offset)
      INTEGER m    ! constituent index
      INTEGER ktp  ! temporary index used to track top 
      INTEGER is   ! nombre de points a ajuster
C
      REAL tmp1, tmp2, tmp3, tmp4
      REAL zx_t, zx_p, zx_q, zx_qs, zx_gam
      REAL zcor, zdelta, zcvm5
C
      REAL qhalf, sh1, sh2, shbs1, shbs2
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
      qhalf(sh1,sh2,shbs1,shbs2) = MIN(MAX(sh1,sh2),
     $                            (shbs2*sh1 + shbs1*sh2)/(shbs1+shbs2))
C
C-----------------------------------------------------------------------
c pas de traceur pour l'instant
      DO m = 1, pcnst
      DO k = 1, klev
      DO i = 1, klon
         cmrb(i,k,m) = 0.0
      ENDDO
      ENDDO
      ENDDO
c
c Les perturbations de la couche limite sont zero pour l'instant
c
      DO m = 1, pcnst
      DO i = 1, klon
         cmrp(i,m) = 0.0
      ENDDO
      ENDDO
      DO i = 1, klon
         thtap(i) = 0.0
         shp(i) = 0.0
         pblh(i) = 1.0
      ENDDO
C
C Ensure that characteristic adjustment time scale (cmftau) assumed
C in estimate of eta isn't smaller than model time scale (deltat)
C
      dt   = deltat
      cats = MAX(dt,cmftau)
      rdt  = 1.0/dt
C
C Compute sb,hb,shbs,hbs
C
      DO k = 1, klev
      DO i = 1, klon
         zx_t = tb(i,k)
         zx_p = p(i,k)
         zx_q = shb(i,k)
           zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
           zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
           zx_qs=MIN(0.5,zx_qs)
           zcor=1./(1.-retv*zx_qs)
           zx_qs=zx_qs*zcor
           zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
         shbs(i,k) = zx_qs
         gam(i,k) = zx_gam
      ENDDO
      ENDDO
C
      DO k=limcnv,klev
         DO i=1,klon
            sb (i,k) = RCPD*tb(i,k) + gz(i,k)
            hb (i,k) = sb(i,k) + RLVTT*shb(i,k)
            hbs(i,k) = sb(i,k) + RLVTT*shbs(i,k)
         ENDDO
      ENDDO
C
C Compute sbh, shbh
C
      DO k=limcnv+1,klev
         km1 = k - 1
         DO i=1,klon
            sbh (i,k) =0.5*(sb(i,km1) + sb(i,k))
            shbh(i,k) =qhalf(shb(i,km1),shb(i,k),shbs(i,km1),shbs(i,k))
            hbh (i,k) =sbh(i,k) + RLVTT*shbh(i,k)
         ENDDO
      ENDDO
C
C Specify properties at top of model (not used, but filling anyway)
C
      DO i=1,klon
         sbh (i,limcnv) = sb(i,limcnv)
         shbh(i,limcnv) = shb(i,limcnv)
         hbh (i,limcnv) = hb(i,limcnv)
      ENDDO
C
C Zero vertically independent control, tendency & diagnostic arrays
C
      DO i=1,klon
         prec(i)  = 0.0
         dzcld(i) = 0.0
         cnb(i)   = 0
         cnt(i)   = klev+1
      ENDDO

      DO k = 1, klev
        DO i = 1,klon
         cmfdt(i,k)  = 0.
         cmfdq(i,k)  = 0.
         cmfdqr(i,k) = 0.
         cmfmc(i,k)  = 0.
         cmfsl(i,k)  = 0.
         cmflq(i,k)  = 0.
        ENDDO
      ENDDO
C
C Begin moist convective mass flux adjustment procedure.
C Formalism ensures that negative cloud liquid water can never occur
C
      DO 70 k=klev-1,limcnv+1,-1
         km1 = k - 1
         kp1 = k + 1
         DO 10 i=1,klon
            eta   (i) = 0.0
            beta  (i) = 0.0
            ds1   (i) = 0.0
            ds2   (i) = 0.0
            ds3   (i) = 0.0
            dq1   (i) = 0.0
            dq2   (i) = 0.0
            dq3   (i) = 0.0
C
C Specification of "cloud base" conditions
C
            qprime    = 0.0
            tprime    = 0.0
C
C Assign tprime within the PBL to be proportional to the quantity
C thtap (which will be bounded by tpmax), passed to this routine by 
C the PBL routine.  Don't allow perturbation to produce a dry 
C adiabatically unstable parcel.  Assign qprime within the PBL to be 
C an appropriately modified value of the quantity shp (which will be 
C bounded by shpmax) passed to this routine by the PBL routine.  The 
C quantity qprime should be less than the local saturation value 
C (qsattp=qsat[t+tprime,p]).  In both cases, thtap and shp are
C linearly reduced toward zero as the PBL top is approached.
C
            pblhgt = MAX(pblh(i),1.0)
            IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.0.0) THEN
               fac1   = MAX(0.0,1.0-gz(i,kp1)/RG/pblhgt)
               tprime = MIN(thtap(i),tpmax)*fac1
               qsattp = shbs(i,kp1) + RCPD/RLVTT*gam(i,kp1)*tprime
               shprme = MIN(MIN(shp(i),shpmax)*fac1,
     $                        MAX(qsattp-shb(i,kp1),0.0))
               qprime = MAX(qprime,shprme)
            ELSE
               tprime = 0.0
               qprime = 0.0
            ENDIF
C
C Specify "updraft" (in-cloud) thermodynamic properties
C
            sc (i)    = sb (i,kp1) + RCPD*tprime
            shc(i)    = shb(i,kp1) + qprime
            hc (i)    = sc (i    ) + RLVTT*shc(i)
            flotab(i) = hc(i) - hbs(i,k)
            dz        = dp(i,k)*RD*tb(i,k)/RG/p(i,k)
            IF (flotab(i).gt.0.0) THEN
               dzcld(i) = dzcld(i) + dz
            ELSE
               dzcld(i) = 0.0
            ENDIF
   10    CONTINUE
C
C Check on moist convective instability
C
         is = 0
         DO i = 1, klon
            IF (flotab(i).GT.0.0) THEN
               ldcum(i) = .TRUE.
               is = is + 1
            ELSE
               ldcum(i) = .FALSE.
            ENDIF
         ENDDO
C
         IF (is.EQ.0) THEN
            DO i=1,klon
               dzcld(i) = 0.0
            ENDDO
            GOTO 70
         ENDIF
C
C Current level just below top level => no overshoot
C
         IF (k.le.limcnv+1) THEN
            DO i=1,klon
            IF (ldcum(i)) THEN
               cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k))
               cldwtr(i) = MAX(0.0,cldwtr(i))
               beta(i)   = 0.0
            ENDIF
            ENDDO
            GOTO 20
         ENDIF
C
C First guess at overshoot parameter using crude buoyancy closure
C 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
C If pre-existing supersaturation in detrainment layer, beta=0
C cldwtr is temporarily equal to RLVTT*l (l=> liquid water)
C
         DO i=1,klon
         IF (ldcum(i)) THEN
            cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k))
            cldwtr(i) = MAX(0.0,cldwtr(i))
            betamx = 1.0 - c0*MAX(0.0,(dzcld(i)-dzmin))
            b1        = (hc(i) - hbs(i,km1))*dp(i,km1)
            b2        = (hc(i) - hbs(i,k  ))*dp(i,k  )
            beta(i)   = MAX(betamn,MIN(betamx, 1.0+b1/b2))
            IF (hbs(i,km1).le.hb(i,km1)) beta(i) = 0.0
         ENDIF
         ENDDO
C
C Bound maximum beta to ensure physically realistic solutions
C
C First check constrains beta so that eta remains positive
C (assuming that eta is already positive for beta equal zero)
c La premiere contrainte de beta est que le flux eta doit etre positif.
C
         DO i=1,klon
         IF (ldcum(i)) THEN
            tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
     $            - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
            tmp2 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))
            IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN
               betamx = 0.99*(tmp1/tmp2)
               beta(i) = MAX(0.0,MIN(betamx,beta(i)))
            ENDIF
C
C Second check involves supersaturation of "detrainment layer"
C small amount of supersaturation acceptable (by ssfac factor)
c La 2e contrainte est que la convection ne doit pas sursaturer
c la "detrainment layer", Neanmoins, une petite sursaturation
c est acceptee (facteur ssfac).
C
            IF (hb(i,km1).lt.hbs(i,km1)) THEN
               tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
     $               - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
               tmp1 = tmp1/dp(i,k)
               tmp2 = gam(i,km1)*(sbh(i,k)-sc(i) + cldwtr(i)) -
     $                 hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
               tmp3 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k)
               tmp4 = (dt/cats)*(hc(i)-hbs(i,k))*tmp2
     $               / (dp(i,km1)*(hbs(i,km1)-hb(i,km1))) + tmp3
               IF ((beta(i)*tmp4-tmp1).GT.0.0) THEN
                  betamx = ssfac*(tmp1/tmp4)
                  beta(i)   = MAX(0.0,MIN(betamx,beta(i)))
               ENDIF
            ELSE 
               beta(i) = 0.0
            ENDIF
C
C Third check to avoid introducing 2 delta x thermodynamic
C noise in the vertical ... constrain adjusted h (or theta e)
C so that the adjustment doesn't contribute to "kinks" in h
C
            g = MIN(0.0,hb(i,k)-hb(i,km1))
            tmp3 = (hb(i,k)-hb(i,km1)-g)*(cats/dt) / (hc(i)-hbs(i,k))
            tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
     $            - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
            tmp1 = tmp1/dp(i,k)
            tmp1 = tmp3*tmp1 + (hc(i) - hbh(i,kp1))/dp(i,k)
            tmp2 = tmp3*(1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k)
     $            + (hc(i)-hbh(i,k)-cldwtr(i))
     $             *(1.0/dp(i,k)+1.0/dp(i,kp1))
            IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN
               betamx = 0.0
               IF (tmp2.NE.0.0) betamx = tmp1/tmp2
               beta(i) = MAX(0.0,MIN(betamx,beta(i)))
            ENDIF
         ENDIF
         ENDDO
C
C Calculate mass flux required for stabilization.
C
C Ensure that the convective mass flux, eta, is positive by
C setting negative values of eta to zero..
C Ensure that estimated mass flux cannot move more than the
C minimum of total mass contained in either layer k or layer k+1.
C Also test for other pathological cases that result in non-
C physical states and adjust eta accordingly.
C
   20    CONTINUE
         DO i=1,klon
         IF (ldcum(i)) THEN
            beta(i) = MAX(0.0,beta(i))
            tmp1 = hc(i) - hbs(i,k)
            tmp2 = ((1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) -
     $               beta(i)*(1.0+gam(i,k))*(sc(i)-sbh(i,k)))/dp(i,k) -
     $              (hbh(i,kp1)-hc(i))/dp(i,kp1)
            eta(i) = tmp1/(tmp2*RG*cats)
            tmass = MIN(dp(i,k),dp(i,kp1))/RG
            IF (eta(i).GT.tmass*rdt .OR. eta(i).LE.0.0) eta(i) = 0.0
C
C Check on negative q in top layer (bound beta)
C
            IF(shc(i)-shbh(i,k).LT.0.0 .AND. beta(i)*eta(i).NE.0.0)THEN
               denom = eta(i)*RG*dt*(shc(i) - shbh(i,k))/dp(i,km1)
               beta(i) = MAX(0.0,MIN(-0.999*shb(i,km1)/denom,beta(i)))
            ENDIF
C
C Check on negative q in middle layer (zero eta)
C
            qtest1 = shb(i,k) + eta(i)*RG*dt*((shc(i) - shbh(i,kp1)) -
     $               (1.0 - beta(i))*cldwtr(i)/RLVTT -
     $               beta(i)*(shc(i) - shbh(i,k)))/dp(i,k)
            IF (qtest1.le.0.0) eta(i) = 0.0
C
C Check on negative q in lower layer (bound eta)
C
            fac1 = -(shbh(i,kp1) - shc(i))/dp(i,kp1)
            qtest2 = shb(i,kp1) - eta(i)*RG*dt*fac1
            IF (qtest2 .lt. 0.0) THEN
               eta(i) = 0.99*shb(i,kp1)/(RG*dt*fac1)
            ENDIF
         ENDIF
         ENDDO
C
C
C Calculate cloud water, rain water, and thermodynamic changes
C
         DO 30 i=1,klon
         IF (ldcum(i)) THEN
            etagdt = eta(i)*RG*dt
            cldwtr(i) = etagdt*cldwtr(i)/RLVTT/RG
            rnwtr(i) = (1.0 - beta(i))*cldwtr(i)
            ds1(i) = etagdt*(sbh(i,kp1) - sc(i))/dp(i,kp1)
            dq1(i) = etagdt*(shbh(i,kp1) - shc(i))/dp(i,kp1)
            ds2(i) = (etagdt*(sc(i) - sbh(i,kp1)) +
     $                RLVTT*RG*cldwtr(i) - beta(i)*etagdt*
     $                (sc(i) - sbh(i,k)))/dp(i,k)
            dq2(i) = (etagdt*(shc(i) - shbh(i,kp1)) -
     $                RG*rnwtr(i) - beta(i)*etagdt*
     $                (shc(i) - shbh(i,k)))/dp(i,k)
            ds3(i) = beta(i)*(etagdt*(sc(i) - sbh(i,k)) -
     $               RLVTT*RG*cldwtr(i))/dp(i,km1)
            dq3(i) = beta(i)*etagdt*(shc(i) - shbh(i,k))/dp(i,km1)
C
C Isolate convective fluxes for later diagnostics
C
            fslkp = eta(i)*(sc(i) - sbh(i,kp1))
            fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) -
     $                       RLVTT*cldwtr(i)*rdt)
            fqlkp = eta(i)*(shc(i) - shbh(i,kp1))
            fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
C
C
C Update thermodynamic profile (update sb, hb, & hbs later)
C
            tb (i,kp1) = tb(i,kp1)  + ds1(i) / RCPD
            tb (i,k  ) = tb(i,k  )  + ds2(i) / RCPD
            tb (i,km1) = tb(i,km1)  + ds3(i) / RCPD
            shb(i,kp1) = shb(i,kp1) + dq1(i)
            shb(i,k  ) = shb(i,k  ) + dq2(i)
            shb(i,km1) = shb(i,km1) + dq3(i)
            prec(i)    = prec(i)    + rnwtr(i)/rhoh2o
C
C Update diagnostic information for final budget
C Tracking temperature & specific humidity tendencies,
C rainout term, convective mass flux, convective liquid
C water static energy flux, and convective total water flux
C
            cmfdt (i,kp1) = cmfdt (i,kp1) + ds1(i)/RCPD*rdt
            cmfdt (i,k  ) = cmfdt (i,k  ) + ds2(i)/RCPD*rdt
            cmfdt (i,km1) = cmfdt (i,km1) + ds3(i)/RCPD*rdt
            cmfdq (i,kp1) = cmfdq (i,kp1) + dq1(i)*rdt
            cmfdq (i,k  ) = cmfdq (i,k  ) + dq2(i)*rdt
            cmfdq (i,km1) = cmfdq (i,km1) + dq3(i)*rdt
            cmfdqr(i,k  ) = cmfdqr(i,k  ) + (RG*rnwtr(i)/dp(i,k))*rdt
            cmfmc (i,kp1) = cmfmc (i,kp1) + eta(i)
            cmfmc (i,k  ) = cmfmc (i,k  ) + beta(i)*eta(i)
            cmfsl (i,kp1) = cmfsl (i,kp1) + fslkp
            cmfsl (i,k  ) = cmfsl (i,k  ) + fslkm
            cmflq (i,kp1) = cmflq (i,kp1) + RLVTT*fqlkp
            cmflq (i,k  ) = cmflq (i,k  ) + RLVTT*fqlkm
            qc    (i,k  ) =                (RG*rnwtr(i)/dp(i,k))*rdt
         ENDIF
   30    CONTINUE
C
C Next, convectively modify passive constituents
C
         DO 50 m=1,pcnst
         DO 40 i=1,klon
         IF (ldcum(i)) THEN
C
C If any of the reported values of the constituent is negative in
C the three adjacent levels, nothing will be done to the profile
C
            IF ((cmrb(i,kp1,m).LT.0.0) .OR.
     $          (cmrb(i,k,m).LT.0.0) .OR.
     $          (cmrb(i,km1,m).LT.0.0)) GOTO 40
C
C Specify constituent interface values (linear interpolation)
C
            cmrh(i,k  ) = 0.5*(cmrb(i,km1,m) + cmrb(i,k  ,m))
            cmrh(i,kp1) = 0.5*(cmrb(i,k  ,m) + cmrb(i,kp1,m))
C
C Specify perturbation properties of constituents in PBL
C
            pblhgt = MAX(pblh(i),1.0)
            IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.0.) THEN
               fac1 = MAX(0.0,1.0-gz(i,kp1)/RG/pblhgt)
               cmrc(i) = cmrb(i,kp1,m) + cmrp(i,m)*fac1
            ELSE
               cmrc(i) = cmrb(i,kp1,m)
            ENDIF
C
C Determine fluxes, flux divergence => changes due to convection
C Logic must be included to avoid producing negative values. A bit
C messy since there are no a priori assumptions about profiles.
C Tendency is modified (reduced) when pending disaster detected.
C
            etagdt = eta(i)*RG*dt
            botflx   = etagdt*(cmrc(i) - cmrh(i,kp1))
            topflx   = beta(i)*etagdt*(cmrc(i)-cmrh(i,k))
            dcmr1(i) = -botflx/dp(i,kp1)
            efac1    = 1.0
            efac2    = 1.0
            efac3    = 1.0
C
            IF (cmrb(i,kp1,m)+dcmr1(i) .LT. 0.0) THEN
               efac1 = MAX(tiny,ABS(cmrb(i,kp1,m)/dcmr1(i)) - eps)
            ENDIF
C
            IF (efac1.EQ.tiny .OR. efac1.GT.1.0) efac1 = 0.0
            dcmr1(i) = -efac1*botflx/dp(i,kp1)
            dcmr2(i) = (efac1*botflx - topflx)/dp(i,k)
C  
            IF (cmrb(i,k,m)+dcmr2(i) .LT. 0.0) THEN
               efac2 = MAX(tiny,ABS(cmrb(i,k  ,m)/dcmr2(i)) - eps)
            ENDIF
C
            IF (efac2.EQ.tiny .OR. efac2.GT.1.0) efac2 = 0.0
            dcmr2(i) = (efac1*botflx - efac2*topflx)/dp(i,k)
            dcmr3(i) = efac2*topflx/dp(i,km1)
C
            IF (cmrb(i,km1,m)+dcmr3(i) .LT. 0.0) THEN
               efac3 = MAX(tiny,ABS(cmrb(i,km1,m)/dcmr3(i)) - eps)
            ENDIF
C
            IF (efac3.EQ.tiny .OR. efac3.GT.1.0) efac3 = 0.0
            efac3    = MIN(efac2,efac3)
            dcmr2(i) = (efac1*botflx - efac3*topflx)/dp(i,k)
            dcmr3(i) = efac3*topflx/dp(i,km1)
C
            cmrb(i,kp1,m) = cmrb(i,kp1,m) + dcmr1(i)
            cmrb(i,k  ,m) = cmrb(i,k  ,m) + dcmr2(i)
            cmrb(i,km1,m) = cmrb(i,km1,m) + dcmr3(i)
         ENDIF
   40    CONTINUE
   50    CONTINUE              ! end of m=1,pcnst loop
C
         IF (k.EQ.limcnv+1) GOTO 60 ! on ne pourra plus glisser
c
c Dans la procedure de glissage ascendant, les variables thermo-
c dynamiques des couches k et km1 servent au calcul des couches
c superieures. Elles ont donc besoin d'une mise-a-jour.
C
         DO i = 1, klon
         IF (ldcum(i)) THEN
            zx_t = tb(i,k)
            zx_p = p(i,k)
            zx_q = shb(i,k)
              zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
              zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
              zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
              zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
              zx_qs=MIN(0.5,zx_qs)
              zcor=1./(1.-retv*zx_qs)
              zx_qs=zx_qs*zcor
              zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
            shbs(i,k) = zx_qs
            gam(i,k) = zx_gam
c
            zx_t = tb(i,km1)
            zx_p = p(i,km1)
            zx_q = shb(i,km1)
              zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
              zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
              zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
              zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
              zx_qs=MIN(0.5,zx_qs)
              zcor=1./(1.-retv*zx_qs)
              zx_qs=zx_qs*zcor
              zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
            shbs(i,km1) = zx_qs
            gam(i,km1) = zx_gam
C
            sb (i,k  ) = sb(i,k  ) + ds2(i)
            sb (i,km1) = sb(i,km1) + ds3(i)
            hb (i,k  ) = sb(i,k  ) + RLVTT*shb(i,k)
            hb (i,km1) = sb(i,km1) + RLVTT*shb(i,km1)
            hbs(i,k  ) = sb(i,k  ) + RLVTT*shbs(i,k  )
            hbs(i,km1) = sb(i,km1) + RLVTT*shbs(i,km1)
C
            sbh (i,k) = 0.5*(sb(i,k) + sb(i,km1))
            shbh(i,k) = qhalf(shb(i,km1),shb(i,k)
     $                       ,shbs(i,km1),shbs(i,k))
            hbh (i,k) = sbh(i,k) + RLVTT*shbh(i,k)
            sbh (i,km1) = 0.5*(sb(i,km1) + sb(i,k-2))
            shbh(i,km1) = qhalf(shb(i,k-2),shb(i,km1),
     $                    shbs(i,k-2),shbs(i,km1))
            hbh (i,km1) = sbh(i,km1) + RLVTT*shbh(i,km1)
         ENDIF
         ENDDO
C
C Ensure that dzcld is reset if convective mass flux zero
C specify the current vertical extent of the convective activity
C top of convective layer determined by size of overshoot param.
C
   60    CONTINUE
         DO i=1,klon
            etagt0 = eta(i).gt.0.0
            IF (.not.etagt0) dzcld(i) = 0.0
            IF (etagt0 .and. beta(i).gt.betamn) THEN
               ktp = km1
            ELSE
               ktp = k
            ENDIF
            IF (etagt0) THEN
               cnt(i) = MIN(cnt(i),ktp)
               cnb(i) = MAX(cnb(i),k)
            ENDIF
         ENDDO
   70 CONTINUE        ! end of k loop
C
C determine whether precipitation, prec, is frozen (snow) or not
C
      DO i=1,klon
         IF (tb(i,klev).LT.tmelt .AND. tb(i,klev-1).lt.tmelt) THEN
             cmfprs(i) = prec(i)*rdt
         ELSE
             cmfprt(i) = prec(i)*rdt
         ENDIF
      ENDDO
C
      RETURN  ! we're all done ... return to calling procedure
      END
