diff --git a/.gitmodules b/.gitmodules index 6fe486cfba..d77aef9a59 100644 --- a/.gitmodules +++ b/.gitmodules @@ -29,7 +29,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_22_003 + fxtag = atmos_phys0_23_000 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 6a2be626e7..ae05948275 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -3268,7 +3268,7 @@ Default: 700.e2 for CAM5, CAM6 and CAM7; all others=> 750.e2 Scheme for ice cloud fraction: 1=wang & sassen, 2=schiller (iciwc), -3=wood & field, 4=Wilson (based on smith), 5=modified slingo (ssat & empyt cloud) +3=wood & field, 4=Wilson (based on smith), 5=modified slingo (ssat & empty cloud) Default: set by build-namelist @@ -3304,7 +3304,7 @@ Default: set by build-namelist -Minimum in-stratus IWC constraint [ kg/kg ] +Minimum in-stratus IWC constraint greater than 0 [ kg/kg ] Default: set by build-namelist diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index 3003033f7a..d5227f0612 100644 --- a/cime_config/testdefs/testlist_cam.xml +++ b/cime_config/testdefs/testlist_cam.xml @@ -334,14 +334,14 @@ - + - + diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/shell_commands new file mode 100644 index 0000000000..eb40ad83e0 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/shell_commands @@ -0,0 +1,2 @@ +./xmlchange ROF_NCPL=\$ATM_NCPL +./xmlchange GLC_NCPL=\$ATM_NCPL diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_cam new file mode 100644 index 0000000000..fc9a71ab44 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_cam @@ -0,0 +1,78 @@ +mfilt=1,1,1,1,1,1 +ndens=1,1,1,1,1,1 +nhtfrq=9,9,9,9,9,9 +write_nstep0=.true. +inithist='ENDOFRUN' + +rad_diag_1 = 'A:Q:H2O', 'N:O2:O2', 'A:CO2:CO2', 'A:O3:O3', 'A:N2O:N2O', 'A:CH4:CH4', + 'N:CFC11:CFC11', 'N:CFC12:CFC12' + +rad_diag_2 = 'A:Q:H2O', 'N:O2:O2', 'A:CO2:CO2', 'A:O3:O3', 'A:N2O:N2O', 'A:CH4:CH4', + 'N:CFC11:CFC11', 'N:CFC12:CFC12', + 'B:MXAER01:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX01_rrtmg.nc', + 'B:MXAER02:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX02_rrtmg.nc', + 'B:MXAER03:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX03_rrtmg.nc', + 'B:MXAER04:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX04_rrtmg.nc', + 'B:MXAER05:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX05_rrtmg.nc', + 'B:MXAER06:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX06_rrtmg.nc', + 'B:MXAER07:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX07_rrtmg.nc', + 'B:MXAER08:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX08_rrtmg.nc', + 'B:MXAER09:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX09_rrtmg.nc', + 'B:MXAER10:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX10_rrtmg.nc', + 'B:MXAER11:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX11_rrtmg.nc', + 'B:MXAER12:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX12_rrtmg.nc', + 'B:MXAER13:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX13_rrtmg.nc', + 'B:MXAER14:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX14_rrtmg.nc', + 'B:MXAER15:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX15_rrtmg.nc', + 'B:MXAER16:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX16_rrtmg.nc', + 'B:MXAER17:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX17_rrtmg.nc', + 'B:MXAER18:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX18_rrtmg.nc', + 'B:MXAER19:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX19_rrtmg.nc', + 'B:MXAER20:$DIN_LOC_ROOT/atm/cam/chem/carma/rrtmg/aerosol_cld_CRMIX20_rrtmg.nc' + +! rad_diag_3 uses the aerosol optics files generated by CARMA in the run +! directory at initialization when carma_do_optics is enabled. +carma_do_optics = .true. + +rad_diag_3 = 'A:Q:H2O', 'N:O2:O2', 'A:CO2:CO2', 'A:O3:O3', 'A:N2O:N2O', 'A:CH4:CH4', + 'N:CFC11:CFC11', 'N:CFC12:CFC12', + 'B:MXAER01:trop_strat_soa5_MXAER01_rrtmg.nc', + 'B:MXAER02:trop_strat_soa5_MXAER02_rrtmg.nc', + 'B:MXAER03:trop_strat_soa5_MXAER03_rrtmg.nc', + 'B:MXAER04:trop_strat_soa5_MXAER04_rrtmg.nc', + 'B:MXAER05:trop_strat_soa5_MXAER05_rrtmg.nc', + 'B:MXAER06:trop_strat_soa5_MXAER06_rrtmg.nc', + 'B:MXAER07:trop_strat_soa5_MXAER07_rrtmg.nc', + 'B:MXAER08:trop_strat_soa5_MXAER08_rrtmg.nc', + 'B:MXAER09:trop_strat_soa5_MXAER09_rrtmg.nc', + 'B:MXAER10:trop_strat_soa5_MXAER10_rrtmg.nc', + 'B:MXAER11:trop_strat_soa5_MXAER11_rrtmg.nc', + 'B:MXAER12:trop_strat_soa5_MXAER12_rrtmg.nc', + 'B:MXAER13:trop_strat_soa5_MXAER13_rrtmg.nc', + 'B:MXAER14:trop_strat_soa5_MXAER14_rrtmg.nc', + 'B:MXAER15:trop_strat_soa5_MXAER15_rrtmg.nc', + 'B:MXAER16:trop_strat_soa5_MXAER16_rrtmg.nc', + 'B:MXAER17:trop_strat_soa5_MXAER17_rrtmg.nc', + 'B:MXAER18:trop_strat_soa5_MXAER18_rrtmg.nc', + 'B:MXAER19:trop_strat_soa5_MXAER19_rrtmg.nc', + 'B:MXAER20:trop_strat_soa5_MXAER20_rrtmg.nc', + 'B:PRSUL01:trop_strat_soa5_PRSUL01_rrtmg.nc', + 'B:PRSUL02:trop_strat_soa5_PRSUL02_rrtmg.nc', + 'B:PRSUL03:trop_strat_soa5_PRSUL03_rrtmg.nc', + 'B:PRSUL04:trop_strat_soa5_PRSUL04_rrtmg.nc', + 'B:PRSUL05:trop_strat_soa5_PRSUL05_rrtmg.nc', + 'B:PRSUL06:trop_strat_soa5_PRSUL06_rrtmg.nc', + 'B:PRSUL07:trop_strat_soa5_PRSUL07_rrtmg.nc', + 'B:PRSUL08:trop_strat_soa5_PRSUL08_rrtmg.nc', + 'B:PRSUL09:trop_strat_soa5_PRSUL09_rrtmg.nc', + 'B:PRSUL10:trop_strat_soa5_PRSUL10_rrtmg.nc', + 'B:PRSUL11:trop_strat_soa5_PRSUL11_rrtmg.nc', + 'B:PRSUL12:trop_strat_soa5_PRSUL12_rrtmg.nc', + 'B:PRSUL13:trop_strat_soa5_PRSUL13_rrtmg.nc', + 'B:PRSUL14:trop_strat_soa5_PRSUL14_rrtmg.nc', + 'B:PRSUL15:trop_strat_soa5_PRSUL15_rrtmg.nc', + 'B:PRSUL16:trop_strat_soa5_PRSUL16_rrtmg.nc', + 'B:PRSUL17:trop_strat_soa5_PRSUL17_rrtmg.nc', + 'B:PRSUL18:trop_strat_soa5_PRSUL18_rrtmg.nc', + 'B:PRSUL19:trop_strat_soa5_PRSUL19_rrtmg.nc', + 'B:PRSUL20:trop_strat_soa5_PRSUL20_rrtmg.nc' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_clm new file mode 100644 index 0000000000..c4cb9d28d6 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_clm @@ -0,0 +1,3 @@ +hist_nhtfrq = 9 +hist_mfilt = 1 +hist_ndens = 1 diff --git a/doc/ChangeLog b/doc/ChangeLog index c3cc2ea0fa..7f502577f7 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,101 @@ =============================================================== +Tag name: cam6_4_183 +Originator(s): jimmielin +Date: 22 Jun 2026 +One-line Summary: Complete CCPPization of cloud fraction for two-moment microphys (cldfrc2m); Add rad_diag CARMA testmod to existing Derecho FCARMA2000climo test +Github PR URL: https://github.com/ESCOMP/CAM/pull/1522 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): +- Closes #1521 - Move two-moment cloud fraction scheme (cldfrc2m) to atmos_phys as compute_cloud_fraction_two_moment. +- cldfrc2m is used broadly across CAM physics (cldwat2m_macro, clubb_intr, CARMA cirrus models) and its namelist parameters are pulled in via USE statements in many places. Moving the computational routines to atmos_phys and making cldfrc2m a thin CAM shim for namelist I/O and parameter distribution facilitates CCPP port of CLUBB, PUMAS, etc. as these namelist parameters will cleanly flow via standard names +- Added keyword arguments calling into cldfrc for rh threshold etc. for clarity +- Closes #1560 - Extend CARMA test on Derecho/Intel with rad_diag testing + +Describe any changes made to build system: N/A + +Describe any changes made to the namelist: N/A + +List any changes to the defaults for the boundary datasets: N/A + +Describe any substantial timing or memory changes: N/A + +Code reviewed by: nusbaume, fvitt (CARMA fix), peverwhee (atmos_phys) + +List all files eliminated: N/A + +List all files added and what they do: +A cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/shell_commands +A cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_cam +A cime_config/testdefs/testmods_dirs/cam/outfrq9s_carma_rad_diag/user_nl_clm + - new testmod that exercises rad_diag for CARMA (prev only Izumi for MAM) + - new testmod that exercises carma_do_optics for CARMA to write out optics files + +List all existing files that have been modified, and describe the changes: + +M .gitmodules +M src/atmos_phys + - update atmos_phys to atmos_phys0_23_000. + +M bld/namelist_files/namelist_definition.xml + - fix typo + +M cime_config/testdefs/testlist_cam.xml + - update Derecho Intel FCARMA test to use the new outfrq9s_carma_rad_diag testmod. + +M src/physics/cam/cldfrc2m.F90 + - Move computational subroutines (astG_PDF, astG_RHU, astG_PDF_single, + astG_RHU_single, aist_single, aist_vector) to atmos_phys + compute_cloud_fraction_two_moment. + - Keep readnl, init, and public namelist-derived constants + (rhmini_const, rhmaxi_const, rhminis_const, rhmaxis_const, + rhminl_const, rhminl_adj_land_const, rhminh_const). + - init now calls compute_cloud_fraction_two_moment_init and + cldfrc_getparams to distribute parameters. + - Newly exposes rhminl_const, rhminl_adj_land_const, rhminh_const + (previously internal) so callers can pass them in explicitly. + +M src/physics/cam/cldwat2m_macro.F90 + - USE compute_cloud_fraction_two_moment directly for subroutines; + USE cldfrc2m shim only for rh constants. + +M src/physics/cam/clubb_intr.F90 + - USE compute_cloud_fraction_two_moment for aist_vector; + USE cldfrc2m shim for rh constants including newly-exposed rhminl_const, + rhminl_adj_land_const, rhminh_const. + +M src/physics/carma/models/cirrus/carma_cloudfraction.F90 +M src/physics/carma/models/cirrus_dust/carma_cloudfraction.F90 + - USE compute_cloud_fraction_two_moment for subroutines; + USE cldfrc2m shim for rh constants. + +If there were any failures reported from running test_driver.sh on any test +platform, and checkin with these failures has been OK'd by the gatekeeper, +then copy the lines from the td.*.status files for the failed tests to the +appropriate machine below. All failed tests must be justified. + +derecho/intel/aux_cam: + SMS_D_Ln9.f19_f19_mg17.FCARMA2000climo.derecho_intel.cam-outfrq9s_carma_rad_diag (Overall: DIFF) details: + - the testmod for this test changed to "outfrq9s_carma_rad_diag". + I cprnc'd the baseline files manually against + /glade/campaign/cesm/community/amwg/cam_baselines/cam6_4_182_intel/ + SMS_D_Ln9.f19_f19_mg17.FCARMA2000climo.derecho_intel.cam-outfrq9s + The cam.h0a file differs only in the field lists due to new variables introduced in rad_diag_2 and 3. + The other files are IDENTICAL. + +derecho/nvhpc/aux_cam: All PASS + +izumi/nag/aux_cam: + ERC_D_Ln9.f10_f10_mt232.FHIST_C5.izumi_nag.cam-outfrq3s_subcol (Overall: FAIL) details: + FAIL ERC_D_Ln9.f10_f10_mt232.FHIST_C5.izumi_nag.cam-outfrq3s_subcol COMPARE_base_rest + - pre-existing failure. + +izumi/gnu/aux_cam: All PASS + +Summarize any changes to answers: All B4B + +=============================================================== + Tag name: cam6_4_182 Originator(s): fvitt Date: 16 Jun 2026 diff --git a/src/atmos_phys b/src/atmos_phys index aa61cd8b28..23904a8784 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit aa61cd8b28678c78b4579c453ca3024de2ca0348 +Subproject commit 23904a8784f676b3191f1391c9f832c3654220f3 diff --git a/src/physics/cam/cldfrc2m.F90 b/src/physics/cam/cldfrc2m.F90 index 77a391fd35..863e26dfa1 100644 --- a/src/physics/cam/cldfrc2m.F90 +++ b/src/physics/cam/cldfrc2m.F90 @@ -1,15 +1,14 @@ module cldfrc2m -! cloud fraction calculations +! two-moment cloud fraction calculations - CAM interface +! Subroutine implementations are now in compute_cloud_fraction_two_moment module. +! this module only exists to provide public constants for the rest of CAM 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, & - svp_water_vect, svp_ice_vect use cam_logfile, only: iulog use cam_abortutils, only: endrun +use compute_cloud_fraction_two_moment, only: CAMstfrac implicit none private @@ -18,17 +17,14 @@ module cldfrc2m public :: & cldfrc2m_readnl, & cldfrc2m_init, & - astG_PDF_single, & - astG_PDF, & - astG_RHU_single, & - astG_RHU, & - aist_single, & - aist_vector, & CAMstfrac, & rhmini_const, & rhmaxi_const, & rhminis_const, & - rhmaxis_const + rhmaxis_const, & + rhminl_const, & + rhminl_adj_land_const, & + rhminh_const ! Namelist variables real(r8) :: cldfrc2m_rhmini ! Minimum rh for ice cloud fraction > 0. @@ -51,20 +47,15 @@ module cldfrc2m ! 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), protected :: rhminl_const ! Critical RH for low-level liquid stratus clouds +real(r8), protected :: rhminl_adj_land_const ! rhminl adjustment for snowfree land +real(r8), protected :: rhminh_const ! Critical RH for high-level liquid stratus clouds + +! Internal variables (used only by init) 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 @@ -123,11 +114,23 @@ end subroutine cldfrc2m_readnl subroutine cldfrc2m_init() use cloud_fraction, only: cldfrc_getparams + use physconst, only: rair + use compute_cloud_fraction_two_moment, only: compute_cloud_fraction_two_moment_init + + character(len=512) :: errmsg + integer :: errflg 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) + call compute_cloud_fraction_two_moment_init( & + masterproc, iulog, & + premit, premib, iceopt, icecrit, & + cldfrc2m_qist_min, cldfrc2m_qist_max, cldfrc2m_do_subgrid_growth, cldfrc2m_do_avg_aist_algs, & + rair, & + errmsg, errflg) + if( masterproc ) then write(iulog,*) 'cldfrc2m parameters:' write(iulog,*) ' rhminl = ', rhminl_const @@ -151,1032 +154,4 @@ 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) - if (rhmaxi .eq. rhmini) then - if (rhi .gt. rhmini) then - rhdif = 1._r8 - else - rhdif = 0._r8 - end if - else - rhdif = (rhi-rhmini) / (rhmaxi - rhmini) - end if - 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.cldfrc2m_qist_min) then - if (cldfrc2m_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._r8,min(1._r8,sqrt(aist*qi/cldfrc2m_qist_min))) - else - ! - ! Default for iceopt=5 - aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_min)) - end if - endif - !maximum - if (icimr.gt.cldfrc2m_qist_max) then - aist = max(0._r8,min(1._r8,qi/cldfrc2m_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(pcols) - 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(pcols) ! Liq sat vapor pressure - real(r8) esi(pcols) ! 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 - - 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), 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) - - if (present(rhmaxi_in)) rhmaxi = rhmaxi_in(i) - 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(i)/esl(i)))) - 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(i)/esi(i)) - if (rhmaxi .eq. rhmini) then - if (rhi .gt. rhmini) then - rhdif = 1._r8 - else - rhdif = 0._r8 - end if - else - rhdif = (rhi-rhmini) / (rhmaxi - rhmini) - end if - 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.cldfrc2m_qist_min) then - if (cldfrc2m_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._r8,min(1._r8,sqrt(aist*qi/cldfrc2m_qist_min))) - else - ! - ! Default for iceopt=5 - aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_min)) - end if - endif - !maximum - if (icimr.gt.cldfrc2m_qist_max) then - aist = max(0._r8,min(1._r8,qi/cldfrc2m_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 diff --git a/src/physics/cam/cldwat2m_macro.F90 b/src/physics/cam/cldwat2m_macro.F90 index f63ccc4e1c..ba2b4cb5e1 100644 --- a/src/physics/cam/cldwat2m_macro.F90 +++ b/src/physics/cam/cldwat2m_macro.F90 @@ -18,9 +18,9 @@ module cldwat2m_macro use wv_saturation, only: qsat_water, svp_water, svp_ice, qsat_ice use cam_logfile, only: iulog use ref_pres, only: top_lev=>trop_cloud_top_lev - use cldfrc2m, only: astG_PDF_single, astG_PDF, astG_RHU_single, & - astG_RHU, aist_single, aist_vector, & - rhmini_const, rhmaxi_const + use compute_cloud_fraction_two_moment, only: astG_PDF_single, astG_PDF, astG_RHU_single, & + astG_RHU, aist_single, aist_vector + use cldfrc2m, only: rhmini_const, rhmaxi_const implicit none private @@ -830,15 +830,29 @@ subroutine mmacro_pcond( lchnk , ncol , dt , p , d U(i,k) = qv(i,k)/qsat_b(i) U_nc(i,k) = U(i,k) enddo + !REMOVECAM: this is no longer needed when CAM is retired and pcols no longer exists + al_st_nc(:,k) = 0._r8 + G_nc(:,k) = 0._r8 + !REMOVECAM_END if( CAMstfrac ) then - call astG_RHU(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol,& - rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k)) + call astG_RHU(U_nc(:ncol,k),p(:ncol,k),qv(:ncol,k),landfrac(:ncol),snowh(:ncol), & + al_st_nc(:ncol,k),G_nc(:ncol,k),ncol, & + rhminl_in=rhminl_arr(:ncol,k), rhminl_adj_land_in=rhminl_adj_land_arr(:ncol,k), & + rhminh_in=rhminh_arr(:ncol,k)) else - call astG_PDF(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol,& - rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k)) + call astG_PDF(U_nc(:ncol,k),p(:ncol,k),qv(:ncol,k),landfrac(:ncol),snowh(:ncol), & + al_st_nc(:ncol,k),G_nc(:ncol,k),ncol, & + rhminl_in=rhminl_arr(:ncol,k), rhminl_adj_land_in=rhminl_adj_land_arr(:ncol,k), & + rhminh_in=rhminh_arr(:ncol,k)) endif - call aist_vector(qv(:,k),T(:,k),p(:,k),qi(:,k),ni(:,k),landfrac(:),snowh(:),ai_st_nc(:,k),ncol,& - rhmaxi_arr(:,k), rhmini_arr(:,k), rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k)) + !REMOVECAM: this is no longer needed when CAM is retired and pcols no longer exists + ai_st_nc(:,k) = 0._r8 + !REMOVECAM_END + call aist_vector(qv(:ncol,k),T(:ncol,k),p(:ncol,k),qi(:ncol,k),ni(:ncol,k), & + landfrac(:ncol),snowh(:ncol),ai_st_nc(:ncol,k),ncol, & + rhmaxi_in=rhmaxi_arr(:ncol,k), rhmini_in=rhmini_arr(:ncol,k), & + rhminl_in=rhminl_arr(:ncol,k), rhminl_adj_land_in=rhminl_adj_land_arr(:ncol,k), & + rhminh_in=rhminh_arr(:ncol,k)) ai_st(:ncol,k) = (1._r8-a_cu(:ncol,k))*ai_st_nc(:ncol,k) al_st(:ncol,k) = (1._r8-a_cu(:ncol,k))*al_st_nc(:ncol,k) @@ -1347,15 +1361,25 @@ subroutine instratus_condensate( lchnk, ncol, k, & call qsat_water(T0_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol) U0_in(:ncol) = qv0_in(:ncol)/qsat_in(:ncol) + al0_st_nc_in(:) = 0._r8 + G0_nc_in(:) = 0._r8 if( CAMstfrac ) then - call astG_RHU(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol,& - rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:)) + call astG_RHU(U0_in(:ncol),p_in(:ncol),qv0_in(:ncol),landfrac(:ncol),snowh(:ncol), & + al0_st_nc_in(:ncol),G0_nc_in(:ncol),ncol, & + rhminl_in=rhminl_in(:ncol), rhminl_adj_land_in=rhminl_adj_land_in(:ncol), & + rhminh_in=rhminh_in(:ncol)) else - call astG_PDF(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol,& - rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:)) + call astG_PDF(U0_in(:ncol),p_in(:ncol),qv0_in(:ncol),landfrac(:ncol),snowh(:ncol), & + al0_st_nc_in(:ncol),G0_nc_in(:ncol),ncol, & + rhminl_in=rhminl_in(:ncol), rhminl_adj_land_in=rhminl_adj_land_in(:ncol), & + rhminh_in=rhminh_in(:ncol)) endif - call aist_vector(qv0_in(:),T0_in(:),p_in(:),qi0_in(:),ni0_in(:),landfrac(:),snowh(:),ai0_st_nc_in(:),ncol,& - rhmaxi_in(:), rhmini_in(:), rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:)) + ai0_st_nc_in(:) = 0._r8 + call aist_vector(qv0_in(:ncol),T0_in(:ncol),p_in(:ncol),qi0_in(:ncol),ni0_in(:ncol), & + landfrac(:ncol),snowh(:ncol),ai0_st_nc_in(:ncol),ncol, & + rhmaxi_in=rhmaxi_in(:ncol), rhmini_in=rhmini_in(:ncol), & + rhminl_in=rhminl_in(:ncol), rhminl_adj_land_in=rhminl_adj_land_in(:ncol), & + rhminh_in=rhminh_in(:ncol)) do i = 1, ncol @@ -1410,13 +1434,14 @@ subroutine instratus_condensate( lchnk, ncol, k, & U0_nc = U0 if( CAMstfrac ) then call astG_RHU_single(U0_nc, p, qv0, landfrac(i), snowh(i), al0_st_nc, G0_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) else call astG_PDF_single(U0_nc, p, qv0, landfrac(i), snowh(i), al0_st_nc, G0_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) endif call aist_single(qv0,T0,p,qi0,landfrac(i),snowh(i),ai0_st_nc,& - rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh) + rhmaxi=rhmaxi, rhmini=rhmini, rhminl=rhminl, & + rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) ai0_st = (1._r8-a_dc-a_sc)*ai0_st_nc al0_st = (1._r8-a_dc-a_sc)*al0_st_nc a0_st = max(ai0_st,al0_st) @@ -1469,10 +1494,10 @@ subroutine instratus_condensate( lchnk, ncol, k, & U_nc = U if( CAMstfrac ) then call astG_RHU_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) else call astG_PDF_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) endif al_st = (1._r8-a_dc-a_sc)*al_st_nc caseid = 0 @@ -1584,17 +1609,18 @@ subroutine instratus_condensate( lchnk, ncol, k, & if( idxmod .eq. 1 ) then call aist_single(qv,T,p,qi,landfrac(i),snowh(i),ai_st_nc,& - rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh) + rhmaxi=rhmaxi, rhmini=rhmini, rhminl=rhminl, & + rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) ai_st = (1._r8-a_dc-a_sc)*ai_st_nc call qsat_water(T, p, es, qs) U = (qv/qs) U_nc = U if( CAMstfrac ) then call astG_RHU_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) else call astG_PDF_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) endif al_st = (1._r8-a_dc-a_sc)*al_st_nc else @@ -1908,10 +1934,10 @@ subroutine funcd_instratus( T, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, & U_nc = U if( CAMstfrac ) then call astG_RHU_single(U_nc, p, qv, landfrac, snowh, al_st_nc, G_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) else call astG_PDF_single(U_nc, p, qv, landfrac, snowh, al_st_nc, G_nc, & - rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh) + rhminl=rhminl, rhminl_adj_land=rhminl_adj_land, rhminh=rhminh) endif al_st = (1._r8-a_dc-a_sc)*al_st_nc dUdt = -(alpha*dqcncdt+beta) diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90 index cbebae4323..c9b05212be 100644 --- a/src/physics/cam/clubb_intr.F90 +++ b/src/physics/cam/clubb_intr.F90 @@ -57,7 +57,7 @@ module clubb_intr clubb_readnl, clubb_init_cnst, clubb_implements_cnst #ifdef CLUBB_SGS - + ! NOTE: the only reason for anything in this section being set to public is for use with SILHS public :: stats_init_clubb, stats_end_timestep_clubb @@ -141,7 +141,7 @@ module clubb_intr rtpthlp_const = 0.01_r8 ! Constant to add to rtpthlp when moments are advected real(r8), parameter :: unset_r8 = huge(1.0_r8) - + integer, parameter :: unset_i = huge(1) ! Commonly used temperature for the melting temp of ice crystals [K] @@ -162,11 +162,11 @@ module clubb_intr clubb_l_intr_sfc_flux_smooth = .false. ! Add a locally calculated roughness to upwp and vpwp sfc fluxes logical :: & - clubb_l_ascending_grid = .false. ! Run clubb in ascending mode, which is opposite of the + clubb_l_ascending_grid = .false. ! Run clubb in ascending mode, which is opposite of the ! cam grid the rest of this code uses, thus it requires ! an expensive array flipping step before calling advance_clubb_core. ! This is mainly for testing, it should not significantly change answers - + logical :: lq(pcnst) logical :: do_rainturb logical :: clubb_do_adv @@ -278,7 +278,7 @@ module clubb_intr clubb_grid_adapt_in_time_method = unset_i, & ! Specifier for how the grid density method should ! be constructed if the grid should be adapted over time ! (set to 0 for no adaptation) - clubb_fill_holes_type = unset_i ! Option for which type of hole filler to use in the + clubb_fill_holes_type = unset_i ! Option for which type of hole filler to use in the ! fill_holes_vertical procedure @@ -410,21 +410,21 @@ module clubb_intr ! Indices for physics buffer (pbuf) ! ! ------------------------------------------------------------ ! integer :: & - wp2_idx, & ! vertical velocity variances - wp3_idx, & ! third moment of vertical velocity - wpthlp_idx, & ! turbulent flux of thetal - wprtp_idx, & ! turbulent flux of total water - rtpthlp_idx, & ! covariance of thetal and rt - rtp2_idx, & ! variance of total water - thlp2_idx, & ! variance of thetal - rtp3_idx, & ! total water 3rd order - thlp3_idx, & ! thetal 3rd order - up2_idx, & ! variance of east-west wind - vp2_idx, & ! variance of north-south wind - up3_idx, & ! east-west wind 3rd order - vp3_idx, & ! north-south wind 3rd order - upwp_idx, & ! east-west momentum flux - vpwp_idx, & ! north-south momentum flux + wp2_idx, & ! vertical velocity variances + wp3_idx, & ! third moment of vertical velocity + wpthlp_idx, & ! turbulent flux of thetal + wprtp_idx, & ! turbulent flux of total water + rtpthlp_idx, & ! covariance of thetal and rt + rtp2_idx, & ! variance of total water + thlp2_idx, & ! variance of thetal + rtp3_idx, & ! total water 3rd order + thlp3_idx, & ! thetal 3rd order + up2_idx, & ! variance of east-west wind + vp2_idx, & ! variance of north-south wind + up3_idx, & ! east-west wind 3rd order + vp3_idx, & ! north-south wind 3rd order + upwp_idx, & ! east-west momentum flux + vpwp_idx, & ! north-south momentum flux wpthvp_idx, & ! buoyancy flux wp2thvp_idx, & ! second order buoyancy term wp2up_idx, & ! w'^2 u' @@ -440,7 +440,7 @@ module clubb_intr wpvp2_idx, & ! w'v'^2 wp2up2_idx, & ! w'^2 u'^2 wp2vp2_idx, & ! w'^2 v'^2 - cld_idx, & ! Cloud fraction + cld_idx, & ! Cloud fraction concld_idx, & ! Convective cloud fraction ast_idx, & ! Stratiform cloud fraction alst_idx, & ! Liquid stratiform cloud fraction @@ -449,9 +449,9 @@ module clubb_intr qist_idx, & ! Physical in-cloud IWC dp_frac_idx, & ! deep convection cloud fraction sh_frac_idx, & ! shallow convection cloud fraction - kvh_idx, & ! CLUBB eddy diffusivity on thermo levels + kvh_idx, & ! CLUBB eddy diffusivity on thermo levels pblh_idx, & ! PBL pbuf - icwmrdp_idx, & ! In cloud mixing ratio for deep convection + icwmrdp_idx, & ! In cloud mixing ratio for deep convection tke_idx, & ! turbulent kinetic energy tpert_idx, & ! temperature perturbation from PBL fice_idx, & ! fice_idx index in physics buffer @@ -487,7 +487,7 @@ module clubb_intr rtpthlp_mc_zt_idx ! added pbuf fields for clubb to have restart bfb when ipdf_call_placement=2 - integer :: & + integer :: & pdf_zm_w_1_idx, & pdf_zm_w_2_idx, & pdf_zm_varnce_w_1_idx, & @@ -546,7 +546,7 @@ subroutine clubb_register_cam( ) call cnst_add(trim(cnst_names(8)),0._r8,0._r8,0._r8,ixup2,longname='CLUBB 2nd moment u wind',cam_outfld=.false.) call cnst_add(trim(cnst_names(9)),0._r8,0._r8,0._r8,ixvp2,longname='CLUBB 2nd moment v wind',cam_outfld=.false.) end if - + ! Determine number of vertical levels used in clubb, thermo variables are nzt_clubb ! and momentum variables are nzm_clubb nzt_clubb = pver + 1 - top_lev @@ -1363,73 +1363,73 @@ subroutine clubb_readnl(nlfile) if ( clubb_detphase_lowtemp == unset_r8 ) call endrun( sub//": FATAL: clubb_detphase_lowtemp not set") if ( clubb_detphase_lowtemp >= meltpt_temp ) call endrun( sub//": ERROR: clubb_detphase_lowtemp must be less than 268.15 K") - call initialize_clubb_config_flags_type_api( clubb_iiPDF_type, & ! In - clubb_ipdf_call_placement, & ! In - clubb_penta_solve_method, & ! In - clubb_tridiag_solve_method, & ! In - clubb_saturation_equation, & ! In - clubb_grid_remap_method, & ! In - clubb_grid_adapt_in_time_method, & ! In - clubb_fill_holes_type, & ! In - clubb_l_use_precip_frac, & ! In - clubb_l_predict_upwp_vpwp, & ! In + call initialize_clubb_config_flags_type_api( clubb_iiPDF_type, & ! In + clubb_ipdf_call_placement, & ! In + clubb_penta_solve_method, & ! In + clubb_tridiag_solve_method, & ! In + clubb_saturation_equation, & ! In + clubb_grid_remap_method, & ! In + clubb_grid_adapt_in_time_method, & ! In + clubb_fill_holes_type, & ! In + clubb_l_use_precip_frac, & ! In + clubb_l_predict_upwp_vpwp, & ! In clubb_l_ho_nontrad_coriolis, & ! In - clubb_l_ho_trad_coriolis, & ! In - clubb_l_min_wp2_from_corr_wx, & ! In - clubb_l_min_xp2_from_corr_wx, & ! In - clubb_l_C2_cloud_frac, & ! In - clubb_l_diffuse_rtm_and_thlm, & ! In - clubb_l_stability_correct_Kh_N2_zm, & ! In - clubb_l_calc_thlp2_rad, & ! In - clubb_l_upwind_xpyp_ta, & ! In - clubb_l_upwind_xm_ma, & ! In - clubb_l_uv_nudge, & ! In - clubb_l_rtm_nudge, & ! In - clubb_l_tke_aniso, & ! In - clubb_l_vert_avg_closure, & ! In - clubb_l_trapezoidal_rule_zt, & ! In - clubb_l_trapezoidal_rule_zm, & ! In - clubb_l_call_pdf_closure_twice, & ! In - clubb_l_standard_term_ta, & ! In - clubb_l_partial_upwind_wp3, & ! In - clubb_l_godunov_upwind_wpxp_ta, & ! In - clubb_l_godunov_upwind_xpyp_ta, & ! In - clubb_l_use_cloud_cover, & ! In - clubb_l_diagnose_correlations, & ! In - clubb_l_calc_w_corr, & ! In - clubb_l_const_Nc_in_cloud, & ! In - clubb_l_fix_w_chi_eta_correlations, & ! In - clubb_l_stability_correct_tau_zm, & ! In - clubb_l_damp_wp2_using_em, & ! In - clubb_l_do_expldiff_rtm_thlm, & ! In - clubb_l_Lscale_plume_centered, & ! In - clubb_l_diag_Lscale_from_tau, & ! In - clubb_l_use_C7_Richardson, & ! In - clubb_l_use_C11_Richardson, & ! In - clubb_l_use_shear_Richardson, & ! In - clubb_l_brunt_vaisala_freq_moist, & ! In - clubb_l_use_thvm_in_bv_freq, & ! In - clubb_l_rcm_supersat_adj, & ! In - clubb_l_damp_wp3_Skw_squared, & ! In - clubb_l_prescribed_avg_deltaz, & ! In - clubb_l_lmm_stepping, & ! In - clubb_l_e3sm_config, & ! In - clubb_l_vary_convect_depth, & ! In - clubb_l_use_tke_in_wp3_pr_turb_term, & ! In - clubb_l_use_tke_in_wp2_wp3_K_dfsn, & ! In - clubb_l_use_wp3_lim_with_smth_Heaviside, & ! In - clubb_l_smooth_Heaviside_tau_wpxp, & ! In - clubb_l_modify_limiters_for_cnvg_test, & ! In - clubb_l_enable_relaxed_clipping, & ! In - clubb_l_linearize_pbl_winds, & ! In - clubb_l_mono_flux_lim_thlm, & ! In - clubb_l_mono_flux_lim_rtm, & ! In - clubb_l_mono_flux_lim_um, & ! In - clubb_l_mono_flux_lim_vm, & ! In - clubb_l_mono_flux_lim_spikefix, & ! In - clubb_l_host_applies_sfc_fluxes, & ! In - clubb_l_wp2_fill_holes_tke, & ! In - clubb_l_add_dycore_grid, & ! In + clubb_l_ho_trad_coriolis, & ! In + clubb_l_min_wp2_from_corr_wx, & ! In + clubb_l_min_xp2_from_corr_wx, & ! In + clubb_l_C2_cloud_frac, & ! In + clubb_l_diffuse_rtm_and_thlm, & ! In + clubb_l_stability_correct_Kh_N2_zm, & ! In + clubb_l_calc_thlp2_rad, & ! In + clubb_l_upwind_xpyp_ta, & ! In + clubb_l_upwind_xm_ma, & ! In + clubb_l_uv_nudge, & ! In + clubb_l_rtm_nudge, & ! In + clubb_l_tke_aniso, & ! In + clubb_l_vert_avg_closure, & ! In + clubb_l_trapezoidal_rule_zt, & ! In + clubb_l_trapezoidal_rule_zm, & ! In + clubb_l_call_pdf_closure_twice, & ! In + clubb_l_standard_term_ta, & ! In + clubb_l_partial_upwind_wp3, & ! In + clubb_l_godunov_upwind_wpxp_ta, & ! In + clubb_l_godunov_upwind_xpyp_ta, & ! In + clubb_l_use_cloud_cover, & ! In + clubb_l_diagnose_correlations, & ! In + clubb_l_calc_w_corr, & ! In + clubb_l_const_Nc_in_cloud, & ! In + clubb_l_fix_w_chi_eta_correlations, & ! In + clubb_l_stability_correct_tau_zm, & ! In + clubb_l_damp_wp2_using_em, & ! In + clubb_l_do_expldiff_rtm_thlm, & ! In + clubb_l_Lscale_plume_centered, & ! In + clubb_l_diag_Lscale_from_tau, & ! In + clubb_l_use_C7_Richardson, & ! In + clubb_l_use_C11_Richardson, & ! In + clubb_l_use_shear_Richardson, & ! In + clubb_l_brunt_vaisala_freq_moist, & ! In + clubb_l_use_thvm_in_bv_freq, & ! In + clubb_l_rcm_supersat_adj, & ! In + clubb_l_damp_wp3_Skw_squared, & ! In + clubb_l_prescribed_avg_deltaz, & ! In + clubb_l_lmm_stepping, & ! In + clubb_l_e3sm_config, & ! In + clubb_l_vary_convect_depth, & ! In + clubb_l_use_tke_in_wp3_pr_turb_term, & ! In + clubb_l_use_tke_in_wp2_wp3_K_dfsn, & ! In + clubb_l_use_wp3_lim_with_smth_Heaviside, & ! In + clubb_l_smooth_Heaviside_tau_wpxp, & ! In + clubb_l_modify_limiters_for_cnvg_test, & ! In + clubb_l_enable_relaxed_clipping, & ! In + clubb_l_linearize_pbl_winds, & ! In + clubb_l_mono_flux_lim_thlm, & ! In + clubb_l_mono_flux_lim_rtm, & ! In + clubb_l_mono_flux_lim_um, & ! In + clubb_l_mono_flux_lim_vm, & ! In + clubb_l_mono_flux_lim_spikefix, & ! In + clubb_l_host_applies_sfc_fluxes, & ! In + clubb_l_wp2_fill_holes_tke, & ! In + clubb_l_add_dycore_grid, & ! In clubb_config_flags ) ! Out #endif @@ -1513,7 +1513,7 @@ subroutine clubb_ini_cam(pbuf_ini) type(err_info_type) :: & err_info ! err_info struct used in CLUBB containing err_code and err_header - + integer :: i, j, k, l ! Indices integer :: nmodes, nspec, m integer :: ixq, ixcldice, ixcldliq, ixnumliq, ixnumice @@ -2123,7 +2123,9 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & init_err_info_api, & cleanup_err_info_api - use cldfrc2m, only: aist_vector, rhmini_const, rhmaxi_const, rhminis_const, rhmaxis_const + use compute_cloud_fraction_two_moment, only: aist_vector + use cldfrc2m, only: rhmini_const, rhmaxi_const, rhminis_const, rhmaxis_const, & + rhminl_const, rhminl_adj_land_const, rhminh_const use cam_history, only: outfld use macrop_driver, only: liquid_macro_tend @@ -2165,26 +2167,26 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & ! Pointers for pbuf ! ! ---------------------------------------------------- ! - real(r8), pointer, dimension(:,:) :: wp2_pbuf ! vertical velocity variance [m^2/s^2] - real(r8), pointer, dimension(:,:) :: wp3_pbuf ! third moment of vertical velocity [m^3/s^3] - real(r8), pointer, dimension(:,:) :: wpthlp_pbuf ! turbulent flux of thetal [m/s K] - real(r8), pointer, dimension(:,:) :: wprtp_pbuf ! turbulent flux of moisture [m/s kg/kg] - real(r8), pointer, dimension(:,:) :: rtpthlp_pbuf ! covariance of thetal and qt [kg/kg K] - real(r8), pointer, dimension(:,:) :: rtp2_pbuf ! moisture variance [kg^2/kg^2] - real(r8), pointer, dimension(:,:) :: thlp2_pbuf ! temperature variance [K^2] - real(r8), pointer, dimension(:,:) :: rtp3_pbuf ! moisture 3rd order [kg^3/kg^3] - real(r8), pointer, dimension(:,:) :: thlp3_pbuf ! temperature 3rd order [K^3] - real(r8), pointer, dimension(:,:) :: up2_pbuf ! east-west wind variance [m^2/s^2] - real(r8), pointer, dimension(:,:) :: vp2_pbuf ! north-south wind variance [m^2/s^2] - real(r8), pointer, dimension(:,:) :: up3_pbuf ! east-west wind 3rd order [m^3/s^3] - real(r8), pointer, dimension(:,:) :: vp3_pbuf ! north-south wind 3rd order [m^3/s^3] - real(r8), pointer, dimension(:,:) :: upwp_pbuf ! east-west momentum flux [m^2/s^2] - real(r8), pointer, dimension(:,:) :: vpwp_pbuf ! north-south momentum flux [m^2/s^2] - real(r8), pointer, dimension(:,:) :: wpthvp_pbuf ! w'th_v' (momentum levels) [m/s K] - real(r8), pointer, dimension(:,:) :: wp2thvp_pbuf ! w'^2 th_v' (thermodynamic levels) [m^2/s^2 K] - real(r8), pointer, dimension(:,:) :: wp2up_pbuf ! w'^2 u' (thermodynamic levels) [m^3/s^3] - real(r8), pointer, dimension(:,:) :: rtpthvp_pbuf ! r_t'th_v' (momentum levels) [kg/kg K] - real(r8), pointer, dimension(:,:) :: thlpthvp_pbuf ! th_l'th_v' (momentum levels) [K^2] + real(r8), pointer, dimension(:,:) :: wp2_pbuf ! vertical velocity variance [m^2/s^2] + real(r8), pointer, dimension(:,:) :: wp3_pbuf ! third moment of vertical velocity [m^3/s^3] + real(r8), pointer, dimension(:,:) :: wpthlp_pbuf ! turbulent flux of thetal [m/s K] + real(r8), pointer, dimension(:,:) :: wprtp_pbuf ! turbulent flux of moisture [m/s kg/kg] + real(r8), pointer, dimension(:,:) :: rtpthlp_pbuf ! covariance of thetal and qt [kg/kg K] + real(r8), pointer, dimension(:,:) :: rtp2_pbuf ! moisture variance [kg^2/kg^2] + real(r8), pointer, dimension(:,:) :: thlp2_pbuf ! temperature variance [K^2] + real(r8), pointer, dimension(:,:) :: rtp3_pbuf ! moisture 3rd order [kg^3/kg^3] + real(r8), pointer, dimension(:,:) :: thlp3_pbuf ! temperature 3rd order [K^3] + real(r8), pointer, dimension(:,:) :: up2_pbuf ! east-west wind variance [m^2/s^2] + real(r8), pointer, dimension(:,:) :: vp2_pbuf ! north-south wind variance [m^2/s^2] + real(r8), pointer, dimension(:,:) :: up3_pbuf ! east-west wind 3rd order [m^3/s^3] + real(r8), pointer, dimension(:,:) :: vp3_pbuf ! north-south wind 3rd order [m^3/s^3] + real(r8), pointer, dimension(:,:) :: upwp_pbuf ! east-west momentum flux [m^2/s^2] + real(r8), pointer, dimension(:,:) :: vpwp_pbuf ! north-south momentum flux [m^2/s^2] + real(r8), pointer, dimension(:,:) :: wpthvp_pbuf ! w'th_v' (momentum levels) [m/s K] + real(r8), pointer, dimension(:,:) :: wp2thvp_pbuf ! w'^2 th_v' (thermodynamic levels) [m^2/s^2 K] + real(r8), pointer, dimension(:,:) :: wp2up_pbuf ! w'^2 u' (thermodynamic levels) [m^3/s^3] + real(r8), pointer, dimension(:,:) :: rtpthvp_pbuf ! r_t'th_v' (momentum levels) [kg/kg K] + real(r8), pointer, dimension(:,:) :: thlpthvp_pbuf ! th_l'th_v' (momentum levels) [K^2] real(r8), pointer, dimension(:,:) :: pdf_zm_w_1_pbuf ! work pointer for pdf_params_zm real(r8), pointer, dimension(:,:) :: pdf_zm_w_2_pbuf ! work pointer for pdf_params_zm real(r8), pointer, dimension(:,:) :: pdf_zm_varnce_w_1_pbuf ! work pointer for pdf_params_zm @@ -2200,15 +2202,15 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & real(r8), pointer, dimension(:,:) :: wpvp2_pbuf ! w'v'^2 (thermodynamic levels) real(r8), pointer, dimension(:,:) :: wp2up2_pbuf ! w'^2 u'^2 (momentum levels) real(r8), pointer, dimension(:,:) :: wp2vp2_pbuf ! w'^2 v'^2 (momentum levels) - real(r8), pointer, dimension(:,:) :: cld_pbuf ! cloud fraction [fraction] - real(r8), pointer, dimension(:,:) :: concld_pbuf ! convective cloud fraction [fraction] - real(r8), pointer, dimension(:,:) :: ast_pbuf ! stratiform cloud fraction [fraction] - real(r8), pointer, dimension(:,:) :: alst_pbuf ! liquid stratiform cloud fraction [fraction] - real(r8), pointer, dimension(:,:) :: aist_pbuf ! ice stratiform cloud fraction [fraction] - real(r8), pointer, dimension(:,:) :: qlst_pbuf ! Physical in-stratus LWC [kg/kg] - real(r8), pointer, dimension(:,:) :: qist_pbuf ! Physical in-stratus IWC [kg/kg] - real(r8), pointer, dimension(:,:) :: deepcu_pbuf ! deep convection cloud fraction [fraction] - real(r8), pointer, dimension(:,:) :: shalcu_pbuf ! shallow convection cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: cld_pbuf ! cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: concld_pbuf ! convective cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: ast_pbuf ! stratiform cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: alst_pbuf ! liquid stratiform cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: aist_pbuf ! ice stratiform cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: qlst_pbuf ! Physical in-stratus LWC [kg/kg] + real(r8), pointer, dimension(:,:) :: qist_pbuf ! Physical in-stratus IWC [kg/kg] + real(r8), pointer, dimension(:,:) :: deepcu_pbuf ! deep convection cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: shalcu_pbuf ! shallow convection cloud fraction [fraction] real(r8), pointer, dimension(:,:) :: khzm_pbuf ! CLUBB's eddy diffusivity of heat/moisture on momentum levels [m^2/s] real(r8), pointer, dimension(:) :: pblh_pbuf ! planetary boundary layer height [m] real(r8), pointer, dimension(:,:) :: tke_pbuf ! turbulent kinetic energy [m^2/s^2] @@ -2274,9 +2276,9 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & real(r8), dimension(state%ncol) :: & deltaz, & - fcor, & ! Coriolis forcing [s^-1] - fcor_y, & ! Non-traditional coriolis forcing [s^-1] - sfc_elevation, & ! Elevation of ground [m AMSL][m] + fcor, & ! Coriolis forcing [s^-1] + fcor_y, & ! Non-traditional coriolis forcing [s^-1] + sfc_elevation, & ! Elevation of ground [m AMSL][m] wpthlp_sfc, & ! w' theta_l' at surface [(m K)/s] wprtp_sfc, & ! w' r_t' at surface [(kg m)/( kg s)] upwp_sfc, & ! u'w' at surface [m^2/s^2] @@ -2293,28 +2295,28 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & wpedsclrp_sfc ! Eddy-scalar flux at surface [{units vary} m/s] real(r8), dimension(state%ncol,nzt_clubb) :: & - rtm, & ! mean moisture mixing ratio [kg/kg] - thlm, & ! mean temperature [K] + rtm, & ! mean moisture mixing ratio [kg/kg] + thlm, & ! mean temperature [K] rcm, & ! CLUBB cloud water mixing ratio [kg/kg] - um, & ! mean east-west wind [m/s] - vm, & ! mean north-south wind [m/s] + um, & ! mean east-west wind [m/s] + vm, & ! mean north-south wind [m/s] thlm_forcing, & ! theta_l forcing (thermodynamic levels) [K/s] rtm_forcing, & ! r_t forcing (thermodynamic levels) [(kg/kg)/s] - um_forcing, & ! u wind forcing (thermodynamic levels) [m/s/s] - vm_forcing, & ! v wind forcing (thermodynamic levels) [m/s/s] - wm_zt, & ! w mean wind component on thermo. levels [m/s] + um_forcing, & ! u wind forcing (thermodynamic levels) [m/s/s] + vm_forcing, & ! v wind forcing (thermodynamic levels) [m/s/s] + wm_zt, & ! w mean wind component on thermo. levels [m/s] rtm_ref, & ! Initial profile of rtm [kg/kg] thlm_ref, & ! Initial profile of thlm [K] um_ref, & ! Initial profile of um [m/s] vm_ref, & ! Initial profile of vm [m/s] ug, & ! U geostrophic wind [m/s] vg, & ! V geostrophic wind [m/s] - p_in_Pa, & ! Air pressure (thermodynamic levels) [Pa] + p_in_Pa, & ! Air pressure (thermodynamic levels) [Pa] rho_zt, & ! Air density on thermo levels [kg/m^3] exner, & ! Exner function (thermodynamic levels) [-] - rho_ds_zt, & ! Dry, static density on thermodynamic levels [kg/m^3] - invrs_rho_ds_zt, & ! Inv. dry, static density on thermo. levels [m^3/kg] - thv_ds_zt, & ! Dry, base-state theta_v on thermo. levels [K] + rho_ds_zt, & ! Dry, static density on thermodynamic levels [kg/m^3] + invrs_rho_ds_zt, & ! Inv. dry, static density on thermo. levels [m^3/kg] + thv_ds_zt, & ! Dry, base-state theta_v on thermo. levels [K] rfrzm, & rvm, & ! water vapor mixing ratio [kg/kg] rtp2_zt, & ! CLUBB R-tot variance on thermo levs @@ -2333,7 +2335,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & pre, & ! input for precip evaporation qrl_clubb, & qclvar, & ! cloud water variance [kg^2/kg^2] - zt_g, & ! Thermodynamic grid of CLUBB [m] + zt_g, & ! Thermodynamic grid of CLUBB [m] Lscale, & dz_g, & ! thickness of layer [m] invrs_dz_g, & ! Inverse of layer thickness [1/m] @@ -2349,11 +2351,11 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & rtp2_forcing, & thlp2_forcing, & rtpthlp_forcing, & - wm_zm, & ! w mean wind component on momentum levels [m/s] + wm_zm, & ! w mean wind component on momentum levels [m/s] rho_zm, & ! Air density on momentum levels [kg/m^3] - rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3] - invrs_rho_ds_zm, & ! Inv. dry, static density on momentum levels [m^3/kg] - thv_ds_zm, & ! Dry, base-state theta_v on momentum levels [K] + rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3] + invrs_rho_ds_zm, & ! Inv. dry, static density on momentum levels [m^3/kg] + thv_ds_zm, & ! Dry, base-state theta_v on momentum levels [K] upwp_pert, & ! Perturbed u'w' [m^2/s^2] vpwp_pert, & ! Perturbed v'w' [m^2/s^2] khzm, & ! Eddy diffusivity of heat/moisture on momentum levels [m^2/s] @@ -2365,7 +2367,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & wprtp_mc, & wpthlp_mc, & rtpthlp_mc, & - zi_g, & ! Momentum grid of CLUBB [m] + zi_g, & ! Momentum grid of CLUBB [m] ! MF Plume mf_dry_a, mf_moist_a, & @@ -2385,12 +2387,12 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & rtm_zm, thlm_zm, & ! momentum grid kappa_zm, p_in_Pa_zm, & ! momentum grid invrs_exner_zm ! momentum grid - + real(r8), dimension(state%ncol,nzt_clubb,sclr_dim) :: & sclrm_forcing, & ! Passive scalar forcing [{units vary}/s] sclrm, & ! Passive scalar mean (thermo. levels) [units vary] sclrp3 ! sclr'^3 (thermo. levels) [{units vary}^3] - + real(r8), dimension(state%ncol,nzm_clubb,sclr_dim) :: & sclrp2, & ! sclr'^2 (momentum levels) [{units vary}^2] sclrprtp, & ! sclr'rt' (momentum levels) [{units vary} (kg/kg)] @@ -2460,8 +2462,8 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & rtp2_zt_output, & ! CLUBB R-tot variance on thermo levs [kg^2/kg^2] wp3_output, & ! wp3 output [m^3/s^3] thl2_zt_output, & ! CLUBB Theta-l variance on thermo levs - wp2_zt_output, & - rcm_in_layer_output, & ! CLUBB in-cloud liquid water mixing ratio [kg/kg] + wp2_zt_output, & + rcm_in_layer_output, & ! CLUBB in-cloud liquid water mixing ratio [kg/kg] pdfp_rtp2_output, & ! Calculated R-tot variance from pdf_params [kg^2/kg^2] wm_zt_output, & ! CLUBB mean W on thermo levs output [m/s] rcm_output, & @@ -2484,6 +2486,11 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & kinwat, & ! Kinematic water vapor flux [m/s] dummy2, & ! dummy variable [units vary] dummy3 ! dummy variable [units vary] + real(r8) :: rhmini_default(pcols) + real(r8) :: rhmaxi_default(pcols) + real(r8) :: rhminl_arr(pcols) + real(r8) :: rhminl_adj_land_arr(pcols) + real(r8) :: rhminh_arr(pcols) real(r8), dimension(pcols,pver) :: & invrs_cpairv, & @@ -2494,7 +2501,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & qvtend, & qctend, & inctend, & - thv, & ! virtual potential temperature [K] + thv, & ! virtual potential temperature [K] th ! potential temperature [K] real(r8), dimension(pcols,nzt_clubb) :: & @@ -2510,7 +2517,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & mixt_frac_max_mag, & dtime, & ! CLUBB time step [s] ubar, & ! surface wind [m/s] - ustar, & ! surface stress [m/s] + ustar, & ! surface stress [m/s] bflx22, & ! Variable for buoyancy flux for pbl [K m/s] zo, & ! roughness height [m] relvarmax, & @@ -2794,7 +2801,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & endif !----------------------------------------- BEGIN GPU SECTION ----------------------------------------- - ! everything within should be functional with the OpenACC code, or be prevented from running + ! everything within should be functional with the OpenACC code, or be prevented from running ! with using OpenACC, see the "ifdef _OPENACC" section above for restriction examples call t_stopf('clubb_tend_cam:non_acc_region') @@ -2843,7 +2850,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & !$acc pdf_params_chnk(lchnk)%covar_chi_eta_1, pdf_params_chnk(lchnk)%covar_chi_eta_2, & !$acc pdf_params_chnk(lchnk)%corr_w_chi_1, pdf_params_chnk(lchnk)%corr_w_chi_2, & !$acc pdf_params_chnk(lchnk)%corr_w_eta_1, pdf_params_chnk(lchnk)%corr_w_eta_2, & - !$acc pdf_params_chnk(lchnk)%corr_chi_eta_1, pdf_params_chnk(lchnk)%corr_chi_eta_2, & + !$acc pdf_params_chnk(lchnk)%corr_chi_eta_1, pdf_params_chnk(lchnk)%corr_chi_eta_2, & !$acc pdf_params_chnk(lchnk)%rsatl_1, pdf_params_chnk(lchnk)%rsatl_2, & !$acc pdf_params_chnk(lchnk)%rc_1, pdf_params_chnk(lchnk)%rc_2, & !$acc pdf_params_chnk(lchnk)%cloud_frac_1, pdf_params_chnk(lchnk)%cloud_frac_2, & @@ -2915,7 +2922,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & vm_pert(i,k) = 0.0_r8 end do end do - + !$acc parallel loop gang vector collapse(2) default(present) do k = 1, nzm_clubb do i = 1, ncol @@ -3091,13 +3098,13 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & if ( clubb_do_adv ) then if (macmic_it == 1) then - + ! Note that some of the moments below can be positive or negative. ! Remove a constant that was added to prevent dynamics from clipping ! them to prevent dynamics from making them positive. do k = 1, nzm_clubb do i = 1, ncol - k_cam = top_lev - 1 + k + k_cam = top_lev - 1 + k rtpthlp_pbuf(i,k) = state_loc%q(i,k_cam,ixrtpthlp) - ( rtpthlp_const * apply_const ) wpthlp_pbuf(i,k) = state_loc%q(i,k_cam, ixwpthlp) - ( wpthlp_const * apply_const ) wprtp_pbuf(i,k) = state_loc%q(i,k_cam, ixwprtp) - ( wprtp_const * apply_const ) @@ -3141,7 +3148,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & do k = 1, nzt_clubb do i = 1, ncol - k_cam = top_lev - 1 + k + k_cam = top_lev - 1 + k ! Define the CLUBB thermodynamic grid (in units of m) zt_g(i,k) = state_loc%zm(i,k_cam) - state_loc%zi(i,pverp) @@ -3161,10 +3168,10 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & do k = 1, nzt_clubb do i = 1, ncol - k_cam = top_lev - 1 + k + k_cam = top_lev - 1 + k p_in_Pa(i,k) = state_loc%pmid(i,k_cam) - + ! Compute inverse exner function consistent with CLUBB's definition, which uses a constant ! surface pressure. CAM's exner (in state) does not. Therefore, for consistent ! treatment with CLUBB code, anytime exner is needed to treat CLUBB variables @@ -3189,7 +3196,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & do k = 1, nzt_clubb do i = 1, ncol - k_cam = top_lev - 1 + k + k_cam = top_lev - 1 + k ! Compute mean w wind on thermo grid, convert from omega to w wm_zt(i,k) = -1._r8 * ( state_loc%omega(i,k_cam) - state_loc%omega(i,pver) ) / ( rho_zt(i,k) * gravit ) @@ -3209,7 +3216,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & deltaz(i) = state_loc%zi(i,pverp-1) - state_loc%zi(i,pverp) - ! Set the surface pressure + ! Set the surface pressure p_sfc(i) = state_loc%pint(i,pverp) ! Set the elevation of the surface @@ -3221,7 +3228,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & !$acc parallel loop gang vector collapse(2) default(present) do k = 1, nzm_clubb do i = 1, ncol - k_cam = top_lev - 1 + k + k_cam = top_lev - 1 + k zi_g(i,k) = state_loc%zi(i,k_cam) - state_loc%zi(i,pverp) end do end do @@ -3298,7 +3305,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & call t_stopf('clubb_tend_cam:acc_copyin') call t_startf('clubb_tend_cam:acc_region') !----------------------------------------- END CLUBB grid initialization ----------------------------------------- - + #ifdef SILHS ! Add forcings for SILHS covariance contributions rtp2_forcing = zt2zm_api( nzm_clubb, nzt_clubb, ncol, gr, rtp2_mc_zt_pbuf(1:ncol,:) ) @@ -3342,7 +3349,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & !$acc parallel loop gang vector default(present) do i = 1, ncol wpthlp_sfc(i) = cam_in%shf(i) / ( cpairv(i,pver,lchnk) * rho_ds_zm(i,nzm_clubb) ) & ! Sensible heat flux - * invrs_exner_zt(i,nzt_clubb) + * invrs_exner_zt(i,nzt_clubb) wprtp_sfc(i) = cam_in%cflx(i,1) / rho_ds_zm(i,nzm_clubb) ! Moisture flux end do @@ -3396,7 +3403,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & vpwp_sfc(1) = -vm(1,nzt_clubb)*ustar**2/ubar end if - + ! Implementation after Thomas Toniazzo (NorESM) and Colin Zarzycki (PSU) ! Other Surface fluxes provided by host model if( (cld_macmic_num_steps > 1) .and. clubb_l_intr_sfc_flux_smooth ) then @@ -3431,7 +3438,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & endif - ! We only need to copy pdf_params from pbuf if this is a restart, we're calling pdf_closure + ! We only need to copy pdf_params from pbuf if this is a restart, we're calling pdf_closure ! at the end of advance_clubb_core, and calling it twice for pdf_params_zm as well if ( is_first_restart_step() & .and. clubb_config_flags%l_call_pdf_closure_twice & @@ -3452,7 +3459,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & if ( edsclr_dim > 0 ) then - ! Copy the cam version of the tracers to the clubb version + ! Copy the cam version of the tracers to the clubb version ! NOTE: if clubb_l_do_expldiff_rtm_thlm=.true., then the last two ! tracers are thlm and rtm, which are added inside clubb icnt=0 @@ -3542,18 +3549,18 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & end if - + if ( clubb_l_ascending_grid ) then - ! CLUBB is to be run in ascending mode, which has the surface at k=1, which is + ! CLUBB is to be run in ascending mode, which has the surface at k=1, which is ! the opposite of the cam grid that the rest of clubb_intr uses, so ! we need to flip the fields (in the vertical dimensions) before calling advance_clubb_core ! ! NOTE: We do not neccesarily flip all arrays, only ones that are used within this - ! subroutine (advance_clubb_core). For example, only the pdf_params fields that + ! subroutine (advance_clubb_core). For example, only the pdf_params fields that ! are used within this subroutine (or used in a subroutine we call) need to - ! be flipped. - + ! be flipped. + call t_startf('clubb_tend_cam:ascending_grid_flip') thlm_forcing = thlm_forcing(:,nzt_clubb:1:-1) @@ -3633,18 +3640,18 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & end if if ( sclr_dim > 0 ) then - + sclrm_forcing = sclrm_forcing(:,nzt_clubb:1:-1,:) sclrm = sclrm(:,nzt_clubb:1:-1,:) sclrp3 = sclrp3(:,nzt_clubb:1:-1,:) - + sclrp2 = sclrp2(:,nzm_clubb:1:-1,:) sclrprtp = sclrprtp(:,nzm_clubb:1:-1,:) sclrpthlp = sclrpthlp(:,nzm_clubb:1:-1,:) wpsclrp = wpsclrp(:,nzm_clubb:1:-1,:) sclrpthvp = sclrpthvp(:,nzm_clubb:1:-1,:) end if - + ! These are flipped, ensuring these are stored in descending mode, regardless of clubb_l_ascending_grid. ! only because these are need to be stored for restarts if ( clubb_config_flags%l_call_pdf_closure_twice ) then @@ -3654,9 +3661,9 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & pdf_params_zm_chnk(lchnk)%varnce_w_2 = pdf_params_zm_chnk(lchnk)%varnce_w_2(:,nzm_clubb:1:-1) pdf_params_zm_chnk(lchnk)%mixt_frac = pdf_params_zm_chnk(lchnk)%mixt_frac (:,nzm_clubb:1:-1) end if - + ! These are flipped, ensuring these are stored in descending mode, regardless of clubb_l_ascending_grid. - ! only for pdfp_rtp2_output calc + ! only for pdfp_rtp2_output calc pdf_params_chnk(lchnk)%mixt_frac = pdf_params_chnk(lchnk)%mixt_frac (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%rt_1 = pdf_params_chnk(lchnk)%rt_1 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%rt_2 = pdf_params_chnk(lchnk)%rt_2 (:,nzt_clubb:1:-1) @@ -3673,7 +3680,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & pdf_params_chnk(lchnk)%thl_2 = pdf_params_chnk(lchnk)%thl_2 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%varnce_thl_1 = pdf_params_chnk(lchnk)%varnce_thl_1(:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%varnce_thl_2 = pdf_params_chnk(lchnk)%varnce_thl_2(:,nzt_clubb:1:-1) - + ! These are flipped for silhs, which uses a cam grid pdf_params_chnk(lchnk)%rc_1 = pdf_params_chnk(lchnk)%rc_1 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%rc_2 = pdf_params_chnk(lchnk)%rc_2 (:,nzt_clubb:1:-1) @@ -3693,7 +3700,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & pdf_params_chnk(lchnk)%corr_chi_eta_2 = pdf_params_chnk(lchnk)%corr_chi_eta_2 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%corr_w_chi_1 = pdf_params_chnk(lchnk)%corr_w_chi_1 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%corr_w_chi_2 = pdf_params_chnk(lchnk)%corr_w_chi_2 (:,nzt_clubb:1:-1) - + call cleanup_grid_api( gr ) @@ -3703,7 +3710,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & deltaz, zi_g(:,1), zi_g(:,nzm_clubb), & ! intent(in) zi_g(:,nzm_clubb:1:-1), zt_g(:,nzt_clubb:1:-1), & ! intent(in) gr, err_info ) ! intent(inout) - + call t_stopf('clubb_tend_cam:ascending_grid_flip') end if @@ -3780,14 +3787,14 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & !$acc wp2up2_pbuf, wp2vp2_pbuf, ice_supersat_frac_pbuf ) call t_stopf('clubb_tend_cam:advance_clubb_core_api') - + if ( clubb_l_ascending_grid ) then call t_startf('clubb_tend_cam:ascending_grid_flip') ! If running in ascending mode, we flip the arrays before calling advance_clubb_core - ! so we need to flip them back. This section should flip every array that was flipped + ! so we need to flip them back. This section should flip every array that was flipped ! before the advance_clubb_core call. thlm_forcing = thlm_forcing(:,nzt_clubb:1:-1) @@ -3880,7 +3887,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & end if if ( sclr_dim > 0 ) then - + sclrm_forcing = sclrm_forcing(:,nzt_clubb:1:-1,:) sclrm = sclrm(:,nzt_clubb:1:-1,:) sclrp3 = sclrp3(:,nzt_clubb:1:-1,:) @@ -3891,7 +3898,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & wpsclrp = wpsclrp(:,nzm_clubb:1:-1,:) sclrpthvp = sclrpthvp(:,nzm_clubb:1:-1,:) end if - + ! These are flipped, ensuring these are stored in descending mode, regardless of clubb_l_ascending_grid ! only because these are need to be stored for restarts if ( clubb_config_flags%l_call_pdf_closure_twice ) then @@ -3901,16 +3908,16 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & pdf_params_zm_chnk(lchnk)%varnce_w_2 = pdf_params_zm_chnk(lchnk)%varnce_w_2(:,nzm_clubb:1:-1) pdf_params_zm_chnk(lchnk)%mixt_frac = pdf_params_zm_chnk(lchnk)%mixt_frac (:,nzm_clubb:1:-1) end if - - ! These are flipped, ensuring these are stored in descending mode, regardless of clubb_l_ascending_grid - ! only for pdfp_rtp2_output calc + + ! These are flipped, ensuring these are stored in descending mode, regardless of clubb_l_ascending_grid + ! only for pdfp_rtp2_output calc pdf_params_chnk(lchnk)%mixt_frac = pdf_params_chnk(lchnk)%mixt_frac (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%rt_1 = pdf_params_chnk(lchnk)%rt_1 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%rt_2 = pdf_params_chnk(lchnk)%rt_2 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%varnce_rt_1 = pdf_params_chnk(lchnk)%varnce_rt_1(:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%varnce_rt_2 = pdf_params_chnk(lchnk)%varnce_rt_2(:,nzt_clubb:1:-1) - ! These are flipped, ensuring these are stored in descending mode, regardless of clubb_l_ascending_grid + ! These are flipped, ensuring these are stored in descending mode, regardless of clubb_l_ascending_grid ! only for update_xp2_mc_api call pdf_params_chnk(lchnk)%w_1 = pdf_params_chnk(lchnk)%w_1 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%w_2 = pdf_params_chnk(lchnk)%w_2 (:,nzt_clubb:1:-1) @@ -3920,7 +3927,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & pdf_params_chnk(lchnk)%thl_2 = pdf_params_chnk(lchnk)%thl_2 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%varnce_thl_1 = pdf_params_chnk(lchnk)%varnce_thl_1(:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%varnce_thl_2 = pdf_params_chnk(lchnk)%varnce_thl_2(:,nzt_clubb:1:-1) - + ! These are flipped for silhs, which uses a cam grid pdf_params_chnk(lchnk)%rc_1 = pdf_params_chnk(lchnk)%rc_1 (:,nzt_clubb:1:-1) pdf_params_chnk(lchnk)%rc_2 = pdf_params_chnk(lchnk)%rc_2 (:,nzt_clubb:1:-1) @@ -3998,7 +4005,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & call t_startf('clubb_tend_cam:do_cldcool') thlp2_rad(:,:) = 0._r8 - + do k = 1, nzt_clubb do i = 1, ncol k_cam = top_lev - 1 + k @@ -4051,7 +4058,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & end do end do - ! pdf_params_zm_chnk is already persistent across calls, but we + ! pdf_params_zm_chnk is already persistent across calls, but we ! save a pbuf version for restarts if ( clubb_config_flags%l_call_pdf_closure_twice ) then !$acc parallel loop gang vector collapse(2) default(present) @@ -4207,7 +4214,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & do k = top_lev, pverp do i = 1, ncol k_clubb = k + 1 - top_lev - tke_pbuf(i,k) = 0.5_r8 * ( up2_pbuf(i,k_clubb) + vp2_pbuf(i,k_clubb) + wp2_pbuf(i,k_clubb) ) + tke_pbuf(i,k) = 0.5_r8 * ( up2_pbuf(i,k_clubb) + vp2_pbuf(i,k_clubb) + wp2_pbuf(i,k_clubb) ) enddo enddo @@ -4337,7 +4344,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & ptend_loc%q(i,k,ixind) = 0._r8 end do end do - + ! Copy CLUBB's edsclr values do k = top_lev, pver do i = 1, ncol @@ -4574,7 +4581,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & ! "cloud_frac"), compute the convective cloud fraction. This follows the formulation ! found in macrophysics code. Assumes that convective cloud is all nonstratiform cloud ! from CLUBB plus the deep convective cloud fraction - ! NOTE: concld_pbuf used to be calculated in the commented-out version below, but since we + ! NOTE: concld_pbuf used to be calculated in the commented-out version below, but since we ! set alst_pbuf=cloud_frac_pbuf, this simplifies to only using deepcu_pbuf. ! This is potentially a bug, but there's not really a "right" way to combine the different ! cloud factions, so it has been left to only use deepcu_pbuf for now @@ -4610,6 +4617,12 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & aist_pbuf(:,:top_lev-1) = 0._r8 qsatfac_pbuf(:, :) = 0._r8 ! Zero out entire profile in case qsatfac is left undefined in aist_vector below + rhmini_default(:) = rhmini_const + rhmaxi_default(:) = rhmaxi_const + rhminl_arr(:) = rhminl_const + rhminl_adj_land_arr(:) = rhminl_adj_land_const + rhminh_arr(:) = rhminh_const + do k = top_lev, pver ! For Type II PSC and for thin cirrus, the clouds can be thin, but @@ -4628,13 +4641,30 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & rhmaxi = rhmaxi_const end where + !REMOVECAM: this is no longer needed when CAM is retired and pcols no longer exists + aist_pbuf(:,k) = 0._r8 + !REMOVECAM_END if ( trim(subcol_scheme) == 'SILHS' ) then - call aist_vector(state_loc%q(:,k,ixq),state_loc%t(:,k),state_loc%pmid(:,k),state_loc%q(:,k,ixcldice), & - state_loc%q(:,k,ixnumice), cam_in%landfrac(:),cam_in%snowhland(:),aist_pbuf(:,k),ncol ) + call aist_vector(state_loc%q(:ncol,k,ixq), state_loc%t(:ncol,k), & + state_loc%pmid(:ncol,k), state_loc%q(:ncol,k,ixcldice), & + state_loc%q(:ncol,k,ixnumice), cam_in%landfrac(:ncol), & + cam_in%snowhland(:ncol), aist_pbuf(:ncol,k), ncol, & + rhmaxi_in=rhmaxi_default(:ncol), & + rhmini_in=rhmini_default(:ncol), & + rhminl_in=rhminl_arr(:ncol), & + rhminl_adj_land_in=rhminl_adj_land_arr(:ncol), & + rhminh_in=rhminh_arr(:ncol)) else - call aist_vector(state_loc%q(:,k,ixq),state_loc%t(:,k),state_loc%pmid(:,k),state_loc%q(:,k,ixcldice), & - state_loc%q(:,k,ixnumice), cam_in%landfrac(:),cam_in%snowhland(:),aist_pbuf(:,k),ncol,& - qsatfac_out=qsatfac_pbuf(:,k), rhmini_in=rhmini, rhmaxi_in=rhmaxi) + call aist_vector(state_loc%q(:ncol,k,ixq), state_loc%t(:ncol,k), & + state_loc%pmid(:ncol,k), state_loc%q(:ncol,k,ixcldice), & + state_loc%q(:ncol,k,ixnumice), cam_in%landfrac(:ncol), & + cam_in%snowhland(:ncol), aist_pbuf(:ncol,k), ncol, & + rhmaxi_in=rhmaxi(:ncol), & + rhmini_in=rhmini(:ncol), & + rhminl_in=rhminl_arr(:ncol), & + rhminl_adj_land_in=rhminl_adj_land_arr(:ncol), & + rhminh_in=rhminh_arr(:ncol), & + qsatfac_out=qsatfac_pbuf(:ncol,k)) endif enddo diff --git a/src/physics/carma/models/cirrus/carma_cloudfraction.F90 b/src/physics/carma/models/cirrus/carma_cloudfraction.F90 index 0ec202041f..5e3614d2d6 100644 --- a/src/physics/carma/models/cirrus/carma_cloudfraction.F90 +++ b/src/physics/carma/models/cirrus/carma_cloudfraction.F90 @@ -30,7 +30,8 @@ subroutine CARMA_CloudFraction(carma, cstate, cam_in, state, icol, cldfrc, rhcri use camsrfexch, only : cam_in_t use ppgrid, only : pcols, pver, pverp - use cldfrc2m, only : astG_RHU_single, astG_PDF_single, aist_single, CAMstfrac + use compute_cloud_fraction_two_moment, only : astG_RHU_single, astG_PDF_single, aist_single, CAMstfrac + use cldfrc2m, only : rhmini_const, rhmaxi_const, rhminl_const, rhminl_adj_land_const, rhminh_const type(carma_type) :: carma !! the carma object type(carmastate_type) :: cstate !! the carma state object @@ -113,15 +114,21 @@ subroutine CARMA_CloudFraction(carma, cstate, cam_in, state, icol, cldfrc, rhcri ! it starts to be used, then a general routine astG_single should be written. if (CAMstfrac) then call astG_RHU_single(ssl + 1._f, state%pmid(icol, iz), state%q(icol, iz, 1), & - cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, rhcrit(iz)) + cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, & + rhminl=rhminl_const, rhminl_adj_land=rhminl_adj_land_const, & + rhminh=rhminh_const, orhmin=rhcrit(iz)) else call astG_PDF_single(ssl + 1._f, state%pmid(icol, iz), state%q(icol, iz, 1), & - cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, rhcrit(iz)) + cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, & + rhminl=rhminl_const, rhminl_adj_land=rhminl_adj_land_const, & + rhminh=rhminh_const, orhmin=rhcrit(iz)) end if ! Now get the ice cloud fraction. call aist_single(state%q(icol, iz, 1), state%t(icol, iz), state%pmid(icol, iz), & - qi(iz), cam_in%landfrac(icol), cam_in%snowhland(icol), icecldf(iz)) + qi(iz), cam_in%landfrac(icol), cam_in%snowhland(icol), icecldf(iz), & + rhmaxi=rhmaxi_const, rhmini=rhmini_const, rhminl=rhminl_const, & + rhminl_adj_land=rhminl_adj_land_const, rhminh=rhminh_const) end do ! Calculate an overall cloud fraction. This may vary depending upon the model, diff --git a/src/physics/carma/models/cirrus_dust/carma_cloudfraction.F90 b/src/physics/carma/models/cirrus_dust/carma_cloudfraction.F90 index 0ec202041f..5e3614d2d6 100644 --- a/src/physics/carma/models/cirrus_dust/carma_cloudfraction.F90 +++ b/src/physics/carma/models/cirrus_dust/carma_cloudfraction.F90 @@ -30,7 +30,8 @@ subroutine CARMA_CloudFraction(carma, cstate, cam_in, state, icol, cldfrc, rhcri use camsrfexch, only : cam_in_t use ppgrid, only : pcols, pver, pverp - use cldfrc2m, only : astG_RHU_single, astG_PDF_single, aist_single, CAMstfrac + use compute_cloud_fraction_two_moment, only : astG_RHU_single, astG_PDF_single, aist_single, CAMstfrac + use cldfrc2m, only : rhmini_const, rhmaxi_const, rhminl_const, rhminl_adj_land_const, rhminh_const type(carma_type) :: carma !! the carma object type(carmastate_type) :: cstate !! the carma state object @@ -113,15 +114,21 @@ subroutine CARMA_CloudFraction(carma, cstate, cam_in, state, icol, cldfrc, rhcri ! it starts to be used, then a general routine astG_single should be written. if (CAMstfrac) then call astG_RHU_single(ssl + 1._f, state%pmid(icol, iz), state%q(icol, iz, 1), & - cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, rhcrit(iz)) + cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, & + rhminl=rhminl_const, rhminl_adj_land=rhminl_adj_land_const, & + rhminh=rhminh_const, orhmin=rhcrit(iz)) else call astG_PDF_single(ssl + 1._f, state%pmid(icol, iz), state%q(icol, iz, 1), & - cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, rhcrit(iz)) + cam_in%landfrac(icol), cam_in%snowhland(icol), liqcldf(iz), Ga, & + rhminl=rhminl_const, rhminl_adj_land=rhminl_adj_land_const, & + rhminh=rhminh_const, orhmin=rhcrit(iz)) end if ! Now get the ice cloud fraction. call aist_single(state%q(icol, iz, 1), state%t(icol, iz), state%pmid(icol, iz), & - qi(iz), cam_in%landfrac(icol), cam_in%snowhland(icol), icecldf(iz)) + qi(iz), cam_in%landfrac(icol), cam_in%snowhland(icol), icecldf(iz), & + rhmaxi=rhmaxi_const, rhmini=rhmini_const, rhminl=rhminl_const, & + rhminl_adj_land=rhminl_adj_land_const, rhminh=rhminh_const) end do ! Calculate an overall cloud fraction. This may vary depending upon the model,