LCOV - code coverage report
Current view: | top level - physics/cam - cloud_fraction.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 88 | 229 | 38.4 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 5 | 6 | 83.3 %

              Line data    Source code
           1             : module cloud_fraction
           2             :
           3             :   ! Cloud fraction parameterization.
           4             :
           5             :
           6             :   use shr_kind_mod,   only: r8 => shr_kind_r8
           7             :   use ppgrid,         only: pcols, pver, pverp
           8             :   use ref_pres,       only: pref_mid
           9             :   use spmd_utils,     only: masterproc
          10             :   use cam_logfile,    only: iulog
          11             :   use cam_abortutils, only: endrun
          12             :   use ref_pres,       only: trop_cloud_top_lev
          13             :
          14             :   implicit none
          15             :   private
          16             :   save
          17             :
          18             :   ! Public interfaces
          19             :   public &
          20             :      cldfrc_readnl,    &! read cldfrc_nl namelist
          21             :      cldfrc_register,  &! add fields to pbuf
          22             :      cldfrc_init,      &! Inititialization of cloud_fraction run-time parameters
          23             :      cldfrc_getparams, &! public access of tuning parameters
          24             :      cldfrc,           &! Computation of cloud fraction
          25             :      cldfrc_fice,      &! Calculate fraction of condensate in ice phase (radiation partitioning)
          26             :      dp1,              &! parameter for deep convection cloud fraction needed in clubb_intr
          27             :      dp2               ! parameter for deep convection cloud fraction needed in clubb_intr
          28             :
          29             :   ! Private data
          30             :   real(r8), parameter :: unset_r8 = huge(1.0_r8)
          31             :
          32             :   ! Top level
          33             :   integer :: top_lev = 1
          34             :
          35             :   ! Physics buffer indices
          36             :   integer :: sh_frac_idx   = 0
          37             :   integer :: dp_frac_idx   = 0
          38             :
          39             :   ! Namelist variables
          40             :   logical  :: cldfrc_freeze_dry           ! switch for Vavrus correction
          41             :   logical  :: cldfrc_ice                  ! switch to compute ice cloud fraction
          42             :   real(r8) :: cldfrc_rhminl = unset_r8    ! minimum rh for low stable clouds
          43             :   real(r8) :: cldfrc_rhminl_adj_land = unset_r8   ! rhminl adjustment for snowfree land
          44             :   real(r8) :: cldfrc_rhminh = unset_r8    ! minimum rh for high stable clouds
          45             :   real(r8) :: cldfrc_sh1    = unset_r8    ! parameter for shallow convection cloud fraction
          46             :   real(r8) :: cldfrc_sh2    = unset_r8    ! parameter for shallow convection cloud fraction
          47             :   real(r8) :: cldfrc_dp1    = unset_r8    ! parameter for deep convection cloud fraction
          48             :   real(r8) :: cldfrc_dp2    = unset_r8    ! parameter for deep convection cloud fraction
          49             :   real(r8) :: cldfrc_premit = unset_r8    ! top pressure bound for mid level cloud
          50             :   real(r8) :: cldfrc_premib  = unset_r8   ! bottom pressure bound for mid level cloud
          51             :   integer  :: cldfrc_iceopt               ! option for ice cloud closure
          52             :                                           ! 1=wang & sassen 2=schiller (iciwc)
          53             :                                           ! 3=wood & field, 4=Wilson (based on smith)
          54             :   real(r8) :: cldfrc_icecrit = unset_r8   ! Critical RH for ice clouds in Wilson & Ballard closure (smaller = more ice clouds)
          55             :
          56             :   real(r8) :: rhminl             ! set from namelist input cldfrc_rhminl
          57             :   real(r8) :: rhminl_adj_land    ! set from namelist input cldfrc_rhminl_adj_land
          58             :   real(r8) :: rhminh             ! set from namelist input cldfrc_rhminh
          59             :   real(r8) :: sh1, sh2           ! set from namelist input cldfrc_sh1, cldfrc_sh2
          60             :   real(r8) :: dp1,dp2            ! set from namelist input cldfrc_dp1, cldfrc_dp2
          61             :   real(r8) :: premit             ! set from namelist input cldfrc_premit
          62             :   real(r8) :: premib             ! set from namelist input cldfrc_premib
          63             :   integer  :: iceopt             ! set from namelist input cldfrc_iceopt
          64             :   real(r8) :: icecrit            ! set from namelist input cldfrc_icecrit
          65             :
          66             :   ! constants
          67             :   real(r8), parameter :: pnot = 1.e5_r8         ! reference pressure
          68             :   real(r8), parameter :: lapse = 6.5e-3_r8      ! U.S. Standard Atmosphere lapse rate
          69             :   real(r8), parameter :: pretop = 1.0e2_r8      ! pressure bounding high cloud
          70             :
          71             :   integer count
          72             :
          73             :   logical :: inversion_cld_off    ! Turns off stratification-based cld frc
          74             :
          75             :   integer :: k700   ! model level nearest 700 mb
          76             :
          77             : !================================================================================================
          78             :   contains
          79             : !================================================================================================
          80             :
          81        1536 : subroutine cldfrc_readnl(nlfile)
          82             :
          83             :    use namelist_utils,  only: find_group_name
          84             :    use units,           only: getunit, freeunit
          85             :    use mpishorthand
          86             :
          87             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          88             :
          89             :    ! Local variables
          90             :    integer :: unitn, ierr
          91             :    character(len=*), parameter :: subname = 'cldfrc_readnl'
          92             :
          93             :    namelist /cldfrc_nl/ cldfrc_freeze_dry,      cldfrc_ice,    cldfrc_rhminl, &
          94             :                         cldfrc_rhminl_adj_land, cldfrc_rhminh, cldfrc_sh1,    &
          95             :                         cldfrc_sh2,             cldfrc_dp1,    cldfrc_dp2,    &
          96             :                         cldfrc_premit,          cldfrc_premib, cldfrc_iceopt, &
          97             :                         cldfrc_icecrit
          98             :    !-----------------------------------------------------------------------------
          99             :
         100        1536 :    if (masterproc) then
         101           2 :       unitn = getunit()
         102           2 :       open( unitn, file=trim(nlfile), status='old' )
         103           2 :       call find_group_name(unitn, 'cldfrc_nl', status=ierr)
         104           2 :       if (ierr == 0) then
         105           2 :          read(unitn, cldfrc_nl, iostat=ierr)
         106           2 :          if (ierr /= 0) then
         107           0 :             call endrun(subname // ':: ERROR reading namelist')
         108             :          end if
         109             :       end if
         110           2 :       close(unitn)
         111           2 :       call freeunit(unitn)
         112             :
         113             :       ! set local variables
         114           2 :       rhminl = cldfrc_rhminl
         115           2 :       rhminl_adj_land = cldfrc_rhminl_adj_land
         116           2 :       rhminh = cldfrc_rhminh
         117           2 :       sh1    = cldfrc_sh1
         118           2 :       sh2    = cldfrc_sh2
         119           2 :       dp1    = cldfrc_dp1
         120           2 :       dp2    = cldfrc_dp2
         121           2 :       premit = cldfrc_premit
         122           2 :       premib  = cldfrc_premib
         123           2 :       iceopt  = cldfrc_iceopt
         124           2 :       icecrit = cldfrc_icecrit
         125             :
         126             :    end if
         127             :
         128             : #ifdef SPMD
         129             :    ! Broadcast namelist variables
         130        1536 :    call mpibcast(cldfrc_freeze_dry, 1, mpilog, 0, mpicom)
         131        1536 :    call mpibcast(cldfrc_ice,        1, mpilog, 0, mpicom)
         132        1536 :    call mpibcast(rhminl,            1, mpir8,  0, mpicom)
         133        1536 :    call mpibcast(rhminl_adj_land,   1, mpir8,  0, mpicom)
         134        1536 :    call mpibcast(rhminh,            1, mpir8,  0, mpicom)
         135        1536 :    call mpibcast(sh1   ,            1, mpir8,  0, mpicom)
         136        1536 :    call mpibcast(sh2   ,            1, mpir8,  0, mpicom)
         137        1536 :    call mpibcast(dp1   ,            1, mpir8,  0, mpicom)
         138        1536 :    call mpibcast(dp2   ,            1, mpir8,  0, mpicom)
         139        1536 :    call mpibcast(premit,            1, mpir8,  0, mpicom)
         140        1536 :    call mpibcast(premib,            1, mpir8,  0, mpicom)
         141        1536 :    call mpibcast(iceopt,            1, mpiint, 0, mpicom)
         142        1536 :    call mpibcast(icecrit,           1, mpir8,  0, mpicom)
         143             : #endif
         144             :
         145        1536 : end subroutine cldfrc_readnl
         146             :
         147             : !================================================================================================
         148             :
         149        1536 : subroutine cldfrc_register
         150             :
         151             :    ! Register fields in the physics buffer.
         152             :
         153             :    use physics_buffer, only : pbuf_add_field, dtype_r8
         154             :
         155             :    !-----------------------------------------------------------------------
         156             :
         157        1536 :    call pbuf_add_field('SH_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), sh_frac_idx)
         158        1536 :    call pbuf_add_field('DP_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), dp_frac_idx)
         159             :
         160        1536 : end subroutine cldfrc_register
         161             :
         162             : !================================================================================================
         163             :
         164        1536 : subroutine cldfrc_getparams(rhminl_out, rhminl_adj_land_out, rhminh_out,  premit_out, &
         165             :                             premib_out, iceopt_out,          icecrit_out)
         166             : !-----------------------------------------------------------------------
         167             : ! Purpose: Return cldfrc tuning parameters
         168             : !-----------------------------------------------------------------------
         169             :
         170             :    real(r8),          intent(out), optional :: rhminl_out
         171             :    real(r8),          intent(out), optional :: rhminl_adj_land_out
         172             :    real(r8),          intent(out), optional :: rhminh_out
         173             :    real(r8),          intent(out), optional :: premit_out
         174             :    real(r8),          intent(out), optional :: premib_out
         175             :    integer,           intent(out), optional :: iceopt_out
         176             :    real(r8),          intent(out), optional :: icecrit_out
         177             :
         178        1536 :    if ( present(rhminl_out) )      rhminl_out = rhminl
         179        1536 :    if ( present(rhminl_adj_land_out) ) rhminl_adj_land_out = rhminl_adj_land
         180        1536 :    if ( present(rhminh_out) )      rhminh_out = rhminh
         181        1536 :    if ( present(premit_out) )      premit_out = premit
         182        1536 :    if ( present(premib_out) )      premib_out  = premib
         183        1536 :    if ( present(iceopt_out) )      iceopt_out  = iceopt
         184        1536 :    if ( present(icecrit_out) )     icecrit_out = icecrit
         185             :
         186        1536 : end subroutine cldfrc_getparams
         187             :
         188             : !===============================================================================
         189             :
         190        1536 : subroutine cldfrc_init
         191             :
         192             :    ! Initialize cloud fraction run-time parameters
         193             :
         194             :    use cam_history,   only:  addfld
         195             :    use phys_control,  only:  phys_getopts
         196             :
         197             :    ! query interfaces for scheme settings
         198             :    character(len=16) :: shallow_scheme, eddy_scheme, macrop_scheme
         199             :
         200             :    integer :: k
         201             :    !-----------------------------------------------------------------------------
         202             :
         203             :    call phys_getopts(shallow_scheme_out = shallow_scheme ,&
         204             :                      eddy_scheme_out    = eddy_scheme    ,&
         205        1536 :                      macrop_scheme_out  = macrop_scheme  )
         206             :
         207             :    ! Limit CAM5 cloud physics to below top cloud level.
         208        1536 :    if ( .not. (macrop_scheme == "rk" .or. macrop_scheme == "SPCAM_sam1mom")) top_lev = trop_cloud_top_lev
         209             :
         210             :    ! Turn off inversion_cld if any UW PBL scheme is being used
         211        1536 :    if ( (eddy_scheme .eq. 'diag_TKE' ) .or. (shallow_scheme .eq.  'UW' ) .or.&
         212             :         (shallow_scheme .eq. 'SPCAM_m2005') ) then
         213           0 :       inversion_cld_off = .true.
         214             :    else
         215        1536 :       inversion_cld_off = .false.
         216             :    endif
         217             :
         218        1536 :    if ( masterproc ) then
         219           2 :       write(iulog,*)'tuning parameters cldfrc_init: inversion_cld_off',inversion_cld_off
         220           2 :       write(iulog,*)'tuning parameters cldfrc_init: dp1',dp1,'dp2',dp2,'sh1',sh1,'sh2',sh2
         221           2 :       if (shallow_scheme .ne. 'UW' .or. shallow_scheme .eq. 'SPCAM_m2005' ) then
         222           2 :          write(iulog,*)'tuning parameters cldfrc_init: rhminl',rhminl,'rhminl_adj_land',rhminl_adj_land, &
         223           4 :                        'rhminh',rhminh,'premit',premit,'premib',premib
         224           2 :          write(iulog,*)'tuning parameters cldfrc_init: iceopt',iceopt,'icecrit',icecrit
         225             :       endif
         226             :    endif
         227             :
         228        1536 :    if (pref_mid(top_lev) > 7.e4_r8) &
         229           0 :         call endrun ('cldfrc_init: model levels bracketing 700 mb not found')
         230             :
         231             :    ! Find vertical level nearest 700 mb.
         232      130560 :    k700 = minloc(abs(pref_mid(top_lev:pver) - 7.e4_r8), 1)
         233             :
         234        1536 :    if (masterproc) then
         235           2 :       write(iulog,*)'cldfrc_init: model level nearest 700 mb is',k700,'which is',pref_mid(k700),'pascals'
         236             :    end if
         237             :
         238        3072 :    call addfld ('SH_CLD', (/ 'lev' /), 'A', 'fraction', 'Shallow convective cloud cover' )
         239        3072 :    call addfld ('DP_CLD', (/ 'lev' /), 'A', 'fraction', 'Deep convective cloud cover'    )
         240             :
         241        1536 : end subroutine cldfrc_init
         242             :
         243             : !===============================================================================
         244             :
         245           0 : subroutine cldfrc(lchnk   ,ncol    , pbuf,  &
         246             :        pmid    ,temp    ,q       ,omga    , phis, &
         247             :        shfrc   ,use_shfrc, &
         248             :        cloud   ,rhcloud, clc     ,pdel    , &
         249             :        cmfmc   ,cmfmc2  ,landfrac,snowh   ,concld  ,cldst   , &
         250             :        ts      ,sst     ,ps      ,zdu     ,ocnfrac ,&
         251             :        rhu00   ,cldice  ,icecldf ,liqcldf ,relhum  ,dindex )
         252             :     !-----------------------------------------------------------------------
         253             :     !
         254             :     ! Purpose:
         255             :     ! Compute cloud fraction
         256             :     !
         257             :     !
         258             :     ! Method:
         259             :     ! This calculate cloud fraction using a relative humidity threshold
         260             :     ! The threshold depends upon pressure, and upon the presence or absence
         261             :     ! of convection as defined by a reasonably large vertical mass flux
         262             :     ! entering that layer from below.
         263             :     !
         264             :     ! Author: Many. Last modified by Jim McCaa
         265             :     !
         266             :     !-----------------------------------------------------------------------
         267        1536 :     use cam_history,   only: outfld
         268             :     use physconst,     only: cappa, gravit, rair, tmelt
         269             :     use wv_saturation, only: qsat, qsat_water, svp_ice_vect
         270             :     use phys_grid,     only: get_rlat_all_p, get_rlon_all_p
         271             :
         272             :
         273             : !RBN - Need this to write shallow,deep fraction to phys buffer.
         274             : !PJR - we should probably make seperate modules for determining convective
         275             : !      clouds and make this one just responsible for relative humidity clouds
         276             :
         277             :     use physics_buffer, only: physics_buffer_desc, pbuf_get_field
         278             :
         279             :     ! Arguments
         280             :     integer, intent(in) :: lchnk                  ! chunk identifier
         281             :     integer, intent(in) :: ncol                   ! number of atmospheric columns
         282             :     integer, intent(in) :: dindex                 ! 0 or 1 to perturb rh
         283             :
         284             :     type(physics_buffer_desc), pointer :: pbuf(:)
         285             :     real(r8), intent(in) :: pmid(pcols,pver)      ! midpoint pressures
         286             :     real(r8), intent(in) :: temp(pcols,pver)      ! temperature
         287             :     real(r8), intent(in) :: q(pcols,pver)         ! specific humidity
         288             :     real(r8), intent(in) :: omga(pcols,pver)      ! vertical pressure velocity
         289             :     real(r8), intent(in) :: cmfmc(pcols,pverp)    ! convective mass flux--m sub c
         290             :     real(r8), intent(in) :: cmfmc2(pcols,pverp)   ! shallow convective mass flux--m sub c
         291             :     real(r8), intent(in) :: snowh(pcols)          ! snow depth (liquid water equivalent)
         292             :     real(r8), intent(in) :: pdel(pcols,pver)      ! pressure depth of layer
         293             :     real(r8), intent(in) :: landfrac(pcols)       ! Land fraction
         294             :     real(r8), intent(in) :: ocnfrac(pcols)        ! Ocean fraction
         295             :     real(r8), intent(in) :: ts(pcols)             ! surface temperature
         296             :     real(r8), intent(in) :: sst(pcols)            ! sea surface temperature
         297             :     real(r8), intent(in) :: ps(pcols)             ! surface pressure
         298             :     real(r8), intent(in) :: zdu(pcols,pver)       ! detrainment rate from deep convection
         299             :     real(r8), intent(in) :: phis(pcols)           ! surface geopotential
         300             :     real(r8), intent(in) :: shfrc(pcols,pver)     ! cloud fraction from convect_shallow
         301             :     real(r8), intent(in) :: cldice(pcols,pver)    ! cloud ice mixing ratio
         302             :     logical,  intent(in)  :: use_shfrc
         303             :
         304             :     ! Output arguments
         305             :     real(r8), intent(out) :: cloud(pcols,pver)     ! cloud fraction
         306             :     real(r8), intent(out) :: rhcloud(pcols,pver)   ! cloud fraction
         307             :     real(r8), intent(out) :: clc(pcols)            ! column convective cloud amount
         308             :     real(r8), intent(out) :: cldst(pcols,pver)     ! cloud fraction
         309             :     real(r8), intent(out) :: rhu00(pcols,pver)     ! RH threshold for cloud
         310             :     real(r8), intent(out) :: relhum(pcols,pver)    ! RH
         311             :     real(r8), intent(out) :: icecldf(pcols,pver)   ! ice cloud fraction
         312             :     real(r8), intent(out) :: liqcldf(pcols,pver)   ! liquid cloud fraction (combined into cloud)
         313             :
         314             :     !---------------------------Local workspace-----------------------------
         315             :     !
         316             :     real(r8) concld(pcols,pver)    ! convective cloud cover
         317             :     real(r8) cld                   ! intermediate scratch variable (low cld)
         318             :     real(r8) dthdpmn(pcols)         ! most stable lapse rate below 750 mb
         319             :     real(r8) dthdp                 ! lapse rate (intermediate variable)
         320             :     real(r8) es(pcols,pver)        ! saturation vapor pressure
         321             :     real(r8) qs(pcols,pver)        ! saturation specific humidity
         322             :     real(r8) rhwght                ! weighting function for rhlim transition
         323             :     real(r8) rh(pcols,pver)        ! relative humidity
         324             :     real(r8) rhdif                 ! intermediate scratch variable
         325             :     real(r8) strat                 ! intermediate scratch variable
         326             :     real(r8) theta(pcols,pver)     ! potential temperature
         327             :     real(r8) rhlim                 ! local rel. humidity threshold estimate
         328             :     real(r8) coef1                 ! coefficient to convert mass flux to mb/d
         329             :     real(r8) clrsky(pcols)         ! temporary used in random overlap calc
         330             :     real(r8) rpdeli(pcols,pver-1) ! 1./(pmid(k+1)-pmid(k))
         331             :     real(r8) rhpert                !the specified perturbation to rh
         332             :
         333           0 :     real(r8), pointer, dimension(:,:) :: deepcu      ! deep convection cloud fraction
         334           0 :     real(r8), pointer, dimension(:,:) :: shallowcu   ! shallow convection cloud fraction
         335             :
         336             :     logical cldbnd(pcols)          ! region below high cloud boundary
         337             :
         338             :     integer i, ierror, k           ! column, level indices
         339             :     integer kp1, ifld
         340             :     integer kdthdp(pcols)
         341             :     integer numkcld                ! number of levels in which to allow clouds
         342             :
         343             :     !  In Cloud Ice Content variables
         344             :     real(r8) :: a,b,c,as,bs,cs        !fit parameters
         345             :     real(r8) :: Kc                    !constant for ice cloud calc (wood & field)
         346             :     real(r8) :: ttmp                  !limited temperature
         347             :     real(r8) :: icicval               !empirical iwc value
         348             :     real(r8) :: rho                   !local air density
         349             :     real(r8) :: esl(pcols,pver)       !liq sat vapor pressure
         350             :     real(r8) :: esi(pcols,pver)       !ice sat vapor pressure
         351             :     real(r8) :: ncf,phi               !Wilson and Ballard parameters
         352             :
         353             :     real(r8) thetas(pcols)                    ! ocean surface potential temperature
         354             :     real(r8) :: clat(pcols)                   ! current latitudes(radians)
         355             :     real(r8) :: clon(pcols)                   ! current longitudes(radians)
         356             :
         357             :     ! Statement functions
         358             :     logical land
         359             :     land(i) = nint(landfrac(i)) == 1
         360             :
         361           0 :     call get_rlat_all_p(lchnk, ncol, clat)
         362           0 :     call get_rlon_all_p(lchnk, ncol, clon)
         363             :
         364           0 :     call pbuf_get_field(pbuf, sh_frac_idx, shallowcu )
         365           0 :     call pbuf_get_field(pbuf, dp_frac_idx, deepcu )
         366             :
         367             :     ! Initialise cloud fraction
         368           0 :     shallowcu = 0._r8
         369           0 :     deepcu    = 0._r8
         370             :
         371             :     !==================================================================================
         372             :     ! PHILOSOPHY OF PRESENT IMPLEMENTATION
         373             :     ! Modification to philosophy for ice supersaturation
         374             :     ! philosophy below is based on RH water only. This is 'liquid condensation'
         375             :     ! or liquid cloud (even though it will freeze immediately to ice)
         376             :     ! The idea is that the RH limits for condensation are strict only for
         377             :     ! water saturation
         378             :     !
         379             :     ! Ice clouds are formed by explicit parameterization of ice nucleation.
         380             :     ! Closure for ice cloud fraction is done on available cloud ice, such that
         381             :     ! the in-cloud ice content matches an empirical fit
         382             :     ! thus, icecldf = min(cldice/icicval,1) where icicval = f(temp,cldice,numice)
         383             :     ! for a first cut, icicval=f(temp) only.
         384             :     ! Combined cloud fraction is maximum overlap  cloud=max(1,max(icecldf,liqcldf))
         385             :     ! No dA/dt term for ice?
         386             :     !
         387             :     ! There are three co-existing cloud types: convective, inversion related low-level
         388             :     ! stratocumulus, and layered cloud (based on relative humidity).  Layered and
         389             :     ! stratocumulus clouds do not compete with convective cloud for which one creates
         390             :     ! the most cloud.  They contribute collectively to the total grid-box average cloud
         391             :     ! amount.  This is reflected in the way in which the total cloud amount is evaluated
         392             :     ! (a sum as opposed to a logical "or" operation)
         393             :     !
         394             :     !==================================================================================
         395             :     ! set defaults for rhu00
         396           0 :     rhu00(:,:) = 2.0_r8
         397             :     ! define rh perturbation in order to estimate rhdfda
         398           0 :     rhpert = 0.01_r8
         399             :
         400             :     !set Wang and Sassen IWC paramters
         401           0 :     a=26.87_r8
         402           0 :     b=0.569_r8
         403           0 :     c=0.002892_r8
         404             :     !set schiller parameters
         405           0 :     as=-68.4202_r8
         406           0 :     bs=0.983917_r8
         407           0 :     cs=2.81795_r8
         408             :     !set wood and field paramters...
         409           0 :     Kc=75._r8
         410             :
         411             :     ! Evaluate potential temperature and relative humidity
         412             :     ! If not computing ice cloud fraction then hybrid RH, if MG then water RH
         413           0 :     if ( cldfrc_ice ) then
         414           0 :        do k = top_lev,pver
         415           0 :           call qsat_water(temp(1:ncol,k), pmid(1:ncol,k), esl(1:ncol,k), qs(1:ncol,k), ncol)
         416           0 :           call svp_ice_vect(temp(1:ncol,k), esi(1:ncol,k), ncol)
         417             :        end do
         418             :     else
         419           0 :        do k = top_lev,pver
         420           0 :           call qsat(temp(1:ncol,k), pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol)
         421             :        end do
         422             :     end if
         423             :
         424           0 :     cloud    = 0._r8
         425           0 :     icecldf  = 0._r8
         426           0 :     liqcldf  = 0._r8
         427           0 :     rhcloud  = 0._r8
         428           0 :     cldst    = 0._r8
         429           0 :     concld   = 0._r8
         430             :
         431           0 :     do k=top_lev,pver
         432           0 :        theta(:ncol,k)    = temp(:ncol,k)*(pnot/pmid(:ncol,k))**cappa
         433             :
         434           0 :        do i=1,ncol
         435           0 :           rh(i,k)     = q(i,k)/qs(i,k)*(1.0_r8+real(dindex,r8)*rhpert)
         436             :           !  record relhum, rh itself will later be modified related with concld
         437           0 :           relhum(i,k) = rh(i,k)
         438             :        end do
         439             :     end do
         440             :
         441             :     ! Initialize other temporary variables
         442           0 :     ierror = 0
         443           0 :     do i=1,ncol
         444             :        ! Adjust thetas(i) in the presence of non-zero ocean heights.
         445             :        ! This reduces the temperature for positive heights according to a standard lapse rate.
         446           0 :        if(ocnfrac(i).gt.0.01_r8) thetas(i)  = &
         447           0 :             ( sst(i) - lapse * phis(i) / gravit) * (pnot/ps(i))**cappa
         448           0 :        if(ocnfrac(i).gt.0.01_r8.and.sst(i).lt.260._r8) ierror = i
         449           0 :        clc(i) = 0.0_r8
         450             :     end do
         451           0 :     coef1 = gravit*864.0_r8    ! conversion to millibars/day
         452             :
         453           0 :     if (ierror > 0) then
         454           0 :        write(iulog,*) 'COLDSST: encountered in cldfrc:', lchnk,ierror,ocnfrac(ierror),sst(ierror)
         455             :     endif
         456             :
         457           0 :     do k=top_lev,pver-1
         458           0 :        rpdeli(:ncol,k) = 1._r8/(pmid(:ncol,k+1) - pmid(:ncol,k))
         459             :     end do
         460             :
         461             :     !
         462             :     ! Estimate of local convective cloud cover based on convective mass flux
         463             :     ! Modify local large-scale relative humidity to account for presence of
         464             :     ! convective cloud when evaluating relative humidity based layered cloud amount
         465             :     !
         466           0 :     concld(:ncol,top_lev:pver) = 0.0_r8
         467             :     !
         468             :     ! cloud mass flux in SI units of kg/m2/s; should produce typical numbers of 20%
         469             :     ! shallow and deep convective cloudiness are evaluated separately (since processes
         470             :     ! are evaluated separately) and summed
         471             :     !
         472             : #ifndef PERGRO
         473           0 :     do k=top_lev,pver
         474           0 :        do i=1,ncol
         475           0 :           if ( .not. use_shfrc ) then
         476           0 :              shallowcu(i,k) = max(0.0_r8,min(sh1*log(1.0_r8+sh2*cmfmc2(i,k+1)),0.30_r8))
         477             :           else
         478           0 :              shallowcu(i,k) = shfrc(i,k)
         479             :           endif
         480           0 :           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))
         481           0 :           concld(i,k) = min(shallowcu(i,k) + deepcu(i,k),0.80_r8)
         482           0 :           rh(i,k) = (rh(i,k) - concld(i,k))/(1.0_r8 - concld(i,k))
         483             :        end do
         484             :     end do
         485             : #endif
         486             :     !==================================================================================
         487             :     !
         488             :     !          ****** Compute layer cloudiness ******
         489             :     !
         490             :     !====================================================================
         491             :     ! Begin the evaluation of layered cloud amount based on (modified) RH
         492             :     !====================================================================
         493             :     !
         494           0 :     numkcld = pver
         495           0 :     do k=top_lev+1,numkcld
         496           0 :        kp1 = min(k + 1,pver)
         497           0 :        do i=1,ncol
         498             :
         499             :           ! This is now designed to apply FOR LIQUID CLOUDS (condensation > RH water)
         500             :
         501           0 :           cldbnd(i) = pmid(i,k).ge.pretop
         502             :
         503           0 :           if ( pmid(i,k).ge.premib ) then
         504             :              !==============================================================
         505             :              ! This is the low cloud (below premib) block
         506             :              !==============================================================
         507             :              ! enhance low cloud activation over land with no snow cover
         508           0 :              if (land(i) .and. (snowh(i) <= 0.000001_r8)) then
         509           0 :                 rhlim = rhminl - rhminl_adj_land
         510             :              else
         511           0 :                 rhlim = rhminl
         512             :              endif
         513             :
         514           0 :              rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
         515           0 :              rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
         516             :
         517             :              ! SJV: decrease cloud amount if very low water vapor content
         518             :              ! (thus very cold): "freeze dry"
         519           0 :              if (cldfrc_freeze_dry) then
         520           0 :                 rhcloud(i,k) = rhcloud(i,k)*max(0.15_r8,min(1.0_r8,q(i,k)/0.0030_r8))
         521             :              endif
         522             :
         523           0 :           else if ( pmid(i,k).lt.premit ) then
         524             :              !==============================================================
         525             :              ! This is the high cloud (above premit) block
         526             :              !==============================================================
         527             :              !
         528           0 :              rhlim = rhminh
         529             :              !
         530           0 :              rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
         531           0 :              rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
         532             :           else
         533             :              !==============================================================
         534             :              ! This is the middle cloud block
         535             :              !==============================================================
         536             :              !
         537             :              !       linear rh threshold transition between thresholds for low & high cloud
         538             :              !
         539           0 :              rhwght = (premib-(max(pmid(i,k),premit)))/(premib-premit)
         540             :
         541           0 :              if (land(i) .and. (snowh(i) <= 0.000001_r8)) then
         542           0 :                 rhlim = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
         543             :              else
         544           0 :                 rhlim = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
         545             :              endif
         546           0 :              rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
         547           0 :              rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
         548             :           end if
         549             :           !==================================================================================
         550             :           ! WE NEED TO DOCUMENT THE PURPOSE OF THIS TYPE OF CODE (ASSOCIATED WITH 2ND CALL)
         551             :           !==================================================================================
         552             :           !      !
         553             :           !      ! save rhlim to rhu00, it handles well by itself for low/high cloud
         554             :           !      !
         555           0 :           rhu00(i,k)=rhlim
         556             :           !==================================================================================
         557             :
         558           0 :           if (cldfrc_ice) then
         559             :
         560             :              ! Evaluate ice cloud fraction based on in-cloud ice content
         561             :
         562             :              !--------ICE CLOUD OPTION 1--------Wang & Sassen 2002
         563             :              !         Evaluate desired in-cloud water content
         564             :              !               icicval = f(temp,cldice,numice)
         565             :              !         Start with a function of temperature.
         566             :              !         Wang & Sassen 2002 (JAS), based on ARM site MMCR (midlat cirrus)
         567             :              !           parameterization valid for 203-253K
         568             :              !           icival > 0 for t>195K
         569           0 :              if (iceopt.lt.3) then
         570           0 :                 if (iceopt.eq.1) then
         571           0 :                    ttmp=max(195._r8,min(temp(i,k),253._r8)) - 273.16_r8
         572           0 :                    icicval=a + b * ttmp + c * ttmp**2._r8
         573             :                    !convert units
         574           0 :                    rho=pmid(i,k)/(rair*temp(i,k))
         575           0 :                    icicval= icicval * 1.e-6_r8 / rho
         576             :                 else
         577             :                    !--------ICE CLOUD OPTION 2--------Schiller 2008 (JGR)
         578             :                    !          Use a curve based on FISH measurements in
         579             :                    !          tropics, mid-lats and arctic. Curve is for 180-250K (raise to 273K?)
         580             :                    !          use median all flights
         581             :
         582           0 :                    ttmp=max(190._r8,min(temp(i,k),273.16_r8))
         583           0 :                    icicval = 10._r8 **(as * bs**ttmp + cs)
         584             :                    !convert units from ppmv to kg/kg
         585           0 :                    icicval= icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
         586             :                 endif
         587             :                 !set icecldfraction  for OPTION 1 or OPTION2
         588           0 :                 icecldf(i,k) =  max(0._r8,min(cldice(i,k)/icicval,1._r8))
         589             :
         590           0 :              else if (iceopt.eq.3) then
         591             :
         592             :                 !--------ICE CLOUD OPTION 3--------Wood & Field 2000 (JAS)
         593             :                 ! eq 6: cloud fraction = 1 - exp (-K * qc/qsati)
         594             :
         595           0 :                 icecldf(i,k)=1._r8 - exp(-Kc*cldice(i,k)/(qs(i,k)*(esi(i,k)/esl(i,k))))
         596           0 :                 icecldf(i,k)=max(0._r8,min(icecldf(i,k),1._r8))
         597             :              else
         598             :                 !--------ICE CLOUD OPTION 4--------Wilson and ballard 1999
         599             :                 ! inversion of smith....
         600             :                 !       ncf = cldice / ((1-RHcrit)*qs)
         601             :                 ! then a function of ncf....
         602           0 :                 ncf =cldice(i,k)/((1._r8 - icecrit)*qs(i,k))
         603           0 :                 if (ncf.le.0._r8) then
         604           0 :                    icecldf(i,k)=0._r8
         605           0 :                 else if (ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8) then
         606           0 :                    icecldf(i,k)=0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
         607           0 :                 else if (ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8) then
         608           0 :                    phi=(acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
         609           0 :                    icecldf(i,k)=(1._r8 - 4._r8 * cos(phi) * cos(phi))
         610             :                 else
         611           0 :                    icecldf(i,k)=1._r8
         612             :                 endif
         613           0 :                 icecldf(i,k)=max(0._r8,min(icecldf(i,k),1._r8))
         614             :              endif
         615             :              !TEST: if ice present, icecldf=1.
         616             :              !          if (cldice(i,k).ge.1.e-8_r8) then
         617             :              !             icecldf(i,k) = 0.99_r8
         618             :              !          endif
         619             :
         620             :              !!          if ((cldice(i,k) .gt. icicval) .or. ((cldice(i,k) .gt. 0._r8) .and. (icecldf(i,k) .eq. 0._r8))) then
         621             :              !          if (cldice(i,k) .gt. 1.e-8_r8) then
         622             :              !             write(iulog,*) 'i,k,pmid,rho,t,cldice,icicval,icecldf,rhcloud: ', &
         623             :              !                i,k,pmid(i,k),rho,temp(i,k),cldice(i,k),icicval,icecldf(i,k),rhcloud(i,k)
         624             :              !          endif
         625             :
         626             :              !         Combine ice and liquid cloud fraction assuming maximum overlap.
         627             :              ! Combined cloud fraction is maximum overlap
         628             :              !          cloud(i,k)=min(1._r8,max(icecldf(i,k),rhcloud(i,k)))
         629             :
         630           0 :              liqcldf(i,k)=(1._r8 - icecldf(i,k))* rhcloud(i,k)
         631           0 :              cloud(i,k)=liqcldf(i,k) + icecldf(i,k)
         632             :           else
         633             :              ! For RK microphysics
         634           0 :              cloud(i,k) = rhcloud(i,k)
         635             :           end if
         636             :        end do
         637             :     end do
         638             :     !
         639             :     ! Add in the marine strat
         640             :     ! MARINE STRATUS SHOULD BE A SPECIAL CASE OF LAYERED CLOUD
         641             :     ! CLOUD CURRENTLY CONTAINS LAYERED CLOUD DETERMINED BY RH CRITERIA
         642             :     ! TAKE THE MAXIMUM OF THE DIAGNOSED LAYERED CLOUD OR STRATOCUMULUS
         643             :     !
         644             :     !===================================================================================
         645             :     !
         646             :     !  SOME OBSERVATIONS ABOUT THE FOLLOWING SECTION OF CODE (missed in earlier look)
         647             :     !  K700 IS SET AS A CONSTANT BASED ON HYBRID COORDINATE: IT DOES NOT DEPEND ON
         648             :     !  LOCAL PRESSURE; THERE IS NO PRESSURE RAMP => LOOKS LEVEL DEPENDENT AND
         649             :     !  DISCONTINUOUS IN SPACE (I.E., STRATUS WILL END SUDDENLY WITH NO TRANSITION)
         650             :     !
         651             :     !  IT APPEARS THAT STRAT IS EVALUATED ACCORDING TO KLEIN AND HARTMANN; HOWEVER,
         652             :     !  THE ACTUAL STRATUS AMOUNT (CLDST) APPEARS TO DEPEND DIRECTLY ON THE RH BELOW
         653             :     !  THE STRONGEST PART OF THE LOW LEVEL INVERSION.
         654             :     !PJR answers: 1) the rh limitation is a physical/mathematical limitation
         655             :     !             cant have more cloud than there is RH
         656             :     !             allowed the cloud to exist two layers below the inversion
         657             :     !             because the numerics frequently make 50% relative humidity
         658             :     !             in level below the inversion which would allow no cloud
         659             :     !             2) since  the cloud is only allowed over ocean, it should
         660             :     !             be very insensitive to surface pressure (except due to
         661             :     !             spectral ringing, which also causes so many other problems
         662             :     !             I didnt worry about it.
         663             :     !
         664             :     !==================================================================================
         665           0 :     if (.not.inversion_cld_off) then
         666             :     !
         667             :     ! Find most stable level below 750 mb for evaluating stratus regimes
         668             :     !
         669           0 :     do i=1,ncol
         670             :        ! Nothing triggers unless a stability greater than this minimum threshold is found
         671           0 :        dthdpmn(i) = -0.125_r8
         672           0 :        kdthdp(i) = 0
         673             :     end do
         674             :     !
         675           0 :     do k=top_lev+1,pver
         676           0 :        do i=1,ncol
         677           0 :           if (pmid(i,k) >= premib .and. ocnfrac(i).gt. 0.01_r8) then
         678             :              ! I think this is done so that dtheta/dp is in units of dg/mb (JJH)
         679           0 :              dthdp = 100.0_r8*(theta(i,k) - theta(i,k-1))*rpdeli(i,k-1)
         680           0 :              if (dthdp < dthdpmn(i)) then
         681           0 :                 dthdpmn(i) = dthdp
         682           0 :                 kdthdp(i) = k     ! index of interface of max inversion
         683             :              end if
         684             :           end if
         685             :        end do
         686             :     end do
         687             :
         688             :     ! Also check between the bottom layer and the surface
         689             :     ! Only perform this check if the criteria were not met above
         690             :
         691           0 :     do i = 1,ncol
         692           0 :        if ( kdthdp(i) .eq. 0 .and. ocnfrac(i).gt.0.01_r8) then
         693           0 :           dthdp = 100.0_r8 * (thetas(i) - theta(i,pver)) / (ps(i)-pmid(i,pver))
         694           0 :           if (dthdp < dthdpmn(i)) then
         695           0 :              dthdpmn(i) = dthdp
         696           0 :              kdthdp(i) = pver     ! index of interface of max inversion
         697             :           endif
         698             :        endif
         699             :     enddo
         700             :
         701           0 :     do i=1,ncol
         702           0 :        if (kdthdp(i) /= 0) then
         703           0 :           k = kdthdp(i)
         704           0 :           kp1 = min(k+1,pver)
         705             :           ! Note: strat will be zero unless ocnfrac > 0.01
         706           0 :           strat = min(1._r8,max(0._r8, ocnfrac(i) * ((theta(i,k700)-thetas(i))*.057_r8-.5573_r8) ) )
         707             :           !
         708             :           ! assign the stratus to the layer just below max inversion
         709             :           ! the relative humidity changes so rapidly across the inversion
         710             :           ! that it is not safe to just look immediately below the inversion
         711             :           ! so limit the stratus cloud by rh in both layers below the inversion
         712             :           !
         713           0 :           cldst(i,k) = min(strat,max(rh(i,k),rh(i,kp1)))
         714             :        end if
         715             :     end do
         716             :     end if  ! .not.inversion_cld_off
         717             :
         718           0 :     do k=top_lev,pver
         719           0 :        do i=1,ncol
         720             :           !
         721             :           !       which is greater; standard layered cloud amount or stratocumulus diagnosis
         722             :           !
         723           0 :           cloud(i,k) = max(rhcloud(i,k),cldst(i,k))
         724             :           !
         725             :           !       add in the contributions of convective cloud (determined separately and accounted
         726             :           !       for by modifications to the large-scale relative humidity.
         727             :           !
         728           0 :           cloud(i,k) = min(cloud(i,k)+concld(i,k), 1.0_r8)
         729             :        end do
         730             :     end do
         731             :
         732           0 :     call outfld( 'SH_CLD  ', shallowcu   , pcols, lchnk )
         733           0 :     call outfld( 'DP_CLD  ', deepcu      , pcols, lchnk )
         734             :
         735             :     !
         736           0 :     return
         737           0 :   end subroutine cldfrc
         738             :
         739             : !================================================================================================
         740             :
         741     1495368 :   subroutine cldfrc_fice(ncol, t, fice, fsnow)
         742             : !
         743             : ! Compute the fraction of the total cloud water which is in ice phase.
         744             : ! The fraction depends on temperature only.
         745             : ! This is the form that was used for radiation, the code came from cldefr originally
         746             : !
         747             : ! Author: B. A. Boville Sept 10, 2002
         748             : !  modified: PJR 3/13/03 (added fsnow to ascribe snow production for convection )
         749             : !-----------------------------------------------------------------------
         750           0 :     use physconst, only: tmelt
         751             :
         752             : ! Arguments
         753             :     integer,  intent(in)  :: ncol                 ! number of active columns
         754             :     real(r8), intent(in)  :: t(:,:)        ! temperature
         755             :
         756             :     real(r8), intent(out) :: fice(:,:)     ! Fractional ice content within cloud
         757             :     real(r8), intent(out) :: fsnow(:,:)    ! Fractional snow content for convection
         758             :
         759             : ! Local variables
         760             :     real(r8) :: tmax_fice                         ! max temperature for cloud ice formation
         761             :     real(r8) :: tmin_fice                         ! min temperature for cloud ice formation
         762             :     real(r8) :: tmax_fsnow                        ! max temperature for transition to convective snow
         763             :     real(r8) :: tmin_fsnow                        ! min temperature for transition to convective snow
         764             :
         765             :     integer :: i,k                                ! loop indexes
         766             :
         767             : !-----------------------------------------------------------------------
         768             :
         769     1495368 :     tmax_fice = tmelt - 10._r8        ! max temperature for cloud ice formation
         770     1495368 :     tmin_fice = tmax_fice - 30._r8    ! min temperature for cloud ice formation
         771     1495368 :     tmax_fsnow = tmelt                ! max temperature for transition to convective snow
         772     1495368 :     tmin_fsnow = tmelt - 5._r8        ! min temperature for transition to convective snow
         773             :
         774   227713248 :     fice(:,:top_lev-1) = 0._r8
         775   227713248 :     fsnow(:,:top_lev-1) = 0._r8
         776             :
         777             : ! Define fractional amount of cloud that is ice
         778   127106280 :     do k=top_lev,pver
         779  2098905480 :        do i=1,ncol
         780             :
         781             : ! If warmer than tmax then water phase
         782  1971799200 :           if (t(i,k) > tmax_fice) then
         783   502589268 :              fice(i,k) = 0.0_r8
         784             :
         785             : ! If colder than tmin then ice phase
         786  1469209932 :           else if (t(i,k) < tmin_fice) then
         787  1082589970 :              fice(i,k) = 1.0_r8
         788             :
         789             : ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
         790             :           else
         791   386619962 :              fice(i,k) =(tmax_fice - t(i,k)) / (tmax_fice - tmin_fice)
         792             :           end if
         793             :
         794             : ! snow fraction partitioning
         795             :
         796             : ! If warmer than tmax then water phase
         797  2097410112 :           if (t(i,k) > tmax_fsnow) then
         798   375834569 :              fsnow(i,k) = 0.0_r8
         799             :
         800             : ! If colder than tmin then ice phase
         801  1595964631 :           else if (t(i,k) < tmin_fsnow) then
         802  1531620515 :              fsnow(i,k) = 1.0_r8
         803             :
         804             : ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
         805             :           else
         806    64344116 :              fsnow(i,k) =(tmax_fsnow - t(i,k)) / (tmax_fsnow - tmin_fsnow)
         807             :           end if
         808             :
         809             :        end do
         810             :     end do
         811             :
         812     1495368 :   end subroutine cldfrc_fice
         813             :
         814             : end module cloud_fraction

Generated by: LCOV version 1.14