#include <misc.h>
#include <params.h>
module cloud_fraction 2,1
!
! Module for cloud fraction routines.
!
! $Id: cloud_fraction.F90,v 1.1.4.7 2004/05/10 21:39:46 jmccaa Exp $
!
use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
private
save
!
! Public interfaces
!
public cldfrc_init ! Inititialization of cloud_fraction run-time parameters
public cldfrc ! Computation of cloud fraction
!
! Private data
!
real(r8) rhminl ! minimum rh for low stable clouds
real(r8) rhminh ! minimum rh for high stable clouds
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
!
contains
subroutine cldfrc_init() 1,1
!
! Purpose:
! Initialize cloud fraction run-time parameters
!
! Author: J. McCaa
!
use dycore, only: dycore_is, get_resolution
if ( dycore_is ('LR') ) then
if ( get_resolution() == '1x1.25' ) then
rhminl = .88
rhminh = .77
sh1 = 0.04
sh2 = 500.0
dp1 = 0.10
dp2 = 500.0
premit = 750.e2 ! top of area defined to be mid-level cloud
else
rhminl = .90
rhminh = .80
sh1 = 0.04
sh2 = 500.0
dp1 = 0.10
dp2 = 500.0
premit = 750.e2 ! top of area defined to be mid-level cloud
endif
else
if ( get_resolution() == 'T85' ) then
rhminl = .91
rhminh = .70
sh1 = 0.07
sh2 = 500.0
dp1 = 0.14
dp2 = 500.0
premit = 250.e2 ! top of area defined to be mid-level cloud
elseif ( get_resolution() == 'T31' ) then
rhminl = .88
rhminh = .80
sh1 = 0.07
sh2 = 500.0
dp1 = 0.14
dp2 = 500.0
premit = 750.e2 ! top of area defined to be mid-level cloud
else
rhminl = .90
rhminh = .80
sh1 = 0.07
sh2 = 500.0
dp1 = 0.14
dp2 = 500.0
premit = 750.e2 ! top of area defined to be mid-level cloud
endif
endif
end subroutine cldfrc_init
subroutine cldfrc(lchnk ,ncol , & 2,9
pmid ,temp ,q ,omga , phis, &
cldtop ,cldbot ,cloud ,clc ,pdel , &
cmfmc ,cmfmc2 ,landfrac,snowh ,concld ,cldst , &
ts ,sst ,ps ,zdu ,ocnfrac ,&
rhu00 ,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
use phys_grid, only: get_rlat_all_p, get_rlon_all_p
use dycore, only: dycore_is, get_resolution
implicit none
real(r8), parameter :: pnot = 1.e5 ! reference pressure
real(r8), parameter :: lapse = 6.5e-3 ! U.S. Standard Atmsophere lapse rate
real(r8), parameter :: premib = 750.e2 ! bottom pressure bound of middle cloud
real(r8), parameter :: pretop = 1.0e2 ! pressure bounding high cloud
!
! 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
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) :: cldtop(pcols) ! top level of convection
real(r8), intent(in) :: cldbot(pcols) ! bottom level of convection
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
!
! Output arguments
!
real(r8), intent(out) :: cloud(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) dmudp ! measure of mass detraining in a layer
!
!---------------------------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) deepcu ! deep convection cloud fraction
real(r8) shallowcu ! shallow convection cloud fraction
logical cldbnd(pcols) ! region below high cloud boundary
integer i, ierror, k ! column, level indices
integer kp1
integer kdthdp(pcols)
integer numkcld ! number of levels in which to allow clouds
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)
!==================================================================================
! PHILOSOPHY OF PRESENT IMPLEMENTATION
!
! 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
! define rh perturbation in order to estimate rhdfda
rhpert = 0.01
!
! Evaluate potential temperature and relative humidity
!
call aqsat(temp ,pmid ,es ,qs ,pcols , &
ncol ,pver ,1 ,pver )
do k=1,pver
do i=1,ncol
theta(i,k) = temp(i,k)*(pnot/pmid(i,k))**cappa
rh(i,k) = q(i,k)/qs(i,k)*(1.0+float(dindex)*rhpert)
!
! record relhum, rh itself will later be modified related with concld
!
relhum(i,k) = rh(i,k)
cloud(i,k) = 0.
cldst(i,k) = 0.
concld(i,k) = 0.
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) thetas(i) = &
( sst(i) - lapse * phis(i) / gravit) * (pnot/ps(i))**cappa
if(ocnfrac(i).gt.0.01.and.sst(i).lt.260.) ierror = i
clc(i) = 0.0
end do
coef1 = gravit*864.0 ! conversion to millibars/day
if (ierror > 0) then
write(6,*) 'COLDSST: encountered in cldfrc:', lchnk,ierror,ocnfrac(ierror),sst(ierror)
endif
do k=1,pver-1
do i=1,ncol
rpdeli(i,k) = 1./(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
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
shallowcu = max(0.0,min(sh1*log(1.0+sh2*cmfmc2(i,k+1)),0.30))
deepcu = max(0.0,min(dp1*log(1.0+dp2*(cmfmc(i,k+1)-cmfmc2(i,k+1))),0.60))
concld(i,k) = min(shallowcu + deepcu,0.80)
rh(i,k) = (rh(i,k) - concld(i,k))/(1.0 - 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
!
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)) then
rhlim = rhminl - 0.10
else
rhlim = rhminl
endif
!
rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
cloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
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)
cloud(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)) then
rhlim = rhminh*rhwght + (rhminl - 0.10)*(1.0-rhwght)
else
rhlim = rhminh*rhwght + rhminl*(1.0-rhwght)
endif
rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
cloud(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
!==================================================================================
end do
!
! Final evaluation of layered cloud fraction
!
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.
!
!==================================================================================
!
! 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
kdthdp(i) = 0
end do
!
do k=2,pver
do i=1,ncol
if (pmid(i,k) >= premib .and. ocnfrac(i).gt. 0.01) then
! I think this is done so that dtheta/dp is in units of dg/mb (JJH)
dthdp = 100.0*(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) then
dthdp = 100.0 * (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-.5573) ) )
!
! 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
!
! AGGREGATE CLOUD CONTRIBUTIONS (cldst should be zero everywhere except at level kdthdp(i))
!
do k=1,pver
do i=1,ncol
!
! which is greater; standard layered cloud amount or stratocumulus diagnosis
!
cloud(i,k) = max(cloud(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