c****************************** plume **********************************
c
      subroutine plume ( ytflx, yqflx, yuflx, yvflx, yt, yq, yu, yv,
     *                   ylwc_c, ylwc_a, ylwc_s, 
     *                   fplumt, fplumq, fplumu, fplumv,
     *                   tracer, ypm, ypthic, ypmcap,
     *                   wplume0, wplume1, plumeprec, heimin, heimax,
     *                   dtplume)

c     Does convection (both pbl and free), by solving a subgrid
c     plume model. For pbl convection, the plume is initiated  
c     at the center of the lowest layer (ie, top of constant-flux layer)
c     with w,t,etc values scaled depending on the surface heat/momentum
c     fluxes based on constant-flux scaling, and is integrated 
c     vertically until it stops.

c     For free plumes, there is an option (commented by "csing" or 
c     "cmult") to do just one vertical integration from the bottom to
c     the top (csing) or to do multiplte integrations stating from
c     each layer in turn (cmult). In both cases, a free plume is 
c     initiated at the center of a model layer if it can rise to 
c     (or beyond) the center of the next layer. Plumes are integrated
c     vertically upwards from layer midpoint to midpoint until the 
c     plume stops. 

c     Condensation can occur within the plumes, forming wthin-plume
c     liquid cloud water, which can precipitate out if it exceeds a
c     maximum value.

c     The plume equations explicitly solve for plume w,theta,q,u,v where
c     all except w are perturbations from the ambient quantities. The
c     plume vertical finite-difference step is upwards from layer
c     midpoint to layer midpoint. These quantities, along with the plume
c     condensation amounts, then imply the large-scale fluxes at layer
c     midpoints and the large-scale precip and latent heat release.
c     These fluxes are accumulated over pbl and free
c     calculations, and then "shifted" downwards to the next lowest
c     layer interface to calculate large-scale convergences.
c     The large-scale flux convergences for heat, water vapor and 
c     momentum are actually put in reservoirs fplum* and released 
c     (see reserv, called from vdif) with a time scale of a few hours.
c     Liquid water cloud amounts and tracers are updated here with
c     no reservoir.

c     Lines beginning with "cloc" are to convert to standalone mode
c     for use with ./Plume/plume.f

c     supplied:
c     ytflx   = upward sensible heat flux from surface (W/m2)
c     yqflx   = upward h2o mass flux from surface (Kg/m2/s)
c     yuflx   = upward u-momentum flux from surface (N/m2)
c     yvflx   = upward v-momentum flux from surface (N/m2)
c     yt      = agcm temperatures (deg K)
c     yq      = agcm specific humidity (Kg/Kg)
c     yu      = agcm eastward  velocity (m/s)
c     yv      = agcm northward velocity  (m/s)
c     ypm     = agcm mid-layer pressures (N/m2)
c     ypthic  = agcm layer pressure-thicknesses (N/m2)
c     ypmcap  = (ypm/surface pressure)**cappa (N/m2)
c     dtplume = plume model timestep (s)

c     modified:
c     ylwc_c  = convective liquid water content (3-D)
c     ylwc_a  = anvil-cirrus liquid water content (3-D) 
c     ylwc_s  = stratiform   liquid water content (3-D)
c     fplum*  = reservoirs of "flux convergence*time" in each agcm layer
c     tracer  = agcm tracer fields. ylwc_[c,a,s] are equivalenced to
c               tracers #1-3, and are changed explicitly below.
c               fplum* are equivalenced to tracers #4-7. Only
c               tracers #8 to #ntrace are mixed as "tracers" below
c               (see local params ntraca,ntracb)
c     wplume0,1 = zonal mean plume vertical velocities (for diagnostics)
c     plumeprec = within-plume precipitation formation rate (3-D)
c                 (for reevap)
c     heimin    = daily min pbl plume heights (for history)
c     heimax    = daily max pbl plume heights (for history)

c ---------------------------------------------------------------------
      INCLUDE '../header.inc'
      INCLUDE 'complume'
c-----------------------------------------------------------------------

      dimension 
     *  ytflx(nlon,norec),        yqflx(nlon,norec),
     *  yuflx(nlon,norec),        yvflx(nlon,norec),

     *  yt(nlon,norec,nlev),      yq(nlon,norec,nlev),
     *  yu(nlon,norec,nlev),      yv(nlon,norec,nlev), 
     *  ylwc_c(nlon,norec,nlev),  ylwc_a(nlon,norec,nlev),
     *  ylwc_s(nlon,norec,nlev), 
     *  fplumt(nlon,norec,nlev),  fplumq(nlon,norec,nlev),
     *  fplumu(nlon,norec,nlev),  fplumv(nlon,norec,nlev), 
     *  tracer(nlon,norec,nlev,ntrace),

     *  ypm(nlon,norec,nlev),
     *  ypthic(nlon,norec,nlev),  ypmcap(nlon,norec,nlev),

     *  wplume0(norec,nlev),      wplume1(norec,nlev),
     *  plumeprec(nlon,norec,nlev),
     *  heimin(nlon,norec),       heimax(nlon,norec)

      dimension
     *  yh(nlon,nlev),         yrho(nlon,nlev),
     *  ylwc_t(nlon,nlev),

     *  dzu(nlon,nlev),        dzl(nlon,nlev),
     *  dztot(nlon,nlev),      dz(nlon,nlev),
     *  dhdz(nlon,nlev),       dqdz(nlon,nlev),
     *  dudz(nlon,nlev),       dvdz(nlon,nlev),
     *  dlwcdz(nlon,nlev),     dadz(nlon,nlev,ntrace),
     *  zelev(nlon,nlev)

      dimension
     *  zh(nlon,nlev),        zq(nlon,nlev), 
     *  zu(nlon,nlev),        zv(nlon,nlev),    
     *  zw(nlon,nlev),        zws(nlon,nlev),
     *  zt(nlon,nlev),        zrho(nlon,nlev),    
     *  zlw(nlon,nlev),       za(nlon,nlev,ntrace),
     *  zcond(nlon,nlev),     zlath(nlon,nlev),
     *  zprec(nlon,nlev),

     *  zconda(nlon,0:nlev),  zheata(nlon,0:nlev),
     *  zpreca(nlon,0:nlev),  zmfa(nlon,0:nlev),
     *  tplu(nlon,0:nlev),    qplu(nlon,0:nlev), 
     *  uplu(nlon,0:nlev),    vplu(nlon,0:nlev),
     *  cpluc(nlon,0:nlev),   cplua(nlon,0:nlev),
     *  cplus(nlon,0:nlev),   aplu(nlon,0:nlev,ntrace)

