cldwat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#undef DEBUG
#include <misc.h>
#include <params.h>
module cldwat
!----------------------------------------------------------------------- 
! 
! Purpose: Prognostic cloud water data and methods.
! 
! Public interfaces:
!
! inimc -- Initialize constants
! pcond -- Calculate prognostic condensate
!
! Author: P. Rasch, with Modifications by Minghua Zhang
! 
!-----------------------------------------------------------------------
   use precision,     only: r8
   use pmgrid,        only: masterproc
   use ppgrid,        only: pcols, pver, pverp
   use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, &
                            cp, epsqs, ttrice

   implicit none

!-----------------------------------------------------------------------
! PUBLIC: Make default data and interfaces private
!-----------------------------------------------------------------------
   private
   public inimc, pcond          ! Public interfaces
   integer, public::  ktop      ! Level above 10 hPa
!-----------------------------------------------------------------------
! PRIVATE: Everything else is private to this module
!-----------------------------------------------------------------------
   real(r8), private:: rhonot   ! air density at surface
   real(r8), private:: t0       ! Freezing temperature
   real(r8), private:: cldmin   ! assumed minimum cloud amount
   real(r8), private:: small    ! small number compared to unity
   real(r8), private:: c        ! constant for graupel like snow cm**(1-d)/s
   real(r8), private:: d        ! constant for graupel like snow
   real(r8), private:: esi      ! collection efficient for ice by snow
   real(r8), private:: esw      ! collection efficient for water by snow
   real(r8), private:: nos      ! particles snow / cm**4
   real(r8), private:: pi       ! Mathematical constant
   real(r8), private:: gravit   ! Gravitational acceleration at surface
   real(r8), private:: rh2o
   real(r8), private:: prhonos
   real(r8), private:: thrpd    ! numerical three added to d
   real(r8), private:: gam3pd   ! gamma function on (3+d)
   real(r8), private:: gam4pd   ! gamma function on (4+d)
   real(r8), private:: rhoi     ! ice density
   real(r8), private:: rhos     ! snow density
   real(r8), private:: rhow     ! water density
   real(r8), private:: mcon01   ! constants used in cloud microphysics
   real(r8), private:: mcon02   ! constants used in cloud microphysics
   real(r8), private:: mcon03   ! constants used in cloud microphysics
   real(r8), private:: mcon04   ! constants used in cloud microphysics
   real(r8), private:: mcon05   ! constants used in cloud microphysics
   real(r8), private:: mcon06   ! constants used in cloud microphysics
   real(r8), private:: mcon07   ! constants used in cloud microphysics
   real(r8), private:: mcon08   ! constants used in cloud microphysics

   integer, private ::  k1mb    ! index of the eta level near 1 mb
#ifdef DEBUG
   integer, private,parameter ::  nlook = 2  ! Number of points to examine
   integer, private ::  ilook(nlook)         ! Longitude index to examine
   integer, private ::  latlook(nlook)       ! Latitude index to examine
   integer, private ::  lchnklook(nlook)     ! Chunk index to examine
   integer, private ::  icollook(nlook)      ! Column index to examine
#endif

contains

subroutine inimc( tmeltx, rhonotx, gravitx, rh2ox )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! initialize constants for the prognostic condensate
! 
! Author: P. Rasch, April 1997
! 
!-----------------------------------------------------------------------
   use phys_grid, only: get_chunk_coord_p
   use pmgrid, only: plev, plevp
   integer k
   real(r8), intent(in) :: tmeltx
   real(r8), intent(in) :: rhonotx
   real(r8), intent(in) :: gravitx
   real(r8), intent(in) :: rh2ox

#include <comhyb.h>

#ifdef CRAY
   real(r8) signgam              ! variable required by cray gamma function
   external gamma
#endif
   rhonot = rhonotx          ! air density at surface (gm/cm3)
   gravit = gravitx
   rh2o   = rh2ox
   rhos = .1                 ! assumed snow density (gm/cm3)
   rhow = 1.                 ! water density
   rhoi = 1.                 ! ice density
   esi = 1.0                 ! collection efficient for ice by snow
   esw = 0.1                 ! collection efficient for water by snow
   t0 = tmeltx               ! approximate freezing temp
   cldmin = 0.02             ! assumed minimum cloud amount
   small = 1.e-22            ! a small number compared to unity
   c = 152.93                ! constant for graupel like snow cm**(1-d)/s
   d = 0.25                  ! constant for graupel like snow
   nos = 3.e-2               ! particles snow / cm**4
   pi = 4.*atan(1.0)
   prhonos = pi*rhos*nos
   thrpd = 3. + d
#ifdef CRAY
   call gamma(3.+d, signgam, gam3pd)
   gam3pd = sign(exp(gam3pd),signgam)
   call gamma(4.+d, signgam, gam4pd)
   gam4pd = sign(exp(gam4pd),signgam)
   write (6,*) ' d, gamma(3+d), gamma(4+d) =', gam3pd, gam4pd
