#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#endifcontainssubroutine 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 = 120 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 '   returnend subroutine inimcsubroutine 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