c     tracers 1-3 are used for liquid water contents (mixed by plume),
c     and tracers 4-7 are used for flux reservoirs (not mixed by plume).
      equivalence 
     *  (dlwcdz,  dadz(1,1,1)), 
     *  (zlw,     za(1,1,1)),
     *  (cpluc,   aplu(1,0,1)),
     *  (cplua,   aplu(1,0,2)),
     *  (cplus,   aplu(1,0,3))
      parameter (ntraca=8, ntracb=ntrace)

      dimension
     *  zarea(nlon,nlev),      zrad(nlon,nlev),
     *  zent(nlon),            zdrag(nlon),

     *  zcfl(nlon),            zmf(nlon,nlev),
     *  ztotc1(nlon),          ztotc2(nlon),
     *  ztota1(nlon),          ztota2(nlon),
     *  ztots1(nlon),          ztots2(nlon),
     *  zlwcold(nlon,nlev)

      dimension
     *  heipbl(nlon),          
     *  laypbl(nlon),           numpbl(nlev),
     *  laytop(nlon),           numfret(nlev),
     *  laybot(nlon),           numfreb(nlev),
     *  wfre1(nlon,nlev),       wfre2(nlon,nlev)

      logical ifpbl,            ifanv(nlon)
      character*3 chaplu(nlev)

c     Diagnostic only:
      dimension qrat(nlev), nqrat(nlev), qneg(nlev), nqneg(nlev)
c     parameter (nprinz=1)
      parameter (nprinz=2)                                         !cloc
      dimension iuprinz(nprinz)

      parameter (pi = 3.14159265358979)

c ---------------------------------------------------------------------
c     statement function for plume radius vs. elevation (not used):
c     plumerad (z) = plumerad1 + min (0.5*z, plumerad2)
c-----------------------------------------------------------------------
c#include "comsatplume"
      INCLUDE 'comsatplume'
c-----------------------------------------------------------------------

      cpair  = 1004.64
      gravit = 9.80616
      latvap = 2.5104e06
      latice = 3.336e5
      latsub = latvap + latice
      rair   = 287.04
      rh2o   = 4.61e2
      rhoh2o = 1.e3
      ch2o = 4.218e3
      cice = 2.106e3                  
      cappa  = rair/cpair
      zvir   = rh2o/rair -1.
      cpwv   = 1.81e3
      cpvir  = cpwv/cpair - 1.

      nstep = 0
      nrstrt = nstep
      twodt = 2.*dtime
      dtbud = dtime
      facbud = 1.

      alw   = .95
      blw_c = 1.
      blw_a = 1.
      blw_s = 1.
      clw_c = .00010
      clw_a = .00003
      clw_s = .00003
      clw_p = .0010
      dlw_l = 20.
      dlw_i = .5
      elw   = .3
      flw_c = 0.2
      flw_a = 0.5
      flw_s = 1.0
      glw = 0.3 ! (Smith, 1990, eq.2.29, c_a)
      tliqa = -15.
      tliqb =  10.
      siganvil = .45








c        Initialize local diagnostics

      zptot = 0.
      zetot = 0.
      zqtot = 0.
      do 10 jk=1,nlev
        numpbl(jk) = 0
        numfret(jk) = 0
        numfreb(jk) = 0
   10 continue

c        Overall loop over latitude

c***********************
      do 1000 jj=1,norec
c***********************

c         Calculate ambient pot.temp, density and total cloud liquid
c         water at layer midpoints

      do 100 jk=1,nlev
        do 102 ji=1,nlon
          yh(ji,jk)     = yt(ji,jj,jk) / ypmcap(ji,jj,jk)
          yrho(ji,jk)   = ypm(ji,jj,jk)
     *                  / (rair*(1.+zvir*yq(ji,jj,jk))*yt(ji,jj,jk))
          ylwc_t(ji,jk) = ylwc_c(ji,jj,jk) + ylwc_a(ji,jj,jk) 
     *                  + ylwc_s(ji,jj,jk)
  102   continue
  100 continue

c        Compute ambient dz, d*/dz for full layers. For bottom layer,
c        just use top half of layer.

      do 120 jk=1,nlev
        do 122 ji=1,nlon
          dztot(ji,jk) = ypthic(ji,jj,jk) / (gravit*yrho(ji,jk))
          dzu(ji,jk) = dztot(ji,jk) * (sig(jk)-sigkmh(jk))/dsigma(jk)
          dzl(ji,jk) = dztot(ji,jk) - dzu(ji,jk)
  122   continue
  120 continue

      do 124 jk=1,nlev-1
        do 126 ji=1,nlon
          dz(ji,jk) = dzl(ji,jk) + dzu(ji,jk+1)
          dhdz(ji,jk)   = (yh(ji,jk)     - yh(ji,jk+1)    ) / dz(ji,jk)
          dqdz(ji,jk)   = (yq(ji,jj,jk)  - yq(ji,jj,jk+1) ) / dz(ji,jk)
          dudz(ji,jk)   = (yu(ji,jj,jk)  - yu(ji,jj,jk+1) ) / dz(ji,jk)
          dvdz(ji,jk)   = (yv(ji,jj,jk)  - yv(ji,jj,jk+1) ) / dz(ji,jk)
          dlwcdz(ji,jk) = (ylwc_t(ji,jk) - ylwc_t(ji,jk+1)) / dz(ji,jk)
  126   continue

        if (ntraca.le.ntracb) then
          do 1260 n=ntraca,ntracb
            do 1262 ji=1,nlon
             dadz(ji,jk,n) = (tracer(ji,jj,jk,n)-tracer(ji,jj,jk+1,n))
     *                       /dz(ji,jk)
 1262       continue
 1260     continue
        endif
  124 continue

c        Compute elevations of layer midpoints (for plume radii calc)

      do 130 ji=1,nlon
        zelev(ji,nlev) = dzl(ji,nlev)
  130 continue
      do 132 jk=nlev-1,1,-1
        do 134 ji=1,nlon
          zelev(ji,jk) = zelev(ji,jk+1) + dz(ji,jk)
  134   continue
  132 continue

c        Zero plume fluxes accumulated over loop 500

      call zero (tplu,  nlon*nlevp)
      call zero (qplu,  nlon*nlevp)
      call zero (uplu,  nlon*nlevp)
      call zero (vplu,  nlon*nlevp)
      call zero (cpluc, nlon*nlevp)
      call zero (cplua, nlon*nlevp)
      call zero (cplus, nlon*nlevp)
      if (ntraca.le.ntracb)
     *  call zero (aplu(1,0,ntraca), nlon*nlevp*(ntracb-ntraca+1))
      call zero (zconda, nlon*nlevp)
      call zero (zheata, nlon*nlevp)
      call zero (zpreca, nlon*nlevp)
      call zero (zmfa,   nlon*nlevp)


c        Loop twice, for "pbl" plumes (ifpbl=.true.) and for "free"
c        plumes (ifpbl=.false., start at second lowest layer)

c*********************
c     do 500 itype=0,1
      do 500 itype=0,0   ! for pbl  only (diagnostics ok only for one)