#else
   if (d==0.25) then
      gam3pd = 2.549256966718531 ! only right for d = 0.25
      gam4pd = 8.285085141835282
   else
      write (6,*) ' can only use d ne 0.25 on a cray '
      stop
   endif
#endif
   mcon01 = pi*nos*c*gam3pd/4.
   mcon02 = 1./(c*gam4pd*sqrt(rhonot)/(6*prhonos**(d/4.)))
   mcon03 = -(0.5+d/4.)
   mcon04 = 4./(4.+d)
   mcon05 = (3+d)/(4+d)
   mcon06 = (3+d)/4.
   mcon07 = mcon01*sqrt(rhonot)*mcon02**mcon05/prhonos**mcon06
   mcon08 = -0.5/(4.+d)

!  find the level about 1mb, we wont do the microphysics above this level
   k1mb = 1
   do k=1,pver-1
      if (hypm(k) < 1.e2 .and. hypm(k+1) >= 1.e2) then
         if (1.e2-hypm(k) < hypm(k+1)-1.e2) then
            k1mb = k
         else
            k1mb = k + 1
         end if
         goto 20
      end if
   end do
   if (masterproc) then
      write(6,*)'inimc: model levels bracketing 1 mb not found'
   end if
!  call endrun
   k1mb = 1
20 if( masterproc ) write(6,*)'inimc: model level nearest 1 mb is',k1mb,'which is',hypm(k1mb),'pascals'

#ifdef DEBUG
!
! Set indicies of the point to examine for debugging
!
   latlook(:) = (/64, 32/)   ! Latitude indices to examine
   ilook(:)   = (/1,   1/)   ! Longitude indicex to examine
   call get_chunk_coord_p( nlook, ilook, latlook, icollook, lchnklook )
#endif
   if( masterproc ) write (6,*) 'cloud water initialization by inimc complete '

   return
end subroutine inimc

subroutine pcond (lchnk   ,ncol    , &
                  tn      ,ttend   ,qn      ,qtend   ,omega   , &
                  cwat    ,p       ,pdel    ,cldn    , &
                  cme     ,evapr   ,prain   ,rmelt   , &     
                  deltat  ,pcflx   ,fwaut   ,fsaut   ,fracw   , &
                  fsacw   ,fsaci   ,lctend  ,rhdfda  ,rhu00, icefrac)  
!----------------------------------------------------------------------- 
! 
! Purpose: 
! The public interface to the cloud water parameterization
! returns tendencies to water vapor, temperature and cloud water variables
! 
! For basic 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.
! 
! For important modifications to improve the method of determining
! condensation/evaporation see Zhang et al (2001, in preparation)
!
! Authors: M. Zhang, W. Lin, P. Rasch and J.E. Kristjansson
!-----------------------------------------------------------------------
   use wv_saturation, only: vqsatd
!
!---------------------------------------------------------------------
!
! Input Arguments
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns

   real(r8), intent(in) :: cldn(pcols,pver)     ! new value of cloud fraction    (fraction)
   real(r8), intent(in) :: cwat(pcols,pver)     ! cloud water (kg/kg)
   real(r8), intent(in) :: omega(pcols,pver)    ! vert pressure vel (Pa/s)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure          (K)
   real(r8), intent(in) :: pcflx(pcols,pverp)   ! convective precip level by level (kg/m2/s)  (DISABLED)
   real(r8), intent(in) :: pdel(pcols,pver)     ! pressure thickness (Pa)
   real(r8), intent(in) :: qn(pcols,pver)       ! new water vapor    (kg/kg)
   real(r8), intent(in) :: qtend(pcols,pver)    ! mixing ratio tend  (kg/kg/s)
   real(r8), intent(in) :: tn(pcols,pver)       ! new temperature    (K)
   real(r8), intent(in) :: ttend(pcols,pver)    ! temp tendencies    (K/s)
   real(r8), intent(in) :: deltat               ! time step to advance solution over
   real(r8), intent(in) :: lctend(pcols,pver)   ! cloud liquid water tendencies   ====wlin
   real(r8), intent(in) :: rhdfda(pcols,pver)   ! dG(a)/da, rh=G(a), when rh>u00  ====wlin
   real(r8), intent(in) :: rhu00 (pcols,pver)   ! Rhlim for cloud                 ====wlin
   real(r8), intent(in) :: icefrac(pcols)       ! sea ice fraction  (fraction)
!
! Output Arguments
!
   real(r8), intent(out) :: cme(pcols,pver)      ! rate of cond-evap within the cloud
   real(r8), intent(out) :: evapr(pcols,pver)    ! rate of evaporation of falling precip (1/s)
   real(r8), intent(out) :: prain(pcols,pver)    ! rate of conversion of condensate to precip (1/s)
   real(r8), intent(out) :: rmelt(pcols,pver)    ! heating rate due to precip phase change (K/s) (DISABLED)

