module cldfrc2m ! cloud fraction calculations use shr_kind_mod, only: r8=>shr_kind_r8 use spmd_utils, only: masterproc use ppgrid, only: pcols use physconst, only: rair use wv_saturation, only: qsat_water, svp_water, svp_ice use cam_logfile, only: iulog use cam_abortutils, only: endrun implicit none private save public :: & cldfrc2m_readnl, & cldfrc2m_init, & astG_PDF_single, & astG_PDF, & astG_RHU_single, & astG_RHU, & aist_single, & aist_vector, & CAMstfrac, & rhmini_const, & rhmaxi_const ! Namelist variables real(r8) :: cldfrc2m_rhmini ! Minimum rh for ice cloud fraction > 0. real(r8) :: cldfrc2m_rhmaxi logical :: cldfrc2m_do_subgrid_growth = .false. ! -------------------------- ! ! Parameters for Ice Stratus ! ! -------------------------- ! real(r8), protected :: rhmini_const ! Minimum rh for ice cloud fraction > 0. real(r8), protected :: rhmaxi_const real(r8), parameter :: qist_min = 1.e-7_r8 ! Minimum in-stratus ice IWC constraint [ kg/kg ] real(r8), parameter :: qist_max = 5.e-3_r8 ! Maximum in-stratus ice IWC constraint [ kg/kg ] ! ----------------------------- ! ! Parameters for Liquid Stratus ! ! ----------------------------- ! 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 real(r8) :: rhminl_const ! Critical RH for low-level liquid stratus clouds real(r8) :: rhminl_adj_land_const ! rhminl adjustment for snowfree land real(r8) :: rhminh_const ! Critical RH for high-level liquid stratus clouds real(r8) :: premit ! Top height for mid-level liquid stratus fraction real(r8) :: premib ! Bottom height for mid-level liquid stratus fraction integer :: iceopt ! option for ice cloud closure ! 1=wang & sassen 2=schiller (iciwc) ! 3=wood & field, 4=Wilson (based on smith) ! 5=modified slingo (ssat & empyt cloud) real(r8) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure ! ( smaller = more ice clouds ) !================================================================================================ contains !================================================================================================ subroutine cldfrc2m_readnl(nlfile) use namelist_utils, only: find_group_name use units, only: getunit, freeunit use spmd_utils, only: mpicom, masterprocid, mpi_logical, mpi_real8 character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input ! Local variables integer :: unitn, ierr character(len=*), parameter :: subname = 'cldfrc2m_readnl' namelist /cldfrc2m_nl/ cldfrc2m_rhmini, cldfrc2m_rhmaxi, cldfrc2m_do_subgrid_growth !----------------------------------------------------------------------------- if (masterproc) then unitn = getunit() open( unitn, file=trim(nlfile), status='old' ) call find_group_name(unitn, 'cldfrc2m_nl', status=ierr) if (ierr == 0) then read(unitn, cldfrc2m_nl, iostat=ierr) if (ierr /= 0) then call endrun(subname // ':: ERROR reading namelist') end if end if close(unitn) call freeunit(unitn) ! set local variables rhmini_const = cldfrc2m_rhmini rhmaxi_const = cldfrc2m_rhmaxi end if ! Broadcast namelist variables call mpi_bcast(rhmini_const, 1, mpi_real8, masterprocid, mpicom, ierr) call mpi_bcast(rhmaxi_const, 1, mpi_real8, masterprocid, mpicom, ierr) call mpi_bcast(cldfrc2m_do_subgrid_growth, 1, mpi_logical,masterprocid, mpicom, ierr) end subroutine cldfrc2m_readnl !================================================================================================ subroutine cldfrc2m_init() use cloud_fraction, only: cldfrc_getparams call cldfrc_getparams(rhminl_out=rhminl_const, rhminl_adj_land_out=rhminl_adj_land_const, & rhminh_out=rhminh_const, premit_out=premit, premib_out=premib, & iceopt_out=iceopt, icecrit_out=icecrit) if( masterproc ) then write(iulog,*) 'cldfrc2m parameters:' write(iulog,*) ' rhminl = ', rhminl_const write(iulog,*) ' rhminl_adj_land = ', rhminl_adj_land_const write(iulog,*) ' rhminh = ', rhminh_const write(iulog,*) ' premit = ', premit write(iulog,*) ' premib = ', premib write(iulog,*) ' iceopt = ', iceopt write(iulog,*) ' icecrit = ', icecrit write(iulog,*) ' rhmini = ', rhmini_const write(iulog,*) ' rhmaxi = ', rhmaxi_const write(iulog,*) ' do_subgrid_growth = ', cldfrc2m_do_subgrid_growth end if end subroutine cldfrc2m_init !================================================================================================ subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, orhmin, & 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(r8), intent(in) :: U ! Relative humidity real(r8), intent(in) :: p ! Pressure [Pa] real(r8), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg] real(r8), intent(in) :: landfrac ! Land fraction real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent) real(r8), intent(out) :: a ! Stratus fraction real(r8), intent(out) :: Ga ! dU/da real(r8), optional, intent(out) :: orhmin ! Critical RH real(r8), optional, intent(in) :: rhminl_in ! Critical relative humidity for low-level liquid stratus real(r8), optional, intent(in) :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land real(r8), optional, intent(in) :: rhminh_in ! Critical relative humidity for high-level liquid stratus ! Local variables integer :: i ! Loop indexes real(r8) dV ! Width of triangular PDF real(r8) cldrh ! RH of stratus cloud real(r8) rhmin ! Critical RH real(r8) rhwght real(r8) :: rhminl real(r8) :: rhminl_adj_land real(r8) :: rhminh ! Statement functions logical land land = nint(landfrac) == 1 ! ---------- ! ! Parameters ! ! ---------- ! cldrh = 1.0_r8 rhminl = rhminl_const if (present(rhminl_in)) rhminl = rhminl_in rhminl_adj_land = rhminl_adj_land_const if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in rhminh = rhminh_const if (present(rhminh_in)) rhminh = rhminh_in ! ---------------- ! ! Main computation ! ! ---------------- ! if( p .ge. premib ) then if( land .and. (snowh.le.0.000001_r8) ) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl endif dV = cldrh - rhmin if( U .ge. 1._r8 ) then a = 1._r8 Ga = 1.e10_r8 elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8) Ga = dV/sqrt(2._r8)*sqrt(1._r8-a) elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8 Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a)) elseif( U .le. (cldrh-dV) ) then a = 0._r8 Ga = 1.e10_r8 endif if( freeze_dry ) then a = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) endif elseif( p .lt. premit ) then rhmin = rhminh dV = cldrh - rhmin if( U .ge. 1._r8 ) then a = 1._r8 Ga = 1.e10_r8 elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8) Ga = dV/sqrt(2._r8)*sqrt(1._r8-a) elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8 Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a)) elseif( U .le. (cldrh-dV) ) then a = 0._r8 Ga = 1.e10_r8 endif else rhwght = (premib-(max(p,premit)))/(premib-premit) ! if( land .and. (snowh.le.0.000001_r8) ) then ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght) ! else rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght) ! endif dV = cldrh - rhmin if( U .ge. 1._r8 ) then a = 1._r8 Ga = 1.e10_r8 elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8) Ga = dV/sqrt(2._r8)*sqrt(1._r8-a) elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8 Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a)) elseif( U .le. (cldrh-dV) ) then a = 0._r8 Ga = 1.e10_r8 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(r8), intent(in) :: U_in(pcols) ! Relative humidity real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa] real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor specific humidity [kg/kg] real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent) real(r8), intent(out) :: a_out(pcols) ! Stratus fraction real(r8), intent(out) :: Ga_out(pcols) ! dU/da integer, intent(in) :: ncol real(r8), optional, intent(in) :: rhminl_in(pcols) ! Critical relative humidity for low-level liquid stratus real(r8), optional, intent(in) :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land real(r8), optional, intent(in) :: rhminh_in(pcols) ! Critical relative humidity for high-level liquid stratus real(r8) :: rhminl ! Critical relative humidity for low-level liquid stratus real(r8) :: rhminl_adj_land ! Adjustment drop of rhminl over the land real(r8) :: rhminh ! Critical relative humidity for high-level liquid stratus real(r8) :: U ! Relative humidity real(r8) :: p ! Pressure [Pa] real(r8) :: qv ! Grid-mean water vapor specific humidity [kg/kg] real(r8) :: landfrac ! Land fraction real(r8) :: snowh ! Snow depth (liquid water equivalent) real(r8) :: a ! Stratus fraction real(r8) :: Ga ! dU/da ! Local variables integer :: i ! Loop indexes real(r8) dV ! Width of triangular PDF real(r8) cldrh ! RH of stratus cloud real(r8) rhmin ! Critical RH real(r8) rhwght ! Statement functions logical land land(i) = nint(landfrac_in(i)) == 1 ! ---------- ! ! Parameters ! ! ---------- ! cldrh = 1.0_r8 rhminl = rhminl_const rhminl_adj_land = rhminl_adj_land_const rhminh = rhminh_const ! ---------------- ! ! Main computation ! ! ---------------- ! a_out(:) = 0._r8 Ga_out(:) = 0._r8 do i = 1, ncol U = U_in(i) p = p_in(i) qv = qv_in(i) landfrac = landfrac_in(i) snowh = snowh_in(i) if (present(rhminl_in)) rhminl = rhminl_in(i) if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i) if (present(rhminh_in)) rhminh = rhminh_in(i) if( p .ge. premib ) then if( land(i) .and. (snowh.le.0.000001_r8) ) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl endif dV = cldrh - rhmin if( U .ge. 1._r8 ) then a = 1._r8 Ga = 1.e10_r8 elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8) Ga = dV/sqrt(2._r8)*sqrt(1._r8-a) elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8 Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a)) elseif( U .le. (cldrh-dV) ) then a = 0._r8 Ga = 1.e10_r8 endif if( freeze_dry ) then a = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) endif elseif( p .lt. premit ) then rhmin = rhminh dV = cldrh - rhmin if( U .ge. 1._r8 ) then a = 1._r8 Ga = 1.e10_r8 elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8) Ga = dV/sqrt(2._r8)*sqrt(1._r8-a) elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8 Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a)) elseif( U .le. (cldrh-dV) ) then a = 0._r8 Ga = 1.e10_r8 endif else rhwght = (premib-(max(p,premit)))/(premib-premit) ! if( land(i) .and. (snowh.le.0.000001_r8) ) then ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght) ! else rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght) ! endif dV = cldrh - rhmin if( U .ge. 1._r8 ) then a = 1._r8 Ga = 1.e10_r8 elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8) Ga = dV/sqrt(2._r8)*sqrt(1._r8-a) elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8 Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a)) elseif( U .le. (cldrh-dV) ) then a = 0._r8 Ga = 1.e10_r8 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, orhmin, & 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(r8), intent(in) :: U ! Relative humidity real(r8), intent(in) :: p ! Pressure [Pa] real(r8), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg] real(r8), intent(in) :: landfrac ! Land fraction real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent) real(r8), intent(out) :: a ! Stratus fraction real(r8), intent(out) :: Ga ! dU/da real(r8), optional, intent(out) :: orhmin ! Critical RH real(r8), optional, intent(in) :: rhminl_in ! Critical relative humidity for low-level liquid stratus real(r8), optional, intent(in) :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land real(r8), optional, intent(in) :: rhminh_in ! Critical relative humidity for high-level liquid stratus ! Local variables real(r8) rhmin ! Critical RH real(r8) rhdif ! Factor for stratus fraction real(r8) rhwght real(r8) :: rhminl real(r8) :: rhminl_adj_land real(r8) :: rhminh ! Statement functions logical land land = nint(landfrac) == 1 rhminl = rhminl_const if (present(rhminl_in)) rhminl = rhminl_in rhminl_adj_land = rhminl_adj_land_const if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in rhminh = rhminh_const if (present(rhminh_in)) rhminh = rhminh_in ! ---------------- ! ! Main computation ! ! ---------------- ! if( p .ge. premib ) then if( land .and. (snowh.le.0.000001_r8) ) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl endif rhdif = (U-rhmin)/(1.0_r8-rhmin) a = min(1._r8,(max(rhdif,0.0_r8))**2) if( (U.ge.1._r8) .or. (U.le.rhmin) ) then Ga = 1.e20_r8 else Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin)) endif if( freeze_dry ) then a = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) endif elseif( p .lt. premit ) then rhmin = rhminh rhdif = (U-rhmin)/(1.0_r8-rhmin) a = min(1._r8,(max(rhdif,0._r8))**2) if( (U.ge.1._r8) .or. (U.le.rhmin) ) then Ga = 1.e20_r8 else Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin)) endif else rhwght = (premib-(max(p,premit)))/(premib-premit) ! if( land .and. (snowh.le.0.000001_r8) ) then ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght) ! else rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght) ! endif rhdif = (U-rhmin)/(1.0_r8-rhmin) a = min(1._r8,(max(rhdif,0._r8))**2) if( (U.ge.1._r8) .or. (U.le.rhmin) ) then Ga = 1.e10_r8 else Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-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(r8), intent(in) :: U_in(pcols) ! Relative humidity real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa] real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor specific humidity [kg/kg] real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent) real(r8), intent(out) :: a_out(pcols) ! Stratus fraction real(r8), intent(out) :: Ga_out(pcols) ! dU/da integer, intent(in) :: ncol real(r8), optional, intent(in) :: rhminl_in(pcols) ! Critical relative humidity for low-level liquid stratus real(r8), optional, intent(in) :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land real(r8), optional, intent(in) :: rhminh_in(pcols) ! Critical relative humidity for high-level liquid stratus real(r8) :: U ! Relative humidity real(r8) :: p ! Pressure [Pa] real(r8) :: qv ! Grid-mean water vapor specific humidity [kg/kg] real(r8) :: landfrac ! Land fraction real(r8) :: snowh ! Snow depth (liquid water equivalent) real(r8) :: rhminl ! Critical relative humidity for low-level liquid stratus real(r8) :: rhminl_adj_land ! Adjustment drop of rhminl over the land real(r8) :: rhminh ! Critical relative humidity for high-level liquid stratus real(r8) :: a ! Stratus fraction real(r8) :: Ga ! dU/da ! Local variables integer i real(r8) rhmin ! Critical RH real(r8) rhdif ! Factor for stratus fraction real(r8) rhwght ! Statement functions logical land land(i) = nint(landfrac_in(i)) == 1 rhminl = rhminl_const rhminl_adj_land = rhminl_adj_land_const rhminh = rhminh_const ! ---------------- ! ! Main computation ! ! ---------------- ! a_out(:) = 0._r8 Ga_out(:) = 0._r8 do i = 1, ncol U = U_in(i) p = p_in(i) qv = qv_in(i) landfrac = landfrac_in(i) snowh = snowh_in(i) if (present(rhminl_in)) rhminl = rhminl_in(i) if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i) if (present(rhminh_in)) rhminh = rhminh_in(i) if( p .ge. premib ) then if( land(i) .and. (snowh.le.0.000001_r8) ) then rhmin = rhminl - rhminl_adj_land else rhmin = rhminl endif rhdif = (U-rhmin)/(1.0_r8-rhmin) a = min(1._r8,(max(rhdif,0.0_r8))**2) if( (U.ge.1._r8) .or. (U.le.rhmin) ) then Ga = 1.e20_r8 else Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin)) endif if( freeze_dry ) then a = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) endif elseif( p .lt. premit ) then rhmin = rhminh rhdif = (U-rhmin)/(1.0_r8-rhmin) a = min(1._r8,(max(rhdif,0._r8))**2) if( (U.ge.1._r8) .or. (U.le.rhmin) ) then Ga = 1.e20_r8 else Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin)) endif else rhwght = (premib-(max(p,premit)))/(premib-premit) ! if( land(i) .and. (snowh.le.0.000001_r8) ) then ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght) ! else rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght) ! endif rhdif = (U-rhmin)/(1.0_r8-rhmin) a = min(1._r8,(max(rhdif,0._r8))**2) if( (U.ge.1._r8) .or. (U.le.rhmin) ) then Ga = 1.e10_r8 else Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-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_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, & qsatfac_out) ! --------------------------------------------------------- ! ! Compute non-physical ice stratus fraction ! ! --------------------------------------------------------- ! real(r8), intent(in) :: qv ! Grid-mean water vapor[kg/kg] real(r8), intent(in) :: T ! Temperature real(r8), intent(in) :: p ! Pressure [Pa] real(r8), intent(in) :: qi ! Grid-mean ice water content [kg/kg] real(r8), intent(in) :: landfrac ! Land fraction real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent) real(r8), intent(out) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) real(r8), optional, intent(in) :: rhmaxi_in real(r8), optional, intent(in) :: rhmini_in ! Critical relative humidity for ice stratus real(r8), optional, intent(in) :: rhminl_in ! Critical relative humidity for low-level liquid stratus real(r8), optional, intent(in) :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land real(r8), optional, intent(in) :: rhminh_in ! Critical relative humidity for high-level liquid stratus real(r8), optional, intent(out) :: qsatfac_out ! Subgrid scaling factor for qsat ! Local variables real(r8) rhmin ! Critical RH real(r8) rhwght real(r8) a,b,c,as,bs,cs ! Fit parameters real(r8) Kc ! Constant for ice cloud calc (wood & field) real(r8) ttmp ! Limited temperature real(r8) icicval ! Empirical IWC value [ kg/kg ] real(r8) rho ! Local air density real(r8) esl ! Liq sat vapor pressure real(r8) esi ! Ice sat vapor pressure real(r8) ncf,phi ! Wilson and Ballard parameters real(r8) es, qs real(r8) rhi ! grid box averaged relative humidity over ice real(r8) minice ! minimum grid box avg ice for having a 'cloud' real(r8) mincld ! minimum ice cloud fraction threshold real(r8) icimr ! in cloud ice mixing ratio real(r8) rhdif ! working variable for slingo scheme real(r8) :: rhmaxi real(r8) :: rhmini real(r8) :: rhminl real(r8) :: rhminl_adj_land real(r8) :: rhminh ! Statement functions logical land land = nint(landfrac) == 1 ! --------- ! ! Constants ! ! --------- ! ! Wang and Sassen IWC paramters ( Option.1 ) a = 26.87_r8 b = 0.569_r8 c = 0.002892_r8 ! Schiller parameters ( Option.2 ) as = -68.4202_r8 bs = 0.983917_r8 cs = 2.81795_r8 ! Wood and Field parameters ( Option.3 ) Kc = 75._r8 ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) ! Slingo modified (option 5) minice = 1.e-12_r8 mincld = 1.e-4_r8 rhmaxi = rhmaxi_const if (present(rhmaxi_in)) rhmaxi = rhmaxi_in rhmini = rhmini_const if (present(rhmini_in)) rhmini = rhmini_in rhminl = rhminl_const if (present(rhminl_in)) rhminl = rhminl_in rhminl_adj_land = rhminl_adj_land_const if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in rhminh = rhminh_const if (present(rhminh_in)) rhminh = rhminh_in if (present(qsatfac_out)) qsatfac_out = 1.0_r8 ! ---------------- ! ! 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._r8,min(T,253._r8)) - 273.16_r8 icicval = a + b * ttmp + c * ttmp**2._r8 rho = p/(rair*T) icicval = icicval * 1.e-6_r8 / rho else ttmp = max(190._r8,min(T,273.16_r8)) icicval = 10._r8 **(as * bs**ttmp + cs) icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8 endif aist = max(0._r8,min(qi/icicval,1._r8)) elseif( iceopt.eq.3 ) then aist = 1._r8 - exp(-Kc*qi/(qs*(esi/esl))) aist = max(0._r8,min(aist,1._r8)) elseif( iceopt.eq.4) then if( p .ge. premib ) then if( land .and. (snowh.le.0.000001_r8) ) 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_r8) ) then ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght) ! else rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght) ! endif endif ncf = qi/((1._r8 - icecrit)*qs) if( ncf.le.0._r8 ) then aist = 0._r8 elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8) elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8 aist = (1._r8 - 4._r8 * cos(phi) * cos(phi)) else aist = 1._r8 endif aist = max(0._r8,min(aist,1._r8)) elseif (iceopt.eq.5) then ! set rh ice cloud fraction rhi= (qv+qi)/qs * (esl/esi) rhdif= (rhi-rhmini) / (rhmaxi - rhmini) aist = min(1.0_r8, max(rhdif,0._r8)**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. cldfrc2m_do_subgrid_growth) then qsatfac_out = max(min(qv / qs, 1._r8), (1._r8 - 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._r8 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 aist = max(0._r8,min(1._r8,qi/qist_min)) endif !maximum if (icimr.gt.qist_max) then aist = max(0._r8,min(1._r8,qi/qist_max)) endif endif endif ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate ! computed after updating 'qi_st'. aist = max(0._r8,min(aist,0.999_r8)) 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 ! ! --------------------------------------------------------- ! real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor[kg/kg] real(r8), intent(in) :: T_in(pcols) ! Temperature real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa] real(r8), intent(in) :: qi_in(pcols) ! Grid-mean ice water content [kg/kg] real(r8), intent(in) :: ni_in(pcols) ! Grid-mean ice water number concentration [#/kg] real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent) real(r8), intent(out) :: aist_out(pcols) ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) integer, intent(in) :: ncol real(r8), optional, intent(in) :: rhmaxi_in real(r8), optional, intent(in) :: rhmini_in(pcols) ! Critical relative humidity for ice stratus real(r8), optional, intent(in) :: rhminl_in(pcols) ! Critical relative humidity for low-level liquid stratus real(r8), optional, intent(in) :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land real(r8), optional, intent(in) :: rhminh_in(pcols) ! Critical relative humidity for high-level liquid stratus real(r8), optional, intent(out) :: qsatfac_out(pcols) ! Subgrid scaling factor for qsat ! Local variables real(r8) qv ! Grid-mean water vapor[kg/kg] real(r8) T ! Temperature real(r8) p ! Pressure [Pa] real(r8) qi ! Grid-mean ice water content [kg/kg] real(r8) ni real(r8) landfrac ! Land fraction real(r8) snowh ! Snow depth (liquid water equivalent) real(r8) rhmaxi ! Critical relative humidity for ice stratus real(r8) rhmini ! Critical relative humidity for ice stratus real(r8) rhminl ! Critical relative humidity for low-level liquid stratus real(r8) rhminl_adj_land ! Adjustment drop of rhminl over the land real(r8) rhminh ! Critical relative humidity for high-level liquid stratus real(r8) aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 ) real(r8) rhmin ! Critical RH real(r8) rhwght real(r8) a,b,c,as,bs,cs,ah,bh,ch ! Fit parameters real(r8) nil real(r8) Kc ! Constant for ice cloud calc (wood & field) real(r8) ttmp ! Limited temperature real(r8) icicval ! Empirical IWC value [ kg/kg ] real(r8) rho ! Local air density real(r8) esl ! Liq sat vapor pressure real(r8) esi ! Ice sat vapor pressure real(r8) ncf,phi ! Wilson and Ballard parameters real(r8) qs real(r8) esat_in(pcols) real(r8) qsat_in(pcols) real(r8) rhi ! grid box averaged relative humidity over ice real(r8) minice ! minimum grid box avg ice for having a 'cloud' real(r8) mincld ! minimum ice cloud fraction threshold real(r8) icimr ! in cloud ice mixing ratio real(r8) 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_r8 b = 0.569_r8 c = 0.002892_r8 ! Schiller parameters ( Option.2 ) as = -68.4202_r8 bs = 0.983917_r8 cs = 2.81795_r8 ! Wood and Field parameters ( Option.3 ) Kc = 75._r8 ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds) ! Slingo modified (option 5) minice = 1.e-12_r8 mincld = 1.e-4_r8 rhmaxi = rhmaxi_const if (present(rhmaxi_in)) rhmaxi = rhmaxi_in rhmini = rhmini_const rhminl = rhminl_const rhminl_adj_land = rhminl_adj_land_const rhminh = rhminh_const if (present(qsatfac_out)) qsatfac_out = 1.0_r8 ! ---------------- ! ! Main computation ! ! ---------------- ! aist_out(:) = 0._r8 esat_in(:) = 0._r8 qsat_in(:) = 0._r8 call qsat_water(T_in(1:ncol), p_in(1:ncol), & esat_in(1:ncol), qsat_in(1: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) esl = svp_water(T) esi = svp_ice(T) if (present(rhmini_in)) rhmini = rhmini_in(i) if (present(rhminl_in)) rhminl = rhminl_in(i) if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i) if (present(rhminh_in)) rhminh = rhminh_in(i) if( iceopt.lt.3 ) then if( iceopt.eq.1 ) then ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8 icicval = a + b * ttmp + c * ttmp**2._r8 rho = p/(rair*T) icicval = icicval * 1.e-6_r8 / rho else ttmp = max(190._r8,min(T,273.16_r8)) icicval = 10._r8 **(as * bs**ttmp + cs) icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8 endif aist = max(0._r8,min(qi/icicval,1._r8)) elseif( iceopt.eq.3 ) then aist = 1._r8 - exp(-Kc*qi/(qs*(esi/esl))) aist = max(0._r8,min(aist,1._r8)) elseif( iceopt.eq.4) then if( p .ge. premib ) then if( land(i) .and. (snowh.le.0.000001_r8) ) 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_r8) ) then ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght) ! else rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght) ! endif endif ncf = qi/((1._r8 - icecrit)*qs) if( ncf.le.0._r8 ) then aist = 0._r8 elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8) elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8 aist = (1._r8 - 4._r8 * cos(phi) * cos(phi)) else aist = 1._r8 endif aist = max(0._r8,min(aist,1._r8)) elseif (iceopt.eq.5) then ! set rh ice cloud fraction rhi= (qv+qi)/qs * (esl/esi) rhdif= (rhi-rhmini) / (rhmaxi - rhmini) aist = min(1.0_r8, max(rhdif,0._r8)**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_r8 bh= 0.0533110_r8 ch= 0.3493813_r8 rho=p/(rair*T) nil=ni*rho/1000._r8 icicval = ah * exp(bh*T) * nil**ch !result is in g m-3, convert to kg H2O / kg air (icimr...) icicval = icicval / rho / 1000._r8 aist = max(0._r8,min(qi/icicval,1._r8)) aist = min(aist,1._r8) 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. cldfrc2m_do_subgrid_growth) then qsatfac_out(i) = max(min(qv / qs, 1._r8), (1._r8 - 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._r8 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 aist = max(0._r8,min(1._r8,qi/qist_min)) endif !maximum if (icimr.gt.qist_max) then aist = max(0._r8,min(1._r8,qi/qist_max)) endif endif endif ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate ! computed after updating 'qi_st'. aist = max(0._r8,min(aist,0.999_r8)) aist_out(i) = aist enddo end subroutine aist_vector !================================================================================================ end module cldfrc2m