c     do 500 itype=1,1   ! for free only (or the other)            !cloc
c*********************

      if (itype.eq.0) then
        ifpbl = .true.
        jkinita = nlev
        jkinitb = nlev
      else
        ifpbl = .false.
        jkinita = nlev  
        jkinitb = nlev                                            !csing
cmult   jkinitb = 2                                               !cmult
      endif

c         Initialize diagnostics for setting of plume velocities
c         (wfre*,wplume*) and indices of bottom,top layers (lay*, num*)

      call zero (wfre1, nlon*nlev)
      call zero (wfre2, nlon*nlev)

      do 140 ji=1,nlon
        heipbl(ji) = dzl(ji,nlev)
        laypbl(ji) = nlevp
        laytop(ji) = nlevp
        laybot(ji) = 0
  140 continue

c************************************
      do 510 jkinit = jkinita,jkinitb
c************************************

c        Assume no new plumes initiated above 150 mb

      if (sig(jkinit).lt.0.150) goto 510

c        Initialize all plume variables to zero or ambient, and
c        initialize diagnostic plume top/bottom locations 

      call zero (zw,   nlon*nlev)
      call zero (zws,  nlon*nlev)
      call zero (zt,   nlon*nlev)
      call zero (zh,   nlon*nlev)
      call zero (zq,   nlon*nlev)
      call zero (zu,   nlon*nlev)
      call zero (zv,   nlon*nlev)
      call zero (zlw,  nlon*nlev)
      call scopy (nlon*nlev, yrho, 1, zrho, 1)
      if (ntraca.le.ntracb)
     *  call zero (za(1,1,ntraca), nlon*nlev*(ntracb-ntraca+1))
      call zero (zarea,nlon*nlev)
      call zero (zcond, nlon*nlev)
      call zero (zlath, nlon*nlev)
      call zero (zprec, nlon*nlev)

c       For "pbl" plumes, initialize plume conditions
c       at midpoint of bottom layer based on surface fluxes from 
c       constant-flux layer. If net buoyancy flux 
c       (ytflx/cp+.622*yqflx) is < 0, leave initial plume z* = zero.

c--------------------
      if (ifpbl) then
c--------------------

        do 150 ji=1,nlon

c         zarea(ji,nlev) = plumearea1  ! .005
c         zrad(ji,nlev) = plumerad1    ! 100.
          zarea(ji,nlev) = xare1                                   !cloc
          zrad (ji,nlev) = xrad1                                   !cloc

          if ( ytflx(ji,jj)/cpair+.622*yqflx(ji,jj) .gt. 0. ) then

            ztau = sqrt (yuflx(ji,jj)**2 + yvflx(ji,jj)**2)
            ztauv = max (.01, sqrt(ztau/yrho(ji,nlev)) )
            zconv = ( (max(ytflx(ji,jj),0.)/(yrho(ji,nlev)*cpair))
     *                * rair * (1.-sig(nlev))
     *              ) ** (1./3.)
c-----
c old:
c           zwt = ytflx(ji,jj) / (cpair*yrho(ji,nlev)*zarea(ji,nlev))
c           zw(ji,nlev) = max (ztauv, (abs(zwt))**(1./3.))
c           zt(ji,nlev) = zwt / zw(ji,nlev)
c           if (abs(zt(ji,nlev)).gt.10.) then
c             zt(ji,nlev) = min (10., max(-10., zt(ji,nlev)))
c             zw(ji,nlev) = zwt / zt(ji,nlev)
c           endif
c           zmfbot = zarea(ji,nlev) * yrho(ji,nlev) * zw(ji,nlev)
c           zq(ji,nlev) = yqflx(ji,jj) / zmfbot
c           zu(ji,nlev) = yuflx(ji,jj) / zmfbot
c           zv(ji,nlev) = yvflx(ji,jj) / zmfbot
c-----
c new:
            zfac = 1. / zarea(ji,nlev)
            zwt = ytflx(ji,jj) / (yrho(ji,nlev)*cpair)
            zw(ji,nlev) = sqrt(zfac) * max(ztauv,zconv) 
            zt(ji,nlev) = zfac*zwt / zw(ji,nlev)
            if (abs(zt(ji,nlev)).gt.10.) then
              zt(ji,nlev) = min (10., max(-10., zt(ji,nlev)))
              zw(ji,nlev) = zfac*zwt / zt(ji,nlev)
            endif
            zq(ji,nlev)= zfac*yqflx(ji,jj)/(yrho(ji,nlev)*zw(ji,nlev))
            zu(ji,nlev)= zfac*yuflx(ji,jj)/(yrho(ji,nlev)*zw(ji,nlev))
            zv(ji,nlev)= zfac*yvflx(ji,jj)/(yrho(ji,nlev)*zw(ji,nlev))

            zh(ji,nlev) = zt(ji,nlev) / ypmcap(ji,jj,nlev)
            zws(ji,nlev)= zw(ji,nlev)**2
            zrho(ji,nlev)= ypm(ji,jj,nlev)
     *                  / ( rair*(1.+zvir*(zq(ji,nlev)+yq(ji,jj,nlev)))
     *                      *(zt(ji,nlev)+yt(ji,jj,nlev)) )

          endif

  150   continue

c---------
      else
c---------

c          For "free" plumes, set initial plume area and radius. 
c          Nb: plume radius affects the physics only via the 
c          entrainment and drag coefficients zent,zdrag. 
c          Plume area does not affect within-plume physics, only the
c          net flux convergences for changing large-scale quantities.
       
        jk = jkinit
        do 160 ji=1,nlon
          zarea(ji,jk) = xare2                                     !cloc
          zrad (ji,jk) = xrad2                                     !cloc
  160   continue

c       Free plume radius dep on column precipitable water vapor:
c       call zero (zpwv, nlon)
c       do 160 jk=1,nlev
c         do 162 ji=1,nlon
c           zpwv(ji) = zpwv(ji) + yq(ji,jj,jk)*ypthic(ji,jj,jk)
c 162     continue
c 160   continue
c       do 164 ji=1,nlon
c         zarea(ji,jkinit) = plumearea2
c         zrad(ji,jkinit) = max (plumerad1, plumerad2*(zpwv(ji)/gravit))
c 164   continue

c       do 160 ji=1,nlon
c         zarea(ji,jkinit) = plumearea2
c         zrad(ji,jkinit) = plumerad2
c 160   continue

c----------
      endif
c----------

c        Integrate plumes upwards from initiating layer to top

c================================
      do 200 jk = jkinit-1, 1, -1
c================================

        jkp = jk+1

c          Step plume from midpoint of layer jkp to midpoint of layer jk

        do 210 ji=1,nlon

c            Force for free plume, or non-zero pbl plume from below

          if( (.not.ifpbl)                   .or.zw(ji,jkp).gt.0.)!csing
