#include <misc.h>#include <params.h>subroutine cldfrc(lchnk   ,ncol    , &                  pmid    ,temp    ,q       ,omga    , &                  cldtop  ,cldbot  ,cloud   ,clc     ,pdel    , &                  cmfmc   ,landfrac,snowh   ,concld  ,cldst   , &                  ts      ,ps      ,zdu     ,ocnfrac ,&                  rhu00   ,relhum  ,dindex )!----------------------------------------------------------------------- ! ! Purpose: ! Compute cloud fraction using scheme of J.M.Slingo,! as modified by J.J.Hack and J.T.Kiehl! ! Method: ! This scheme is based on the operational scheme used in the ECMWF model! A full description of its development can be found in Slingo (1987),! which appears in the QJRMS July issue.  A number of modifications have! been introduced to the original scheme in the following implementation! ! Author: J. Hack! !-----------------------------------------------------------------------   use precision   use ppgrid   use physconst, only: cappa, gravit, rair   use cldconst   use wv_saturation, only: aqsat   use dycore   implicit none   real(r8), parameter :: pnot = 1.e5                 ! reference pressure!! 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) :: 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) :: ps(pcols)             ! surface pressure   real(r8), intent(in) :: zdu(pcols,pver)       ! detrainment rate from deep convection!! 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) cld8(pcols)           ! low cloud fraction estimate   real(r8) cld9(pcols)           ! mid and high cloud fraction estimate#ifdef STDCONCLD   real(r8) cck(pcols)            ! convective cloud per level (assuming!                                  random overlap in convective layer)   real(r8) zrth                  ! reciprocal of no. of convective layers   real(r8) ccldt(pcols)          ! estimate of total convective cloud#endif   real(r8) dthtdp(pcols,pver)    ! lapse rate (d theta/dp) below 750 mb   real(r8) dtdpmn(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) premib                ! bottom pressure bound of middle cloud   real(r8) pretop                ! pressure bounding high cloud   real(r8) rh(pcols,pver)        ! relative humidity#ifdef OLDLOWCLD   real(r8) rhb                   ! intermediate scratch variable   real(r8) pdepth                ! intermediate scratch variable   real(r8) stratfac              ! intermediate scratch variable#endif   real(r8) rhdif                 ! intermediate scratch variable   real(r8) strat                 ! intermediate scratch variable   real(r8) theta(pcols,pver)     ! potential temperature   real(r8) bvf                   ! brunt-vaisalla frequency   real(r8) rbvflim               ! bound on inverse of bvf   real(r8) rho                   ! local density (used to calculate bvf)   real(r8) rhlim                 ! local rel. humidity threshold estimate   real(r8) rhden                 ! intermediate scratch variable   real(r8) rhdif2                ! intermediate scratch variable   real(r8) rhminl                ! minimum rh for low stable clouds   real(r8) rhminh                ! minimum rh for high stable clouds   real(r8) mcbar(pcols)          ! mean convective scale motion in column   real(r8) dpsum(pcols)          ! vertical sum of delta-p (k-1 levels)   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   logical lol(pcols)             ! region of low level cloud   logical cldbnd(pcols)          ! region below high cloud boundary   integer i,k                    ! longitude, level indices   integer kp1   integer kdthdp(pcols)   integer numkcld                ! number of levels in which to allow clouds   real(r8) thetas(pcols)!! Statement functions!   logical land   logical ocean   land(i) = nint(landfrac(i)) == 1   ocean(i) = nint(ocnfrac(i)) == 1!! Set bound for inverse of brunt-vaisalla frequency and minimum relative! humidity thresholds for stable clouds.  These are the principal! "disposable" parameters for the cloud fraction scheme!   rbvflim = 1./0.00035! set defaults for rhu00   rhu00(:,:) = 2.0   if ( dycore_is ('LR') ) then        rhminl = .90   else        rhminl = .85   endif   rhminh = .90!! 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!   do i=1,ncol      thetas(i)  = ts(i)*(pnot/ps(i))**cappa      clc(i) = 0.0   end do   coef1 = gravit*864.0    ! conversion to millibars/day   do i=1,ncol      mcbar(i) = 0.0      dpsum(i) = 0.0   end do   do k=1,pver-1      do i=1,ncol         rpdeli(i,k) = 1./(pmid(i,k+1) - pmid(i,k))      end do   end do!! Calculate mean convective motion throughout column (in units of mb/day)!   do k=1,pver-1      do i=1,ncol         mcbar(i) = mcbar(i) + max(cmfmc(i,k+1)*coef1,0._r8)*pdel(i,k)         dpsum(i) = dpsum(i) + pdel(i,k)      end do   end do!! Estimate of total convective cloud cover based on mean convective motion!#ifdef STDCONCLD   do i=1,ncol      cck(i) = 0.0      mcbar(i) = max(mcbar(i)/dpsum(i),1.0e-15_r8)      ccldt(i) = min(0.035*log(1.0+mcbar(i)),0.80_r8)      if ((cldbot(i) - cldtop(i)) >= 1.0) then!! Inverse of depth of convection (depth is expressed in model levels)!         zrth = 1.0/(cldbot(i) - cldtop(i))!! Compute amount of convective cloud at each level so that! after random overlap, the total convective cloud cover is ccldt!         cck(i) = 1.0 - (1.0 - ccldt(i))**zrth      end if   end do!! Vertically distribute cloud in convective layer!   do k=1,pver-1      do i=1,ncol         if (k <= cldbot(i) .and. k >= cldtop(i)) then            concld(i,k) = cck(i)            rh(i,k) = (rh(i,k) - concld(i,k))/(1.0 - concld(i,k))         end if      end do   end do#else!     make the convective cloud depend on the conv. mass detraining!     for upper levels only (above 500mb), since Xu and Kreuger showed!     rh is a very poor predictor of those clouds   do k = 1,pver-1      do i = 1,ncol         if (pmid(i,k) < 5.e4) then!               dmudp = (cmfmc(i,k+1)-cmfmc(i,k))/pdel(i,k)