module cloud_fraction 5,4
! Cloud fraction parameterization.
use shr_kind_mod
, only: r8 => shr_kind_r8
use spmd_utils, only: masterproc
use cam_logfile, only: iulog
use abortutils, only: endrun
implicit none
private
save
! Public interfaces
public &
cldfrc_readnl, &! read cldfrc_nl namelist
cldfrc_init, &! Inititialization of cloud_fraction run-time parameters
cldfrc_getparams, &! public access of tuning parameters
cldfrc ! Computation of cloud fraction
! Private data
real(r8), parameter :: unset_r8 = huge(1.0_r8)
! Namelist variables
logical :: cldfrc_freeze_dry ! switch for Vavrus correction
logical :: cldfrc_ice ! switch to compute ice cloud fraction
real(r8) :: cldfrc_rhminl = unset_r8 ! minimum rh for low stable clouds
real(r8) :: cldfrc_rhminh = unset_r8 ! minimum rh for high stable clouds
real(r8) :: rhminl ! set from namelist input cldfrc_rhminl
real(r8) :: rhminh ! set from namelist input cldfrc_rhminh
real(r8) :: sh1, sh2 ! parameters for shallow convection cloud fraction
real(r8) :: dp1,dp2 ! parameters for deep convection cloud fraction
real(r8) :: premit ! top pressure bound for mid level cloud
real(r8), parameter :: pnot = 1.e5_r8 ! reference pressure
real(r8), parameter :: lapse = 6.5e-3_r8 ! U.S. Standard Atmsophere lapse rate
real(r8), parameter :: premib_uw = 750.e2_r8 ! bottom pressure bound of middle cloud for UW
real(r8), parameter :: premib_cam = 750.e2_r8 ! bottom pressure bound of middle cloud for CAM
real(r8) :: premib ! bottom pressure bound of middle cloud
real(r8), parameter :: pretop = 1.0e2_r8 ! pressure bounding high cloud
integer :: iceopt = 4 ! option for ice cloud closure
! 1=wang & sassen 2=schiller (iciwc)
! 3=wood & field, 4=Wilson (based on smith)
real(r8), parameter :: icecrit = 0.95_r8 ! Critical RH for ice clouds in Wilson & Ballard closure (smaller = more ice clouds)
integer count
character(len=16) :: microp_scheme
logical :: inversion_cld_off ! Turns off stratification-based cld frc
!================================================================================================
contains
!================================================================================================
subroutine cldfrc_readnl(nlfile) 1,11
use namelist_utils, only: find_group_name
use units, only: getunit, freeunit
use mpishorthand
character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
! Local variables
integer :: unitn, ierr
character(len=*), parameter :: subname = 'cldfrc_readnl'
namelist /cldfrc_nl/ cldfrc_freeze_dry, cldfrc_ice, cldfrc_rhminl, cldfrc_rhminh
!-----------------------------------------------------------------------------
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'cldfrc_nl', status=ierr)
if (ierr == 0) then
read(unitn, cldfrc_nl, iostat=ierr)
if (ierr /= 0) then
call endrun(subname // ':: ERROR reading namelist')
end if
end if
close(unitn)
call freeunit(unitn)
! set local variables
rhminl = cldfrc_rhminl
rhminh = cldfrc_rhminh
end if
#ifdef SPMD
! Broadcast namelist variables
call mpibcast(cldfrc_freeze_dry, 1, mpilog, 0, mpicom)
call mpibcast(cldfrc_ice, 1, mpilog, 0, mpicom)
call mpibcast(rhminl, 1, mpir8, 0, mpicom)
call mpibcast(rhminh, 1, mpir8, 0, mpicom)
#endif
end subroutine cldfrc_readnl
!================================================================================================
subroutine cldfrc_getparams(rhminl_out, rhminh_out ) 1
!-----------------------------------------------------------------------
! Purpose: Return cldfrc tuning parameters
!-----------------------------------------------------------------------
real(r8), intent(out), optional :: rhminl_out
real(r8), intent(out), optional :: rhminh_out
if ( present(rhminl_out) ) rhminl_out = rhminl
if ( present(rhminh_out) ) rhminh_out = rhminh
end subroutine cldfrc_getparams
!===============================================================================
subroutine cldfrc_init 1,7
!
! Purpose:
! Initialize cloud fraction run-time parameters
!
! Author: J. McCaa
!
use dycore, only: dycore_is, get_resolution
use ppgrid, only: pver
use chemistry, only: chem_is
use phys_control, only: phys_getopts
! horizontal grid specifier
character(len=32) :: hgrid
! query interfaces for scheme settings
character(len=16) :: shallow_scheme, eddy_scheme
call phys_getopts(shallow_scheme_out = shallow_scheme ,&
eddy_scheme_out = eddy_scheme, &
microp_scheme_out = microp_scheme)
hgrid = get_resolution()
if ( shallow_scheme .eq. 'UW' ) then
premib = premib_uw
else
premib = premib_cam
endif
sh2 = 500.0_r8
dp2 = 500.0_r8
premit = 750.e2_r8
iceopt = 1
if ( dycore_is ('LR') ) then
sh1 = 0.04_r8
dp1 = 0.10_r8
if ( trim(hgrid) == '0.23x0.31' ) then
premit = 250.e2_r8
elseif ( trim(hgrid) == '0.47x0.63' .or. trim(hgrid) == '0.5x0.625' ) then
premit = 250.e2_r8
elseif ( trim(hgrid) == '1.9x2.5' .or. trim(hgrid) == '2x2.5' ) then
if ( microp_scheme .eq. 'MG' .and. shallow_scheme .eq. 'UW') then
dp1 = 0.04_r8
dp2 = 675.0_r8
endif
endif
else
sh1 = 0.07_r8
dp1 = 0.14_r8
if ( trim(hgrid) == 'T341' ) then
premit = 250.e2_r8 ! top of area defined to be mid-level cloud
elseif ( trim(hgrid) == 'T170' ) then
premit = 250.e2_r8
elseif ( trim(hgrid) == 'T85' ) then
premit = 250.e2_r8
endif
endif
! Turn off inversion_cld if any UW PBL scheme is being used
if ( (eddy_scheme .eq. 'diag_TKE' ) .or. (shallow_scheme .eq. 'UW' )) then
inversion_cld_off = .true.
else
inversion_cld_off = .false.
endif
if ( masterproc ) then
write(iulog,*)'tuning parameters cldfrc_init: inversion_cld_off',inversion_cld_off
write(iulog,*)'tuning parameters cldfrc_init: dp1',dp1,'dp2',dp2,'sh1',sh1,'sh2',sh2
if (shallow_scheme .ne. 'UW' ) then
write(iulog,*)'tuning parameters cldfrc_init: rhminl',rhminl,'rhminh',rhminh,'premit',premit,'premib',premib
endif
endif
end subroutine cldfrc_init
subroutine cldfrc(lchnk ,ncol , pbuf, & 4,15
pmid ,temp ,q ,omga , phis, &
shfrc ,use_shfrc, &
cloud ,rhcloud, clc ,pdel , &
cmfmc ,cmfmc2 ,landfrac,snowh ,concld ,cldst , &
ts ,sst ,ps ,zdu ,ocnfrac ,&
rhu00 ,cldice ,icecldf ,liqcldf ,relhum ,dindex )
!-----------------------------------------------------------------------
!
! Purpose:
! Compute cloud fraction
!
!
! Method:
! This calculate cloud fraction using a relative humidity threshold
! The threshold depends upon pressure, and upon the presence or absence
! of convection as defined by a reasonably large vertical mass flux
! entering that layer from below.
!
! Author: Many. Last modified by Jim McCaa
!
!-----------------------------------------------------------------------
use ppgrid
use physconst, only: cappa, gravit, rair, tmelt
use cldconst
use wv_saturation, only: aqsat, aqsat_water, polysvp
use phys_grid, only: get_rlat_all_p, get_rlon_all_p
use dycore, only: dycore_is, get_resolution
!RBN - Need this to write shallow,deep fraction to phys buffer.
!PJR - we should probably make seperate modules for determining convective
! clouds and make this one just responsible for relative humidity clouds
use phys_buffer, only: pbuf_size_max, pbuf_fld, pbuf_old_tim_idx, pbuf_get_fld_idx
! Arguments
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
integer, intent(in) :: dindex ! 0 or 1 to perturb rh
type(pbuf_fld), intent(inout) :: pbuf(pbuf_size_max)
real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures
real(r8), intent(in) :: temp(pcols,pver) ! temperature
real(r8), intent(in) :: q(pcols,pver) ! specific humidity
real(r8), intent(in) :: omga(pcols,pver) ! vertical pressure velocity
real(r8), intent(in) :: cmfmc(pcols,pverp) ! convective mass flux--m sub c
real(r8), intent(in) :: cmfmc2(pcols,pverp) ! shallow convective mass flux--m sub c
real(r8), intent(in) :: snowh(pcols) ! snow depth (liquid water equivalent)
real(r8), intent(in) :: pdel(pcols,pver) ! pressure depth of layer
real(r8), intent(in) :: landfrac(pcols) ! Land fraction
real(r8), intent(in) :: ocnfrac(pcols) ! Ocean fraction
real(r8), intent(in) :: ts(pcols) ! surface temperature
real(r8), intent(in) :: sst(pcols) ! sea surface temperature
real(r8), intent(in) :: ps(pcols) ! surface pressure
real(r8), intent(in) :: zdu(pcols,pver) ! detrainment rate from deep convection
real(r8), intent(in) :: phis(pcols) ! surface geopotential
real(r8), intent(in) :: shfrc(pcols,pver) ! cloud fraction from convect_shallow
real(r8), intent(in) :: cldice(pcols,pver) ! cloud ice mixing ratio
logical, intent(in) :: use_shfrc
! Output arguments
real(r8), intent(out) :: cloud(pcols,pver) ! cloud fraction
real(r8), intent(out) :: rhcloud(pcols,pver) ! cloud fraction
real(r8), intent(out) :: clc(pcols) ! column convective cloud amount
real(r8), intent(out) :: cldst(pcols,pver) ! cloud fraction
real(r8), intent(out) :: rhu00(pcols,pver) ! RH threshold for cloud
real(r8), intent(out) :: relhum(pcols,pver) ! RH
real(r8), intent(out) :: icecldf(pcols,pver) ! ice cloud fraction
real(r8), intent(out) :: liqcldf(pcols,pver) ! liquid cloud fraction (combined into cloud)
!---------------------------Local workspace-----------------------------
!
real(r8) concld(pcols,pver) ! convective cloud cover
real(r8) cld ! intermediate scratch variable (low cld)
real(r8) dthdpmn(pcols) ! most stable lapse rate below 750 mb
real(r8) dthdp ! lapse rate (intermediate variable)
real(r8) es(pcols,pver) ! saturation vapor pressure
real(r8) qs(pcols,pver) ! saturation specific humidity
real(r8) rhwght ! weighting function for rhlim transition
real(r8) rh(pcols,pver) ! relative humidity
real(r8) rhdif ! intermediate scratch variable
real(r8) strat ! intermediate scratch variable
real(r8) theta(pcols,pver) ! potential temperature
real(r8) rhlim ! local rel. humidity threshold estimate
real(r8) coef1 ! coefficient to convert mass flux to mb/d
real(r8) clrsky(pcols) ! temporary used in random overlap calc
real(r8) rpdeli(pcols,pver-1) ! 1./(pmid(k+1)-pmid(k))
real(r8) rhpert !the specified perturbation to rh
real(r8), pointer, dimension(:,:) :: deepcu ! deep convection cloud fraction
real(r8), pointer, dimension(:,:) :: shallowcu ! shallow convection cloud fraction
logical cldbnd(pcols) ! region below high cloud boundary
integer i, ierror, k ! column, level indices
integer kp1, ifld
integer kdthdp(pcols)
integer numkcld ! number of levels in which to allow clouds
! In Cloud Ice Content variables
real(r8) :: a,b,c,as,bs,cs !fit parameters
real(r8) :: Kc !constant for ice cloud calc (wood & field)
real(r8) :: ttmp !limited temperature
real(r8) :: icicval !empirical iwc value
real(r8) :: rho !local air density
real(r8) :: esl(pcols,pver) !liq sat vapor pressure
real(r8) :: esi(pcols,pver) !ice sat vapor pressure
real(r8) :: ncf,phi !Wilson and Ballard parameters
real(r8) thetas(pcols) ! ocean surface potential temperature
real(r8) :: clat(pcols) ! current latitudes(radians)
real(r8) :: clon(pcols) ! current longitudes(radians)
! Statement functions
logical land
land(i) = nint(landfrac(i)) == 1
call get_rlat_all_p(lchnk, ncol, clat)
call get_rlon_all_p(lchnk, ncol, clon)
ifld = pbuf_get_fld_idx('SH_FRAC')
shallowcu => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
ifld = pbuf_get_fld_idx('DP_FRAC')
deepcu => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
! Initialise cloud fraction
shallowcu = 0._r8
deepcu = 0._r8
!==================================================================================
! PHILOSOPHY OF PRESENT IMPLEMENTATION
!++ag ice3
! Modification to philosophy for ice supersaturation
! philosophy below is based on RH water only. This is 'liquid condensation'
! or liquid cloud (even though it will freeze immediately to ice)
! The idea is that the RH limits for condensation are strict only for
! water saturation
!
! Ice clouds are formed by explicit parameterization of ice nucleation.
! Closure for ice cloud fraction is done on available cloud ice, such that
! the in-cloud ice content matches an empirical fit
! thus, icecldf = min(cldice/icicval,1) where icicval = f(temp,cldice,numice)
! for a first cut, icicval=f(temp) only.
! Combined cloud fraction is maximum overlap cloud=max(1,max(icecldf,liqcldf))
! No dA/dt term for ice?
!--ag
!
! There are three co-existing cloud types: convective, inversion related low-level
! stratocumulus, and layered cloud (based on relative humidity). Layered and
! stratocumulus clouds do not compete with convective cloud for which one creates
! the most cloud. They contribute collectively to the total grid-box average cloud
! amount. This is reflected in the way in which the total cloud amount is evaluated
! (a sum as opposed to a logical "or" operation)
!
!==================================================================================
! set defaults for rhu00
rhu00(:,:) = 2.0_r8
! define rh perturbation in order to estimate rhdfda
rhpert = 0.01_r8
!set Wang and Sassen IWC paramters
a=26.87_r8
b=0.569_r8
c=0.002892_r8
!set schiller parameters
as=-68.4202_r8
bs=0.983917_r8
cs=2.81795_r8
!set wood and field paramters...
Kc=75._r8
! Evaluate potential temperature and relative humidity
! If not computing ice cloud fraction then hybrid RH, if MG then water RH
if ( .not. cldfrc_ice ) then
call aqsat(temp, pmid, es, qs, pcols, &
ncol, pver, 1, pver)
else
call aqsat_water(temp, pmid, es, qs, pcols, &
ncol, pver, 1, pver)
endif
do k=1,pver
do i=1,ncol
if (cldfrc_ice) then
!++ag calculate qsat ice from qsatw
esl(i,k)=polysvp(temp(i,k),0)
esi(i,k)=polysvp(temp(i,k),1)
end if
theta(i,k) = temp(i,k)*(pnot/pmid(i,k))**cappa
rh(i,k) = q(i,k)/qs(i,k)*(1.0_r8+real(dindex,r8)*rhpert)
! record relhum, rh itself will later be modified related with concld
relhum(i,k) = rh(i,k)
cloud(i,k) = 0._r8
icecldf(i,k) = 0._r8
liqcldf(i,k) = 0._r8
rhcloud(i,k) = 0._r8
cldst(i,k) = 0._r8
concld(i,k) = 0._r8
end do
end do
! Initialize other temporary variables
ierror = 0
do i=1,ncol
! Adjust thetas(i) in the presence of non-zero ocean heights.
! This reduces the temperature for positive heights according to a standard lapse rate.
if(ocnfrac(i).gt.0.01_r8) thetas(i) = &
( sst(i) - lapse * phis(i) / gravit) * (pnot/ps(i))**cappa
if(ocnfrac(i).gt.0.01_r8.and.sst(i).lt.260._r8) ierror = i
clc(i) = 0.0_r8
end do
coef1 = gravit*864.0_r8 ! conversion to millibars/day
if (ierror > 0) then
write(iulog,*) 'COLDSST: encountered in cldfrc:', lchnk,ierror,ocnfrac(ierror),sst(ierror)
endif
do k=1,pver-1
do i=1,ncol
rpdeli(i,k) = 1._r8/(pmid(i,k+1) - pmid(i,k))
end do
end do
!
! Estimate of local convective cloud cover based on convective mass flux
! Modify local large-scale relative humidity to account for presence of
! convective cloud when evaluating relative humidity based layered cloud amount
!
do k=1,pver
do i=1,ncol
concld(i,k) = 0.0_r8
end do
end do
!
! cloud mass flux in SI units of kg/m2/s; should produce typical numbers of 20%
! shallow and deep convective cloudiness are evaluated separately (since processes
! are evaluated separately) and summed
!
#ifndef PERGRO
do k=1,pver-1
do i=1,ncol
if ( .not. use_shfrc ) then
shallowcu(i,k) = max(0.0_r8,min(sh1*log(1.0_r8+sh2*cmfmc2(i,k+1)),0.30_r8))
else
shallowcu(i,k) = shfrc(i,k)
endif
deepcu(i,k) = max(0.0_r8,min(dp1*log(1.0_r8+dp2*(cmfmc(i,k+1)-cmfmc2(i,k+1))),0.60_r8))
concld(i,k) = min(shallowcu(i,k) + deepcu(i,k),0.80_r8)
rh(i,k) = (rh(i,k) - concld(i,k))/(1.0_r8 - concld(i,k))
end do
end do
#endif
!==================================================================================
!
! ****** Compute layer cloudiness ******
!
!====================================================================
! Begin the evaluation of layered cloud amount based on (modified) RH
!====================================================================
!
numkcld = pver
do k=2,numkcld
kp1 = min(k + 1,pver)
do i=1,ncol
!++ag This is now designed to apply FOR LIQUID CLOUDS (condensation > RH water)
cldbnd(i) = pmid(i,k).ge.pretop
if ( pmid(i,k).ge.premib ) then
!==============================================================
! This is the low cloud (below premib) block
!==============================================================
! enhance low cloud activation over land with no snow cover
if (land(i) .and. (snowh(i) <= 0.000001_r8)) then
rhlim = rhminl - 0.10_r8
! rhlim = rhminl
else
rhlim = rhminl
endif
rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
! SJV: decrease cloud amount if very low water vapor content
! (thus very cold): "freeze dry"
if (cldfrc_freeze_dry) then
rhcloud(i,k) = rhcloud(i,k)*max(0.15_r8,min(1.0_r8,q(i,k)/0.0030_r8))
endif
else if ( pmid(i,k).lt.premit ) then
!==============================================================
! This is the high cloud (above premit) block
!==============================================================
!
rhlim = rhminh
!
rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
else
!==============================================================
! This is the middle cloud block
!==============================================================
!
! linear rh threshold transition between thresholds for low & high cloud
!
rhwght = (premib-(max(pmid(i,k),premit)))/(premib-premit)
if (land(i) .and. (snowh(i) <= 0.000001_r8)) then
rhlim = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
else
rhlim = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
endif
rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
end if
!==================================================================================
! WE NEED TO DOCUMENT THE PURPOSE OF THIS TYPE OF CODE (ASSOCIATED WITH 2ND CALL)
!==================================================================================
! !
! ! save rhlim to rhu00, it handles well by itself for low/high cloud
! !
rhu00(i,k)=rhlim
!==================================================================================
if (cldfrc_ice) then
! Evaluate ice cloud fraction based on in-cloud ice content
!--------ICE CLOUD OPTION 1--------Wang & Sassen 2002
! Evaluate desired in-cloud water content
! icicval = f(temp,cldice,numice)
! Start with a function of temperature.
! Wang & Sassen 2002 (JAS), based on ARM site MMCR (midlat cirrus)
! parameterization valid for 203-253K
! icival > 0 for t>195K
if (iceopt.lt.3) then
if (iceopt.eq.1) then
ttmp=max(195._r8,min(temp(i,k),253._r8)) - 273.16_r8
icicval=a + b * ttmp + c * ttmp**2._r8
!convert units
rho=pmid(i,k)/(rair*temp(i,k))
icicval= icicval * 1.e-6_r8 / rho
else
!--------ICE CLOUD OPTION 2--------Schiller 2008 (JGR)
! Use a curve based on FISH measurements in
! tropics, mid-lats and arctic. Curve is for 180-250K (raise to 273K?)
! use median all flights
ttmp=max(190._r8,min(temp(i,k),273.16_r8))
icicval = 10._r8 **(as * bs**ttmp + cs)
!convert units from ppmv to kg/kg
icicval= icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
endif
!set icecldfraction for OPTION 1 or OPTION2
icecldf(i,k) = max(0._r8,min(cldice(i,k)/icicval,1._r8))
else if (iceopt.eq.3) then
!--------ICE CLOUD OPTION 3--------Wood & Field 2000 (JAS)
! eq 6: cloud fraction = 1 - exp (-K * qc/qsati)
icecldf(i,k)=1._r8 - exp(-Kc*cldice(i,k)/(qs(i,k)*(esi(i,k)/esl(i,k))))
icecldf(i,k)=max(0._r8,min(icecldf(i,k),1._r8))
else
!--------ICE CLOUD OPTION 4--------Wilson and ballard 1999
! inversion of smith....
! ncf = cldice / ((1-RHcrit)*qs)
! then a function of ncf....
ncf =cldice(i,k)/((1._r8 - icecrit)*qs(i,k))
if (ncf.le.0._r8) then
icecldf(i,k)=0._r8
else if (ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8) then
icecldf(i,k)=0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
else if (ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8) then
phi=(acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
icecldf(i,k)=(1._r8 - 4._r8 * cos(phi) * cos(phi))
else
icecldf(i,k)=1._r8
endif
icecldf(i,k)=max(0._r8,min(icecldf(i,k),1._r8))
endif
!TEST: if ice present, icecldf=1.
! if (cldice(i,k).ge.1.e-8_r8) then
! icecldf(i,k) = 0.99_r8
! endif
!! if ((cldice(i,k) .gt. icicval) .or. ((cldice(i,k) .gt. 0._r8) .and. (icecldf(i,k) .eq. 0._r8))) then
! if (cldice(i,k) .gt. 1.e-8_r8) then
! write(iulog,*) 'i,k,pmid,rho,t,cldice,icicval,icecldf,rhcloud: ', &
! i,k,pmid(i,k),rho,temp(i,k),cldice(i,k),icicval,icecldf(i,k),rhcloud(i,k)
! endif
! Combine ice and liquid cloud fraction assuming maximum overlap.
! Combined cloud fraction is maximum overlap
! cloud(i,k)=min(1._r8,max(icecldf(i,k),rhcloud(i,k)))
liqcldf(i,k)=(1._r8 - icecldf(i,k))* rhcloud(i,k)
cloud(i,k)=liqcldf(i,k) + icecldf(i,k)
else
! For RK microphysics
cloud(i,k) = rhcloud(i,k)
end if
end do
end do
!
! Add in the marine strat
! MARINE STRATUS SHOULD BE A SPECIAL CASE OF LAYERED CLOUD
! CLOUD CURRENTLY CONTAINS LAYERED CLOUD DETERMINED BY RH CRITERIA
! TAKE THE MAXIMUM OF THE DIAGNOSED LAYERED CLOUD OR STRATOCUMULUS
!
!===================================================================================
!
! SOME OBSERVATIONS ABOUT THE FOLLOWING SECTION OF CODE (missed in earlier look)
! K700 IS SET AS A CONSTANT BASED ON HYBRID COORDINATE: IT DOES NOT DEPEND ON
! LOCAL PRESSURE; THERE IS NO PRESSURE RAMP => LOOKS LEVEL DEPENDENT AND
! DISCONTINUOUS IN SPACE (I.E., STRATUS WILL END SUDDENLY WITH NO TRANSITION)
!
! IT APPEARS THAT STRAT IS EVALUATED ACCORDING TO KLEIN AND HARTMANN; HOWEVER,
! THE ACTUAL STRATUS AMOUNT (CLDST) APPEARS TO DEPEND DIRECTLY ON THE RH BELOW
! THE STRONGEST PART OF THE LOW LEVEL INVERSION.
!PJR answers: 1) the rh limitation is a physical/mathematical limitation
! cant have more cloud than there is RH
! allowed the cloud to exist two layers below the inversion
! because the numerics frequently make 50% relative humidity
! in level below the inversion which would allow no cloud
! 2) since the cloud is only allowed over ocean, it should
! be very insensitive to surface pressure (except due to
! spectral ringing, which also causes so many other problems
! I didnt worry about it.
!
!==================================================================================
if (.not.inversion_cld_off) then
!
! Find most stable level below 750 mb for evaluating stratus regimes
!
do i=1,ncol
! Nothing triggers unless a stability greater than this minimum threshold is found
dthdpmn(i) = -0.125_r8
kdthdp(i) = 0
end do
!
do k=2,pver
do i=1,ncol
if (pmid(i,k) >= premib .and. ocnfrac(i).gt. 0.01_r8) then
! I think this is done so that dtheta/dp is in units of dg/mb (JJH)
dthdp = 100.0_r8*(theta(i,k) - theta(i,k-1))*rpdeli(i,k-1)
if (dthdp < dthdpmn(i)) then
dthdpmn(i) = dthdp
kdthdp(i) = k ! index of interface of max inversion
end if
end if
end do
end do
! Also check between the bottom layer and the surface
! Only perform this check if the criteria were not met above
do i = 1,ncol
if ( kdthdp(i) .eq. 0 .and. ocnfrac(i).gt.0.01_r8) then
dthdp = 100.0_r8 * (thetas(i) - theta(i,pver)) / (ps(i)-pmid(i,pver))
if (dthdp < dthdpmn(i)) then
dthdpmn(i) = dthdp
kdthdp(i) = pver ! index of interface of max inversion
endif
endif
enddo
do i=1,ncol
if (kdthdp(i) /= 0) then
k = kdthdp(i)
kp1 = min(k+1,pver)
! Note: strat will be zero unless ocnfrac > 0.01
strat = min(1._r8,max(0._r8, ocnfrac(i) * ((theta(i,k700)-thetas(i))*.057_r8-.5573_r8) ) )
!
! assign the stratus to the layer just below max inversion
! the relative humidity changes so rapidly across the inversion
! that it is not safe to just look immediately below the inversion
! so limit the stratus cloud by rh in both layers below the inversion
!
cldst(i,k) = min(strat,max(rh(i,k),rh(i,kp1)))
end if
end do
end if ! .not.inversion_cld_off
do k=1,pver
do i=1,ncol
!
! which is greater; standard layered cloud amount or stratocumulus diagnosis
!
cloud(i,k) = max(rhcloud(i,k),cldst(i,k))
!
! add in the contributions of convective cloud (determined separately and accounted
! for by modifications to the large-scale relative humidity.
!
cloud(i,k) = min(cloud(i,k)+concld(i,k), 1.0_r8)
end do
end do
!
return
end subroutine cldfrc
end module cloud_fraction