!
! Local workspace
!
   integer i                 ! work variable
   integer iter              ! #iterations for precipitation calculation
   integer k                 ! work variable
   integer l                 ! work variable

   real(r8) cldm(pcols)          ! mean cloud fraction over the time step
   real(r8) cldmax(pcols)        ! max cloud fraction above
   real(r8) coef(pcols)          ! conversion time scale for condensate to rain
   real(r8) conke                ! rate of evaporation of precipitation:
   real(r8) cwm(pcols)           ! cwat mixing ratio at midpoint of time step
   real(r8) cwn(pcols)           ! cwat mixing ratio at end
   real(r8) denom                ! work variable
   real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
   real(r8) es(pcols)            ! sat. vapor pressure
   real(r8) fice(pcols)          ! fraction of cwat that is ice
   real(r8) fracw(pcols,pver)    ! relative importance of collection of liquid by rain
   real(r8) fsaci(pcols,pver)    ! relative importance of collection of ice by snow
   real(r8) fsacw(pcols,pver)    ! relative importance of collection of liquid by snow
   real(r8) fsaut(pcols,pver)    ! relative importance of ice auto conversion
   real(r8) fwaut(pcols,pver)    ! relative importance of warm cloud autoconversion
   real(r8) gamma(pcols)         ! d qs / dT
   real(r8) iceab(pcols)         ! rate of ice only from above
   real(r8) icwc(pcols)          ! in-cloud water content (kg/kg)
   real(r8) mincld               ! a small cloud fraction to avoid / zero
   real(r8) omeps                ! 1 minus epsilon
   real(r8) omsm                 ! a number just less than unity (for rounding)
   real(r8) precab(pcols)        ! rate of precipitation (kg / (m**2 * s))
   real(r8) prect(pcols)         ! rate of precipitation including convection (kg / (m**2 * s))
   real(r8) prprov(pcols)        ! provisional value of precip at btm of layer
   real(r8) prtmp                ! work variable
   real(r8) q(pcols,pver)        ! mixing ratio before time step ignoring condensate
   real(r8) qs(pcols)            ! spec. hum. of water vapor
   real(r8) qsn, esn             ! work variable
   real(r8) qsp(pcols,pver)      ! sat pt mixing ratio
   real(r8) qtl(pcols)           ! tendency which would saturate the grid box in deltat
   real(r8) qtmp, ttmp           ! work variable
   real(r8) relhum1(pcols)        ! relative humidity
   real(r8) relhum(pcols)        ! relative humidity
   real(r8) tc                   ! crit temp of transition to ice
   real(r8) t(pcols,pver)        ! temp before time step ignoring condensate
   real(r8) tsp(pcols,pver)      ! sat pt temperature
   real(r8) pol                  ! work variable
   real(r8) cdt                  ! work variable

! Extra local work space for cloud scheme modification       

   real(r8) cpohl                !Cp/Hlatv
   real(r8) hlocp                !Hlatv/Cp
   real(r8) dto2                 !0.5*deltat (delta=2.0*dt)
   real(r8) calpha(pcols)        !alpha of new C - E scheme formulation
   real(r8) cbeta (pcols)        !beta  of new C - E scheme formulation
   real(r8) cbetah(pcols)        !beta_hat at saturation portion 
   real(r8) cgamma(pcols)        !gamma of new C - E scheme formulation
   real(r8) cgamah(pcols)        !gamma_hat at saturation portion
   real(r8) rcgama(pcols)        !gamma/gamma_hat
   real(r8) csigma(pcols)        !sigma of new C - E scheme formulation
   real(r8) cmec1 (pcols)        !c1    of new C - E scheme formulation
   real(r8) cmec2 (pcols)        !c2    of new C - E scheme formulation
   real(r8) cmec3 (pcols)        !c3    of new C - E scheme formulation
   real(r8) cmec4 (pcols)        !c4    of new C - E scheme formulation
   real(r8) cmeres(pcols)        !residual cond of over-sat after cme and evapr
   real(r8) ctmp                 !a scalar representation of cmeres
   real(r8) clrh2o               ! Ratio of latvap to water vapor gas const
!
!------------------------------------------------------------
#include <comadj.h>              
!------------------------------------------------------------
!
   clrh2o = hlatv/rh2o   ! Ratio of latvap to water vapor gas const
   omeps = 1.0 - epsqs
#ifdef PERGRO
   mincld = 1.e-4
   iter = 1   ! number of times to iterate the precipitation calculation
#else
   mincld = 1.e-4
   iter = 2
#endif
   omsm = 0.99999
   cpohl = cp/hlatv
   hlocp = hlatv/cp
   dto2=0.5*deltat
!
! Constant for computing rate of evaporation of precipitation:
!
   conke = 1.e-5
!
! initialize a few single level fields
!
   do i = 1,ncol
      precab(i) = 0.0
      prect(i) = 0.0
      iceab(i) = 0.0                ! latent heat of precip above
      cldmax(i) = 0.0
   end do
!
! initialize multi-level fields 
!
   do k = 1,pver
      do i = 1,ncol
         q(i,k) = qn(i,k) 
         t(i,k) = tn(i,k)
      end do
   end do
   cme  (:ncol,:) = 0._r8
   evapr(:ncol,:) = 0._r8
   prain(:ncol,:) = 0._r8
   rmelt(:ncol,:) = 0._r8
   fwaut(:ncol,:) = 0._r8