moistconvection.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#if ( defined DIAGNS )
! Determine chunk latitudes and longitudes
   call get_lat_all_p(lchnk, ncol, lat)
   call get_lon_all_p(lchnk, ncol, lon)
#endif
!
! Ensure that characteristic adjustment time scale (cmftau) assumed
! in estimate of eta isn't smaller than model time scale (ztodt)
! The time over which the convection is assumed to act (the adjustment
! time scale) can be applied with each application of the three-level
! cloud model, or applied to the column tendencies after a "hard"
! adjustment (i.e., on a 2-delta t time scale) is evaluated
!
   if (rlxclm) then
      cats   = ztodt             ! relaxation applied to column
      adjfac = ztodt/(max(ztodt,cmftau))
   else
      cats   = max(ztodt,cmftau) ! relaxation applied to triplet
      adjfac = 1.0
   endif
   rtdt = 1.0/ztodt
!
! Move temperature and moisture into working storage
!
   do k=limcnv,pver
      do i=1,ncol
         tb (i,k) = t(i,k)
         shb(i,k) = q(i,k,1)
      end do
   end do
   do k=1,pver
      do i=1,ncol
         icwmr(i,k) = 0.
      end do
   end do
!
! Compute sb,hb,shbs,hbs
!
   call aqsatd(tb      ,pmid    ,estemp ,shbs    ,gam     , &
               pcols   ,ncol    ,pver   ,limcnv  ,pver    )
!
   do k=limcnv,pver
      do i=1,ncol
         sb (i,k) = cp*tb(i,k) + zm(i,k)*grav + phis(i)
         hb (i,k) = sb(i,k) + hlat*shb(i,k)
         hbs(i,k) = sb(i,k) + hlat*shbs(i,k)
      end do
   end do
!
! Compute sbh, shbh
!
   do k=limcnv+1,pver
      do i=1,ncol
         sbh (i,k) = 0.5*(sb(i,k-1) + sb(i,k))
         shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
         hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
      end do
   end do
!
! Specify properties at top of model (not used, but filling anyway)
!
   do i=1,ncol
      sbh (i,limcnv) = sb(i,limcnv)
      shbh(i,limcnv) = shb(i,limcnv)
      hbh (i,limcnv) = hb(i,limcnv)
   end do
!
! Zero vertically independent control, tendency & diagnostic arrays
!
   do i=1,ncol
      prec(i)  = 0.0
      dzcld(i) = 0.0
      cnb(i)   = 0.0
      cnt(i)   = float(pver+1)
   end do
#if ( defined DIAGNS )
!#######################################################################
!#                                                                     #
!#    output initial thermodynamic profile if debug diagnostics        #
!#                                                                     #
   do i=1,ncol
     if ((lat(i).eq.jloc) .and. (lon(i).eq.iloc) &
         .and. (nstep.ge.nsloc)) then
!#                                                                     #
!#       approximate vertical integral of moist static energy          #
!#       and total preciptable water                                   #
!#                                                                     #
      hsum1 = 0.0
      qsum1 = 0.0
      do k=limcnv,pver
         hsum1 = hsum1 + pdel(i,k)*rgrav*hb(i,k)
         qsum1 = qsum1 + pdel(i,k)*rgrav*shb(i,k)
      end do
!#                                                                     #
      write (6,8010)
      fac = grav*864.
      do k=limcnv,pver
         rh = shb(i,k)/shbs(i,k)
         write(6,8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k),cmfsl(i,k), cmflq(i,k)
         write(6,8040) tb(i,k),shb(i,k),rh,sb(i,k),hb(i,k),hbs(i,k),ztodt*cmfdt(i,k), &
                       ztodt*cmfdq(i,k),ztodt*cmfdqr(i,k)
      end do
      write(6, 8000) prec(i)
     end if
   enddo
#endif
!#                                                                     #
!#                                                                     #
!#######################################################################
!
! Begin moist convective mass flux adjustment procedure.
! Formalism ensures that negative cloud liquid water can never occur
!
   do 70 k=pver-1,limcnv+1,-1
      do 10 i=1,ncol
         etagdt(i) = 0.0
         eta   (i) = 0.0
         beta  (i) = 0.0
         ds1   (i) = 0.0
         ds2   (i) = 0.0
         ds3   (i) = 0.0
         dq1   (i) = 0.0
         dq2   (i) = 0.0
         dq3   (i) = 0.0
!
! Specification of "cloud base" conditions
!
         qprime    = 0.0
         tprime    = 0.0
