convection.F90

Go to the documentation of this file.

module convection

IMPLICIT NONE

PRIVATE
PUBLIC :: rase, rase0, acritn, sundq3_ice, dqsats_ras, dqsat_ras

CONTAINS

 SUBROUTINE rase( IDIM, IRUN, K0, ICMIN, DT,                  &
                  CONS_CP, CONS_ALHL, CONS_GRAV, CONS_RGAS,   &
                  CONS_H2OMW, CONS_AIRMW, CONS_VIREPS,        &
                  SEEDRAS, SIGE,                              &
                  KCBL, WGT0, WGT1, FRLAND, TS,               &
                  THO, QHO, UHO, VHO,                         &
                  CO_AUTO, PLE,                               &
                  CLW, FLXD, CNV_PRC3, CNV_UPDFRC,            &
                  RASPARAMS, ESTBLX                           )

 IMPLICIT NONE

 !INPUTS
 INTEGER,                        INTENT(IN   ) ::  idim, irun, k0, icmin
 REAL(8), DIMENSION (IDIM,K0+1), INTENT(IN   ) ::  ple
 REAL(8), DIMENSION (     K0+1), INTENT(IN   ) ::  sige
 REAL(8),                        INTENT(IN   ) ::  dt, cons_cp, cons_alhl, cons_grav, cons_rgas
 REAL(8),                        INTENT(IN   ) ::  cons_h2omw,cons_airmw,cons_vireps
 INTEGER, DIMENSION (IDIM)  ,    INTENT(IN   ) ::  seedras
 INTEGER, DIMENSION (IDIM  ),    INTENT(IN   ) ::  kcbl
 REAL(8), DIMENSION (IDIM  ),    INTENT(IN   ) ::  ts, frland
 REAL(8), DIMENSION (IDIM     ), INTENT(IN   ) ::  co_auto
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(IN   ) ::  wgt0,wgt1
 REAL(8), DIMENSION(:),          INTENT(IN   ) ::  rasparams
 REAL(8), DIMENSION(:),          INTENT(IN   ) ::  estblx

 !OUTPUTS
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(  OUT) ::  clw , flxd
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(  OUT) ::  cnv_prc3
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(  OUT) ::  cnv_updfrc

 !PROGNOSTIC
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(INOUT) ::  tho, qho, uho, vho

 !LOCALS
 INTEGER :: i, ic, l, kk, k

 !Parameters
 REAL(8), PARAMETER :: onepkap = 1.+ 2./7., daylen = 86400.0
 REAL(8), PARAMETER :: rhmax   = 0.9999
 REAL(8), PARAMETER :: cbl_qpert = 0.0, cbl_tpert = 1.0
 REAL(8), PARAMETER :: cbl_tpert_mxocn = 2.0, cbl_tpert_mxlnd = 4.0

 !Constants
 REAL(8) :: grav, cp, alhl, cpbg, alhi, cpi, gravi, ddt, lbcp 
   
 !Rasparams
 REAL(8) :: fricfac, cli_crit, rasal1, rasal2 
 REAL(8) :: friclambda
 REAL(8) :: sdqv2, sdqv3, sdqvt1
 REAL(8) :: acritfac,  pblfrac, autorampb
 REAL(8) :: maxdallowed, rhmn, rhmx

 REAL(8) :: mxdiam
 REAL(8) :: tx2, tx3, akm, acr, alm, tth, qqh, dqx
 REAL(8) :: wfn, tem, trg, trgexp, evp, wlq, qcc
 REAL(8) :: cli, te_a, c00_x, cli_crit_x, toki
 REAL(8) :: dt_lyr, rate, cvw_x, closs, f2, f3, f4
 REAL(8) :: wght0, prcbl, rndu
 REAL(8) :: lambda_min, lambda_max
 REAL(8) :: tpert, qpert
 REAL(8) :: uht, vht

 REAL(8), DIMENSION(K0) :: poi_sv, qoi_sv, uoi_sv, voi_sv
 REAL(8), DIMENSION(K0) :: poi, qoi, uoi, voi, dqq, bet, gam, cll
 REAL(8), DIMENSION(K0) :: poi_c, qoi_c
 REAL(8), DIMENSION(K0) :: prh,  pri,  ght, dpt, dpb, pki
 REAL(8), DIMENSION(K0) :: ucu, vcu
 REAL(8), DIMENSION(K0) :: cln, rns, pol
 REAL(8), DIMENSION(K0) :: qst, ssl,  rmf, rnn, rn1, rmfc, rmfp
 REAL(8), DIMENSION(K0) :: gms, eta, gmh, eht,  gm1, hcc, rmfd
 REAL(8), DIMENSION(K0) :: hol, hst, qol, zol, hcld, cll0,cllx,clli
 REAL(8), DIMENSION(K0) :: bke , cvw, updfrc
 REAL(8), DIMENSION(K0) :: rasal, updfrp,bk2,dll0,dllx
 REAL(8), DIMENSION(K0) :: wght, massf
 REAL(8), DIMENSION(K0) :: qss, dqs, pf, pk, tempf, zlo
 REAL(8), DIMENSION(K0+1) :: prj, prs, qht, sht ,zet, zle, pke

 !Initialize Local Arrays
 poi = 0.0
 qoi = 0.0
 uoi = 0.0
 voi = 0.0
 dqq = 0.0
 bet = 0.0
 gam = 0.0
 poi_c = 0.0
 qoi_c = 0.0
 prh = 0.0
 pri = 0.0
 ght = 0.0
 dpt = 0.0
 dpb = 0.0
 pki = 0.0
 ucu = 0.0
 vcu = 0.0
 cln = 0.0
 pol = 0.0
 qst = 0.0
 ssl = 0.0
 rmf = 0.0
 rnn = 0.0
 rn1 = 0.0
 gms = 0.0
 eta = 0.0
 gmh = 0.0
 eht = 0.0
 gm1 = 0.0
 hcc = 0.0
 hol = 0.0
 hst = 0.0
 qol = 0.0
 zol = 0.0
 hcld = 0.0
 bke = 0.0
 cvw = 0.0
 updfrc = 0.0
 rasal = 0.0
 updfrp = 0.0
 bk2 = 0.0
 wght = 0.0
 massf = 0.0
 qss = 0.0
 dqs = 0.0
 pf = 0.0
 pk = 0.0
 tempf = 0.0
 zlo = 0.0
 prj = 0.0
 prs = 0.0
 qht = 0.0
 sht = 0.0
 zet = 0.0
 zle = 0.0
 pke = 0.0

 !Initialize Outputs
 cnv_prc3  =0.
 cnv_updfrc=0. 
 flxd = 0.
 clw = 0.

 fricfac      = rasparams(1)     !  ---  1
 cli_crit     = rasparams(4)     !  ---  4
 rasal1       = rasparams(5)     !  ---  5
 rasal2       = rasparams(6)     !  ---  6
 friclambda   = rasparams(11)    !  --- 11
 sdqv2        = rasparams(14)    !  --- 14
 sdqv3        = rasparams(15)    !  --- 15
 sdqvt1       = rasparams(16)    !  --- 16
 acritfac     = rasparams(17)    !  --- 17
 pblfrac      = rasparams(20)    !  --- 20
 autorampb    = rasparams(21)    !  --- 21
 rhmn         = rasparams(24)    !  --- 24
 maxdallowed  = rasparams(23)    !  --- 24
 rhmx         = rasparams(25)    !  --- 25

 grav  = cons_grav
 alhl  = cons_alhl
 cp    = cons_cp
 cpi   = 1.0/cp      
 alhi  = 1.0/alhl
 gravi = 1.0/grav
 cpbg  = cp*gravi
 ddt   = daylen/dt
 lbcp  = alhl*cpi
 