cmult     if(((.not.ifpbl).and.jkp.eq.jkinit).or.zw(ji,jkp).gt.0.)!cmult
     *      then

            zent(ji)  = 0.20 / zrad(ji,jkp)
            zdrag(ji) = 0.5  / (pi*zrad(ji,jkp))

            zexp = exp(-zent(ji)*dz(ji,jk))
            zh(ji,jk) = zh(ji,jkp)*zexp
     *                  - (dhdz(ji,jk)/zent(ji))  *(1.-zexp)
            zq(ji,jk) = zq(ji,jkp)*zexp
     *                  - (dqdz(ji,jk)/zent(ji))  *(1.-zexp)
            zu(ji,jk) = zu(ji,jkp)*zexp
     *                  - (dudz(ji,jk)/zent(ji))  *(1.-zexp)
            zv(ji,jk) = zv(ji,jkp)*zexp
     *                  - (dvdz(ji,jk)/zent(ji))  *(1.-zexp)
            zlw(ji,jk)= zlw(ji,jkp)*zexp
     *                  - (dlwcdz(ji,jk)/zent(ji))*(1.-zexp)
            zt(ji,jk) = zh(ji,jk)*ypmcap(ji,jj,jk)

c              Test for parcel saturation, and condense amount that 
c              results in exactly saturated air (zcond, zlath).
c              Add to plume liquid water content, converting excess 
c              over maximum (namelist param clw_p) to within-plume
c              precip (zprec).
c              zlath is latent heat of condensation, determined by
c              ambient level temperature (ie, rain or snow), plus
c              correction for specific heats to be consistent with lsx
c              (imagining all phase changes occur at tmelt).

            ztabs = zt(ji,jk) + yt(ji,jj,jk)
            zqabs = zq(ji,jk) + yq(ji,jj,jk) 
            zqsat = qstblf (ztabs, ypm(ji,jj,jk))
            if (zqabs.gt.zqsat) then
c             zlath(ji,jk) = cvmgt ( latvap + cpwv*(ztabs-tmelt)
c    *                               - ch2o*(yt(ji,jj,jk)-tmelt),
c    *                               latsub + cpwv*(ztabs-tmelt)
c    *                               - cice*(yt(ji,jj,jk)-tmelt),
c    *                               yt(ji,jj,jk).ge.tmelt )
c             For now, all cloud water is liquid with no specific heat
              zlath(ji,jk) = latvap + cpwv*(ztabs-tmelt)
              dqz = (zlath(ji,jk)/rh2o) * zqsat / (ztabs**2)
              cpz = cpair*(1.+cpvir*zqabs)
              delq = (zqabs-zqsat) / (1.+dqz*zlath(ji,jk)/cpz)
              zqnew = zqabs - delq
              delt = delq*zlath(ji,jk) / (cpair*(1.+cpvir*zqnew))
              ztnew = ztabs + delt
              zq(ji,jk) = zqnew - yq(ji,jj,jk)
              zt(ji,jk) = ztnew - yt(ji,jj,jk)
              zh(ji,jk) = zt(ji,jk)/ypmcap(ji,jj,jk)

              zlwnew = zlw(ji,jk) + ylwc_t(ji,jk) + delq
c             zclw_p = clw_p                               ! (kg/kg dep)
              zclw_p = clw_p / yrho(ji,jk)                 ! (kg/m3 dep)
              delw = max (0., zlwnew-zclw_p)
              zlw(ji,jk) = min (zlwnew, zclw_p) - ylwc_t(ji,jk)
            else
              delq = 0.
              delt = 0.
              delw = 0.
            endif

            zrho(ji,jk) = ypm(ji,jj,jk)
     *                  / ( rair*(1.+zvir*(zq(ji,jk)+yq(ji,jj,jk)))
     *                      *(zt(ji,jk)+yt(ji,jj,jk)) )

            drhog0= gravit * (yrho(ji,jkp)-zrho(ji,jkp))/zrho(ji,jkp)
            drhog1= gravit * (yrho(ji,jk) -zrho(ji,jk) )/zrho(ji,jk) 
            a_rho =  (drhog0 + (drhog1-drhog0)/(1.-zexp)) / zent(ji)
            b_rho = -(         (drhog1-drhog0)/(1.-zexp)) / zent(ji)

            zexp2 = zexp*zexp
c           or: include vertical drag
c           zexp2 = exp(-2.*(zent(ji)+zdrag(ji))*dz(ji,jk))
c           a_rho = a_rho* zent(ji) / (zent(ji)+zdrag(ji))
c           b_rho = b_rho* zent(ji) / (2.*(zent(ji)+zdrag(ji))-zent(ji))

            zws(ji,jk) = zws(ji,jkp)*zexp2
     *                 + a_rho*(1.-zexp2) + 2.*b_rho*(zexp-zexp2)
            zw(ji,jk) = sqrt ( max(zws(ji,jk),0.) )

c             Test if detrainment >=0 is violated. If so, do corrections

            if (zw(ji,jk).gt.0.) then
              drho0 = yrho(ji,jkp)-zrho(ji,jkp)
              drho1 = yrho(ji,jk) -zrho(ji,jk)
              zwsmax = (   zrho(ji,jkp)*zws(ji,jkp)
     *                   + gravit*0.5*(drho0+drho1)*dz(ji,jk)
     *                 ) / zrho(ji,jk)
              if (zwsmax.gt.0. .and. zws(ji,jk).gt.zwsmax) then
                zws(ji,jk) = zwsmax
                zw(ji,jk) = sqrt(zws(ji,jk))
                zma = zrho(ji,jkp)*zw(ji,jkp)
                zmb = zrho(ji,jk) *zw(ji,jk)
                zent(ji) = (zmb-zma)/(0.5*(zma+zmb)*dz(ji,jk))
                zexp = exp(-zent(ji)*dz(ji,jk))
                zh(ji,jk) = zh(ji,jkp)*zexp 
     *                      - (dhdz(ji,jk)/zent(ji))  *(1.-zexp)
     *                      + delt/ypmcap(ji,jj,jk)
                zq(ji,jk) = zq(ji,jkp)*zexp
     *                      - (dqdz(ji,jk)/zent(ji))  *(1.-zexp)
     *                      - delq
                zu(ji,jk) = zu(ji,jkp)*zexp
     *                      - (dudz(ji,jk)/zent(ji))  *(1.-zexp)
                zv(ji,jk) = zv(ji,jkp)*zexp
     *                      - (dvdz(ji,jk)/zent(ji))  *(1.-zexp)
                zlw(ji,jk)= zlw(ji,jkp)*zexp
     *                      - (dlwcdz(ji,jk)/zent(ji))*(1.-zexp)
     *                      + delq - delw
                zt(ji,jk) = zh(ji,jk)*ypmcap(ji,jj,jk)
              endif
            endif

