SUBROUTINE cucnvc
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  ! subprogram:    cucnvc      convective precipitation parameterization
  !   prgrmmr: janjic          org: w/np2      date: 93-11-02       
  !     
  ! abstract:
  !     cucnvc calculates the sub-grid scale convection including 
  !     deep and shallow convective clouds following the scheme described
  !     by betts (1986) but with significant modifications.  in addition,
  !     the latent heat release and moisture change due to precipitating
  !     and non-precipitating clouds are computed.
  !     
  !     
  ! program history log:
  !   87-09-??  janjic     - originator
  !   90-11-21  janjic     - two sets of dsp profiles (fast and slow)
  !                          replace the original one set
  !   95-03-25  black      - conversion from 1-d to 2-d in horizontal
  !   96-03-28  black      - added external edge
  !     
  ! usage: call cucnvc from main program ebu
  !
  !   input argument list:
  !     none     
  !  
  !   output argument list: 
  !     none
  !     
  !   output files:
  !     none
  !     
  !   subprograms called:
  !  
  !     unique:
  !        ttblex
  !  
  !     library:
  !        none
  !  
  !   common blocks: ctlblk
  !                  loops
  !                  masks
  !                  phys
  !                  vrbls
  !                  cnvcld
  !                  pvrbls
  !                  acmclh
  !                  indx
  !     ******************************************************************
  !     *                                                                *
  !     *  references:                                                   *
  !     *                                                                *
  !     *  betts, a.k., 1986:  a new convective adjustment scheme.       *
  !     *    part i: observational and theoretical basis.  quart. j. r.  *
  !     *    met. soc., 112, 677-691.                                    *
  !     *                                                                *
  !     *  betts, a.k., and m.j. miller, 1986:  a new convective         *
  !     *    adjustment scheme.  part ii: single column tests using      *
  !     *    gate wave, bomex, atex and arctic air mass data sets.       *
  !     *    quart. j. r. met. soc., 112, 693-709.                       *
  !     *                                                                *
  !     *                                                                *
  !     ******************************************************************
  ! *** warning: this subroutine will not work if lm.lt.12;
  !              must be called in the same step with profq2 because profq
  !              defines ape;
  IMPLICIT NONE
  INCLUDE "cuparm.f90"
  INCLUDE "parmeta.f90"
  INCLUDE "parm.tbl.f90"
  INTEGER, PARAMETER :: la=13,ksmud=0,nrow= 0,imjm=im*jm-jm/2
  INTEGER, PARAMETER :: jam=6+2*(jm-10),lp1=lm+1,lm1=lm-1,lscrch=4*lm+1+la+1 
  INTEGER, PARAMETER :: l1=la+1,l2=la+lm+1,l3=la+2*lm+1,l4=la+3*lm+1
  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "phys.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "cnvcld.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "acmclh.f90"
  INCLUDE "indx.f90"
  !
  INTEGER, DIMENSION(lm) :: ntopd,nbotd,ntops,nbots,ndpthd,ndpths
  REAL   , DIMENSION(lm) :: trefk,qrefk,pk,apek,tk,thsk,psk,apesk,qk,therk
  REAL   , DIMENSION(lm) :: thvref,dift,difq,qsatk,fpk
  INTEGER, DIMENSION(im,jm) :: ltop,lbot,iptb,ithtb
  INTEGER, DIMENSION(imjm ) :: ishal,jshal,ideep,jdeep,ilres,jlres,ihres,jhres
  LOGICAL, DIMENSION(im,jm) :: correctp
  REAL   , DIMENSION(im,jm) :: ptop,pbot,pdsl,pp,qq,psp,thbt,thesp,p
  REAL   , DIMENSION(im,jm) :: wf1,wf2,wf3,wf4,precol,vpsfck
  REAL   , DIMENSION(im,jm,lm) :: ape,tref,tmod,qmod
  REAL   , DIMENSION(0:im+1,0:jm+1) :: dspb,dsp0,dspt,tl,ql,tne,tse,qne,qse
  REAL    :: aetal,apekl,apekxx,apekxy,apesp,avrgt,avrgtl,bq,bqk,bqs00k,bqs10k
  REAL    :: cthrs,den,dentpy,dhdt,difql,diftl,dpkl,dpmix,dqref,drhdp,drheat
  REAL    :: dsp,dsp0k,dspbk,dsptk,dst,dstq,dtcnvc,dtdeta,dthem,efi,fefi,fptk
  REAL    :: hcorr,otsum,pbtk,pk0,pkb,pkl,pkt,pkxxxx,pkxxxy,potsum,pp1,ppk
  REAL    :: preck,presk,psum,pthrs,ptpk,pz0,qkl,qnew,qotsum,qq1,qqk
  REAL    :: qrfkl,qrftp,qsum,rdp0t,rdpsum,rdtcnvc,rhh,rhl,rhmax,rotsum
  REAL    :: rsmk,rtbar,smix,smk,sq,sqk,sqs00k,sqs10k,stabdl,sumde,sumdp
  REAL    :: sumdt,tauk,tcorr,therkx,therky,thskl,thtpk,thvmkl,tkl,tpsp,tq,tqk
  REAL    :: trefkx,trfkl,tskl,tth,tthbt,tthes,tthk
  INTEGER :: i,iq,iqtb,it,iter,ittb,ivi,j,kb,khdeep,khshal,knumh,knuml,kount
  INTEGER :: krow,ks,l,l0,l0m1,lb,lbm1,lbtk,lcor,ll,lmhij,lqm,lshu,ltp1,ltp2
  INTEGER :: ltpk,ltsh,n,ndeep,ndepth,ndstn,ndstp,negds,nneg,nshal,nsmud

  dspb = 0.0
  dsp0 = 0.0
  dspt = 0.0
  tl = 0.0
  ql = 0.0
  tne = 0.0
  tse = 0.0
  qne = 0.0
  qse = 0.0

  avcnvc=avcnvc+1.
  acutim=acutim+1.
  dtcnvc=ncnvc*dt
  rdtcnvc=1./dtcnvc
  tauk=dtcnvc/trel
  cthrs=(0.006350/86400.)*dtcnvc/cprlg

  IF(ntsd.EQ.1)THEN
     cldefi=avgefi*sm+stefi*(h1-sm)
  ENDIF

  !*** diagnostics

  ntopd=0
  nbotd =0
  ntops =0
  nbots =0
  ndpths=0
  ndpthd=0

  !*** preparations

  DO j=1,jm
     DO i=1,im
        lbot (i,j)=lmh(i,j)
        thesp(i,j)=d00
        pdsl (i,j)=res(i,j)*pd(i,j)
        precol(i,j)=d00
        tref(i,j,1)=t(i,j,1)
     END DO
  END DO

  !*** padding specific humidity if too small
  !    restore ape to scratch array

  DO l=1,lm
     DO j=1,jm
        DO i=1,im
           ape(i,j,l)=(h10e5/(pdsl(i,j)*aeta(l)+pt))**capa
        END DO
     END DO
  END DO
  DO l=1,lm
     DO j=1,jm
        DO i=1,im
           IF(q(i,j,l).LT.epsq)q(i,j,l)=htm(i,j,l)*epsq
        END DO
     END DO
  END DO

  !*** search for maximum buoyancy level

  DO kb=1,lm
     DO j=3,jm-2
        DO i=1,im
           IF (kb <= lmh(i,j)) THEN
              IF (aeta(kb)*pdsl(i,j)+pt >=  &
                   aeta(lmh(i,j))*pdsl(i,j)+pt-pbm ) THEN

                 !*** trial maximum buoyancy level variables

                 tthbt=t(i,j,kb)*ape(i,j,kb)
                 tth=(tthbt-thl)*rdth
                 ittb=INT(tth)+1

                 !*** keeping indices within the table

                 IF(ittb.LT.1)THEN
                    ittb=1
                    qq1=d00
                 ELSE IF(ittb.GE.jtb)THEN
                    ittb=jtb-1
                    qq1=d00
                 ELSE
                    qq1=tth-AINT(tth)
                 ENDIF

                 !*** scaling spec. humidity & table index

                 bq  = (qs0(ittb+1)-qs0(ittb))*qq1+qs0(ittb)
                 sq  = (sqs(ittb+1)-sqs(ittb))*qq1+sqs(ittb)
                 tq  = (q(i,j,kb)-bq)/sq*rdq
                 iqtb=INT(tq)+1

                 !*** keeping indices within the table

                 IF(iqtb.LT.1)THEN
                    iqtb=1
                    pp1=d00
                 ELSEIF(iqtb.GE.itb)THEN
                    iqtb=itb-1
                    pp1=d00
                 ELSE
                    pp1 = tq - AINT(tq)
                 ENDIF

                 !*** saturation point variables at the bottom

                 tpsp=ptbl(iqtb,ittb)+&
                      (ptbl(iqtb+1,ittb)-ptbl(iqtb,ittb))*pp1+&
                      (ptbl(iqtb,ittb+1)-ptbl(iqtb,ittb))*qq1+&
                      (ptbl(iqtb,ittb)-ptbl(iqtb+1,ittb)-&
                      ptbl(iqtb,ittb+1)+ptbl(iqtb+1,ittb+1))*pp1*qq1
                 apesp=(h10e5/tpsp)**capa
                 tthes=tthbt*EXP(elocp*q(i,j,kb)*apesp/tthbt)

                 !*** check for maximum buoyancy

                 IF(tthes.GT.thesp(i,j))THEN
                    psp  (i,j)=tpsp
                    thbt (i,j)=tthbt
                    thesp(i,j)=tthes
                 END IF
              ENDIF
           END IF
        END DO
     END DO
  END DO

  !*** choose cloud base as model level just below psp

  DO l=1,lm1
     aetal=aeta(l)
     DO j=3,jm-2
        DO i=1,im
           p(i,j)=pdsl(i,j)*aetal+pt
           IF(p(i,j).LT.psp(i,j).AND.p(i,j).GE.pqm)lbot(i,j)=l+1
        END DO
     END DO
  END DO

  !*** warning: lbot must not be gt lmh(i,j)-1 in shallow convection
  !*** make sure cloud base is at least pone above the surface

  DO j=3,jm-2
     DO i=1,im
        lmhij=lmh(i,j)
        pbot(i,j)=aeta(lbot(i,j))*pdsl(i,j)+pt
        vpsfck(i,j)=aeta(lmhij)*pdsl(i,j)+pt
        correctp(i,j) = pbot(i,j).GE.vpsfck(i,j)-pone.OR.lbot(i,j).GE.lmhij
     END DO
  END DO

  DO j=3,jm-2
     DO i=1,im
        IF (correctp(i,j)) THEN
           DO l=1,lmh(i,j)-1
              p(i,j)=aeta(l)*pdsl(i,j)+pt
              IF(p(i,j).LT.vpsfck(i,j)-pone)lbot(i,j)=l
           END DO
           pbot(i,j)=aeta(lbot(i,j))*pdsl(i,j)+pt
        ENDIF
     END DO
  END DO

  !*** cloud top computation

  DO j=1,jm
     DO i=1,im
        ltop(i,j)=lbot(i,j)
        ptop(i,j)=pbot(i,j)
     END DO
  END DO

  !*** scaling pressure & tt table index

  DO l=lm,1,-1
     knuml=0
     knumh=0
     DO j=3,jm-2
        DO i=1,im
           presk=pdsl(i,j)*aeta(l)+pt
           IF(presk.LT.plq)THEN
              knuml=knuml+1
              ilres(knuml)=i
              jlres(knuml)=j
           ELSE
              knumh=knumh+1
              ihres(knumh)=i
              jhres(knumh)=j
           ENDIF
        END DO
     END DO

     !***  compute parcel temperature along moist adiabat for pressure<pl

     IF(knuml.GT.0)THEN
        CALL ttblex(tref(1,1,l),ttbl,itb,jtb,knuml,ilres,jlres,pdsl, &
             aeta(l),htm(0,0,l),pt,pl,qq,pp,rdp,the0,sthe, &
             rdthe,thesp,iptb,ithtb)
     ENDIF

     !***  compute parcel temperature along moist adiabat for pressure>pl

     IF(knumh.GT.0)THEN
        CALL ttblex(tref(1,1,l),ttblq,itbq,jtbq,knumh,ihres,jhres,pdsl, &
             aeta(l),htm(0,0,l),pt,plq,qq,pp,rdpq,the0q,stheq, &
             rdtheq,thesp,iptb,ithtb)
     ENDIF
  END DO

  !*** buoyancy check

  DO l=lm,1,-1
     DO j=3,jm-2
        DO i=1,im
           IF(tref(i,j,l).GT.t(i,j,l)-dttop)ltop(i,j)=l
        ENDDO
     ENDDO
  END DO

  !*** cloud top pressure

  DO j=3,jm-2
     DO i=1,im
        ptop(i,j)=aeta(ltop(i,j))*pdsl(i,j)+pt
     END DO
  END DO

  !*** define and smooth dsps and cldefi
  !*** unified or separate land/sea conv

  IF(unis)THEN
     DO j=1,jm
        DO i=1,im
           efi=cldefi(i,j)
           dspb(i,j)=(efi-efimn)*slopbs+dspbss
           dsp0(i,j)=(efi-efimn)*slop0s+dsp0ss
           dspt(i,j)=(efi-efimn)*slopts+dsptss
        END DO
     END DO
  ELSEIF(.NOT.unil)THEN
     DO j=1,jm
        DO i=1,im
           efi=cldefi(i,j)
           dspb(i,j)=((efi-efimn)*slopbs+dspbss)*sm(i,j) &
                +((efi-efimn)*slopbl+dspbsl)*(h1-sm(i,j))
           dsp0(i,j)=((efi-efimn)*slop0s+dsp0ss)*sm(i,j) &
                +((efi-efimn)*slop0l+dsp0sl)*(h1-sm(i,j))
           dspt(i,j)=((efi-efimn)*slopts+dsptss)*sm(i,j) &
                +((efi-efimn)*sloptl+dsptsl)*(h1-sm(i,j))
        END DO
     END DO
  ELSE
     DO j=1,jm
        DO i=1,im
           efi=cldefi(i,j)
           dspb(i,j)=((efi-efimn)*slopbl+dspbsl)
           dsp0(i,j)=((efi-efimn)*slop0l+dsp0sl)
           dspt(i,j)=((efi-efimn)*sloptl+dsptsl) 
        END DO
     END DO
  ENDIF

  !*** extending sea structures inland along coastline

  IF(nrow.GT.0.AND..NOT.unis.AND..NOT.unil)THEN
     DO j=1,jm
        DO i=1,im
           wf1(i,j)=0.
           wf2(i,j)=0.
           wf3(i,j)=0.
           wf4(i,j)=0.
        END DO
     END DO

     krow=nrow

     DO kount=1,krow
        DO j=3,jm-2
           DO i=1,im
              wf1(i,j)=(dspb(i+ihw(j),j-1)+dspb(i+ihe(j),j-1) &
                   +dspb(i+ihw(j),j+1)+dspb(i+ihe(j),j+1)+4.*dspb(i,j)) &
                   *hbm2(i,j)*0.125
              wf2(i,j)=(dsp0(i+ihw(j),j-1)+dsp0(i+ihe(j),j-1) &
                   +dsp0(i+ihw(j),j+1)+dsp0(i+ihe(j),j+1)+4.*dsp0(i,j)) &
                   *hbm2(i,j)*0.125
              wf3(i,j)=(dspt(i+ihw(j),j-1)+dspt(i+ihe(j),j-1) &
                   +dspt(i+ihw(j),j+1)+dspt(i+ihe(j),j+1)+4.*dspt(i,j)) &
                   *hbm2(i,j)*0.125
              wf4(i,j)=(cldefi(i+ihw(j),j-1)+cldefi(i+ihe(j),j-1) &
                   +cldefi(i+ihw(j),j+1)+cldefi(i+ihe(j),j+1)+4.*cldefi(i,j)) &
                   *hbm2(i,j)*0.125
           END DO
        END DO

        DO j=1,jm
           DO i=1,im
              smk =sm(i,j)
              rsmk=1.-smk
              dspb  (i,j)=dspb  (i,j)*smk+wf1(i,j)*rsmk
              dsp0  (i,j)=dsp0  (i,j)*smk+wf2(i,j)*rsmk
              dspt  (i,j)=dspt  (i,j)*smk+wf3(i,j)*rsmk 
              cldefi(i,j)=cldefi(i,j)*smk+wf4(i,j)*rsmk
           END DO
        END DO
     END DO
  ENDIF

  !*** initialize changes of t and q due to convection

  DO l=1,lm
     DO j=1,jm
        DO i=1,im
           tmod(i,j,l)=d00
           qmod(i,j,l)=d00
        END DO
     END DO
  END DO

  !*** clean up and gather deep convection points

  DO j=3,jm-2
     DO i=1,im
        IF(ltop(i,j).GT.lbot(i,j))THEN
           ltop(i,j)=lbot(i,j)
           ptop(i,j)=pbot(i,j)
        ENDIF
        IF(hbm2(i,j).LT.0.90)THEN
           ltop(i,j)=lbot(i,j)
           ptop(i,j)=pbot(i,j)
        ENDIF
        IF(ptop(i,j).GT.pbot(i,j)-pno.OR.ltop(i,j).GT.lbot(i,j)-2) &
             cldefi(i,j)=avgefi*sm(i,j)+stefi*(h1-sm(i,j))
     END DO
  END DO

  khdeep=0
  DO j=3,jm-2
     DO i=1,im
        IF(ptop(i,j).LT.pbot(i,j)-psh)THEN
           khdeep=khdeep+1
           ideep(khdeep)=i
           jdeep(khdeep)=j
        ENDIF
     END DO
  END DO

  !*** horizontal loop for deep convection

  DO n=1,khdeep
     i=ideep(n)
     j=jdeep(n)
     ltpk=ltop(i,j)
     lbtk=lbot(i,j)

     !***  deep convection

     lb   =lbtk
     efi  =cldefi(i,j)
     dspbk=dspb  (i,j)
     dsp0k=dsp0  (i,j)
     dsptk=dspt  (i,j)

     !*** initialize variables in the convective column

     !***  one should note that the values assigned to the array trefk
     !***  in the next loop are really only relevant in anchoring the
     !***  reference temperature profile at level lb.  when building the
     !***  reference profile from cloud base, then assigning the
     !***  ambient temperature to trefk is acceptable.  however, when
     !***  building the reference profile from some other level (such as
     !***  one level above the ground), then trefk should be filled with
     !***  the temperatures in tref(i,j,l) which are the temperatures of
     !***  the moist adiabat through cloud base.  by the time the line 
     !***  numbered 450 has been reached, trefk actually does hold the
     !***  reference temperature profile.

     DO l=1,lm
        dift  (l)=d00
        difq  (l)=d00
        tkl      =t(i,j,l)
        tk    (l)=tkl
        trefk (l)=tkl
        qkl      =q(i,j,l)
        qk    (l)=qkl
        qrefk (l)=qkl
        pkl      =aeta(l)*pdsl(i,j)+pt
        pk    (l)=pkl
        psk   (l)=pkl
        apekl    =ape(i,j,l)
        apek  (l)=apekl
        therk (l)=tref(i,j,l)*apekl
     END DO

     !*** deep convection reference temperature profile

     ltp1=ltpk+1
     lbm1=lb-1
     pkb=pk(lb)
     pkt=pk(ltpk)

     !*** temperature reference profile below freezing level

     l0=lb
     pk0=pk(lb)
     trefkx=trefk(lb)
     therkx=therk(lb)
     apekxx=apek(lb)
     therky=therk(lbm1)
     apekxy=apek(lbm1)
     DO l=ltpk,lbm1
        ivi=ltpk+lbm1-l

        !*** temperature reference profile above freezing level

        IF(t(i,j,ivi+1).LT.t0)THEN
           l0m1=l0-1
           rdp0t=h1/(pk0-pkt)
           dthem=therk(l0)-trefk(l0)*apek(l0)
           DO ll=ltpk,l0m1
              trefk(ll)=(therk(ll)-(pk(ll)-pkt)*dthem*rdp0t)/apek(ll)
           END DO
           EXIT
        END IF

        stabdl=stabd
        trefkx=((therky-therkx)*stabdl+trefkx*apekxx)/apekxy
        trefk(ivi)=trefkx
        apekxx=apekxy
        therkx=therky
        apekxy=apek(ivi-1)
        therky=therk(ivi-1)
        l0=ivi
        pk0=pk(l0)
     END DO

     !*** freezing level at or above the cloud top

     l0m1=l0-1

     !*** deep convection reference humidity profile

     DO l=ltpk,lb

        !*** saturation pressure difference

        IF(pkb-pk0.GE.pfrz)THEN

           IF(l.LT.l0)THEN
              dsp=((pk0-pk(l))*dsptk+(pk(l)-pkt)*dsp0k)/(pk0-pkt)
           ELSE
              dsp=((pkb-pk(l))*dsp0k+(pk(l)-pk0)*dspbk)/(pkb-pk0)
           ENDIF
        ELSE
           dsp=dspc
        ENDIF

        !*** humidity profile

        IF(pk(l).GT.pqm)THEN
           psk(l)=pk(l)+dsp
           apesk(l)=(h10e5/psk(l))**capa
           thsk(l)=trefk(l)*apek(l)
           qrefk(l)=pq0/psk(l)*EXP(a2*(thsk(l)-a3*apesk(l)) &
                /(thsk(l)-a4*apesk(l)))
        ELSE
           qrefk(l)=q(i,j,l)
        ENDIF
     END DO

     !*** enthalpy conservation integral

     DO iter=1,2
        sumde=d00
        sumdp=d00
        DO l=ltpk,lb
           sumde=((tk(l)-trefk(l))*cp+(qk(l)-qrefk(l))*elwv)*deta(l) &
                +sumde
           sumdp=sumdp+deta(l)
        END DO
        hcorr=sumde/(sumdp-deta(ltpk))
        lcor=ltpk+1

        !*** find lqm
        DO l=1,lb
           IF(pk(l).LE.pqm)lqm=l
        END DO

        !*** above lqm correct temperature only

        IF(lcor.LE.lqm)THEN
           DO l=lcor,lqm
              trefk(l)=trefk(l)+hcorr*rcp
           END DO
           lcor=lqm+1
        ENDIF

        !*** below lqm correct both temperature and moisture

        DO l=lcor,lb
           tskl=trefk(l)*apek(l)/apesk(l)
           dhdt=qrefk(l)*a23m4l/(tskl-a4)**2+cp
           trefk(l)=hcorr/dhdt+trefk(l)
           thskl=trefk(l)*apek(l)
           qrefk(l)=pq0/psk(l)*EXP(a2*(thskl-a3*apesk(l)) &
                /(thskl-a4*apesk(l)))
        END DO
     END DO

     !*** heating, moistening, precipitation

     dentpy=d00
     avrgt =d00
     preck =d00
     DO l=ltpk,lb
        tkl    =tk(l)
        diftl  =(trefk(l)-tkl  )*tauk
        difql  =(qrefk(l)-qk(l))*tauk
        avrgtl =(tkl+tkl+diftl)
        dentpy =(diftl*cp+difql*elwv)*deta(l)/avrgtl+dentpy
        avrgt  =avrgtl*deta(l)+avrgt
        preck  =deta(l)*diftl+preck
        dift(l)=diftl
        difq(l)=difql
     END DO
     dentpy=dentpy+dentpy
     avrgt =avrgt/(sumdp+sumdp)

     !*** swap if entropy and/or precip .lt. 0 ..

     IF(dentpy.LT.epsntp.OR.preck.LT.d00)THEN
        IF(oct90)THEN
           cldefi(i,j)=efimn
        ELSE
           cldefi(i,j)=efimn*sm(i,j)+stefi*(h1-sm(i,j))
        ENDIF

        !*** search for shallow cloud top

        lbtk=lbot(i,j)
        ltsh=lbtk
        lbm1=lbtk-1
        pbtk=pk(lbtk)
        ptpk=pbtk-psh

        !*** cloud top is the level just below pbtk-psh
        DO l=1,lm
           IF(pk(l).LE.ptpk)ltpk=l+1
        END DO
        ptpk=pk(ltpk)

        !*** highest level allowed is level just below pshu

        IF(ptpk.LE.pshu)THEN
           DO l=1,lm
              IF(pk(l).LE.pshu)lshu=l+1
           END DO
           ltpk=lshu
           ptpk=pk(ltpk)
        ENDIF
        ltp1=ltpk+1
        ltp2=ltpk+2

        DO l=ltpk,lbtk
           qsatk(l)=pq0/pk(l)*EXP(a2*(tk(l)-a3)/(tk(l)-a4))
        END DO

        rhh=qk(ltpk)/qsatk(ltpk)
        rhmax=0.
        DO l=ltp1,lbm1
           rhl=qk(l)/qsatk(l)
           drhdp=(rhh-rhl)/(pk(l-1)-pk(l))
           IF(drhdp.GT.rhmax)THEN
              ltsh=l-1
              rhmax=drhdp
           ENDIF
           rhh=rhl
        END DO

        ltop(i,j)=ltsh

        !*** cloud must be at least two layers thick

        IF(lbot(i,j)-ltop(i,j).LT.2)ltop(i,j)=lbot(i,j)-2
        ptop(i,j)=pk(ltop(i,j))
        go to 600
     ENDIF

     !*** deep convection otherwise

     !***  keep the land value of efi equal to 1 until precip surpasses
     !***  a threshold value, currently set to 0.25 in per 24 hrs

     pthrs=cthrs/pdsl(i,j)
     drheat=(preck*sm(i,j)+amax1(epsp,preck-pthrs)*(h1-sm(i,j))) &
          *cp/avrgt
     efi=efifc*dentpy/drheat

     !*** unified or separate land/sea conv. 

     IF(oct90)THEN
        IF(unis)THEN
           efi=cldefi(i,j)*fcp+efi*fcc
        ELSEIF(.NOT.unil)THEN
           efi=(cldefi(i,j)*fcp+efi*fcc)*sm(i,j)+h1-sm(i,j)
        ELSE
           efi=h1
        ENDIF
     ELSE
        efi=cldefi(i,j)*fcp+efi*fcc
     ENDIF

     IF(efi.GT.h1   )    efi=h1
     IF(efi.LT.efimn)    efi=efimn
     IF(preck.EQ.0.)     efi=h1
     cldefi(i,j)=efi

     fefi=efmnt+slope*(efi-efimn)
     !     fefi=amax1(efi,efmnt)

     preck=preck*fefi

     !*** update precipitation, temperature & moisture

     precol(i,j)=pdsl(i,j)*preck*cprlg
     prec  (i,j)=pdsl(i,j)*preck*cprlg+prec  (i,j)
     cuprec(i,j)=pdsl(i,j)*preck*cprlg+cuprec(i,j)
     acprec(i,j)=pdsl(i,j)*preck*cprlg+acprec(i,j)
     cuppt(i,j) =pdsl(i,j)*preck*cprlg+cuppt(i,j)

     DO l=ltpk,lb
        tmod(i,j,l)=dift(l)*fefi
        qmod(i,j,l)=difq(l)*fefi
     END DO

     !*** end of deep convection
