cldwat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

   use comsrf,        only: landm
!
! input args
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns
   integer, intent(in) :: k                     ! level index

   real(r8), intent(in) :: precab(pcols)        ! rate of precipitation from above (kg / (m**2 * s))
   real(r8), intent(in) :: t(pcols,pver)        ! temperature       (K)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure          (Pa)
   real(r8), intent(in) :: cldm(pcols)          ! cloud fraction
   real(r8), intent(in) :: cldmax(pcols)        ! max cloud fraction above this level
   real(r8), intent(in) :: cwm(pcols)           ! condensate mixing ratio (kg/kg)
   real(r8), intent(in) :: fice(pcols)          ! fraction of cwat that is ice
   real(r8), intent(in) :: icefrac(pcols)       ! sea ice fraction 

! output arguments
   real(r8), intent(out) :: coef(pcols)          ! conversion rate (1/s)
   real(r8), intent(out) :: fwaut(pcols)         ! relative importance of liquid autoconversion (a diagnostic)
   real(r8), intent(out) :: fsaut(pcols)         ! relative importance of ice autoconversion    (a diagnostic)
   real(r8), intent(out) :: fracw(pcols)         ! relative importance of rain accreting liquid (a diagnostic)
   real(r8), intent(out) :: fsacw(pcols)         ! relative importance of snow accreting liquid (a diagnostic)
   real(r8), intent(out) :: fsaci(pcols)         ! relative importance of snow accreting ice    (a diagnostic)

! work variables

   integer i
   integer ii
   integer ind(pcols)
   integer ncols

   real(r8), parameter :: degrad = 57.296 ! divide by this to convert degrees to radians
   real(r8) alpha                ! ratio of 3rd moment radius to 2nd
   real(r8) capc                 ! constant for autoconversion
   real(r8) capn                 ! local cloud particles / cm3
   real(r8) capnoice             ! local cloud particles when not over sea ice / cm3
   real(r8) capnsi               ! sea ice cloud particles / cm3
   real(r8) capnc                ! cold and oceanic cloud particles / cm3
   real(r8) capnw                ! warm continental cloud particles / cm3
   real(r8) ciaut                ! coefficient of autoconversion of ice (1/s)
   real(r8) ciautb               ! coefficient of autoconversion of ice (1/s)
   real(r8) cldloc(pcols)        ! non-zero amount of cloud
   real(r8) cldpr(pcols)         ! assumed cloudy volume occupied by rain and cloud
   real(r8) con1                 ! work constant
   real(r8) con2                 ! work constant
   real(r8) convfw               ! constant used for fall velocity calculation
   real(r8) cracw                ! constant used for rain accreting water
   real(r8) critpr               ! critical precip rate collection efficiency changes
   real(r8) csacx                ! constant used for snow accreting liquid or ice
   real(r8) dtice                ! interval for transition from liquid to ice
   real(r8) effc                 ! collection efficiency
   real(r8) icemr(pcols)         ! in-cloud ice mixing ratio
   real(r8) icrit                ! threshold for autoconversion of ice
   real(r8) icritc               ! threshold for autoconversion of cold ice
   real(r8) icritw               ! threshold for autoconversion of warm ice
   real(r8) kconst               ! const for terminal velocity (stokes regime)
   real(r8) liqmr(pcols)         ! in-cloud liquid water mixing ratio
   real(r8) pracw                ! rate of rain accreting water
   real(r8) prlloc(pcols)        ! local rain flux in mm/day
   real(r8) prscgs(pcols)        ! local snow amount in cgs units
   real(r8) psaci                ! rate of collection of ice by snow (lin et al 1983)
   real(r8) psacw                ! rate of collection of liquid by snow (lin et al 1983)
   real(r8) psaut                ! rate of autoconversion of ice condensate
   real(r8) ptot                 ! total rate of conversion
   real(r8) pwaut                ! rate of autoconversion of liquid condensate
   real(r8) r3l                  ! volume radius
   real(r8) r3lcrit              ! critical radius at which autoconversion become efficient
   real(r8) rainmr(pcols)        ! in-cloud rain mixing ratio
   real(r8) rat1                 ! work constant
   real(r8) rat2                 ! work constant
   real(r8) rdtice               ! recipricol of dtice
   real(r8) rho(pcols)           ! density (mks units)
   real(r8) rhocgs               ! density (cgs units)
   real(r8) rlat(pcols)          ! latitude (radians)
   real(r8) snowfr               ! fraction of precipate existing as snow
   real(r8) totmr(pcols)         ! in-cloud total condensate mixing ratio
   real(r8) vfallw               ! fall speed of precipitate as liquid
   real(r8) wp                   ! weight factor used in calculating pressure dep of autoconversion
   real(r8) wsi                  ! weight factor for sea ice
   real(r8) wt                   ! fraction of ice
   real(r8) wland                ! fraction of land