c              If vertical velocity = 0, end of plume;
c              else, re-compute zu and zv with (non-linear) form drag

            if (zw(ji,jk).eq.0.) then

              zh(ji,jk) = 0.
              zq(ji,jk) = 0.
              zu(ji,jk) = 0.
              zv(ji,jk) = 0.
              zlw(ji,jk) = 0.
              zt(ji,jk) = 0.
              zrho(ji,jk) = yrho(ji,jk)
              zcond(ji,jk) = 0.
              zlath(ji,jk) = 0.
              zprec(ji,jk) = 0.

c             Either:
c             heipbl(ji)= heipbl(ji) + dz(ji,jk) * zws(ji,jkp)
c    *                / (max(zws(ji,jkp),1.e-10) - min(zws(ji,jk),0.))
c             Or:
c             Solve quadratic eqn for exp(-zent*z),ie,for z where w=0,
c             for accurate increment to pbl height (diagnostic only).
c             Only valid for no vertical drag(if commented out above).
c             cvmgt(zz1) mod is fudge to avoid quadratic singularity.
              if (ifpbl) then
                zz1 = zws(ji,jkp) - a_rho - 2.*b_rho
                zz1 = cvmgt (zz1, zz1+1.e-9, abs(zz1).gt.1.e-10)
                zz2 = (-b_rho + sqrt(max(0.,b_rho**2-zz1*a_rho))) / zz1
                heipbl(ji) = heipbl(ji) - alog(zz2)/zent(ji)
                laypbl(ji) = jkp
              else
                if (zw(ji,jkp).gt.0.) laytop(ji) = min (laytop(ji),jkp)
              endif

              zws(ji,jk) = 0.

            else

              zenu= zent(ji) + zdrag(ji)*abs(0.5*(zu(ji,jkp)+zu(ji,jk)))
     *                         / (0.5*(zw(ji,jkp)+zw(ji,jk)))
              zenv= zent(ji) + zdrag(ji)*abs(0.5*(zv(ji,jkp)+zv(ji,jk)))
     *                         / (0.5*(zw(ji,jkp)+zw(ji,jk)))
              zexu  = exp(-zenu*dz(ji,jk))
              zexv  = exp(-zenv*dz(ji,jk))
              zu(ji,jk) = zu(ji,jkp)*zexu - (dudz(ji,jk)/zenu)*(1.-zexu)
              zv(ji,jk) = zv(ji,jkp)*zexv - (dvdz(ji,jk)/zenv)*(1.-zexv)

              zcond(ji,jk) = zarea(ji,jkp)*zrho(ji,jk)*zw(ji,jk)*delq
              zprec(ji,jk) = zarea(ji,jkp)*zrho(ji,jk)*zw(ji,jk)*delw

              if (ifpbl) then
                heipbl(ji) = heipbl(ji) + dz(ji,jk)
              else
                laybot(ji) = max (laybot(ji),jkp)
              endif

            endif

          endif

  210   continue

c          In same way, integrate plume tracers through curr layer

        if (ntraca.le.ntracb) then
          do 2100 n=ntraca,ntracb
            do 2102 ji=1,nlon
              if (zw(ji,jk).gt.0.) then
                zexp = exp(-zent(ji)*dz(ji,jk))
                za(ji,jk,n) = za(ji,jkp,n)*zexp
     *                      - (dadz(ji,jk,n)/zent(ji))*(1.-zexp)
              else
                za(ji,jk,n) = 0.
              endif
 2102       continue
 2100     continue
        endif

c          Set new plume radius and area.
c          Also count number of plumes, skip out of vertical 
c          integration for pbl or multiple-initiation plumes
c          if none left for this strip.

        nplu = 0
        do 230 ji=1,nlon
          zrad(ji,jk)  = zrad(ji,jkp)
          zarea(ji,jk) = zarea(ji,jkp)
c         zw(ji,jk)  = zw(ji,jk) * zarea(ji,jkp)/zarea(ji,jk)
c         zws(ji,jk) = zw(ji,jk)**2
          if (zw(ji,jk).gt.0.) nplu = nplu + 1
  230   continue
        if (ifpbl .and.  nplu.eq.0) goto 202                      !csing
cmult   if (             nplu.eq.0) goto 202                      !cmult

c=============
  200 continue
  202 continue
c=============
  
c        Accumulate large-scale fluxes at layer midpoints. 
c        Also accumulate mass flux zmfa for cfl check below

      do 300 jk=1,nlev

        do 302 ji=1,nlon
          zmf(ji,jk)  = zarea(ji,jk)*zrho(ji,jk)*zw(ji,jk)
          zmfa(ji,jk) = zmfa(ji,jk) + zmf(ji,jk)
  302   continue

        if (ifpbl .and. jk.eq.nlev) then
          do 304 ji=1,nlon
            tplu(ji,jk) = tplu(ji,jk) + ytflx(ji,jj)
            qplu(ji,jk) = qplu(ji,jk) + yqflx(ji,jj)
            uplu(ji,jk) = uplu(ji,jk) + yuflx(ji,jj)
            vplu(ji,jk) = vplu(ji,jk) + yvflx(ji,jj)
  304     continue
        else
          do 306 ji=1,nlon
            tplu(ji,jk) = tplu(ji,jk) + zt(ji,jk)*zmf(ji,jk)*cpair
            qplu(ji,jk) = qplu(ji,jk) + zq(ji,jk)*zmf(ji,jk)
c           comment out next two lines to turn off "cumulus friction"
            uplu(ji,jk) = uplu(ji,jk) + zu(ji,jk)*zmf(ji,jk)
            vplu(ji,jk) = vplu(ji,jk) + zv(ji,jk)*zmf(ji,jk)

            cpluc(ji,jk) = cpluc(ji,jk)
     *                   + ( zlw(ji,jk) + ylwc_t(ji,jk)
     *                                  - ylwc_c(ji,jj,jk) )*zmf(ji,jk)
            cplua(ji,jk) = cplua(ji,jk) - ylwc_a(ji,jj,jk)  *zmf(ji,jk)
            cplus(ji,jk) = cplus(ji,jk) - ylwc_s(ji,jj,jk)  *zmf(ji,jk)

            zconda(ji,jk)= zconda(ji,jk) + zcond(ji,jk)
            zheata(ji,jk)= zheata(ji,jk) + zcond(ji,jk)*zlath(ji,jk)
            zpreca(ji,jk)= zpreca(ji,jk) + zprec(ji,jk)
  306     continue
        endif

        if (ntraca.le.ntracb) then
          do 3060 n=ntraca,ntracb
            do 3062 ji=1,nlon
              aplu(ji,jk,n) = aplu(ji,jk,n) + za(ji,jk,n)*zmf(ji,jk)
 3062       continue
 3060     continue
        endif

  300 continue

c        For "pbl" plumes,set history daily min/max pbl heights,
c        accumulate diagnostic number of pbl-tops in each layer(numpbl),
c        and save diagnostic zonal mean plume vertical velocities in
c        wplume0. For "free" plumes, accumulate wfre1,2 for later
c        setting (after all initiating layers done) of wplume1.