!
! Assign tprime within the PBL to be proportional to the quantity
! tpert (which will be bounded by tpmax), passed to this routine by
! the PBL routine.  Don't allow perturbation to produce a dry
! adiabatically unstable parcel.  Assign qprime within the PBL to be
! an appropriately modified value of the quantity qpert (which will be
! bounded by shpmax) passed to this routine by the PBL routine.  The
! quantity qprime should be less than the local saturation value
! (qsattp=qsat[t+tprime,p]).  In both cases, tpert and qpert are
! linearly reduced toward zero as the PBL top is approached.
!
         pblhgt = max(pblht(i),1.0_r8)
         if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0 ) then
            fac1   = max(0.0_r8,1.0-zm(i,k+1)/pblhgt)
            tprime = min(tpert(i),tpmax)*fac1
            qsattp = shbs(i,k+1) + cp*rhlat*gam(i,k+1)*tprime
            shprme = min(min(qpert(i,1),shpmax)*fac1,max(qsattp-shb(i,k+1),0.0_r8))
            qprime = max(qprime,shprme)
         else
            tprime = 0.0
            qprime = 0.0
         end if
!
! Specify "updraft" (in-cloud) thermodynamic properties
!
         sc (i)    = sb (i,k+1) + cp*tprime
         shc(i)    = shb(i,k+1) + qprime
         hc (i)    = sc (i    ) + hlat*shc(i)
         vtemp4(i) = hc(i) - hbs(i,k)
         dz        = pdel(i,k)*rgas*tb(i,k)*rgrav/pmid(i,k)
         if (vtemp4(i) > 0.0) then
            dzcld(i) = dzcld(i) + dz
         else
            dzcld(i) = 0.0
         end if
10       continue
#if ( defined DIAGNS )
!#######################################################################
!#                                                                     #
!#    output thermodynamic perturbation information                    #
!#                                                                     #
         do i=1,ncol
           if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
            write (6,8090) k+1,sc(iloc),shc(iloc),hc(iloc)
           end if
         enddo
!#                                                                     #
!#######################################################################
#endif
!
! Check on moist convective instability
! Build index vector of points where instability exists
!
         len1 = 0
         do i=1,ncol
            if (vtemp4(i) > 0.0) then
               len1 = len1 + 1
               indx1(len1) = i
            end if
         end do
         if (len1 <= 0) go to 70
!
! Current level just below top level => no overshoot
!
         if (k <= limcnv+1) then
            do ii=1,len1
               i = indx1(ii)
               temp1     = vtemp4(i)/(1.0 + gam(i,k))
               cldwtr(i) = max(0.0_r8,(sb(i,k) - sc(i) + temp1))
               beta(i)   = 0.0
               vtemp3(i) = (1.0 + gam(i,k))*(sc(i) - sbh(i,k))
            end do
         else
!
! First guess at overshoot parameter using crude buoyancy closure
! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
! If pre-existing supersaturation in detrainment layer, beta=0
! cldwtr is temporarily equal to hlat*l (l=> liquid water)
!
!CDIR$ IVDEP
            do ii=1,len1
               i = indx1(ii)
               temp1     = vtemp4(i)/(1.0 + gam(i,k))
               cldwtr(i) = max(0.0_r8,(sb(i,k)-sc(i)+temp1))
               betamx(i) = 1.0 - c0*max(0.0_r8,(dzcld(i)-dzmin))
               b1        = (hc(i) - hbs(i,k-1))*pdel(i,k-1)
               b2        = (hc(i) - hbs(i,k  ))*pdel(i,k  )
               beta(i)   = max(betamn,min(betamx(i), 1.0 + b1/b2))
               if (hbs(i,k-1) <= hb(i,k-1)) beta(i) = 0.0
!
! Bound maximum beta to ensure physically realistic solutions
!
! First check constrains beta so that eta remains positive
! (assuming that eta is already positive for beta equal zero)
!
               vtemp1(i) = -(hbh(i,k+1) - hc(i))*pdel(i,k)*rpdel(i,k+1)+ &
                           (1.0 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i))
               vtemp2(i) = (1.0 + gam(i,k))*(sc(i) - sbh(i,k))
               vtemp3(i) = vtemp2(i)
               if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0.) then
                  betamx(i) = 0.99*(vtemp1(i)/vtemp2(i))
                  beta(i)   = max(0.0_r8,min(betamx(i),beta(i)))
               end if
            end do