!      real(r8) csaci
!      real(r8) csacw
!      real(r8) cwaut
!      real(r8) efact
!      real(r8) lamdas
!      real(r8) lcrit
!      real(r8) rcwm
!      real(r8) r3lc2
!      real(r8) snowmr(pcols)
!      real(r8) vfalls


!     inline statement functions
   real(r8) heavy, heavym, a1, a2, heavyp, heavymp
   heavy(a1,a2) = max(0._r8,sign(1._r8,a1-a2))  ! heavyside function
   heavym(a1,a2) = max(0.01_r8,sign(1._r8,a1-a2))  ! modified heavyside function
!
! New heavyside functions to perhaps address error growth problems
!
   heavyp(a1,a2) = a1/(a2+a1+1.e-36)
   heavymp(a1,a2) = (a1+0.01*a2)/(a2+a1+1.e-36)

! critical precip rate at which we assume the collector drops can change the
! drop size enough to enhance the auto-conversion process (mm/day)
   critpr = 0.5

   convfw = 1.94*2.13*sqrt(rhow*1000.*9.81*2.7e-4)

! liquid microphysics
!      cracw = 6                 ! beheng
   cracw = .884*sqrt(9.81/(rhow*1000.*2.7e-4)) ! tripoli and cotton

! ice microphysics
!      ciautb = 6.e-4
!      ciautb = 1.e-3
   ciautb = 5.e-4
!      icritw = 1.e-5
!      icritw = 5.e-5
   icritw = 4.e-4
!      icritc = 4.e-6
!      icritc = 6.e-6
   icritc = 5.e-6

   dtice = 20.
   rdtice = 1./dtice

   capnw = 400.              ! warm continental cloud particles / cm3
   capnc = 150.              ! cold and oceanic cloud particles / cm3
!  capnsi = 40.              ! sea ice cloud particles density  / cm3
   capnsi =  5.              ! sea ice cloud particles density  / cm3

   kconst = 1.18e6           ! const for terminal velocity

!      effc = 1.                 ! autoconv collection efficiency following boucher 96
!      effc = .55*0.05           ! autoconv collection efficiency following baker 93
   effc = 0.55                ! autoconv collection efficiency following tripoli and cotton
!   effc = 0.    ! turn off warm-cloud autoconv
   alpha = 1.1**4
   capc = pi**(-.333)*kconst*effc *(0.75)**(1.333)*alpha  ! constant for autoconversion

   r3lcrit = 15.0e-6         ! 15.0u  crit radius where liq conversion begins
!
! find all the points where we need to do the microphysics
! and set the output variables to zero
!
   ncols = 0
   do i = 1,ncol
      coef(i) = 0.
      fwaut(i) = 0.
      fsaut(i) = 0.
      fracw(i) = 0.
      fsacw(i) = 0.
      fsaci(i) = 0.
      liqmr(i) = 0.
      rainmr(i) = 0.
      if (cwm(i) > 1.e-20) then
         ncols = ncols + 1
         ind(ncols) = i
      endif
   end do

!cdir$ ivdep
   do ii = 1,ncols
      i = ind(ii)
!
! the local cloudiness at this level
!
      cldloc(i) = max(cldmin,cldm(i))
!
! a weighted mean between max cloudiness above, and this layer
!
      cldpr(i) = max(cldmin,(cldmax(i)+cldm(i))*0.5)
!
! decompose the suspended condensate into
! an incloud liquid and ice phase component
!
      totmr(i) = cwm(i)/cldloc(i)
      icemr(i) = totmr(i)*fice(i)
      liqmr(i) = totmr(i)*(1.-fice(i))
!
! density
!
      rho(i) = p(i,k)/(287.*t(i,k))
      rhocgs = rho(i)*1.e-3     ! density in cgs units
!
! decompose the precipitate into a liquid and ice phase
!
      if (t(i,k) > t0) then
         vfallw = convfw/sqrt(rho(i))
         rainmr(i) = precab(i)/(rho(i)*vfallw*cldpr(i))
         snowfr = 0
      else
         snowfr = 1
         rainmr(i) = 0.
      endif
!
! local snow amount in cgs units
!
      prscgs(i) = precab(i)/cldpr(i)*0.1*snowfr
!
! local rain amount in mm/day
!
      prlloc(i) = precab(i)*86400./cldpr(i)
   end do

   con1 = 1./(1.333*pi)**0.333 * 0.01 ! meters
!
! calculate the conversion terms
!
   call get_rlat_all_p(lchnk, ncol, rlat)

