From a0f02d09e4c5e08580091cac3c9a9638d74493a5 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 31 Mar 2026 13:10:39 -0400 Subject: [PATCH 1/9] Bring in cldfrc2m (compute_cloud_fraction_two_moment) cherry-picks of 51c81a7e part-f6e0bda6 part-be911bdff from uwshcu dev --- .../compute_cloud_fraction_two_moment.F90 | 1086 +++++++++++++++++ .../compute_cloud_fraction_two_moment.meta | 86 ++ ...ute_cloud_fraction_two_moment_namelist.xml | 113 ++ 3 files changed, 1285 insertions(+) create mode 100644 schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 create mode 100644 schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta create mode 100644 schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 new file mode 100644 index 00000000..f5964886 --- /dev/null +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -0,0 +1,1086 @@ +! Compute cloud fractions. +! This was used in combination with the two-moment MG microphysics in CAM, +! hence the "two moment" distinction. Formerly called "cldfrc2m". +! +! Provides liquid stratus (PDF and RHU methods) and ice stratus cloud fraction. +! Original authors: Sungsu Park, Andrew Gettelman +module compute_cloud_fraction_two_moment + + use ccpp_kinds, only: kind_phys + implicit none + private + + ! CCPP-compliant subroutines: + public :: compute_cloud_fraction_two_moment_init + + ! Underlying subroutines that are not CCPP schemes + ! but are used by other schemes and internally use the parameters + ! read from the above init phase. + public :: astG_PDF_single, astG_PDF + public :: astG_RHU_single, astG_RHU + public :: aist_single, aist_vector + + !REMOVECAM: this public parameter is used by CAM code and this public declaration can be removed + ! after CAM is retired and all code uses this parameter is migrated to receive this via + ! the CCPP framework. + public :: CAMstfrac + !REMOVECAM_END + + ! Hardcoded parameters + logical, parameter :: CAMstfrac = .false. ! If .true. (.false.), + ! use Slingo (triangular PDF-based) liquid stratus fraction + logical, parameter :: freeze_dry = .false. ! If .true., use 'freeze dry' in liquid stratus fraction formula + + ! Module variables set by init (used internally by subroutines) + real(kind_phys) :: premit ! Top pressure bound for mid-level clouds [Pa] + real(kind_phys) :: premib ! Bottom pressure bound for mid-level clouds [Pa] + integer :: iceopt ! Ice cloud closure option + ! 1=wang & sassen 2=schiller (iciwc) + ! 3=wood & field, 4=Wilson (based on smith) + ! 5=modified slingo (ssat & empty cloud) + real(kind_phys) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure + ! ( smaller = more ice clouds ) + real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [ kg/kg ] + real(kind_phys) :: qist_max ! Maximum in-stratus IWC constraint [ kg/kg ] + logical :: do_subgrid_growth + logical :: do_avg_aist_algs + real(kind_phys) :: rair ! Gas constant of dry air [J kg-1 K-1] + +contains + +!================================================================================================ + +!> \section arg_table_compute_cloud_fraction_two_moment_init Argument Table +!! \htmlinclude compute_cloud_fraction_two_moment_init.html +subroutine compute_cloud_fraction_two_moment_init( & + amIRoot, iulog, & + premit_in, premib_in, iceopt_in, icecrit_in, & + qist_min_in, qist_max_in, do_subgrid_growth_in, do_avg_aist_algs_in, & + rair_in, & + errmsg, errflg) + + ! Input arguments + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + real(kind_phys), intent(in) :: premit_in + real(kind_phys), intent(in) :: premib_in + integer, intent(in) :: iceopt_in + real(kind_phys), intent(in) :: icecrit_in + real(kind_phys), intent(in) :: qist_min_in + real(kind_phys), intent(in) :: qist_max_in + logical, intent(in) :: do_subgrid_growth_in + logical, intent(in) :: do_avg_aist_algs_in + real(kind_phys), intent(in) :: rair_in + + ! Output arguments + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errflg = 0 + errmsg = '' + + ! Set module variables from input arguments + premit = premit_in + premib = premib_in + iceopt = iceopt_in + icecrit = icecrit_in + qist_min = qist_min_in + qist_max = qist_max_in + do_subgrid_growth = do_subgrid_growth_in + do_avg_aist_algs = do_avg_aist_algs_in + rair = rair_in + + if(amIRoot) then + write(iulog,*) 'compute_cloud_fraction_two_moment_init parameters:' + write(iulog,*) ' premit = ', premit + write(iulog,*) ' premib = ', premib + write(iulog,*) ' iceopt = ', iceopt + write(iulog,*) ' icecrit = ', icecrit + write(iulog,*) ' qist_min = ', qist_min + write(iulog,*) ' qist_max = ', qist_max + write(iulog,*) ' do_subgrid_growth = ', do_subgrid_growth + write(iulog,*) ' do_avg_aist_algs = ', do_avg_aist_algs + end if + +end subroutine compute_cloud_fraction_two_moment_init + +!================================================================================================ + +subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! analytical formulation of triangular PDF. ! + ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! + ! so using constant 'dV' assume that width is proportional ! + ! to the saturation specific humidity. ! + ! dV ~ 0.1. ! + ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. In fact, it does not ! + ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! + ! they will produce the same results. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U ! Relative humidity + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a ! Stratus fraction + real(kind_phys), intent(out) :: Ga ! dU/da + + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys), optional, intent(out) :: orhmin ! Critical RH + + ! Local variables + real(kind_phys) dV ! Width of triangular PDF + real(kind_phys) cldrh ! RH of stratus cloud + real(kind_phys) rhmin ! Critical RH + real(kind_phys) rhwght + + ! Statement functions + logical land + land = nint(landfrac) == 1 + + ! ---------- ! + ! Parameters ! + ! ---------- ! + + cldrh = 1.0_kind_phys + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + if( p .ge. premib ) then + + if( land .and. (snowh.le.0.000001_kind_phys) ) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + endif + + dV = cldrh - rhmin + + if( U .ge. 1._kind_phys ) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) + elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) + elseif( U .le. (cldrh-dV) ) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + endif + + if( freeze_dry ) then + a = a *max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + endif + + elseif( p .lt. premit ) then + + rhmin = rhminh + dV = cldrh - rhmin + + if( U .ge. 1._kind_phys ) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) + elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) + elseif( U .le. (cldrh-dV) ) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + endif + + else + + rhwght = (premib-(max(p,premit)))/(premib-premit) + + ! if( land .and. (snowh.le.0.000001_kind_phys) ) then + ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) + ! else + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) + ! endif + + dV = cldrh - rhmin + + if( U .ge. 1._kind_phys ) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) + elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) + elseif( U .le. (cldrh-dV) ) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + endif + + endif + + if (present(orhmin)) orhmin = rhmin + +end subroutine astG_PDF_single + +!================================================================================================ + +subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & + rhminl_in, rhminl_adj_land_in, rhminh_in ) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! analytical formulation of triangular PDF. ! + ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! + ! so using constant 'dV' assume that width is proportional ! + ! to the saturation specific humidity. ! + ! dV ~ 0.1. ! + ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. In fact, it does not ! + ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! + ! they will produce the same results. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U_in(:) ! Relative humidity + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg/kg] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction + real(kind_phys), intent(out) :: Ga_out(:) ! dU/da + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: U ! Relative humidity + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg/kg] + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys) :: a ! Stratus fraction + real(kind_phys) :: Ga ! dU/da + + ! Local variables + integer :: i ! Loop indexes + real(kind_phys) :: dV ! Width of triangular PDF + real(kind_phys) :: cldrh ! RH of stratus cloud + real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhwght + + ! Statement functions + logical land + land(i) = nint(landfrac_in(i)) == 1 + + ! ---------- ! + ! Parameters ! + ! ---------- ! + + cldrh = 1.0_kind_phys + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + a_out(:) = 0._kind_phys + Ga_out(:) = 0._kind_phys + + do i = 1, ncol + + U = U_in(i) + p = p_in(i) + qv = qv_in(i) + landfrac = landfrac_in(i) + snowh = snowh_in(i) + + rhminl = rhminl_in(i) + rhminl_adj_land = rhminl_adj_land_in(i) + rhminh = rhminh_in(i) + + if( p .ge. premib ) then + + if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + endif + + dV = cldrh - rhmin + + if( U .ge. 1._kind_phys ) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) + elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) + elseif( U .le. (cldrh-dV) ) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + endif + + if( freeze_dry ) then + a = a *max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + endif + + elseif( p .lt. premit ) then + + rhmin = rhminh + dV = cldrh - rhmin + + if( U .ge. 1._kind_phys ) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) + elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) + elseif( U .le. (cldrh-dV) ) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + endif + + else + + rhwght = (premib-(max(p,premit)))/(premib-premit) + + ! if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then + ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) + ! else + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) + ! endif + + dV = cldrh - rhmin + + if( U .ge. 1._kind_phys ) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) + elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) + elseif( U .le. (cldrh-dV) ) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + endif + + endif + + a_out(i) = a + Ga_out(i) = Ga + + enddo + +end subroutine astG_PDF +!================================================================================================ + +subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! CAM35 cloud fraction formula. ! + ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! + ! For the other cases, I should re-define 'rhminl,rhminh' & ! + ! 'premib,premit'. ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U ! Relative humidity + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a ! Stratus fraction + real(kind_phys), intent(out) :: Ga ! dU/da + + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys), optional, intent(out) :: orhmin ! Critical RH + + ! Local variables + real(kind_phys) rhmin ! Critical RH + real(kind_phys) rhdif ! Factor for stratus fraction + real(kind_phys) rhwght + + ! Statement functions + logical land + land = nint(landfrac) == 1 + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + if( p .ge. premib ) then + + if( land .and. (snowh.le.0.000001_kind_phys) ) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + endif + rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) + a = min(1._kind_phys,(max(rhdif,0.0_kind_phys))**2) + if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) + endif + if( freeze_dry ) then + a = a*max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + endif + + elseif( p .lt. premit ) then + + rhmin = rhminh + rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) + a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) + if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) + endif + + else + + rhwght = (premib-(max(p,premit)))/(premib-premit) + + ! if( land .and. (snowh.le.0.000001_kind_phys) ) then + ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) + ! else + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) + ! endif + + rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) + a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) + if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then + Ga = 1.e10_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) + endif + + endif + + if (present(orhmin)) orhmin = rhmin + +end subroutine astG_RHU_single + +!================================================================================================ + +subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & + rhminl_in, rhminl_adj_land_in, rhminh_in ) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! CAM35 cloud fraction formula. ! + ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! + ! For the other cases, I should re-define 'rhminl,rhminh' & ! + ! 'premib,premit'. ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U_in(:) ! Relative humidity + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg/kg] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction + real(kind_phys), intent(out) :: Ga_out(:) ! dU/da + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: U ! Relative humidity + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg/kg] + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: a ! Stratus fraction + real(kind_phys) :: Ga ! dU/da + + ! Local variables + integer i + real(kind_phys) rhmin ! Critical RH + real(kind_phys) rhdif ! Factor for stratus fraction + real(kind_phys) rhwght + + ! Statement functions + logical land + land(i) = nint(landfrac_in(i)) == 1 + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + a_out(:) = 0._kind_phys + Ga_out(:) = 0._kind_phys + + do i = 1, ncol + + U = U_in(i) + p = p_in(i) + qv = qv_in(i) + landfrac = landfrac_in(i) + snowh = snowh_in(i) + + rhminl = rhminl_in(i) + rhminl_adj_land = rhminl_adj_land_in(i) + rhminh = rhminh_in(i) + + if( p .ge. premib ) then + + if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + endif + rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) + a = min(1._kind_phys,(max(rhdif,0.0_kind_phys))**2) + if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) + endif + if( freeze_dry ) then + a = a*max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) + endif + + elseif( p .lt. premit ) then + + rhmin = rhminh + rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) + a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) + if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) + endif + + else + + rhwght = (premib-(max(p,premit)))/(premib-premit) + + ! if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then + ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) + ! else + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) + ! endif + + rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) + a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) + if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then + Ga = 1.e10_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) + endif + + endif + + a_out(i) = a + Ga_out(i) = Ga + + enddo + +end subroutine astG_RHU + +!================================================================================================ + +subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & + rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh, & + qsatfac_out) + + ! --------------------------------------------------------- ! + ! Compute non-physical ice stratus fraction ! + ! --------------------------------------------------------- ! + + use wv_saturation, only: qsat_water, svp_water, svp_ice + + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor[kg/kg] + real(kind_phys), intent(in) :: T ! Temperature + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qi ! Grid-mean ice water content [kg/kg] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + + real(kind_phys), intent(in) :: rhmaxi + real(kind_phys), intent(in) :: rhmini ! Critical relative humidity for ice stratus + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + real(kind_phys), optional, intent(out) :: qsatfac_out ! Subgrid scaling factor for qsat + + ! Local variables + real(kind_phys) rhmin ! Critical RH + real(kind_phys) rhwght + + real(kind_phys) a,b,c,as,bs,cs ! Fit parameters + real(kind_phys) Kc ! Constant for ice cloud calc (wood & field) + real(kind_phys) ttmp ! Limited temperature + real(kind_phys) icicval ! Empirical IWC value [ kg/kg ] + real(kind_phys) rho ! Local air density + real(kind_phys) esl ! Liq sat vapor pressure + real(kind_phys) esi ! Ice sat vapor pressure + real(kind_phys) ncf,phi ! Wilson and Ballard parameters + real(kind_phys) es, qs + + real(kind_phys) rhi ! grid box averaged relative humidity over ice + real(kind_phys) minice ! minimum grid box avg ice for having a 'cloud' + real(kind_phys) mincld ! minimum ice cloud fraction threshold + real(kind_phys) icimr ! in cloud ice mixing ratio + real(kind_phys) rhdif ! working variable for slingo scheme + + ! Statement functions + logical land + land = nint(landfrac) == 1 + + ! --------- ! + ! Constants ! + ! --------- ! + + ! Wang and Sassen IWC paramters ( Option.1 ) + a = 26.87_kind_phys + b = 0.569_kind_phys + c = 0.002892_kind_phys + ! Schiller parameters ( Option.2 ) + as = -68.4202_kind_phys + bs = 0.983917_kind_phys + cs = 2.81795_kind_phys + ! Wood and Field parameters ( Option.3 ) + Kc = 75._kind_phys + ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) + ! Slingo modified (option 5) + minice = 1.e-12_kind_phys + mincld = 1.e-4_kind_phys + + if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys + + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + call qsat_water(T, p, es, qs) + esl = svp_water(T) + esi = svp_ice(T) + + if( iceopt.lt.3 ) then + if( iceopt.eq.1 ) then + ttmp = max(195._kind_phys,min(T,253._kind_phys)) - 273.16_kind_phys + icicval = a + b * ttmp + c * ttmp**2._kind_phys + rho = p/(rair*T) + icicval = icicval * 1.e-6_kind_phys / rho + else + ttmp = max(190._kind_phys,min(T,273.16_kind_phys)) + icicval = 10._kind_phys **(as * bs**ttmp + cs) + icicval = icicval * 1.e-6_kind_phys * 18._kind_phys / 28.97_kind_phys + endif + aist = max(0._kind_phys,min(qi/icicval,1._kind_phys)) + elseif( iceopt.eq.3 ) then + aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi/esl))) + aist = max(0._kind_phys,min(aist,1._kind_phys)) + elseif( iceopt.eq.4) then + if( p .ge. premib ) then + if( land .and. (snowh.le.0.000001_kind_phys) ) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + endif + elseif( p .lt. premit ) then + rhmin = rhminh + else + rhwght = (premib-(max(p,premit)))/(premib-premit) + ! if( land .and. (snowh.le.0.000001_kind_phys) ) then + ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) + ! else + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) + ! endif + endif + ncf = qi/((1._kind_phys - icecrit)*qs) + if( ncf.le.0._kind_phys ) then + aist = 0._kind_phys + elseif( ncf.gt.0._kind_phys .and. ncf.le.1._kind_phys/6._kind_phys ) then + aist = 0.5_kind_phys*(6._kind_phys * ncf)**(2._kind_phys/3._kind_phys) + elseif( ncf.gt.1._kind_phys/6._kind_phys .and. ncf.lt.1._kind_phys ) then + phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+4._kind_phys*3.1415927_kind_phys)/3._kind_phys + aist = (1._kind_phys - 4._kind_phys * cos(phi) * cos(phi)) + else + aist = 1._kind_phys + endif + aist = max(0._kind_phys,min(aist,1._kind_phys)) + elseif (iceopt.eq.5) then + ! set rh ice cloud fraction + rhi= (qv+qi)/qs * (esl/esi) + if (rhmaxi .eq. rhmini) then + if (rhi .gt. rhmini) then + rhdif = 1._kind_phys + else + rhdif = 0._kind_phys + end if + else + rhdif = (rhi-rhmini) / (rhmaxi - rhmini) + end if + aist = min(1.0_kind_phys, max(rhdif,0._kind_phys)**2) + + ! Similar to alpha in Wilson & Ballard (1999), determine a + ! scaling factor for saturation vapor pressure that reflects + ! the cloud fraction, rhmini, and rhmaxi. + ! + ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less. + if (present(qsatfac_out) .and. do_subgrid_growth) then + qsatfac_out = max(min(qv / qs, 1._kind_phys), (1._kind_phys - aist) * rhmini + aist * rhmaxi) + end if + + ! limiter to remove empty cloud and ice with no cloud + ! and set icecld fraction to mincld if ice exists + + if (qi.lt.minice) then + aist=0._kind_phys + else + aist=max(mincld,aist) + endif + + ! enforce limits on icimr + if (qi.ge.minice) then + icimr=qi/aist + + !minimum + if (icimr.lt.qist_min) then + if (do_avg_aist_algs) then + ! + ! Take the geometric mean of the iceopt=4 and iceopt=5 values. + ! Mods developed by Thomas Toniazzo for NorESM. + aist = max(0._kind_phys,min(1._kind_phys,sqrt(aist*qi/qist_min))) + else + ! + ! Default for iceopt=5 + aist = max(0._kind_phys,min(1._kind_phys,qi/qist_min)) + end if + endif + !maximum + if (icimr.gt.qist_max) then + aist = max(0._kind_phys,min(1._kind_phys,qi/qist_max)) + endif + + endif + endif + + ! 0.999_kind_phys is added to prevent infinite 'ql_st' at the end of instratus_condensate + ! computed after updating 'qi_st'. + + aist = max(0._kind_phys,min(aist,0.999_kind_phys)) + +end subroutine aist_single + +!================================================================================================ + +subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, aist_out, ncol, & + rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, & + qsatfac_out ) + + ! --------------------------------------------------------- ! + ! Compute non-physical ice stratus fraction ! + ! --------------------------------------------------------- ! + + use wv_saturation, only: qsat_water, svp_water_vect, svp_ice_vect + + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor[kg/kg] + real(kind_phys), intent(in) :: T_in(:) ! Temperature + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qi_in(:) ! Grid-mean ice water content [kg/kg] + real(kind_phys), intent(in) :: ni_in(:) ! Grid-mean ice water number concentration [#/kg] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: aist_out(:) ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhmaxi_in(:) + real(kind_phys), intent(in) :: rhmini_in(:) ! Critical relative humidity for ice stratus + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + real(kind_phys), optional, intent(out) :: qsatfac_out(:) ! Subgrid scaling factor for qsat + + ! Local variables + + real(kind_phys) qv ! Grid-mean water vapor[kg/kg] + real(kind_phys) T ! Temperature + real(kind_phys) p ! Pressure [Pa] + real(kind_phys) qi ! Grid-mean ice water content [kg/kg] + real(kind_phys) ni + real(kind_phys) landfrac ! Land fraction + real(kind_phys) snowh ! Snow depth (liquid water equivalent) + + real(kind_phys) rhmaxi ! Critical relative humidity for ice stratus + real(kind_phys) rhmini ! Critical relative humidity for ice stratus + real(kind_phys) rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + + real(kind_phys) rhmin ! Critical RH + real(kind_phys) rhwght + + real(kind_phys) a,b,c,as,bs,cs,ah,bh,ch ! Fit parameters + real(kind_phys) nil + real(kind_phys) Kc ! Constant for ice cloud calc (wood & field) + real(kind_phys) ttmp ! Limited temperature + real(kind_phys) icicval ! Empirical IWC value [ kg/kg ] + real(kind_phys) rho ! Local air density + real(kind_phys) esl(ncol) ! Liq sat vapor pressure + real(kind_phys) esi(ncol) ! Ice sat vapor pressure + real(kind_phys) ncf,phi ! Wilson and Ballard parameters + real(kind_phys) qs + real(kind_phys) esat_in(ncol) + real(kind_phys) qsat_in(ncol) + + real(kind_phys) rhi ! grid box averaged relative humidity over ice + real(kind_phys) minice ! minimum grid box avg ice for having a 'cloud' + real(kind_phys) mincld ! minimum ice cloud fraction threshold + real(kind_phys) icimr ! in cloud ice mixing ratio + real(kind_phys) rhdif ! working variable for slingo scheme + + integer :: i + + + ! Statement functions + logical land + land(i) = nint(landfrac_in(i)) == 1 + + ! --------- ! + ! Constants ! + ! --------- ! + + ! Wang and Sassen IWC paramters ( Option.1 ) + a = 26.87_kind_phys + b = 0.569_kind_phys + c = 0.002892_kind_phys + ! Schiller parameters ( Option.2 ) + as = -68.4202_kind_phys + bs = 0.983917_kind_phys + cs = 2.81795_kind_phys + ! Wood and Field parameters ( Option.3 ) + Kc = 75._kind_phys + ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) + ! Slingo modified (option 5) + minice = 1.e-12_kind_phys + mincld = 1.e-4_kind_phys + + if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + aist_out(:) = 0._kind_phys + esat_in(:) = 0._kind_phys + qsat_in(:) = 0._kind_phys + + call qsat_water(T_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol) + call svp_water_vect(T_in(1:ncol), esl(1:ncol), ncol) + call svp_ice_vect(T_in(1:ncol), esi(1:ncol), ncol) + + do i = 1, ncol + + landfrac = landfrac_in(i) + snowh = snowh_in(i) + T = T_in(i) + qv = qv_in(i) + p = p_in(i) + qi = qi_in(i) + ni = ni_in(i) + qs = qsat_in(i) + + rhmaxi = rhmaxi_in(i) + rhmini = rhmini_in(i) + rhminl = rhminl_in(i) + rhminl_adj_land = rhminl_adj_land_in(i) + rhminh = rhminh_in(i) + + if( iceopt.lt.3 ) then + if( iceopt.eq.1 ) then + ttmp = max(195._kind_phys,min(T,253._kind_phys)) - 273.16_kind_phys + icicval = a + b * ttmp + c * ttmp**2._kind_phys + rho = p/(rair*T) + icicval = icicval * 1.e-6_kind_phys / rho + else + ttmp = max(190._kind_phys,min(T,273.16_kind_phys)) + icicval = 10._kind_phys **(as * bs**ttmp + cs) + icicval = icicval * 1.e-6_kind_phys * 18._kind_phys / 28.97_kind_phys + endif + aist = max(0._kind_phys,min(qi/icicval,1._kind_phys)) + elseif( iceopt.eq.3 ) then + aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi(i)/esl(i)))) + aist = max(0._kind_phys,min(aist,1._kind_phys)) + elseif( iceopt.eq.4) then + if( p .ge. premib ) then + if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + endif + elseif( p .lt. premit ) then + rhmin = rhminh + else + rhwght = (premib-(max(p,premit)))/(premib-premit) + ! if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then + ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) + ! else + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) + ! endif + endif + ncf = qi/((1._kind_phys - icecrit)*qs) + if( ncf.le.0._kind_phys ) then + aist = 0._kind_phys + elseif( ncf.gt.0._kind_phys .and. ncf.le.1._kind_phys/6._kind_phys ) then + aist = 0.5_kind_phys*(6._kind_phys * ncf)**(2._kind_phys/3._kind_phys) + elseif( ncf.gt.1._kind_phys/6._kind_phys .and. ncf.lt.1._kind_phys ) then + phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+4._kind_phys*3.1415927_kind_phys)/3._kind_phys + aist = (1._kind_phys - 4._kind_phys * cos(phi) * cos(phi)) + else + aist = 1._kind_phys + endif + aist = max(0._kind_phys,min(aist,1._kind_phys)) + elseif (iceopt.eq.5) then + ! set rh ice cloud fraction + rhi= (qv+qi)/qs * (esl(i)/esi(i)) + if (rhmaxi .eq. rhmini) then + if (rhi .gt. rhmini) then + rhdif = 1._kind_phys + else + rhdif = 0._kind_phys + end if + else + rhdif = (rhi-rhmini) / (rhmaxi - rhmini) + end if + aist = min(1.0_kind_phys, max(rhdif,0._kind_phys)**2) + + elseif (iceopt.eq.6) then + !----- ICE CLOUD OPTION 6: fit based on T and Number (Gettelman: based on Heymsfield obs) + ! Use observations from Heymsfield et al 2012 of IWC and Ni v. Temp + ! Multivariate fit follows form of Boudala 2002: ICIWC = a * exp(b*T) * N^c + ! a=6.73e-8, b=0.05, c=0.349 + ! N is #/L, so need to convert Ni_L=N*rhoa/1000. + ah= 6.73834e-08_kind_phys + bh= 0.0533110_kind_phys + ch= 0.3493813_kind_phys + rho=p/(rair*T) + nil=ni*rho/1000._kind_phys + icicval = ah * exp(bh*T) * nil**ch + !result is in g m-3, convert to kg H2O / kg air (icimr...) + icicval = icicval / rho / 1000._kind_phys + aist = max(0._kind_phys,min(qi/icicval,1._kind_phys)) + aist = min(aist,1._kind_phys) + + endif + + if (iceopt.eq.5 .or. iceopt.eq.6) then + + ! Similar to alpha in Wilson & Ballard (1999), determine a + ! scaling factor for saturation vapor pressure that reflects + ! the cloud fraction, rhmini, and rhmaxi. + ! + ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less. + if (present(qsatfac_out) .and. do_subgrid_growth) then + qsatfac_out(i) = max(min(qv / qs, 1._kind_phys), (1._kind_phys - aist) * rhmini + aist * rhmaxi) + end if + + ! limiter to remove empty cloud and ice with no cloud + ! and set icecld fraction to mincld if ice exists + + if (qi.lt.minice) then + aist=0._kind_phys + else + aist=max(mincld,aist) + endif + + ! enforce limits on icimr + if (qi.ge.minice) then + icimr=qi/aist + + !minimum + if (icimr.lt.qist_min) then + if (do_avg_aist_algs) then + ! + ! Take the geometric mean of the iceopt=4 and iceopt=5 values. + ! Mods developed by Thomas Toniazzo for NorESM. + aist = max(0._kind_phys,min(1._kind_phys,sqrt(aist*qi/qist_min))) + else + ! + ! Default for iceopt=5 + aist = max(0._kind_phys,min(1._kind_phys,qi/qist_min)) + end if + endif + !maximum + if (icimr.gt.qist_max) then + aist = max(0._kind_phys,min(1._kind_phys,qi/qist_max)) + endif + + endif + endif + + ! 0.999_kind_phys is added to prevent infinite 'ql_st' at the end of instratus_condensate + ! computed after updating 'qi_st'. + + aist = max(0._kind_phys,min(aist,0.999_kind_phys)) + + aist_out(i) = aist + + enddo + +end subroutine aist_vector + +!================================================================================================ + +end module compute_cloud_fraction_two_moment diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta new file mode 100644 index 00000000..23808fc3 --- /dev/null +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta @@ -0,0 +1,86 @@ +[ccpp-table-properties] + name = compute_cloud_fraction_two_moment + type = scheme + dependencies = ../../to_be_ccppized/wv_saturation.F90 + +[ccpp-arg-table] + name = compute_cloud_fraction_two_moment_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ premit_in ] + standard_name = tunable_parameter_for_top_pressure_bound_for_mid_level_clouds_for_cloud_fraction + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = in +[ premib_in ] + standard_name = tunable_parameter_for_bottom_pressure_bound_for_mid_level_liquid_stratus_for_cloud_fraction + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = in +[ iceopt_in ] + standard_name = control_for_ice_cloud_fraction + units = 1 + type = integer + dimensions = () + intent = in +[ icecrit_in ] + standard_name = tunable_parameter_for_critical_relative_humidity_for_ice_clouds_for_cloud_fraction_using_wilson_and_ballard_scheme + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ qist_min_in ] + standard_name = tunable_parameter_for_minimum_in_stratus_ice_water_content_for_two_moment_cloud_fraction + units = kg kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ qist_max_in ] + standard_name = tunable_parameter_for_maximum_in_stratus_ice_water_content_for_two_moment_cloud_fraction + units = kg kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ do_subgrid_growth_in ] + standard_name = do_subgrid_growth_for_two_moment_cloud_fraction + units = flag + type = logical + dimensions = () + intent = in +[ do_avg_aist_algs_in ] + standard_name = do_average_ice_stratus_algorithms_for_two_moment_cloud_fraction + units = flag + type = logical + dimensions = () + intent = in +[ rair_in ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml b/schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml new file mode 100644 index 00000000..978fb7a7 --- /dev/null +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml @@ -0,0 +1,113 @@ + + + + + + + + real + cldfrc2m + cldfrc2m_nl + tunable_parameter_for_minimum_relative_humidity_for_ice_stratus_for_two_moment_cloud_fraction + 1 + + Minimum relative humidity for ice cloud fraction greater than 0. + + + 0.80D0 + + + + real + cldfrc2m + cldfrc2m_nl + tunable_parameter_for_maximum_relative_humidity_for_ice_stratus_for_two_moment_cloud_fraction + 1 + + Maximum relative humidity for ice cloud fraction. + + + 1.10D0 + + + + real + cldfrc2m + cldfrc2m_nl + tunable_parameter_for_minimum_relative_humidity_for_ice_stratus_in_stratosphere_for_two_moment_cloud_fraction + 1 + + Minimum relative humidity for ice cloud fraction greater than 0 in the stratosphere. + + + 1.0D0 + + + + real + cldfrc2m + cldfrc2m_nl + tunable_parameter_for_maximum_relative_humidity_for_ice_stratus_in_stratosphere_for_two_moment_cloud_fraction + 1 + + Maximum relative humidity for ice cloud fraction in the stratosphere. + + + 1.0D0 + + + + real + cldfrc2m + cldfrc2m_nl + tunable_parameter_for_minimum_in_stratus_ice_water_content_for_two_moment_cloud_fraction + kg kg-1 + + Minimum in-stratus IWC constraint for iceopt=5. + + + 1.0D-7 + + + + real + cldfrc2m + cldfrc2m_nl + tunable_parameter_for_maximum_in_stratus_ice_water_content_for_two_moment_cloud_fraction + kg kg-1 + + Maximum in-stratus IWC constraint for iceopt=5. + + + 1.0D-3 + + + + logical + cldfrc2m + cldfrc2m_nl + do_subgrid_growth_for_two_moment_cloud_fraction + flag + + Enable subgrid-scale cloud growth for ice stratus (iceopt=5). + Default: FALSE + + + .false. + + + + logical + cldfrc2m + cldfrc2m_nl + do_average_ice_stratus_algorithms_for_two_moment_cloud_fraction + flag + + Use geometric mean of iceopt=4 and iceopt=5 algorithms for ice stratus. + Default: FALSE + + + .false. + + + From 3e5f92535c91dd1102f470220aa623c24fe75d63 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 12 Jun 2026 13:41:29 -0400 Subject: [PATCH 2/9] Apply suggestions from code review Co-authored-by: Jesse Nusbaumer --- .../compute_cloud_fraction_two_moment.F90 | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 index f5964886..3a362975 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -21,7 +21,7 @@ module compute_cloud_fraction_two_moment public :: aist_single, aist_vector !REMOVECAM: this public parameter is used by CAM code and this public declaration can be removed - ! after CAM is retired and all code uses this parameter is migrated to receive this via + ! after CAM is retired and all code that uses this parameter is migrated to receive this via ! the CCPP framework. public :: CAMstfrac !REMOVECAM_END @@ -40,8 +40,8 @@ module compute_cloud_fraction_two_moment ! 5=modified slingo (ssat & empty cloud) real(kind_phys) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure ! ( smaller = more ice clouds ) - real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [ kg/kg ] - real(kind_phys) :: qist_max ! Maximum in-stratus IWC constraint [ kg/kg ] + real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [ kg kg-1 ] + real(kind_phys) :: qist_max ! Maximum in-stratus IWC constraint [ kg kg-1 ] logical :: do_subgrid_growth logical :: do_avg_aist_algs real(kind_phys) :: rair ! Gas constant of dry air [J kg-1 K-1] @@ -157,7 +157,7 @@ subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ ! Main computation ! ! ---------------- ! - if( p .ge. premib ) then + if( p >= premib ) then if( land .and. (snowh.le.0.000001_kind_phys) ) then rhmin = rhminl - rhminl_adj_land @@ -510,7 +510,7 @@ subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, nco ! --------------------------------------------------------- ! ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! CAM35 cloud fraction formula. ! + ! CAM3.5 cloud fraction formula. ! ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! ! For the other cases, I should re-define 'rhminl,rhminh' & ! ! 'premib,premit'. ! @@ -688,16 +688,20 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & ! --------- ! ! Wang and Sassen IWC paramters ( Option.1 ) + ! DOI: 10.1175/1520-0469(2002)059<2291:CCMPRU>2.0.CO;2 a = 26.87_kind_phys b = 0.569_kind_phys c = 0.002892_kind_phys ! Schiller parameters ( Option.2 ) + ! DOI: 10.1029/2008JD010342 as = -68.4202_kind_phys bs = 0.983917_kind_phys cs = 2.81795_kind_phys ! Wood and Field parameters ( Option.3 ) + ! DOI: 10.1175/1520-0469(2000)057<1888:RBTWCW>2.0.CO;2 Kc = 75._kind_phys ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) + ! DOI: 10.1002/qj.49712555707 ! Slingo modified (option 5) minice = 1.e-12_kind_phys mincld = 1.e-4_kind_phys From 16c075aa8238aa7cd440c2f21d1b5c109cf1ae11 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 12 Jun 2026 13:54:12 -0400 Subject: [PATCH 3/9] Fix some formatting (take 1) --- .../compute_cloud_fraction_two_moment.F90 | 1962 ++++++++--------- 1 file changed, 953 insertions(+), 1009 deletions(-) diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 index 3a362975..f0fad5db 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -7,6 +7,7 @@ module compute_cloud_fraction_two_moment use ccpp_kinds, only: kind_phys + implicit none private @@ -27,24 +28,24 @@ module compute_cloud_fraction_two_moment !REMOVECAM_END ! Hardcoded parameters - logical, parameter :: CAMstfrac = .false. ! If .true. (.false.), - ! use Slingo (triangular PDF-based) liquid stratus fraction - logical, parameter :: freeze_dry = .false. ! If .true., use 'freeze dry' in liquid stratus fraction formula + logical, parameter :: CAMstfrac = .false. ! If .true. (.false.), + ! use Slingo (triangular PDF-based) liquid stratus fraction + logical, parameter :: freeze_dry = .false. ! If .true., use 'freeze dry' in liquid stratus fraction formula ! Module variables set by init (used internally by subroutines) - real(kind_phys) :: premit ! Top pressure bound for mid-level clouds [Pa] - real(kind_phys) :: premib ! Bottom pressure bound for mid-level clouds [Pa] - integer :: iceopt ! Ice cloud closure option - ! 1=wang & sassen 2=schiller (iciwc) - ! 3=wood & field, 4=Wilson (based on smith) - ! 5=modified slingo (ssat & empty cloud) - real(kind_phys) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure - ! ( smaller = more ice clouds ) - real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [ kg kg-1 ] - real(kind_phys) :: qist_max ! Maximum in-stratus IWC constraint [ kg kg-1 ] - logical :: do_subgrid_growth - logical :: do_avg_aist_algs - real(kind_phys) :: rair ! Gas constant of dry air [J kg-1 K-1] + real(kind_phys) :: premit ! Top pressure bound for mid-level clouds [Pa] + real(kind_phys) :: premib ! Bottom pressure bound for mid-level clouds [Pa] + integer :: iceopt ! Ice cloud closure option + ! 1=wang & sassen 2=schiller (iciwc) + ! 3=wood & field, 4=Wilson (based on smith) + ! 5=modified slingo (ssat & empty cloud) + real(kind_phys) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure + ! ( smaller = more ice clouds ) + real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [ kg kg-1 ] + real(kind_phys) :: qist_max ! Maximum in-stratus IWC constraint [ kg kg-1 ] + logical :: do_subgrid_growth + logical :: do_avg_aist_algs + real(kind_phys) :: rair ! Gas constant of dry air [J kg-1 K-1] contains @@ -52,728 +53,932 @@ module compute_cloud_fraction_two_moment !> \section arg_table_compute_cloud_fraction_two_moment_init Argument Table !! \htmlinclude compute_cloud_fraction_two_moment_init.html -subroutine compute_cloud_fraction_two_moment_init( & + subroutine compute_cloud_fraction_two_moment_init( & amIRoot, iulog, & premit_in, premib_in, iceopt_in, icecrit_in, & qist_min_in, qist_max_in, do_subgrid_growth_in, do_avg_aist_algs_in, & rair_in, & errmsg, errflg) - ! Input arguments - logical, intent(in) :: amIRoot - integer, intent(in) :: iulog - real(kind_phys), intent(in) :: premit_in - real(kind_phys), intent(in) :: premib_in - integer, intent(in) :: iceopt_in - real(kind_phys), intent(in) :: icecrit_in - real(kind_phys), intent(in) :: qist_min_in - real(kind_phys), intent(in) :: qist_max_in - logical, intent(in) :: do_subgrid_growth_in - logical, intent(in) :: do_avg_aist_algs_in - real(kind_phys), intent(in) :: rair_in - - ! Output arguments - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errflg - - errflg = 0 - errmsg = '' - - ! Set module variables from input arguments - premit = premit_in - premib = premib_in - iceopt = iceopt_in - icecrit = icecrit_in - qist_min = qist_min_in - qist_max = qist_max_in - do_subgrid_growth = do_subgrid_growth_in - do_avg_aist_algs = do_avg_aist_algs_in - rair = rair_in - - if(amIRoot) then - write(iulog,*) 'compute_cloud_fraction_two_moment_init parameters:' - write(iulog,*) ' premit = ', premit - write(iulog,*) ' premib = ', premib - write(iulog,*) ' iceopt = ', iceopt - write(iulog,*) ' icecrit = ', icecrit - write(iulog,*) ' qist_min = ', qist_min - write(iulog,*) ' qist_max = ', qist_max - write(iulog,*) ' do_subgrid_growth = ', do_subgrid_growth - write(iulog,*) ' do_avg_aist_algs = ', do_avg_aist_algs - end if - -end subroutine compute_cloud_fraction_two_moment_init + ! Input arguments + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + real(kind_phys), intent(in) :: premit_in + real(kind_phys), intent(in) :: premib_in + integer, intent(in) :: iceopt_in + real(kind_phys), intent(in) :: icecrit_in + real(kind_phys), intent(in) :: qist_min_in + real(kind_phys), intent(in) :: qist_max_in + logical, intent(in) :: do_subgrid_growth_in + logical, intent(in) :: do_avg_aist_algs_in + real(kind_phys), intent(in) :: rair_in + + ! Output arguments + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errflg = 0 + errmsg = '' + + ! Set module variables from input arguments + premit = premit_in + premib = premib_in + iceopt = iceopt_in + icecrit = icecrit_in + qist_min = qist_min_in + qist_max = qist_max_in + do_subgrid_growth = do_subgrid_growth_in + do_avg_aist_algs = do_avg_aist_algs_in + rair = rair_in + + if (amIRoot) then + write (iulog, *) 'compute_cloud_fraction_two_moment_init parameters:' + write (iulog, *) ' premit = ', premit + write (iulog, *) ' premib = ', premib + write (iulog, *) ' iceopt = ', iceopt + write (iulog, *) ' icecrit = ', icecrit + write (iulog, *) ' qist_min = ', qist_min + write (iulog, *) ' qist_max = ', qist_max + write (iulog, *) ' do_subgrid_growth = ', do_subgrid_growth + write (iulog, *) ' do_avg_aist_algs = ', do_avg_aist_algs + end if + + end subroutine compute_cloud_fraction_two_moment_init !================================================================================================ -subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) - - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! analytical formulation of triangular PDF. ! - ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! - ! so using constant 'dV' assume that width is proportional ! - ! to the saturation specific humidity. ! - ! dV ~ 0.1. ! - ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. In fact, it does not ! - ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! - ! they will produce the same results. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U ! Relative humidity - real(kind_phys), intent(in) :: p ! Pressure [Pa] - real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg] - real(kind_phys), intent(in) :: landfrac ! Land fraction - real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: a ! Stratus fraction - real(kind_phys), intent(out) :: Ga ! dU/da - - real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus - - real(kind_phys), optional, intent(out) :: orhmin ! Critical RH - - ! Local variables - real(kind_phys) dV ! Width of triangular PDF - real(kind_phys) cldrh ! RH of stratus cloud - real(kind_phys) rhmin ! Critical RH - real(kind_phys) rhwght - - ! Statement functions - logical land - land = nint(landfrac) == 1 - - ! ---------- ! - ! Parameters ! - ! ---------- ! - - cldrh = 1.0_kind_phys - - ! ---------------- ! - ! Main computation ! - ! ---------------- ! - - if( p >= premib ) then - - if( land .and. (snowh.le.0.000001_kind_phys) ) then - rhmin = rhminl - rhminl_adj_land - else - rhmin = rhminl - endif - - dV = cldrh - rhmin - - if( U .ge. 1._kind_phys ) then - a = 1._kind_phys - Ga = 1.e10_kind_phys - elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then - a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) - Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) - elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys - Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) - elseif( U .le. (cldrh-dV) ) then - a = 0._kind_phys - Ga = 1.e10_kind_phys - endif - - if( freeze_dry ) then - a = a *max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - endif - - elseif( p .lt. premit ) then - - rhmin = rhminh - dV = cldrh - rhmin - - if( U .ge. 1._kind_phys ) then - a = 1._kind_phys - Ga = 1.e10_kind_phys - elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then - a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) - Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) - elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys - Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) - elseif( U .le. (cldrh-dV) ) then - a = 0._kind_phys - Ga = 1.e10_kind_phys - endif - - else - - rhwght = (premib-(max(p,premit)))/(premib-premit) - - ! if( land .and. (snowh.le.0.000001_kind_phys) ) then - ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) - ! else - rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) - ! endif - - dV = cldrh - rhmin - - if( U .ge. 1._kind_phys ) then - a = 1._kind_phys - Ga = 1.e10_kind_phys - elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then - a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) - Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) - elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys - Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) - elseif( U .le. (cldrh-dV) ) then - a = 0._kind_phys - Ga = 1.e10_kind_phys - endif - - endif - - if (present(orhmin)) orhmin = rhmin - -end subroutine astG_PDF_single + subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! analytical formulation of triangular PDF. ! + ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! + ! so using constant 'dV' assume that width is proportional ! + ! to the saturation specific humidity. ! + ! dV ~ 0.1. ! + ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. In fact, it does not ! + ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! + ! they will produce the same results. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U ! Relative humidity + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a ! Stratus fraction + real(kind_phys), intent(out) :: Ga ! dU/da + + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys), optional, intent(out) :: orhmin ! Critical RH + + ! Local variables + real(kind_phys) :: dV ! Width of triangular PDF + real(kind_phys) :: cldrh ! RH of stratus cloud + real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhwght + + ! Statement functions + logical land + land = nint(landfrac) == 1 + + ! ---------- ! + ! Parameters ! + ! ---------- ! + + cldrh = 1.0_kind_phys + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + if (p >= premib) then + + if (land .and. (snowh <= 0.000001_kind_phys)) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + end if + + dV = cldrh - rhmin + + if (U >= 1._kind_phys) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif (U > (cldrh - dV/6._kind_phys) .and. U < 1._kind_phys) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) + elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) + elseif (U <= (cldrh - dV)) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + end if + + if (freeze_dry) then + a = a*max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + end if + + elseif (p < premit) then + + rhmin = rhminh + dV = cldrh - rhmin + + if (U >= 1._kind_phys) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif (U > (cldrh - dV/6._kind_phys) .and. U < 1._kind_phys) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) + elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) + elseif (U <= (cldrh - dV)) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + end if + + else + + rhwght = (premib - (max(p, premit)))/(premib - premit) + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) + + dV = cldrh - rhmin + + if (U >= 1._kind_phys) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif (U > (cldrh - dV/6._kind_phys) .and. U < 1._kind_phys) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) + elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) + elseif (U <= (cldrh - dV)) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + end if + + end if + + if (present(orhmin)) orhmin = rhmin + + end subroutine astG_PDF_single !================================================================================================ -subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & - rhminl_in, rhminl_adj_land_in, rhminh_in ) - - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! analytical formulation of triangular PDF. ! - ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! - ! so using constant 'dV' assume that width is proportional ! - ! to the saturation specific humidity. ! - ! dV ~ 0.1. ! - ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. In fact, it does not ! - ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! - ! they will produce the same results. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U_in(:) ! Relative humidity - real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] - real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg/kg] - real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction - real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction - real(kind_phys), intent(out) :: Ga_out(:) ! dU/da - integer, intent(in) :: ncol - - real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus - - real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus - - real(kind_phys) :: U ! Relative humidity - real(kind_phys) :: p ! Pressure [Pa] - real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg/kg] - real(kind_phys) :: landfrac ! Land fraction - real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) - - real(kind_phys) :: a ! Stratus fraction - real(kind_phys) :: Ga ! dU/da - - ! Local variables - integer :: i ! Loop indexes - real(kind_phys) :: dV ! Width of triangular PDF - real(kind_phys) :: cldrh ! RH of stratus cloud - real(kind_phys) :: rhmin ! Critical RH - real(kind_phys) :: rhwght - - ! Statement functions - logical land - land(i) = nint(landfrac_in(i)) == 1 - - ! ---------- ! - ! Parameters ! - ! ---------- ! - - cldrh = 1.0_kind_phys - - ! ---------------- ! - ! Main computation ! - ! ---------------- ! - - a_out(:) = 0._kind_phys - Ga_out(:) = 0._kind_phys - - do i = 1, ncol - - U = U_in(i) - p = p_in(i) - qv = qv_in(i) - landfrac = landfrac_in(i) - snowh = snowh_in(i) - - rhminl = rhminl_in(i) - rhminl_adj_land = rhminl_adj_land_in(i) - rhminh = rhminh_in(i) - - if( p .ge. premib ) then - - if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then - rhmin = rhminl - rhminl_adj_land - else - rhmin = rhminl - endif - - dV = cldrh - rhmin - - if( U .ge. 1._kind_phys ) then - a = 1._kind_phys - Ga = 1.e10_kind_phys - elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then - a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) - Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) - elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys - Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) - elseif( U .le. (cldrh-dV) ) then - a = 0._kind_phys - Ga = 1.e10_kind_phys - endif - - if( freeze_dry ) then - a = a *max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - endif - - elseif( p .lt. premit ) then - - rhmin = rhminh - dV = cldrh - rhmin - - if( U .ge. 1._kind_phys ) then - a = 1._kind_phys - Ga = 1.e10_kind_phys - elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then - a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) - Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) - elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys - Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) - elseif( U .le. (cldrh-dV) ) then - a = 0._kind_phys - Ga = 1.e10_kind_phys - endif - - else - - rhwght = (premib-(max(p,premit)))/(premib-premit) - - ! if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then - ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) - ! else - rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) - ! endif - - dV = cldrh - rhmin - - if( U .ge. 1._kind_phys ) then - a = 1._kind_phys - Ga = 1.e10_kind_phys - elseif( U .gt. (cldrh-dV/6._kind_phys) .and. U .lt. 1._kind_phys ) then - a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U-cldrh)/dV)**(2._kind_phys/3._kind_phys) - Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys-a) - elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._kind_phys) ) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys+(U-cldrh)/dV))-2._kind_phys*3.141592_kind_phys)))**2._kind_phys - Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a)-sqrt(a)) - elseif( U .le. (cldrh-dV) ) then - a = 0._kind_phys - Ga = 1.e10_kind_phys - endif - - endif - - a_out(i) = a - Ga_out(i) = Ga - - enddo - -end subroutine astG_PDF + subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & + rhminl_in, rhminl_adj_land_in, rhminh_in) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! analytical formulation of triangular PDF. ! + ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! + ! so using constant 'dV' assume that width is proportional ! + ! to the saturation specific humidity. ! + ! dV ~ 0.1. ! + ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. In fact, it does not ! + ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! + ! they will produce the same results. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U_in(:) ! Relative humidity + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction + real(kind_phys), intent(out) :: Ga_out(:) ! dU/da + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: U ! Relative humidity + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys) :: a ! Stratus fraction + real(kind_phys) :: Ga ! dU/da + + ! Local variables + integer :: i ! Loop indexes + real(kind_phys) :: dV ! Width of triangular PDF + real(kind_phys) :: cldrh ! RH of stratus cloud + real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhwght + + ! Statement functions + logical land + land(i) = nint(landfrac_in(i)) == 1 + + ! ---------- ! + ! Parameters ! + ! ---------- ! + + cldrh = 1.0_kind_phys + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + a_out(:) = 0._kind_phys + Ga_out(:) = 0._kind_phys + + do i = 1, ncol + U = U_in(i) + p = p_in(i) + qv = qv_in(i) + landfrac = landfrac_in(i) + snowh = snowh_in(i) + + rhminl = rhminl_in(i) + rhminl_adj_land = rhminl_adj_land_in(i) + rhminh = rhminh_in(i) + + if (p >= premib) then + if (land(i) .and. (snowh <= 0.000001_kind_phys)) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + end if + + dV = cldrh - rhmin + + if (U >= 1._kind_phys) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif (U > (cldrh - dV/6._kind_phys) .and. U < 1._kind_phys) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) + elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) + elseif (U <= (cldrh - dV)) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + end if + + if (freeze_dry) then + a = a*max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + end if + elseif (p < premit) then + rhmin = rhminh + dV = cldrh - rhmin + + if (U >= 1._kind_phys) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif (U > (cldrh - dV/6._kind_phys) .and. U < 1._kind_phys) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) + elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) + elseif (U <= (cldrh - dV)) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + end if + else + rhwght = (premib - (max(p, premit)))/(premib - premit) + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) + dV = cldrh - rhmin + + if (U >= 1._kind_phys) then + a = 1._kind_phys + Ga = 1.e10_kind_phys + elseif (U > (cldrh - dV/6._kind_phys) .and. U < 1._kind_phys) then + a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) + Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) + elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) + elseif (U <= (cldrh - dV)) then + a = 0._kind_phys + Ga = 1.e10_kind_phys + end if + end if + + a_out(i) = a + Ga_out(i) = Ga + end do + + end subroutine astG_PDF !================================================================================================ -subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) - - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! CAM35 cloud fraction formula. ! - ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! - ! For the other cases, I should re-define 'rhminl,rhminh' & ! - ! 'premib,premit'. ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U ! Relative humidity - real(kind_phys), intent(in) :: p ! Pressure [Pa] - real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg] - real(kind_phys), intent(in) :: landfrac ! Land fraction - real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: a ! Stratus fraction - real(kind_phys), intent(out) :: Ga ! dU/da - - real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus - - real(kind_phys), optional, intent(out) :: orhmin ! Critical RH - - ! Local variables - real(kind_phys) rhmin ! Critical RH - real(kind_phys) rhdif ! Factor for stratus fraction - real(kind_phys) rhwght - - ! Statement functions - logical land - land = nint(landfrac) == 1 - - ! ---------------- ! - ! Main computation ! - ! ---------------- ! - - if( p .ge. premib ) then - - if( land .and. (snowh.le.0.000001_kind_phys) ) then - rhmin = rhminl - rhminl_adj_land - else - rhmin = rhminl - endif - rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) - a = min(1._kind_phys,(max(rhdif,0.0_kind_phys))**2) - if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then - Ga = 1.e20_kind_phys - else - Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) - endif - if( freeze_dry ) then - a = a*max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - endif - - elseif( p .lt. premit ) then - - rhmin = rhminh - rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) - a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) - if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then - Ga = 1.e20_kind_phys - else - Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) - endif - - else - - rhwght = (premib-(max(p,premit)))/(premib-premit) - - ! if( land .and. (snowh.le.0.000001_kind_phys) ) then - ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) - ! else - rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) - ! endif - - rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) - a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) - if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then - Ga = 1.e10_kind_phys - else - Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) - endif - - endif - - if (present(orhmin)) orhmin = rhmin - -end subroutine astG_RHU_single + subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! CAM35 cloud fraction formula. ! + ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! + ! For the other cases, I should re-define 'rhminl,rhminh' & ! + ! 'premib,premit'. ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U ! Relative humidity + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a ! Stratus fraction + real(kind_phys), intent(out) :: Ga ! dU/da + + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys), optional, intent(out) :: orhmin ! Critical RH + + ! Local variables + real(kind_phys) rhmin ! Critical RH + real(kind_phys) rhdif ! Factor for stratus fraction + real(kind_phys) rhwght + + ! Statement functions + logical land + land = nint(landfrac) == 1 + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + if (p >= premib) then + + if (land .and. (snowh <= 0.000001_kind_phys)) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + end if + rhdif = (U - rhmin)/(1.0_kind_phys - rhmin) + a = min(1._kind_phys, (max(rhdif, 0.0_kind_phys))**2) + if ((U >= 1._kind_phys) .or. (U <= rhmin)) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys - rhmin)*((1._kind_phys - rhmin)/(U - rhmin)) + end if + if (freeze_dry) then + a = a*max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + end if + + else if (p < premit) then + rhmin = rhminh + rhdif = (U - rhmin)/(1.0_kind_phys - rhmin) + a = min(1._kind_phys, (max(rhdif, 0._kind_phys))**2) + if ((U >= 1._kind_phys) .or. (U <= rhmin)) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys - rhmin)*((1._kind_phys - rhmin)/(U - rhmin)) + end if + else + rhwght = (premib - (max(p, premit)))/(premib - premit) + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) + + rhdif = (U - rhmin)/(1.0_kind_phys - rhmin) + a = min(1._kind_phys, (max(rhdif, 0._kind_phys))**2) + if ((U >= 1._kind_phys) .or. (U <= rhmin)) then + Ga = 1.e10_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys - rhmin)*((1._kind_phys - rhmin)/(U - rhmin)) + end if + end if + + if (present(orhmin)) orhmin = rhmin + + end subroutine astG_RHU_single !================================================================================================ -subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & - rhminl_in, rhminl_adj_land_in, rhminh_in ) - - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! CAM3.5 cloud fraction formula. ! - ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! - ! For the other cases, I should re-define 'rhminl,rhminh' & ! - ! 'premib,premit'. ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U_in(:) ! Relative humidity - real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] - real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg/kg] - real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction - real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction - real(kind_phys), intent(out) :: Ga_out(:) ! dU/da - integer, intent(in) :: ncol - - real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus - - real(kind_phys) :: U ! Relative humidity - real(kind_phys) :: p ! Pressure [Pa] - real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg/kg] - real(kind_phys) :: landfrac ! Land fraction - real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) - - real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus - - real(kind_phys) :: a ! Stratus fraction - real(kind_phys) :: Ga ! dU/da - - ! Local variables - integer i - real(kind_phys) rhmin ! Critical RH - real(kind_phys) rhdif ! Factor for stratus fraction - real(kind_phys) rhwght - - ! Statement functions - logical land - land(i) = nint(landfrac_in(i)) == 1 - - ! ---------------- ! - ! Main computation ! - ! ---------------- ! - - a_out(:) = 0._kind_phys - Ga_out(:) = 0._kind_phys - - do i = 1, ncol - - U = U_in(i) - p = p_in(i) - qv = qv_in(i) - landfrac = landfrac_in(i) - snowh = snowh_in(i) - - rhminl = rhminl_in(i) - rhminl_adj_land = rhminl_adj_land_in(i) - rhminh = rhminh_in(i) - - if( p .ge. premib ) then - - if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then - rhmin = rhminl - rhminl_adj_land - else - rhmin = rhminl - endif - rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) - a = min(1._kind_phys,(max(rhdif,0.0_kind_phys))**2) - if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then - Ga = 1.e20_kind_phys - else - Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) - endif - if( freeze_dry ) then - a = a*max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - Ga = Ga/max(0.15_kind_phys,min(1.0_kind_phys,qv/0.0030_kind_phys)) - endif - - elseif( p .lt. premit ) then - - rhmin = rhminh - rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) - a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) - if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then - Ga = 1.e20_kind_phys - else - Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) - endif - - else - - rhwght = (premib-(max(p,premit)))/(premib-premit) - - ! if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then - ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) - ! else - rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) - ! endif - - rhdif = (U-rhmin)/(1.0_kind_phys-rhmin) - a = min(1._kind_phys,(max(rhdif,0._kind_phys))**2) - if( (U.ge.1._kind_phys) .or. (U.le.rhmin) ) then - Ga = 1.e10_kind_phys - else - Ga = 0.5_kind_phys*(1._kind_phys-rhmin)*((1._kind_phys-rhmin)/(U-rhmin)) - endif - - endif - - a_out(i) = a - Ga_out(i) = Ga - - enddo - -end subroutine astG_RHU + subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & + rhminl_in, rhminl_adj_land_in, rhminh_in) + + ! --------------------------------------------------------- ! + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! + ! CAM3.5 cloud fraction formula. ! + ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! + ! For the other cases, I should re-define 'rhminl,rhminh' & ! + ! 'premib,premit'. ! + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! + ! G is discontinuous across U = 1. ! + ! --------------------------------------------------------- ! + + real(kind_phys), intent(in) :: U_in(:) ! Relative humidity + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction + real(kind_phys), intent(out) :: Ga_out(:) ! dU/da + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: U ! Relative humidity + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: a ! Stratus fraction + real(kind_phys) :: Ga ! dU/da + + ! Local variables + integer i + real(kind_phys) rhmin ! Critical RH + real(kind_phys) rhdif ! Factor for stratus fraction + real(kind_phys) rhwght + + ! Statement functions + logical land + land(i) = nint(landfrac_in(i)) == 1 + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + a_out(:) = 0._kind_phys + Ga_out(:) = 0._kind_phys + + do i = 1, ncol + U = U_in(i) + p = p_in(i) + qv = qv_in(i) + landfrac = landfrac_in(i) + snowh = snowh_in(i) + + rhminl = rhminl_in(i) + rhminl_adj_land = rhminl_adj_land_in(i) + rhminh = rhminh_in(i) + + if (p >= premib) then + if (land(i) .and. (snowh <= 0.000001_kind_phys)) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + end if + rhdif = (U - rhmin)/(1.0_kind_phys - rhmin) + a = min(1._kind_phys, (max(rhdif, 0.0_kind_phys))**2) + if ((U >= 1._kind_phys) .or. (U <= rhmin)) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys - rhmin)*((1._kind_phys - rhmin)/(U - rhmin)) + end if + if (freeze_dry) then + a = a*max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + Ga = Ga/max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) + end if + else if (p < premit) then + rhmin = rhminh + rhdif = (U - rhmin)/(1.0_kind_phys - rhmin) + a = min(1._kind_phys, (max(rhdif, 0._kind_phys))**2) + if ((U >= 1._kind_phys) .or. (U <= rhmin)) then + Ga = 1.e20_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys - rhmin)*((1._kind_phys - rhmin)/(U - rhmin)) + end if + else + rhwght = (premib - (max(p, premit)))/(premib - premit) + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) + + rhdif = (U - rhmin)/(1.0_kind_phys - rhmin) + a = min(1._kind_phys, (max(rhdif, 0._kind_phys))**2) + if ((U >= 1._kind_phys) .or. (U <= rhmin)) then + Ga = 1.e10_kind_phys + else + Ga = 0.5_kind_phys*(1._kind_phys - rhmin)*((1._kind_phys - rhmin)/(U - rhmin)) + end if + end if + + a_out(i) = a + Ga_out(i) = Ga + end do + + end subroutine astG_RHU !================================================================================================ -subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & - rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh, & - qsatfac_out) - - ! --------------------------------------------------------- ! - ! Compute non-physical ice stratus fraction ! - ! --------------------------------------------------------- ! - - use wv_saturation, only: qsat_water, svp_water, svp_ice - - real(kind_phys), intent(in) :: qv ! Grid-mean water vapor[kg/kg] - real(kind_phys), intent(in) :: T ! Temperature - real(kind_phys), intent(in) :: p ! Pressure [Pa] - real(kind_phys), intent(in) :: qi ! Grid-mean ice water content [kg/kg] - real(kind_phys), intent(in) :: landfrac ! Land fraction - real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) - - real(kind_phys), intent(in) :: rhmaxi - real(kind_phys), intent(in) :: rhmini ! Critical relative humidity for ice stratus - real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus - real(kind_phys), optional, intent(out) :: qsatfac_out ! Subgrid scaling factor for qsat - - ! Local variables - real(kind_phys) rhmin ! Critical RH - real(kind_phys) rhwght - - real(kind_phys) a,b,c,as,bs,cs ! Fit parameters - real(kind_phys) Kc ! Constant for ice cloud calc (wood & field) - real(kind_phys) ttmp ! Limited temperature - real(kind_phys) icicval ! Empirical IWC value [ kg/kg ] - real(kind_phys) rho ! Local air density - real(kind_phys) esl ! Liq sat vapor pressure - real(kind_phys) esi ! Ice sat vapor pressure - real(kind_phys) ncf,phi ! Wilson and Ballard parameters - real(kind_phys) es, qs - - real(kind_phys) rhi ! grid box averaged relative humidity over ice - real(kind_phys) minice ! minimum grid box avg ice for having a 'cloud' - real(kind_phys) mincld ! minimum ice cloud fraction threshold - real(kind_phys) icimr ! in cloud ice mixing ratio - real(kind_phys) rhdif ! working variable for slingo scheme - - ! Statement functions - logical land - land = nint(landfrac) == 1 - - ! --------- ! - ! Constants ! - ! --------- ! - - ! Wang and Sassen IWC paramters ( Option.1 ) - ! DOI: 10.1175/1520-0469(2002)059<2291:CCMPRU>2.0.CO;2 - a = 26.87_kind_phys - b = 0.569_kind_phys - c = 0.002892_kind_phys - ! Schiller parameters ( Option.2 ) - ! DOI: 10.1029/2008JD010342 - as = -68.4202_kind_phys - bs = 0.983917_kind_phys - cs = 2.81795_kind_phys - ! Wood and Field parameters ( Option.3 ) - ! DOI: 10.1175/1520-0469(2000)057<1888:RBTWCW>2.0.CO;2 - Kc = 75._kind_phys - ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) - ! DOI: 10.1002/qj.49712555707 - ! Slingo modified (option 5) - minice = 1.e-12_kind_phys - mincld = 1.e-4_kind_phys - - if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys - - - ! ---------------- ! - ! Main computation ! - ! ---------------- ! - - call qsat_water(T, p, es, qs) - esl = svp_water(T) - esi = svp_ice(T) - - if( iceopt.lt.3 ) then - if( iceopt.eq.1 ) then - ttmp = max(195._kind_phys,min(T,253._kind_phys)) - 273.16_kind_phys - icicval = a + b * ttmp + c * ttmp**2._kind_phys - rho = p/(rair*T) - icicval = icicval * 1.e-6_kind_phys / rho - else - ttmp = max(190._kind_phys,min(T,273.16_kind_phys)) - icicval = 10._kind_phys **(as * bs**ttmp + cs) - icicval = icicval * 1.e-6_kind_phys * 18._kind_phys / 28.97_kind_phys - endif - aist = max(0._kind_phys,min(qi/icicval,1._kind_phys)) - elseif( iceopt.eq.3 ) then - aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi/esl))) - aist = max(0._kind_phys,min(aist,1._kind_phys)) - elseif( iceopt.eq.4) then - if( p .ge. premib ) then - if( land .and. (snowh.le.0.000001_kind_phys) ) then - rhmin = rhminl - rhminl_adj_land - else - rhmin = rhminl - endif - elseif( p .lt. premit ) then - rhmin = rhminh - else - rhwght = (premib-(max(p,premit)))/(premib-premit) - ! if( land .and. (snowh.le.0.000001_kind_phys) ) then - ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) - ! else - rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) - ! endif - endif - ncf = qi/((1._kind_phys - icecrit)*qs) - if( ncf.le.0._kind_phys ) then - aist = 0._kind_phys - elseif( ncf.gt.0._kind_phys .and. ncf.le.1._kind_phys/6._kind_phys ) then - aist = 0.5_kind_phys*(6._kind_phys * ncf)**(2._kind_phys/3._kind_phys) - elseif( ncf.gt.1._kind_phys/6._kind_phys .and. ncf.lt.1._kind_phys ) then - phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+4._kind_phys*3.1415927_kind_phys)/3._kind_phys - aist = (1._kind_phys - 4._kind_phys * cos(phi) * cos(phi)) - else - aist = 1._kind_phys - endif - aist = max(0._kind_phys,min(aist,1._kind_phys)) - elseif (iceopt.eq.5) then + subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & + rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh, & + qsatfac_out) + + ! --------------------------------------------------------- ! + ! Compute non-physical ice stratus fraction ! + ! --------------------------------------------------------- ! + + use wv_saturation, only: qsat_water, svp_water, svp_ice + + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor[kg kg-1] + real(kind_phys), intent(in) :: T ! Temperature + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qi ! Grid-mean ice water content [kg kg-1] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + + real(kind_phys), intent(in) :: rhmaxi + real(kind_phys), intent(in) :: rhmini ! Critical relative humidity for ice stratus + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + real(kind_phys), optional, intent(out) :: qsatfac_out ! Subgrid scaling factor for qsat + + ! Local variables + real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhwght + + real(kind_phys) :: a, b, c, as, bs, cs ! Fit parameters + real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) + real(kind_phys) :: ttmp ! Limited temperature + real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] + real(kind_phys) :: rho ! Local air density + real(kind_phys) :: esl ! Liq sat vapor pressure + real(kind_phys) :: esi ! Ice sat vapor pressure + real(kind_phys) :: ncf, phi ! Wilson and Ballard parameters + real(kind_phys) :: es, qs + + real(kind_phys) :: rhi ! grid box averaged relative humidity over ice + real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' + real(kind_phys) :: mincld ! minimum ice cloud fraction threshold + real(kind_phys) :: icimr ! in cloud ice mixing ratio + real(kind_phys) :: rhdif ! working variable for slingo scheme + + ! Statement functions + logical land + land = nint(landfrac) == 1 + + ! --------- ! + ! Constants ! + ! --------- ! + + ! Wang and Sassen IWC paramters ( Option.1 ) + ! DOI: 10.1175/1520-0469(2002)059<2291:CCMPRU>2.0.CO;2 + a = 26.87_kind_phys + b = 0.569_kind_phys + c = 0.002892_kind_phys + ! Schiller parameters ( Option.2 ) + ! DOI: 10.1029/2008JD010342 + as = -68.4202_kind_phys + bs = 0.983917_kind_phys + cs = 2.81795_kind_phys + ! Wood and Field parameters ( Option.3 ) + ! DOI: 10.1175/1520-0469(2000)057<1888:RBTWCW>2.0.CO;2 + Kc = 75._kind_phys + ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) + ! DOI: 10.1002/qj.49712555707 + ! Slingo modified (option 5) + minice = 1.e-12_kind_phys + mincld = 1.e-4_kind_phys + + if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + call qsat_water(T, p, es, qs) + esl = svp_water(T) + esi = svp_ice(T) + + if (iceopt < 3) then + if (iceopt == 1) then + ttmp = max(195._kind_phys, min(T, 253._kind_phys)) - 273.16_kind_phys + icicval = a + b*ttmp + c*ttmp**2._kind_phys + rho = p/(rair*T) + icicval = icicval*1.e-6_kind_phys/rho + else + ttmp = max(190._kind_phys, min(T, 273.16_kind_phys)) + icicval = 10._kind_phys**(as*bs**ttmp + cs) + icicval = icicval*1.e-6_kind_phys*18._kind_phys/28.97_kind_phys + end if + aist = max(0._kind_phys, min(qi/icicval, 1._kind_phys)) + elseif (iceopt == 3) then + aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi/esl))) + aist = max(0._kind_phys, min(aist, 1._kind_phys)) + elseif (iceopt == 4) then + if (p >= premib) then + if (land .and. (snowh <= 0.000001_kind_phys)) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + end if + else if (p < premit) then + rhmin = rhminh + else + rhwght = (premib - (max(p, premit)))/(premib - premit) + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) + end if + ncf = qi/((1._kind_phys - icecrit)*qs) + if (ncf <= 0._kind_phys) then + aist = 0._kind_phys + elseif (ncf > 0._kind_phys .and. ncf <= 1._kind_phys/6._kind_phys) then + aist = 0.5_kind_phys*(6._kind_phys*ncf)**(2._kind_phys/3._kind_phys) + elseif (ncf > 1._kind_phys/6._kind_phys .and. ncf < 1._kind_phys) then + phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+& + 4._kind_phys*3.1415927_kind_phys)/3._kind_phys + aist = (1._kind_phys - 4._kind_phys*cos(phi)*cos(phi)) + else + aist = 1._kind_phys + end if + aist = max(0._kind_phys, min(aist, 1._kind_phys)) + else if (iceopt == 5) then + ! set rh ice cloud fraction + rhi = (qv + qi)/qs*(esl/esi) + if (rhmaxi == rhmini) then + if (rhi > rhmini) then + rhdif = 1._kind_phys + else + rhdif = 0._kind_phys + end if + else + rhdif = (rhi - rhmini)/(rhmaxi - rhmini) + end if + aist = min(1.0_kind_phys, max(rhdif, 0._kind_phys)**2) + + ! Similar to alpha in Wilson & Ballard (1999), determine a + ! scaling factor for saturation vapor pressure that reflects + ! the cloud fraction, rhmini, and rhmaxi. + ! + ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less. + if (present(qsatfac_out) .and. do_subgrid_growth) then + qsatfac_out = max(min(qv/qs, 1._kind_phys), (1._kind_phys - aist)*rhmini + aist*rhmaxi) + end if + + ! limiter to remove empty cloud and ice with no cloud + ! and set icecld fraction to mincld if ice exists + + if (qi < minice) then + aist = 0._kind_phys + else + aist = max(mincld, aist) + end if + + ! enforce limits on icimr + if (qi >= minice) then + icimr = qi/aist + + !minimum + if (icimr < qist_min) then + if (do_avg_aist_algs) then + ! + ! Take the geometric mean of the iceopt=4 and iceopt=5 values. + ! Mods developed by Thomas Toniazzo for NorESM. + aist = max(0._kind_phys, min(1._kind_phys, sqrt(aist*qi/qist_min))) + else + ! + ! Default for iceopt=5 + aist = max(0._kind_phys, min(1._kind_phys, qi/qist_min)) + end if + end if + !maximum + if (icimr > qist_max) then + aist = max(0._kind_phys, min(1._kind_phys, qi/qist_max)) + end if + end if + end if + + ! 0.999_kind_phys is added to prevent infinite 'ql_st' at the end of instratus_condensate + ! computed after updating 'qi_st'. + aist = max(0._kind_phys, min(aist, 0.999_kind_phys)) + + end subroutine aist_single + +!================================================================================================ + + subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, aist_out, ncol, & + rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, & + qsatfac_out) + + ! --------------------------------------------------------- ! + ! Compute non-physical ice stratus fraction ! + ! --------------------------------------------------------- ! + + use wv_saturation, only: qsat_water, svp_water_vect, svp_ice_vect + + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor[kg kg-1] + real(kind_phys), intent(in) :: T_in(:) ! Temperature + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qi_in(:) ! Grid-mean ice water content [kg kg-1] + real(kind_phys), intent(in) :: ni_in(:) ! Grid-mean ice water number concentration [#/kg] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: aist_out(:) ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhmaxi_in(:) + real(kind_phys), intent(in) :: rhmini_in(:) ! Critical relative humidity for ice stratus + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + real(kind_phys), optional, intent(out) :: qsatfac_out(:) ! Subgrid scaling factor for qsat + + ! Local variables + + real(kind_phys) :: qv ! Grid-mean water vapor [kg kg-1] + real(kind_phys) :: T ! Temperature + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qi ! Grid-mean ice water content [kg kg-1] + real(kind_phys) :: ni + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys) :: rhmaxi ! Critical relative humidity for ice stratus + real(kind_phys) :: rhmini ! Critical relative humidity for ice stratus + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + + real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhwght + + real(kind_phys) :: a, b, c, as, bs, cs, ah, bh, ch! Fit parameters + real(kind_phys) :: nil + real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) + real(kind_phys) :: ttmp ! Limited temperature + real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] + real(kind_phys) :: rho ! Local air density + real(kind_phys) :: esl(ncol) ! Liq sat vapor pressure + real(kind_phys) :: esi(ncol) ! Ice sat vapor pressure + real(kind_phys) :: ncf, phi ! Wilson and Ballard parameters + real(kind_phys) :: qs + real(kind_phys) :: esat_in(ncol) + real(kind_phys) :: qsat_in(ncol) + + real(kind_phys) :: rhi ! grid box averaged relative humidity over ice + real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' + real(kind_phys) :: mincld ! minimum ice cloud fraction threshold + real(kind_phys) :: icimr ! in cloud ice mixing ratio + real(kind_phys) :: rhdif ! working variable for slingo scheme + + integer :: i + + ! Statement functions + logical land + land(i) = nint(landfrac_in(i)) == 1 + + ! --------- ! + ! Constants ! + ! --------- ! + + ! Wang and Sassen IWC paramters ( Option.1 ) + a = 26.87_kind_phys + b = 0.569_kind_phys + c = 0.002892_kind_phys + ! Schiller parameters ( Option.2 ) + as = -68.4202_kind_phys + bs = 0.983917_kind_phys + cs = 2.81795_kind_phys + ! Wood and Field parameters ( Option.3 ) + Kc = 75._kind_phys + ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) + ! Slingo modified (option 5) + minice = 1.e-12_kind_phys + mincld = 1.e-4_kind_phys + + if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys + + ! ---------------- ! + ! Main computation ! + ! ---------------- ! + + aist_out(:) = 0._kind_phys + esat_in(:) = 0._kind_phys + qsat_in(:) = 0._kind_phys + + call qsat_water(T_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol) + call svp_water_vect(T_in(1:ncol), esl(1:ncol), ncol) + call svp_ice_vect(T_in(1:ncol), esi(1:ncol), ncol) + + do i = 1, ncol + + landfrac = landfrac_in(i) + snowh = snowh_in(i) + T = T_in(i) + qv = qv_in(i) + p = p_in(i) + qi = qi_in(i) + ni = ni_in(i) + qs = qsat_in(i) + + rhmaxi = rhmaxi_in(i) + rhmini = rhmini_in(i) + rhminl = rhminl_in(i) + rhminl_adj_land = rhminl_adj_land_in(i) + rhminh = rhminh_in(i) + + if (iceopt < 3) then + if (iceopt == 1) then + ttmp = max(195._kind_phys, min(T, 253._kind_phys)) - 273.16_kind_phys + icicval = a + b*ttmp + c*ttmp**2._kind_phys + rho = p/(rair*T) + icicval = icicval*1.e-6_kind_phys/rho + else + ttmp = max(190._kind_phys, min(T, 273.16_kind_phys)) + icicval = 10._kind_phys**(as*bs**ttmp + cs) + icicval = icicval*1.e-6_kind_phys*18._kind_phys/28.97_kind_phys + end if + aist = max(0._kind_phys, min(qi/icicval, 1._kind_phys)) + elseif (iceopt == 3) then + aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi(i)/esl(i)))) + aist = max(0._kind_phys, min(aist, 1._kind_phys)) + elseif (iceopt == 4) then + if (p >= premib) then + if (land(i) .and. (snowh <= 0.000001_kind_phys)) then + rhmin = rhminl - rhminl_adj_land + else + rhmin = rhminl + end if + elseif (p < premit) then + rhmin = rhminh + else + rhwght = (premib - (max(p, premit)))/(premib - premit) + rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) + end if + ncf = qi/((1._kind_phys - icecrit)*qs) + if (ncf <= 0._kind_phys) then + aist = 0._kind_phys + elseif (ncf > 0._kind_phys .and. ncf <= 1._kind_phys/6._kind_phys) then + aist = 0.5_kind_phys*(6._kind_phys*ncf)**(2._kind_phys/3._kind_phys) + elseif (ncf > 1._kind_phys/6._kind_phys .and. ncf < 1._kind_phys) then + phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+4._kind_phys*3.1415927_kind_phys)/3._kind_phys + aist = (1._kind_phys - 4._kind_phys*cos(phi)*cos(phi)) + else + aist = 1._kind_phys + end if + aist = max(0._kind_phys, min(aist, 1._kind_phys)) + elseif (iceopt == 5) then ! set rh ice cloud fraction - rhi= (qv+qi)/qs * (esl/esi) - if (rhmaxi .eq. rhmini) then - if (rhi .gt. rhmini) then - rhdif = 1._kind_phys - else - rhdif = 0._kind_phys - end if + rhi = (qv + qi)/qs*(esl(i)/esi(i)) + if (rhmaxi == rhmini) then + if (rhi > rhmini) then + rhdif = 1._kind_phys + else + rhdif = 0._kind_phys + end if else - rhdif = (rhi-rhmini) / (rhmaxi - rhmini) + rhdif = (rhi - rhmini)/(rhmaxi - rhmini) end if - aist = min(1.0_kind_phys, max(rhdif,0._kind_phys)**2) + aist = min(1.0_kind_phys, max(rhdif, 0._kind_phys)**2) + + elseif (iceopt == 6) then + !----- ICE CLOUD OPTION 6: fit based on T and Number (Gettelman: based on Heymsfield obs) + ! Use observations from Heymsfield et al 2012 of IWC and Ni v. Temp + ! Multivariate fit follows form of Boudala 2002: ICIWC = a * exp(b*T) * N^c + ! a=6.73e-8, b=0.05, c=0.349 + ! N is #/L, so need to convert Ni_L=N*rhoa/1000. + ah = 6.73834e-08_kind_phys + bh = 0.0533110_kind_phys + ch = 0.3493813_kind_phys + rho = p/(rair*T) + nil = ni*rho/1000._kind_phys + icicval = ah*exp(bh*T)*nil**ch + !result is in g m-3, convert to kg H2O / kg air (icimr...) + icicval = icicval/rho/1000._kind_phys + aist = max(0._kind_phys, min(qi/icicval, 1._kind_phys)) + aist = min(aist, 1._kind_phys) + + end if + + if (iceopt == 5 .or. iceopt == 6) then ! Similar to alpha in Wilson & Ballard (1999), determine a ! scaling factor for saturation vapor pressure that reflects @@ -781,310 +986,49 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & ! ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less. if (present(qsatfac_out) .and. do_subgrid_growth) then - qsatfac_out = max(min(qv / qs, 1._kind_phys), (1._kind_phys - aist) * rhmini + aist * rhmaxi) + qsatfac_out(i) = max(min(qv/qs, 1._kind_phys), (1._kind_phys - aist)*rhmini + aist*rhmaxi) end if ! limiter to remove empty cloud and ice with no cloud ! and set icecld fraction to mincld if ice exists - if (qi.lt.minice) then - aist=0._kind_phys + if (qi < minice) then + aist = 0._kind_phys else - aist=max(mincld,aist) - endif + aist = max(mincld, aist) + end if ! enforce limits on icimr - if (qi.ge.minice) then - icimr=qi/aist - - !minimum - if (icimr.lt.qist_min) then - if (do_avg_aist_algs) then - ! - ! Take the geometric mean of the iceopt=4 and iceopt=5 values. - ! Mods developed by Thomas Toniazzo for NorESM. - aist = max(0._kind_phys,min(1._kind_phys,sqrt(aist*qi/qist_min))) - else - ! - ! Default for iceopt=5 - aist = max(0._kind_phys,min(1._kind_phys,qi/qist_min)) - end if - endif - !maximum - if (icimr.gt.qist_max) then - aist = max(0._kind_phys,min(1._kind_phys,qi/qist_max)) - endif - - endif - endif - - ! 0.999_kind_phys is added to prevent infinite 'ql_st' at the end of instratus_condensate - ! computed after updating 'qi_st'. - - aist = max(0._kind_phys,min(aist,0.999_kind_phys)) - -end subroutine aist_single + if (qi >= minice) then + icimr = qi/aist + + !minimum + if (icimr < qist_min) then + if (do_avg_aist_algs) then + ! + ! Take the geometric mean of the iceopt=4 and iceopt=5 values. + ! Mods developed by Thomas Toniazzo for NorESM. + aist = max(0._kind_phys, min(1._kind_phys, sqrt(aist*qi/qist_min))) + else + ! + ! Default for iceopt=5 + aist = max(0._kind_phys, min(1._kind_phys, qi/qist_min)) + end if + end if + !maximum + if (icimr > qist_max) then + aist = max(0._kind_phys, min(1._kind_phys, qi/qist_max)) + end if -!================================================================================================ + end if + end if -subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, aist_out, ncol, & - rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, & - qsatfac_out ) - - ! --------------------------------------------------------- ! - ! Compute non-physical ice stratus fraction ! - ! --------------------------------------------------------- ! - - use wv_saturation, only: qsat_water, svp_water_vect, svp_ice_vect - - real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor[kg/kg] - real(kind_phys), intent(in) :: T_in(:) ! Temperature - real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] - real(kind_phys), intent(in) :: qi_in(:) ! Grid-mean ice water content [kg/kg] - real(kind_phys), intent(in) :: ni_in(:) ! Grid-mean ice water number concentration [#/kg] - real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction - real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: aist_out(:) ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) - integer, intent(in) :: ncol - - real(kind_phys), intent(in) :: rhmaxi_in(:) - real(kind_phys), intent(in) :: rhmini_in(:) ! Critical relative humidity for ice stratus - real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus - real(kind_phys), optional, intent(out) :: qsatfac_out(:) ! Subgrid scaling factor for qsat - - ! Local variables - - real(kind_phys) qv ! Grid-mean water vapor[kg/kg] - real(kind_phys) T ! Temperature - real(kind_phys) p ! Pressure [Pa] - real(kind_phys) qi ! Grid-mean ice water content [kg/kg] - real(kind_phys) ni - real(kind_phys) landfrac ! Land fraction - real(kind_phys) snowh ! Snow depth (liquid water equivalent) - - real(kind_phys) rhmaxi ! Critical relative humidity for ice stratus - real(kind_phys) rhmini ! Critical relative humidity for ice stratus - real(kind_phys) rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys) rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys) rhminh ! Critical relative humidity for high-level liquid stratus - - real(kind_phys) aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) - - real(kind_phys) rhmin ! Critical RH - real(kind_phys) rhwght - - real(kind_phys) a,b,c,as,bs,cs,ah,bh,ch ! Fit parameters - real(kind_phys) nil - real(kind_phys) Kc ! Constant for ice cloud calc (wood & field) - real(kind_phys) ttmp ! Limited temperature - real(kind_phys) icicval ! Empirical IWC value [ kg/kg ] - real(kind_phys) rho ! Local air density - real(kind_phys) esl(ncol) ! Liq sat vapor pressure - real(kind_phys) esi(ncol) ! Ice sat vapor pressure - real(kind_phys) ncf,phi ! Wilson and Ballard parameters - real(kind_phys) qs - real(kind_phys) esat_in(ncol) - real(kind_phys) qsat_in(ncol) - - real(kind_phys) rhi ! grid box averaged relative humidity over ice - real(kind_phys) minice ! minimum grid box avg ice for having a 'cloud' - real(kind_phys) mincld ! minimum ice cloud fraction threshold - real(kind_phys) icimr ! in cloud ice mixing ratio - real(kind_phys) rhdif ! working variable for slingo scheme - - integer :: i - - - ! Statement functions - logical land - land(i) = nint(landfrac_in(i)) == 1 - - ! --------- ! - ! Constants ! - ! --------- ! - - ! Wang and Sassen IWC paramters ( Option.1 ) - a = 26.87_kind_phys - b = 0.569_kind_phys - c = 0.002892_kind_phys - ! Schiller parameters ( Option.2 ) - as = -68.4202_kind_phys - bs = 0.983917_kind_phys - cs = 2.81795_kind_phys - ! Wood and Field parameters ( Option.3 ) - Kc = 75._kind_phys - ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) - ! Slingo modified (option 5) - minice = 1.e-12_kind_phys - mincld = 1.e-4_kind_phys - - if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys - - ! ---------------- ! - ! Main computation ! - ! ---------------- ! - - aist_out(:) = 0._kind_phys - esat_in(:) = 0._kind_phys - qsat_in(:) = 0._kind_phys - - call qsat_water(T_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol) - call svp_water_vect(T_in(1:ncol), esl(1:ncol), ncol) - call svp_ice_vect(T_in(1:ncol), esi(1:ncol), ncol) - - do i = 1, ncol - - landfrac = landfrac_in(i) - snowh = snowh_in(i) - T = T_in(i) - qv = qv_in(i) - p = p_in(i) - qi = qi_in(i) - ni = ni_in(i) - qs = qsat_in(i) - - rhmaxi = rhmaxi_in(i) - rhmini = rhmini_in(i) - rhminl = rhminl_in(i) - rhminl_adj_land = rhminl_adj_land_in(i) - rhminh = rhminh_in(i) - - if( iceopt.lt.3 ) then - if( iceopt.eq.1 ) then - ttmp = max(195._kind_phys,min(T,253._kind_phys)) - 273.16_kind_phys - icicval = a + b * ttmp + c * ttmp**2._kind_phys - rho = p/(rair*T) - icicval = icicval * 1.e-6_kind_phys / rho - else - ttmp = max(190._kind_phys,min(T,273.16_kind_phys)) - icicval = 10._kind_phys **(as * bs**ttmp + cs) - icicval = icicval * 1.e-6_kind_phys * 18._kind_phys / 28.97_kind_phys - endif - aist = max(0._kind_phys,min(qi/icicval,1._kind_phys)) - elseif( iceopt.eq.3 ) then - aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi(i)/esl(i)))) - aist = max(0._kind_phys,min(aist,1._kind_phys)) - elseif( iceopt.eq.4) then - if( p .ge. premib ) then - if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then - rhmin = rhminl - rhminl_adj_land - else - rhmin = rhminl - endif - elseif( p .lt. premit ) then - rhmin = rhminh - else - rhwght = (premib-(max(p,premit)))/(premib-premit) - ! if( land(i) .and. (snowh.le.0.000001_kind_phys) ) then - ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys-rhwght) - ! else - rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys-rhwght) - ! endif - endif - ncf = qi/((1._kind_phys - icecrit)*qs) - if( ncf.le.0._kind_phys ) then - aist = 0._kind_phys - elseif( ncf.gt.0._kind_phys .and. ncf.le.1._kind_phys/6._kind_phys ) then - aist = 0.5_kind_phys*(6._kind_phys * ncf)**(2._kind_phys/3._kind_phys) - elseif( ncf.gt.1._kind_phys/6._kind_phys .and. ncf.lt.1._kind_phys ) then - phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+4._kind_phys*3.1415927_kind_phys)/3._kind_phys - aist = (1._kind_phys - 4._kind_phys * cos(phi) * cos(phi)) - else - aist = 1._kind_phys - endif - aist = max(0._kind_phys,min(aist,1._kind_phys)) - elseif (iceopt.eq.5) then - ! set rh ice cloud fraction - rhi= (qv+qi)/qs * (esl(i)/esi(i)) - if (rhmaxi .eq. rhmini) then - if (rhi .gt. rhmini) then - rhdif = 1._kind_phys - else - rhdif = 0._kind_phys - end if - else - rhdif = (rhi-rhmini) / (rhmaxi - rhmini) - end if - aist = min(1.0_kind_phys, max(rhdif,0._kind_phys)**2) - - elseif (iceopt.eq.6) then - !----- ICE CLOUD OPTION 6: fit based on T and Number (Gettelman: based on Heymsfield obs) - ! Use observations from Heymsfield et al 2012 of IWC and Ni v. Temp - ! Multivariate fit follows form of Boudala 2002: ICIWC = a * exp(b*T) * N^c - ! a=6.73e-8, b=0.05, c=0.349 - ! N is #/L, so need to convert Ni_L=N*rhoa/1000. - ah= 6.73834e-08_kind_phys - bh= 0.0533110_kind_phys - ch= 0.3493813_kind_phys - rho=p/(rair*T) - nil=ni*rho/1000._kind_phys - icicval = ah * exp(bh*T) * nil**ch - !result is in g m-3, convert to kg H2O / kg air (icimr...) - icicval = icicval / rho / 1000._kind_phys - aist = max(0._kind_phys,min(qi/icicval,1._kind_phys)) - aist = min(aist,1._kind_phys) - - endif - - if (iceopt.eq.5 .or. iceopt.eq.6) then - - ! Similar to alpha in Wilson & Ballard (1999), determine a - ! scaling factor for saturation vapor pressure that reflects - ! the cloud fraction, rhmini, and rhmaxi. - ! - ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less. - if (present(qsatfac_out) .and. do_subgrid_growth) then - qsatfac_out(i) = max(min(qv / qs, 1._kind_phys), (1._kind_phys - aist) * rhmini + aist * rhmaxi) - end if - - ! limiter to remove empty cloud and ice with no cloud - ! and set icecld fraction to mincld if ice exists - - if (qi.lt.minice) then - aist=0._kind_phys - else - aist=max(mincld,aist) - endif - - ! enforce limits on icimr - if (qi.ge.minice) then - icimr=qi/aist - - !minimum - if (icimr.lt.qist_min) then - if (do_avg_aist_algs) then - ! - ! Take the geometric mean of the iceopt=4 and iceopt=5 values. - ! Mods developed by Thomas Toniazzo for NorESM. - aist = max(0._kind_phys,min(1._kind_phys,sqrt(aist*qi/qist_min))) - else - ! - ! Default for iceopt=5 - aist = max(0._kind_phys,min(1._kind_phys,qi/qist_min)) - end if - endif - !maximum - if (icimr.gt.qist_max) then - aist = max(0._kind_phys,min(1._kind_phys,qi/qist_max)) - endif - - endif - endif - - ! 0.999_kind_phys is added to prevent infinite 'ql_st' at the end of instratus_condensate - ! computed after updating 'qi_st'. - - aist = max(0._kind_phys,min(aist,0.999_kind_phys)) - - aist_out(i) = aist - - enddo - -end subroutine aist_vector + ! 0.999_kind_phys is added to prevent infinite 'ql_st' at the end of instratus_condensate + ! computed after updating 'qi_st'. + aist = max(0._kind_phys, min(aist, 0.999_kind_phys)) + aist_out(i) = aist + end do -!================================================================================================ + end subroutine aist_vector end module compute_cloud_fraction_two_moment From bd71c062c1cfc3a9367c05a25111933960afc690 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 12 Jun 2026 14:42:42 -0400 Subject: [PATCH 4/9] Continue to fix formatting (take 2) --- .../compute_cloud_fraction_two_moment.F90 | 516 +++++++++--------- 1 file changed, 246 insertions(+), 270 deletions(-) diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 index f0fad5db..8305de2b 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -5,7 +5,6 @@ ! Provides liquid stratus (PDF and RHU methods) and ice stratus cloud fraction. ! Original authors: Sungsu Park, Andrew Gettelman module compute_cloud_fraction_two_moment - use ccpp_kinds, only: kind_phys implicit none @@ -28,29 +27,29 @@ module compute_cloud_fraction_two_moment !REMOVECAM_END ! Hardcoded parameters - logical, parameter :: CAMstfrac = .false. ! If .true. (.false.), + logical, parameter :: CAMstfrac = .false. ! If .true. (.false.), ! use Slingo (triangular PDF-based) liquid stratus fraction - logical, parameter :: freeze_dry = .false. ! If .true., use 'freeze dry' in liquid stratus fraction formula + logical, parameter :: freeze_dry = .false. ! If .true., use 'freeze dry' in liquid stratus fraction formula ! Module variables set by init (used internally by subroutines) - real(kind_phys) :: premit ! Top pressure bound for mid-level clouds [Pa] - real(kind_phys) :: premib ! Bottom pressure bound for mid-level clouds [Pa] - integer :: iceopt ! Ice cloud closure option + real(kind_phys) :: premit ! Top pressure bound for mid-level clouds [Pa] + real(kind_phys) :: premib ! Bottom pressure bound for mid-level clouds [Pa] + + ! Ice cloud closure option ! 1=wang & sassen 2=schiller (iciwc) ! 3=wood & field, 4=Wilson (based on smith) ! 5=modified slingo (ssat & empty cloud) - real(kind_phys) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure - ! ( smaller = more ice clouds ) - real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [ kg kg-1 ] - real(kind_phys) :: qist_max ! Maximum in-stratus IWC constraint [ kg kg-1 ] - logical :: do_subgrid_growth - logical :: do_avg_aist_algs - real(kind_phys) :: rair ! Gas constant of dry air [J kg-1 K-1] + integer :: iceopt -contains - -!================================================================================================ + real(kind_phys) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure + ! (smaller = more ice clouds) + real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [kg kg-1] + real(kind_phys) :: qist_max ! Maximum in-stratus IWC constraint [kg kg-1] + logical :: do_subgrid_growth + logical :: do_avg_aist_algs + real(kind_phys) :: rair ! Gas constant of dry air [J kg-1 K-1] +contains !> \section arg_table_compute_cloud_fraction_two_moment_init Argument Table !! \htmlinclude compute_cloud_fraction_two_moment_init.html subroutine compute_cloud_fraction_two_moment_init( & @@ -61,21 +60,21 @@ subroutine compute_cloud_fraction_two_moment_init( & errmsg, errflg) ! Input arguments - logical, intent(in) :: amIRoot - integer, intent(in) :: iulog - real(kind_phys), intent(in) :: premit_in - real(kind_phys), intent(in) :: premib_in - integer, intent(in) :: iceopt_in - real(kind_phys), intent(in) :: icecrit_in - real(kind_phys), intent(in) :: qist_min_in - real(kind_phys), intent(in) :: qist_max_in - logical, intent(in) :: do_subgrid_growth_in - logical, intent(in) :: do_avg_aist_algs_in - real(kind_phys), intent(in) :: rair_in + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + real(kind_phys), intent(in) :: premit_in + real(kind_phys), intent(in) :: premib_in + integer, intent(in) :: iceopt_in + real(kind_phys), intent(in) :: icecrit_in + real(kind_phys), intent(in) :: qist_min_in + real(kind_phys), intent(in) :: qist_max_in + logical, intent(in) :: do_subgrid_growth_in + logical, intent(in) :: do_avg_aist_algs_in + real(kind_phys), intent(in) :: rair_in ! Output arguments character(len=512), intent(out) :: errmsg - integer, intent(out) :: errflg + integer, intent(out) :: errflg errflg = 0 errmsg = '' @@ -105,61 +104,52 @@ subroutine compute_cloud_fraction_two_moment_init( & end subroutine compute_cloud_fraction_two_moment_init -!================================================================================================ - subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) - - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! analytical formulation of triangular PDF. ! - ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! - ! so using constant 'dV' assume that width is proportional ! - ! to the saturation specific humidity. ! - ! dV ~ 0.1. ! - ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. In fact, it does not ! - ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! - ! they will produce the same results. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U ! Relative humidity - real(kind_phys), intent(in) :: p ! Pressure [Pa] + ! --------------------------------------------------------- + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the + ! analytical formulation of triangular PDF. + ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', + ! so using constant 'dV' assume that width is proportional + ! to the saturation specific humidity. + ! dV ~ 0.1. + ! cldrh : RH of in-stratus( = 1 if no supersaturation) + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is + ! G is discontinuous across U = 1. In fact, it does not + ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that + ! they will produce the same results. + ! --------------------------------------------------------- + + real(kind_phys), intent(in) :: U ! Relative humidity + real(kind_phys), intent(in) :: p ! Pressure [Pa] real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] - real(kind_phys), intent(in) :: landfrac ! Land fraction - real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) - real(kind_phys), intent(out) :: a ! Stratus fraction + real(kind_phys), intent(out) :: a ! Stratus fraction real(kind_phys), intent(out) :: Ga ! dU/da - real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus - real(kind_phys), optional, intent(out) :: orhmin ! Critical RH + real(kind_phys), optional, intent(out) :: orhmin ! Critical RH ! Local variables - real(kind_phys) :: dV ! Width of triangular PDF - real(kind_phys) :: cldrh ! RH of stratus cloud - real(kind_phys) :: rhmin ! Critical RH + real(kind_phys), parameter :: cldrh = 1.0_kind_phys ! RH of stratus cloud + + real(kind_phys) :: dV ! Width of triangular PDF + real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght ! Statement functions logical land land = nint(landfrac) == 1 - ! ---------- ! - ! Parameters ! - ! ---------- ! - - cldrh = 1.0_kind_phys - ! ---------------- ! ! Main computation ! ! ---------------- ! if (p >= premib) then - if (land .and. (snowh <= 0.000001_kind_phys)) then rhmin = rhminl - rhminl_adj_land else @@ -175,8 +165,9 @@ subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*& + (acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) elseif (U <= (cldrh - dV)) then a = 0._kind_phys @@ -187,9 +178,7 @@ subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ a = a*max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) Ga = Ga/max(0.15_kind_phys, min(1.0_kind_phys, qv/0.0030_kind_phys)) end if - elseif (p < premit) then - rhmin = rhminh dV = cldrh - rhmin @@ -200,16 +189,15 @@ subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*& + (acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) elseif (U <= (cldrh - dV)) then a = 0._kind_phys Ga = 1.e10_kind_phys end if - else - rhwght = (premib - (max(p, premit)))/(premib - premit) rhmin = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) @@ -222,83 +210,74 @@ subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*& + (acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) elseif (U <= (cldrh - dV)) then a = 0._kind_phys Ga = 1.e10_kind_phys end if - end if if (present(orhmin)) orhmin = rhmin end subroutine astG_PDF_single -!================================================================================================ - subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & rhminl_in, rhminl_adj_land_in, rhminh_in) + ! --------------------------------------------------------- + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the + ! analytical formulation of triangular PDF. + ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', + ! so using constant 'dV' assume that width is proportional + ! to the saturation specific humidity. + ! dV ~ 0.1. + ! cldrh : RH of in-stratus( = 1 if no supersaturation) + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is + ! G is discontinuous across U = 1. In fact, it does not + ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that + ! they will produce the same results. + ! --------------------------------------------------------- + + real(kind_phys), intent(in) :: U_in(:) ! Relative humidity + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction + real(kind_phys), intent(out) :: Ga_out(:) ! dU/da + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! analytical formulation of triangular PDF. ! - ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', ! - ! so using constant 'dV' assume that width is proportional ! - ! to the saturation specific humidity. ! - ! dV ~ 0.1. ! - ! cldrh : RH of in-stratus( = 1 if no supersaturation) ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. In fact, it does not ! - ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that ! - ! they will produce the same results. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U_in(:) ! Relative humidity - real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] - real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg kg-1] - real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction - real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction - real(kind_phys), intent(out) :: Ga_out(:) ! dU/da - integer, intent(in) :: ncol - - real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus - real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus - - real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus - - real(kind_phys) :: U ! Relative humidity - real(kind_phys) :: p ! Pressure [Pa] - real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] - real(kind_phys) :: landfrac ! Land fraction - real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) - - real(kind_phys) :: a ! Stratus fraction - real(kind_phys) :: Ga ! dU/da + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + + real(kind_phys) :: U ! Relative humidity + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys) :: a ! Stratus fraction + real(kind_phys) :: Ga ! dU/da ! Local variables - integer :: i ! Loop indexes - real(kind_phys) :: dV ! Width of triangular PDF - real(kind_phys) :: cldrh ! RH of stratus cloud - real(kind_phys) :: rhmin ! Critical RH + real(kind_phys), parameter :: cldrh = 1.0_kind_phys ! RH of stratus cloud + integer :: i ! Loop indexes + real(kind_phys) :: dV ! Width of triangular PDF + real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght ! Statement functions logical land land(i) = nint(landfrac_in(i)) == 1 - ! ---------- ! - ! Parameters ! - ! ---------- ! - - cldrh = 1.0_kind_phys - ! ---------------- ! ! Main computation ! ! ---------------- ! @@ -333,8 +312,9 @@ subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, nco a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*& + (acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) elseif (U <= (cldrh - dV)) then a = 0._kind_phys @@ -356,8 +336,9 @@ subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, nco a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*& + (acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) elseif (U <= (cldrh - dV)) then a = 0._kind_phys @@ -375,8 +356,9 @@ subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, nco a = 1._kind_phys - (-3._kind_phys/sqrt(2._kind_phys)*(U - cldrh)/dV)**(2._kind_phys/3._kind_phys) Ga = dV/sqrt(2._kind_phys)*sqrt(1._kind_phys - a) elseif (U > (cldrh - dV) .and. U <= (cldrh - dV/6._kind_phys)) then - a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*(acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & - (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys + a = 4._kind_phys*(cos((1._kind_phys/3._kind_phys)*& + (acos((3._kind_phys/2._kind_phys/sqrt(2._kind_phys))* & + (1._kind_phys + (U - cldrh)/dV)) - 2._kind_phys*3.141592_kind_phys)))**2._kind_phys Ga = dV/sqrt(2._kind_phys)*(1._kind_phys/sqrt(a) - sqrt(a)) elseif (U <= (cldrh - dV)) then a = 0._kind_phys @@ -389,32 +371,31 @@ subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, nco end do end subroutine astG_PDF -!================================================================================================ subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! CAM35 cloud fraction formula. ! - ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! - ! For the other cases, I should re-define 'rhminl,rhminh' & ! - ! 'premib,premit'. ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U ! Relative humidity - real(kind_phys), intent(in) :: p ! Pressure [Pa] - real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] - real(kind_phys), intent(in) :: landfrac ! Land fraction - real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: a ! Stratus fraction - real(kind_phys), intent(out) :: Ga ! dU/da - - real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + ! --------------------------------------------------------- + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the + ! CAM35 cloud fraction formula. + ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core + ! For the other cases, I should re-define 'rhminl,rhminh' & + ! 'premib,premit'. + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is + ! G is discontinuous across U = 1. + ! --------------------------------------------------------- + + real(kind_phys), intent(in) :: U ! Relative humidity + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a ! Stratus fraction + real(kind_phys), intent(out) :: Ga ! dU/da + + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus real(kind_phys), optional, intent(out) :: orhmin ! Critical RH @@ -476,53 +457,51 @@ subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ end subroutine astG_RHU_single -!================================================================================================ - subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & rhminl_in, rhminl_adj_land_in, rhminh_in) - ! --------------------------------------------------------- ! - ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! - ! CAM3.5 cloud fraction formula. ! - ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! - ! For the other cases, I should re-define 'rhminl,rhminh' & ! - ! 'premib,premit'. ! - ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is ! - ! G is discontinuous across U = 1. ! - ! --------------------------------------------------------- ! - - real(kind_phys), intent(in) :: U_in(:) ! Relative humidity - real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] - real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg kg-1] - real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction - real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) - - real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction - real(kind_phys), intent(out) :: Ga_out(:) ! dU/da - integer, intent(in) :: ncol - - real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + ! --------------------------------------------------------- + ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the + ! CAM3.5 cloud fraction formula. + ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core + ! For the other cases, I should re-define 'rhminl,rhminh' & + ! 'premib,premit'. + ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is + ! G is discontinuous across U = 1. + ! --------------------------------------------------------- + + real(kind_phys), intent(in) :: U_in(:) ! Relative humidity + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + + real(kind_phys), intent(out) :: a_out(:) ! Stratus fraction + real(kind_phys), intent(out) :: Ga_out(:) ! dU/da + integer, intent(in) :: ncol + + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus - real(kind_phys) :: U ! Relative humidity - real(kind_phys) :: p ! Pressure [Pa] - real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] - real(kind_phys) :: landfrac ! Land fraction - real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + real(kind_phys) :: U ! Relative humidity + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qv ! Grid-mean water vapor specific humidity [kg kg-1] + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) - real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus - real(kind_phys) :: a ! Stratus fraction - real(kind_phys) :: Ga ! dU/da + real(kind_phys) :: a ! Stratus fraction + real(kind_phys) :: Ga ! dU/da ! Local variables - integer i - real(kind_phys) rhmin ! Critical RH - real(kind_phys) rhdif ! Factor for stratus fraction - real(kind_phys) rhwght + integer :: i + real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhdif ! Factor for stratus fraction + real(kind_phys) :: rhwght ! Statement functions logical land @@ -591,53 +570,51 @@ subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, nco end subroutine astG_RHU -!================================================================================================ - subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh, & qsatfac_out) - ! --------------------------------------------------------- ! - ! Compute non-physical ice stratus fraction ! - ! --------------------------------------------------------- ! + ! ------------------------------------------ + ! Compute non-physical ice stratus fraction + ! ------------------------------------------ use wv_saturation, only: qsat_water, svp_water, svp_ice - real(kind_phys), intent(in) :: qv ! Grid-mean water vapor[kg kg-1] - real(kind_phys), intent(in) :: T ! Temperature - real(kind_phys), intent(in) :: p ! Pressure [Pa] - real(kind_phys), intent(in) :: qi ! Grid-mean ice water content [kg kg-1] - real(kind_phys), intent(in) :: landfrac ! Land fraction - real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) + real(kind_phys), intent(in) :: qv ! Grid-mean water vapor[kg kg-1] + real(kind_phys), intent(in) :: T ! Temperature + real(kind_phys), intent(in) :: p ! Pressure [Pa] + real(kind_phys), intent(in) :: qi ! Grid-mean ice water content [kg kg-1] + real(kind_phys), intent(in) :: landfrac ! Land fraction + real(kind_phys), intent(in) :: snowh ! Snow depth (liquid water equivalent) - real(kind_phys), intent(out) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + real(kind_phys), intent(out) :: aist ! Non-physical ice stratus fraction (0<=aist<=1) real(kind_phys), intent(in) :: rhmaxi - real(kind_phys), intent(in) :: rhmini ! Critical relative humidity for ice stratus - real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhmini ! Critical relative humidity for ice stratus + real(kind_phys), intent(in) :: rhminl ! Critical relative humidity for low-level liquid stratus real(kind_phys), intent(in) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus - real(kind_phys), optional, intent(out) :: qsatfac_out ! Subgrid scaling factor for qsat + real(kind_phys), intent(in) :: rhminh ! Critical relative humidity for high-level liquid stratus + real(kind_phys), optional, intent(out) :: qsatfac_out ! Subgrid scaling factor for qsat ! Local variables - real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght - real(kind_phys) :: a, b, c, as, bs, cs ! Fit parameters - real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) - real(kind_phys) :: ttmp ! Limited temperature - real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] - real(kind_phys) :: rho ! Local air density - real(kind_phys) :: esl ! Liq sat vapor pressure - real(kind_phys) :: esi ! Ice sat vapor pressure - real(kind_phys) :: ncf, phi ! Wilson and Ballard parameters + real(kind_phys) :: a, b, c, as, bs, cs ! Fit parameters + real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) + real(kind_phys) :: ttmp ! Limited temperature + real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] + real(kind_phys) :: rho ! Local air density + real(kind_phys) :: esl ! Liq sat vapor pressure + real(kind_phys) :: esi ! Ice sat vapor pressure + real(kind_phys) :: ncf, phi ! Wilson and Ballard parameters real(kind_phys) :: es, qs - real(kind_phys) :: rhi ! grid box averaged relative humidity over ice - real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' - real(kind_phys) :: mincld ! minimum ice cloud fraction threshold - real(kind_phys) :: icimr ! in cloud ice mixing ratio - real(kind_phys) :: rhdif ! working variable for slingo scheme + real(kind_phys) :: rhi ! grid box averaged relative humidity over ice + real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' + real(kind_phys) :: mincld ! minimum ice cloud fraction threshold + real(kind_phys) :: icimr ! in cloud ice mixing ratio + real(kind_phys) :: rhdif ! working variable for slingo scheme ! Statement functions logical land @@ -710,7 +687,7 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & elseif (ncf > 0._kind_phys .and. ncf <= 1._kind_phys/6._kind_phys) then aist = 0.5_kind_phys*(6._kind_phys*ncf)**(2._kind_phys/3._kind_phys) elseif (ncf > 1._kind_phys/6._kind_phys .and. ncf < 1._kind_phys) then - phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+& + phi = (acos(3._kind_phys*(1._kind_phys - ncf)/2._kind_phys**(3._kind_phys/2._kind_phys)) + & 4._kind_phys*3.1415927_kind_phys)/3._kind_phys aist = (1._kind_phys - 4._kind_phys*cos(phi)*cos(phi)) else @@ -779,77 +756,75 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & end subroutine aist_single -!================================================================================================ - subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, aist_out, ncol, & rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, & qsatfac_out) - ! --------------------------------------------------------- ! - ! Compute non-physical ice stratus fraction ! - ! --------------------------------------------------------- ! + ! ------------------------------------------ + ! Compute non-physical ice stratus fraction + ! ------------------------------------------ use wv_saturation, only: qsat_water, svp_water_vect, svp_ice_vect - real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor[kg kg-1] - real(kind_phys), intent(in) :: T_in(:) ! Temperature - real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] - real(kind_phys), intent(in) :: qi_in(:) ! Grid-mean ice water content [kg kg-1] - real(kind_phys), intent(in) :: ni_in(:) ! Grid-mean ice water number concentration [#/kg] - real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction - real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) + real(kind_phys), intent(in) :: qv_in(:) ! Grid-mean water vapor [kg kg-1] + real(kind_phys), intent(in) :: T_in(:) ! Temperature + real(kind_phys), intent(in) :: p_in(:) ! Pressure [Pa] + real(kind_phys), intent(in) :: qi_in(:) ! Grid-mean ice water content [kg kg-1] + real(kind_phys), intent(in) :: ni_in(:) ! Grid-mean ice water number concentration [# kg-1] + real(kind_phys), intent(in) :: landfrac_in(:) ! Land fraction + real(kind_phys), intent(in) :: snowh_in(:) ! Snow depth (liquid water equivalent) - real(kind_phys), intent(out) :: aist_out(:) ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) - integer, intent(in) :: ncol + real(kind_phys), intent(out) :: aist_out(:) ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + integer, intent(in) :: ncol real(kind_phys), intent(in) :: rhmaxi_in(:) - real(kind_phys), intent(in) :: rhmini_in(:) ! Critical relative humidity for ice stratus - real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus + real(kind_phys), intent(in) :: rhmini_in(:) ! Critical relative humidity for ice stratus + real(kind_phys), intent(in) :: rhminl_in(:) ! Critical relative humidity for low-level liquid stratus real(kind_phys), intent(in) :: rhminl_adj_land_in(:) ! Adjustment drop of rhminl over the land - real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus - real(kind_phys), optional, intent(out) :: qsatfac_out(:) ! Subgrid scaling factor for qsat + real(kind_phys), intent(in) :: rhminh_in(:) ! Critical relative humidity for high-level liquid stratus + real(kind_phys), optional, intent(out) :: qsatfac_out(:) ! Subgrid scaling factor for qsat ! Local variables - real(kind_phys) :: qv ! Grid-mean water vapor [kg kg-1] - real(kind_phys) :: T ! Temperature - real(kind_phys) :: p ! Pressure [Pa] - real(kind_phys) :: qi ! Grid-mean ice water content [kg kg-1] + real(kind_phys) :: qv ! Grid-mean water vapor [kg kg-1] + real(kind_phys) :: T ! Temperature + real(kind_phys) :: p ! Pressure [Pa] + real(kind_phys) :: qi ! Grid-mean ice water content [kg kg-1] real(kind_phys) :: ni - real(kind_phys) :: landfrac ! Land fraction - real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) + real(kind_phys) :: landfrac ! Land fraction + real(kind_phys) :: snowh ! Snow depth (liquid water equivalent) - real(kind_phys) :: rhmaxi ! Critical relative humidity for ice stratus - real(kind_phys) :: rhmini ! Critical relative humidity for ice stratus - real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus - real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land - real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus + real(kind_phys) :: rhmaxi ! Critical relative humidity for ice stratus + real(kind_phys) :: rhmini ! Critical relative humidity for ice stratus + real(kind_phys) :: rhminl ! Critical relative humidity for low-level liquid stratus + real(kind_phys) :: rhminl_adj_land ! Adjustment drop of rhminl over the land + real(kind_phys) :: rhminh ! Critical relative humidity for high-level liquid stratus - real(kind_phys) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) + real(kind_phys) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) - real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght - real(kind_phys) :: a, b, c, as, bs, cs, ah, bh, ch! Fit parameters + real(kind_phys) :: a, b, c, as, bs, cs, ah, bh, ch ! Fit parameters real(kind_phys) :: nil - real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) - real(kind_phys) :: ttmp ! Limited temperature - real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] - real(kind_phys) :: rho ! Local air density - real(kind_phys) :: esl(ncol) ! Liq sat vapor pressure - real(kind_phys) :: esi(ncol) ! Ice sat vapor pressure - real(kind_phys) :: ncf, phi ! Wilson and Ballard parameters + real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) + real(kind_phys) :: ttmp ! Limited temperature + real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] + real(kind_phys) :: rho ! Local air density + real(kind_phys) :: esl(ncol) ! Liq sat vapor pressure + real(kind_phys) :: esi(ncol) ! Ice sat vapor pressure + real(kind_phys) :: ncf, phi ! Wilson and Ballard parameters real(kind_phys) :: qs real(kind_phys) :: esat_in(ncol) real(kind_phys) :: qsat_in(ncol) - real(kind_phys) :: rhi ! grid box averaged relative humidity over ice - real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' - real(kind_phys) :: mincld ! minimum ice cloud fraction threshold - real(kind_phys) :: icimr ! in cloud ice mixing ratio - real(kind_phys) :: rhdif ! working variable for slingo scheme + real(kind_phys) :: rhi ! grid box averaged relative humidity over ice + real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' + real(kind_phys) :: mincld ! minimum ice cloud fraction threshold + real(kind_phys) :: icimr ! in cloud ice mixing ratio + real(kind_phys) :: rhdif ! working variable for slingo scheme - integer :: i + integer :: i ! Statement functions logical land @@ -939,7 +914,8 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a elseif (ncf > 0._kind_phys .and. ncf <= 1._kind_phys/6._kind_phys) then aist = 0.5_kind_phys*(6._kind_phys*ncf)**(2._kind_phys/3._kind_phys) elseif (ncf > 1._kind_phys/6._kind_phys .and. ncf < 1._kind_phys) then - phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+4._kind_phys*3.1415927_kind_phys)/3._kind_phys + phi = (acos(3._kind_phys*(1._kind_phys-ncf)/2._kind_phys**(3._kind_phys/2._kind_phys))+& + 4._kind_phys*3.1415927_kind_phys)/3._kind_phys aist = (1._kind_phys - 4._kind_phys*cos(phi)*cos(phi)) else aist = 1._kind_phys From 3ea33c4112674b90494e03bee89f94e0710154c0 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 12 Jun 2026 14:54:35 -0400 Subject: [PATCH 5/9] Cleanup parameters; fix standard name. --- .../compute_cloud_fraction.meta | 2 +- .../compute_cloud_fraction_namelist.xml | 2 +- .../compute_cloud_fraction_two_moment.F90 | 129 +++++++----------- .../compute_cloud_fraction_two_moment.meta | 2 +- ...ute_cloud_fraction_two_moment_namelist.xml | 72 ++++++++++ 5 files changed, 126 insertions(+), 81 deletions(-) diff --git a/schemes/cloud_fraction/compute_cloud_fraction.meta b/schemes/cloud_fraction/compute_cloud_fraction.meta index c202af13..a5085f56 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction.meta +++ b/schemes/cloud_fraction/compute_cloud_fraction.meta @@ -79,7 +79,7 @@ dimensions = () intent = in [ premib_in ] - standard_name = tunable_parameter_for_bottom_pressure_bound_for_mid_level_liquid_stratus_for_cloud_fraction + standard_name = tunable_parameter_for_bottom_pressure_bound_for_mid_level_clouds_for_cloud_fraction units = Pa type = real | kind = kind_phys dimensions = () diff --git a/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml b/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml index f91074fe..c9794bce 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml +++ b/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml @@ -181,7 +181,7 @@ real cldfrc cldfrc_nl - tunable_parameter_for_bottom_pressure_bound_for_mid_level_liquid_stratus_for_cloud_fraction + tunable_parameter_for_bottom_pressure_bound_for_mid_level_clouds_for_cloud_fraction Pa Bottom pressure bound for mid-level liquid stratus fraction. diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 index 8305de2b..41c0cf59 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -36,11 +36,32 @@ module compute_cloud_fraction_two_moment real(kind_phys) :: premib ! Bottom pressure bound for mid-level clouds [Pa] ! Ice cloud closure option - ! 1=wang & sassen 2=schiller (iciwc) - ! 3=wood & field, 4=Wilson (based on smith) - ! 5=modified slingo (ssat & empty cloud) + ! 1 = Wang and Sassen + ! 2 = Schiller (iciwc) + ! 3 = Wood and Field + ! 4 = Wilson (based on Smith) + ! 5 = modified Slingo (ssat & empty cloud) integer :: iceopt + ! Ice cloud fraction fit parameters (used by aist_single and aist_vector) + ! Wang and Sassen IWC parameters (iceopt=1) + ! DOI: 10.1175/1520-0469(2002)059<2291:CCMPRU>2.0.CO;2 + real(kind_phys), parameter :: wang_sassen_a = 26.87_kind_phys + real(kind_phys), parameter :: wang_sassen_b = 0.569_kind_phys + real(kind_phys), parameter :: wang_sassen_c = 0.002892_kind_phys + ! Schiller parameters (iceopt=2) + ! DOI: 10.1029/2008JD010342 + real(kind_phys), parameter :: schiller_a = -68.4202_kind_phys + real(kind_phys), parameter :: schiller_b = 0.983917_kind_phys + real(kind_phys), parameter :: schiller_c = 2.81795_kind_phys + ! Wood and Field parameters (iceopt=3) + ! DOI: 10.1175/1520-0469(2000)057<1888:RBTWCW>2.0.CO;2 + real(kind_phys), parameter :: wood_field_Kc = 75._kind_phys + ! Minimum grid box avg ice for having a 'cloud' + real(kind_phys), parameter :: minice = 1.e-12_kind_phys + ! Minimum ice cloud fraction threshold + real(kind_phys), parameter :: mincld = 1.e-4_kind_phys + real(kind_phys) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure ! (smaller = more ice clouds) real(kind_phys) :: qist_min ! Minimum in-stratus IWC constraint [kg kg-1] @@ -104,7 +125,7 @@ subroutine compute_cloud_fraction_two_moment_init( & end subroutine compute_cloud_fraction_two_moment_init - subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) + pure subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) ! --------------------------------------------------------- ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! analytical formulation of triangular PDF. @@ -224,8 +245,8 @@ subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ end subroutine astG_PDF_single - subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & - rhminl_in, rhminl_adj_land_in, rhminh_in) + pure subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & + rhminl_in, rhminl_adj_land_in, rhminh_in) ! --------------------------------------------------------- ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the ! analytical formulation of triangular PDF. @@ -372,7 +393,7 @@ subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, nco end subroutine astG_PDF - subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) + pure subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) ! --------------------------------------------------------- ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the @@ -457,8 +478,8 @@ subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_ end subroutine astG_RHU_single - subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & - rhminl_in, rhminl_adj_land_in, rhminh_in) + pure subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, & + rhminl_in, rhminl_adj_land_in, rhminh_in) ! --------------------------------------------------------- ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the @@ -600,8 +621,6 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght - real(kind_phys) :: a, b, c, as, bs, cs ! Fit parameters - real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) real(kind_phys) :: ttmp ! Limited temperature real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] real(kind_phys) :: rho ! Local air density @@ -611,8 +630,6 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & real(kind_phys) :: es, qs real(kind_phys) :: rhi ! grid box averaged relative humidity over ice - real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' - real(kind_phys) :: mincld ! minimum ice cloud fraction threshold real(kind_phys) :: icimr ! in cloud ice mixing ratio real(kind_phys) :: rhdif ! working variable for slingo scheme @@ -620,29 +637,6 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & logical land land = nint(landfrac) == 1 - ! --------- ! - ! Constants ! - ! --------- ! - - ! Wang and Sassen IWC paramters ( Option.1 ) - ! DOI: 10.1175/1520-0469(2002)059<2291:CCMPRU>2.0.CO;2 - a = 26.87_kind_phys - b = 0.569_kind_phys - c = 0.002892_kind_phys - ! Schiller parameters ( Option.2 ) - ! DOI: 10.1029/2008JD010342 - as = -68.4202_kind_phys - bs = 0.983917_kind_phys - cs = 2.81795_kind_phys - ! Wood and Field parameters ( Option.3 ) - ! DOI: 10.1175/1520-0469(2000)057<1888:RBTWCW>2.0.CO;2 - Kc = 75._kind_phys - ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) - ! DOI: 10.1002/qj.49712555707 - ! Slingo modified (option 5) - minice = 1.e-12_kind_phys - mincld = 1.e-4_kind_phys - if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys ! ---------------- ! @@ -656,17 +650,17 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & if (iceopt < 3) then if (iceopt == 1) then ttmp = max(195._kind_phys, min(T, 253._kind_phys)) - 273.16_kind_phys - icicval = a + b*ttmp + c*ttmp**2._kind_phys + icicval = wang_sassen_a + wang_sassen_b*ttmp + wang_sassen_c*ttmp**2._kind_phys rho = p/(rair*T) icicval = icicval*1.e-6_kind_phys/rho else ttmp = max(190._kind_phys, min(T, 273.16_kind_phys)) - icicval = 10._kind_phys**(as*bs**ttmp + cs) + icicval = 10._kind_phys**(schiller_a*schiller_b**ttmp + schiller_c) icicval = icicval*1.e-6_kind_phys*18._kind_phys/28.97_kind_phys end if aist = max(0._kind_phys, min(qi/icicval, 1._kind_phys)) elseif (iceopt == 3) then - aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi/esl))) + aist = 1._kind_phys - exp(-wood_field_Kc*qi/(qs*(esi/esl))) aist = max(0._kind_phys, min(aist, 1._kind_phys)) elseif (iceopt == 4) then if (p >= premib) then @@ -805,9 +799,6 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght - real(kind_phys) :: a, b, c, as, bs, cs, ah, bh, ch ! Fit parameters - real(kind_phys) :: nil - real(kind_phys) :: Kc ! Constant for ice cloud calc (wood & field) real(kind_phys) :: ttmp ! Limited temperature real(kind_phys) :: icicval ! Empirical IWC value [ kg kg-1 ] real(kind_phys) :: rho ! Local air density @@ -819,10 +810,15 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a real(kind_phys) :: qsat_in(ncol) real(kind_phys) :: rhi ! grid box averaged relative humidity over ice - real(kind_phys) :: minice ! minimum grid box avg ice for having a 'cloud' - real(kind_phys) :: mincld ! minimum ice cloud fraction threshold real(kind_phys) :: icimr ! in cloud ice mixing ratio real(kind_phys) :: rhdif ! working variable for slingo scheme + real(kind_phys) :: nil ! Ice number concentration [#/L] + + ! Heymsfield/Boudala fit parameters (iceopt=6) + ! Boudala et al. (2002), Heymsfield et al. (2012) + real(kind_phys), parameter :: heymsfield_a = 6.73834e-08_kind_phys + real(kind_phys), parameter :: heymsfield_b = 0.0533110_kind_phys + real(kind_phys), parameter :: heymsfield_c = 0.3493813_kind_phys integer :: i @@ -830,25 +826,6 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a logical land land(i) = nint(landfrac_in(i)) == 1 - ! --------- ! - ! Constants ! - ! --------- ! - - ! Wang and Sassen IWC paramters ( Option.1 ) - a = 26.87_kind_phys - b = 0.569_kind_phys - c = 0.002892_kind_phys - ! Schiller parameters ( Option.2 ) - as = -68.4202_kind_phys - bs = 0.983917_kind_phys - cs = 2.81795_kind_phys - ! Wood and Field parameters ( Option.3 ) - Kc = 75._kind_phys - ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) - ! Slingo modified (option 5) - minice = 1.e-12_kind_phys - mincld = 1.e-4_kind_phys - if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys ! ---------------- ! @@ -883,17 +860,17 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a if (iceopt < 3) then if (iceopt == 1) then ttmp = max(195._kind_phys, min(T, 253._kind_phys)) - 273.16_kind_phys - icicval = a + b*ttmp + c*ttmp**2._kind_phys + icicval = wang_sassen_a + wang_sassen_b*ttmp + wang_sassen_c*ttmp**2._kind_phys rho = p/(rair*T) icicval = icicval*1.e-6_kind_phys/rho else ttmp = max(190._kind_phys, min(T, 273.16_kind_phys)) - icicval = 10._kind_phys**(as*bs**ttmp + cs) + icicval = 10._kind_phys**(schiller_a*schiller_b**ttmp + schiller_c) icicval = icicval*1.e-6_kind_phys*18._kind_phys/28.97_kind_phys end if aist = max(0._kind_phys, min(qi/icicval, 1._kind_phys)) elseif (iceopt == 3) then - aist = 1._kind_phys - exp(-Kc*qi/(qs*(esi(i)/esl(i)))) + aist = 1._kind_phys - exp(-wood_field_Kc*qi/(qs*(esi(i)/esl(i)))) aist = max(0._kind_phys, min(aist, 1._kind_phys)) elseif (iceopt == 4) then if (p >= premib) then @@ -937,17 +914,16 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a elseif (iceopt == 6) then !----- ICE CLOUD OPTION 6: fit based on T and Number (Gettelman: based on Heymsfield obs) - ! Use observations from Heymsfield et al 2012 of IWC and Ni v. Temp - ! Multivariate fit follows form of Boudala 2002: ICIWC = a * exp(b*T) * N^c + ! Use observations from + ! Heymsfield et al. 2013 (https://doi.org/10.1175/JAS-D-12-0124.1) of IWC and Ni v. Temp + ! Multivariate fit follows form of + ! Boudala 2002 (https://doi.org/10.1002/joc.774): ICIWC = a * exp(b*T) * N^c ! a=6.73e-8, b=0.05, c=0.349 ! N is #/L, so need to convert Ni_L=N*rhoa/1000. - ah = 6.73834e-08_kind_phys - bh = 0.0533110_kind_phys - ch = 0.3493813_kind_phys rho = p/(rair*T) nil = ni*rho/1000._kind_phys - icicval = ah*exp(bh*T)*nil**ch - !result is in g m-3, convert to kg H2O / kg air (icimr...) + icicval = heymsfield_a*exp(heymsfield_b*T)*nil**heymsfield_c + ! result is in g m-3, convert to kg H2O / kg air (icimr...) icicval = icicval/rho/1000._kind_phys aist = max(0._kind_phys, min(qi/icicval, 1._kind_phys)) aist = min(aist, 1._kind_phys) @@ -955,10 +931,10 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a end if if (iceopt == 5 .or. iceopt == 6) then - ! Similar to alpha in Wilson & Ballard (1999), determine a ! scaling factor for saturation vapor pressure that reflects ! the cloud fraction, rhmini, and rhmaxi. + ! https://doi.org/10.1002/qj.49712555707 ! ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less. if (present(qsatfac_out) .and. do_subgrid_growth) then @@ -967,7 +943,6 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a ! limiter to remove empty cloud and ice with no cloud ! and set icecld fraction to mincld if ice exists - if (qi < minice) then aist = 0._kind_phys else @@ -978,15 +953,13 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a if (qi >= minice) then icimr = qi/aist - !minimum + ! minimum: if (icimr < qist_min) then if (do_avg_aist_algs) then - ! ! Take the geometric mean of the iceopt=4 and iceopt=5 values. ! Mods developed by Thomas Toniazzo for NorESM. aist = max(0._kind_phys, min(1._kind_phys, sqrt(aist*qi/qist_min))) else - ! ! Default for iceopt=5 aist = max(0._kind_phys, min(1._kind_phys, qi/qist_min)) end if diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta index 23808fc3..20360008 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.meta @@ -25,7 +25,7 @@ dimensions = () intent = in [ premib_in ] - standard_name = tunable_parameter_for_bottom_pressure_bound_for_mid_level_liquid_stratus_for_cloud_fraction + standard_name = tunable_parameter_for_bottom_pressure_bound_for_mid_level_clouds_for_cloud_fraction units = Pa type = real | kind = kind_phys dimensions = () diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml b/schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml index 978fb7a7..90a9a0b1 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment_namelist.xml @@ -4,6 +4,78 @@ + + real cldfrc2m From f06b8ea8f07699e54f84f44c4fd7595f6799e296 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 12 Jun 2026 14:59:45 -0400 Subject: [PATCH 6/9] CAM35 -> CAM3.5 --- schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 index 41c0cf59..31260950 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -397,7 +397,7 @@ pure subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl ! --------------------------------------------------------- ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the - ! CAM35 cloud fraction formula. + ! CAM3.5 cloud fraction formula. ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core ! For the other cases, I should re-define 'rhminl,rhminh' & ! 'premib,premit'. From e110fce0ae2986645644c834a02d45d95d99f3e1 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 12 Jun 2026 15:30:46 -0400 Subject: [PATCH 7/9] Extract land statement function (impure) to pure elemental function. Assisted-by: claude-opus:4.6[1m] --- .../compute_cloud_fraction_two_moment.F90 | 36 +++++++------------ 1 file changed, 12 insertions(+), 24 deletions(-) diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 index 31260950..18966980 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -71,6 +71,12 @@ module compute_cloud_fraction_two_moment real(kind_phys) :: rair ! Gas constant of dry air [J kg-1 K-1] contains + + pure elemental logical function is_land(landfrac) + real(kind_phys), intent(in) :: landfrac + is_land = nint(landfrac) == 1 + end function is_land + !> \section arg_table_compute_cloud_fraction_two_moment_init Argument Table !! \htmlinclude compute_cloud_fraction_two_moment_init.html subroutine compute_cloud_fraction_two_moment_init( & @@ -162,16 +168,13 @@ pure subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght - ! Statement functions - logical land - land = nint(landfrac) == 1 ! ---------------- ! ! Main computation ! ! ---------------- ! if (p >= premib) then - if (land .and. (snowh <= 0.000001_kind_phys)) then + if (is_land(landfrac) .and. (snowh <= 0.000001_kind_phys)) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl @@ -295,9 +298,6 @@ pure subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out real(kind_phys) :: rhmin ! Critical RH real(kind_phys) :: rhwght - ! Statement functions - logical land - land(i) = nint(landfrac_in(i)) == 1 ! ---------------- ! ! Main computation ! @@ -318,7 +318,7 @@ pure subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out rhminh = rhminh_in(i) if (p >= premib) then - if (land(i) .and. (snowh <= 0.000001_kind_phys)) then + if (is_land(landfrac_in(i)) .and. (snowh <= 0.000001_kind_phys)) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl @@ -425,9 +425,6 @@ pure subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl real(kind_phys) rhdif ! Factor for stratus fraction real(kind_phys) rhwght - ! Statement functions - logical land - land = nint(landfrac) == 1 ! ---------------- ! ! Main computation ! @@ -435,7 +432,7 @@ pure subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl if (p >= premib) then - if (land .and. (snowh <= 0.000001_kind_phys)) then + if (is_land(landfrac) .and. (snowh <= 0.000001_kind_phys)) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl @@ -524,9 +521,6 @@ pure subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out real(kind_phys) :: rhdif ! Factor for stratus fraction real(kind_phys) :: rhwght - ! Statement functions - logical land - land(i) = nint(landfrac_in(i)) == 1 ! ---------------- ! ! Main computation ! @@ -547,7 +541,7 @@ pure subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out rhminh = rhminh_in(i) if (p >= premib) then - if (land(i) .and. (snowh <= 0.000001_kind_phys)) then + if (is_land(landfrac_in(i)) .and. (snowh <= 0.000001_kind_phys)) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl @@ -633,9 +627,6 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & real(kind_phys) :: icimr ! in cloud ice mixing ratio real(kind_phys) :: rhdif ! working variable for slingo scheme - ! Statement functions - logical land - land = nint(landfrac) == 1 if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys @@ -664,7 +655,7 @@ subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, & aist = max(0._kind_phys, min(aist, 1._kind_phys)) elseif (iceopt == 4) then if (p >= premib) then - if (land .and. (snowh <= 0.000001_kind_phys)) then + if (is_land(landfrac) .and. (snowh <= 0.000001_kind_phys)) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl @@ -822,9 +813,6 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a integer :: i - ! Statement functions - logical land - land(i) = nint(landfrac_in(i)) == 1 if (present(qsatfac_out)) qsatfac_out = 1.0_kind_phys @@ -874,7 +862,7 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a aist = max(0._kind_phys, min(aist, 1._kind_phys)) elseif (iceopt == 4) then if (p >= premib) then - if (land(i) .and. (snowh <= 0.000001_kind_phys)) then + if (is_land(landfrac_in(i)) .and. (snowh <= 0.000001_kind_phys)) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl From b64314a38edaa324952287ee72808085adf7f6cd Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 15 Jun 2026 12:55:04 -0400 Subject: [PATCH 8/9] Address review comments --- .../compute_cloud_fraction_two_moment.F90 | 22 +++++++++++-------- suites/suite_cam5.xml | 9 ++++++++ suites/suite_cam7.xml | 11 +++++++++- 3 files changed, 32 insertions(+), 10 deletions(-) diff --git a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 index 18966980..3f4fe54c 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 +++ b/schemes/cloud_fraction/compute_cloud_fraction_two_moment.F90 @@ -41,6 +41,8 @@ module compute_cloud_fraction_two_moment ! 3 = Wood and Field ! 4 = Wilson (based on Smith) ! 5 = modified Slingo (ssat & empty cloud) + ! 6 = fit based on T and Number + ! using observations from Heymsfield et al. (2013), Boudala (2002) integer :: iceopt ! Ice cloud fraction fit parameters (used by aist_single and aist_vector) @@ -72,11 +74,6 @@ module compute_cloud_fraction_two_moment contains - pure elemental logical function is_land(landfrac) - real(kind_phys), intent(in) :: landfrac - is_land = nint(landfrac) == 1 - end function is_land - !> \section arg_table_compute_cloud_fraction_two_moment_init Argument Table !! \htmlinclude compute_cloud_fraction_two_moment_init.html subroutine compute_cloud_fraction_two_moment_init( & @@ -131,6 +128,11 @@ subroutine compute_cloud_fraction_two_moment_init( & end subroutine compute_cloud_fraction_two_moment_init + pure elemental logical function is_land(landfrac) + real(kind_phys), intent(in) :: landfrac + is_land = nint(landfrac) == 1 + end function is_land + pure subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl_adj_land, rhminh, orhmin) ! --------------------------------------------------------- ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the @@ -421,9 +423,9 @@ pure subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, rhminl, rhminl real(kind_phys), optional, intent(out) :: orhmin ! Critical RH ! Local variables - real(kind_phys) rhmin ! Critical RH - real(kind_phys) rhdif ! Factor for stratus fraction - real(kind_phys) rhwght + real(kind_phys) :: rhmin ! Critical RH + real(kind_phys) :: rhdif ! Factor for stratus fraction + real(kind_phys) :: rhwght ! ---------------- ! @@ -806,7 +808,9 @@ subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, a real(kind_phys) :: nil ! Ice number concentration [#/L] ! Heymsfield/Boudala fit parameters (iceopt=6) - ! Boudala et al. (2002), Heymsfield et al. (2012) + ! Boudala et al. (2002), Heymsfield et al. (2013) + ! https://doi.org/10.1002/joc.774 + ! https://doi.org/10.1175/JAS-D-12-0124.1 real(kind_phys), parameter :: heymsfield_a = 6.73834e-08_kind_phys real(kind_phys), parameter :: heymsfield_b = 0.0533110_kind_phys real(kind_phys), parameter :: heymsfield_c = 0.3493813_kind_phys diff --git a/suites/suite_cam5.xml b/suites/suite_cam5.xml index 1de13ff5..51444771 100644 --- a/suites/suite_cam5.xml +++ b/suites/suite_cam5.xml @@ -23,6 +23,9 @@ zm_conv_options + + compute_cloud_fraction_two_moment + check_energy_gmean @@ -85,6 +88,12 @@ Not implemented yet. --> + + compute_cloud_fraction + sima_state_diagnostics diff --git a/suites/suite_cam7.xml b/suites/suite_cam7.xml index 76fb3e13..84d45b02 100644 --- a/suites/suite_cam7.xml +++ b/suites/suite_cam7.xml @@ -99,7 +99,16 @@ tropopause_find tropopause_diagnostics - + + compute_cloud_fraction + compute_cloud_fraction_two_moment + + + +