c--------------------
      if (ifpbl) then
c--------------------

        do 400 ji=1,nlon
          heimin(ji,jj) = min (heipbl(ji), heimin(ji,jj))
          heimax(ji,jj) = max (heipbl(ji), heimax(ji,jj))
          if (laypbl(ji).ne.nlevp)
     *      numpbl(laypbl(ji)) = numpbl(laypbl(ji)) + 1
  400   continue

        do 402 jk=1,nlev
          wplume0(jj,jk) = 0.
          nzw = 0
          do 404 ji=1,nlon
            if (zw(ji,jk).gt.0.) then
              wplume0(jj,jk) = wplume0(jj,jk) + zw(ji,jk)
              nzw = nzw + 1
            endif
  404     continue
          if (nzw.gt.0) wplume0(jj,jk) = wplume0(jj,jk)/nzw
  402   continue

c---------
      else
c---------

        do 412 jk=1,nlev
          do 414 ji=1,nlon
c           wfre1(ji,jk) = wfre1(ji,jk) + zarea(ji,jk)*zrho(ji,jk)
c    *                                    *zw(ji,jk)
c           wfre2(ji,jk) = wfre2(ji,jk) + zarea(ji,jk)*zrho(ji,jk)
            wfre1(ji,jk) = max (wfre1(ji,jk),zw(ji,jk))
            wfre2(ji,jk) = 1.
  414     continue
  412   continue

c----------
      endif
c----------

c        End of loop over initiating layers

c*************
  510 continue
c*************

c        Set diagnostic zonal "free" plume bottom and top locations
c        and vertical velocity wplume1

      if (.not.ifpbl) then
  
        do 520 ji=1,nlon
          if (laytop(ji).ne.nlevp)
     *      numfret(laytop(ji)) = numfret(laytop(ji)) + 1
          if (laybot(ji).ne.0) 
     *      numfreb(laybot(ji)) = numfreb(laybot(ji)) + 1
  520   continue

        do 530 jk=1,nlev
          wplume1(jj,jk) = 0.
          nzw = 0
          do 532 ji=1,nlon
            if (wfre1(ji,jk).gt.0. .and. wfre2(ji,jk).gt.0.) then
              nzw = nzw + 1
              wplume1(jj,jk)= wplume1(jj,jk) + wfre1(ji,jk)/wfre2(ji,jk)
            endif
  532     continue
          if (nzw.gt.0) wplume1(jj,jk) = wplume1(jj,jk)/nzw
  530   continue

      endif

c        End of "pbl" or "free" loop

c*************
  500 continue
c*************

c        Check CFL mass-flux  criterion, find worst violation for each
c        column, and reduce fluxes uniformly for a column to make CFL
c        ok. (But not for bottom layer, otherwise wouldn't conserve
c        surface fluxes).
 
      call zero (zcfl, nlon)

      do 600 jk=1,nlev-1
        do 602 ji=1,nlon
          zcfl(ji) = max ( zcfl(ji), 
     *                     zmfa(ji,jk)*dtplume*gravit
     *                     / min (ypthic(ji,jj,jk),ypthic(ji,jj,jk+1))
     *                   ) 
  602   continue
  600 continue
 
c     cflmax = 1.e6
      cflmax = 0.75

      do 610 jk=1,nlev-1
        do 612 ji=1,nlon
          zfac = cflmax / max (cflmax, zcfl(ji))
          zmfa(ji,jk)  = zmfa(ji,jk)  * zfac
          tplu(ji,jk)  = tplu(ji,jk)  * zfac
          qplu(ji,jk)  = qplu(ji,jk)  * zfac
          uplu(ji,jk)  = uplu(ji,jk)  * zfac
          vplu(ji,jk)  = vplu(ji,jk)  * zfac
          cpluc(ji,jk) = cpluc(ji,jk) * zfac
          cplua(ji,jk) = cplua(ji,jk) * zfac
          cplus(ji,jk) = cplus(ji,jk) * zfac
          zconda(ji,jk)= zconda(ji,jk)* zfac
          zheata(ji,jk)= zheata(ji,jk)* zfac
          zpreca(ji,jk)= zpreca(ji,jk)* zfac
  612   continue

        if (ntraca.le.ntracb) then
          do 6120 n=ntraca,ntracb
            do 6122 ji=1,nlon
              zfac = cflmax / max (cflmax, zcfl(ji))
              aplu(ji,jk,n) = aplu(ji,jk,n) * zfac
 6122       continue
 6120     continue
        endif
  610 continue 

c     zz = 0.
c     do 8880 ji=1,nlon
c       zz = zz +  cflmax / max (cflmax, zcfl(ji))
c8880 continue
c     write(*,8881) nstep, jj, ji, cflmax, zz/nlon
c8881 format('plume: zonal mean cfl: nstep,jj,cflmax,zfac=',2i4,2f6.3)

c        Compute flux convergence into each layer and add to reservoirs
c        fplum* ("flux*time per layer", in commun) for use in vdif.
c        Also increment tracer fields themselves (no reservoir).
c        Use dtplume since reservoirs are not leapfrog variables
c        (see comments in reserv).

c        For simplicity, shift fluxes at layer midpoints down to next
c        lowest interface, (so give zconda,zheata,zpreca to layer below)
 
      do 700 jk=1,nlev
        do 702 ji=1,nlon
          fplumt(ji,jj,jk) = fplumt(ji,jj,jk)
     *                     + (  tplu(ji,jk)-tplu(ji,jk-1) 
     *                        + zheata(ji,jk-1)           ) * dtplume
          fplumq(ji,jj,jk) = fplumq(ji,jj,jk)
     *                     + (  qplu(ji,jk)-qplu(ji,jk-1)
     *                        - zconda(ji,jk-1)           ) * dtplume
          fplumu(ji,jj,jk) = fplumu(ji,jj,jk)
     *                     + (  uplu(ji,jk)-uplu(ji,jk-1) ) * dtplume
          fplumv(ji,jj,jk) = fplumv(ji,jj,jk)
     *                     + (  vplu(ji,jk)-vplu(ji,jk-1) ) * dtplume

          zqtot = zqtot
     *          + (qplu(ji,jk)-qplu(ji,jk-1))*dtplume*cosbud(jj)/nlon
          zptot = zptot + zconda(ji,jk-1)*dtplume*cosbud(jj)/nlon
          if (ifpbl .and. jk.eq.nlev)
     *      zetot = zetot + yqflx(ji,jj)*dtplume*cosbud(jj)/nlon
  702   continue

        if (ntraca.le.ntracb) then
          do 7020 n=ntraca,ntracb
            do 7022 ji=1,nlon
              zm = ypthic(ji,jj,jk)/gravit
              tracer(ji,jj,jk,n) = tracer(ji,jj,jk,n)
     *                      + (aplu(ji,jk,n)-aplu(ji,jk-1,n))*dtplume/zm
 7022       continue
 7020     continue
        endif
  700 continue