i = 1

    !CALL FINDBASE
    k = kcbl(i)

    IF ( k > 0 ) THEN 

       !Get saturation specific humidity and gradient wrt to T
       pke  = (ple(i,:)/1000.)**(cons_rgas/cons_cp)
       pf   = 0.5*(ple(i,1:k0) +  ple(i,2:k0+1  ) )
       pk   = (pf/1000.)**(cons_rgas/cons_cp)
       tempf = tho(i,:)*pk

       zle = 0.0
       zlo = 0.0
       zle(k0+1) = 0.
       do l=k0,1,-1
          zle(l) = tho(i,l)   * (1.+cons_vireps*qho(i,l))
          zlo(l) = zle(l+1) + (cons_cp/cons_grav)*( pke(l+1)-pk(l  ) ) * zle(l)
          zle(l) = zlo(l)   + (cons_cp/cons_grav)*( pk(l)  -pke(l)   ) * zle(l)
       end do

       tpert  = cbl_tpert * ( ts(i) - ( tempf(k0)+ cons_grav*zlo(k0)/cons_cp )  ) 
       qpert  = cbl_qpert !* ( QSSFC - Q(:,:,K0) ) [CBL_QPERT = 0.0]
       tpert  = max( tpert , 0.0 )
       qpert  = max( qpert , 0.0 )
      
       if (frland(i) < 0.1) then
          tpert = min( tpert , cbl_tpert_mxocn ) ! ocean
       else
          tpert = min( tpert , cbl_tpert_mxlnd ) ! land
       endif

       call dqsat_ras(dqs,qss,tempf,pf,k0,estblx,cons_h2omw,cons_airmw)

       do kk=icmin,k+1
          prj(kk) = pke(kk)
       enddo

       prs(icmin:k0+1) = ple(i,icmin:k0+1)
       poi(icmin:k)   = tho(i,icmin:k)
       qoi(icmin:k)   = qho(i,icmin:k)
       uoi(icmin:k)   = uho(i,icmin:k)
       voi(icmin:k)   = vho(i,icmin:k)

       qst(icmin:k) = qss(icmin:k)
       dqq(icmin:k) = dqs(icmin:k)

       !Mass fraction of each layer below cloud base
       massf(:) = wgt0(i,:)

       !RESET PRESSURE at bottom edge of CBL 
       prcbl = prs(k)
       do l= k,k0
          prcbl = prcbl + massf(l)*( prs(l+1)-prs(l) )
       end do
       prs(k+1) = prcbl
       prj(k+1) = (prs(k+1)/1000.)**(cons_rgas/cons_cp)

       DO l=k,icmin,-1
          pol(l)  = 0.5*(prs(l)+prs(l+1))
          prh(l)  = (prs(l+1)*prj(l+1)-prs(l)*prj(l)) / (onepkap*(prs(l+1)-prs(l)))
          pki(l)  = 1.0 / prh(l)
          dpt(l)  = prh(l  ) - prj(l)
          dpb(l)  = prj(l+1) - prh(l)
          pri(l)  = .01 / (prs(l+1)-prs(l))
       ENDDO

       !RECALCULATE PROFILE QUAN. IN LOWEST STRAPPED LAYER
       if ( k<=k0) then
          poi(k) = 0.
          qoi(k) = 0.
          uoi(k) = 0.
          voi(k) = 0.

          !SPECIFY WEIGHTS GIVEN TO EACH LAYER WITHIN SUBCLOUD "SUPERLAYER"
          wght = 0.
          DO l=k,k0
             wght(l)   = massf(l) * ( ple(i,l+1) - ple(i,l) ) / ( prs(k+1)   - prs(k)  )
          END DO      

          DO l=k,k0
             poi(k) = poi(k) + wght(l)*tho(i,l)
             qoi(k) = qoi(k) + wght(l)*qho(i,l)
             uoi(k) = uoi(k) + wght(l)*uho(i,l)
             voi(k) = voi(k) + wght(l)*vho(i,l)
          ENDDO

          call dqsats_ras(dqq(k), qst(k), poi(k)*prh(k),pol(k),estblx,cons_h2omw,cons_airmw)

       endif

       rndu = max( seedras(i)/1000000., 1e-6 )
       mxdiam = maxdallowed*( rndu**(-1./2.) )

       DO l=k,icmin,-1
          bet(l)  = dqq(l)*pki(l)  !*
          gam(l)  = pki(l)/(1.0+lbcp*dqq(l)) !*
          IF (l<k) THEN
             ght(l+1) = gam(l)*dpb(l) + gam(l+1)*dpt(l+1)
             gm1(l+1) = 0.5*lbcp*(dqq(l  )/(alhl*(1.0+lbcp*dqq(l  ))) + dqq(l+1)/(alhl*(1.0+lbcp*dqq(l+1))) )
          ENDIF
       ENDDO

       rns  = 0.
       cll  = 0.
       rmf  = 0.
       rmfd = 0.
       rmfc = 0.
       rmfp = 0.
       cll0 = 0.
       dll0 = 0.
       cllx = 0.
       dllx = 0.
       clli = 0.

       poi_sv = poi
       qoi_sv = qoi
       uoi_sv = uoi
       voi_sv = voi

       cvw     = 0.0
       updfrc  = 0.0
       updfrp  = 0.0

       hol=0.            ! HOL initialized here in order not to confuse Valgrind debugger
       zet(k+1) = 0
       sht(k+1) = cp*poi(k)*prj(k+1)
       DO l=k,icmin,-1
          qol(l)  = min(qst(l)*rhmax,qoi(l))
          qol(l)  = max( 0.000, qol(l) )     ! GUARDIAN/NEG.PREC.
          ssl(l)  = cp*prj(l+1)*poi(l) + grav*zet(l+1)
          hol(l)  = ssl(l) + qol(l)*alhl
          hst(l)  = ssl(l) + qst(l)*alhl
          tem     = poi(l)*(prj(l+1)-prj(l))*cpbg
          zet(l)  = zet(l+1) + tem
          zol(l)  = zet(l+1) + (prj(l+1)-prh(l))*poi(l)*cpbg
       ENDDO

       DO ic =  k,icmin+1,-1

          ucu(icmin:) = 0.
          vcu(icmin:) = 0.
              
          alm   = 0.
          trg   = min(1.,(qoi(k)/qst(k)-rhmn)/(rhmx-rhmn))

          f4  = min(   1.0,  max( 0.0 , (autorampb-sige(ic))/0.2 )  )  
     
          IF (trg <= 1.0e-5) THEN
             cycle !================>>
          ENDIF

          !RECOMPUTE SOUNDING UP TO DETRAINMENT LEVEL
          poi_c = poi
          qoi_c = qoi
          poi_c(k) =  poi_c(k) + tpert
          qoi_c(k) =  qoi_c(k) + qpert

          zet(k+1) = 0.
          sht(k+1) = cp*poi_c(k)*prj(k+1)
          DO l=k,ic,-1
             qol(l)  = min(qst(l)*rhmax,qoi_c(l))
             qol(l)  = max( 0.000, qol(l) )     ! GUARDIAN/NEG.PREC.
             ssl(l)  = cp*prj(l+1)*poi_c(l) + grav*zet(l+1)
             hol(l)  = ssl(l) + qol(l)*alhl
             hst(l)  = ssl(l) + qst(l)*alhl
             tem     = poi_c(l)*(prj(l+1)-prj(l))*cpbg
             zet(l)  = zet(l+1) + tem
             zol(l)  = zet(l+1) + (prj(l+1)-prh(l))*poi_c(l)*cpbg
          ENDDO

          DO l=ic+1,k
             tem  = (prj(l)-prh(l-1))/(prh(l)-prh(l-1))
             sht(l)  = ssl(l-1) + tem*(ssl(l)-ssl(l-1)) 
             qht(l)  = .5*(qol(l)+qol(l-1))
          ENDDO

          !CALCULATE LAMBDA, ETA, AND WORKFUNCTION
          lambda_min = .2/mxdiam
          lambda_max = .2/  200. 

          IF (hol(k) <= hst(ic)) THEN
             cycle !================>>
          ENDIF

          !LAMBDA CALCULATION: MS-A18

          tem  =       (hst(ic)-hol(ic))*(zol(ic)-zet(ic+1)) 
          DO l=ic+1,k-1
             tem = tem + (hst(ic)-hol(l ))*(zet(l )-zet(l +1))
          ENDDO

          IF (tem <= 0.0) THEN
             cycle !================>>
          ENDIF

          alm = (hol(k)-hst(ic)) / tem

          IF (alm > lambda_max) THEN
             cycle !================>>
          ENDIF

          toki=1.0

          IF (alm < lambda_min) THEN
             toki = ( alm/lambda_min )**2
          ENDIF

          !ETA CALCULATION: MS-A2
          DO l=ic+1,k
             eta(l) = 1.0 + alm * (zet(l )-zet(k))
          ENDDO
          eta(ic) = 1.0 + alm * (zol(ic)-zet(k))

          !WORKFUNCTION CALCULATION:  MS-A22
          wfn     = 0.0
          hcc(k)  = hol(k)
          DO l=k-1,ic+1,-1
             hcc(l) = hcc(l+1) + (eta(l) - eta(l+1))*hol(l)
             tem    = hcc(l+1)*dpb(l) + hcc(l)*dpt(l)
             eht(l) = eta(l+1)*dpb(l) + eta(l)*dpt(l)
             wfn    = wfn + (tem - eht(l)*hst(l))*gam(l)
          ENDDO
          hcc(ic) = hst(ic)*eta(ic)
          wfn     = wfn + (hcc(ic+1)-hst(ic)*eta(ic+1))*gam(ic)*dpb(ic)

          !VERTICAL VELOCITY/KE CALCULATION (ADDED 12/2001 JTB)
          bk2(k)   = 0.0
          bke(k)   = 0.0
          hcld(k)  = hol(k)
          DO l=k-1,ic,-1
             hcld(l) = ( eta(l+1)*hcld(l+1) + (eta(l) - eta(l+1))*hol(l) ) / eta(l)
             tem     = (hcld(l)-hst(l) )*(zet(l)-zet(l+1))/ (1.0+lbcp*dqq(l))
             bke(l)  = bke(l+1) + grav * tem / ( cp*prj(l+1)*poi(l) )
             bk2(l)  = bk2(l+1) + grav * max(tem,0.0) / ( cp*prj(l+1)*poi(l) )
             cvw(l) = sqrt(  2.0* max( bk2(l) , 0.0 )  )    
          ENDDO

          !ALPHA CALCULATION 
          IF ( zet(ic) <  2000. ) THEN
             rasal(ic) = rasal1
          ENDIF
          IF ( zet(ic)  >= 2000. ) THEN 
             rasal(ic) = rasal1 + (rasal2-rasal1)*(zet(ic) - 2000.)/8000.
          ENDIF
          rasal(ic) = min( rasal(ic) , 1.0e5 )
    
          rasal(ic) = dt / rasal(ic)

          DO l = ic,k 
             cvw(l) = max(  cvw(l) , 1.00 )
          ENDDO

          CALL acritn(pol(ic), prs(k), acr, acritfac)

          IF (wfn <= acr) THEN
             cycle !================>>
          ENDIF

          wlq     = qol(k)
          uht     = uoi(k)
          vht     = voi(k)
          rnn(k)  = 0.
          cll0(k) = 0.

          DO l=k-1,ic,-1
             tem   = eta(l) - eta(l+1)
             wlq   = wlq + tem * qol(l)
             uht   = uht + tem * uoi(l)
             vht   = vht + tem * voi(l)

             IF (l>ic) THEN
                tx2   = 0.5*(qst(l)+qst(l-1))*eta(l)
                tx3   = 0.5*(hst(l)+hst(l-1))*eta(l)
                qcc   = tx2 + gm1(l)*(hcc(l)-tx3)
                cll0(l)  = (wlq-qcc)
             ELSE
                cll0(l)   = (wlq-qst(ic)*eta(ic))
             ENDIF

             cll0(l)    = max( cll0(l) , 0.00 )

             cli  = cll0(l) / eta(l)
             te_a = poi(l)*prh(l)

             CALL sundq3_ice( te_a,  sdqv2, sdqv3, sdqvt1, f2 , f3 )

             c00_x  =  co_auto(i) * f2 * f3  * f4 
             cli_crit_x =  cli_crit / ( f2 * f3 )
       
             rate = c00_x * ( 1.0 - exp( -(cli)**2 / cli_crit_x**2 ) )

             cvw_x     = max( cvw(l) , 1.00 )  

             dt_lyr  = ( zet(l)-zet(l+1) )/cvw_x 

             closs   = cll0(l) * rate * dt_lyr
             closs   = min( closs , cll0(l) )

             cll0(l) = cll0(l) - closs
             dll0(l) = closs

             IF (closs>0.) then
                wlq = wlq - closs
                rnn(l) = closs 
             else
                rnn(l) = 0.
             ENDIF

          ENDDO

          wlq = wlq - qst(ic)*eta(ic)

          !CALCULATE GAMMAS AND KERNEL
          gms(k) =          (sht(k)-ssl(k))*pri(k)
          gmh(k) = gms(k) + (qht(k)-qol(k))*pri(k)*alhl
          akm    = gmh(k)*gam(k-1)*dpb(k-1)

          tx2     = gmh(k)
          DO l=k-1,ic+1,-1
             gms(l) = ( eta(l  )*(sht(l)-ssl(l  )) + eta(l+1)*(ssl(l)-sht(l+1)) )     *pri(l)
             gmh(l) = gms(l) + ( eta(l  )*(qht(l)-qol(l  )) + eta(l+1)*(qol(l)-qht(l+1)) )*alhl*pri(l)
             tx2 = tx2 + (eta(l) - eta(l+1)) * gmh(l)
             akm = akm - gms(l)*eht(l)*pki(l) + tx2*ght(l)
          ENDDO

          gms(ic) = eta(ic+1)*(ssl(ic)-sht(ic+1))*pri(ic)
          akm     = akm - gms(ic)*eta(ic+1)*dpb(ic)*pki(ic)

          gmh(ic) =   gms(ic) + ( eta(ic+1)*(qol(ic)-qht(ic+1))*alhl + eta(ic)*(hst(ic)-hol(ic)))*pri(ic)

          !CLOUD BASE MASS FLUX
          IF (akm >= 0.0 .OR. wlq < 0.0)  THEN
             cycle !================>>
          ENDIF

          wfn = - (wfn-acr)/akm
          wfn = min( ( rasal(ic)*trg*toki )*wfn  ,   (prs(k+1)-prs(k) )*(100.*pblfrac))
     
          !CUMULATIVE PRECIP AND CLOUD-BASE MASS FLUX FOR OUTPUT
          tem      = wfn*gravi
          cll(ic) = cll(ic) + wlq*tem
          rmf(ic) = rmf(ic) +     tem
          rmfd(ic) = rmfd(ic) +     tem * eta(ic)

          DO l=ic+1,k
             rmfp(l) = tem * eta(l)
             rmfc(l) = rmfc(l)  +  rmfp(l)
           
             dllx(l) = dllx(l)  +  tem*dll0(l)        

             IF ( cvw(l) > 0.0 ) THEN
                updfrp(l) = rmfp(l) * (ddt/daylen) * 1000. / ( cvw(l) * prs(l) )
             ELSE 
                updfrp(l) = 0.0       
             ENDIF

             clli(l) = cll0(l)/eta(l)

             updfrc(l) =  updfrc(l) +  updfrp(l)      
 
          ENDDO

          !THETA AND Q CHANGE DUE TO CLOUD TYPE IC
          DO l=ic,k
             rns(l) = rns(l) + rnn(l)*tem
             gmh(l) = gmh(l) * wfn
             gms(l) = gms(l) * wfn
             qoi(l) = qoi(l) + (gmh(l) - gms(l)) * alhi
             poi(l) = poi(l) + gms(l)*pki(l)*cpi
             qst(l) = qst(l) + gms(l)*bet(l)*cpi
          ENDDO

          wfn     = wfn*0.5 *1.0           !*FRICFAC*0.5

          !CUMULUS FRICTION
          IF (fricfac <= 0.0) THEN
             cycle !================>>
          ENDIF

          wfn     = wfn*fricfac*exp( -alm / friclambda )
          tem     = wfn*pri(k)

          ucu(k)  = ucu(k) + tem * (uoi(k-1) - uoi(k))
          vcu(k)  = vcu(k) + tem * (voi(k-1) - voi(k))

          DO l=k-1,ic+1,-1
           tem    = wfn*pri(l)
           ucu(l) = ucu(l) + tem * ( (uoi(l-1) - uoi(l)) * eta(l) + (uoi(l) - uoi(l+1)) * eta(l+1) )
           vcu(l) = vcu(l) + tem * ( (voi(l-1) - voi(l)) * eta(l) + (voi(l) - voi(l+1)) * eta(l+1) )
          ENDDO

          tem     = wfn*pri(ic)
          ucu(ic) = ucu(ic) + (2.*(uht-uoi(ic)*(eta(ic)-eta(ic+1))) - (uoi(ic)+uoi(ic+1))*eta(ic+1))*tem
          vcu(ic) = vcu(ic) + (2.*(vht-voi(ic)*(eta(ic)-eta(ic+1))) - (voi(ic)+voi(ic+1))*eta(ic+1))*tem

          DO l=ic,k
             uoi(l) = uoi(l) + ucu(l)
             voi(l) = voi(l) + vcu(l)
          ENDDO

       ENDDO !CLOUD LOOP

       IF ( sum( rmf(icmin:k) ) > 0.0 ) THEN

          DO l=icmin,k
             tem    = pri(l)*grav
             cnv_prc3(i,l) = rns(l)*tem     
          ENDDO

          tho(i,icmin:k-1) = poi(icmin:k-1)
          qho(i,icmin:k-1) = qoi(icmin:k-1)
          uho(i,icmin:k-1) = uoi(icmin:k-1)
          vho(i,icmin:k-1) = voi(icmin:k-1)
          cnv_updfrc(i,icmin:k-1)   =  updfrc(icmin:k-1)

          !De-strap tendencies from RAS
          wght   = wgt1(i,:)

          !Scale properly by layer masses
          wght0 = 0.
          DO l=k,k0 
             wght0 = wght0 + wght(l)* ( ple(i,l+1) - ple(i,l) )
          END DO
         
          wght0 = ( prs(k+1)   - prs(k)  )/wght0
          wght  = wght0 * wght

          DO l=k,k0 
             tho(i,l) =  tho(i,l) + wght(l)*(poi(k) - poi_sv(k))
             qho(i,l) =  qho(i,l) + wght(l)*(qoi(k) - qoi_sv(k))
             uho(i,l) =  uho(i,l) + wght(l)*(uoi(k) - uoi_sv(k))
             vho(i,l) =  vho(i,l) + wght(l)*(voi(k) - voi_sv(k))
          END DO

          flxd(i,icmin:k) = rmfd(icmin:k) * ddt/daylen
          clw(i,icmin:k) = cll(icmin:k) * ddt/daylen

          flxd(i,1:icmin-1) = 0.
          clw(i,1:icmin-1) = 0.

          IF ( k < k0 ) THEN 
             flxd(i,k:k0) = 0.
             clw(i,k:k0) = 0.
          END IF

       ELSE

          flxd(i,:) = 0.
          clw(i,:) = 0.

       ENDIF 

    ELSE 
     
       flxd(i,:) = 0.
       clw(i,:) = 0.

    ENDIF

