LCOV - code coverage report
Current view: | top level - physics/cam - cldfrc2m.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 44 | 401 | 11.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 2 | 8 | 25.0 %

              Line data    Source code
           1             : module cldfrc2m
           2             :
           3             : ! cloud fraction calculations
           4             :
           5             : use shr_kind_mod,     only: r8=>shr_kind_r8
           6             : use spmd_utils,       only: masterproc
           7             : use ppgrid,           only: pcols
           8             : use physconst,        only: rair
           9             : use wv_saturation,    only: qsat_water, svp_water, svp_ice, &
          10             :                             svp_water_vect, svp_ice_vect
          11             : use cam_logfile,      only: iulog
          12             : use cam_abortutils,   only: endrun
          13             :
          14             : implicit none
          15             : private
          16             : save
          17             :
          18             : public ::           &
          19             :    cldfrc2m_readnl, &
          20             :    cldfrc2m_init,   &
          21             :    astG_PDF_single, &
          22             :    astG_PDF,        &
          23             :    astG_RHU_single, &
          24             :    astG_RHU,        &
          25             :    aist_single,     &
          26             :    aist_vector,     &
          27             :    CAMstfrac,       &
          28             :    rhmini_const,    &
          29             :    rhmaxi_const,    &
          30             :    rhminis_const,   &
          31             :    rhmaxis_const
          32             :
          33             : ! Namelist variables
          34             : real(r8) :: cldfrc2m_rhmini            ! Minimum rh for ice cloud fraction > 0.
          35             : real(r8) :: cldfrc2m_rhmaxi
          36             : real(r8) :: cldfrc2m_rhminis           ! Minimum rh for ice cloud fraction > 0 in the stratsophere.
          37             : real(r8) :: cldfrc2m_rhmaxis
          38             : real(r8) :: cldfrc2m_qist_min          ! Minimum in-stratus IWC constraint [ kg/kg ]
          39             : real(r8) :: cldfrc2m_qist_max          ! Maximum in-stratus IWC constraint [ kg/kg ]
          40             : logical  :: cldfrc2m_do_subgrid_growth = .false.
          41             : logical  :: cldfrc2m_do_avg_aist_algs = .false.
          42             : ! -------------------------- !
          43             : ! Parameters for Ice Stratus !
          44             : ! -------------------------- !
          45             : real(r8),  protected :: rhmini_const                 ! Minimum rh for ice cloud fraction > 0.
          46             : real(r8),  protected :: rhmaxi_const
          47             : real(r8),  protected :: rhminis_const                ! Minimum rh for ice cloud fraction > 0.
          48             : real(r8),  protected :: rhmaxis_const
          49             :
          50             : ! ----------------------------- !
          51             : ! Parameters for Liquid Stratus !
          52             : ! ----------------------------- !
          53             :
          54             : logical,  parameter  :: CAMstfrac    = .false.    ! If .true. (.false.),
          55             :                                                   ! use Slingo (triangular PDF-based) liquid stratus fraction
          56             : logical,  parameter  :: freeze_dry   = .false.    ! If .true., use 'freeze dry' in liquid stratus fraction formula
          57             : real(r8)             :: rhminl_const              ! Critical RH for low-level  liquid stratus clouds
          58             : real(r8)             :: rhminl_adj_land_const     ! rhminl adjustment for snowfree land
          59             : real(r8)             :: rhminh_const              ! Critical RH for high-level liquid stratus clouds
          60             : real(r8)             :: premit                    ! Top    height for mid-level liquid stratus fraction
          61             : real(r8)             :: premib                    ! Bottom height for mid-level liquid stratus fraction
          62             : integer              :: iceopt                    ! option for ice cloud closure
          63             :                                                   ! 1=wang & sassen 2=schiller (iciwc)
          64             :                                                   ! 3=wood & field, 4=Wilson (based on smith)
          65             :                                                   ! 5=modified slingo (ssat & empyt cloud)
          66             : real(r8)             :: icecrit                   ! Critical RH for ice clouds in Wilson & Ballard closure
          67             :                                                   ! ( smaller = more ice clouds )
          68             :
          69             : !================================================================================================
          70             : contains
          71             : !================================================================================================
          72             :
          73        1536 : subroutine cldfrc2m_readnl(nlfile)
          74             :
          75             :    use namelist_utils,  only: find_group_name
          76             :    use units,           only: getunit, freeunit
          77             :    use spmd_utils,      only: mpicom, masterprocid, mpi_logical, mpi_real8
          78             :
          79             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          80             :
          81             :    ! Local variables
          82             :    integer :: unitn, ierr
          83             :    character(len=*), parameter :: subname = 'cldfrc2m_readnl'
          84             :
          85             :    namelist /cldfrc2m_nl/ cldfrc2m_rhmini, cldfrc2m_rhmaxi, cldfrc2m_rhminis, cldfrc2m_rhmaxis, cldfrc2m_do_subgrid_growth, &
          86             :                           cldfrc2m_qist_min, cldfrc2m_qist_max, cldfrc2m_do_avg_aist_algs
          87             :    !-----------------------------------------------------------------------------
          88             :
          89        1536 :    if (masterproc) then
          90           2 :       unitn = getunit()
          91           2 :       open( unitn, file=trim(nlfile), status='old' )
          92           2 :       call find_group_name(unitn, 'cldfrc2m_nl', status=ierr)
          93           2 :       if (ierr == 0) then
          94           2 :          read(unitn, cldfrc2m_nl, iostat=ierr)
          95           2 :          if (ierr /= 0) then
          96           0 :             call endrun(subname // ':: ERROR reading namelist')
          97             :          end if
          98             :       end if
          99           2 :       close(unitn)
         100           2 :       call freeunit(unitn)
         101             :
         102             :       ! set local variables
         103           2 :       rhmini_const  = cldfrc2m_rhmini
         104           2 :       rhmaxi_const  = cldfrc2m_rhmaxi
         105           2 :       rhminis_const = cldfrc2m_rhminis
         106           2 :       rhmaxis_const = cldfrc2m_rhmaxis
         107             :    end if
         108             :
         109             :    ! Broadcast namelist variables
         110        1536 :    call mpi_bcast(rhmini_const,               1, mpi_real8,  masterprocid, mpicom, ierr)
         111        1536 :    call mpi_bcast(rhmaxi_const,               1, mpi_real8,  masterprocid, mpicom, ierr)
         112        1536 :    call mpi_bcast(rhminis_const,              1, mpi_real8,  masterprocid, mpicom, ierr)
         113        1536 :    call mpi_bcast(rhmaxis_const,              1, mpi_real8,  masterprocid, mpicom, ierr)
         114        1536 :    call mpi_bcast(cldfrc2m_qist_min,          1, mpi_real8,  masterprocid, mpicom, ierr)
         115        1536 :    call mpi_bcast(cldfrc2m_qist_max,          1, mpi_real8,  masterprocid, mpicom, ierr)
         116        1536 :    call mpi_bcast(cldfrc2m_do_subgrid_growth, 1, mpi_logical,masterprocid, mpicom, ierr)
         117        1536 :    call mpi_bcast(cldfrc2m_do_avg_aist_algs,  1, mpi_logical,masterprocid, mpicom, ierr)
         118             :
         119        1536 : end subroutine cldfrc2m_readnl
         120             :
         121             : !================================================================================================
         122             :
         123        1536 : subroutine cldfrc2m_init()
         124             :
         125             :    use cloud_fraction, only: cldfrc_getparams
         126             :
         127             :    call cldfrc_getparams(rhminl_out=rhminl_const, rhminl_adj_land_out=rhminl_adj_land_const,  &
         128             :                          rhminh_out=rhminh_const, premit_out=premit, premib_out=premib, &
         129        1536 :                          iceopt_out=iceopt, icecrit_out=icecrit)
         130             :
         131        1536 :    if( masterproc ) then
         132           2 :       write(iulog,*) 'cldfrc2m parameters:'
         133           2 :       write(iulog,*) '  rhminl            = ', rhminl_const
         134           2 :       write(iulog,*) '  rhminl_adj_land   = ', rhminl_adj_land_const
         135           2 :       write(iulog,*) '  rhminh            = ', rhminh_const
         136           2 :       write(iulog,*) '  premit            = ', premit
         137           2 :       write(iulog,*) '  premib            = ', premib
         138           2 :       write(iulog,*) '  iceopt            = ', iceopt
         139           2 :       write(iulog,*) '  icecrit           = ', icecrit
         140           2 :       write(iulog,*) '  rhmini            = ', rhmini_const
         141           2 :       write(iulog,*) '  rhmaxi            = ', rhmaxi_const
         142           2 :       write(iulog,*) '  rhminis           = ', rhminis_const
         143           2 :       write(iulog,*) '  rhmaxis           = ', rhmaxis_const
         144           2 :       write(iulog,*) '  do_subgrid_growth = ', cldfrc2m_do_subgrid_growth
         145           2 :       write(iulog,*) '  do_avg_aist_algs  = ', cldfrc2m_do_avg_aist_algs
         146           2 :       write(iulog,*) '  cldfrc2m_qist_min = ', cldfrc2m_qist_min
         147           2 :       write(iulog,*) '  cldfrc2m_qist_max = ', cldfrc2m_qist_max
         148             :    end if
         149             :
         150        1536 : end subroutine cldfrc2m_init
         151             :
         152             : !================================================================================================
         153             :
         154             :
         155           0 : subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, orhmin, &
         156             :                            rhminl_in, rhminl_adj_land_in, rhminh_in )
         157             :
         158             :    ! --------------------------------------------------------- !
         159             :    ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
         160             :    ! analytical formulation of triangular PDF.                 !
         161             :    ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
         162             :    ! so using constant 'dV' assume that width is proportional  !
         163             :    ! to the saturation specific humidity.                      !
         164             :    !    dV ~ 0.1.                                              !
         165             :    !    cldrh : RH of in-stratus( = 1 if no supersaturation)   !
         166             :    ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
         167             :    ! G is discontinuous across U = 1.  In fact, it does not    !
         168             :    ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that   !
         169             :    ! they will produce the same results.                       !
         170             :    ! --------------------------------------------------------- !
         171             :
         172             :    real(r8), intent(in)  :: U                     ! Relative humidity
         173             :    real(r8), intent(in)  :: p                     ! Pressure [Pa]
         174             :    real(r8), intent(in)  :: qv                    ! Grid-mean water vapor specific humidity [kg/kg]
         175             :    real(r8), intent(in)  :: landfrac              ! Land fraction
         176             :    real(r8), intent(in)  :: snowh                 ! Snow depth (liquid water equivalent)
         177             :
         178             :    real(r8), intent(out) :: a                     ! Stratus fraction
         179             :    real(r8), intent(out) :: Ga                    ! dU/da
         180             :    real(r8), optional, intent(out) :: orhmin      ! Critical RH
         181             :
         182             :    real(r8), optional, intent(in)  :: rhminl_in          ! Critical relative humidity for low-level  liquid stratus
         183             :    real(r8), optional, intent(in)  :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land
         184             :    real(r8), optional, intent(in)  :: rhminh_in          ! Critical relative humidity for high-level liquid stratus
         185             :
         186             :    ! Local variables
         187             :    integer :: i                                   ! Loop indexes
         188             :    real(r8) dV                                    ! Width of triangular PDF
         189             :    real(r8) cldrh                                 ! RH of stratus cloud
         190             :    real(r8) rhmin                                 ! Critical RH
         191             :    real(r8) rhwght
         192             :
         193             :    real(r8) :: rhminl
         194             :    real(r8) :: rhminl_adj_land
         195             :    real(r8) :: rhminh
         196             :
         197             :    ! Statement functions
         198             :    logical land
         199           0 :    land = nint(landfrac) == 1
         200             :
         201             :    ! ---------- !
         202             :    ! Parameters !
         203             :    ! ---------- !
         204             :
         205           0 :    cldrh  = 1.0_r8
         206             :
         207           0 :    rhminl = rhminl_const
         208           0 :    if (present(rhminl_in)) rhminl = rhminl_in
         209           0 :    rhminl_adj_land = rhminl_adj_land_const
         210           0 :    if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in
         211           0 :    rhminh = rhminh_const
         212           0 :    if (present(rhminh_in)) rhminh = rhminh_in
         213             :
         214             :    ! ---------------- !
         215             :    ! Main computation !
         216             :    ! ---------------- !
         217             :
         218           0 :    if( p .ge. premib ) then
         219             :
         220           0 :        if( land .and. (snowh.le.0.000001_r8) ) then
         221           0 :            rhmin = rhminl - rhminl_adj_land
         222             :        else
         223             :            rhmin = rhminl
         224             :        endif
         225             :
         226           0 :        dV = cldrh - rhmin
         227             :
         228           0 :        if( U .ge. 1._r8 ) then
         229           0 :            a  = 1._r8
         230           0 :            Ga = 1.e10_r8
         231           0 :        elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
         232           0 :            a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
         233           0 :            Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
         234           0 :        elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
         235             :            a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
         236           0 :                       (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
         237           0 :            Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
         238           0 :        elseif( U .le. (cldrh-dV) ) then
         239           0 :            a  = 0._r8
         240           0 :            Ga = 1.e10_r8
         241             :        endif
         242             :
         243             :        if( freeze_dry ) then
         244             :            a  = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         245             :            Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         246             :        endif
         247             :
         248           0 :    elseif( p .lt. premit ) then
         249             :
         250           0 :        rhmin = rhminh
         251           0 :        dV    = cldrh - rhmin
         252             :
         253           0 :        if( U .ge. 1._r8 ) then
         254           0 :            a  = 1._r8
         255           0 :            Ga = 1.e10_r8
         256           0 :        elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
         257           0 :            a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
         258           0 :            Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
         259           0 :        elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
         260             :            a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
         261           0 :                       (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
         262           0 :            Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
         263           0 :        elseif( U .le. (cldrh-dV) ) then
         264           0 :            a  = 0._r8
         265           0 :            Ga = 1.e10_r8
         266             :        endif
         267             :
         268             :    else
         269             :
         270           0 :        rhwght = (premib-(max(p,premit)))/(premib-premit)
         271             :
         272             :      ! if( land .and. (snowh.le.0.000001_r8) ) then
         273             :      !     rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
         274             :      ! else
         275           0 :            rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
         276             :      ! endif
         277             :
         278           0 :        dV    = cldrh - rhmin
         279             :
         280           0 :        if( U .ge. 1._r8 ) then
         281           0 :            a  = 1._r8
         282           0 :            Ga = 1.e10_r8
         283           0 :        elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
         284           0 :            a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
         285           0 :            Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
         286           0 :        elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
         287             :            a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
         288           0 :                          (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
         289           0 :            Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
         290           0 :        elseif( U .le. (cldrh-dV) ) then
         291           0 :            a  = 0._r8
         292           0 :            Ga = 1.e10_r8
         293             :        endif
         294             :
         295             :    endif
         296             :
         297           0 :    if (present(orhmin)) orhmin = rhmin
         298             :
         299        1536 : end subroutine astG_PDF_single
         300             :
         301             : !================================================================================================
         302             :
         303           0 : subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, &
         304             :                     rhminl_in, rhminl_adj_land_in, rhminh_in )
         305             :
         306             :    ! --------------------------------------------------------- !
         307             :    ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
         308             :    ! analytical formulation of triangular PDF.                 !
         309             :    ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
         310             :    ! so using constant 'dV' assume that width is proportional  !
         311             :    ! to the saturation specific humidity.                      !
         312             :    !    dV ~ 0.1.                                              !
         313             :    !    cldrh : RH of in-stratus( = 1 if no supersaturation)   !
         314             :    ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
         315             :    ! G is discontinuous across U = 1.  In fact, it does not    !
         316             :    ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that   !
         317             :    ! they will produce the same results.                       !
         318             :    ! --------------------------------------------------------- !
         319             :
         320             :    real(r8), intent(in)  :: U_in(pcols)           ! Relative humidity
         321             :    real(r8), intent(in)  :: p_in(pcols)           ! Pressure [Pa]
         322             :    real(r8), intent(in)  :: qv_in(pcols)          ! Grid-mean water vapor specific humidity [kg/kg]
         323             :    real(r8), intent(in)  :: landfrac_in(pcols)    ! Land fraction
         324             :    real(r8), intent(in)  :: snowh_in(pcols)       ! Snow depth (liquid water equivalent)
         325             :
         326             :    real(r8), intent(out) :: a_out(pcols)          ! Stratus fraction
         327             :    real(r8), intent(out) :: Ga_out(pcols)         ! dU/da
         328             :    integer,  intent(in)  :: ncol
         329             :
         330             :    real(r8), optional, intent(in)  :: rhminl_in(pcols)                ! Critical relative humidity for low-level  liquid stratus
         331             :    real(r8), optional, intent(in)  :: rhminl_adj_land_in(pcols)       ! Adjustment drop of rhminl over the land
         332             :    real(r8), optional, intent(in)  :: rhminh_in(pcols)                ! Critical relative humidity for high-level liquid stratus
         333             :
         334             :    real(r8)              :: rhminl                ! Critical relative humidity for low-level  liquid stratus
         335             :    real(r8)              :: rhminl_adj_land       ! Adjustment drop of rhminl over the land
         336             :    real(r8)              :: rhminh                ! Critical relative humidity for high-level liquid stratus
         337             :
         338             :    real(r8)              :: U                     ! Relative humidity
         339             :    real(r8)              :: p                     ! Pressure [Pa]
         340             :    real(r8)              :: qv                    ! Grid-mean water vapor specific humidity [kg/kg]
         341             :    real(r8)              :: landfrac              ! Land fraction
         342             :    real(r8)              :: snowh                 ! Snow depth (liquid water equivalent)
         343             :
         344             :    real(r8)              :: a                     ! Stratus fraction
         345             :    real(r8)              :: Ga                    ! dU/da
         346             :
         347             :    ! Local variables
         348             :    integer :: i                                   ! Loop indexes
         349             :    real(r8) dV                                    ! Width of triangular PDF
         350             :    real(r8) cldrh                                 ! RH of stratus cloud
         351             :    real(r8) rhmin                                 ! Critical RH
         352             :    real(r8) rhwght
         353             :
         354             :    ! Statement functions
         355             :    logical land
         356             :    land(i) = nint(landfrac_in(i)) == 1
         357             :
         358             :    ! ---------- !
         359             :    ! Parameters !
         360             :    ! ---------- !
         361             :
         362           0 :    cldrh  = 1.0_r8
         363             :
         364           0 :    rhminl          = rhminl_const
         365           0 :    rhminl_adj_land = rhminl_adj_land_const
         366           0 :    rhminh          = rhminh_const
         367             :
         368             :    ! ---------------- !
         369             :    ! Main computation !
         370             :    ! ---------------- !
         371             :
         372           0 :    a_out(:)  = 0._r8
         373           0 :    Ga_out(:) = 0._r8
         374             :
         375           0 :    do i = 1, ncol
         376             :
         377           0 :    U        = U_in(i)
         378           0 :    p        = p_in(i)
         379           0 :    qv       = qv_in(i)
         380           0 :    landfrac = landfrac_in(i)
         381           0 :    snowh    = snowh_in(i)
         382             :
         383           0 :    if (present(rhminl_in))          rhminl          = rhminl_in(i)
         384           0 :    if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i)
         385           0 :    if (present(rhminh_in))          rhminh          = rhminh_in(i)
         386             :
         387           0 :    if( p .ge. premib ) then
         388             :
         389           0 :        if( land(i) .and. (snowh.le.0.000001_r8) ) then
         390           0 :            rhmin = rhminl - rhminl_adj_land
         391             :        else
         392             :            rhmin = rhminl
         393             :        endif
         394             :
         395           0 :        dV = cldrh - rhmin
         396             :
         397           0 :        if( U .ge. 1._r8 ) then
         398             :            a  = 1._r8
         399             :            Ga = 1.e10_r8
         400           0 :        elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
         401           0 :            a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
         402           0 :            Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
         403           0 :        elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
         404             :            a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
         405           0 :                       (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
         406           0 :            Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
         407           0 :        elseif( U .le. (cldrh-dV) ) then
         408           0 :            a  = 0._r8
         409           0 :            Ga = 1.e10_r8
         410             :        endif
         411             :
         412             :        if( freeze_dry ) then
         413             :            a  = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         414             :            Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         415             :        endif
         416             :
         417           0 :    elseif( p .lt. premit ) then
         418             :
         419           0 :        rhmin = rhminh
         420           0 :        dV    = cldrh - rhmin
         421             :
         422           0 :        if( U .ge. 1._r8 ) then
         423             :            a  = 1._r8
         424             :            Ga = 1.e10_r8
         425           0 :        elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
         426           0 :            a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
         427           0 :            Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
         428           0 :        elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
         429             :            a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
         430           0 :                       (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
         431           0 :            Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
         432           0 :        elseif( U .le. (cldrh-dV) ) then
         433           0 :            a  = 0._r8
         434           0 :            Ga = 1.e10_r8
         435             :        endif
         436             :
         437             :    else
         438             :
         439           0 :        rhwght = (premib-(max(p,premit)))/(premib-premit)
         440             :
         441             :      ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
         442             :      !     rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
         443             :      ! else
         444           0 :            rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
         445             :      ! endif
         446             :
         447           0 :        dV    = cldrh - rhmin
         448             :
         449           0 :        if( U .ge. 1._r8 ) then
         450             :            a  = 1._r8
         451             :            Ga = 1.e10_r8
         452           0 :        elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
         453           0 :            a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
         454           0 :            Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
         455           0 :        elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
         456             :            a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
         457           0 :                          (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
         458           0 :            Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
         459           0 :        elseif( U .le. (cldrh-dV) ) then
         460           0 :            a  = 0._r8
         461           0 :            Ga = 1.e10_r8
         462             :        endif
         463             :
         464             :    endif
         465             :
         466           0 :    a_out(i)  = a
         467           0 :    Ga_out(i) = Ga
         468             :
         469             :    enddo
         470             :
         471           0 : end subroutine astG_PDF
         472             : !================================================================================================
         473             :
         474           0 : subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, orhmin, &
         475             :                               rhminl_in, rhminl_adj_land_in, rhminh_in )
         476             :
         477             :    ! --------------------------------------------------------- !
         478             :    ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
         479             :    ! CAM35 cloud fraction formula.                             !
         480             :    ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !
         481             :    ! For the other cases, I should re-define 'rhminl,rhminh' & !
         482             :    ! 'premib,premit'.                                          !
         483             :    ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
         484             :    ! G is discontinuous across U = 1.                          !
         485             :    ! --------------------------------------------------------- !
         486             :
         487             :    real(r8), intent(in)  :: U               ! Relative humidity
         488             :    real(r8), intent(in)  :: p               ! Pressure [Pa]
         489             :    real(r8), intent(in)  :: qv              ! Grid-mean water vapor specific humidity [kg/kg]
         490             :    real(r8), intent(in)  :: landfrac        ! Land fraction
         491             :    real(r8), intent(in)  :: snowh           ! Snow depth (liquid water equivalent)
         492             :
         493             :    real(r8), intent(out) :: a               ! Stratus fraction
         494             :    real(r8), intent(out) :: Ga              ! dU/da
         495             :    real(r8), optional, intent(out) :: orhmin ! Critical RH
         496             :
         497             :    real(r8), optional, intent(in)  :: rhminl_in          ! Critical relative humidity for low-level  liquid stratus
         498             :    real(r8), optional, intent(in)  :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land
         499             :    real(r8), optional, intent(in)  :: rhminh_in          ! Critical relative humidity for high-level liquid stratus
         500             :
         501             :    ! Local variables
         502             :    real(r8) rhmin                                 ! Critical RH
         503             :    real(r8) rhdif                                 ! Factor for stratus fraction
         504             :    real(r8) rhwght
         505             :
         506             :    real(r8) :: rhminl
         507             :    real(r8) :: rhminl_adj_land
         508             :    real(r8) :: rhminh
         509             :
         510             :    ! Statement functions
         511             :    logical land
         512           0 :    land = nint(landfrac) == 1
         513             :
         514           0 :    rhminl = rhminl_const
         515           0 :    if (present(rhminl_in)) rhminl = rhminl_in
         516           0 :    rhminl_adj_land = rhminl_adj_land_const
         517           0 :    if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in
         518           0 :    rhminh = rhminh_const
         519           0 :    if (present(rhminh_in)) rhminh = rhminh_in
         520             :
         521             :    ! ---------------- !
         522             :    ! Main computation !
         523             :    ! ---------------- !
         524             :
         525           0 :    if( p .ge. premib ) then
         526             :
         527           0 :        if( land .and. (snowh.le.0.000001_r8) ) then
         528           0 :            rhmin = rhminl - rhminl_adj_land
         529             :        else
         530             :            rhmin = rhminl
         531             :        endif
         532           0 :        rhdif = (U-rhmin)/(1.0_r8-rhmin)
         533           0 :        a  = min(1._r8,(max(rhdif,0.0_r8))**2)
         534           0 :        if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
         535           0 :             Ga = 1.e20_r8
         536             :        else
         537           0 :             Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
         538             :        endif
         539             :        if( freeze_dry ) then
         540             :            a  = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         541             :            Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         542             :        endif
         543             :
         544           0 :    elseif( p .lt. premit ) then
         545             :
         546           0 :        rhmin = rhminh
         547           0 :        rhdif = (U-rhmin)/(1.0_r8-rhmin)
         548           0 :        a  = min(1._r8,(max(rhdif,0._r8))**2)
         549           0 :        if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
         550           0 :             Ga = 1.e20_r8
         551             :        else
         552           0 :             Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
         553             :        endif
         554             :
         555             :    else
         556             :
         557           0 :        rhwght = (premib-(max(p,premit)))/(premib-premit)
         558             :
         559             :      ! if( land .and. (snowh.le.0.000001_r8) ) then
         560             :      !     rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
         561             :      ! else
         562           0 :            rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
         563             :      ! endif
         564             :
         565           0 :        rhdif = (U-rhmin)/(1.0_r8-rhmin)
         566           0 :        a  = min(1._r8,(max(rhdif,0._r8))**2)
         567           0 :        if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
         568           0 :             Ga = 1.e10_r8
         569             :        else
         570           0 :             Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
         571             :        endif
         572             :
         573             :    endif
         574             :
         575           0 :    if (present(orhmin)) orhmin = rhmin
         576             :
         577           0 : end subroutine astG_RHU_single
         578             :
         579             : !================================================================================================
         580             :
         581           0 : subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, &
         582             :                     rhminl_in, rhminl_adj_land_in, rhminh_in )
         583             :
         584             :    ! --------------------------------------------------------- !
         585             :    ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
         586             :    ! CAM35 cloud fraction formula.                             !
         587             :    ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !
         588             :    ! For the other cases, I should re-define 'rhminl,rhminh' & !
         589             :    ! 'premib,premit'.                                          !
         590             :    ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
         591             :    ! G is discontinuous across U = 1.                          !
         592             :    ! --------------------------------------------------------- !
         593             :
         594             :    real(r8), intent(in)  :: U_in(pcols)           ! Relative humidity
         595             :    real(r8), intent(in)  :: p_in(pcols)           ! Pressure [Pa]
         596             :    real(r8), intent(in)  :: qv_in(pcols)          ! Grid-mean water vapor specific humidity [kg/kg]
         597             :    real(r8), intent(in)  :: landfrac_in(pcols)    ! Land fraction
         598             :    real(r8), intent(in)  :: snowh_in(pcols)       ! Snow depth (liquid water equivalent)
         599             :
         600             :    real(r8), intent(out) :: a_out(pcols)          ! Stratus fraction
         601             :    real(r8), intent(out) :: Ga_out(pcols)         ! dU/da
         602             :    integer,  intent(in)  :: ncol
         603             :
         604             :    real(r8), optional, intent(in)  :: rhminl_in(pcols)          ! Critical relative humidity for low-level  liquid stratus
         605             :    real(r8), optional, intent(in)  :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land
         606             :    real(r8), optional, intent(in)  :: rhminh_in(pcols)          ! Critical relative humidity for high-level liquid stratus
         607             :
         608             :    real(r8)              :: U                     ! Relative humidity
         609             :    real(r8)              :: p                     ! Pressure [Pa]
         610             :    real(r8)              :: qv                    ! Grid-mean water vapor specific humidity [kg/kg]
         611             :    real(r8)              :: landfrac              ! Land fraction
         612             :    real(r8)              :: snowh                 ! Snow depth (liquid water equivalent)
         613             :
         614             :    real(r8)              :: rhminl                ! Critical relative humidity for low-level  liquid stratus
         615             :    real(r8)              :: rhminl_adj_land       ! Adjustment drop of rhminl over the land
         616             :    real(r8)              :: rhminh                ! Critical relative humidity for high-level liquid stratus
         617             :
         618             :    real(r8)              :: a                     ! Stratus fraction
         619             :    real(r8)              :: Ga                    ! dU/da
         620             :
         621             :    ! Local variables
         622             :    integer  i
         623             :    real(r8) rhmin                                 ! Critical RH
         624             :    real(r8) rhdif                                 ! Factor for stratus fraction
         625             :    real(r8) rhwght
         626             :
         627             :    ! Statement functions
         628             :    logical land
         629             :    land(i) = nint(landfrac_in(i)) == 1
         630             :
         631           0 :    rhminl          = rhminl_const
         632           0 :    rhminl_adj_land = rhminl_adj_land_const
         633           0 :    rhminh          = rhminh_const
         634             :
         635             :    ! ---------------- !
         636             :    ! Main computation !
         637             :    ! ---------------- !
         638             :
         639           0 :    a_out(:) = 0._r8
         640           0 :    Ga_out(:) = 0._r8
         641             :
         642           0 :    do i = 1, ncol
         643             :
         644           0 :    U        = U_in(i)
         645           0 :    p        = p_in(i)
         646           0 :    qv       = qv_in(i)
         647           0 :    landfrac = landfrac_in(i)
         648           0 :    snowh    = snowh_in(i)
         649             :
         650           0 :    if (present(rhminl_in))          rhminl          = rhminl_in(i)
         651           0 :    if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i)
         652           0 :    if (present(rhminh_in))          rhminh          = rhminh_in(i)
         653             :
         654           0 :    if( p .ge. premib ) then
         655             :
         656           0 :        if( land(i) .and. (snowh.le.0.000001_r8) ) then
         657           0 :            rhmin = rhminl - rhminl_adj_land
         658             :        else
         659             :            rhmin = rhminl
         660             :        endif
         661           0 :        rhdif = (U-rhmin)/(1.0_r8-rhmin)
         662           0 :        a  = min(1._r8,(max(rhdif,0.0_r8))**2)
         663           0 :        if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
         664             :             Ga = 1.e20_r8
         665             :        else
         666           0 :             Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
         667             :        endif
         668             :        if( freeze_dry ) then
         669             :            a  = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         670             :            Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
         671             :        endif
         672             :
         673           0 :    elseif( p .lt. premit ) then
         674             :
         675           0 :        rhmin = rhminh
         676           0 :        rhdif = (U-rhmin)/(1.0_r8-rhmin)
         677           0 :        a  = min(1._r8,(max(rhdif,0._r8))**2)
         678           0 :        if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
         679             :             Ga = 1.e20_r8
         680             :        else
         681           0 :             Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
         682             :        endif
         683             :
         684             :    else
         685             :
         686           0 :        rhwght = (premib-(max(p,premit)))/(premib-premit)
         687             :
         688             :      ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
         689             :      !     rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
         690             :      ! else
         691           0 :            rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
         692             :      ! endif
         693             :
         694           0 :        rhdif = (U-rhmin)/(1.0_r8-rhmin)
         695           0 :        a  = min(1._r8,(max(rhdif,0._r8))**2)
         696           0 :        if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
         697             :             Ga = 1.e10_r8
         698             :        else
         699           0 :             Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
         700             :        endif
         701             :
         702             :    endif
         703             :
         704           0 :    a_out(i)  = a
         705           0 :    Ga_out(i) = Ga
         706             :
         707             :    enddo
         708             :
         709           0 : end subroutine astG_RHU
         710             :
         711             : !================================================================================================
         712             :
         713           0 : subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, &
         714             :                        rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, &
         715             :                        qsatfac_out)
         716             :
         717             :    ! --------------------------------------------------------- !
         718             :    ! Compute non-physical ice stratus fraction                 !
         719             :    ! --------------------------------------------------------- !
         720             :
         721             :    real(r8), intent(in)  :: qv              ! Grid-mean water vapor[kg/kg]
         722             :    real(r8), intent(in)  :: T               ! Temperature
         723             :    real(r8), intent(in)  :: p               ! Pressure [Pa]
         724             :    real(r8), intent(in)  :: qi              ! Grid-mean ice water content [kg/kg]
         725             :    real(r8), intent(in)  :: landfrac        ! Land fraction
         726             :    real(r8), intent(in)  :: snowh           ! Snow depth (liquid water equivalent)
         727             :
         728             :    real(r8), intent(out) :: aist            ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
         729             :
         730             :    real(r8), optional, intent(in)  :: rhmaxi_in
         731             :    real(r8), optional, intent(in)  :: rhmini_in          ! Critical relative humidity for               ice stratus
         732             :    real(r8), optional, intent(in)  :: rhminl_in          ! Critical relative humidity for low-level  liquid stratus
         733             :    real(r8), optional, intent(in)  :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land
         734             :    real(r8), optional, intent(in)  :: rhminh_in          ! Critical relative humidity for high-level liquid stratus
         735             :    real(r8), optional, intent(out) :: qsatfac_out        ! Subgrid scaling factor for qsat
         736             :
         737             :    ! Local variables
         738             :    real(r8) rhmin                           ! Critical RH
         739             :    real(r8) rhwght
         740             :
         741             :    real(r8) a,b,c,as,bs,cs                  ! Fit parameters
         742             :    real(r8) Kc                              ! Constant for ice cloud calc (wood & field)
         743             :    real(r8) ttmp                            ! Limited temperature
         744             :    real(r8) icicval                         ! Empirical IWC value [ kg/kg ]
         745             :    real(r8) rho                             ! Local air density
         746             :    real(r8) esl                             ! Liq sat vapor pressure
         747             :    real(r8) esi                             ! Ice sat vapor pressure
         748             :    real(r8) ncf,phi                         ! Wilson and Ballard parameters
         749             :    real(r8) es, qs
         750             :
         751             :    real(r8) rhi                             ! grid box averaged relative humidity over ice
         752             :    real(r8) minice                          ! minimum grid box avg ice for having a 'cloud'
         753             :    real(r8) mincld                          ! minimum ice cloud fraction threshold
         754             :    real(r8) icimr                           ! in cloud ice mixing ratio
         755             :    real(r8) rhdif                           ! working variable for slingo scheme
         756             :
         757             :    real(r8) :: rhmaxi
         758             :    real(r8) :: rhmini
         759             :    real(r8) :: rhminl
         760             :    real(r8) :: rhminl_adj_land
         761             :    real(r8) :: rhminh
         762             :
         763             :    ! Statement functions
         764             :    logical land
         765           0 :    land = nint(landfrac) == 1
         766             :
         767             :    ! --------- !
         768             :    ! Constants !
         769             :    ! --------- !
         770             :
         771             :    ! Wang and Sassen IWC paramters ( Option.1 )
         772           0 :      a = 26.87_r8
         773           0 :      b = 0.569_r8
         774           0 :      c = 0.002892_r8
         775             :    ! Schiller parameters ( Option.2 )
         776           0 :      as = -68.4202_r8
         777           0 :      bs = 0.983917_r8
         778           0 :      cs = 2.81795_r8
         779             :    ! Wood and Field parameters ( Option.3 )
         780           0 :      Kc = 75._r8
         781             :    ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
         782             :    ! Slingo modified (option 5)
         783           0 :      minice = 1.e-12_r8
         784           0 :      mincld = 1.e-4_r8
         785             :
         786           0 :    rhmaxi = rhmaxi_const
         787           0 :    if (present(rhmaxi_in)) rhmaxi = rhmaxi_in
         788           0 :    rhmini = rhmini_const
         789           0 :    if (present(rhmini_in)) rhmini = rhmini_in
         790           0 :    rhminl = rhminl_const
         791           0 :    if (present(rhminl_in)) rhminl = rhminl_in
         792           0 :    rhminl_adj_land = rhminl_adj_land_const
         793           0 :    if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in
         794           0 :    rhminh = rhminh_const
         795           0 :    if (present(rhminh_in)) rhminh = rhminh_in
         796           0 :    if (present(qsatfac_out)) qsatfac_out = 1.0_r8
         797             :
         798             :
         799             :    ! ---------------- !
         800             :    ! Main computation !
         801             :    ! ---------------- !
         802             :
         803           0 :      call qsat_water(T, p, es, qs)
         804           0 :      esl = svp_water(T)
         805           0 :      esi = svp_ice(T)
         806             :
         807           0 :      if( iceopt.lt.3 ) then
         808           0 :          if( iceopt.eq.1 ) then
         809           0 :              ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
         810           0 :              icicval = a + b * ttmp + c * ttmp**2._r8
         811           0 :              rho = p/(rair*T)
         812           0 :              icicval = icicval * 1.e-6_r8 / rho
         813             :          else
         814           0 :              ttmp = max(190._r8,min(T,273.16_r8))
         815           0 :              icicval = 10._r8 **(as * bs**ttmp + cs)
         816           0 :              icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
         817             :          endif
         818           0 :          aist =  max(0._r8,min(qi/icicval,1._r8))
         819           0 :      elseif( iceopt.eq.3 ) then
         820           0 :          aist = 1._r8 - exp(-Kc*qi/(qs*(esi/esl)))
         821           0 :          aist = max(0._r8,min(aist,1._r8))
         822           0 :      elseif( iceopt.eq.4) then
         823           0 :          if( p .ge. premib ) then
         824           0 :              if( land .and. (snowh.le.0.000001_r8) ) then
         825             :                  rhmin = rhminl - rhminl_adj_land
         826             :              else
         827             :                  rhmin = rhminl
         828             :              endif
         829           0 :          elseif( p .lt. premit ) then
         830             :              rhmin = rhminh
         831             :          else
         832           0 :              rhwght = (premib-(max(p,premit)))/(premib-premit)
         833             :            ! if( land .and. (snowh.le.0.000001_r8) ) then
         834             :            !     rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
         835             :            ! else
         836           0 :                  rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
         837             :            ! endif
         838             :          endif
         839           0 :          ncf = qi/((1._r8 - icecrit)*qs)
         840           0 :          if( ncf.le.0._r8 ) then
         841           0 :              aist = 0._r8
         842           0 :          elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then
         843           0 :              aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
         844           0 :          elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
         845           0 :              phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
         846           0 :              aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
         847             :          else
         848           0 :              aist = 1._r8
         849             :          endif
         850           0 :              aist = max(0._r8,min(aist,1._r8))
         851           0 :      elseif (iceopt.eq.5) then
         852             :         ! set rh ice cloud fraction
         853           0 :         rhi= (qv+qi)/qs * (esl/esi)
         854           0 :         if (rhmaxi .eq. rhmini) then
         855           0 :            if (rhi .gt. rhmini) then
         856             :               rhdif = 1._r8
         857             :            else
         858           0 :               rhdif = 0._r8
         859             :            end if
         860             :         else
         861           0 :            rhdif = (rhi-rhmini) / (rhmaxi - rhmini)
         862             :         end if
         863           0 :         aist = min(1.0_r8, max(rhdif,0._r8)**2)
         864             :
         865             :         ! Similar to alpha in Wilson & Ballard (1999), determine a
         866             :         ! scaling factor for saturation vapor pressure that reflects
         867             :         ! the cloud fraction, rhmini, and rhmaxi.
         868             :         !
         869             :         ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less.
         870           0 :         if (present(qsatfac_out) .and. cldfrc2m_do_subgrid_growth) then
         871           0 :            qsatfac_out = max(min(qv / qs, 1._r8), (1._r8 - aist) * rhmini + aist * rhmaxi)
         872             :         end if
         873             :
         874             :         ! limiter to remove empty cloud and ice with no cloud
         875             :         ! and set icecld fraction to mincld if ice exists
         876             :
         877           0 :         if (qi.lt.minice) then
         878           0 :            aist=0._r8
         879             :         else
         880           0 :            aist=max(mincld,aist)
         881             :         endif
         882             :
         883             :         ! enforce limits on icimr
         884           0 :         if (qi.ge.minice) then
         885           0 :            icimr=qi/aist
         886             :
         887             :            !minimum
         888           0 :            if (icimr.lt.cldfrc2m_qist_min) then
         889           0 :              if (cldfrc2m_do_avg_aist_algs) then
         890             :                !
         891             :                ! Take the geometric mean of the iceopt=4 and iceopt=5 values.
         892             :                ! Mods developed by Thomas Toniazzo for NorESM.
         893           0 :                aist = max(0._r8,min(1._r8,sqrt(aist*qi/cldfrc2m_qist_min)))
         894             :              else
         895             :                !
         896             :                ! Default for iceopt=5
         897           0 :                aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_min))
         898             :              end if
         899             :            endif
         900             :            !maximum
         901           0 :            if (icimr.gt.cldfrc2m_qist_max) then
         902           0 :               aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_max))
         903             :            endif
         904             :
         905             :         endif
         906             :      endif
         907             :
         908             :    ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
         909             :    ! computed after updating 'qi_st'.
         910             :
         911           0 :      aist = max(0._r8,min(aist,0.999_r8))
         912             :
         913           0 : end subroutine aist_single
         914             :
         915             : !================================================================================================
         916             :
         917           0 : subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, aist_out, ncol, &
         918             :                        rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, &
         919             :                        qsatfac_out )
         920             :
         921             :    ! --------------------------------------------------------- !
         922             :    ! Compute non-physical ice stratus fraction                 !
         923             :    ! --------------------------------------------------------- !
         924             :
         925             :    real(r8), intent(in)  :: qv_in(pcols)       ! Grid-mean water vapor[kg/kg]
         926             :    real(r8), intent(in)  :: T_in(pcols)        ! Temperature
         927             :    real(r8), intent(in)  :: p_in(pcols)        ! Pressure [Pa]
         928             :    real(r8), intent(in)  :: qi_in(pcols)       ! Grid-mean ice water content [kg/kg]
         929             :    real(r8), intent(in)  :: ni_in(pcols)       ! Grid-mean ice water number concentration [#/kg]
         930             :    real(r8), intent(in)  :: landfrac_in(pcols) ! Land fraction
         931             :    real(r8), intent(in)  :: snowh_in(pcols)    ! Snow depth (liquid water equivalent)
         932             :
         933             :    real(r8), intent(out) :: aist_out(pcols)    ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
         934             :    integer,  intent(in)  :: ncol
         935             :
         936             :    real(r8), optional, intent(in)  :: rhmaxi_in(pcols)
         937             :    real(r8), optional, intent(in)  :: rhmini_in(pcols)          ! Critical relative humidity for               ice stratus
         938             :    real(r8), optional, intent(in)  :: rhminl_in(pcols)          ! Critical relative humidity for low-level  liquid stratus
         939             :    real(r8), optional, intent(in)  :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land
         940             :    real(r8), optional, intent(in)  :: rhminh_in(pcols)          ! Critical relative humidity for high-level liquid stratus
         941             :    real(r8), optional, intent(out) :: qsatfac_out(pcols)        ! Subgrid scaling factor for qsat
         942             :
         943             :    ! Local variables
         944             :
         945             :    real(r8) qv                              ! Grid-mean water vapor[kg/kg]
         946             :    real(r8) T                               ! Temperature
         947             :    real(r8) p                               ! Pressure [Pa]
         948             :    real(r8) qi                              ! Grid-mean ice water content [kg/kg]
         949             :    real(r8) ni
         950             :    real(r8) landfrac                        ! Land fraction
         951             :    real(r8) snowh                           ! Snow depth (liquid water equivalent)
         952             :
         953             :    real(r8) rhmaxi                          ! Critical relative humidity for               ice stratus
         954             :    real(r8) rhmini                          ! Critical relative humidity for               ice stratus
         955             :    real(r8) rhminl                          ! Critical relative humidity for low-level  liquid stratus
         956             :    real(r8) rhminl_adj_land                 ! Adjustment drop of rhminl over the land
         957             :    real(r8) rhminh                          ! Critical relative humidity for high-level liquid stratus
         958             :
         959             :    real(r8) aist                            ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
         960             :
         961             :    real(r8) rhmin                           ! Critical RH
         962             :    real(r8) rhwght
         963             :
         964             :    real(r8) a,b,c,as,bs,cs,ah,bh,ch         ! Fit parameters
         965             :    real(r8) nil
         966             :    real(r8) Kc                              ! Constant for ice cloud calc (wood & field)
         967             :    real(r8) ttmp                            ! Limited temperature
         968             :    real(r8) icicval                         ! Empirical IWC value [ kg/kg ]
         969             :    real(r8) rho                             ! Local air density
         970             :    real(r8) esl(pcols)                      ! Liq sat vapor pressure
         971             :    real(r8) esi(pcols)                      ! Ice sat vapor pressure
         972             :    real(r8) ncf,phi                         ! Wilson and Ballard parameters
         973             :    real(r8) qs
         974             :    real(r8) esat_in(pcols)
         975             :    real(r8) qsat_in(pcols)
         976             :
         977             :    real(r8) rhi                             ! grid box averaged relative humidity over ice
         978             :    real(r8) minice                          ! minimum grid box avg ice for having a 'cloud'
         979             :    real(r8) mincld                          ! minimum ice cloud fraction threshold
         980             :    real(r8) icimr                           ! in cloud ice mixing ratio
         981             :    real(r8) rhdif                           ! working variable for slingo scheme
         982             :
         983             :    integer i
         984             :
         985             :
         986             :    ! Statement functions
         987             :    logical land
         988             :    land(i) = nint(landfrac_in(i)) == 1
         989             :
         990             :    ! --------- !
         991             :    ! Constants !
         992             :    ! --------- !
         993             :
         994             :    ! Wang and Sassen IWC paramters ( Option.1 )
         995           0 :      a = 26.87_r8
         996           0 :      b = 0.569_r8
         997           0 :      c = 0.002892_r8
         998             :    ! Schiller parameters ( Option.2 )
         999           0 :      as = -68.4202_r8
        1000           0 :      bs = 0.983917_r8
        1001           0 :      cs = 2.81795_r8
        1002             :    ! Wood and Field parameters ( Option.3 )
        1003           0 :      Kc = 75._r8
        1004             :    ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
        1005             :    ! Slingo modified (option 5)
        1006           0 :      minice = 1.e-12_r8
        1007           0 :      mincld = 1.e-4_r8
        1008             :
        1009           0 :      rhmaxi          = rhmaxi_const
        1010             :
        1011           0 :      rhmini          = rhmini_const
        1012           0 :      rhminl          = rhminl_const
        1013           0 :      rhminl_adj_land = rhminl_adj_land_const
        1014           0 :      rhminh          = rhminh_const
        1015             :
        1016           0 :      if (present(qsatfac_out)) qsatfac_out = 1.0_r8
        1017             :
        1018             :    ! ---------------- !
        1019             :    ! Main computation !
        1020             :    ! ---------------- !
        1021             :
        1022           0 :      aist_out(:) = 0._r8
        1023           0 :      esat_in(:)  = 0._r8
        1024           0 :      qsat_in(:)  = 0._r8
        1025             :
        1026           0 :      call qsat_water(T_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol)
        1027           0 :      call svp_water_vect(T_in(1:ncol), esl(1:ncol), ncol)
        1028           0 :      call svp_ice_vect(T_in(1:ncol), esi(1:ncol), ncol)
        1029             :
        1030           0 :      do i = 1, ncol
        1031             :
        1032           0 :        landfrac = landfrac_in(i)
        1033           0 :        snowh = snowh_in(i)
        1034           0 :        T = T_in(i)
        1035           0 :        qv = qv_in(i)
        1036           0 :        p = p_in(i)
        1037           0 :        qi = qi_in(i)
        1038           0 :        ni = ni_in(i)
        1039           0 :        qs = qsat_in(i)
        1040             :
        1041           0 :        if (present(rhmaxi_in))          rhmaxi          = rhmaxi_in(i)
        1042           0 :        if (present(rhmini_in))          rhmini          = rhmini_in(i)
        1043           0 :        if (present(rhminl_in))          rhminl          = rhminl_in(i)
        1044           0 :        if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i)
        1045           0 :        if (present(rhminh_in))          rhminh          = rhminh_in(i)
        1046             :
        1047           0 :        if( iceopt.lt.3 ) then
        1048           0 :          if( iceopt.eq.1 ) then
        1049           0 :              ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
        1050           0 :              icicval = a + b * ttmp + c * ttmp**2._r8
        1051           0 :              rho = p/(rair*T)
        1052           0 :              icicval = icicval * 1.e-6_r8 / rho
        1053             :          else
        1054           0 :              ttmp = max(190._r8,min(T,273.16_r8))
        1055           0 :              icicval = 10._r8 **(as * bs**ttmp + cs)
        1056           0 :              icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
        1057             :          endif
        1058           0 :          aist =  max(0._r8,min(qi/icicval,1._r8))
        1059           0 :        elseif( iceopt.eq.3 ) then
        1060           0 :          aist = 1._r8 - exp(-Kc*qi/(qs*(esi(i)/esl(i))))
        1061           0 :          aist = max(0._r8,min(aist,1._r8))
        1062           0 :        elseif( iceopt.eq.4) then
        1063           0 :          if( p .ge. premib ) then
        1064           0 :              if( land(i) .and. (snowh.le.0.000001_r8) ) then
        1065             :                  rhmin = rhminl - rhminl_adj_land
        1066             :              else
        1067             :                  rhmin = rhminl
        1068             :              endif
        1069           0 :          elseif( p .lt. premit ) then
        1070             :              rhmin = rhminh
        1071             :          else
        1072           0 :              rhwght = (premib-(max(p,premit)))/(premib-premit)
        1073             :            ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
        1074             :            !     rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
        1075             :            ! else
        1076           0 :                  rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
        1077             :            ! endif
        1078             :          endif
        1079           0 :          ncf = qi/((1._r8 - icecrit)*qs)
        1080           0 :          if( ncf.le.0._r8 ) then
        1081             :              aist = 0._r8
        1082           0 :          elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then
        1083           0 :              aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
        1084           0 :          elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
        1085           0 :              phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
        1086           0 :              aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
        1087             :          else
        1088             :              aist = 1._r8
        1089             :          endif
        1090           0 :              aist = max(0._r8,min(aist,1._r8))
        1091           0 :        elseif (iceopt.eq.5) then
        1092             :          ! set rh ice cloud fraction
        1093           0 :          rhi= (qv+qi)/qs * (esl(i)/esi(i))
        1094           0 :          if (rhmaxi .eq. rhmini) then
        1095           0 :            if (rhi .gt. rhmini) then
        1096             :               rhdif = 1._r8
        1097             :            else
        1098           0 :               rhdif = 0._r8
        1099             :            end if
        1100             :          else
        1101           0 :            rhdif = (rhi-rhmini) / (rhmaxi - rhmini)
        1102             :          end if
        1103           0 :          aist = min(1.0_r8, max(rhdif,0._r8)**2)
        1104             :
        1105           0 :        elseif (iceopt.eq.6) then
        1106             :          !----- ICE CLOUD OPTION 6: fit based on T and Number (Gettelman: based on Heymsfield obs)
        1107             :          ! Use observations from Heymsfield et al 2012 of IWC and Ni v. Temp
        1108             :          ! Multivariate fit follows form of Boudala 2002: ICIWC = a * exp(b*T) * N^c
        1109             :          ! a=6.73e-8, b=0.05, c=0.349
        1110             :          ! N is #/L, so need to convert Ni_L=N*rhoa/1000.
        1111           0 :          ah= 6.73834e-08_r8
        1112           0 :          bh= 0.0533110_r8
        1113           0 :          ch= 0.3493813_r8
        1114           0 :          rho=p/(rair*T)
        1115           0 :          nil=ni*rho/1000._r8
        1116           0 :          icicval = ah * exp(bh*T) * nil**ch
        1117             :          !result is in g m-3, convert to kg H2O / kg air (icimr...)
        1118           0 :          icicval = icicval / rho / 1000._r8
        1119           0 :          aist =  max(0._r8,min(qi/icicval,1._r8))
        1120           0 :          aist =  min(aist,1._r8)
        1121             :
        1122             :        endif
        1123             :
        1124           0 :        if (iceopt.eq.5 .or. iceopt.eq.6) then
        1125             :
        1126             :          ! Similar to alpha in Wilson & Ballard (1999), determine a
        1127             :          ! scaling factor for saturation vapor pressure that reflects
        1128             :          ! the cloud fraction, rhmini, and rhmaxi.
        1129             :          !
        1130             :          ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less.
        1131           0 :          if (present(qsatfac_out) .and. cldfrc2m_do_subgrid_growth) then
        1132           0 :            qsatfac_out(i) = max(min(qv / qs, 1._r8), (1._r8 - aist) * rhmini + aist * rhmaxi)
        1133             :          end if
        1134             :
        1135             :          ! limiter to remove empty cloud and ice with no cloud
        1136             :          ! and set icecld fraction to mincld if ice exists
        1137             :
        1138           0 :          if (qi.lt.minice) then
        1139             :            aist=0._r8
        1140             :          else
        1141           0 :            aist=max(mincld,aist)
        1142             :          endif
        1143             :
        1144             :          ! enforce limits on icimr
        1145           0 :          if (qi.ge.minice) then
        1146           0 :            icimr=qi/aist
        1147             :
        1148             :            !minimum
        1149           0 :            if (icimr.lt.cldfrc2m_qist_min) then
        1150           0 :              if (cldfrc2m_do_avg_aist_algs) then
        1151             :                !
        1152             :                ! Take the geometric mean of the iceopt=4 and iceopt=5 values.
        1153             :                ! Mods developed by Thomas Toniazzo for NorESM.
        1154           0 :                aist = max(0._r8,min(1._r8,sqrt(aist*qi/cldfrc2m_qist_min)))
        1155             :              else
        1156             :                !
        1157             :                ! Default for iceopt=5
        1158           0 :                aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_min))
        1159             :              end if
        1160             :            endif
        1161             :            !maximum
        1162           0 :            if (icimr.gt.cldfrc2m_qist_max) then
        1163           0 :               aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_max))
        1164             :            endif
        1165             :
        1166             :          endif
        1167             :        endif
        1168             :
        1169             :        ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
        1170             :        ! computed after updating 'qi_st'.
        1171             :
        1172           0 :        aist = max(0._r8,min(aist,0.999_r8))
        1173             :
        1174           0 :        aist_out(i) = aist
        1175             :
        1176             :      enddo
        1177             :
        1178           0 : end subroutine aist_vector
        1179             :
        1180             : !================================================================================================
        1181             :
        1182             : end module cldfrc2m

Generated by: LCOV version 1.14