c        Increment large-scale liquid water contents (ylwc_*, kg/kg)
c        set within-plume precipitation (plumeprec, kg/m2/s), increment
c        global budget terms (d[m,h]clouf, kg/m2)

      do 720 ji=1,nlon
        ifanv(ji) = .true.
  720 continue

      do 730 jk=1,nlev                                   ! top to bottom
        do 732 ji=1,nlon

          zm = ypthic(ji,jj,jk)/gravit
          zc = cpluc (ji,jk)   - cpluc(ji,jk-1)
     *       + zconda(ji,jk-1) - zpreca(ji,jk-1)
          zlwcold(ji,jk) = ylwc_c(ji,jj,jk)

c         Give to anvil if plume top is high enough, else to conv
          if (zc.ne.0.) then
            if (ifanv(ji) .and. sig(jk).lt.siganvil) then 
              ylwc_a(ji,jj,jk)= ylwc_a(ji,jj,jk) + dtplume*zc/zm
            else
              ylwc_c(ji,jj,jk)= ylwc_c(ji,jj,jk) + dtplume*zc/zm
            endif
            ifanv(ji) = .false.
          endif

          ylwc_a(ji,jj,jk) = ylwc_a(ji,jj,jk)
     *                     + (cplua(ji,jk)-cplua(ji,jk-1)) * dtplume/zm
          ylwc_s(ji,jj,jk) = ylwc_s(ji,jj,jk)
     *                     + (cplus(ji,jk)-cplus(ji,jk-1)) * dtplume/zm

          plumeprec(ji,jj,jk) = zpreca(ji,jk-1)

          dmclouf = dmclouf + zconda(ji,jk-1)*(cosbud(jj)/nlon)*dtplume
          dhcloud = dhcloud + zheata(ji,jk-1)*(cosbud(jj)/nlon)*dtplume
  732   continue
  730 continue

c        Crudely fix very slightly negative cloud liquid water contents
c        (otherwise causes blowups in cldcmp/radctl)

      call zero (ztotc1, nlon)
      call zero (ztota1, nlon)
      call zero (ztots1, nlon)
      call zero (ztotc2, nlon)
      call zero (ztota2, nlon)
      call zero (ztots2, nlon)
      do 750 jk=1,nlev
        do 752 ji=1,nlon
          zm = ypthic(ji,jj,jk)/gravit
          ztotc1(ji) = ztotc1(ji) + ylwc_c(ji,jj,jk)*zm
          ztota1(ji) = ztota1(ji) + ylwc_a(ji,jj,jk)*zm
          ztots1(ji) = ztots1(ji) + ylwc_s(ji,jj,jk)*zm
          ylwc_c(ji,jj,jk) = max (ylwc_c(ji,jj,jk), 0.) 
          ylwc_a(ji,jj,jk) = max (ylwc_a(ji,jj,jk), 0.) 
          ylwc_s(ji,jj,jk) = max (ylwc_s(ji,jj,jk), 0.) 
          ztotc2(ji) = ztotc2(ji) + ylwc_c(ji,jj,jk)*zm
          ztota2(ji) = ztota2(ji) + ylwc_a(ji,jj,jk)*zm
          ztots2(ji) = ztots2(ji) + ylwc_s(ji,jj,jk)*zm
  752   continue
  750 continue

      do 754 ji=1,nlon
c       Since z*2 > z*1, if z*1>0, so is z*2:
        ztotc1(ji) = cvmgt ( ztotc1(ji)/max(ztotc2(ji),1.e-20), 0.,
     *                       ztotc1(ji).gt.0. )
        ztota1(ji) = cvmgt ( ztota1(ji)/max(ztota2(ji),1.e-20), 0.,
     *                       ztota1(ji).gt.0.)
        ztots1(ji) = cvmgt ( ztots1(ji)/max(ztots2(ji),1.e-20), 0.,
     *                       ztots1(ji).gt.0.)
  754 continue

      do 756 jk=1,nlev
        do 758 ji=1,nlon
          ylwc_c(ji,jj,jk) = ztotc1(ji) * ylwc_c(ji,jj,jk)
          ylwc_a(ji,jj,jk) = ztota1(ji) * ylwc_a(ji,jj,jk)
          ylwc_s(ji,jj,jk) = ztots1(ji) * ylwc_s(ji,jj,jk)
  758   continue
  756 continue

c*****
c     Diagnostic printout: cloud liquid water
c     do 8882 ji=1,nlon
c       do 8883 jk=1,nlev
c         if (ylwc_c(ji,jj,jk).lt.0.) then
c           write (*,*) 'plume: ji=',ji,' jj=',jj,' jk=',jk
c           write (*,8887)
c           do 8884 k=1,nlev 
c             zm = ypthic(ji,jj,k)/gravit
c             zca = (cpluc(ji,k)-cpluc(ji,k-1))*dtplume/zm
c             zcb = zconda(ji,k-1)*dtplume/zm
c             zcc = -zpreca(ji,k-1)*dtplume/zm
c             write(*,8888) k, 
c    *                      zw(ji,k),
c    *                      zmf(ji,k)*dtplume*gravit/ypthic(ji,jj,k),
c    *                      ylwc_t(ji,k)*1.e3,
c    *                      (zlw(ji,k)+ylwc_t(ji,k))*1.e3,
c    *                      zlwcold(ji,k)*1.e3,
c    *                      ylwc_c(ji,jj,k)*1.e3,
c    *                      zca*1.e3,
c    *                      zcb*1.e3,
c    *                      zcc*1.e3
c8887         format(' k', 2x, '         zw', '        zmf',
c    *                        '      ylwc_t', '        zlw',
c    *                        '    ylw_cold', '    ylw_cnew',
c    *                        '     d(cplu)',
c    *                        '   d(zconda)', '   d(zpreca)')
c8888         format(i2, 2x, 9f12.8)
c8884       continue
c           call endrun
c         endif
c8883   continue
c8882 continue
c*****

c     iuprinz(1) = nout
      iuprinz(1) = 6
      iuprinz(2) = iuout

c*****
c     Diagnostic printout: vertical profile at one point
c     nz = nint(86400./dtime)
c     if (nstep.gt.nstop-nz .and. jj.eq.norec/2) then
c       zzmax = -1.e20
c       do 8000 ji=1,nlon
c         zz = 0.
c         do 8002 jk=1,nlev
c           if (zconda(ji,jk-1).gt.0.) zz = zz + 1.
c8002     continue
c         if (zz.gt.zzmax) then
c           jim =ji
c           zzmax = zz
c         endif
c8000   continue

c       if (jj.eq.14) then
c         jim = 18
        if (jj.eq.1) then                                          !cloc
          jim = 1                                                  !cloc
          do 8010 ipz=1,nprinz