600  CONTINUE
  END DO

  ndeep=0
  DO j=3,jm-2
     DO i=1,im
        ltpk=ltop(i,j)
        lbtk=lbot(i,j)
        lb  =lmh(i,j)-1
        IF(ptop(i,j).LT.pbot(i,j)-psh)THEN
           ndeep=ndeep+1
           ndepth=lb-ltpk
           ntopd (ltpk  )=ntopd (ltpk  )+1
           nbotd (lb    )=nbotd (lb    )+1
           ndpthd(ndepth)=ndpthd(ndepth)+1
        ENDIF
     END DO
  END DO
  nneg=khdeep-ndeep

  !*** gather shallow convection points

  khshal=0
  ndstn =0
  ndstp =0
  DO j=3,jm-2
     DO i=1,im
        IF(ptop(i,j).GT.pbot(i,j)-pno.OR. &
             ltop(i,j).GT.lbot(i,j)-2) CYCLE
        IF(ptop(i,j)+h1.GE.pbot(i,j)-psh)THEN
           khshal=khshal+1
           ishal(khshal)=i
           jshal(khshal)=j
        ENDIF
     END DO
  END DO

  !*** horizontal loop for shallow convection

  DO n=1,khshal
     i=ishal(n)
     j=jshal(n)

     !*** shallow convection

     pz0=pd(i,j)+pt
     DO l=1,lm
        tkl      =t  (i,j,l)
        tk   (l) =tkl
        trefk(l) =tkl
        qkl      =q  (i,j,l)
        qk   (l) =qkl
        qrefk(l) =qkl
        pkl      =aeta(l)*pdsl(i,j)+pt
        pk   (l) =pkl
        qsatk(l) =pq0/pk(l)*EXP(a2*(tk(l)-a3)/(tk(l)-a4))
        apekl    =ape(i,j,l)

        !*** choose the pressure function
        !         fpk  (l) =alog(pkl)
        !         fpk  (l) =pkl
        !         fpk  (l) =-1./pkl

        apek (l) =apekl
        thvmkl   =tkl*apekl*(qkl*d608+h1)
        thvref(l)=thvmkl
     END DO

     !*** shallow cloud top

     lbtk=lbot(i,j)
     lbm1=lbtk-1
     ptpk=ptop(i,j)
     ltp1=ltop(i,j)
     ltpk=ltop(i,j)-1

     IF(ptop(i,j).GT.pbot(i,j)-pno.OR.ltop(i,j).GT.lbot(i,j)-2)THEN
        ltop(i,j)=lbot(i,j)
        ptop(i,j)=pbot(i,j)
        go to 800
     ENDIF

     !*** scaling potential temperature & table index at top

     !     thtpk=t(i,j,ltp1)*ape(i,j,ltp1)
     thtpk=t(i,j,ltpk)*ape(i,j,ltpk)

     tthk =(thtpk-thl)*rdth
     qqk  =tthk-AINT(tthk)
     it   =INT(tthk)+1
     IF(it.LT.1)THEN
        it=1
        qqk=d00
     ENDIF
     IF(it.GE.jtb)THEN
        it=jtb-1
        qqk=d00
     ENDIF

     !*** base and scaling factor for spec. humidity at top

     bqs00k=qs0(it)
     sqs00k=sqs(it)
     bqs10k=qs0(it+1)
     sqs10k=sqs(it+1)

     !*** scaling spec. humidity & table index at top

     bqk=(bqs10k-bqs00k)*qqk+bqs00k
     sqk=(sqs10k-sqs00k)*qqk+sqs00k
     !     tqk=(q(i,j,ltp1)-bqk)/sqk*rdq
     tqk=(q(i,j,ltpk)-bqk)/sqk*rdq

     ppk=tqk-AINT(tqk)
     iq =INT(tqk)+1
     IF(iq.LT.1)THEN
        iq=1
        ppk=d00
     ENDIF
     IF(iq.GE.itb)THEN
        iq=itb-1
        ppk=d00
     ENDIF

     !*** cloud top saturation point pressure

     ptpk= ptbl(iq,it  ) &
          +(ptbl(iq+1,it)-ptbl(iq,it))*ppk &
          +(ptbl(iq,it+1)-ptbl(iq,it))*qqk &
          +(ptbl(iq  ,it  )-ptbl(iq+1,it  ) &
          -ptbl(iq  ,it+1)+ptbl(iq+1,it+1))*ppk*qqk

     dpmix=ptpk-psp(i,j)
     IF(ABS(dpmix).LT.h3000)dpmix=-h3000

     !*** temperature profile slope

     smix=(thtpk-thbt(i,j))/dpmix*stabs

     trefkx=trefk(lbtk+1)
     pkxxxx=pk(lbtk+1)
     pkxxxy=pk(lbtk)
     apekxx=apek(lbtk+1)
     apekxy=apek(lbtk)
     DO l=ltp1,lbtk
        ivi=ltp1+lbtk-l
        trefkx=((pkxxxy-pkxxxx)*smix &
             +trefkx*apekxx)/apekxy
        trefk(ivi)=trefkx
        apekxx=apekxy
        pkxxxx=pkxxxy
        apekxy=apek(ivi-1)
        pkxxxy=pk(ivi-1)
     END DO

     !*** temperature reference profile correction

     sumdt=d00
     sumdp=d00

     DO l=ltp1,lbtk
        sumdt=(tk(l)-trefk(l))*deta(l)+sumdt
        sumdp=sumdp+deta(l)
     END DO

     rdpsum=1./sumdp
     fpk(lbtk)=trefk(lbtk)

     tcorr=sumdt*rdpsum

     DO l=ltp1,lbtk
        trfkl   =trefk(l)+tcorr
        trefk(l)=trfkl
        fpk  (l)=trfkl
     END DO

     !*** humidity profile equations

     psum  =0.
     qsum  =0.
     potsum=0.
     qotsum=0.
     otsum =0.
     dst   =0.
     fptk  =fpk(ltp1)

     DO l=ltp1,lbtk
        dpkl  =fpk(l)-fptk
        psum  =dpkl *deta(l)+psum
        qsum  =qk(l)*deta(l)+qsum
        rtbar =2./(trefk(l)+tk(l))
        otsum =deta(l)*rtbar+otsum
        potsum=dpkl   *rtbar*deta(l)+potsum
        qotsum=qk(l)  *rtbar*deta(l)+qotsum
        dst   =(trefk(l)-tk(l))*rtbar*deta(l)+dst
     END DO

     psum  =psum*rdpsum
     qsum  =qsum*rdpsum
     rotsum=1./otsum
     potsum=potsum*rotsum
     qotsum=qotsum*rotsum
     dst   =dst   *rotsum*cp/elwv

     !*** ensure positive entropy change

     IF(dst.GT.0.)THEN
        ltop(i,j)=lbot(i,j)
        ptop(i,j)=pbot(i,j)
        ndstp=ndstp+1
        go to 800
     ELSE
        dstq=dst*epsdn
     ENDIF

     !*** check for isothermal atmosphere

     den=potsum-psum
     IF(-den/psum.LT.5.e-5)THEN
        ltop(i,j)=lbot(i,j)
        ptop(i,j)=pbot(i,j)
        go to 800

        !*** slope of the reference humidity profile

     ELSE
        dqref=(qotsum-dstq-qsum)/den
     ENDIF
     !***  humidity does not increase with height

     IF(dqref.LT.0.)THEN
        ltop(i,j)=lbot(i,j)
        ptop(i,j)=pbot(i,j)
        go to 800
     ENDIF

     !*** humidity at the cloud top

     qrftp=qsum-dqref*psum

     !*** humidity profile

     DO l=ltp1,lbtk
        qrfkl=(fpk(l)-fptk)*dqref+qrftp

        !*** supersaturation or negative q not allowed

        qnew =(qrfkl-qk(l))*tauk+qk(l)
        IF(qnew.GT.qsatk(l)*stresh.OR.qnew.LT.d00)THEN
           ltop(i,j)=lbot(i,j)
           ptop(i,j)=pbot(i,j)
           go to 800
        ENDIF
        thvref(l)=trefk(l)*apek(l)*(qrfkl*d608+h1)
        qrefk(l)=qrfkl
     END DO

     !*** eliminate impossible slopes (betts, dtheta/dq)

     DO l=ltp1,lbtk
        dtdeta=(thvref(l-1)-thvref(l))/(aeta(l)-aeta(l-1))
        IF(dtdeta.LT.epsth)THEN
           ltop(i,j)=lbot(i,j)
           ptop(i,j)=pbot(i,j)
           go to 800
        ENDIF
     END DO

     !*** diagnostics

     IF(dst.GT.0.)THEN
        ndstp=ndstp+1
     ELSE
        ndstn=ndstn+1
     ENDIF

     dentpy=d00
     DO l=ltp1,lbtk
        dentpy=dentpy+((trefk(l)-tk(l))*cp+(qrefk(l)-qk(l))*elwv) &
             /(tk(l)+trefk(l))*deta(l)
     END DO

     !*** relaxation towards reference profiles

     DO l=ltp1,lbtk
        tmod(i,j,l)=(trefk(l)-tk(l))*tauk
        qmod(i,j,l)=(qrefk(l)-qk(l))*tauk
     END DO

     !*** end of shallow convection
