cldwat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

   fsaut(:ncol,:) = 0._r8
   fracw(:ncol,:) = 0._r8
   fsacw(:ncol,:) = 0._r8
   fsaci(:ncol,:) = 0._r8
!
! find the wet bulb temp and saturation value
! for the provisional t and q without condensation
!
   call findsp (lchnk, ncol, qn, tn, p, tsp, qsp)
   do 800 k = k1mb,pver
      call vqsatd (t(1,k), p(1,k), es, qs, gamma, ncol)
      do i = 1,ncol
         relhum(i) = q(i,k)/qs(i)
!
         cldm(i) = max(cldn(i,k),mincld)
!
! the max cloud fraction above this level
!
         cldmax(i) = max(cldmax(i), cldm(i))

! define the coefficients for C - E calculation

         calpha(i) = 1.0/qs(i)
         cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl
         cbetah(i) = 1.0/qs(i)*gamma(i)*cpohl
         cgamma(i) = calpha(i)+hlatv*cbeta(i)/cp
         cgamah(i) = calpha(i)+hlatv*cbetah(i)/cp
         rcgama(i) = cgamma(i)/cgamah(i)

         if(cldm(i) > mincld) then
            icwc(i) = max(0._r8,cwat(i,k)/cldm(i))
         else
            icwc(i) = 0.0
         endif

!
! initial guess of evaporation, will be updated within iteration
!
         evapr(i,k) = conke*(1. - cldm(i))*sqrt(precab(i)) &
                        *(1. - min(relhum(i),1._r8))

!
! zero cmeres before iteration for each level
!
         cmeres(i)=0.0

      end do
      do i = 1,ncol
!
! fractions of ice at this level
!
         tc = t(i,k) - t0
         fice(i) = max(0._r8,min(-tc*0.05,1.0_r8))
!
! calculate the cooling due to a phase change of the rainwater
! from above
!
         if (t(i,k) >= t0) then
!              rmelt(i,k) =  -hlatf/cp*iceab(i)*gravit/pdel(i,k)
            rmelt(i,k) = 0.
            iceab(i) = 0.
         else
            rmelt(i,k) = 0.
         endif
      end do

!
! calculate cme and formation of precip. 
!
! The cloud microphysics is highly nonlinear and coupled with cme
! Both rain processes and cme are calculated iteratively.
! 
      do 100 l = 1,iter

        do i = 1,ncol

