diff --git a/bld/build-namelist b/bld/build-namelist index efd42d4949..90bd503820 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -3380,6 +3380,9 @@ if ($clubb_sgs =~ /$TRUE/io) { add_default($nl, 'clubb_mf_L0'); add_default($nl, 'clubb_mf_ent0'); add_default($nl, 'clubb_mf_nup'); + + #Turn on HB scheme where CLUBB not active + add_default($nl, 'do_hb_above_clubb'); } # Force exit if running cam_dev and CLUBB is off @@ -3960,6 +3963,7 @@ if ($dyn =~ /se/) { se_kmin_jet se_kmax_jet se_molecular_diff + se_pgf_formulation ); my %opts; diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index 75cc4f4677..751ab84c6d 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -266,7 +266,7 @@ atm/waccm/ic/wa3_ne5np4_1950_spinup.cam2.i.1960-01-01-00000_c150810.nc atm/waccm/ic/FW2000_ne30_L70_01-01-0001_c200602.nc atm/waccm/ic/FWsc2000climo_ne30pg3_L70_0002-01-01_c221103.nc -atm/waccm/ic/FW2000_ne30pg3_L70_01-01-0001_c200602.nc +atm/waccm/ic/FW2000.ne30pg3_ne30pg3_nlev70_c230906.nc atm/waccm/ic/FWsc2000_ne30pg3_L110_01-01-0001_c200521.nc atm/waccm/ic/FW2000_CONUS_30x8_L70_01-01-0001_c200602.nc @@ -2117,7 +2117,9 @@ 1 2.0 60.0 - +.false. +.true. +.true. .true. 0.2 @@ -3011,6 +3013,8 @@ + 1 +2 2 .true. diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index e3f4333864..dfe47f84e5 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -5169,6 +5169,13 @@ Bretherton; 'HBR' for Rasch modified version of 'HB'. Default: set by build-namelist + +Logical: If True activate Holtslag and Boville vertical diffusion scheme where CLUBB is not active + (note that CLUBB top is dynamic in each column) +Default: Set by build-namelist. + + + + + 1: Exner version of pressure gradient force (PGF) + see Appendix A in https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2022MS003192 + + 2: Traditional pressure gradient formulation (grad p) + + Default: Set by build-namelist. + + diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam index f695957589..41eded30db 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam @@ -3,6 +3,5 @@ ndens=1,1,1,1,1,1,1,1,1 nhtfrq=9,9,9,9,9,9,9,9,9 inithist='ENDOFRUN' -se_nsplit = 30 state_debug_checks = .true. diff --git a/doc/ChangeLog b/doc/ChangeLog index 6ee1ed2f62..af188b509a 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,198 @@ =============================================================== +Tag name: cam6_3_127 +Originator(s): pel, jet +Date: Sept 12, 2023 +One-line Summary: Option to turn on HB where CLUBB is not active +Github PR URL: https://github.com/ESCOMP/CAM/pull/849 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): + - Free atmosphere Richardson number based mixing (using the free atmosphere + part of the Holtslag-Boville PBL scheme) for those layers where CLUBB is + not active (turned on in cam6 or cam_dev). Also helps stabilize less diffusive + dynamical cores. The issue for this PR is #846 and it will also address issue #772 + +Describe any changes made to build system: NA + +Describe any changes made to the namelist: + - Introduced namelist variable do_hb_above_clubb. Default False unless using cam6 or cam_dev then + Holtslag-Boville PBL scheme will be used above CLUBB top) + - new namelist variable se_pgf_formulation to allow use of older pgf calculation to help with waccm regression errors + Can have an interger value of 1 or 2. + 1: Exner version of pressure gradient force (PGF) + 2: Traditional pressure gradient formulation (grad p) + +List any changes to the defaults for the boundary datasets: N/A + - New Boundary data for ERP and SMS WACCM runs - these runs require a spun up IC. + inputdata/atm/waccm/ic/FW2000.ne30pg3_ne30pg3_nlev70_c230906.nc + +Describe any substantial timing or memory changes: N/A + +Code reviewed by:fvitt, pel, eaton, cacraig, jesse + +List all files eliminated: N/A + +List all files added and what they do: N/A + +List all existing files that have been modified, and describe the changes: + +M bld/build-namelist + - Turn on HB scheme where CLUBB not active + +M bld/namelist_files/namelist_defaults_cam.xml +M bld/namelist_files/namelist_definition.xml + - add new namelist parameter for do_hb_above_clubb + - add new ic file to fix existing waccm regression test failures + +M cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam + - use default se_nsplit for waccm runs using new ic file + +M src/dynamics/mpas/dp_coupling.F90 + - bug fix, looping should only be over wet species + +M src/dynamics/se/dp_coupling.F90 + - bug fix, looping should only be over wet species + +M src/dynamics/se/dycore/global_norms_mod.F90 +M src/physics/cam/clubb_intr.F90 + - set pbuf field so that HB scheme is only applied above CLUBB top + +M src/physics/cam/hb_diff.F90 +M src/physics/cam/vertical_diffusion.F90 + - add call to run HB scheme where CLUBB not active + +M src/physics/cam/geopotential.F90 + - bug fix for converting wet to dry mixing ratio + +M src/physics/cam/phys_control.F90 + - do_hb_above_clubb defined here and provided via getopts. + +M src/dynamics/se/dycore/prim_advance_mod.F90 +M src/dynamics/se/dyn_comp.F90 +M src/dynamics/se/dycore/control_mod.F90 + - remove duplicate timing call + - allow use of older pgf formulation to help with waccm regression errors + +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. + +cheyenne/intel/aux_cam: all BFB except: + + ERP_Ln9_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq9s (Overall: FAIL) details: + - pre-existing failure + + ERC_D_Ln9_Vnuopc.ne16_ne16_mg17.FADIAB.cheyenne_intel.cam-terminator (Overall: NLFAIL) details: + ERC_D_Ln9_Vnuopc.ne16_ne16_mg17.QPC5HIST.cheyenne_intel.cam-outfrq3s_usecase (Overall: NLFAIL) details: + ERS_Ln9_Vnuopc.ne0TESTONLYne5x4_ne0TESTONLYne5x4_mg37.FADIAB.cheyenne_intel.cam-outfrq3s_refined (Overall: NLFAIL) details: + - expected NLFAIL due to new se_pgf_formulation namelist variable. + + SMS_Ld5_Vnuopc.f09_f09_mg17.PC6.cheyenne_intel.cam-cam6_port_f09 (Overall: NLFAIL) details: + - expected NLFAIL due to new do_hb_above_clubb namelist variable. + + ERC_D_Ln9_P144x1_Vnuopc.ne16pg3_ne16pg3_mg17.QPC6HIST.cheyenne_intel.cam-outfrq3s_ttrac_usecase (Overall: DIFF) details: + ERC_D_Ln9_Vnuopc.f19_f19_mg17.QPC6.cheyenne_intel.cam-outfrq3s_cosp (Overall: DIFF) details: + ERP_D_Ln9.ne30pg3_ne30pg3_mg17.FLTHIST.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_D_Ln9_Vnuopc.f09_f09_mg17.QSC6.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_D_Ln9_Vnuopc.f19_f19_mg17.QPC6.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_D_Ln9_Vnuopc.ne30pg3_ne30pg3_mg17.F2000dev.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_Ld3_Vnuopc.f09_f09_mg17.FWHIST.cheyenne_intel.cam-reduced_hist1d (Overall: DIFF) details: + ERP_Ln9_P24x3_Vnuopc.f45_f45_mg37.QPWmaC6.cheyenne_intel.cam-outfrq9s_mee_fluxes (Overall: DIFF) details: + ERP_Ln9_Vnuopc.C96_C96_mg17.F2000climo.cheyenne_intel.cam-outfrq9s_mg3 (Overall: DIFF) details: + ERP_Ln9_Vnuopc.f09_f09_mg17.F1850.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_Ln9_Vnuopc.f09_f09_mg17.F2000climo.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_Ln9_Vnuopc.f09_f09_mg17.F2000dev.cheyenne_intel.cam-outfrq9s_mg3 (Overall: DIFF) details: + ERP_Ln9_Vnuopc.f09_f09_mg17.F2010climo.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_Ln9_Vnuopc.f09_f09_mg17.FHIST_BDRD.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_Ln9_Vnuopc.f19_f19_mg17.FWsc1850.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_Ln9_Vnuopc.ne30_ne30_mg17.FCnudged.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERP_Ln9_Vnuopc.ne30pg3_ne30pg3_mg17.FW2000climo.cheyenne_intel.cam-outfrq9s_wcm_ne30 (Overall: DIFF) details: + ERS_Ld3_Vnuopc.f10_f10_mg37.F1850.cheyenne_intel.cam-outfrq1d_14dec_ghg_cam_dev (Overall: DIFF) details: + ERS_Ln9_P288x1_Vnuopc.mpasa120_mpasa120.F2000climo.cheyenne_intel.cam-outfrq9s_mpasa120 (Overall: DIFF) details: + ERS_Ln9_P36x1_Vnuopc.mpasa480_mpasa480.F2000climo.cheyenne_intel.cam-outfrq9s_mpasa480 (Overall: DIFF) details: + ERS_Ln9_Vnuopc.f09_f09_mg17.FX2000.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + ERS_Ln9_Vnuopc.f19_f19_mg17.FXSD.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9.ne30pg3_ne30pg3_mg17.FMTHIST.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f09_f09_mg17.FCts2nudged.cheyenne_intel.cam-outfrq9s_leapday (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f09_f09_mg17.FCvbsxHIST.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f09_f09_mg17.FSD.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f19_f19_mg17.FWma2000climo.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f19_f19_mg17.FWma2000climo.cheyenne_intel.cam-outfrq9s_waccm_ma_mam4 (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f19_f19_mg17.FXHIST.cheyenne_intel.cam-outfrq9s_amie (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f19_f19_mg17.QPC2000climo.cheyenne_intel.cam-outfrq3s_usecase (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.ne16_ne16_mg17.FX2000.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.ne16_ne16_mg17.QPX2000.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc_P720x1.ne0ARCTICne30x4_ne0ARCTICne30x4_mt12.FHIST.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc_P720x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc_P720x1.ne30pg3_ne30pg3_mg17.FCLTHIST.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.T42_T42.FSCAM.cheyenne_intel.cam-outfrq9s (Overall: DIFF) details: + SMS_Ld1_Vnuopc.f09_f09_mg17.FW2000climo.cheyenne_intel.cam-outfrq1d (Overall: DIFF) details: + SMS_Ld1_Vnuopc.f19_f19.F2000dev.cheyenne_intel.cam-outfrq1d (Overall: DIFF) details: + SMS_Ld1_Vnuopc.ne30pg3_ne30pg3_mg17.FC2010climo.cheyenne_intel.cam-outfrq1d (Overall: DIFF) details: + SMS_Lh12_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq3h (Overall: DIFF) details: + SMS_Lm13_Vnuopc.f10_f10_mg37.F2000climo.cheyenne_intel.cam-outfrq1m (Overall: DIFF) details: + SMS_Ln9_Vnuopc.f09_f09_mg17.F2010climo.cheyenne_intel.cam-nudging (Overall: DIFF) details: + SMS_Ln9_Vnuopc.f09_f09_mg17.FW1850.cheyenne_intel.cam-reduced_hist3s (Overall: DIFF) details: + SMS_Ln9_Vnuopc.f19_f19.F2000climo.cheyenne_intel.cam-silhs (Overall: DIFF) details: + SMS_Ln9_Vnuopc.f19_f19_mg17.FHIST.cheyenne_intel.cam-outfrq9s_nochem (Overall: DIFF) details: + - expected diff for cam_dev and cam6 run using do_hb_above_clubb=.true. + + +izumi/nag/aux_cam: all BFB except: + DAE_Vnuopc.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) + - pre-existing failure + + ERC_D_Ln9_Vnuopc.ne16_ne16_mg17.QPC4.izumi_nag.cam-outfrq3s_usecase (Overall: NLFAIL) details: + ERC_D_Ln9_Vnuopc.ne5_ne5_mg37.QPC5.izumi_nag.cam-outfrq3s_ttrac (Overall: NLFAIL) details: + ERC_D_Ln9_Vnuopc.ne16pg3_ne16pg3_mg17.QPC4.izumi_nag.cam-outfrq3s_usecase (Overall: NLFAIL) details: + ERI_D_Ln18_Vnuopc.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq3s_bwic (Overall: NLFAIL) details: + ERI_D_Ln18_Vnuopc.ne5pg3_ne5pg3_mg37.FADIAB.izumi_nag.cam-outfrq3s_bwic (Overall: NLFAIL) details: + ERS_Ln27_Vnuopc.ne5pg3_ne5pg3_mg37.FKESSLER.izumi_nag.cam-outfrq9s (Overall: NLFAIL) details: + ERS_Ln9_Vnuopc.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq9s (Overall: NLFAIL) details: + PEM_D_Ln9_Vnuopc.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq3s (Overall: NLFAIL) details: + PLB_D_Ln9_Vnuopc.ne5_ne5_mg37.QPC5.izumi_nag.cam-ttrac_loadbal0 (Overall: NLFAIL) details: + PLB_D_Ln9_Vnuopc.ne5_ne5_mg37.QPC5.izumi_nag.cam-ttrac_loadbal1 (Overall: NLFAIL) details: + PLB_D_Ln9_Vnuopc.ne5_ne5_mg37.QPC5.izumi_nag.cam-ttrac_loadbal3 (Overall: NLFAIL) details: + SMS_D_Ln6_Vnuopc.ne5_ne5_mg37.QPWmaC4.izumi_nag.cam-outfrq3s_physgrid_tem (Overall: NLFAIL) details: + SMS_D_Ln9_P1x1_Vnuopc.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq3s (Overall: NLFAIL) details: + - expected NLFAIL due to new se_pgf_formulation namelist variable. + + ERC_D_Ln9_Vnuopc.f10_f10_mg37.QPC6.izumi_nag.cam-outfrq3s_am (Overall: DIFF) details: + ERC_D_Ln9_Vnuopc.f10_f10_mg37.QPC6.izumi_nag.cam-outfrq3s_convmic (Overall: DIFF) details: + ERC_D_Ln9_Vnuopc.f10_f10_mg37.QPC6.izumi_nag.cam-outfrq3s_cospsathist (Overall: DIFF) details: + ERC_D_Ln9_Vnuopc.f10_f10_mg37.QPC6.izumi_nag.cam-outfrq3s (Overall: DIFF) details: + ERC_D_Ln9_Vnuopc.f10_f10_mg37.QPWmaC6.izumi_nag.cam-outfrq3s (Overall: DIFF) details: + ERI_D_Ln18_Vnuopc.f19_f19_mg17.QPC6.izumi_nag.cam-ghgrmp_e8 (Overall: DIFF) details: + ERP_Ln9_Vnuopc.ne5pg3_ne5pg3_mg37.QPC6.izumi_nag.cam-outfrq9s_clubbmf (Overall: DIFF) details: + SMS_D_Ln3_Vnuopc.ne5pg3_ne5pg3_mg37.QPX2000.izumi_nag.cam-outfrq3s (Overall: DIFF) details: + SMS_D_Ln9_Vnuopc.f10_f10_mg37.QPC6.izumi_nag.cam-outfrq3s_ba (Overall: DIFF) details: + SMS_P48x1_D_Ln3_Vnuopc.f09_f09_mg17.QPC6HIST.izumi_nag.cam-outfrq3s_co2cycle_usecase (Overall: DIFF) details: + - expected diff for cam_dev/cam6 runs using do_hb_above_clubb=.true. + +izumi/gnu/aux_cam: all BFB except: + + ERC_D_Ln9_Vnuopc.ne5_ne5_mg37.QPC4.izumi_gnu.cam-outfrq3s_nudging_ne5_L26 (Overall: NLFAIL) details: + ERC_D_Ln9_Vnuopc.ne5_ne5_mg37.QPC5.izumi_gnu.cam-outfrq3s_ba (Overall: NLFAIL) details: + ERC_D_Ln9_Vnuopc.ne5pg2_ne5pg2_mg37.FADIAB.izumi_gnu.cam-outfrq3s (Overall: NLFAIL) details: + ERC_D_Ln9_Vnuopc.ne5pg4_ne5pg4_mg37.FADIAB.izumi_gnu.cam-outfrq3s (Overall: NLFAIL) details: + ERP_Ln9_Vnuopc.ne5_ne5_mg37.FHS94.izumi_gnu.cam-outfrq9s (Overall: NLFAIL) details: + ERP_Ln9_Vnuopc.ne5_ne5_mg37.QPC5.izumi_gnu.cam-outfrq9s (Overall: NLFAIL) details: + PEM_D_Ln9_Vnuopc.ne5pg3_ne5pg3_mg37.FADIAB.izumi_gnu.cam-outfrq3s (Overall: NLFAIL) details: + PLB_D_Ln9_Vnuopc.ne5pg3_ne5pg3_mg37.QPC5.izumi_gnu.cam-ttrac_loadbal0 (Overall: NLFAIL) details: + PLB_D_Ln9_Vnuopc.ne5pg3_ne5pg3_mg37.QPC5.izumi_gnu.cam-ttrac_loadbal1 (Overall: NLFAIL) details: + PLB_D_Ln9_Vnuopc.ne5pg3_ne5pg3_mg37.QPC5.izumi_gnu.cam-ttrac_loadbal3 (Overall: NLFAIL) details: + SMS_D_Ln9_Vnuopc.ne5pg3_ne5pg3_mg37.QPC5.izumi_gnu.cam-outfrq3s_ttrac (Overall: NLFAIL) details: + - expected NLFAIL due to new se_pgf_formulation namelist variable. + + SMS_D_Ln9.f10_f10_mg37.2000_CAM%DEV%GHGMAM4_CLM50%SP_CICE%PRES_DOCN%DOM_MOSART_SGLC_SWAV_SESP.izumi_gnu.cam-outfrq9s (Overall: DIFF) + - expected diff for cam_dev/cam6 runs using do_hb_above_clubb=.true. + +Summarize any changes to answers: climate changing for cam6/cam-dev runs using CLUBB + +=============================================================== +=============================================================== + Tag name: cam6_3_126 Originator(s): gdicker Date: 6 Sep 2023 diff --git a/src/dynamics/mpas/dp_coupling.F90 b/src/dynamics/mpas/dp_coupling.F90 index 792a7d54b0..9366c6fac3 100644 --- a/src/dynamics/mpas/dp_coupling.F90 +++ b/src/dynamics/mpas/dp_coupling.F90 @@ -10,7 +10,7 @@ module dp_coupling use constituents, only: pcnst, cnst_type use physconst, only: gravit, cappa, rh2o, zvir use air_composition,only: cpairv, rairv - +use air_composition,only: dry_air_species_num use spmd_dyn, only: local_dp_map, block_buf_nrecs, chunk_buf_nrecs use spmd_utils, only: mpicom, iam, masterproc @@ -243,7 +243,7 @@ subroutine p_d_coupling(phys_state, phys_tend, dyn_in) allocate( t_tend(pver,nCellsSolve), stat=ierr) if( ierr /= 0 ) call endrun(subname//':failed to allocate t_tend array') - allocate( q_tend(thermodynamic_active_species_num,pver,nCellsSolve), stat=ierr) + allocate( q_tend(thermodynamic_active_species_num-dry_air_species_num,pver,nCellsSolve), stat=ierr) if( ierr /= 0 ) call endrun(subname//':failed to allocate q_tend array') nullify(tend_physics) @@ -285,7 +285,7 @@ subroutine p_d_coupling(phys_state, phys_tend, dyn_in) ! ! compute tendencies for thermodynamic active species ! - do m=1,thermodynamic_active_species_num + do m=dry_air_species_num + 1,thermodynamic_active_species_num idx_phys = thermodynamic_active_species_idx(m) idx_dycore = thermodynamic_active_species_idx_dycore(m) if (idx_dycore==index_qv) index_qv_phys = m @@ -388,7 +388,7 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d) ! To be consistent with total energy formula in physic's check_energy module only ! include water vapor in moist pdel. factor(:ncol,k) = 1.0_r8 - do m_cnst=1,thermodynamic_active_species_num + do m_cnst=dry_air_species_num + 1,thermodynamic_active_species_num m = thermodynamic_active_species_idx(m_cnst) ! at this point all q's are dry factor(:ncol,k) = factor(:ncol,k)+phys_state(lchnk)%q(:ncol,k,m) @@ -600,7 +600,7 @@ subroutine derived_tend(nCellsSolve, nCells, t_tend, u_tend, v_tend, q_tend, dyn ! ! Rnew and Rold are only needed for diagnostics purposes ! - do m=1,thermodynamic_active_species_num + do m=dry_air_species_num+1,thermodynamic_active_species_num idx_thermo(m) = m idx_dycore = thermodynamic_active_species_idx_dycore(m) do iCell = 1, nCellsSolve @@ -612,7 +612,7 @@ subroutine derived_tend(nCellsSolve, nCells, t_tend, u_tend, v_tend, q_tend, dyn call get_R(qktmp,idx_thermo,Rnew) Rnew = Rnew*cv/Rgas - do m=1,thermodynamic_active_species_num + do m=dry_air_species_num+1,thermodynamic_active_species_num idx_dycore = thermodynamic_active_species_idx_dycore(m) do iCell = 1, nCellsSolve do k = 1, pver @@ -673,7 +673,7 @@ subroutine derived_tend(nCellsSolve, nCells, t_tend, u_tend, v_tend, q_tend, dyn ! ! temporarily save thermodynamic active species (n+1) ! - do m=1,thermodynamic_active_species_num + do m=dry_air_species_num+1,thermodynamic_active_species_num idx_dycore = thermodynamic_active_species_idx_dycore(m) qk(m,:,: ) = tracers(idx_dycore,:,1:nCellsSolve) tracers(idx_dycore,:,1:nCellsSolve)= qk(m,:,: )-dtime*q_tend(m,:,1:nCellsSolve) @@ -685,7 +685,7 @@ subroutine derived_tend(nCellsSolve, nCells, t_tend, u_tend, v_tend, q_tend, dyn ux(:,1:nCellsSolve)+dtime*u_tend(:,1:nCellsSolve)/rho_zz(:,1:nCellsSolve), & uy(:,1:nCellsSolve)+dtime*v_tend(:,1:nCellsSolve)/rho_zz(:,1:nCellsSolve),'dAP') ! revert - do m=1,thermodynamic_active_species_num + do m=dry_air_species_num+1,thermodynamic_active_species_num idx_dycore = thermodynamic_active_species_idx_dycore(m) tracers(idx_dycore,:,1:nCellsSolve)= qk(m,:,: ) end do @@ -762,7 +762,7 @@ subroutine hydrostatic_pressure(nCells, nVertLevels, qsize, index_qv, zz, zgrid, do k = nVertLevels, 1, -1 rhodryk = zz(k,iCell)* rho_zz(k,iCell) !full CAM physics density rhok = 1.0_r8 - do idx=1,thermodynamic_active_species_num + do idx=dry_air_species_num+1,thermodynamic_active_species_num rhok = rhok+q(thermodynamic_active_species_idx_dycore(idx),k,iCell) end do rhok = rhok*rhodryk @@ -772,7 +772,7 @@ subroutine hydrostatic_pressure(nCells, nVertLevels, qsize, index_qv, zz, zgrid, k = nVertLevels sum_water = 1.0_r8 - do idx=1,thermodynamic_active_species_num + do idx=dry_air_species_num+1,thermodynamic_active_species_num sum_water = sum_water+q(thermodynamic_active_species_idx_dycore(idx),k,iCell) end do rhok = sum_water*zz(k,iCell) * rho_zz(k,iCell) @@ -790,7 +790,7 @@ subroutine hydrostatic_pressure(nCells, nVertLevels, qsize, index_qv, zz, zgrid, ! compute hydrostatic dry interface pressure so that (pintdry(k+1)-pintdry(k))/g is pseudo density ! sum_water = 1.0_r8 - do idx=1,thermodynamic_active_species_num + do idx=dry_air_species_num+1,thermodynamic_active_species_num sum_water = sum_water+q(thermodynamic_active_species_idx_dycore(idx),k,iCell) end do thetavk = theta_m(k,iCell)/sum_water!convert modified theta to virtual theta @@ -881,7 +881,7 @@ subroutine tot_energy_dyn(nCells, nVertLevels, qsize, index_qv, zz, zgrid, rho_z u(iCell,k) = ux(k,iCell) v(iCell,k) = uy(k,iCell) phis(iCell) = zgrid(1,iCell)*gravit - do idx=1,thermodynamic_active_species_num + do idx=dry_air_species_num+1,thermodynamic_active_species_num idx_tmp = thermodynamic_active_species_idx_dycore(idx) tracers(iCell,k,idx_tmp) = q(idx_tmp,k,iCell) end do diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 7dae784315..056eb45f93 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -605,7 +605,7 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) ! wet pressure variables (should be removed from physics!) factor_array(:,:) = 1.0_r8 - do m_cnst=1,thermodynamic_active_species_num + do m_cnst=dry_air_species_num + 1,thermodynamic_active_species_num m = thermodynamic_active_species_idx(m_cnst) do k=1,nlev do i=1,ncol diff --git a/src/dynamics/se/dycore/control_mod.F90 b/src/dynamics/se/dycore/control_mod.F90 index 4c1127c45b..053f478c6a 100644 --- a/src/dynamics/se/dycore/control_mod.F90 +++ b/src/dynamics/se/dycore/control_mod.F90 @@ -124,5 +124,7 @@ module control_mod real(r8), public :: molecular_diff = -1.0_r8 integer, public :: vert_remap_uvTq_alg, vert_remap_tracer_alg - + + + integer, public :: pgf_formulation = -1 !PGF formulation - see prim_advance_mod.F90 end module control_mod diff --git a/src/dynamics/se/dycore/global_norms_mod.F90 b/src/dynamics/se/dycore/global_norms_mod.F90 index 843fd88bc7..2b53ef95c7 100644 --- a/src/dynamics/se/dycore/global_norms_mod.F90 +++ b/src/dynamics/se/dycore/global_norms_mod.F90 @@ -681,8 +681,10 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& do k=1,nlev ! Vertical profile from FV dycore (see Lauritzen et al. 2012 DOI:10.1177/1094342011410088) scale1 = 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(sponge_del4_lev)/pmid(k)))) - nu_div_lev(k) = (1.0_r8-scale1)*nu_div+scale1*nu_div_max - if (sponge_del4_nu_fac.ne.1.0_r8) then + if (sponge_del4_nu_div_fac /= 1.0_r8) then + nu_div_lev(k) = (1.0_r8-scale1)*nu_div+scale1*nu_div_max + end if + if (sponge_del4_nu_fac /= 1.0_r8) then nu_lev(k) = (1.0_r8-scale1)*nu +scale1*nu_max nu_t_lev(k) = (1.0_r8-scale1)*nu_p +scale1*nu_max end if diff --git a/src/dynamics/se/dycore/prim_advance_mod.F90 b/src/dynamics/se/dycore/prim_advance_mod.F90 index c9f1ac194b..f9199bf7ce 100644 --- a/src/dynamics/se/dycore/prim_advance_mod.F90 +++ b/src/dynamics/se/dycore/prim_advance_mod.F90 @@ -78,7 +78,7 @@ subroutine prim_advance_exp(elem, fvm, deriv, hvcoord, hybrid,dt, tl, nets, net real (kind=r8) :: qwater(np,np,nlev,thermodynamic_active_species_num,nets:nete) integer :: qidx(thermodynamic_active_species_num) real (kind=r8) :: kappa(np,np,nlev,nets:nete) - call t_startf('prim_advance_exp') + nm1 = tl%nm1 n0 = tl%n0 np1 = tl%np1 @@ -257,7 +257,6 @@ subroutine prim_advance_exp(elem, fvm, deriv, hvcoord, hybrid,dt, tl, nets, net call omp_set_nested(.false.) - call t_stopf('prim_advance_exp') end subroutine prim_advance_exp @@ -976,6 +975,7 @@ subroutine compute_and_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& ! ! =================================== use dimensions_mod, only: np, nc, nlev, use_cslam + use control_mod, only: pgf_formulation use hybrid_mod, only: hybrid_t use element_mod, only: element_t use derivative_mod, only: derivative_t, divergence_sphere, gradient_sphere, vorticity_sphere @@ -1030,9 +1030,8 @@ subroutine compute_and_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& real (kind=r8) :: vtens1(np,np,nlev),vtens2(np,np,nlev),ttens(np,np,nlev) real (kind=r8) :: stashdp3d (np,np,nlev),tempdp3d(np,np), tempflux(nc,nc,4) real (kind=r8) :: ckk, term, T_v(np,np,nlev) - real (kind=r8), dimension(np,np,2) :: grad_exner - real (kind=r8), dimension(np,np,2) :: grad_exner_term - real (kind=r8), dimension(np,np,2) :: grad_logexner + real (kind=r8), dimension(np,np,2) :: pgf_term + real (kind=r8), dimension(np,np,2) :: grad_exner,grad_logexner real (kind=r8) :: T0,T1 real (kind=r8), dimension(np,np) :: theta_v @@ -1177,51 +1176,52 @@ subroutine compute_and_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& call gradient_sphere(Ephi(:,:),deriv,elem(ie)%Dinv,vtemp) density_inv(:,:) = R_dry(:,:,k)*T_v(:,:,k)/p_full(:,:,k) - if (dry_air_species_num==0) then - exner(:,:)=(p_full(:,:,k)/hvcoord%ps0)**kappa(:,:,k,ie) - theta_v(:,:)=T_v(:,:,k)/exner(:,:) - call gradient_sphere(exner(:,:),deriv,elem(ie)%Dinv,grad_exner) - - grad_exner_term(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,1) - grad_exner_term(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,2) + if (pgf_formulation==1) then + if (dry_air_species_num==0) then + exner(:,:)=(p_full(:,:,k)/hvcoord%ps0)**kappa(:,:,k,ie) + theta_v(:,:)=T_v(:,:,k)/exner(:,:) + call gradient_sphere(exner(:,:),deriv,elem(ie)%Dinv,grad_exner) + pgf_term(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,1) + pgf_term(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,2) + else + exner(:,:)=(p_full(:,:,k)/hvcoord%ps0)**kappa(:,:,k,ie) + theta_v(:,:)=T_v(:,:,k)/exner(:,:) + call gradient_sphere(exner(:,:),deriv,elem(ie)%Dinv,grad_exner) + call gradient_sphere(kappa(:,:,k,ie),deriv,elem(ie)%Dinv,grad_kappa_term) + suml = exner(:,:)*LOG(p_full(:,:,k)/hvcoord%ps0) + grad_kappa_term(:,:,1)=-suml*grad_kappa_term(:,:,1) + grad_kappa_term(:,:,2)=-suml*grad_kappa_term(:,:,2) + pgf_term(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,1)+grad_kappa_term(:,:,1)) + pgf_term(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,2)+grad_kappa_term(:,:,2)) + end if + ! balanced ref profile correction: + ! reference temperature profile (Simmons and Jiabin, 1991, QJRMS, Section 2a) + ! + ! Tref = T0+T1*Exner + ! T1 = .0065*Tref*Cp/g ! = ~191 + ! T0 = Tref-T1 ! = ~97 + ! + T1 = lapse_rate*Tref*cpair*rga + T0 = Tref-T1 + if (hvcoord%hybm(k)>0) then + !only apply away from constant pressure levels + call gradient_sphere(log(exner(:,:)),deriv,elem(ie)%Dinv,grad_logexner) + pgf_term(:,:,1)=pgf_term(:,:,1) + & + cpair*T0*(grad_logexner(:,:,1)-grad_exner(:,:,1)/exner(:,:)) + pgf_term(:,:,2)=pgf_term(:,:,2) + & + cpair*T0*(grad_logexner(:,:,2)-grad_exner(:,:,2)/exner(:,:)) + end if + elseif (pgf_formulation==2) then + pgf_term(:,:,1) = density_inv(:,:)*grad_p_full(:,:,1,k) + pgf_term(:,:,2) = density_inv(:,:)*grad_p_full(:,:,2,k) else - exner(:,:)=(p_full(:,:,k)/hvcoord%ps0)**kappa(:,:,k,ie) - theta_v(:,:)=T_v(:,:,k)/exner(:,:) - call gradient_sphere(exner(:,:),deriv,elem(ie)%Dinv,grad_exner) - - call gradient_sphere(kappa(:,:,k,ie),deriv,elem(ie)%Dinv,grad_kappa_term) - suml = exner(:,:)*LOG(p_full(:,:,k)/hvcoord%ps0) - grad_kappa_term(:,:,1)=-suml*grad_kappa_term(:,:,1) - grad_kappa_term(:,:,2)=-suml*grad_kappa_term(:,:,2) - - grad_exner_term(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,1)+grad_kappa_term(:,:,1)) - grad_exner_term(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,2)+grad_kappa_term(:,:,2)) + call endrun('ERROR: bad choice of pgf_formulation (must be 1 or 2)') end if - ! balanced ref profile correction: - ! reference temperature profile (Simmons and Jiabin, 1991, QJRMS, Section 2a) - ! - ! Tref = T0+T1*Exner - ! T1 = .0065*Tref*Cp/g ! = ~191 - ! T0 = Tref-T1 ! = ~97 - ! - T1 = lapse_rate*Tref*cpair*rga - T0 = Tref-T1 - - if (hvcoord%hybm(k)>0) then - call gradient_sphere(log(exner(:,:)),deriv,elem(ie)%Dinv,grad_logexner) - grad_exner_term(:,:,1)=grad_exner_term(:,:,1) + & - cpair*T0*(grad_logexner(:,:,1)-grad_exner(:,:,1)/exner(:,:)) - grad_exner_term(:,:,2)=grad_exner_term(:,:,2) + & - cpair*T0*(grad_logexner(:,:,2)-grad_exner(:,:,2)/exner(:,:)) - end if - - - do j=1,np do i=1,np - glnps1 = grad_exner_term(i,j,1) - glnps2 = grad_exner_term(i,j,2) + glnps1 = pgf_term(i,j,1) + glnps2 = pgf_term(i,j,2) v1 = elem(ie)%state%v(i,j,1,k,n0) v2 = elem(ie)%state%v(i,j,2,k,n0) diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 683565267c..ac46b8cc46 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -110,7 +110,7 @@ subroutine dyn_readnl(NLFileName) use control_mod, only: topology, variable_nsplit use control_mod, only: fine_ne, hypervis_power, hypervis_scaling use control_mod, only: max_hypervis_courant, statediag_numtrac,refined_mesh - use control_mod, only: molecular_diff + use control_mod, only: molecular_diff, pgf_formulation use control_mod, only: sponge_del4_nu_div_fac, sponge_del4_nu_fac, sponge_del4_lev use dimensions_mod, only: ne, npart use dimensions_mod, only: large_Courant_incr @@ -167,6 +167,7 @@ subroutine dyn_readnl(NLFileName) integer :: se_kmin_jet integer :: se_kmax_jet real(r8) :: se_molecular_diff + integer :: se_pgf_formulation namelist /dyn_se_inparm/ & se_fine_ne, & ! For refined meshes @@ -211,7 +212,8 @@ subroutine dyn_readnl(NLFileName) se_fvm_supercycling_jet, & se_kmin_jet, & se_kmax_jet, & - se_molecular_diff + se_molecular_diff, & + se_pgf_formulation !-------------------------------------------------------------------------- @@ -285,6 +287,7 @@ subroutine dyn_readnl(NLFileName) call MPI_bcast(se_kmin_jet, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(se_kmax_jet, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(se_molecular_diff, 1, mpi_real8, masterprocid, mpicom, ierr) + call MPI_bcast(se_pgf_formulation, 1, mpi_integer, masterprocid, mpicom, ierr) if (se_npes <= 0) then call endrun('dyn_readnl: ERROR: se_npes must be > 0') @@ -352,6 +355,7 @@ subroutine dyn_readnl(NLFileName) kmax_jet = se_kmax_jet variable_nsplit = .false. molecular_diff = se_molecular_diff + pgf_formulation = se_pgf_formulation if (fv_nphys > 0) then ! Use finite volume physics grid and CSLAM for tracer advection diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90 index 07444cbfa3..bb38481e39 100644 --- a/src/physics/cam/clubb_intr.F90 +++ b/src/physics/cam/clubb_intr.F90 @@ -17,18 +17,18 @@ module clubb_intr ! ! !----------------------------------------------------------------------------------------------------- ! - use shr_kind_mod, only: r8=>shr_kind_r8 + use shr_kind_mod, only: r8=>shr_kind_r8 use ppgrid, only: pver, pverp, pcols, begchunk, endchunk use phys_control, only: phys_getopts - use physconst, only: cpair, gravit, rga, latvap, latice, zvir, rh2o, karman + use physconst, only: cpair, gravit, rga, latvap, latice, zvir, rh2o use air_composition, only: rairv, cpairv use cam_history_support, only: max_fieldname_len - use spmd_utils, only: masterproc - use constituents, only: pcnst, cnst_add - use pbl_utils, only: calc_ustar, calc_obklen - use ref_pres, only: top_lev => trop_cloud_top_lev - use zm_conv_intr, only: zmconv_microp + use spmd_utils, only: masterproc + use constituents, only: pcnst, cnst_add + use pbl_utils, only: calc_ustar, calc_obklen + use ref_pres, only: top_lev => trop_cloud_top_lev + use zm_conv_intr, only: zmconv_microp #ifdef CLUBB_SGS use clubb_api_module, only: pdf_parameter, implicit_coefs_terms use clubb_api_module, only: clubb_config_flags_type, grid, stats, nu_vertical_res_dep @@ -77,7 +77,6 @@ module clubb_intr logical, public :: do_cldcool logical :: clubb_do_icesuper - #ifdef CLUBB_SGS type(clubb_config_flags_type), public :: clubb_config_flags real(r8), dimension(nparams), public :: clubb_params ! Adjustable CLUBB parameters (C1, C2 ...) @@ -306,7 +305,7 @@ module clubb_intr clubb_l_mono_flux_lim_vm, & ! Flag to turn on monotonic flux limiter for vm clubb_l_mono_flux_lim_spikefix, & ! Flag to implement monotonic flux limiter code that ! eliminates spurious drying tendencies at model top - clubb_l_intr_sfc_flux_smooth = .false. ! Add a locally calculated roughness to upwp and vpwp sfc fluxes + clubb_l_intr_sfc_flux_smooth = .false.! Add a locally calculated roughness to upwp and vpwp sfc fluxes ! Constant parameters logical, parameter, private :: & @@ -323,6 +322,7 @@ module clubb_intr logical :: clubb_do_liqsupersat = .false. logical :: clubb_do_energyfix = .true. logical :: history_budget + logical :: do_hb_above_clubb = .false. integer :: history_budget_histfile_num integer :: edsclr_dim ! Number of scalars to transport in CLUBB integer :: offset @@ -390,7 +390,8 @@ module clubb_intr qsatfac_idx, & ! subgrid cloud water saturation scaling factor ice_supersat_idx, & ! ice cloud fraction for SILHS rcm_idx, & ! Cloud water mixing ratio for SILHS - ztodt_idx ! physics timestep for SILHS + ztodt_idx,& ! physics timestep for SILHS + clubbtop_idx ! level index for CLUBB top ! Indices for microphysical covariance tendencies integer :: & @@ -467,13 +468,15 @@ subroutine clubb_register_cam( ) !------------------------------------------------ ! ! Add CLUBB fields to pbuf - use physics_buffer, only: pbuf_add_field, dtype_r8, dyn_time_lvls + use physics_buffer, only: pbuf_add_field, dtype_r8, dtype_i4, dyn_time_lvls use subcol_utils, only: subcol_get_scheme - + + !----- Begin Code ----- call phys_getopts( eddy_scheme_out = eddy_scheme, & deep_scheme_out = deep_scheme, & history_budget_out = history_budget, & - history_budget_histfile_num_out = history_budget_histfile_num ) + history_budget_histfile_num_out = history_budget_histfile_num, & + do_hb_above_clubb_out = do_hb_above_clubb) subcol_scheme = subcol_get_scheme() if (trim(subcol_scheme) == 'SILHS') then @@ -498,6 +501,9 @@ 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 + if (do_hb_above_clubb) then + call pbuf_add_field('clubbtop', 'physpkg', dtype_i4, (/pcols/), clubbtop_idx) + endif ! put pbuf_add calls here (see macrop_driver.F90 for sample) use indicies defined at top call pbuf_add_field('pblh', 'global', dtype_r8, (/pcols/), pblh_idx) @@ -806,7 +812,7 @@ subroutine clubb_readnl(nlfile) clubb_tridiag_solve_method, & clubb_up2_sfc_coef, & clubb_wpxp_L_thresh - + !----- Begin Code ----- ! Determine if we want clubb_history to be output @@ -1273,8 +1279,6 @@ subroutine clubb_ini_cam(pbuf2d) ! From CAM libraries use cam_history, only: addfld, add_default, horiz_only - use ref_pres, only: pref_mid - use hb_diff, only: init_hb_diff use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_mode_num_idx, rad_cnst_get_mam_mmr_idx use cam_abortutils, only: endrun @@ -1338,8 +1342,6 @@ subroutine clubb_ini_cam(pbuf2d) integer :: err_code ! Code for when CLUBB fails integer :: i, j, k, l ! Indices - integer :: ntop_eddy ! Top interface level to which eddy vertical diffusion is applied ( = 1 ) - integer :: nbot_eddy ! Bottom interface level to which eddy vertical diffusion is applied ( = pver ) integer :: nmodes, nspec, m integer :: ixq, ixcldice, ixcldliq, ixnumliq, ixnumice integer :: lptr @@ -1395,7 +1397,8 @@ subroutine clubb_ini_cam(pbuf2d) call phys_getopts(prog_modal_aero_out=prog_modal_aero, & history_amwg_out=history_amwg, & - history_clubb_out=history_clubb) + history_clubb_out=history_clubb, & + do_hb_above_clubb_out=do_hb_above_clubb) ! Select variables to apply tendencies back to CAM @@ -1644,17 +1647,6 @@ subroutine clubb_ini_cam(pbuf2d) write(iulog,'(a,i0,a)') " CLUBB configurable flags " call print_clubb_config_flags_api( iulog, clubb_config_flags ) ! Intent(in) end if - - ! ----------------------------------------------------------------- ! - ! Set-up HB diffusion. Only initialized to diagnose PBL depth ! - ! ----------------------------------------------------------------- ! - - ! Initialize eddy diffusivity module - - ntop_eddy = 1 ! if >1, must be <= nbot_molec - nbot_eddy = pver ! currently always pver - - call init_hb_diff( gravit, cpair, ntop_eddy, nbot_eddy, pref_mid, karman, eddy_scheme ) ! ----------------------------------------------------------------- ! ! Add output fields for the history files @@ -1958,6 +1950,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & physics_ptend_sum, physics_update, set_wet_to_dry use physics_buffer, only: pbuf_old_tim_idx, pbuf_get_field, physics_buffer_desc + use physics_buffer, only: pbuf_set_field use constituents, only: cnst_get_ind, cnst_type use camsrfexch, only: cam_in_t @@ -1966,7 +1959,6 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & use cam_logfile, only: iulog use tropopause, only: tropopause_findChemTrop use time_manager, only: get_nstep, is_first_restart_step - #ifdef CLUBB_SGS use hb_diff, only: pblintd use scamMOD, only: single_column,scm_clubb_iop_name @@ -3642,7 +3634,14 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & clubbtop(i) = clubbtop(i) + 1 end do end do - + ! + ! set pbuf field so that HB scheme is only applied above CLUBB top + ! + if (do_hb_above_clubb) then + call pbuf_set_field(pbuf, clubbtop_idx, clubbtop) + endif + + ! Compute integrals for static energy, kinetic energy, water vapor, and liquid water ! after CLUBB is called. This is for energy conservation purposes. se_a(:) = 0._r8 diff --git a/src/physics/cam/geopotential.F90 b/src/physics/cam/geopotential.F90 index 93e99644ac..7672fc92e0 100644 --- a/src/physics/cam/geopotential.F90 +++ b/src/physics/cam/geopotential.F90 @@ -1,4 +1,3 @@ - module geopotential !--------------------------------------------------------------------------------- @@ -39,7 +38,8 @@ subroutine geopotential_t( & !----------------------------------------------------------------------- use ppgrid, only : pcols -use air_composition, only: thermodynamic_active_species_num,thermodynamic_active_species_idx +use air_composition, only: thermodynamic_active_species_num +use air_composition, only: thermodynamic_active_species_idx, dry_air_species_num !------------------------------Arguments-------------------------------- ! ! Input arguments @@ -126,7 +126,7 @@ subroutine geopotential_t( & ! Compute factor for converting wet to dry mixing ratio (eq.7) ! qfac = 1.0_r8 - do idx = 1,thermodynamic_active_species_num + do idx = dry_air_species_num + 1,thermodynamic_active_species_num do k = 1,pver do i = 1,ncol qfac(i,k) = qfac(i,k)-q(i,k,thermodynamic_active_species_idx(idx)) @@ -137,7 +137,7 @@ subroutine geopotential_t( & ! Compute sum of dry water mixing ratios sum_dry_mixing_ratio = 1.0_r8 - do idx = 1,thermodynamic_active_species_num + do idx = dry_air_species_num + 1,thermodynamic_active_species_num do k = 1,pver do i = 1,ncol sum_dry_mixing_ratio(i,k) = sum_dry_mixing_ratio(i,k)& diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90 index f1b67d68a0..ba97978e72 100644 --- a/src/physics/cam/hb_diff.F90 +++ b/src/physics/cam/hb_diff.F90 @@ -35,6 +35,7 @@ module hb_diff ! Public interfaces public init_hb_diff public compute_hb_diff + public compute_hb_free_atm_diff public pblintd ! ! PBL limits @@ -131,7 +132,7 @@ end subroutine init_hb_diff !=============================================================================== - subroutine compute_hb_diff(lchnk, ncol, & + subroutine compute_hb_diff(ncol , & th ,t ,q ,z ,zi , & pmid ,u ,v ,taux ,tauy , & shflx ,qflx ,obklen ,ustar ,pblh , & @@ -139,8 +140,8 @@ subroutine compute_hb_diff(lchnk, ncol, & tpert ,qpert ,cldn ,ocnfrac ,tke , & ri , & eddy_scheme) - !----------------------------------------------------------------------- - ! + !----------------------------------------------------------------------- + ! ! Purpose: ! Interface routines for calcualtion and diatnostics of turbulence related ! coefficients @@ -155,7 +156,6 @@ subroutine compute_hb_diff(lchnk, ncol, & ! ! Input arguments ! - integer, intent(in) :: lchnk ! chunk index (for debug only) integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: th(pcols,pver) ! potential temperature [K] @@ -196,8 +196,8 @@ subroutine compute_hb_diff(lchnk, ncol, & real(r8) :: rrho(pcols) ! 1./bottom level density real(r8) :: wstar(pcols) ! convective velocity scale [m/s] real(r8) :: kqfs(pcols) ! kinematic surf constituent flux (kg/m2/s) - real(r8) :: khfs(pcols) ! kinimatic surface heat flux - real(r8) :: kbfs(pcols) ! surface buoyancy flux + real(r8) :: khfs(pcols) ! kinimatic surface heat flux + real(r8) :: kbfs(pcols) ! surface buoyancy flux real(r8) :: kvf(pcols,pverp) ! free atmospheric eddy diffsvty [m2/s] real(r8) :: s2(pcols,pver) ! shear squared real(r8) :: n2(pcols,pver) ! brunt vaisaila frequency @@ -236,7 +236,7 @@ subroutine compute_hb_diff(lchnk, ncol, & ! ! Get pbl exchange coefficients ! - call austausch_pbl(lchnk, ncol, & + call austausch_pbl(ncol , & z ,kvf ,kqfs ,khfs ,kbfs , & obklen ,ustar ,wstar ,pblh ,kvm , & kvh ,cgh ,cgs ,tpert ,qpert , & @@ -244,9 +244,95 @@ subroutine compute_hb_diff(lchnk, ncol, & ! kvq(:ncol,:) = kvh(:ncol,:) - - return end subroutine compute_hb_diff + + subroutine compute_hb_free_atm_diff(ncol, & + th ,t ,q ,z , & + pmid ,u ,v ,taux ,tauy , & + shflx ,qflx ,obklen ,ustar , & + kvm ,kvh ,kvq ,cgh ,cgs , & + ri ) + !----------------------------------------------------------------------- + ! + ! This is a version of compute_hb_diff that only computes free + ! atmosphere exchange (no PBL computations) + ! + ! Author: B. Stevens (rewrite August 2000) + ! Modified by Thomas Toniazzo and Peter H. Lauritzen (June 2023) + ! + !----------------------------------------------------------------------- + + use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm + + !------------------------------Arguments-------------------------------- + ! + ! Input arguments + ! + integer, intent(in) :: ncol ! number of atmospheric columns + + real(r8), intent(in) :: th(pcols,pver) ! potential temperature [K] + real(r8), intent(in) :: t(pcols,pver) ! temperature (used for density) + real(r8), intent(in) :: q(pcols,pver) ! specific humidity [kg/kg] + real(r8), intent(in) :: z(pcols,pver) ! height above surface [m] + real(r8), intent(in) :: u(pcols,pver) ! zonal velocity + real(r8), intent(in) :: v(pcols,pver) ! meridional velocity + real(r8), intent(in) :: taux(pcols) ! zonal stress [N/m2] + real(r8), intent(in) :: tauy(pcols) ! meridional stress [N/m2] + real(r8), intent(in) :: shflx(pcols) ! sensible heat flux + real(r8), intent(in) :: qflx(pcols) ! water vapor flux + real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures + ! + ! Output arguments + ! + real(r8), intent(out) :: kvm(pcols,pverp) ! eddy diffusivity for momentum [m2/s] + real(r8), intent(out) :: kvh(pcols,pverp) ! eddy diffusivity for heat [m2/s] + real(r8), intent(out) :: kvq(pcols,pverp) ! eddy diffusivity for constituents [m2/s] + real(r8), intent(out) :: cgh(pcols,pverp) ! counter-gradient term for heat [J/kg/m] + real(r8), intent(out) :: cgs(pcols,pverp) ! counter-gradient star (cg/flux) + real(r8), intent(out) :: ustar(pcols) ! surface friction velocity [m/s] + real(r8), intent(out) :: obklen(pcols) ! Obukhov length + real(r8), intent(out) :: ri(pcols,pver) ! richardson number: n2/s2 + ! + !---------------------------Local workspace----------------------------- + ! + real(r8) :: thv(pcols,pver) ! virtual potential temperature + real(r8) :: rrho(pcols) ! 1./bottom level density + real(r8) :: kqfs(pcols) ! kinematic surf constituent flux (kg/m2/s) + real(r8) :: khfs(pcols) ! kinimatic surface heat flux + real(r8) :: kbfs(pcols) ! surface buoyancy flux + real(r8) :: kvf(pcols,pverp) ! free atmospheric eddy diffsvty [m2/s] + real(r8) :: s2(pcols,pver) ! shear squared + real(r8) :: n2(pcols,pver) ! brunt vaisaila frequency + + ! virtual potential temperature + call virtem(ncol, (pver-ntop_turb+1), th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), thv(:ncol,ntop_turb:)) + + ! Compute ustar, Obukhov length, and kinematic surface fluxes. + call calc_ustar(ncol, t(:ncol,pver),pmid(:ncol,pver),taux(:ncol),tauy(:ncol), & + rrho(:ncol),ustar(:ncol)) + call calc_obklen(ncol, th(:ncol,pver), thv(:ncol,pver), qflx(:ncol), & + shflx(:ncol), rrho(:ncol), ustar(:ncol), & + khfs(:ncol), kqfs(:ncol), kbfs(:ncol), & + obklen(:ncol)) + ! Calculate s2, n2, and Richardson number. + call trbintd(ncol , & + thv ,z ,u ,v , & + s2 ,n2 ,ri ) + ! + ! Get free atmosphere exchange coefficients + ! + call austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, & + ml2, ri, s2, kvf) + + kvq(:ncol,:) = kvf(:ncol,:) + kvm(:ncol,:) = kvf(:ncol,:) + kvh(:ncol,:) = kvf(:ncol,:) + cgh(:ncol,:) = 0._r8 + cgs(:ncol,:) = 0._r8 + + end subroutine compute_hb_free_atm_diff + + ! !=============================================================================== subroutine trbintd(ncol , & @@ -488,7 +574,7 @@ subroutine pblintd(ncol , & end subroutine pblintd ! !=============================================================================== - subroutine austausch_pbl(lchnk ,ncol , & + subroutine austausch_pbl(ncol , & z ,kvf ,kqfs ,khfs ,kbfs , & obklen ,ustar ,wstar ,pblh ,kvm , & kvh ,cgh ,cgs ,tpert ,qpert , & @@ -524,7 +610,6 @@ subroutine austausch_pbl(lchnk ,ncol , & ! ! Input arguments ! - integer, intent(in) :: lchnk ! local chunk index (for debug only) integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: z(pcols,pver) ! height above surface [m] diff --git a/src/physics/cam/phys_control.F90 b/src/physics/cam/phys_control.F90 index d301c39948..92ccac1335 100644 --- a/src/physics/cam/phys_control.F90 +++ b/src/physics/cam/phys_control.F90 @@ -72,6 +72,8 @@ module phys_control logical :: history_chemspecies_srf = .false. logical :: do_clubb_sgs +logical :: do_hb_above_clubb = .false. ! enable HB vertical mixing above clubb top + ! Check validity of physics_state objects in physics_update. logical :: state_debug_checks = .false. @@ -136,7 +138,7 @@ subroutine phys_ctl_readnl(nlfile) do_clubb_sgs, state_debug_checks, use_hetfrz_classnuc, use_gw_oro, use_gw_front, & use_gw_front_igw, use_gw_convect_dp, use_gw_convect_sh, cld_macmic_num_steps, & offline_driver, convproc_do_aer, cam_snapshot_before_num, cam_snapshot_after_num, & - cam_take_snapshot_before, cam_take_snapshot_after, cam_physics_mesh, use_hemco + cam_take_snapshot_before, cam_take_snapshot_after, cam_physics_mesh, use_hemco, do_hb_above_clubb !----------------------------------------------------------------------------- if (masterproc) then @@ -199,6 +201,7 @@ subroutine phys_ctl_readnl(nlfile) call mpi_bcast(cam_take_snapshot_before, len(cam_take_snapshot_before), mpi_character, masterprocid, mpicom, ierr) call mpi_bcast(cam_take_snapshot_after, len(cam_take_snapshot_after), mpi_character, masterprocid, mpicom, ierr) call mpi_bcast(cam_physics_mesh, len(cam_physics_mesh), mpi_character, masterprocid, mpicom, ierr) + call mpi_bcast(do_hb_above_clubb, 1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(use_hemco, 1, mpi_logical, masterprocid, mpicom, ierr) use_spcam = ( cam_physpkg_is('spcam_sam1mom') & @@ -257,6 +260,12 @@ subroutine phys_ctl_readnl(nlfile) end if end if + ! do_hb_above_clubb requires that CLUBB is being used + if (do_hb_above_clubb .and. .not. do_clubb_sgs) then + write(iulog,*)'do_hb_above_clubb requires CLUBB to be active' + call endrun('do_hb_above_clubb incompatible with do_clubb_sgs = .false.') + endif + ! Macro/micro co-substepping support. if (cld_macmic_num_steps > 1) then if (microp_scheme /= "MG" .or. (macrop_scheme /= "park" .and. macrop_scheme /= "CLUBB_SGS")) then @@ -315,7 +324,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi cam_chempkg_out, prog_modal_aero_out, macrop_scheme_out, & do_clubb_sgs_out, use_spcam_out, state_debug_checks_out, cld_macmic_num_steps_out, & offline_driver_out, convproc_do_aer_out, cam_snapshot_before_num_out, cam_snapshot_after_num_out,& - cam_take_snapshot_before_out, cam_take_snapshot_after_out, physics_grid_out) + cam_take_snapshot_before_out, cam_take_snapshot_after_out, physics_grid_out, do_hb_above_clubb_out) !----------------------------------------------------------------------- ! Purpose: Return runtime settings ! deep_scheme_out : deep convection scheme @@ -363,6 +372,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi character(len=32), intent(out), optional :: cam_take_snapshot_before_out character(len=32), intent(out), optional :: cam_take_snapshot_after_out character(len=cl), intent(out), optional :: physics_grid_out + logical, intent(out), optional :: do_hb_above_clubb_out if ( present(deep_scheme_out ) ) deep_scheme_out = deep_scheme if ( present(shallow_scheme_out ) ) shallow_scheme_out = shallow_scheme @@ -402,6 +412,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi if ( present(cam_take_snapshot_before_out) ) cam_take_snapshot_before_out = cam_take_snapshot_before if ( present(cam_take_snapshot_after_out ) ) cam_take_snapshot_after_out = cam_take_snapshot_after if ( present(physics_grid_out ) ) physics_grid_out = cam_physics_mesh + if ( present(do_hb_above_clubb_out ) ) do_hb_above_clubb_out = do_hb_above_clubb end subroutine phys_getopts diff --git a/src/physics/cam/vertical_diffusion.F90 b/src/physics/cam/vertical_diffusion.F90 index dd2cc721f3..de47c75fb1 100644 --- a/src/physics/cam/vertical_diffusion.F90 +++ b/src/physics/cam/vertical_diffusion.F90 @@ -72,7 +72,6 @@ module vertical_diffusion use ref_pres, only : do_molec_diff, nbot_molec use phys_control, only : phys_getopts use time_manager, only : is_first_step - implicit none private save @@ -129,6 +128,9 @@ module vertical_diffusion integer :: taubljx_idx = -1 integer :: taubljy_idx = -1 +! pbuf field for clubb top above which HB (Holtslag Boville) scheme may be enabled +integer :: clubbtop_idx = -1 + logical :: diff_cnsrv_mass_check ! do mass conservation check logical :: do_iss ! switch for implicit turbulent surface stress logical :: prog_modal_aero = .false. ! set true if prognostic modal aerosols are present @@ -137,6 +139,7 @@ module vertical_diffusion logical :: do_pbl_diags = .false. logical :: waccmx_mode = .false. +logical :: do_hb_above_clubb = .false. contains @@ -389,6 +392,8 @@ subroutine vertical_diffusion_init(pbuf2d) if (masterproc) write(iulog, fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY =', ntop_eddy, 'NBOT_EDDY =', nbot_eddy + call phys_getopts(do_hb_above_clubb_out=do_hb_above_clubb) + select case ( eddy_scheme ) case ( 'diag_TKE', 'SPCAM_m2005' ) if( masterproc ) write(iulog,*) & @@ -401,6 +406,18 @@ subroutine vertical_diffusion_init(pbuf2d) call addfld('HB_ri', (/ 'lev' /), 'A', 'no', 'Richardson Number (HB Scheme), I' ) case ( 'CLUBB_SGS' ) do_pbl_diags = .true. + call init_hb_diff(gravit, cpair, ntop_eddy, nbot_eddy, pref_mid, karman, eddy_scheme) + ! + ! run HB scheme where CLUBB is not active when running cam_dev or cam6 physics + ! else init_hb_diff is called just for diagnostic purposes + ! + if (do_hb_above_clubb) then + if( masterproc ) then + write(iulog,*) 'vertical_diffusion_init: ' + write(iulog,*) 'eddy_diffusivity scheme where CLUBB is not active: Holtslag and Boville' + end if + call addfld('HB_ri', (/ 'lev' /), 'A', 'no', 'Richardson Number (HB Scheme), I' ) + end if end select ! ------------------------------------------- ! @@ -599,6 +616,11 @@ subroutine vertical_diffusion_init(pbuf2d) kvh_idx = pbuf_get_index('kvh') end if + if (do_hb_above_clubb) then + ! pbuf field denoting top of clubb + clubbtop_idx = pbuf_get_index('clubbtop') + end if + ! Initialization of some pbuf fields if (is_first_step()) then ! Initialization of pbuf fields tke, kvh, kvm are done in phys_inidat @@ -658,7 +680,7 @@ subroutine vertical_diffusion_tend( & use trb_mtn_stress_cam, only : trb_mtn_stress_tend use beljaars_drag_cam, only : beljaars_drag_tend use eddy_diff_cam, only : eddy_diff_tend - use hb_diff, only : compute_hb_diff + use hb_diff, only : compute_hb_diff, compute_hb_free_atm_diff use wv_saturation, only : qsat use molec_diff, only : compute_molec_diff, vd_lu_qdecomp use constituents, only : qmincg, qmin, cnst_type @@ -839,6 +861,7 @@ subroutine vertical_diffusion_tend( & real(r8) :: tauy(pcols) real(r8) :: shflux(pcols) real(r8) :: cflux(pcols,pcnst) + integer, pointer :: clubbtop(:) ! (pcols) logical :: lq(pcnst) @@ -979,42 +1002,67 @@ subroutine vertical_diffusion_tend( & ! Modification : We may need to use 'taux' instead of 'tautotx' here, for ! consistency with the previous HB scheme. - call compute_hb_diff( lchnk , ncol , & - th , state%t , state%q , state%zm , state%zi, & - state%pmid, state%u , state%v , tautotx , tautoty , & - cam_in%shf, cam_in%cflx(:,1), obklen , ustar , pblh , & - kvm , kvh , kvq , cgh , cgs , & - tpert , qpert , cldn , cam_in%ocnfrac , tke , & + + call compute_hb_diff(ncol , & + th , state%t , state%q , state%zm , state%zi , & + state%pmid, state%u , state%v , tautotx , tautoty , & + cam_in%shf, cam_in%cflx(:,1), obklen , ustar , pblh , & + kvm , kvh , kvq , cgh , cgs , & + tpert , qpert , cldn , cam_in%ocnfrac , tke , & ri , & - eddy_scheme ) + eddy_scheme) call outfld( 'HB_ri', ri, pcols, lchnk ) case ( 'CLUBB_SGS' ) - - ! CLUBB has only a bare-bones placeholder here. If using CLUBB, the - ! PBL diffusion will happen before coupling, so vertical_diffusion - ! is only handling other things, e.g. some boundary conditions, tms, - ! and molecular diffusion. - - call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol)) - - call calc_ustar( ncol, state%t(:ncol,pver), state%pmid(:ncol,pver), & - cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol), ustar(:ncol)) - ! Use actual qflux, not lhf/latvap as was done previously - call calc_obklen( ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), & - cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), & - khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol)) - - ! These tendencies all applied elsewhere. - kvm = 0._r8 - kvh = 0._r8 - kvq = 0._r8 - - ! Not defined since PBL is not actually running here. - cgh = 0._r8 - cgs = 0._r8 - + ! + ! run HB scheme where CLUBB is not active when running cam_dev + ! + if (do_hb_above_clubb) then + call compute_hb_free_atm_diff( ncol , & + th , state%t , state%q , state%zm , & + state%pmid, state%u , state%v , tautotx , tautoty , & + cam_in%shf, cam_in%cflx(:,1), obklen , ustar , & + kvm , kvh , kvq , cgh , cgs , & + ri ) + + call pbuf_get_field(pbuf, clubbtop_idx, clubbtop) + ! + ! zero out HB where CLUBB is active + ! + do i=1,ncol + do k=clubbtop(i),pverp + kvm(i,k) = 0.0_r8 + kvh(i,k) = 0.0_r8 + kvq(i,k) = 0.0_r8 + cgs(i,k) = 0.0_r8 + cgh(i,k) = 0.0_r8 + end do + end do + + call outfld( 'HB_ri', ri, pcols, lchnk ) + else + ! CLUBB has only a bare-bones placeholder here. If using CLUBB, the + ! PBL diffusion will happen before coupling, so vertical_diffusion + ! is only handling other things, e.g. some boundary conditions, tms, + ! and molecular diffusion. + + call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol)) + + call calc_ustar( ncol, state%t(:ncol,pver), state%pmid(:ncol,pver), & + cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol), ustar(:ncol)) + ! Use actual qflux, not lhf/latvap as was done previously + call calc_obklen( ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), & + cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), & + khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol)) + ! These tendencies all applied elsewhere. + kvm = 0._r8 + kvh = 0._r8 + kvq = 0._r8 + ! Not defined since PBL is not actually running here. + cgh = 0._r8 + cgs = 0._r8 + end if end select call outfld( 'ustar', ustar(:), pcols, lchnk )