800  CONTINUE
  END DO

  !*** diagnostics

  nshal=0
  DO j=3,jm-2
     DO i=1,im
        ltpk=ltop(i,j)
        lbtk=lbot(i,j)
        ptpk=ptop(i,j)
        pbtk=pbot(i,j)
        IF ( (ptpk <= pbtk-pno)  .AND. &
             (ltpk <= lbtk-2  )  .AND. &
             (ptpk >= pbtk-psh)          ) THEN
           nshal=nshal+1
           ntops(ltpk)=ntops(ltpk)+1
           nbots(lbtk)=nbots(lbtk)+1
           ndepth=lbtk-ltpk
           ndpths(ndepth)=ndpths(ndepth)+1
        ENDIF
     END DO
  END DO
  negds=khshal-nshal

  !*** smoothing temperature & humidity corrections

  IF(ksmud.EQ.0)THEN
     DO l=1,lm

        !***  update the fundamental prognostic arrays

        DO j=1,jm
           DO i=1,im
              t(i,j,l)=t(i,j,l)+tmod(i,j,l)
              q(i,j,l)=q(i,j,l)+qmod(i,j,l)

              !***  accumulate latent heating due to convection.
              !***  scale by the reciprocal of the period at which this routine
              !***  is called.  this period is the convection timestep.

              tcucn(i,j,l)=tcucn(i,j,l)+tmod(i,j,l)*rdtcnvc
           END DO
        END DO
     END DO
  ELSE
     DO l=1,lm
        DO j=1,jm
           DO i=1,im
              tl(i,j)=tmod(i,j,l)
              ql(i,j)=qmod(i,j,l)
           END DO
        END DO


        nsmud=ksmud
        DO ks=1,nsmud
           DO j=1,jm-1
              DO i=1,im
                 tne(i,j)=(tl(i+ihe(j),j+1)-tl(i,j)) &
                      *htm(i,j,l)*htm(i+ihe(j),j+1,l)
                 qne(i,j)=(ql(i+ihe(j),j+1)-ql(i,j)) &
                      *htm(i,j,l)*htm(i+ihe(j),j+1,l)
              END DO
           END DO
           DO j=2,jm
              DO i=1,im
                 tse(i,j)=(tl(i+ihe(j),j-1)-tl(i,j)) &
                      *htm(i+ihe(j),j-1,l)*htm(i,j,l)
                 qse(i,j)=(ql(i+ihe(j),j-1)-ql(i,j)) &
                      *htm(i+ihe(j),j-1,l)*htm(i,j,l)
              END DO
           END DO
           DO j=3,jm-2
              DO i=1,im
                 tl(i,j)=(tne(i,j)-tne(i+ihw(j),j-1)+tse(i,j)-tse(i+ihw(j),j+1)) &
                      *hbm2(i,j)*d125+tl(i,j)
                 ql(i,j)=(qne(i,j)-qne(i+ihw(j),j-1)+qse(i,j)-qse(i+ihw(j),j+1)) &
                      *hbm2(i,j)*d125+ql(i,j)
              END DO
           END DO
        END DO

        !***  update the fundamental prognostic arrays

        DO j=1,jm
           DO i=1,im
              t(i,j,l)=t(i,j,l)+tl(i,j)
              q(i,j,l)=q(i,j,l)+ql(i,j)
           END DO
        END DO

        !***  accumulate latent heating due to convection.
        !***  scale by the reciprocal of the period at which this routine
        !***  is called.  this period is the convection timestep.

        DO j=1,jm
           DO i=1,im
              tcucn(i,j,l)=tcucn(i,j,l)+tl(i,j)*rdtcnvc
           ENDDO
        ENDDO
     END DO
  ENDIF

  !*** save cloud top and bottom for radiation

  DO  j=1,jm
     DO  i=1,im
        htop(i,j)=MIN(float(ltop(i,j)),htop(i,j))
        hbot(i,j)=MAX(float(lbot(i,j)),hbot(i,j))
     END DO
  END DO

  !*** diagnostics 

  !     write(list,950)ntsd,nshal,ndeep,nneg,negds,ndstn,ndstp
  !         do 940 l=1,lm
  !     write(list,952)l
  !     write(list,954)nbots(l),ntops(l),ndpths(l), &
  !                  nbotd(l),ntopd(l),ndpthd(l)
  ! 940 continue
  ! 950 format(' ntsd=',i3,i8,' shallow, ',i4,' deep, ', &
  !     i4,' neg., ',i4,' neg. shall.,',i4,' dst.lt.0, ',i4,' dst.gt.0')
  ! 952 format(' layer (from top),',i2)
  ! 954 format('     nbots=',i4,'     ntops=',i4,'     ndpths=',i4,&
  !           '     nbotd=',i4,'     ntopd=',i4,'     ndpthd=',i4)
END SUBROUTINE cucnvc