END SUBROUTINE rase

SUBROUTINE acritn( PL, PLB, ACR, ACRITFAC )
      
 IMPLICIT NONE

 REAL(8), INTENT(IN ) :: pl, plb, acritfac
 REAL(8), INTENT(OUT) :: acr

 INTEGER :: iwk
      
 REAL(8), PARAMETER :: ph(15)=(/150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, &
                             550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0/)

 REAL(8), PARAMETER :: a(15)=(/ 1.6851, 1.1686, 0.7663, 0.5255, 0.4100, 0.3677, &
                             0.3151, 0.2216, 0.1521, 0.1082, 0.0750, 0.0664, &
                             0.0553, 0.0445, 0.0633  /) 

 iwk = int(pl * 0.02 - 0.999999999)

 IF (iwk .GT. 1 .AND. iwk .LE. 15) THEN
    acr = a(iwk-1) + (pl-ph(iwk-1))*.02*(a(iwk)-a(iwk-1))
 ELSEIF(iwk > 15) THEN
    acr = a(15)
 ELSE
    acr = a(1)
 ENDIF

 acr = acritfac  * acr * (plb - pl)

ENDSUBROUTINE acritn

SUBROUTINE sundq3_ice( TEMP,RATE2,RATE3,TE1, F2, F3)