!cdir$ ivdep
   do ii = 1,ncols
      i = ind(ii)
      rhocgs = rho(i)*1.e-3     ! density in cgs units
!
! exponential temperature factor
!
!        efact = exp(0.025*(t(i,k)-t0))
!
! some temperature dependent constants
!
      wt = min(1._r8,max(0._r8,(t0-t(i,k))*rdtice))
      icrit = icritc*wt + icritw*(1-wt)
!
! linear weight factor in pressure (1 near sfc, 0 at .8 of sfc)
!
      wp = min(1._r8,max(0._r8,(p(i,k)-0.8*p(i,pver))/(0.2*p(i,pver))))
!
! near land near sfc raise the number concentration
! except south of 60S where land is so clean we treat it as ocean
!
      if (rlat(i) < -60./degrad) then
         wland = 0.
      else
         wland = landm(i,lchnk)
      endif
      
      capnoice =  wland*(capnw*wp + capnc*(1-wp))+ &
                  (1.-wland)*capnc
!
!     modify the estimated value to acknowledge cloud water number changes over sea ice
!
      wsi = icefrac(i) ! (1 all sea ice, 0, all ocean or land)
      capn = (capnsi*wsi + capnoice*(1-wsi))*wp + capnoice*(1-wp)

!!!!  enable the following line to get old cloudwater behavior
!!!!  capn = capnoice

      if (icefrac(i) > 0.001) then
         capn = capnsi
      else
         capn = capnoice
      endif

#ifdef DEBUG
      if ( (lat(i) == latlook(1)) .or. (lat(i) == latlook(2)) ) then
         if (i == ilook(1)) then
            write (6,*) ' findmcnew: lat, k, icefrac, landm, wp, capnoice, capn ', &
                 lat(i), k, icefrac(i), landm(i,lat(i)), wp, capnoice, capn
         endif
      endif
#endif

!
! useful terms in following calculations
!
      rat1 = rhocgs/rhow
      rat2 = liqmr(i)/capn
      con2 = (rat1*rat2)**0.333
!
! volume radius
!
!        r3l = (rhocgs*liqmr(i)/(1.333*pi*capn*rhow))**0.333 * 0.01 ! meters
      r3l = con1*con2
!
! critical threshold for autoconversion if modified for mixed phase
! clouds to mimic a bergeron findeisen process
! r3lc2 = r3lcrit*(1.-0.5*fice(i)*(1-fice(i)))
!
! autoconversion of liquid
!
!        cwaut = 2.e-4
!        cwaut = 1.e-3
!        lcrit = 2.e-4
!        lcrit = 5.e-4
!        pwaut = max(0._r8,liqmr(i)-lcrit)*cwaut
!
! pwaut is following tripoli and cotton (and many others)
! we reduce the autoconversion below critpr, because these are regions where
! the drop size distribution is likely to imply much smaller collector drops than
! those relevant for a cloud distribution corresponding to the value of effc = 0.55
! suggested by cotton (see austin 1995 JAS, baker 1993)

! easy to follow form
!        pwaut = capc*liqmr(i)**2*rhocgs/rhow
!    $           *(liqmr(i)*rhocgs/(rhow*capn))**(.333)
!    $           *heavy(r3l,r3lcrit)
!    $           *max(0.10_r8,min(1._r8,prlloc(i)/critpr))
! somewhat faster form
#define HEAVYNEW
#ifdef HEAVYNEW
!#ifdef PERGRO
      pwaut = capc*liqmr(i)**2*rat1*con2*heavymp(r3l,r3lcrit) * &
              max(0.10_r8,min(1._r8,prlloc(i)/critpr))
#else
      pwaut = capc*liqmr(i)**2*rat1*con2*heavym(r3l,r3lcrit)* &
              max(0.10_r8,min(1._r8,prlloc(i)/critpr))
#endif
!
! autoconversion of ice
!
!        ciaut = ciautb*efact
      ciaut = ciautb
!        psaut = capc*totmr(i)**2*rhocgs/rhoi
!     $           *(totmr(i)*rhocgs/(rhoi*capn))**(.333)
!
! autoconversion of ice condensate
!
#ifdef PERGRO
      psaut = heavyp(icemr(i),icrit)*icemr(i)*ciaut
#else
      psaut = max(0._r8,icemr(i)-icrit)*ciaut
#endif
!
! collection of liquid by rain
!
!        pracw = cracw*rho(i)*liqmr(i)*rainmr(i) !(beheng 1994)
      pracw = cracw*rho(i)*sqrt(rho(i))*liqmr(i)*rainmr(i) !(tripoli and cotton)
!!      pracw = 0.
!