!
! Second check involves supersaturation of "detrainment layer"
! small amount of supersaturation acceptable (by ssfac factor)
!
!CDIR$ IVDEP
            do ii=1,len1
               i = indx1(ii)
               if (hb(i,k-1) < hbs(i,k-1)) then
                  vtemp1(i) = vtemp1(i)*rpdel(i,k)
                  temp2 = gam(i,k-1)*(sbh(i,k) - sc(i) + cldwtr(i)) -  &
                          hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
                  temp3 = vtemp3(i)*rpdel(i,k)
                  vtemp2(i) = (ztodt/cats)*(hc(i) - hbs(i,k))*temp2/ &
                              (pdel(i,k-1)*(hbs(i,k-1) - hb(i,k-1))) + temp3
                  if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0.) then
                     betamx(i) = ssfac*(vtemp1(i)/vtemp2(i))
                     beta(i)   = max(0.0_r8,min(betamx(i),beta(i)))
                  end if
               else
                  beta(i) = 0.0
               end if
            end do
!
! Third check to avoid introducing 2 delta x thermodynamic
! noise in the vertical ... constrain adjusted h (or theta e)
! so that the adjustment doesn't contribute to "kinks" in h
!
!CDIR$ IVDEP
            do ii=1,len1
               i = indx1(ii)
               g = min(0.0_r8,hb(i,k) - hb(i,k-1))
               temp1 = (hb(i,k) - hb(i,k-1) - g)*(cats/ztodt)/(hc(i) - hbs(i,k))
               vtemp1(i) = temp1*vtemp1(i) + (hc(i) - hbh(i,k+1))*rpdel(i,k)
               vtemp2(i) = temp1*vtemp3(i)*rpdel(i,k) + (hc(i) - hbh(i,k) - cldwtr(i))* &
                           (rpdel(i,k) + rpdel(i,k+1))
               if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0.) then
                  if (vtemp2(i) /= 0.0) then
                     betamx(i) = vtemp1(i)/vtemp2(i)
                  else
                     betamx(i) = 0.0
                  end if
                  beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
               end if
            end do
         end if
!
! Calculate mass flux required for stabilization.
!
! Ensure that the convective mass flux, eta, is positive by
! setting negative values of eta to zero..
! Ensure that estimated mass flux cannot move more than the
! minimum of total mass contained in either layer k or layer k+1.
! Also test for other pathological cases that result in non-
! physical states and adjust eta accordingly.
!
!CDIR$ IVDEP
         do ii=1,len1
            i = indx1(ii)
            beta(i) = max(0.0_r8,beta(i))
            temp1 = hc(i) - hbs(i,k)
            temp2 = ((1.0 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i)) - &
                      beta(i)*vtemp3(i))*rpdel(i,k) - (hbh(i,k+1) - hc(i))*rpdel(i,k+1)
            eta(i) = temp1/(temp2*grav*cats)
            tmass = min(pdel(i,k),pdel(i,k+1))*rgrav
            if (eta(i) > tmass*rtdt .or. eta(i) <= 0.0) eta(i) = 0.0
!
! Check on negative q in top layer (bound beta)
!
            if (shc(i)-shbh(i,k) < 0.0 .and. beta(i)*eta(i) /= 0.0) then
               denom = eta(i)*grav*ztodt*(shc(i) - shbh(i,k))*rpdel(i,k-1)
               beta(i) = max(0.0_r8,min(-0.999*shb(i,k-1)/denom,beta(i)))
            end if
!
! Check on negative q in middle layer (zero eta)
!
            qtest1 = shb(i,k) + eta(i)*grav*ztodt*((shc(i) - shbh(i,k+1)) - &
                     (1.0 - beta(i))*cldwtr(i)*rhlat - beta(i)*(shc(i) - shbh(i,k)))* &
	             rpdel(i,k)
            if (qtest1 <= 0.0) eta(i) = 0.0
!
! Check on negative q in lower layer (bound eta)
!
            fac1 = -(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
            qtest2 = shb(i,k+1) - eta(i)*grav*ztodt*fac1
            if (qtest2 < 0.0) then
               eta(i) = 0.99*shb(i,k+1)/(grav*ztodt*fac1)
            end if
            etagdt(i) = eta(i)*grav*ztodt
         end do
!
#if ( defined DIAGNS )
!#######################################################################
!#                                                                     #
         do i=1,ncol
           if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep >= nsloc)) then
            write(6,8080) beta(iloc), eta(iloc)
           end if
         enddo
!#                                                                     #
!#######################################################################
#endif
!
! Calculate cloud water, rain water, and thermodynamic changes
!
!CDIR$ IVDEP
         do 30 ii=1,len1
            i = indx1(ii)
            icwmr(i,k) = cldwtr(i)*rhlat
            cldwtr(i) = etagdt(i)*cldwtr(i)*rhlat*rgrav
            rnwtr(i) = (1.0 - beta(i))*cldwtr(i)