IMPLICIT NONE

REAL(8), INTENT( IN) :: temp,rate2,rate3,te1
REAL(8), INTENT(OUT) :: f2, f3

REAL(8) :: xx, yy,te0,te2,jump1  !,RATE2,RATE3,TE1

 te0=273.
 te2=200.
 jump1=  (rate2-1.0) / ( ( te0-te1 )**0.333 ) 

 ! Ice - phase treatment
 IF ( temp .GE. te0 ) THEN 
    f2   = 1.0
    f3   = 1.0
 ENDIF

 IF ( ( temp .GE. te1 ) .AND. ( temp .LT. te0 ) )THEN 
    f2   = 1.0 + jump1 * (( te0 - temp )**0.3333)
    f3   = 1.0
 ENDIF

 IF ( temp .LT. te1 ) THEN 
    f2   = rate2 + (rate3-rate2)*(te1-temp)/(te1-te2)
    f3   = 1.0
 ENDIF

 IF ( f2 .GT. 27.0 ) THEN
    f2 = 27.0
 ENDIF

endsubroutine sundq3_ice

subroutine dqsat_ras(DQSi,QSSi,TEMP,PLO,lm,ESTBLX,CONS_H2OMW,CONS_AIRMW)
!COMPUTES SATURATION VAPOUR PRESSURE QSSi AND GRADIENT w.r.t TEMPERATURE DQSi.
!INPUTS ARE TEMPERATURE AND PLO (PRESSURE AT T-LEVELS)
!VALES ARE COMPUTED FROM LOOK-UP TALBE (PIECEWISE LINEAR)

 IMPLICIT NONE

 !Inputs
 INTEGER :: lm
 REAL(8), dimension(lm) :: temp, plo
 REAL(8) :: estblx(:)
 REAL(8) :: cons_h2omw, cons_airmw

 !Outputs
 REAL(8), dimension(lm) :: dqsi, qssi

 !Locals
 REAL(8), parameter :: max_mixing_ratio = 1.0
 REAL(8) :: esfac

 INTEGER :: k

 REAL(8) :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
 INTEGER :: it

 INTEGER, parameter :: degsubs    =  100
 REAL(8), parameter :: tmintbl    =  150.0, tmaxtbl = 333.0
 INTEGER, parameter :: tablesize  =  nint(tmaxtbl-tmintbl)*degsubs + 1

 esfac = cons_h2omw/cons_airmw

 do k=1,lm

    tl = temp(k)
    pl = plo(k)

    pp = pl*100.0

    if (tl<=tmintbl) then
       ti = tmintbl
    elseif(tl>=tmaxtbl-.001) then
       ti = tmaxtbl-.001
    else
       ti = tl
    end if

    tt = (ti - tmintbl)*degsubs+1
    it = int(tt)

    dqq =  estblx(it+1) - estblx(it)
    qq  =  (tt-it)*dqq + estblx(it)

    if (pp <= qq) then
       qsat = max_mixing_ratio
       dqsat = 0.0
    else
       dd = 1.0/(pp - (1.0-esfac)*qq)
       qsat = esfac*qq*dd
       dqsat = (esfac*degsubs)*dqq*pp*(dd*dd)
    end if

    dqsi(k) = dqsat
    qssi(k) = qsat

 end do