c           write(iuprinz(ipz),8012) nstep*dtime/86400.,jim,jj
c8012       format(/' day=',f8.3,'  ji=',i3,'  jj=',i3)
            write(iuprinz(ipz),*)
            write(iuprinz(ipz),8014) 
     *       (jk, sigkmh(jk),
     *            zarea(jim,jk), zw(jim,jk),
     *            zt(jim,jk),    zq(jim,jk), zrho(jim,jk),
     *            yh(jim,jk),    yt(jim,jj,jk),
     *            yq(jim,jj,jk), yrho(jim,jk),
     *            tplu(jim,jk),  qplu(jim,jk)*.864e5,
     *            zconda(jim,jk)*.864e5, 
     *            fplumt(jim,jj,jk)/(cpair*ypthic(jim,jj,jk)/gravit),
     *            fplumq(jim,jj,jk)/(ypthic(jim,jj,jk)/gravit),
c    *      jk=nlev,1,-1)
     *      jk=1,nlev)
 8014       format(' jk',4x,'sig',
     *             5x,'zarea',8x,'zw',
     *             8x,'zt',8x,'zq',6x,'zrho',
     *             8x,'yh',8x,'yt',
     *             8x,'yq',6x,'yrho',
     *             6x,'tplu',6x,'qplu',
     *             4x,'zconda',
     *             4x,'fplumt',4x,'fplumq'
     *           /(i3,f7.3,14f10.5))
 8010     continue
        endif
c     endif
c*****

c*****
c     Diagnostic printout: one line for one point
c     nz = nint(86400./dtime)
c     if (nstep.gt.nstop-2*nz .and. nstep.le.nstop-nz.and.jj.eq.norec/2)
c       then
c       if (mod(nstep,nz).eq.1) write(iuout,8020) ji,jj
c       ji = nlon/2
      if (jj.eq.1) then                                            !cloc
        ji = 1                                                     !cloc
 8020   format(/' ji,jj=',2i3)
        do 8022 jk=1,nlev
          chaplu(jk) = '   '
 8022   continue
        chaplu(laytop(ji)) = '***'
 
        do 8030 ipz=1,nprinz
c         write(iuprinz(ipz),8032)
c8032     format(5x,'day',2x,'heipbl',3x,'ytflx',
c    *           5x,5(6x,'yh',1x),5x,5(6x,'yq',1x) )
c         write(iuprinz(ipz),8034) nstep*dtime/86400., 
c    *      heipbl(ji), ytflx(ji,jj),
c    *      (yh(ji,jk),chaplu(jk)(1:1), 
c    *      jk=nlev,nlev-4,-1),
c    *      (yq(ji,jj,jk),chaplu(jk)(1:1), jk=nlev,nlev-4,-1)
c8034     format(f8.3,f8.1,f8.1,5x,5(f8.2,a1),5x,5(f8.5,a1))
          write(iuprinz(ipz),'(a,f8.1)') 'heipbl=',heipbl(ji)
 8030    continue
      endif
c*****

c        End of overall loop over latitude

c*************
 1000 continue
c*************

c*****
c     Diagnostic printout: number of negative q's globally vs level
      zz = dtime / (3600.)
      qnegmax = 1.e20
      nqnegtot = 0
      do 8040 jk=1,nlev
        qrat(jk) = 0.
        qneg(jk) = 0.
        nqrat(jk) = 0
        nqneg(jk) = 0
        do 8042 jj=1,norec
        do 8042 ji=1,nlon
          dqnew = zz * fplumq(ji,jj,jk) / (ypthic(ji,jj,jk)/gravit)
          if (abs(yq(ji,jj,jk)).gt.0.00005) then
            qrat(jk) = qrat(jk) + abs(dqnew) / abs(yq(ji,jj,jk))
            nqrat(jk) = nqrat(jk) + 1
          endif
          qnew = yq(ji,jj,jk) + dqnew
          if (qnew.lt.0.) then
            qneg(jk) = qneg(jk) + qnew
            nqneg(jk) = nqneg(jk) + 1
            nqnegtot = nqnegtot + 1
            if (qnew.lt.qnegmax) then
              qnegmax = qnew
              jkmax = jk
              jimax = ji
              jjmax = jj
            endif
          endif
 8042   continue
        if (nqrat(jk).gt.0) qrat(jk) = qrat(jk)/nqrat(jk)
        if( nqneg(jk).gt.0) qneg(jk) = qneg(jk)/nqneg(jk)
 8040 continue
 
      if (nqnegtot.gt.0) then
        do 8050 ipz=1,nprinz
          write(iuprinz(ipz),8052)
     *       nstep*dtime/86400., zptot, zetot, zqtot,
     *       qnegmax, jimax,jjmax,jkmax,
     *       nqrat(jkmax),qrat(jkmax),nqneg(jkmax),qneg(jkmax)
 8052     format(/' day=',f8.3,
     *            '   zptot=',f10.5,'   zetot=',f10.5,'   zqtot=',f10.5
     *           /' qnegmax=',f10.5,'  at=',3i4
     *           /' nqrat=',i4,'   qrat=',f9.6,'  nqneg=',i4,
     *            '  qneg=',f9.6)
 8050   continue
      endif
c*****

c*****
c     Diagnostic printout: max heimax, numfret
c     nz = nint(86400./dtime)
c     if (nstep.gt.nstop-nz) then
c       zz = -1.e20
c       do 8060 jj=1,norec
c       do 8060 ji=1,nlon
c         if (heimax(ji,jj).gt.zz) then
c           izz = ji
c           jzz = jj
c           zz = heimax(ji,jj)
c         endif
c8060   continue
c       write(*,8062) nstep*dtime/86400., zz, izz, jzz,
c    *                (numfret(jk),jk=nlev,1,-1) 
c8062   format(' day=',f8.3,'  heimax=',f8.1,'  at ji,jj=',2i3,
c    *         '  numfret=',(18i5))
c     endif
c*****
      
c*****
c     Diagnostic printout: number of pbl and plume tops in each layer
c     (written to fort.87) 
      if (mod(nstep,nint(1.*86400./dtime)) .eq. 0) then
        iu = 87
        write(iu,8070) nstep*dtime/86400.
        write(iu,8072) 'pbl  ', (numpbl (jk),jk=nlev,1,-1)
        write(iu,8072) 'fret ', (numfret(jk),jk=nlev,1,-1)
        write(iu,8072) 'freb ', (numfreb(jk),jk=nlev,1,-1)
        ier = flush (iu)
 8070   format('day: ',f8.3)
 8072   format(4x,'num',a,'=',(18i5))
      endif
c*****

c     DO J=1,ND
c       IF (ZMF(1,J) .NE. 0.0) THEN
c         YT(1,1,J) = YT(1,1,J) + TPLU(1,J) / (ZMF(1,J)*CPAIR)
c         YQ(1,1,J) = YQ(1,1,J) + QPLU(1,J) / (ZMF(1,J))
c       END IF
c     END DO


      return
      end