!
! calculation of cme has 4 scenarios
! ==================================
!
           if(relhum(i) > rhu00(i,k)) then
    
           ! 1. whole grid saturation
           ! ========================
              if(relhum(i) >= 0.999_r8 .or. cldm(i) >= 0.999_r8 ) then
                 cme(i,k)=(calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i)

           ! 2. fractional saturation
           ! ========================
              else
                  csigma(i) = 1.0/(rhdfda(i,k)+cgamma(i)*icwc(i))
                  cmec1(i) = (1.0-cldm(i))*csigma(i)*rhdfda(i,k)
                  cmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0-rcgama(i)*cldm(i))*   &
                             csigma(i)*calpha(i)*icwc(i)
                  cmec3(i) = cldm(i)*cbetah(i)/cgamah(i) +  &
                           (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i)
                  cmec4(i) = csigma(i)*cgamma(i)*icwc(i)

                  ! Q=C-E=-C1*Al + C2*Aq - C3* At + C4*Er
  
                  cme(i,k) = -cmec1(i)*lctend(i,k) + cmec2(i)*qtend(i,k)  &
                             -cmec3(i)*ttend(i,k) + cmec4(i)*evapr(i,k)
               endif

           ! 3. when rh < rhu00, evaporate existing cloud water
           ! ================================================== 
           else if(cwat(i,k) > 0.0)then
              ! liquid water should be evaporated but not to exceed 
              ! saturation point. if qn > qsp, not to evaporate cwat
              cme(i,k)=-min(max(0._r8,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat 

           ! 4. no condensation nor evaporation
           ! ==================================
           else
              cme(i,k)=0.0
           endif

  
        end do    !end loop for cme update

! Because of the finite time step, 
! place a bound here not to exceed wet bulb point
! and not to evaporate more than available water
!
         do i = 1, ncol
            qtmp = qn(i,k) - cme(i,k)*deltat

! possibilities to have qtmp > qsp
!
!   1. if qn > qs(tn), it condenses; 
!      if after applying cme,  qtmp > qsp,  more condensation is applied. 
!      
!   2. if qn < qs, evaporation should not exceed qsp,
    
            if(qtmp > qsp(i,k)) then
              cme(i,k) = cme(i,k) + (qtmp-qsp(i,k))/deltat
            endif

!
! if net evaporation, it should not exceed available cwat
!
            if(cme(i,k) < -cwat(i,k)/deltat)  &
               cme(i,k) = -cwat(i,k)/deltat
!
! addition of residual condensation from previous step of iteration
!
            cme(i,k) = cme(i,k) + cmeres(i)

         end do

         do i = 1,ncol
!
! as a safe limit, condensation should not reduce grid mean rh below rhu00
! 
           if(cme(i,k) > 0.0 .and. relhum(i) > rhu00(i,k) )  &
              cme(i,k) = min(cme(i,k), (qn(i,k)-qs(i)*rhu00(i,k))/deltat)
!
! initial guess for cwm (mean cloud water over time step) if 1st iteration
!
           if(l < 2) then
             cwm(i) = max(cwat(i,k)+cme(i,k)*dto2,  0._r8)
           endif

         enddo

! provisional precipitation falling through model layer
         do i = 1,ncol
            prprov(i) = prect(i) + prain(i,k)*pdel(i,k)/gravit
         end do

! calculate conversion of condensate to precipitation by cloud microphysics 
         call findmcnew (lchnk   ,ncol    , &
                         k       ,prprov  ,t       ,p       , &
                         cwm     ,cldm    ,cldmax  ,fice    ,coef    , &
                         fwaut(1,k),fsaut(1,k),fracw(1,k),fsacw(1,k),fsaci(1,k), &
                         icefrac)
!
! calculate the precip rate
!
         do i = 1,ncol
            if (cldm(i) > 0) then  
!
! first predict the cloud water
!
               cdt = coef(i)*deltat
               if (cdt > 0.01) then
                  pol = cme(i,k)/coef(i) ! production over loss
                  cwn(i) = max(0._r8,(cwat(i,k)-pol)*exp(-cdt)+ pol)
               else
                  cwn(i) = max(0._r8,(cwat(i,k) + cme(i,k)*deltat)/(1+cdt))
               endif
!
! now back out the tendency of net rain production
!
               prain(i,k) =  max(0._r8,cme(i,k)-(cwn(i)-cwat(i,k))/deltat)
            else
               prain(i,k) = 0.0
               cwn(i) = 0.
            endif
!
! update any remaining  provisional values
!
            cwm(i) = (cwn(i) + cwat(i,k))*0.5
!
! update in cloud water
!
            if(cldm(i) > mincld) then
               icwc(i) = cwm(i)/cldm(i)
            else
               icwc(i) = 0.0
            endif

         end do              ! end of do i = 1,ncol

!
! calculate provisional value of cloud water for
! evaporation of rain (evapr) calculation
!
      do i = 1,ncol
         qtmp = qn(i,k) - cme(i,k)*deltat
         ttmp = tn(i,k) + rmelt(i,k)*deltat + hlocp*deltat*cme(i,k)
         esn = estblf(ttmp)
         qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)
         qtl(i) = max((qsn - qtmp)/deltat,0._r8)
         relhum1(i) = qtmp/qsn
      end do
!
      do i = 1,ncol
#ifdef PERGRO
         evapr(i,k) = conke*(1. - max(cldm(i),mincld))* &
                      sqrt(precab(i))*(1. - min(relhum1(i),1._r8))
#else
         evapr(i,k) = conke*(1. - cldm(i))*sqrt(precab(i)) &
                      *(1. - min(relhum1(i),1._r8))
#endif
!
! limit the evaporation to the amount which is entering the box
! or saturates the box
!
         prtmp = precab(i)*gravit/pdel(i,k)
         evapr(i,k) = min(evapr(i,k), prtmp, qtl(i))*omsm
#ifdef PERGRO
!           zeroing needed for pert growth
         evapr(i,k) = 0.
#endif
      end do

! now remove the residual of any over-saturation. Normally,
! the oversaturated water vapor should have been removed by 
! cme formulation plus constraints by wet bulb tsp/qsp
! as computed above. However, because of non-linearity,
! addition of (cme-evapr) to update t and q may still cause
! a very small amount of over saturation. It is called a
! residual of over-saturation because theoretically, cme
! should have taken care of all of large scale condensation.
! 

       do i = 1,ncol
          qtmp = qn(i,k)-(cme(i,k)-evapr(i,k))*deltat
          ttmp = tn(i,k)+(rmelt(i,k)+hlocp*(cme(i,k)-evapr(i,k)) )  &
                      *deltat
          esn = estblf(ttmp)
          qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)
          !
          !Upper stratosphere and mesosphere, qsn calculated
          !above may be negative. Here just to skip it instead
          !of resetting it to 1 as in aqsat
          !
          if(qtmp > qsn .and. qsn > 0) then
             !calculate dqsdt, a more precise calculation
             !which taking into account different range of T 
             !can be found in aqsatd.F. Here follows
             !cond.F to calculate it.
             !
             denom = (p(i,k)-omeps*esn)*ttmp*ttmp
             dqsdt = clrh2o*qsn*p(i,k)/denom
             !
             !now extra condensation to bring air to just saturation
             !
             ctmp = (qtmp-qsn)/(1.+hlocp*dqsdt)/deltat
             cme(i,k) = cme(i,k)+ctmp
!
! save residual on cmeres to addtion to cme on entering next iteration
! cme exit here contain the residual but overrided if back to iteration
!
             cmeres(i) = ctmp
          else
             cmeres(i) = 0.0
          endif
       end do
              
 100 continue              ! end of do l = 1,iter
!
! precipitation
!
      do i = 1,ncol
         prtmp = pdel(i,k) / gravit *(prain(i,k) - evapr(i,k))
         iceab(i) = iceab(i) + fice(i)*prtmp
         precab(i) = precab(i) + prtmp
         prect(i) = prect(i) + prtmp + pcflx(i,k+1)
         if ((precab(i)) < 1.e-10) then      
            precab(i) = 0.
         endif
         if ((prect(i)) < 1.e-10) then      
            prect(i) = 0.
         endif
      end do
 800 continue                ! level loop (k=1,pver)

   return
end subroutine pcond

!##############################################################################

subroutine findmcnew (lchnk   ,ncol    , &
                      k       ,precab  ,t       ,p       , &
                      cwm     ,cldm    ,cldmax  ,fice    ,coef    , &
                      fwaut   ,fsaut   ,fracw   ,fsacw   ,fsaci   , &
                      icefrac )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! calculate the conversion of condensate to precipitate
! 
! Method: 
! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
!  model climate using diagnosed and 
!  predicted condensate parameterizations, 1998, J. Clim., 11,
!  pp1587---1614.
! 
! Author: P. Rasch
! 
!-----------------------------------------------------------------------
   use phys_grid, only: get_rlat_all_p