end subroutine dqsat_ras

subroutine dqsats_ras(DQSi,QSSi,TEMP,PLO,ESTBLX,CONS_H2OMW,CONS_AIRMW)
!COMPUTES SATURATION VAPOUR PRESSURE QSSi AND GRADIENT w.r.t TEMPERATURE DQSi.
!INPUTS ARE TEMPERATURE AND PLO (PRESSURE AT T-LEVELS)
!VALES ARE COMPUTED FROM LOOK-UP TALBE (PIECEWISE LINEAR)

 IMPLICIT NONE

 !Inputs
 REAL(8) :: temp, plo
 REAL(8) :: estblx(:)
 REAL(8) :: cons_h2omw, cons_airmw

 !Outputs
 REAL(8) :: dqsi, qssi

 !Locals
 REAL(8), parameter :: max_mixing_ratio = 1.0
 REAL(8) :: esfac

 REAL(8) :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
 INTEGER :: it

 INTEGER, parameter :: degsubs    =  100
 REAL(8), parameter :: tmintbl    =  150.0, tmaxtbl = 333.0
 INTEGER, parameter :: tablesize  =  nint(tmaxtbl-tmintbl)*degsubs + 1

 esfac = cons_h2omw/cons_airmw

 tl = temp
 pl = plo

 pp = pl*100.0
 
 if (tl<=tmintbl) then
    ti = tmintbl
 elseif(tl>=tmaxtbl-.001) then
    ti = tmaxtbl-.001
 else
    ti = tl
 end if

 tt = (ti - tmintbl)*degsubs+1
 it = int(tt)

 dqq =  estblx(it+1) - estblx(it)
 qq  =  (tt-it)*dqq + estblx(it)

 if (pp <= qq) then
    qsat = max_mixing_ratio
    dqsat = 0.0
 else
    dd = 1.0/(pp - (1.0-esfac)*qq)
    qsat = esfac*qq*dd
    dqsat = (esfac*degsubs)*dqq*pp*(dd*dd)
 end if

 dqsi = dqsat
 qssi = qsat

end subroutine dqsats_ras

 SUBROUTINE rase0(IDIM, IRUN, K0, ICMIN, DT,                  &
                  CONS_CP, CONS_ALHL, CONS_GRAV, CONS_RGAS,   &
                  CONS_H2OMW, CONS_AIRMW, CONS_VIREPS,        &
                  SEEDRAS, SIGE,                              &
                  KCBL, WGT0, WGT1, FRLAND, TS,               &
                  THO, QHO,                                   &
                  CO_AUTO, PLE,                               &
                  CLW, FLXD, CNV_PRC3, CNV_UPDFRC,            &
                  RASPARAMS, ESTBLX                           )

 IMPLICIT NONE

 !INPUTS
 INTEGER,                        INTENT(IN   ) ::  idim, irun, k0, icmin
 REAL(8), DIMENSION (IDIM,K0+1), INTENT(IN   ) ::  ple
 REAL(8), DIMENSION (     K0+1), INTENT(IN   ) ::  sige
 REAL(8),                        INTENT(IN   ) ::  dt, cons_cp, cons_alhl, cons_grav, cons_rgas
 REAL(8),                        INTENT(IN   ) ::  cons_h2omw,cons_airmw,cons_vireps
 INTEGER, DIMENSION (IDIM)  ,    INTENT(IN   ) ::  seedras
 INTEGER, DIMENSION (IDIM  ),    INTENT(IN   ) ::  kcbl
 REAL(8), DIMENSION (IDIM  ),    INTENT(IN   ) ::  ts, frland
 REAL(8), DIMENSION (IDIM     ), INTENT(IN   ) ::  co_auto
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(IN   ) ::  wgt0,wgt1
 REAL(8), DIMENSION(:),          INTENT(IN   ) ::  rasparams
 REAL(8), DIMENSION(:),          INTENT(IN   ) ::  estblx

 !OUTPUTS
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(  OUT) ::  clw , flxd
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(  OUT) ::  cnv_prc3
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(  OUT) ::  cnv_updfrc

 !PROGNOSTIC
 REAL(8), DIMENSION (IDIM,K0  ), INTENT(INOUT) ::  tho, qho

 !LOCALS
 INTEGER :: i, ic, l, kk, k

 !Parameters
 REAL(8), PARAMETER :: onepkap = 1.+ 2./7., daylen = 86400.0
 REAL(8), PARAMETER :: rhmax   = 0.9999
 REAL(8), PARAMETER :: cbl_qpert = 0.0, cbl_tpert = 1.0
 REAL(8), PARAMETER :: cbl_tpert_mxocn = 2.0, cbl_tpert_mxlnd = 4.0

 !Constants
 REAL(8) :: grav, cp, alhl, cpbg, alhi, cpi, gravi, ddt, lbcp 
   
 !Rasparams
 REAL(8) :: fricfac, cli_crit, rasal1, rasal2 
 REAL(8) :: friclambda
 REAL(8) :: sdqv2, sdqv3, sdqvt1
 REAL(8) :: acritfac,  pblfrac, autorampb
 REAL(8) :: maxdallowed, rhmn, rhmx

 REAL(8) :: mxdiam
 REAL(8) :: tx2, tx3, akm, acr, alm, tth, qqh, dqx
 REAL(8) :: wfn, tem, trg, trgexp, evp, wlq, qcc
 REAL(8) :: cli, te_a, c00_x, cli_crit_x, toki
 REAL(8) :: dt_lyr, rate, cvw_x, closs, f2, f3, f4
 REAL(8) :: wght0, prcbl, rndu
 REAL(8) :: lambda_min, lambda_max
 REAL(8) :: tpert,qpert

 REAL(8), DIMENSION(K0) :: poi_sv, qoi_sv
 REAL(8), DIMENSION(K0) :: poi, qoi, dqq, bet, gam, cll
 REAL(8), DIMENSION(K0) :: poi_c, qoi_c
 REAL(8), DIMENSION(K0) :: prh,  pri,  ght, dpt, dpb, pki
 REAL(8), DIMENSION(K0) :: cln, rns, pol,dm
 REAL(8), DIMENSION(K0) :: qst, ssl,  rmf, rnn, rn1, rmfc, rmfp
 REAL(8), DIMENSION(K0) :: gms, eta, gmh, eht,  gm1, hcc, rmfd
 REAL(8), DIMENSION(K0) :: hol, hst, qol, zol, hcld, cll0,cllx,clli
 REAL(8), DIMENSION(K0) :: bke , cvw, updfrc
 REAL(8), DIMENSION(K0) :: rasal, updfrp,bk2,dll0,dllx
 REAL(8), DIMENSION(K0) :: wght, massf
 REAL(8), DIMENSION(K0) :: qss, dqs, pf, pk, tempf, zlo
 REAL(8), DIMENSION(K0+1) :: prj, prs, qht, sht ,zet, zle, pke

 cnv_prc3  =0.
 cnv_updfrc=0. 
 flxd = 0.
 clw = 0.

 fricfac      = rasparams(1)     !  ---  1
 cli_crit     = rasparams(4)     !  ---  4
 rasal1       = rasparams(5)     !  ---  5
 rasal2       = rasparams(6)     !  ---  6
 friclambda   = rasparams(11)    !  --- 11
 sdqv2        = rasparams(14)    !  --- 14
 sdqv3        = rasparams(15)    !  --- 15
 sdqvt1       = rasparams(16)    !  --- 16
 acritfac     = rasparams(17)    !  --- 17
 pblfrac      = rasparams(20)    !  --- 20
 autorampb    = rasparams(21)    !  --- 21
 rhmn         = rasparams(24)    !  --- 24
 maxdallowed  = rasparams(23)    !  --- 24
 rhmx         = rasparams(25)    !  --- 25

 grav  = cons_grav
 alhl  = cons_alhl
 cp    = cons_cp
 cpi   = 1.0/cp      
 alhi  = 1.0/alhl
 gravi = 1.0/grav
 cpbg  = cp*gravi
 ddt   = daylen/dt
 lbcp  = alhl*cpi
 
 DO i = 1,irun

    !CALL FINDBASE
    k = kcbl(i)

    IF ( k > 0 ) THEN 

       !Get saturation specific humidity and gradient wrt to T
       pke  = (ple(i,:)/1000.)**(cons_rgas/cons_cp)
       pf   = 0.5*(ple(i,1:k0) +  ple(i,2:k0+1  ) )
       pk   = (pf/1000.)**(cons_rgas/cons_cp)
       tempf = tho(i,:)*pk

       zle = 0.0
       zlo = 0.0
       zle(k0+1) = 0.
       do l=k0,1,-1
          zle(l) = tho(i,l)   * (1.+cons_vireps*qho(i,l))
          zlo(l) = zle(l+1) + (cons_cp/cons_grav)*( pke(l+1)-pk(l  ) ) * zle(l)
          zle(l) = zlo(l)   + (cons_cp/cons_grav)*( pk(l)  -pke(l)   ) * zle(l)
       end do

       tpert  = cbl_tpert * ( ts(i) - ( tempf(k0)+ cons_grav*zlo(k0)/cons_cp )  ) 
       qpert  = cbl_qpert !* ( QSSFC - Q(:,:,K0) ) [CBL_QPERT = 0.0]
       tpert  = max( tpert , 0.0 )
       qpert  = max( qpert , 0.0 )
      
       if (frland(i) < 0.1) then
          tpert = min( tpert , cbl_tpert_mxocn ) ! ocean
       else
          tpert = min( tpert , cbl_tpert_mxlnd ) ! land
       endif

       call dqsat_ras(dqs,qss,tempf,pf,k0,estblx,cons_h2omw,cons_airmw)

       do kk=icmin,k+1
          prj(kk) = pke(kk)
       enddo

       poi=0.        ! These initialized here in order not to confuse Valgrind debugger
       qoi=0.        ! Do not believe it actually makes any difference.

       prs(icmin:k0+1) = ple(i,icmin:k0+1)
       poi(icmin:k)   = tho(i,icmin:k)
       qoi(icmin:k)   = qho(i,icmin:k)

       qst(icmin:k) = qss(icmin:k)
       dqq(icmin:k) = dqs(icmin:k)

       !Mass fraction of each layer below cloud base
       massf(:) = wgt0(i,:)

       !RESET PRESSURE at bottom edge of CBL 
       prcbl = prs(k)
       do l= k,k0
          prcbl = prcbl + massf(l)*( prs(l+1)-prs(l) )
       end do
       prs(k+1) = prcbl
       prj(k+1) = (prs(k+1)/1000.)**(cons_rgas/cons_cp)

       DO l=k,icmin,-1
          pol(l)  = 0.5*(prs(l)+prs(l+1))
          prh(l)  = (prs(l+1)*prj(l+1)-prs(l)*prj(l)) / (onepkap*(prs(l+1)-prs(l)))
          pki(l)  = 1.0 / prh(l)
          dpt(l)  = prh(l  ) - prj(l)
          dpb(l)  = prj(l+1) - prh(l)
          pri(l)  = .01 / (prs(l+1)-prs(l))
       ENDDO

       !RECALCULATE PROFILE QUAN. IN LOWEST STRAPPED LAYER
       if ( k<=k0) then
          poi(k) = 0.
          qoi(k) = 0.

          !SPECIFY WEIGHTS GIVEN TO EACH LAYER WITHIN SUBCLOUD "SUPERLAYER"
          wght = 0.
          DO l=k,k0
             wght(l)   = massf(l) * ( ple(i,l+1) - ple(i,l) ) / ( prs(k+1)   - prs(k)  )
          END DO      

          DO l=k,k0
             poi(k) = poi(k) + wght(l)*tho(i,l)
             qoi(k) = qoi(k) + wght(l)*qho(i,l)
          ENDDO

          call dqsats_ras(dqq(k), qst(k), poi(k)*prh(k),pol(k),estblx,cons_h2omw,cons_airmw)

       endif

       rndu = max( seedras(i)/1000000., 1e-6 )
       mxdiam = maxdallowed*( rndu**(-1./2.) )

       DO l=k,icmin,-1
          bet(l)  = dqq(l)*pki(l)  !*
          gam(l)  = pki(l)/(1.0+lbcp*dqq(l)) !*
          IF (l<k) THEN
             ght(l+1) = gam(l)*dpb(l) + gam(l+1)*dpt(l+1)
             gm1(l+1) = 0.5*lbcp*(dqq(l  )/(alhl*(1.0+lbcp*dqq(l  ))) + dqq(l+1)/(alhl*(1.0+lbcp*dqq(l+1))) )
          ENDIF
       ENDDO

       rns  = 0.
       cll  = 0.
       rmf  = 0.
       rmfd = 0.
       rmfc = 0.
       rmfp = 0.
       cll0 = 0.
       dll0 = 0.
       cllx = 0.
       dllx = 0.
       clli = 0.

       poi_sv = poi
       qoi_sv = qoi

       cvw     = 0.0
       updfrc  = 0.0
       updfrp  = 0.0

       hol=0.            ! HOL initialized here in order not to confuse Valgrind debugger
       zet(k+1) = 0
       sht(k+1) = cp*poi(k)*prj(k+1)
       DO l=k,icmin,-1
          qol(l)  = min(qst(l)*rhmax,qoi(l))
          qol(l)  = max( 0.000, qol(l) )     ! GUARDIAN/NEG.PREC.
          ssl(l)  = cp*prj(l+1)*poi(l) + grav*zet(l+1)
          hol(l)  = ssl(l) + qol(l)*alhl
          hst(l)  = ssl(l) + qst(l)*alhl
          tem     = poi(l)*(prj(l+1)-prj(l))*cpbg
          zet(l)  = zet(l+1) + tem
          zol(l)  = zet(l+1) + (prj(l+1)-prh(l))*poi(l)*cpbg
       ENDDO

       DO ic =  k,icmin+1,-1
              
          alm   = 0.
          trg   = min(1.,(qoi(k)/qst(k)-rhmn)/(rhmx-rhmn))

          f4  = min(   1.0,  max( 0.0 , (autorampb-sige(ic))/0.2 )  )  
     
          IF (trg <= 1.0e-5) THEN
             cycle !================>>
          ENDIF

          !RECOMPUTE SOUNDING UP TO DETRAINMENT LEVEL
          poi_c = poi
          qoi_c = qoi
          poi_c(k) =  poi_c(k) + tpert
          qoi_c(k) =  qoi_c(k) + qpert

          zet(k+1) = 0.
          sht(k+1) = cp*poi_c(k)*prj(k+1)
          DO l=k,ic,-1
             qol(l)  = min(qst(l)*rhmax,qoi_c(l))
             qol(l)  = max( 0.000, qol(l) )     ! GUARDIAN/NEG.PREC.
             ssl(l)  = cp*prj(l+1)*poi_c(l) + grav*zet(l+1)
             hol(l)  = ssl(l) + qol(l)*alhl
             hst(l)  = ssl(l) + qst(l)*alhl
             tem     = poi_c(l)*(prj(l+1)-prj(l))*cpbg
             zet(l)  = zet(l+1) + tem
             zol(l)  = zet(l+1) + (prj(l+1)-prh(l))*poi_c(l)*cpbg
          ENDDO

          DO l=ic+1,k
             tem  = (prj(l)-prh(l-1))/(prh(l)-prh(l-1))
             sht(l)  = ssl(l-1) + tem*(ssl(l)-ssl(l-1)) 
             qht(l)  = .5*(qol(l)+qol(l-1))
          ENDDO

          !CALCULATE LAMBDA, ETA, AND WORKFUNCTION
          lambda_min = .2/mxdiam
          lambda_max = .2/  200. 

          IF (hol(k) <= hst(ic)) THEN
             cycle !================>>
          ENDIF

          !LAMBDA CALCULATION: MS-A18

          tem  =       (hst(ic)-hol(ic))*(zol(ic)-zet(ic+1)) 
          DO l=ic+1,k-1
             tem = tem + (hst(ic)-hol(l ))*(zet(l )-zet(l +1))
          ENDDO

          IF (tem <= 0.0) THEN
             cycle !================>>
          ENDIF

          alm = (hol(k)-hst(ic)) / tem

          IF (alm > lambda_max) THEN
             cycle !================>>
          ENDIF

          toki=1.0

          IF (alm < lambda_min) THEN
             toki = ( alm/lambda_min )**2
          ENDIF

          !ETA CALCULATION: MS-A2
          DO l=ic+1,k
             eta(l) = 1.0 + alm * (zet(l )-zet(k))
          ENDDO
          eta(ic) = 1.0 + alm * (zol(ic)-zet(k))

          !WORKFUNCTION CALCULATION:  MS-A22
          wfn     = 0.0
          hcc(k)  = hol(k)
          DO l=k-1,ic+1,-1
             hcc(l) = hcc(l+1) + (eta(l) - eta(l+1))*hol(l)
             tem    = hcc(l+1)*dpb(l) + hcc(l)*dpt(l)
             eht(l) = eta(l+1)*dpb(l) + eta(l)*dpt(l)
             wfn    = wfn + (tem - eht(l)*hst(l))*gam(l)
          ENDDO
          hcc(ic) = hst(ic)*eta(ic)
          wfn     = wfn + (hcc(ic+1)-hst(ic)*eta(ic+1))*gam(ic)*dpb(ic)

          !VERTICAL VELOCITY/KE CALCULATION (ADDED 12/2001 JTB)
          bk2(k)   = 0.0
          bke(k)   = 0.0
          hcld(k)  = hol(k)
          DO l=k-1,ic,-1
             hcld(l) = ( eta(l+1)*hcld(l+1) + (eta(l) - eta(l+1))*hol(l) ) / eta(l)
             tem     = (hcld(l)-hst(l) )*(zet(l)-zet(l+1))/ (1.0+lbcp*dqq(l))
             bke(l)  = bke(l+1) + grav * tem / ( cp*prj(l+1)*poi(l) )
             bk2(l)  = bk2(l+1) + grav * max(tem,0.0) / ( cp*prj(l+1)*poi(l) )
             cvw(l) = sqrt(  2.0* max( bk2(l) , 0.0 )  )    
          ENDDO

          !ALPHA CALCULATION 
          IF ( zet(ic) <  2000. ) THEN
             rasal(ic) = rasal1
          ENDIF
          IF ( zet(ic)  >= 2000. ) THEN 
             rasal(ic) = rasal1 + (rasal2-rasal1)*(zet(ic) - 2000.)/8000.
          ENDIF
          rasal(ic) = min( rasal(ic) , 1.0e5 )
    
          rasal(ic) = dt / rasal(ic)

          cvw(ic:k) = max(  cvw(ic:k) , 1.00 )

          CALL acritn(pol(ic), prs(k), acr, acritfac)

          IF (wfn <= acr) THEN
             cycle !================>>
          ENDIF

          wlq     = qol(k)
          rnn(k)  = 0.
          cll0(k) = 0.

          DO l=k-1,ic,-1
             tem   = eta(l) - eta(l+1)
             wlq   = wlq + tem * qol(l)

             IF (l>ic) THEN
                tx2   = 0.5*(qst(l)+qst(l-1))*eta(l)
                tx3   = 0.5*(hst(l)+hst(l-1))*eta(l)
                qcc   = tx2 + gm1(l)*(hcc(l)-tx3)
                cll0(l)  = (wlq-qcc)
             ELSE
                cll0(l)   = (wlq-qst(ic)*eta(ic))
             ENDIF

             cll0(l)    = max( cll0(l) , 0.00 )

             cli  = cll0(l) / eta(l)
             te_a = poi(l)*prh(l)

             CALL sundq3_ice( te_a,  sdqv2, sdqv3, sdqvt1, f2 , f3 )

             c00_x  =  co_auto(i) * f2 * f3  * f4 
             cli_crit_x =  cli_crit / ( f2 * f3 )
       
             rate = c00_x * ( 1.0 - exp( -(cli)**2 / cli_crit_x**2 ) )

             cvw_x     = max( cvw(l) , 1.00 )  

             dt_lyr  = ( zet(l)-zet(l+1) )/cvw_x 

             closs   = cll0(l) * rate * dt_lyr
             closs   = min( closs , cll0(l) )

             cll0(l) = cll0(l) - closs
             dll0(l) = closs

             IF (closs>0.) then
                wlq = wlq - closs
                rnn(l) = closs 
             else
                rnn(l) = 0.
             ENDIF

          ENDDO

          wlq = wlq - qst(ic)*eta(ic)

          !CALCULATE GAMMAS AND KERNEL
          gms(k) =          (sht(k)-ssl(k))*pri(k)
          gmh(k) = gms(k) + (qht(k)-qol(k))*pri(k)*alhl
          akm    = gmh(k)*gam(k-1)*dpb(k-1)

          tx2     = gmh(k)
          DO l=k-1,ic+1,-1
             gms(l) = ( eta(l  )*(sht(l)-ssl(l  )) + eta(l+1)*(ssl(l)-sht(l+1)) )     *pri(l)
             gmh(l) = gms(l) + ( eta(l  )*(qht(l)-qol(l  )) + eta(l+1)*(qol(l)-qht(l+1)) )*alhl*pri(l)
             tx2 = tx2 + (eta(l) - eta(l+1)) * gmh(l)
             akm = akm - gms(l)*eht(l)*pki(l) + tx2*ght(l)
          ENDDO

          gms(ic) = eta(ic+1)*(ssl(ic)-sht(ic+1))*pri(ic)
          akm     = akm - gms(ic)*eta(ic+1)*dpb(ic)*pki(ic)

          gmh(ic) =   gms(ic) + ( eta(ic+1)*(qol(ic)-qht(ic+1))*alhl + eta(ic)*(hst(ic)-hol(ic)))*pri(ic)

          !CLOUD BASE MASS FLUX
          IF (akm >= 0.0 .OR. wlq < 0.0)  THEN
             cycle !================>>
          ENDIF

          wfn = - (wfn-acr)/akm
          wfn = min( ( rasal(ic)*trg*toki )*wfn  ,   (prs(k+1)-prs(k) )*(100.*pblfrac))
     
          !CUMULATIVE PRECIP AND CLOUD-BASE MASS FLUX FOR OUTPUT
          tem      = wfn*gravi
          cll(ic) = cll(ic) + wlq*tem
          rmf(ic) = rmf(ic) +     tem
          rmfd(ic) = rmfd(ic) +     tem * eta(ic)

          DO l=ic+1,k
             rmfp(l) = tem * eta(l)
             rmfc(l) = rmfc(l)  +  rmfp(l)
           
             dllx(l) = dllx(l)  +  tem*dll0(l)        

             IF ( cvw(l) > 0.0 ) THEN
                updfrp(l) = rmfp(l) * (ddt/daylen) * 1000. / ( cvw(l) * prs(l) )
             ELSE 
                updfrp(l) = 0.0       
             ENDIF

             clli(l) = cll0(l)/eta(l)

             updfrc(l) =  updfrc(l) +  updfrp(l)      
 
          ENDDO

          !THETA AND Q CHANGE DUE TO CLOUD TYPE IC
          DO l=ic,k
             rns(l) = rns(l) + rnn(l)*tem
             gmh(l) = gmh(l) * wfn
             gms(l) = gms(l) * wfn
             qoi(l) = qoi(l) + (gmh(l) - gms(l)) * alhi
             poi(l) = poi(l) + gms(l)*pki(l)*cpi
             qst(l) = qst(l) + gms(l)*bet(l)*cpi
          ENDDO

       ENDDO !CLOUD LOOP

       IF ( sum( rmf(icmin:k) ) > 0.0 ) THEN

          DO l=icmin,k
             tem    = pri(l)*grav
             cnv_prc3(i,l) = rns(l)*tem     
          ENDDO

          tho(i,icmin:k-1) = poi(icmin:k-1)
          qho(i,icmin:k-1) = qoi(icmin:k-1)
          cnv_updfrc(i,icmin:k-1)   =  updfrc(icmin:k-1)

          !De-strap tendencies from RAS
          wght   = wgt1(i,:)

          !Scale properly by layer masses
          wght0 = 0.
          DO l=k,k0 
             wght0 = wght0 + wght(l)* ( ple(i,l+1) - ple(i,l) )
          END DO
         
          wght0 = ( prs(k+1)   - prs(k)  )/wght0
          wght  = wght0 * wght

          DO l=k,k0 
             tho(i,l) =  tho(i,l) + wght(l)*(poi(k) - poi_sv(k))
             qho(i,l) =  qho(i,l) + wght(l)*(qoi(k) - qoi_sv(k))
          END DO

          flxd(i,icmin:k) = rmfd(icmin:k) * ddt/daylen
          clw(i,icmin:k) = cll(icmin:k) * ddt/daylen

          flxd(i,1:icmin-1) = 0.
          clw(i,1:icmin-1) = 0.

          IF ( k < k0 ) THEN 
             flxd(i,k:k0) = 0.
             clw(i,k:k0) = 0.
          END IF

       ELSE

          flxd(i,:) = 0.
          clw(i,:) = 0.

       ENDIF 

    ELSE 
     
       flxd(i,:) = 0.
       clw(i,:) = 0.

    ENDIF

 ENDDO

END SUBROUTINE rase0

END MODULE convection