From 4210a61acbdfd197bf00a48728cd575fb1e0dcda Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 15 Feb 2024 17:08:59 -0500 Subject: [PATCH 01/45] misc updates for the pam_debug routines --- .../eam/src/physics/crm/pam/pam_debug.h | 117 ++++++++++++------ 1 file changed, 82 insertions(+), 35 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_debug.h b/components/eam/src/physics/crm/pam/pam_debug.h index 5186d51c6a9d..f4d8cdb1eda2 100644 --- a/components/eam/src/physics/crm/pam/pam_debug.h +++ b/components/eam/src/physics/crm/pam/pam_debug.h @@ -16,29 +16,33 @@ inline void pam_debug_init( pam::PamCoupler &coupler ) { auto nens = coupler.get_option("ncrms"); //------------------------------------------------------------------------------------------------ dm_device.register_and_allocate("debug_save_temp", "saved temp for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); - dm_device.register_and_allocate("debug_save_rhod", "saved rhod for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_rhov", "saved rhov for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_rhoc", "saved rhoc for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_rhoi", "saved rhoi for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); + dm_device.register_and_allocate("debug_save_uvel", "saved uvel for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); + dm_device.register_and_allocate("debug_save_wvel", "saved wvel for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); auto debug_save_temp = dm_device.get("debug_save_temp"); - auto debug_save_rhod = dm_device.get("debug_save_rhod"); auto debug_save_rhov = dm_device.get("debug_save_rhov"); auto debug_save_rhoc = dm_device.get("debug_save_rhoc"); auto debug_save_rhoi = dm_device.get("debug_save_rhoi"); + auto debug_save_uvel = dm_device.get("debug_save_uvel"); + auto debug_save_wvel = dm_device.get("debug_save_wvel"); //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp"); - auto rhod = dm_device.get("density_dry"); auto rhov = dm_device.get("water_vapor"); auto rhoc = dm_device.get("cloud_water"); auto rhoi = dm_device.get("ice"); + auto uvel = dm_device.get("uvel"); + auto wvel = dm_device.get("wvel"); //------------------------------------------------------------------------------------------------ parallel_for("copy data to saved debug variables", SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { debug_save_temp(k,j,i,iens) = temp(k,j,i,iens); - debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens); debug_save_rhov(k,j,i,iens) = rhov(k,j,i,iens); debug_save_rhoc(k,j,i,iens) = rhoc(k,j,i,iens); debug_save_rhoi(k,j,i,iens) = rhoi(k,j,i,iens); + debug_save_uvel(k,j,i,iens) = uvel(k,j,i,iens); + debug_save_wvel(k,j,i,iens) = wvel(k,j,i,iens); }); //------------------------------------------------------------------------------------------------ } @@ -48,81 +52,124 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; auto &dm_device = coupler.get_data_manager_device_readwrite(); - auto &dm_host = coupler.get_data_manager_host_readwrite(); + auto &dm_host = coupler.get_data_manager_host_readonly(); auto nx = coupler.get_option("crm_nx"); auto ny = coupler.get_option("crm_ny"); auto nz = coupler.get_option("crm_nz"); auto nens = coupler.get_option("ncrms"); - auto temp = dm_device.get("temp"); - auto rhod = dm_device.get("density_dry"); - auto rhov = dm_device.get("water_vapor"); - auto rhoc = dm_device.get("cloud_water"); - auto rhoi = dm_device.get("ice"); + auto temp = dm_device.get("temp"); + auto rhov = dm_device.get("water_vapor"); + auto rhoc = dm_device.get("cloud_water"); + auto rhoi = dm_device.get("ice"); + auto uvel = dm_device.get("uvel"); + auto wvel = dm_device.get("wvel"); auto debug_save_temp = dm_device.get("debug_save_temp"); - auto debug_save_rhod = dm_device.get("debug_save_rhod"); auto debug_save_rhov = dm_device.get("debug_save_rhov"); auto debug_save_rhoc = dm_device.get("debug_save_rhoc"); auto debug_save_rhoi = dm_device.get("debug_save_rhoi"); - real grav = 9.80616; - auto input_phis = dm_host.get("input_phis").createDeviceCopy(); + auto debug_save_uvel = dm_device.get("debug_save_uvel"); + auto debug_save_wvel = dm_device.get("debug_save_wvel"); auto lat = dm_host.get("latitude" ).createDeviceCopy(); auto lon = dm_host.get("longitude" ).createDeviceCopy(); + auto input_phis = dm_host.get("input_phis").createDeviceCopy(); + real grav = 9.80616; //------------------------------------------------------------------------------------------------ // Check for invalid values parallel_for("", SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { + auto phis = input_phis(iens)/grav; // Check for NaNs const auto is_nan_t_atm = isnan( temp(k,j,i,iens) ); - const auto is_nan_r_atm = isnan( rhod(k,j,i,iens) ); const auto is_nan_q_atm = isnan( rhov(k,j,i,iens) ); - if ( is_nan_t_atm || is_nan_r_atm || is_nan_q_atm ) { - auto phis = input_phis(iens)/grav; - printf("PAM-DEBUG nan-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g \n", + if ( is_nan_t_atm || is_nan_q_atm ) { + printf("PAM-DEBUG nan-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, temp(k,j,i,iens), - rhod(k,j,i,iens), rhov(k,j,i,iens), rhoc(k,j,i,iens), rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), debug_save_temp(k,j,i,iens), - debug_save_rhod(k,j,i,iens), debug_save_rhov(k,j,i,iens), debug_save_rhoc(k,j,i,iens), - debug_save_rhoi(k,j,i,iens) + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) ); } // Check for negative values const auto is_neg_t_atm = temp(k,j,i,iens)<0; - const auto is_neg_r_atm = rhod(k,j,i,iens)<0; const auto is_neg_q_atm = rhov(k,j,i,iens)<0; - if ( is_neg_t_atm || is_neg_r_atm || is_neg_q_atm ) { - auto phis = input_phis(iens)/grav; - printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g \n", + if ( is_neg_t_atm || is_neg_q_atm ) { + printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, temp(k,j,i,iens), - rhod(k,j,i,iens), rhov(k,j,i,iens), rhoc(k,j,i,iens), rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), debug_save_temp(k,j,i,iens), - debug_save_rhod(k,j,i,iens), debug_save_rhov(k,j,i,iens), debug_save_rhoc(k,j,i,iens), - debug_save_rhoi(k,j,i,iens) + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) + ); + } + // Check for low temperature + const auto is_low_t = temp(k,j,i,iens)<100; + if ( is_low_t ) { + printf("PAM-DEBUG low-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + nstep,id,k,i,iens,lat(iens),lon(iens),phis, + temp(k,j,i,iens), + rhov(k,j,i,iens), + rhoc(k,j,i,iens), + rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), + debug_save_temp(k,j,i,iens), + debug_save_rhov(k,j,i,iens), + debug_save_rhoc(k,j,i,iens), + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) + ); + } + // Check for large vertical velocity + const auto is_large_pos_w = wvel(k,j,i,iens)> 40; + const auto is_large_neg_w = wvel(k,j,i,iens)<-40; + if ( is_large_pos_w || is_large_neg_w ) { + printf("PAM-DEBUG large-W - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + nstep,id,k,i,iens,lat(iens),lon(iens),phis, + temp(k,j,i,iens), + rhov(k,j,i,iens), + rhoc(k,j,i,iens), + rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), + debug_save_temp(k,j,i,iens), + debug_save_rhov(k,j,i,iens), + debug_save_rhoc(k,j,i,iens), + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) ); } // update saved previous values debug_save_temp(k,j,i,iens) = temp(k,j,i,iens); - debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens); debug_save_rhov(k,j,i,iens) = rhov(k,j,i,iens); debug_save_rhoc(k,j,i,iens) = rhoc(k,j,i,iens); debug_save_rhoi(k,j,i,iens) = rhoi(k,j,i,iens); + debug_save_uvel(k,j,i,iens) = uvel(k,j,i,iens); + debug_save_wvel(k,j,i,iens) = wvel(k,j,i,iens); }); //------------------------------------------------------------------------------------------------ } // print the min and max of PAM state variables to help look for problems // void print_state_min_max( pam::PamCoupler &coupler, std::string id ) { - // auto &dm_device = coupler.get_data_manager_device_readwrite(); + // auto &dm_device = coupler.get_data_manager_device_readonly(); // auto nz = coupler.get_option("crm_nz"); // auto nx = coupler.get_option("crm_nx"); // auto ny = coupler.get_option("crm_ny"); @@ -184,18 +231,18 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { void pam_debug_print_state( pam::PamCoupler &coupler, int id ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; - auto &dm_device = coupler.get_data_manager_device_readwrite(); + auto &dm_device = coupler.get_data_manager_device_readonly(); auto nz = coupler.get_option("crm_nz"); auto nx = coupler.get_option("crm_nx"); auto ny = coupler.get_option("crm_ny"); auto nens = coupler.get_option("ncrms"); // auto pmid = coupler.compute_pressure_array(); // auto zmid = dm_device.get("vertical_midpoint_height" ); - auto temp = dm_device.get("temp"); - auto rho_v = dm_device.get("water_vapor"); - auto rho_c = dm_device.get("cloud_water"); - auto rho_d = dm_device.get("density_dry"); - auto rho_i = dm_device.get("ice"); + auto temp = dm_device.get("temp"); + auto rho_v = dm_device.get("water_vapor"); + auto rho_c = dm_device.get("cloud_water"); + auto rho_d = dm_device.get("density_dry"); + auto rho_i = dm_device.get("ice"); // parallel_for("", SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { parallel_for("pam_debug_print_state", SimpleBounds<2>(nz,nx), YAKL_LAMBDA (int k, int i) { int j = 0; From 0def774fe8d97cc5a890f48604062dcabea96635 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 15 Feb 2024 17:19:56 -0500 Subject: [PATCH 02/45] add MMF_PAM_dyn_per_phys namelist parameter --- components/eam/bld/build-namelist | 3 +++ .../eam/bld/namelist_files/namelist_defaults_eam.xml | 2 +- components/eam/bld/namelist_files/namelist_definition.xml | 6 ++++++ components/eam/src/physics/cam/phys_control.F90 | 8 ++++++-- components/eam/src/physics/crm/crm_physics.F90 | 7 +++++++ 5 files changed, 23 insertions(+), 3 deletions(-) diff --git a/components/eam/bld/build-namelist b/components/eam/bld/build-namelist index 5a29d02bba48..c94bb5054d57 100755 --- a/components/eam/bld/build-namelist +++ b/components/eam/bld/build-namelist @@ -3885,6 +3885,9 @@ if ($cfg->get('use_MMF')) { # MMF CRM domain orientation add_default($nl, 'MMF_orientation_angle'); + + # PAM dynamics sub-stepping parameter + add_default($nl, 'MMF_PAM_dyn_per_phys'); } add_default($nl, 'do_aerocom_ind3'); diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index c5a61b500d21..4768f1799a4e 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -879,7 +879,7 @@ 0 .false. .true. - + 2 90.0 0.0 diff --git a/components/eam/bld/namelist_files/namelist_definition.xml b/components/eam/bld/namelist_files/namelist_definition.xml index fb5049df2110..dc63a488a6b6 100644 --- a/components/eam/bld/namelist_files/namelist_definition.xml +++ b/components/eam/bld/namelist_files/namelist_definition.xml @@ -4769,6 +4769,12 @@ Defaults: 3D CRM: 0 + +Number of PAM dycor calls per physics time step. +Default: 2 + + Wavenumber cutoff for filtered MMF variance transport. diff --git a/components/eam/src/physics/cam/phys_control.F90 b/components/eam/src/physics/cam/phys_control.F90 index c4d4b7fba7fa..41cd05e2714e 100644 --- a/components/eam/src/physics/cam/phys_control.F90 +++ b/components/eam/src/physics/cam/phys_control.F90 @@ -72,6 +72,7 @@ module phys_control logical :: use_crm_accel = .false. ! true => use MMF CRM mean-state acceleration (MSA) real(r8) :: crm_accel_factor = 2.D0 ! CRM acceleration factor logical :: crm_accel_uv = .true. ! true => apply MMF CRM MSA to momentum fields +integer :: MMF_PAM_dyn_per_phys = 2 ! PAM CRM dynamics steps per CRM physics steps logical :: use_subcol_microp = .false. ! if .true. then use sub-columns in microphysics @@ -229,7 +230,7 @@ subroutine phys_ctl_readnl(nlfile) eddy_scheme, microp_scheme, macrop_scheme, radiation_scheme, srf_flux_avg, & MMF_microphysics_scheme, MMF_orientation_angle, use_MMF, use_ECPP, & use_MMF_VT, MMF_VT_wn_max, use_MMF_ESMT, & - use_crm_accel, crm_accel_factor, crm_accel_uv, & + use_crm_accel, crm_accel_factor, crm_accel_uv, MMF_PAM_dyn_per_phys, & use_subcol_microp, atm_dep_flux, history_amwg, history_verbose, history_vdiag, & get_presc_aero_data,history_aerosol, history_aero_optics, & is_output_interactive_volc, & @@ -314,6 +315,7 @@ subroutine phys_ctl_readnl(nlfile) call mpibcast(use_crm_accel, 1 , mpilog, 0, mpicom) call mpibcast(crm_accel_factor, 1 , mpir8, 0, mpicom) call mpibcast(crm_accel_uv, 1 , mpilog, 0, mpicom) + call mpibcast(MMF_PAM_dyn_per_phys, 1 , mpiint, 0, mpicom) call mpibcast(use_subcol_microp, 1 , mpilog, 0, mpicom) call mpibcast(atm_dep_flux, 1 , mpilog, 0, mpicom) call mpibcast(history_amwg, 1 , mpilog, 0, mpicom) @@ -598,7 +600,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, & prog_modal_aero_out, macrop_scheme_out, ideal_phys_option_out, & use_MMF_out, use_ECPP_out, MMF_microphysics_scheme_out, & MMF_orientation_angle_out, use_MMF_VT_out, MMF_VT_wn_max_out, use_MMF_ESMT_out, & - use_crm_accel_out, crm_accel_factor_out, crm_accel_uv_out, & + use_crm_accel_out, crm_accel_factor_out, crm_accel_uv_out, MMF_PAM_dyn_per_phys_out, & do_clubb_sgs_out, do_shoc_sgs_out, do_tms_out, state_debug_checks_out, & linearize_pbl_winds_out, export_gustiness_out, & do_aerocom_ind3_out, & @@ -641,6 +643,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, & logical, intent(out), optional :: use_crm_accel_out real(r8), intent(out), optional :: crm_accel_factor_out logical, intent(out), optional :: crm_accel_uv_out + integer, intent(out), optional :: MMF_PAM_dyn_per_phys_out logical, intent(out), optional :: use_subcol_microp_out logical, intent(out), optional :: atm_dep_flux_out logical, intent(out), optional :: history_amwg_out @@ -745,6 +748,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, & if ( present(use_crm_accel_out ) ) use_crm_accel_out = use_crm_accel if ( present(crm_accel_factor_out ) ) crm_accel_factor_out = crm_accel_factor if ( present(crm_accel_uv_out ) ) crm_accel_uv_out = crm_accel_uv + if ( present(MMF_PAM_dyn_per_phys ) ) MMF_PAM_dyn_per_phys_out = MMF_PAM_dyn_per_phys if ( present(use_subcol_microp_out ) ) use_subcol_microp_out = use_subcol_microp if ( present(macrop_scheme_out ) ) macrop_scheme_out = macrop_scheme diff --git a/components/eam/src/physics/crm/crm_physics.F90 b/components/eam/src/physics/crm/crm_physics.F90 index 66e5b57e28eb..31654e2bb7ac 100644 --- a/components/eam/src/physics/crm/crm_physics.F90 +++ b/components/eam/src/physics/crm/crm_physics.F90 @@ -620,6 +620,8 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, logical :: use_MMF_VT_tmp ! flag for MMF variance transport (for Fortran CRM) logical :: use_MMF_ESMT_tmp ! flag for MMF scalar momentum transport (for Fortran CRM) integer :: MMF_VT_wn_max ! wavenumber cutoff for filtered variance transport + + integer :: MMF_PAM_dyn_per_phys ! PAM CRM dynamics steps per CRM physics steps real(r8) :: tmp_e_sat ! temporary saturation vapor pressure real(r8) :: tmp_q_sat ! temporary saturation specific humidity @@ -732,6 +734,9 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, use_crm_accel = use_crm_accel_tmp crm_accel_uv = crm_accel_uv_tmp + ! PAM namelist options + call phys_getopts(MMF_PAM_dyn_per_phys_out = MMF_PAM_dyn_per_phys) + nstep = get_nstep() itim = pbuf_old_tim_idx() ! "Old" pbuf time index (what does all this mean?) @@ -1459,6 +1464,8 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, call pam_set_option('enable_physics_tend_stats', .false. ) + call pam_set_option('crm_dyn_per_phys', MMF_PAM_dyn_per_phys ) + call pam_set_option('is_first_step', (nstep<=1) ) call pam_set_option('is_restart', (nsrest>0) ) call pam_set_option('am_i_root', masterproc ) From 9a25fc689c29fcbcf6725f09859fcb91c3c46a2d Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 15 Feb 2024 17:58:22 -0500 Subject: [PATCH 03/45] pam driver updates --- .../eam/src/physics/crm/pam/pam_driver.cpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index f131bbd44f09..128974a856e4 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -24,7 +24,7 @@ #include "p3_f90.hpp" #include "pam_debug.h" -bool constexpr enable_check_state = false; +bool constexpr enable_check_state = true; extern "C" void pam_driver() { //------------------------------------------------------------------------------------------------ @@ -50,12 +50,6 @@ extern "C" void pam_driver() { //------------------------------------------------------------------------------------------------ // set various coupler options coupler.set_option("gcm_physics_dt",gcm_dt); - #ifdef MMF_PAM_DPP - // this is leftover from debugging, but it might still be useful for testing values of crm_per_phys - coupler.set_option("crm_per_phys",MMF_PAM_DPP); - #else - coupler.set_option("crm_per_phys",2); // # of PAM-C dynamics steps per physics - #endif coupler.set_option("sponge_num_layers",crm_nz*0.3); // depth of sponge layer coupler.set_option("sponge_time_scale",60); // minimum damping timescale at top coupler.set_option("crm_acceleration_ceaseflag",false); @@ -143,6 +137,9 @@ extern "C" void pam_driver() { dycore.declare_current_profile_as_hydrostatic(coupler,/*use_gcm_data=*/true); #endif + #ifdef MMF_DISABLE_DENSITY_RECALL + do_density_save_recall = false; + #endif //------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------ @@ -168,8 +165,9 @@ extern "C" void pam_driver() { // run a PAM time step coupler.run_module( "apply_gcm_forcing_tendencies" , modules::apply_gcm_forcing_tendencies ); - coupler.run_module( "radiation" , [&] (pam::PamCoupler &coupler) {rad .timeStep(coupler);} ); if (enable_check_state) { pam_debug_check_state(coupler, 2, nstep); } + coupler.run_module( "radiation" , [&] (pam::PamCoupler &coupler) {rad .timeStep(coupler);} ); + if (enable_check_state) { pam_debug_check_state(coupler, 3, nstep); } // Dynamics if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } @@ -177,13 +175,13 @@ extern "C" void pam_driver() { coupler.run_module( "dycore", [&] (pam::PamCoupler &coupler) {dycore.timeStep(coupler);} ); if (do_density_save_recall) { pam_state_recall_dry_density(coupler); } if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"dycor"); } - if (enable_check_state) { pam_debug_check_state(coupler, 3, nstep); } + if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } // Sponge layer damping if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } coupler.run_module( "sponge_layer", modules::sponge_layer ); if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sponge"); } - if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } + // if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } // Apply hyperdiffusion to account for lack of horizontal mixing in SHOC pam_hyperdiffusion(coupler); From 0e93989e805d8ae7ecd59d2ffb955d0c10bcba30 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 15 Feb 2024 17:59:05 -0500 Subject: [PATCH 04/45] add temporary MMF_PAM_HDT for debugging --- components/eam/src/physics/crm/pam/pam_hyperdiffusion.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h index 83dc83252c68..4452c583afe3 100644 --- a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h +++ b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h @@ -25,7 +25,11 @@ inline void pam_hyperdiffusion( pam::PamCoupler &coupler ) { auto uvel = dm_device.get("uvel" ); auto vvel = dm_device.get("vvel" ); //------------------------------------------------------------------------------------------------ + #ifdef MMF_PAM_HDT + real constexpr hd_timescale = MMF_PAM_HDT; // damping time scale [sec] + #else real constexpr hd_timescale = 10.0; // damping time scale [sec] + #endif //------------------------------------------------------------------------------------------------ real4d hd_temp("hd_temp",nz,ny,nx,nens); real4d hd_rhod("hd_rhod",nz,ny,nx,nens); From 047f01af0a3025d1b64d63afce65a90df2bededb Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 15 Feb 2024 17:59:38 -0500 Subject: [PATCH 05/45] Add temporary MMF_ALT_DENSITY_RECALL for debugging --- .../eam/src/physics/crm/pam/pam_state.h | 29 +++++++++++++++---- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_state.h b/components/eam/src/physics/crm/pam/pam_state.h index e8f6e4ab668a..34cc42a4a788 100644 --- a/components/eam/src/physics/crm/pam/pam_state.h +++ b/components/eam/src/physics/crm/pam/pam_state.h @@ -261,6 +261,7 @@ inline void pam_state_save_dry_density( pam::PamCoupler &coupler ) { inline void pam_state_recall_dry_density( pam::PamCoupler &coupler ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; + using yakl::atomicAdd; auto &dm_device = coupler.get_data_manager_device_readwrite(); auto nens = coupler.get_option("ncrms"); auto nz = coupler.get_option("crm_nz"); @@ -269,11 +270,29 @@ inline void pam_state_recall_dry_density( pam::PamCoupler &coupler ) { auto crm_rho_d = dm_device.get("density_dry"); auto tmp_rho_d = dm_device.get("density_dry_save"); //------------------------------------------------------------------------------------------------ - parallel_for( "Recall CRM dry density", - SimpleBounds<4>(nz,ny,nx,nens), - YAKL_LAMBDA (int k_crm, int j, int i, int iens) { - crm_rho_d(k_crm,j,i,iens) = tmp_rho_d(k_crm,j,i,iens); - }); + #ifdef MMF_ALT_DENSITY_RECALL + // initialize horizontal mean of previous CRM dry density + real2d old_hmean_rho_d("old_hmean_rho_d" ,nz,nens); + real2d new_hmean_rho_d("new_hmean_rho_d" ,nz,nens); + parallel_for(SimpleBounds<2>(nz,nens), YAKL_LAMBDA (int k, int iens) { + old_hmean_rho_d(k,iens) = 0; + new_hmean_rho_d(k,iens) = 0; + }); + real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions + // calculate horizontal mean of previous CRM dry density + parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { + atomicAdd( old_hmean_rho_d(k,iens), tmp_rho_d(k,j,i,iens) * r_nx_ny ); + atomicAdd( new_hmean_rho_d(k,iens), crm_rho_d(k,j,i,iens) * r_nx_ny ); + }); + // replace horizontal mean dry density with prevoius value + parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { + crm_rho_d(k,j,i,iens) = crm_rho_d(k,j,i,iens) - new_hmean_rho_d(k,iens) + old_hmean_rho_d(k,iens); + }); + #else + parallel_for( "Recall CRM dry density",SimpleBounds<4>(nz,ny,nx,nens),YAKL_LAMBDA (int k_crm, int j, int i, int iens) { + crm_rho_d(k_crm,j,i,iens) = tmp_rho_d(k_crm,j,i,iens); + }); + #endif //------------------------------------------------------------------------------------------------ } From 6df3ae38c9cbb25897ef9341639c8d7a220ef957 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 22 Feb 2024 08:01:34 -0800 Subject: [PATCH 06/45] fix typo --- components/eam/src/physics/cam/phys_control.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/cam/phys_control.F90 b/components/eam/src/physics/cam/phys_control.F90 index 543fd29753f2..b87ea442c4e7 100644 --- a/components/eam/src/physics/cam/phys_control.F90 +++ b/components/eam/src/physics/cam/phys_control.F90 @@ -745,7 +745,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, & if ( present(use_crm_accel_out ) ) use_crm_accel_out = use_crm_accel if ( present(crm_accel_factor_out ) ) crm_accel_factor_out = crm_accel_factor if ( present(crm_accel_uv_out ) ) crm_accel_uv_out = crm_accel_uv - if ( present(MMF_PAM_dyn_per_phys ) ) MMF_PAM_dyn_per_phys_out = MMF_PAM_dyn_per_phys + if ( present(MMF_PAM_dyn_per_phys_out) ) MMF_PAM_dyn_per_phys_out = MMF_PAM_dyn_per_phys if ( present(use_subcol_microp_out ) ) use_subcol_microp_out = use_subcol_microp if ( present(macrop_scheme_out ) ) macrop_scheme_out = macrop_scheme From 8b89b6fc87e9a0dbbd272d786e18f550b314a151 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Mon, 26 Feb 2024 16:05:17 -0800 Subject: [PATCH 07/45] adjust default MMF2 configuration --- components/eam/cime_config/config_component.xml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/components/eam/cime_config/config_component.xml b/components/eam/cime_config/config_component.xml index c7585072ac94..16ff080ae839 100755 --- a/components/eam/cime_config/config_component.xml +++ b/components/eam/cime_config/config_component.xml @@ -75,10 +75,11 @@ -crm pam -pam_dycor spam -crm_dt 8 -crm pam -pam_dycor awfl -crm_dt 8 -use_MMF -nlev 60 -crm_nz 50 - -crm_dx 2000 -crm_nx 64 -crm_ny 1 + -crm_dx 2000 -crm_nx 64 -crm_ny 1 -crm_nx_rad 4 -crm_ny_rad 1 + -crm_dx 2000 -crm_nx 45 -crm_ny 1 -crm_nx_rad 5 -crm_ny_rad 1 -MMF_microphysics_scheme sam1mom -chem none -MMF_microphysics_scheme p3 -chem none - -crm_nx_rad 4 -crm_ny_rad 1 -rad rrtmgp -rrtmgpxx + -rad rrtmgp -rrtmgpxx -use_MMF_VT -use_MMF_ESMT -aquaplanet From ef9d3281c9051dadc2627f6abcf9d4e0308647bd Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 27 Feb 2024 08:18:08 -0800 Subject: [PATCH 08/45] add PAM coupler options to override SPAM parameter defaults --- components/eam/src/physics/crm/pam/pam_driver.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index 128974a856e4..c325945f8803 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -54,6 +54,14 @@ extern "C" void pam_driver() { coupler.set_option("sponge_time_scale",60); // minimum damping timescale at top coupler.set_option("crm_acceleration_ceaseflag",false); //------------------------------------------------------------------------------------------------ + // coupler options for SPAM dycor + coupler.set_option("spam_couple_wind_exact_inverse",true); + coupler.set_option("spam_clip_negative_densities",true); + coupler.set_option("spam_clip_vertical_velocities",true); + coupler.set_option("spam_adjust_crm_per_phys_using_vert_cfl",true); + coupler.set_option("spam_target_cfl",0.7); + coupler.set_option("spam_max_w",50.0); + //------------------------------------------------------------------------------------------------ // Allocate the coupler state and retrieve host/device data managers coupler.allocate_coupler_state( crm_nz , crm_ny , crm_nx , nens ); From 644bf4c00c6dc3295fc3996455318e0fdd6f9c51 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 27 Feb 2024 08:24:59 -0800 Subject: [PATCH 09/45] update pam_state_recall_dry_density --- .../eam/src/physics/crm/pam/pam_state.h | 54 +++++++++++-------- 1 file changed, 31 insertions(+), 23 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_state.h b/components/eam/src/physics/crm/pam/pam_state.h index 34cc42a4a788..3266899ba217 100644 --- a/components/eam/src/physics/crm/pam/pam_state.h +++ b/components/eam/src/physics/crm/pam/pam_state.h @@ -258,6 +258,13 @@ inline void pam_state_save_dry_density( pam::PamCoupler &coupler ) { // recall CRM dry density saved previously - only for anelastic case +// Note - The need for this arises because the SPAM dycor diagnoses the dry +// density in a way that preserves the total density to match the reference +// density. However, this is problematic in the presence of the GCM forcing +// which naturally changes the total density. Therefore, we can recall the +// previous dry density and use it to replace the horizontal mean returned from +// the dycor, which ensures that the impact of GCM forcing is preserved while +// also retaining the dry density perturbations created by the dycor. inline void pam_state_recall_dry_density( pam::PamCoupler &coupler ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; @@ -270,29 +277,30 @@ inline void pam_state_recall_dry_density( pam::PamCoupler &coupler ) { auto crm_rho_d = dm_device.get("density_dry"); auto tmp_rho_d = dm_device.get("density_dry_save"); //------------------------------------------------------------------------------------------------ - #ifdef MMF_ALT_DENSITY_RECALL - // initialize horizontal mean of previous CRM dry density - real2d old_hmean_rho_d("old_hmean_rho_d" ,nz,nens); - real2d new_hmean_rho_d("new_hmean_rho_d" ,nz,nens); - parallel_for(SimpleBounds<2>(nz,nens), YAKL_LAMBDA (int k, int iens) { - old_hmean_rho_d(k,iens) = 0; - new_hmean_rho_d(k,iens) = 0; - }); - real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions - // calculate horizontal mean of previous CRM dry density - parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { - atomicAdd( old_hmean_rho_d(k,iens), tmp_rho_d(k,j,i,iens) * r_nx_ny ); - atomicAdd( new_hmean_rho_d(k,iens), crm_rho_d(k,j,i,iens) * r_nx_ny ); - }); - // replace horizontal mean dry density with prevoius value - parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { - crm_rho_d(k,j,i,iens) = crm_rho_d(k,j,i,iens) - new_hmean_rho_d(k,iens) + old_hmean_rho_d(k,iens); - }); - #else - parallel_for( "Recall CRM dry density",SimpleBounds<4>(nz,ny,nx,nens),YAKL_LAMBDA (int k_crm, int j, int i, int iens) { - crm_rho_d(k_crm,j,i,iens) = tmp_rho_d(k_crm,j,i,iens); - }); - #endif + // initialize horizontal mean of previous CRM dry density + real2d old_hmean_rho_d("old_hmean_rho_d" ,nz,nens); + real2d new_hmean_rho_d("new_hmean_rho_d" ,nz,nens); + parallel_for(SimpleBounds<2>(nz,nens), YAKL_LAMBDA (int k, int iens) { + old_hmean_rho_d(k,iens) = 0; + new_hmean_rho_d(k,iens) = 0; + }); + real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions + // calculate horizontal mean of previous CRM dry density + parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { + atomicAdd( old_hmean_rho_d(k,iens), tmp_rho_d(k,j,i,iens) * r_nx_ny ); + atomicAdd( new_hmean_rho_d(k,iens), crm_rho_d(k,j,i,iens) * r_nx_ny ); + }); + // replace horizontal mean dry density with previous value + parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { + crm_rho_d(k,j,i,iens) = crm_rho_d(k,j,i,iens) - new_hmean_rho_d(k,iens) + old_hmean_rho_d(k,iens); + }); + //------------------------------------------------------------------------------------------------ + // Original simple appraoch to completely reinstate the previous dry density field + // (this was shown to negatively impact the model stability) + // This was used initially, but the above + // parallel_for( "Recall CRM dry density",SimpleBounds<4>(nz,ny,nx,nens),YAKL_LAMBDA (int k_crm, int j, int i, int iens) { + // crm_rho_d(k_crm,j,i,iens) = tmp_rho_d(k_crm,j,i,iens); + // }); //------------------------------------------------------------------------------------------------ } From 6b7dcbeae6ac486a59041eb78e882ae06e1a8ea0 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 27 Feb 2024 08:37:57 -0800 Subject: [PATCH 10/45] update PAM submodule --- components/eam/src/physics/crm/pam/external | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/external b/components/eam/src/physics/crm/pam/external index 87731d56aeee..9d692dccdd67 160000 --- a/components/eam/src/physics/crm/pam/external +++ b/components/eam/src/physics/crm/pam/external @@ -1 +1 @@ -Subproject commit 87731d56aeee4bfae4750e17732828ba186367a7 +Subproject commit 9d692dccdd67d9c4dc278c76290504a2af8cd187 From 0af5931b90c565976fdcf94472ab187903ae6d50 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 26 Mar 2024 12:38:59 -0700 Subject: [PATCH 11/45] PAM VT update --- .../eam/src/physics/crm/pam/pam_driver.cpp | 14 ++++++------- .../physics/crm/pam/pam_variance_transport.h | 20 +++++++++++-------- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index 2453f2ad61b2..e6f85ed0d5ac 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -25,7 +25,7 @@ #include "p3_f90.hpp" #include "pam_debug.h" -bool constexpr enable_check_state = true; +bool constexpr enable_check_state = false; extern "C" void pam_driver() { //------------------------------------------------------------------------------------------------ @@ -180,10 +180,11 @@ extern "C" void pam_driver() { // Apply forcing tendencies if (use_MMF_VT) { pam_variance_transport_apply_forcing(coupler); } - coupler.run_module( "apply_gcm_forcing_tendencies" , modules::apply_gcm_forcing_tendencies ); if (enable_check_state) { pam_debug_check_state(coupler, 2, nstep); } - coupler.run_module( "radiation" , [&] (pam::PamCoupler &coupler) {rad .timeStep(coupler);} ); + coupler.run_module( "apply_gcm_forcing_tendencies" , modules::apply_gcm_forcing_tendencies ); if (enable_check_state) { pam_debug_check_state(coupler, 3, nstep); } + coupler.run_module( "radiation" , [&] (pam::PamCoupler &coupler) {rad .timeStep(coupler);} ); + if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } // Dynamics if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } @@ -191,13 +192,12 @@ extern "C" void pam_driver() { coupler.run_module( "dycore", [&] (pam::PamCoupler &coupler) {dycore.timeStep(coupler);} ); if (do_density_save_recall) { pam_state_recall_dry_density(coupler); } if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"dycor"); } - if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } + if (enable_check_state) { pam_debug_check_state(coupler, 5, nstep); } // Sponge layer damping if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } coupler.run_module( "sponge_layer", modules::sponge_layer ); if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sponge"); } - // if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } // Apply hyperdiffusion to account for lack of horizontal mixing in SHOC pam_hyperdiffusion(coupler); @@ -207,13 +207,13 @@ extern "C" void pam_driver() { if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } coupler.run_module( "sgs", [&] (pam::PamCoupler &coupler) {sgs .timeStep(coupler);} ); if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sgs"); } - if (enable_check_state) { pam_debug_check_state(coupler, 5, nstep); } + if (enable_check_state) { pam_debug_check_state(coupler, 6, nstep); } // Microphysics - P3 if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } coupler.run_module( "micro", [&] (pam::PamCoupler &coupler) {micro .timeStep(coupler);} ); if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"micro"); } - if (enable_check_state) { pam_debug_check_state(coupler, 6, nstep); } + if (enable_check_state) { pam_debug_check_state(coupler, 7, nstep); } // CRM mean state acceleration if (use_crm_accel && !coupler.get_option("crm_acceleration_ceaseflag")) { diff --git a/components/eam/src/physics/crm/pam/pam_variance_transport.h b/components/eam/src/physics/crm/pam/pam_variance_transport.h index b7067631203f..1033df4280d4 100644 --- a/components/eam/src/physics/crm/pam/pam_variance_transport.h +++ b/components/eam/src/physics/crm/pam/pam_variance_transport.h @@ -2,6 +2,7 @@ #include "pam_coupler.h" + inline void pam_variance_transport_init( pam::PamCoupler &coupler ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; @@ -24,6 +25,7 @@ inline void pam_variance_transport_init( pam::PamCoupler &coupler ) { //------------------------------------------------------------------------------------------------ } + inline void pam_variance_transport_diagnose( pam::PamCoupler &coupler ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; @@ -83,6 +85,7 @@ inline void pam_variance_transport_diagnose( pam::PamCoupler &coupler ) { //------------------------------------------------------------------------------------------------ } + inline void pam_variance_transport_compute_forcing( pam::PamCoupler &coupler ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; @@ -135,9 +138,10 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { // large-scale forcing from variance transport. This is meant to // protect against creating unstable situations, although // problematic scenarios were extremely rare in testing. - // A scaling limit of +/- 10% was found to be adequate. - real constexpr pert_scale_min = 1.0 - 0.005; - real constexpr pert_scale_max = 1.0 + 0.005; + // A scaling limit of +/- 10% was found to be adequate in SAM, + // but PAM seems much more sensitive (not sure why), so we use 0.1% here + real constexpr pert_scale_min = 1.0 - 0.001; + real constexpr pert_scale_max = 1.0 + 0.001; //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); auto rhov = dm_device.get("water_vapor" ); @@ -162,14 +166,15 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { //------------------------------------------------------------------------------------------------ // calculate scaling factor for local perturbations parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { + real tmp; // initialize scaling factors to 1.0 temp_pert_scale(k,n) = 1.0; rhov_pert_scale(k,n) = 1.0; uvel_pert_scale(k,n) = 1.0; - // set scaling factors as long as there are perturbations to scale - if (vt_temp(k,n)>0.0) { temp_pert_scale(k,n) = sqrt( 1.0 + crm_dt * vt_temp_forcing_tend(k,n) / vt_temp(k,n) ); } - if (vt_rhov(k,n)>0.0) { rhov_pert_scale(k,n) = sqrt( 1.0 + crm_dt * vt_rhov_forcing_tend(k,n) / vt_rhov(k,n) ); } - if (vt_uvel(k,n)>0.0) { uvel_pert_scale(k,n) = sqrt( 1.0 + crm_dt * vt_uvel_forcing_tend(k,n) / vt_uvel(k,n) ); } + // calculate variance scaling factor + tmp = 1.+crm_dt*vt_temp_forcing_tend(k,n)/vt_temp(k,n); if (tmp>0){ temp_pert_scale(k,n) = sqrt(tmp); } + tmp = 1.+crm_dt*vt_rhov_forcing_tend(k,n)/vt_rhov(k,n); if (tmp>0){ rhov_pert_scale(k,n) = sqrt(tmp); } + tmp = 1.+crm_dt*vt_uvel_forcing_tend(k,n)/vt_uvel(k,n); if (tmp>0){ uvel_pert_scale(k,n) = sqrt(tmp); } // enforce minimum scaling temp_pert_scale(k,n) = std::max( temp_pert_scale(k,n), pert_scale_min ); rhov_pert_scale(k,n) = std::max( rhov_pert_scale(k,n), pert_scale_min ); @@ -261,7 +266,6 @@ inline void pam_variance_transport_copy_to_host( pam::PamCoupler &coupler ) { vt_temp_feedback_tend.deep_copy_to(output_vt_temp_tend_host); vt_rhov_feedback_tend.deep_copy_to(output_vt_rhov_tend_host); vt_uvel_feedback_tend.deep_copy_to(output_vt_uvel_tend_host); - yakl::fence(); //------------------------------------------------------------------------------------------------ } From 6fe86f61e08038fb95f964b0cfb5907d5a53423a Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 26 Mar 2024 12:47:40 -0700 Subject: [PATCH 12/45] update PAM submodule --- components/eam/src/physics/crm/pam/external | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/external b/components/eam/src/physics/crm/pam/external index 9d692dccdd67..ff652f3f2a6c 160000 --- a/components/eam/src/physics/crm/pam/external +++ b/components/eam/src/physics/crm/pam/external @@ -1 +1 @@ -Subproject commit 9d692dccdd67d9c4dc278c76290504a2af8cd187 +Subproject commit ff652f3f2a6c2d675cd957883245aa6c036ef110 From 59b9476bbc96ed060ff12e4a7ba3480795b67872 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 26 Mar 2024 12:48:49 -0700 Subject: [PATCH 13/45] update PAM hyperdiffusion timescale default --- components/eam/src/physics/crm/pam/pam_hyperdiffusion.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h index 4452c583afe3..b632b0e49093 100644 --- a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h +++ b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h @@ -28,7 +28,7 @@ inline void pam_hyperdiffusion( pam::PamCoupler &coupler ) { #ifdef MMF_PAM_HDT real constexpr hd_timescale = MMF_PAM_HDT; // damping time scale [sec] #else - real constexpr hd_timescale = 10.0; // damping time scale [sec] + real constexpr hd_timescale = 60.0; // damping time scale [sec] #endif //------------------------------------------------------------------------------------------------ real4d hd_temp("hd_temp",nz,ny,nx,nens); From ebdb83028c0c026a227c038ffdf9279c575bab42 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 29 Mar 2024 14:05:15 -0700 Subject: [PATCH 14/45] revert CRM timestep to 5 sec for MMF2 --- components/eam/cime_config/config_component.xml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/eam/cime_config/config_component.xml b/components/eam/cime_config/config_component.xml index 577c3256db42..784fe44195ed 100755 --- a/components/eam/cime_config/config_component.xml +++ b/components/eam/cime_config/config_component.xml @@ -72,8 +72,8 @@ -crm samxx -crm_dt 10 - -crm pam -pam_dycor spam -crm_dt 8 - -crm pam -pam_dycor awfl -crm_dt 8 + -crm pam -pam_dycor spam -crm_dt 5 + -crm pam -pam_dycor awfl -crm_dt 5 -use_MMF -nlev 60 -crm_nz 50 -crm_dx 2000 -crm_nx 64 -crm_ny 1 -crm_nx_rad 4 -crm_ny_rad 1 -crm_dx 2000 -crm_nx 45 -crm_ny 1 -crm_nx_rad 5 -crm_ny_rad 1 From 0194a2b1ea565e95765fd741d95def019414c3b6 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 5 Apr 2024 08:18:19 -0700 Subject: [PATCH 15/45] pam_debug update --- .../eam/src/physics/crm/pam/pam_debug.h | 41 +++++++++++++------ 1 file changed, 29 insertions(+), 12 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_debug.h b/components/eam/src/physics/crm/pam/pam_debug.h index 937d18cf154e..2d3dc1efb894 100644 --- a/components/eam/src/physics/crm/pam/pam_debug.h +++ b/components/eam/src/physics/crm/pam/pam_debug.h @@ -27,28 +27,32 @@ inline void pam_debug_init( pam::PamCoupler &coupler ) { auto nens = coupler.get_option("ncrms"); //------------------------------------------------------------------------------------------------ dm_device.register_and_allocate("debug_save_temp", "saved temp for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); + dm_device.register_and_allocate("debug_save_rhod", "saved rhod for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_rhov", "saved rhov for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_rhoc", "saved rhoc for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_rhoi", "saved rhoi for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_uvel", "saved uvel for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); dm_device.register_and_allocate("debug_save_wvel", "saved wvel for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} ); auto debug_save_temp = dm_device.get("debug_save_temp"); + auto debug_save_rhod = dm_device.get("debug_save_rhod"); auto debug_save_rhov = dm_device.get("debug_save_rhov"); auto debug_save_rhoc = dm_device.get("debug_save_rhoc"); auto debug_save_rhoi = dm_device.get("debug_save_rhoi"); auto debug_save_uvel = dm_device.get("debug_save_uvel"); auto debug_save_wvel = dm_device.get("debug_save_wvel"); //------------------------------------------------------------------------------------------------ - auto temp = dm_device.get("temp"); - auto rhov = dm_device.get("water_vapor"); - auto rhoc = dm_device.get("cloud_water"); - auto rhoi = dm_device.get("ice"); - auto uvel = dm_device.get("uvel"); - auto wvel = dm_device.get("wvel"); + auto temp = dm_device.get("temp"); + auto rhod = dm_device.get("density_dry"); + auto rhov = dm_device.get("water_vapor"); + auto rhoc = dm_device.get("cloud_water"); + auto rhoi = dm_device.get("ice"); + auto uvel = dm_device.get("uvel"); + auto wvel = dm_device.get("wvel"); //------------------------------------------------------------------------------------------------ parallel_for("copy data to saved debug variables", SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { debug_save_temp(k,j,i,iens) = temp(k,j,i,iens); + debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens); debug_save_rhov(k,j,i,iens) = rhov(k,j,i,iens); debug_save_rhoc(k,j,i,iens) = rhoc(k,j,i,iens); debug_save_rhoi(k,j,i,iens) = rhoi(k,j,i,iens); @@ -69,12 +73,14 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { auto nz = coupler.get_option("crm_nz"); auto nens = coupler.get_option("ncrms"); auto temp = dm_device.get("temp"); + auto rhod = dm_device.get("density_dry"); auto rhov = dm_device.get("water_vapor"); auto rhoc = dm_device.get("cloud_water"); auto rhoi = dm_device.get("ice"); auto uvel = dm_device.get("uvel"); auto wvel = dm_device.get("wvel"); auto debug_save_temp = dm_device.get("debug_save_temp"); + auto debug_save_rhod = dm_device.get("debug_save_rhod"); auto debug_save_rhov = dm_device.get("debug_save_rhov"); auto debug_save_rhoc = dm_device.get("debug_save_rhoc"); auto debug_save_rhoi = dm_device.get("debug_save_rhoi"); @@ -90,18 +96,21 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { auto phis = input_phis(iens)/grav; // Check for NaNs const auto is_nan_t_atm = isnan( temp(k,j,i,iens) ); + const auto is_nan_d_atm = isnan( rhod(k,j,i,iens) ); const auto is_nan_q_atm = isnan( rhov(k,j,i,iens) ); - if ( is_nan_t_atm || is_nan_q_atm ) { + if ( is_nan_t_atm || is_nan_q_atm || is_nan_d_atm ) { auto phis = input_phis(iens)/grav; - printf("PAM-DEBUG nan-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + printf("PAM-DEBUG nan-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, temp(k,j,i,iens), + rhod(k,j,i,iens), rhov(k,j,i,iens), rhoc(k,j,i,iens), rhoi(k,j,i,iens), uvel(k,j,i,iens), wvel(k,j,i,iens), debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), debug_save_rhov(k,j,i,iens), debug_save_rhoc(k,j,i,iens), debug_save_rhoi(k,j,i,iens), @@ -111,18 +120,21 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { } // Check for negative values const auto is_neg_t_atm = temp(k,j,i,iens)<0; + const auto is_neg_d_atm = rhod(k,j,i,iens)<0; const auto is_neg_q_atm = rhov(k,j,i,iens)<0; - if ( is_neg_t_atm || is_neg_q_atm ) { + if ( is_neg_t_atm || is_neg_q_atm || is_neg_d_atm ) { auto phis = input_phis(iens)/grav; - printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, temp(k,j,i,iens), + rhod(k,j,i,iens), rhov(k,j,i,iens), rhoc(k,j,i,iens), rhoi(k,j,i,iens), uvel(k,j,i,iens), wvel(k,j,i,iens), debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), debug_save_rhov(k,j,i,iens), debug_save_rhoc(k,j,i,iens), debug_save_rhoi(k,j,i,iens), @@ -133,15 +145,17 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { // Check for low temperature const auto is_low_t = temp(k,j,i,iens)<100; if ( is_low_t ) { - printf("PAM-DEBUG low-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + printf("PAM-DEBUG low-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, temp(k,j,i,iens), + rhod(k,j,i,iens), rhov(k,j,i,iens), rhoc(k,j,i,iens), rhoi(k,j,i,iens), uvel(k,j,i,iens), wvel(k,j,i,iens), debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), debug_save_rhov(k,j,i,iens), debug_save_rhoc(k,j,i,iens), debug_save_rhoi(k,j,i,iens), @@ -153,15 +167,17 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { const auto is_large_pos_w = wvel(k,j,i,iens)> 40; const auto is_large_neg_w = wvel(k,j,i,iens)<-40; if ( is_large_pos_w || is_large_neg_w ) { - printf("PAM-DEBUG large-W - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + printf("PAM-DEBUG large-W - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, temp(k,j,i,iens), + rhod(k,j,i,iens), rhov(k,j,i,iens), rhoc(k,j,i,iens), rhoi(k,j,i,iens), uvel(k,j,i,iens), wvel(k,j,i,iens), debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), debug_save_rhov(k,j,i,iens), debug_save_rhoc(k,j,i,iens), debug_save_rhoi(k,j,i,iens), @@ -171,6 +187,7 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { } // update saved previous values debug_save_temp(k,j,i,iens) = temp(k,j,i,iens); + debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens); debug_save_rhov(k,j,i,iens) = rhov(k,j,i,iens); debug_save_rhoc(k,j,i,iens) = rhoc(k,j,i,iens); debug_save_rhoi(k,j,i,iens) = rhoi(k,j,i,iens); From 730931f76f3a17a9c2da0bad2407190707f377d3 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 9 Apr 2024 11:17:12 -0700 Subject: [PATCH 16/45] bug fix in pam_state --- .../eam/src/physics/crm/pam/pam_state.h | 34 ++++++++++++++----- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_state.h b/components/eam/src/physics/crm/pam/pam_state.h index 3266899ba217..a34a867eb6e8 100644 --- a/components/eam/src/physics/crm/pam/pam_state.h +++ b/components/eam/src/physics/crm/pam/pam_state.h @@ -144,7 +144,7 @@ inline void pam_state_set_reference_state( pam::PamCoupler &coupler ) { real2d hmean_temp ("hmean_temp" ,nz ,nens); real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions // initialize horizontal means - parallel_for(SimpleBounds<2>(nz,nens), YAKL_LAMBDA (int k, int iens) { + parallel_for(SimpleBounds<2>(nz+1,nens), YAKL_LAMBDA (int k, int iens) { hmean_pint(k,iens) = 0; if (k < nz) { hmean_pmid (k,iens) = 0; @@ -166,24 +166,42 @@ inline void pam_state_set_reference_state( pam::PamCoupler &coupler ) { atomicAdd( hmean_temp (k,iens), temp (k,j,i,iens) * r_nx_ny ); }); // calculate interface pressure from mid-level pressure - parallel_for(SimpleBounds<4>(nz+1,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int iens) { + // parallel_for(SimpleBounds<4>(nz+1,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int iens) { + parallel_for(SimpleBounds<2>(nz+1,nens) , YAKL_LAMBDA (int k, int iens) { if (k == 0 ) { - real rho = rho_d(k ,j,i,iens)+rho_v(k ,j,i,iens); - real dz = zint(k+1,iens)-zint(k ,iens); + real rho = hmean_rho_d(k,iens)+hmean_rho_v(k,iens); + real dz = zint(k+1,iens)-zint(k,iens); hmean_pint(k,iens) = hmean_pmid(k ,iens) + grav*rho*dz/2; } else if (k == nz) { - real rho = rho_d(k-1,j,i,iens)+rho_v(k-1,j,i,iens); - real dz = zint(k ,iens)-zint(k-1,iens); + real rho = hmean_rho_d(k-1,iens)+hmean_rho_v(k-1,iens); + real dz = zint(k,iens)-zint(k-1,iens); hmean_pint(k,iens) = hmean_pmid(k-1,iens) - grav*rho*dz/2; } else { - real rhokm1 = rho_d(k-1,j,i,iens)+rho_v(k-1,j,i,iens); - real rhokm0 = rho_d(k ,j,i,iens)+rho_v(k ,j,i,iens); + real rhokm1 = hmean_rho_d(k-1,iens)+hmean_rho_v(k-1,iens);; + real rhokm0 = hmean_rho_d(k ,iens)+hmean_rho_v(k ,iens); real dzkm1 = zint(k ,iens)-zint(k-1,iens); real dzkm0 = zint(k+1,iens)-zint(k ,iens); hmean_pint(k,iens) = 0.5_fp * ( hmean_pmid(k-1,iens) - grav*rhokm1*dzkm1/2 + hmean_pmid(k ,iens) + grav*rhokm0*dzkm0/2 ); } }); + + auto &dm_host = coupler.get_data_manager_host_readonly(); + auto lat = dm_host.get("latitude" ).createDeviceCopy(); + auto lon = dm_host.get("longitude" ).createDeviceCopy(); + auto input_phis = dm_host.get("input_phis").createDeviceCopy(); + // check that interface pressure is reasonable + parallel_for(SimpleBounds<2>(nz+1,nens) , YAKL_LAMBDA (int k, int iens) { + if ( hmean_pint(k,iens) < hmean_pmid(k,iens) ) { + auto phis = input_phis(iens)/grav; + printf("PAM-DEBUG bad-pint - k:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- pint:%8.2g pmid:%8.2g \n", + k,iens,lat(iens),lon(iens),phis, + hmean_pint(k,iens), + hmean_pmid(k,iens) + ); + } + }); + // set anelastic reference state from CRM horizontal mean parallel_for(SimpleBounds<2>(nz+1,nens) , YAKL_LAMBDA (int k, int iens) { ref_presi(k,iens) = hmean_pint(k,iens); From b094fc0db0534982a3f52aadd2ba5466e09be4e7 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 9 Apr 2024 11:20:48 -0700 Subject: [PATCH 17/45] pam_debug updates --- .../eam/src/physics/crm/pam/pam_debug.h | 73 ++++++++++++------- 1 file changed, 47 insertions(+), 26 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_debug.h b/components/eam/src/physics/crm/pam/pam_debug.h index 2d3dc1efb894..43d2b5bf4f5d 100644 --- a/components/eam/src/physics/crm/pam/pam_debug.h +++ b/components/eam/src/physics/crm/pam/pam_debug.h @@ -142,32 +142,31 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { debug_save_wvel(k,j,i,iens) ); } - // Check for low temperature - const auto is_low_t = temp(k,j,i,iens)<100; - if ( is_low_t ) { - printf("PAM-DEBUG low-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", - nstep,id,k,i,iens,lat(iens),lon(iens),phis, - temp(k,j,i,iens), - rhod(k,j,i,iens), - rhov(k,j,i,iens), - rhoc(k,j,i,iens), - rhoi(k,j,i,iens), - uvel(k,j,i,iens), - wvel(k,j,i,iens), - debug_save_temp(k,j,i,iens), - debug_save_rhod(k,j,i,iens), - debug_save_rhov(k,j,i,iens), - debug_save_rhoc(k,j,i,iens), - debug_save_rhoi(k,j,i,iens), - debug_save_uvel(k,j,i,iens), - debug_save_wvel(k,j,i,iens) - ); - } - // Check for large vertical velocity - const auto is_large_pos_w = wvel(k,j,i,iens)> 40; - const auto is_large_neg_w = wvel(k,j,i,iens)<-40; - if ( is_large_pos_w || is_large_neg_w ) { - printf("PAM-DEBUG large-W - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + // // Check for low temperature + // const auto is_low_t = temp(k,j,i,iens)<100; + // if ( is_low_t ) { + // printf("PAM-DEBUG low-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + // nstep,id,k,i,iens,lat(iens),lon(iens),phis, + // temp(k,j,i,iens), + // rhod(k,j,i,iens), + // rhov(k,j,i,iens), + // rhoc(k,j,i,iens), + // rhoi(k,j,i,iens), + // uvel(k,j,i,iens), + // wvel(k,j,i,iens), + // debug_save_temp(k,j,i,iens), + // debug_save_rhod(k,j,i,iens), + // debug_save_rhov(k,j,i,iens), + // debug_save_rhoc(k,j,i,iens), + // debug_save_rhoi(k,j,i,iens), + // debug_save_uvel(k,j,i,iens), + // debug_save_wvel(k,j,i,iens) + // ); + // } + // Check for large temperature drops + const auto is_drop_t = (temp(k,j,i,iens)-debug_save_temp(k,j,i,iens))<-50; + if ( is_drop_t ) { + printf("PAM-DEBUG drop-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, temp(k,j,i,iens), rhod(k,j,i,iens), @@ -185,6 +184,28 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { debug_save_wvel(k,j,i,iens) ); } + // // Check for large vertical velocity + // const auto is_large_pos_w = wvel(k,j,i,iens)> 40; + // const auto is_large_neg_w = wvel(k,j,i,iens)<-40; + // if ( is_large_pos_w || is_large_neg_w ) { + // printf("PAM-DEBUG large-W - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + // nstep,id,k,i,iens,lat(iens),lon(iens),phis, + // temp(k,j,i,iens), + // rhod(k,j,i,iens), + // rhov(k,j,i,iens), + // rhoc(k,j,i,iens), + // rhoi(k,j,i,iens), + // uvel(k,j,i,iens), + // wvel(k,j,i,iens), + // debug_save_temp(k,j,i,iens), + // debug_save_rhod(k,j,i,iens), + // debug_save_rhov(k,j,i,iens), + // debug_save_rhoc(k,j,i,iens), + // debug_save_rhoi(k,j,i,iens), + // debug_save_uvel(k,j,i,iens), + // debug_save_wvel(k,j,i,iens) + // ); + // } // update saved previous values debug_save_temp(k,j,i,iens) = temp(k,j,i,iens); debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens); From 0720d0dca888c047638492a6c5312fb5932d88aa Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 9 Apr 2024 11:24:16 -0700 Subject: [PATCH 18/45] disable mean-state acceleration for dry density --- .../eam/src/physics/crm/pam/pam_accelerate.h | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_accelerate.h b/components/eam/src/physics/crm/pam/pam_accelerate.h index a48ca628f691..a7b179a85c9d 100644 --- a/components/eam/src/physics/crm/pam/pam_accelerate.h +++ b/components/eam/src/physics/crm/pam/pam_accelerate.h @@ -25,7 +25,7 @@ inline void pam_accelerate_init( pam::PamCoupler &coupler ) { auto crm_accel_uv = coupler.get_option("crm_accel_uv"); //------------------------------------------------------------------------------------------------ dm_device.register_and_allocate("accel_save_t", "saved temperature for MSA", {nz,nens}, {"z","nens"} ); - dm_device.register_and_allocate("accel_save_r", "saved dry density for MSA", {nz,nens}, {"z","nens"} ); + // dm_device.register_and_allocate("accel_save_r", "saved dry density for MSA", {nz,nens}, {"z","nens"} ); dm_device.register_and_allocate("accel_save_q", "saved total water for MSA", {nz,nens}, {"z","nens"} ); dm_device.register_and_allocate("accel_save_u", "saved uvel for MSA", {nz,nens}, {"z","nens"} ); dm_device.register_and_allocate("accel_save_v", "saved vvel for MSA", {nz,nens}, {"z","nens"} ); @@ -45,14 +45,14 @@ inline void pam_accelerate_diagnose( pam::PamCoupler &coupler ) { auto crm_accel_uv = coupler.get_option("crm_accel_uv"); //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); - auto rhod = dm_device.get("density_dry"); + // auto rhod = dm_device.get("density_dry"); auto rhov = dm_device.get("water_vapor"); auto rhol = dm_device.get("cloud_water"); auto rhoi = dm_device.get("ice" ); auto uvel = dm_device.get("uvel" ); auto vvel = dm_device.get("vvel" ); auto accel_save_t = dm_device.get("accel_save_t"); - auto accel_save_r = dm_device.get("accel_save_r"); + // auto accel_save_r = dm_device.get("accel_save_r"); auto accel_save_q = dm_device.get("accel_save_q"); auto accel_save_u = dm_device.get("accel_save_u"); auto accel_save_v = dm_device.get("accel_save_v"); @@ -60,7 +60,7 @@ inline void pam_accelerate_diagnose( pam::PamCoupler &coupler ) { // compute horizontal means needed later for mean-state acceleration parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { accel_save_t(k,n) = 0.0; - accel_save_r(k,n) = 0.0; + // accel_save_r(k,n) = 0.0; accel_save_q(k,n) = 0.0; if (crm_accel_uv) { accel_save_u(k,n) = 0.0; @@ -70,7 +70,7 @@ inline void pam_accelerate_diagnose( pam::PamCoupler &coupler ) { real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { yakl::atomicAdd( accel_save_t(k,n), temp(k,j,i,n) * r_nx_ny ); - yakl::atomicAdd( accel_save_r(k,n), rhod(k,j,i,n) * r_nx_ny ); + // yakl::atomicAdd( accel_save_r(k,n), rhod(k,j,i,n) * r_nx_ny ); yakl::atomicAdd( accel_save_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny ); if (crm_accel_uv) { yakl::atomicAdd( accel_save_u(k,n), uvel(k,j,i,n) * r_nx_ny ); @@ -93,14 +93,14 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { auto nx = coupler.get_option("crm_nx"); //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); - auto rhod = dm_device.get("density_dry"); + // auto rhod = dm_device.get("density_dry"); auto rhov = dm_device.get("water_vapor"); auto rhol = dm_device.get("cloud_water"); auto rhoi = dm_device.get("ice" ); auto uvel = dm_device.get("uvel" ); auto vvel = dm_device.get("vvel" ); auto accel_save_t = dm_device.get("accel_save_t"); - auto accel_save_r = dm_device.get("accel_save_r"); + // auto accel_save_r = dm_device.get("accel_save_r"); auto accel_save_q = dm_device.get("accel_save_q"); auto accel_save_u = dm_device.get("accel_save_u"); auto accel_save_v = dm_device.get("accel_save_v"); @@ -109,12 +109,12 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { real crm_accel_factor = coupler.get_option("crm_accel_factor"); //------------------------------------------------------------------------------------------------ real2d hmean_t ("hmean_t", nz,nens); - real2d hmean_r ("hmean_r", nz,nens); + // real2d hmean_r ("hmean_r", nz,nens); real2d hmean_q ("hmean_q", nz,nens); real2d hmean_u ("hmean_u", nz,nens); real2d hmean_v ("hmean_v", nz,nens); real2d ttend_acc("ttend_acc", nz,nens); - real2d rtend_acc("rtend_acc", nz,nens); + // real2d rtend_acc("rtend_acc", nz,nens); real2d qtend_acc("qtend_acc", nz,nens); real2d utend_acc("utend_acc", nz,nens); real2d vtend_acc("vtend_acc", nz,nens); @@ -127,7 +127,7 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { // Compute the horizontal mean for each variable parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { hmean_t(k,n) = 0.0; - hmean_r(k,n) = 0.0; + // hmean_r(k,n) = 0.0; hmean_q(k,n) = 0.0; if (crm_accel_uv) { hmean_u(k,n) = 0.0; @@ -137,7 +137,7 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { yakl::atomicAdd( hmean_t(k,n), temp(k,j,i,n) * r_nx_ny ); - yakl::atomicAdd( hmean_r(k,n), rhod(k,j,i,n) * r_nx_ny ); + // yakl::atomicAdd( hmean_r(k,n), rhod(k,j,i,n) * r_nx_ny ); yakl::atomicAdd( hmean_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny ); if (crm_accel_uv) { yakl::atomicAdd( hmean_u(k,n), uvel(k,j,i,n) * r_nx_ny ); @@ -150,7 +150,7 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { ScalarLiveOut ceaseflag_liveout(false); parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { ttend_acc(k,n) = hmean_t(k,n) - accel_save_t(k,n); - rtend_acc(k,n) = hmean_r(k,n) - accel_save_r(k,n); + // rtend_acc(k,n) = hmean_r(k,n) - accel_save_r(k,n); qtend_acc(k,n) = hmean_q(k,n) - accel_save_q(k,n); if (crm_accel_uv) { utend_acc(k,n) = hmean_u(k,n) - accel_save_u(k,n); @@ -193,7 +193,7 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { temp(k,j,i,n) = temp(k,j,i,n) + crm_accel_factor * ttend_acc(k,n); - rhod(k,j,i,n) = rhod(k,j,i,n) + crm_accel_factor * rtend_acc(k,n); + // rhod(k,j,i,n) = rhod(k,j,i,n) + crm_accel_factor * rtend_acc(k,n); rhov(k,j,i,n) = rhov(k,j,i,n) + crm_accel_factor * qtend_acc(k,n); if (crm_accel_uv) { uvel(k,j,i,n) = uvel(k,j,i,n) + crm_accel_factor * utend_acc(k,n); From 081b18326623709d8f7019d4fcaedec551c51852 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 9 Apr 2024 11:29:11 -0700 Subject: [PATCH 19/45] update pam_driver --- components/eam/src/physics/crm/pam/pam_driver.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index e6f85ed0d5ac..69e67a90894c 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -25,7 +25,7 @@ #include "p3_f90.hpp" #include "pam_debug.h" -bool constexpr enable_check_state = false; +bool constexpr enable_check_state = true; extern "C" void pam_driver() { //------------------------------------------------------------------------------------------------ @@ -62,7 +62,7 @@ extern "C" void pam_driver() { coupler.set_option("spam_clip_vertical_velocities",true); coupler.set_option("spam_adjust_crm_per_phys_using_vert_cfl",true); coupler.set_option("spam_target_cfl",0.7); - coupler.set_option("spam_max_w",50.0); + coupler.set_option("spam_max_w",30.0); //------------------------------------------------------------------------------------------------ // Allocate the coupler state and retrieve host/device data managers coupler.allocate_coupler_state( crm_nz , crm_ny , crm_nx , nens ); @@ -91,8 +91,10 @@ extern "C" void pam_driver() { // Copy input CRM state (saved by the GCM) to coupler pam_state_copy_input_to_coupler(coupler); - // // update CRM dry density to match GCM and disable dry density forcing - // pam_state_update_dry_density(coupler); + #ifdef MMF_DISABLE_DENSITY_FORCING + // update CRM dry density to match GCM and disable dry density forcing + pam_state_update_dry_density(coupler); + #endif // if debugging - initialize saved state variables and check initial CRM state if (enable_check_state) { From fcb801bf90e7bd379b510e2c2ff3ee75e5992186 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 12:06:32 -0400 Subject: [PATCH 20/45] adjust white space --- components/eam/src/physics/crm/crm_physics.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/eam/src/physics/crm/crm_physics.F90 b/components/eam/src/physics/crm/crm_physics.F90 index 2beeded06bb3..51628d5ea74e 100644 --- a/components/eam/src/physics/crm/crm_physics.F90 +++ b/components/eam/src/physics/crm/crm_physics.F90 @@ -1384,8 +1384,8 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, call pam_mirror_array_readwrite( 'output_ni_mean', crm_output%ni_mean, '' ) call pam_mirror_array_readwrite( 'output_qm_mean', crm_output%qm_mean, '' ) call pam_mirror_array_readwrite( 'output_bm_mean', crm_output%bm_mean, '' ) - call pam_mirror_array_readwrite( 'output_rho_d_mean', crm_output%rho_d_mean, '' ) - call pam_mirror_array_readwrite( 'output_rho_v_mean', crm_output%rho_v_mean, '' ) + call pam_mirror_array_readwrite( 'output_rho_d_mean', crm_output%rho_d_mean, '' ) + call pam_mirror_array_readwrite( 'output_rho_v_mean', crm_output%rho_v_mean, '' ) call pam_mirror_array_readwrite( 'output_qt_ls', crm_output%qt_ls, '' ) call pam_mirror_array_readwrite( 'output_t_ls', crm_output%t_ls, '' ) From 44ba5145df639f554be9e931740848f39a63464f Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 12:13:13 -0400 Subject: [PATCH 21/45] change how PAM reference state pressure is set --- .../eam/src/physics/crm/pam/pam_state.h | 57 +++++++++---------- 1 file changed, 26 insertions(+), 31 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_state.h b/components/eam/src/physics/crm/pam/pam_state.h index a34a867eb6e8..fb6553df07ef 100644 --- a/components/eam/src/physics/crm/pam/pam_state.h +++ b/components/eam/src/physics/crm/pam/pam_state.h @@ -108,6 +108,7 @@ inline void pam_state_set_reference_state( pam::PamCoupler &coupler ) { auto &dm_device = coupler.get_data_manager_device_readwrite(); auto nens = coupler.get_option("ncrms"); auto nz = coupler.get_option("crm_nz"); + auto gcm_nlev = coupler.get_option("gcm_nlev"); auto ref_presi = dm_device.get("ref_presi"); auto ref_pres = dm_device.get("ref_pres"); auto ref_rho_d = dm_device.get("ref_density_dry"); @@ -131,6 +132,14 @@ inline void pam_state_set_reference_state( pam::PamCoupler &coupler ) { auto rho_v = dm_device.get("water_vapor"); auto rho_c = dm_device.get("cloud_water"); auto rho_i = dm_device.get("ice"); + + auto &dm_host = coupler.get_data_manager_host_readonly(); + auto input_pmid = dm_host.get("input_pmid").createDeviceCopy(); + auto input_pint = dm_host.get("input_pint").createDeviceCopy(); + + auto lat = dm_host.get("latitude" ).createDeviceCopy(); + auto lon = dm_host.get("longitude" ).createDeviceCopy(); + auto input_phis = dm_host.get("input_phis").createDeviceCopy(); //------------------------------------------------------------------------------------------------ // Set anelastic reference state with current CRM mean state @@ -157,48 +166,34 @@ inline void pam_state_set_reference_state( pam::PamCoupler &coupler ) { }); // calculate horizontal means parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { - real pmid_tmp = rho_d(k,j,i,iens)*R_d*temp(k,j,i,iens) + rho_v(k,j,i,iens)*R_v*temp(k,j,i,iens); - atomicAdd( hmean_pmid (k,iens), pmid_tmp * r_nx_ny ); atomicAdd( hmean_rho_d(k,iens), rho_d (k,j,i,iens) * r_nx_ny ); atomicAdd( hmean_rho_v(k,iens), rho_v (k,j,i,iens) * r_nx_ny ); atomicAdd( hmean_rho_c(k,iens), rho_c (k,j,i,iens) * r_nx_ny ); atomicAdd( hmean_rho_i(k,iens), rho_i (k,j,i,iens) * r_nx_ny ); atomicAdd( hmean_temp (k,iens), temp (k,j,i,iens) * r_nx_ny ); }); - // calculate interface pressure from mid-level pressure - // parallel_for(SimpleBounds<4>(nz+1,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int iens) { - parallel_for(SimpleBounds<2>(nz+1,nens) , YAKL_LAMBDA (int k, int iens) { - if (k == 0 ) { - real rho = hmean_rho_d(k,iens)+hmean_rho_v(k,iens); - real dz = zint(k+1,iens)-zint(k,iens); - hmean_pint(k,iens) = hmean_pmid(k ,iens) + grav*rho*dz/2; - } else if (k == nz) { - real rho = hmean_rho_d(k-1,iens)+hmean_rho_v(k-1,iens); - real dz = zint(k,iens)-zint(k-1,iens); - hmean_pint(k,iens) = hmean_pmid(k-1,iens) - grav*rho*dz/2; - } else { - real rhokm1 = hmean_rho_d(k-1,iens)+hmean_rho_v(k-1,iens);; - real rhokm0 = hmean_rho_d(k ,iens)+hmean_rho_v(k ,iens); - real dzkm1 = zint(k ,iens)-zint(k-1,iens); - real dzkm0 = zint(k+1,iens)-zint(k ,iens); - hmean_pint(k,iens) = 0.5_fp * ( hmean_pmid(k-1,iens) - grav*rhokm1*dzkm1/2 + - hmean_pmid(k ,iens) + grav*rhokm0*dzkm0/2 ); - } + + // Use GCM state for reference pressure - previously, the current CRM pressure + // was used, but the way the interface pressure was calculated led to problems + // in edge cases (ex. over topography) - switching to GCM pressure works well + parallel_for( SimpleBounds<2>(nz+1,nens) , YAKL_LAMBDA (int k_crm, int iens) { + int k_gcm = (gcm_nlev+1)-1-k_crm; + hmean_pint(k_crm,iens) = input_pint(k_gcm,iens); + }); + parallel_for(SimpleBounds<2>(nz,nens), YAKL_LAMBDA (int k_crm, int iens) { + int k_gcm = gcm_nlev-1-k_crm; + hmean_pmid(k_crm,iens) = input_pmid(k_gcm,iens); }); - auto &dm_host = coupler.get_data_manager_host_readonly(); - auto lat = dm_host.get("latitude" ).createDeviceCopy(); - auto lon = dm_host.get("longitude" ).createDeviceCopy(); - auto input_phis = dm_host.get("input_phis").createDeviceCopy(); // check that interface pressure is reasonable - parallel_for(SimpleBounds<2>(nz+1,nens) , YAKL_LAMBDA (int k, int iens) { + parallel_for(SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int iens) { if ( hmean_pint(k,iens) < hmean_pmid(k,iens) ) { auto phis = input_phis(iens)/grav; - printf("PAM-DEBUG bad-pint - k:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- pint:%8.2g pmid:%8.2g \n", - k,iens,lat(iens),lon(iens),phis, - hmean_pint(k,iens), - hmean_pmid(k,iens) - ); + printf("PAM-STATE - bad interface pressure for reference state - "+\ + "k:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- "+\ + "pint:%12.4f pmid:%12.4f \n", + k,iens,lat(iens),lon(iens),phis, + hmean_pint(k,iens),hmean_pmid(k,iens)); } }); From 8c8cb1c9e358119154f9ee4f974fc2d219a15c3b Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 12:14:45 -0400 Subject: [PATCH 22/45] update PAM variance transport --- .../physics/crm/pam/pam_variance_transport.h | 39 ++++++++++--------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_variance_transport.h b/components/eam/src/physics/crm/pam/pam_variance_transport.h index 1033df4280d4..ae9f0747ee31 100644 --- a/components/eam/src/physics/crm/pam/pam_variance_transport.h +++ b/components/eam/src/physics/crm/pam/pam_variance_transport.h @@ -38,6 +38,8 @@ inline void pam_variance_transport_diagnose( pam::PamCoupler &coupler ) { //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); auto rhov = dm_device.get("water_vapor"); + auto rhoc = dm_device.get("cloud_water"); + auto rhoi = dm_device.get("ice" ); auto uvel = dm_device.get("uvel" ); auto vt_temp = dm_device.get("vt_temp" ); auto vt_rhov = dm_device.get("vt_rhov" ); @@ -64,15 +66,17 @@ inline void pam_variance_transport_diagnose( pam::PamCoupler &coupler ) { // calculate horizontal mean real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int n) { + real rhot = rhov(k,j,i,n) + rhoc(k,j,i,n) + rhoi(k,j,i,n); yakl::atomicAdd( temp_mean(k,n), temp(k,j,i,n)*r_nx_ny ); - yakl::atomicAdd( rhov_mean(k,n), rhov(k,j,i,n)*r_nx_ny ); + yakl::atomicAdd( rhov_mean(k,n), rhot *r_nx_ny ); yakl::atomicAdd( uvel_mean(k,n), uvel(k,j,i,n)*r_nx_ny ); }); //------------------------------------------------------------------------------------------------ // calculate fluctuations from horz mean parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int n) { + real rhot = rhov(k,j,i,n) + rhoc(k,j,i,n) + rhoi(k,j,i,n); vt_temp_pert(k,j,i,n) = temp(k,j,i,n) - temp_mean(k,n); - vt_rhov_pert(k,j,i,n) = rhov(k,j,i,n) - rhov_mean(k,n); + vt_rhov_pert(k,j,i,n) = rhot - rhov_mean(k,n); vt_uvel_pert(k,j,i,n) = uvel(k,j,i,n) - uvel_mean(k,n); }); //------------------------------------------------------------------------------------------------ @@ -101,9 +105,6 @@ inline void pam_variance_transport_compute_forcing( pam::PamCoupler &coupler ) { // update CRM variance values pam_variance_transport_diagnose(coupler); //------------------------------------------------------------------------------------------------ - auto temp = dm_device.get("temp" ); - auto rhov = dm_device.get("water_vapor" ); - auto uvel = dm_device.get("uvel" ); auto vt_temp = dm_device.get("vt_temp" ); auto vt_rhov = dm_device.get("vt_rhov" ); auto vt_uvel = dm_device.get("vt_uvel" ); @@ -134,14 +135,16 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { auto nx = coupler.get_option("crm_nx"); auto crm_dt = coupler.get_option("crm_dt"); //------------------------------------------------------------------------------------------------ + // update CRM variance values + pam_variance_transport_diagnose(coupler); + //------------------------------------------------------------------------------------------------ // min and max perturbation scaling values are used to limit the // large-scale forcing from variance transport. This is meant to // protect against creating unstable situations, although // problematic scenarios were extremely rare in testing. - // A scaling limit of +/- 10% was found to be adequate in SAM, - // but PAM seems much more sensitive (not sure why), so we use 0.1% here - real constexpr pert_scale_min = 1.0 - 0.001; - real constexpr pert_scale_max = 1.0 + 0.001; + // A scaling limit of +/- 10% was found to be adequate. + real constexpr pert_scale_min = 1.0 - 0.1; + real constexpr pert_scale_max = 1.0 + 0.1; //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); auto rhov = dm_device.get("water_vapor" ); @@ -161,9 +164,6 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { real2d rhov_pert_scale("rhov_pert_scale", nz, nens); real2d uvel_pert_scale("uvel_pert_scale", nz, nens); //------------------------------------------------------------------------------------------------ - // update CRM variance values - pam_variance_transport_diagnose(coupler); - //------------------------------------------------------------------------------------------------ // calculate scaling factor for local perturbations parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { real tmp; @@ -172,9 +172,15 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { rhov_pert_scale(k,n) = 1.0; uvel_pert_scale(k,n) = 1.0; // calculate variance scaling factor - tmp = 1.+crm_dt*vt_temp_forcing_tend(k,n)/vt_temp(k,n); if (tmp>0){ temp_pert_scale(k,n) = sqrt(tmp); } - tmp = 1.+crm_dt*vt_rhov_forcing_tend(k,n)/vt_rhov(k,n); if (tmp>0){ rhov_pert_scale(k,n) = sqrt(tmp); } - tmp = 1.+crm_dt*vt_uvel_forcing_tend(k,n)/vt_uvel(k,n); if (tmp>0){ uvel_pert_scale(k,n) = sqrt(tmp); } + real tmp_t_scale = -1.0; + real tmp_q_scale = -1.0; + real tmp_u_scale = -1.0; + if (vt_temp(k,n)>0.0) { tmp_t_scale = 1. + crm_dt*vt_temp_forcing_tend(k,n) / vt_temp(k,n); } + if (vt_rhov(k,n)>0.0) { tmp_q_scale = 1. + crm_dt*vt_rhov_forcing_tend(k,n) / vt_rhov(k,n); } + if (vt_uvel(k,n)>0.0) { tmp_u_scale = 1. + crm_dt*vt_uvel_forcing_tend(k,n) / vt_uvel(k,n); } + if (tmp>0.0){ temp_pert_scale(k,n) = sqrt(tmp_t_scale); } + if (tmp>0.0){ rhov_pert_scale(k,n) = sqrt(tmp_q_scale); } + if (tmp>0.0){ uvel_pert_scale(k,n) = sqrt(tmp_u_scale); } // enforce minimum scaling temp_pert_scale(k,n) = std::max( temp_pert_scale(k,n), pert_scale_min ); rhov_pert_scale(k,n) = std::max( rhov_pert_scale(k,n), pert_scale_min ); @@ -213,9 +219,6 @@ inline void pam_variance_transport_compute_feedback( pam::PamCoupler &coupler ) // update CRM variance values pam_variance_transport_diagnose(coupler); //------------------------------------------------------------------------------------------------ - auto temp = dm_device.get("temp" ); - auto rhov = dm_device.get("water_vapor"); - auto uvel = dm_device.get("uvel" ); auto vt_temp = dm_device.get("vt_temp" ); auto vt_rhov = dm_device.get("vt_rhov" ); auto vt_uvel = dm_device.get("vt_uvel" ); From 98726502a74a6067dd2f4e96c04db227e5f3a790 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 12:16:11 -0400 Subject: [PATCH 23/45] update MMF2 default config --- components/eam/cime_config/config_component.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/components/eam/cime_config/config_component.xml b/components/eam/cime_config/config_component.xml index 784fe44195ed..056b607953f9 100755 --- a/components/eam/cime_config/config_component.xml +++ b/components/eam/cime_config/config_component.xml @@ -72,15 +72,15 @@ -crm samxx -crm_dt 10 - -crm pam -pam_dycor spam -crm_dt 5 - -crm pam -pam_dycor awfl -crm_dt 5 + -crm pam -pam_dycor spam -crm_dt 10 + -crm pam -pam_dycor awfl -crm_dt 10 -use_MMF -nlev 60 -crm_nz 50 -crm_dx 2000 -crm_nx 64 -crm_ny 1 -crm_nx_rad 4 -crm_ny_rad 1 -crm_dx 2000 -crm_nx 45 -crm_ny 1 -crm_nx_rad 5 -crm_ny_rad 1 -MMF_microphysics_scheme sam1mom -chem none -MMF_microphysics_scheme p3 -chem none -rad rrtmgp -rrtmgpxx - -use_MMF_VT + -use_MMF_VT -use_MMF_ESMT -aquaplanet -aquaplanet -rce From 2380d339ca04b8f6c878d2fcfdf76ecaa408663f Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 12:29:53 -0400 Subject: [PATCH 24/45] bug fix in PAM variance transport --- .../eam/src/physics/crm/pam/pam_variance_transport.h | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_variance_transport.h b/components/eam/src/physics/crm/pam/pam_variance_transport.h index ae9f0747ee31..e258b3b7f3ad 100644 --- a/components/eam/src/physics/crm/pam/pam_variance_transport.h +++ b/components/eam/src/physics/crm/pam/pam_variance_transport.h @@ -143,8 +143,8 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { // protect against creating unstable situations, although // problematic scenarios were extremely rare in testing. // A scaling limit of +/- 10% was found to be adequate. - real constexpr pert_scale_min = 1.0 - 0.1; - real constexpr pert_scale_max = 1.0 + 0.1; + real constexpr pert_scale_min = 1.0 - 0.05; + real constexpr pert_scale_max = 1.0 + 0.05; //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); auto rhov = dm_device.get("water_vapor" ); @@ -166,7 +166,6 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { //------------------------------------------------------------------------------------------------ // calculate scaling factor for local perturbations parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { - real tmp; // initialize scaling factors to 1.0 temp_pert_scale(k,n) = 1.0; rhov_pert_scale(k,n) = 1.0; @@ -178,9 +177,9 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { if (vt_temp(k,n)>0.0) { tmp_t_scale = 1. + crm_dt*vt_temp_forcing_tend(k,n) / vt_temp(k,n); } if (vt_rhov(k,n)>0.0) { tmp_q_scale = 1. + crm_dt*vt_rhov_forcing_tend(k,n) / vt_rhov(k,n); } if (vt_uvel(k,n)>0.0) { tmp_u_scale = 1. + crm_dt*vt_uvel_forcing_tend(k,n) / vt_uvel(k,n); } - if (tmp>0.0){ temp_pert_scale(k,n) = sqrt(tmp_t_scale); } - if (tmp>0.0){ rhov_pert_scale(k,n) = sqrt(tmp_q_scale); } - if (tmp>0.0){ uvel_pert_scale(k,n) = sqrt(tmp_u_scale); } + if (tmp_t_scale>0.0){ temp_pert_scale(k,n) = sqrt(tmp_t_scale); } + if (tmp_q_scale>0.0){ rhov_pert_scale(k,n) = sqrt(tmp_q_scale); } + if (tmp_u_scale>0.0){ uvel_pert_scale(k,n) = sqrt(tmp_u_scale); } // enforce minimum scaling temp_pert_scale(k,n) = std::max( temp_pert_scale(k,n), pert_scale_min ); rhov_pert_scale(k,n) = std::max( rhov_pert_scale(k,n), pert_scale_min ); From 63df74c0bb17552d825f55cb3b04edc5f5b83014 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 12:30:50 -0400 Subject: [PATCH 25/45] implement PAM driver subcycling --- .../eam/src/physics/crm/pam/pam_driver.cpp | 173 ++++++++++++------ 1 file changed, 122 insertions(+), 51 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index 69e67a90894c..0833d2beff8d 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -1,5 +1,4 @@ #include "pam_coupler.h" -// #include "params.h" #include "Dycore.h" #include "Microphysics.h" #include "SGS.h" @@ -14,26 +13,90 @@ #include "pam_output.h" #include "pam_accelerate.h" #include "pam_variance_transport.h" +#include "pam_hyperdiffusion.h" #include "sponge_layer.h" #include "surface_friction.h" #include "scream_cxx_interface_finalize.h" -#include "pam_hyperdiffusion.h" - // Needed for p3_init #include "p3_functions.hpp" #include "p3_f90.hpp" #include "pam_debug.h" -bool constexpr enable_check_state = true; +bool constexpr enable_check_state = false; -extern "C" void pam_driver() { + +inline int pam_driver_set_subcycle_timestep( pam::PamCoupler &coupler, real crm_dt_fixed ) { + // calculate the CFL condition and adjust the PAM time loop subcylcing //------------------------------------------------------------------------------------------------ - using yakl::intrinsics::abs; - using yakl::intrinsics::maxval; - using yakl::atomicAdd; using yakl::c::parallel_for; using yakl::c::SimpleBounds; + using yakl::atomicMax; + //------------------------------------------------------------------------------------------------ + auto nens = coupler.get_option("ncrms"); + auto gcm_nlev = coupler.get_option("gcm_nlev"); + auto crm_nz = coupler.get_option("crm_nz"); + auto crm_nx = coupler.get_option("crm_nx"); + auto crm_ny = coupler.get_option("crm_ny"); + auto crm_dx = coupler.get_option("crm_dx"); + auto crm_dy = coupler.get_option("crm_dy"); + auto &dm_device = coupler.get_data_manager_device_readonly(); + auto &dm_host = coupler.get_data_manager_host_readonly(); + auto uvel = dm_device.get("uvel"); + auto wvel = dm_device.get("wvel"); + auto input_zint = dm_host.get("input_zint").createDeviceCopy(); + //------------------------------------------------------------------------------------------------ + yakl::ParallelMax pmax( crm_nz*nens ); + real cfl = 0; + int num_subcycle = 1; + int constexpr max_num_subcycle = 10; + real2d wvel_max("wvel_max",crm_nz,nens); + real2d uvel_max("uvel_max",crm_nz,nens); + real2d cfl_max("cfl_max", crm_nz,nens); + //------------------------------------------------------------------------------------------------ + // initialize max U and W arrays + parallel_for( SimpleBounds<2>(crm_nz,nens) , YAKL_LAMBDA (int k, int n) { + wvel_max(k,n) = 0.0; + uvel_max(k,n) = 0.0; + }); + // calculate max U and W + parallel_for( SimpleBounds<4>(crm_nz,crm_ny,crm_nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { + yakl::atomicMax(uvel_max(k,n), sqrt(uvel(k,j,i,n)*uvel(k,j,i,n)) ); + yakl::atomicMax(wvel_max(k,n), fabs(wvel(k,j,i,n)) ); + }); + // find max CFL between horizontal and vertical CFL values + parallel_for( SimpleBounds<2>(crm_nz,nens) , YAKL_LAMBDA (int k, int n) { + int k_gcm = gcm_nlev-1-k; + real crm_dz = input_zint(k_gcm,n) - input_zint(k_gcm+1,n); + real cfl_u = uvel_max(k,n)*crm_dt_fixed/crm_dx; + real cfl_w = wvel_max(k,n)*crm_dt_fixed/crm_dz; + cfl_max(k,n) = max(cfl_u,cfl_w); + }); + // calculate final CFL across ensemble + real cfl_loc = pmax(cfl_max.data()); + cfl = max(cfl,cfl_loc); + // update number of subcycles and time step + num_subcycle = max(num_subcycle,max(1,static_cast(ceil(cfl/0.7)))); + real crm_dt_subcycle = crm_dt_fixed / num_subcycle; + coupler.set_option("crm_dt",crm_dt_subcycle); + // check for excessive subcylcing - don't exit, just print + if(num_subcycle > max_num_subcycle) { + real umax = pmax(uvel_max.data()); + real wmax = pmax(wvel_max.data()); + printf("PAM_DRIVER - WARNING: excessive subcycling!"+\ + " - num_subcycle: %3.3d dt: %8.4f cfl: %8.4f umax: %8.2f wmax: %8.2f \n", + num_subcycle,crm_dt_subcycle,cfl,umax,wmax); + // exit(-1); + } + + return num_subcycle; + //------------------------------------------------------------------------------------------------ +} + + +extern "C" void pam_driver() { + // This is the primary method for running the PAM CRM in E3SM-MMF. + //------------------------------------------------------------------------------------------------ auto &coupler = pam_interface::get_coupler(); //------------------------------------------------------------------------------------------------ // retreive coupler options @@ -42,8 +105,10 @@ extern "C" void pam_driver() { auto crm_nz = coupler.get_option("crm_nz"); auto crm_nx = coupler.get_option("crm_nx"); auto crm_ny = coupler.get_option("crm_ny"); + // auto crm_dx = coupler.get_option("crm_dx"); + // auto crm_dy = coupler.get_option("crm_dy"); auto gcm_dt = coupler.get_option("gcm_dt"); - auto crm_dt = coupler.get_option("crm_dt"); + auto crm_dt_fixed = coupler.get_option("crm_dt"); auto is_first_step = coupler.get_option("is_first_step"); auto is_restart = coupler.get_option("is_restart"); bool use_crm_accel = coupler.get_option("use_crm_accel"); @@ -162,7 +227,7 @@ extern "C" void pam_driver() { //------------------------------------------------------------------------------------------------ // set number of CRM steps - int nstop = int(gcm_dt/crm_dt); + int nstop = int(gcm_dt/crm_dt_fixed); // for mean-state acceleration adjust nstop and diagnose horizontal means if (use_crm_accel) { @@ -171,51 +236,58 @@ extern "C" void pam_driver() { }; // Run the CRM - real etime_crm = 0; int nstep = 0; - // while (etime_crm < gcm_dt) { while (nstep < nstop) { - if (crm_dt == 0.) { crm_dt = dycore.compute_time_step(coupler); } - if (etime_crm + crm_dt > gcm_dt) { crm_dt = gcm_dt - etime_crm; } + //-------------------------------------------------------------------------- + //-------------------------------------------------------------------------- + auto num_subcycle = pam_driver_set_subcycle_timestep(coupler,crm_dt_fixed); + #if defined(MMF_PAM_DYCOR_SPAM) + dycore.update_dt(coupler); + #endif if (enable_check_state) { pam_debug_check_state(coupler, 1, nstep); } - // Apply forcing tendencies - if (use_MMF_VT) { pam_variance_transport_apply_forcing(coupler); } - if (enable_check_state) { pam_debug_check_state(coupler, 2, nstep); } - coupler.run_module( "apply_gcm_forcing_tendencies" , modules::apply_gcm_forcing_tendencies ); - if (enable_check_state) { pam_debug_check_state(coupler, 3, nstep); } - coupler.run_module( "radiation" , [&] (pam::PamCoupler &coupler) {rad .timeStep(coupler);} ); - if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } - - // Dynamics - if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } - if (do_density_save_recall) { pam_state_save_dry_density(coupler); } - coupler.run_module( "dycore", [&] (pam::PamCoupler &coupler) {dycore.timeStep(coupler);} ); - if (do_density_save_recall) { pam_state_recall_dry_density(coupler); } - if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"dycor"); } - if (enable_check_state) { pam_debug_check_state(coupler, 5, nstep); } - - // Sponge layer damping - if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } - coupler.run_module( "sponge_layer", modules::sponge_layer ); - if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sponge"); } - - // Apply hyperdiffusion to account for lack of horizontal mixing in SHOC - pam_hyperdiffusion(coupler); - - // Turbulence - SHOC - coupler.run_module( "compute_surface_friction", modules::compute_surface_friction ); - if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } - coupler.run_module( "sgs", [&] (pam::PamCoupler &coupler) {sgs .timeStep(coupler);} ); - if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sgs"); } - if (enable_check_state) { pam_debug_check_state(coupler, 6, nstep); } - - // Microphysics - P3 - if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } - coupler.run_module( "micro", [&] (pam::PamCoupler &coupler) {micro .timeStep(coupler);} ); - if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"micro"); } - if (enable_check_state) { pam_debug_check_state(coupler, 7, nstep); } + // loop for adaptive subcyling based on CFL + for(int icycle=1; icycle<=num_subcycle; icycle++) { + + // Apply forcing tendencies + if (use_MMF_VT) { pam_variance_transport_apply_forcing(coupler); } + if (enable_check_state) { pam_debug_check_state(coupler, 2, nstep); } + coupler.run_module( "apply_gcm_forcing_tendencies" , modules::apply_gcm_forcing_tendencies ); + if (enable_check_state) { pam_debug_check_state(coupler, 3, nstep); } + coupler.run_module( "radiation" , [&] (pam::PamCoupler &coupler) {rad .timeStep(coupler);} ); + if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); } + + // Dynamics + if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } + if (do_density_save_recall) { pam_state_save_dry_density(coupler); } + coupler.run_module( "dycore", [&] (pam::PamCoupler &coupler) {dycore.timeStep(coupler);} ); + if (do_density_save_recall) { pam_state_recall_dry_density(coupler); } + if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"dycor"); } + if (enable_check_state) { pam_debug_check_state(coupler, 5, nstep); } + + // Sponge layer damping + if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } + coupler.run_module( "sponge_layer", modules::sponge_layer ); + if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sponge"); } + + // Apply hyperdiffusion to account for lack of horizontal mixing in SHOC + pam_hyperdiffusion(coupler); + + // Turbulence - SHOC + coupler.run_module( "compute_surface_friction", modules::compute_surface_friction ); + if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } + coupler.run_module( "sgs", [&] (pam::PamCoupler &coupler) {sgs .timeStep(coupler);} ); + if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sgs"); } + if (enable_check_state) { pam_debug_check_state(coupler, 6, nstep); } + + // Microphysics - P3 + if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); } + coupler.run_module( "micro", [&] (pam::PamCoupler &coupler) {micro .timeStep(coupler);} ); + if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"micro"); } + if (enable_check_state) { pam_debug_check_state(coupler, 7, nstep); } + + } // num_subcycle // CRM mean state acceleration if (use_crm_accel && !coupler.get_option("crm_acceleration_ceaseflag")) { @@ -227,7 +299,6 @@ extern "C" void pam_driver() { pam_radiation_timestep_aggregation(coupler); pam_statistics_timestep_aggregation(coupler); - etime_crm += crm_dt; nstep += 1; } //------------------------------------------------------------------------------------------------ From 7cc042a07306b2370931de60e03cc077f70f8454 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 13:47:26 -0400 Subject: [PATCH 26/45] bug fix --- components/eam/src/physics/crm/pam/pam_driver.cpp | 2 +- components/eam/src/physics/crm/pam/pam_state.h | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index 0833d2beff8d..a05b5cf8b037 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -83,7 +83,7 @@ inline int pam_driver_set_subcycle_timestep( pam::PamCoupler &coupler, real crm_ if(num_subcycle > max_num_subcycle) { real umax = pmax(uvel_max.data()); real wmax = pmax(wvel_max.data()); - printf("PAM_DRIVER - WARNING: excessive subcycling!"+\ + printf("PAM_DRIVER - WARNING: excessive subcycling!" " - num_subcycle: %3.3d dt: %8.4f cfl: %8.4f umax: %8.2f wmax: %8.2f \n", num_subcycle,crm_dt_subcycle,cfl,umax,wmax); // exit(-1); diff --git a/components/eam/src/physics/crm/pam/pam_state.h b/components/eam/src/physics/crm/pam/pam_state.h index fb6553df07ef..270fef7f6c23 100644 --- a/components/eam/src/physics/crm/pam/pam_state.h +++ b/components/eam/src/physics/crm/pam/pam_state.h @@ -189,8 +189,8 @@ inline void pam_state_set_reference_state( pam::PamCoupler &coupler ) { parallel_for(SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int iens) { if ( hmean_pint(k,iens) < hmean_pmid(k,iens) ) { auto phis = input_phis(iens)/grav; - printf("PAM-STATE - bad interface pressure for reference state - "+\ - "k:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- "+\ + printf("PAM-STATE - bad interface pressure for reference state - " + "k:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- " "pint:%12.4f pmid:%12.4f \n", k,iens,lat(iens),lon(iens),phis, hmean_pint(k,iens),hmean_pmid(k,iens)); From f714cf4e2d6612210692eecdf5b2e4eda3a82933 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 13:48:07 -0400 Subject: [PATCH 27/45] bug fix in PAM variance transport --- .../physics/crm/pam/pam_variance_transport.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_variance_transport.h b/components/eam/src/physics/crm/pam/pam_variance_transport.h index e258b3b7f3ad..1eb7779f3469 100644 --- a/components/eam/src/physics/crm/pam/pam_variance_transport.h +++ b/components/eam/src/physics/crm/pam/pam_variance_transport.h @@ -118,9 +118,9 @@ inline void pam_variance_transport_compute_forcing( pam::PamCoupler &coupler ) { // calculate variance transport forcing parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k_crm, int n) { int k_gcm = gcm_nlev-1-k_crm; - vt_temp_forcing_tend(k_crm,n) = ( gcm_vt_temp(k_gcm,n) - vt_temp(k_crm,n) ) * gcm_dt ; - vt_rhov_forcing_tend(k_crm,n) = ( gcm_vt_rhov(k_gcm,n) - vt_rhov(k_crm,n) ) * gcm_dt ; - vt_uvel_forcing_tend(k_crm,n) = ( gcm_vt_uvel(k_gcm,n) - vt_uvel(k_crm,n) ) * gcm_dt ; + vt_temp_forcing_tend(k_crm,n) = ( gcm_vt_temp(k_gcm,n) - vt_temp(k_crm,n) ) / gcm_dt ; + vt_rhov_forcing_tend(k_crm,n) = ( gcm_vt_rhov(k_gcm,n) - vt_rhov(k_crm,n) ) / gcm_dt ; + vt_uvel_forcing_tend(k_crm,n) = ( gcm_vt_uvel(k_gcm,n) - vt_uvel(k_crm,n) ) / gcm_dt ; }); //------------------------------------------------------------------------------------------------ } @@ -174,9 +174,9 @@ inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { real tmp_t_scale = -1.0; real tmp_q_scale = -1.0; real tmp_u_scale = -1.0; - if (vt_temp(k,n)>0.0) { tmp_t_scale = 1. + crm_dt*vt_temp_forcing_tend(k,n) / vt_temp(k,n); } - if (vt_rhov(k,n)>0.0) { tmp_q_scale = 1. + crm_dt*vt_rhov_forcing_tend(k,n) / vt_rhov(k,n); } - if (vt_uvel(k,n)>0.0) { tmp_u_scale = 1. + crm_dt*vt_uvel_forcing_tend(k,n) / vt_uvel(k,n); } + if (vt_temp(k,n)>0.0) { tmp_t_scale = 1. + crm_dt * vt_temp_forcing_tend(k,n) / vt_temp(k,n); } + if (vt_rhov(k,n)>0.0) { tmp_q_scale = 1. + crm_dt * vt_rhov_forcing_tend(k,n) / vt_rhov(k,n); } + if (vt_uvel(k,n)>0.0) { tmp_u_scale = 1. + crm_dt * vt_uvel_forcing_tend(k,n) / vt_uvel(k,n); } if (tmp_t_scale>0.0){ temp_pert_scale(k,n) = sqrt(tmp_t_scale); } if (tmp_q_scale>0.0){ rhov_pert_scale(k,n) = sqrt(tmp_q_scale); } if (tmp_u_scale>0.0){ uvel_pert_scale(k,n) = sqrt(tmp_u_scale); } @@ -236,9 +236,9 @@ inline void pam_variance_transport_compute_feedback( pam::PamCoupler &coupler ) parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k_crm, int n) { int k_gcm = gcm_nlev-1-k_crm; if (k_crm Date: Fri, 12 Apr 2024 13:48:58 -0400 Subject: [PATCH 28/45] update default MMF_PAM_dyn_per_phys --- components/eam/bld/namelist_files/namelist_defaults_eam.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index 74f2572ea040..d7e9afecf20e 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -871,7 +871,7 @@ 0 .false. .true. - 2 + 1 90.0 0.0 From 98a6e4025e9617bfdd0ba27699fa69bbbd959772 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 15:22:03 -0400 Subject: [PATCH 29/45] update default PAM hyperdiffusion timescale --- components/eam/src/physics/crm/pam/pam_hyperdiffusion.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h index b632b0e49093..e75d80265445 100644 --- a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h +++ b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h @@ -28,7 +28,7 @@ inline void pam_hyperdiffusion( pam::PamCoupler &coupler ) { #ifdef MMF_PAM_HDT real constexpr hd_timescale = MMF_PAM_HDT; // damping time scale [sec] #else - real constexpr hd_timescale = 60.0; // damping time scale [sec] + real constexpr hd_timescale = 120.0; // damping time scale [sec] #endif //------------------------------------------------------------------------------------------------ real4d hd_temp("hd_temp",nz,ny,nx,nens); From 2c54c4d9aa7db1659374433517d7a2636357c690 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Fri, 12 Apr 2024 15:22:46 -0400 Subject: [PATCH 30/45] update PAM default crm_dx to 3km --- components/eam/cime_config/config_component.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/cime_config/config_component.xml b/components/eam/cime_config/config_component.xml index 056b607953f9..ef9433ef9664 100755 --- a/components/eam/cime_config/config_component.xml +++ b/components/eam/cime_config/config_component.xml @@ -76,7 +76,7 @@ -crm pam -pam_dycor awfl -crm_dt 10 -use_MMF -nlev 60 -crm_nz 50 -crm_dx 2000 -crm_nx 64 -crm_ny 1 -crm_nx_rad 4 -crm_ny_rad 1 - -crm_dx 2000 -crm_nx 45 -crm_ny 1 -crm_nx_rad 5 -crm_ny_rad 1 + -crm_dx 3000 -crm_nx 45 -crm_ny 1 -crm_nx_rad 5 -crm_ny_rad 1 -MMF_microphysics_scheme sam1mom -chem none -MMF_microphysics_scheme p3 -chem none -rad rrtmgp -rrtmgpxx From 748ad0a808b404ebb3cb43b5032c7f663c3cf177 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Wed, 17 Apr 2024 18:58:56 -0400 Subject: [PATCH 31/45] reduce hd_timescale back to 60 sec fr PAM --- components/eam/src/physics/crm/pam/pam_hyperdiffusion.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h index e75d80265445..b632b0e49093 100644 --- a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h +++ b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h @@ -28,7 +28,7 @@ inline void pam_hyperdiffusion( pam::PamCoupler &coupler ) { #ifdef MMF_PAM_HDT real constexpr hd_timescale = MMF_PAM_HDT; // damping time scale [sec] #else - real constexpr hd_timescale = 120.0; // damping time scale [sec] + real constexpr hd_timescale = 60.0; // damping time scale [sec] #endif //------------------------------------------------------------------------------------------------ real4d hd_temp("hd_temp",nz,ny,nx,nens); From 6b92cabe43740e3a36df9c7aa718d9a5a4bc0301 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Wed, 17 Apr 2024 18:59:54 -0400 Subject: [PATCH 32/45] update PAM submodule --- components/eam/src/physics/crm/pam/external | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/external b/components/eam/src/physics/crm/pam/external index ff652f3f2a6c..8f9538db7138 160000 --- a/components/eam/src/physics/crm/pam/external +++ b/components/eam/src/physics/crm/pam/external @@ -1 +1 @@ -Subproject commit ff652f3f2a6c2d675cd957883245aa6c036ef110 +Subproject commit 8f9538db7138baf26fa2b0f163fae3050d3fe29d From f750a1d9caca3e853a215fbf012888d07bba3176 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 18 Apr 2024 16:23:30 -0400 Subject: [PATCH 33/45] make PAM dry density acceleration optional --- .../eam/src/physics/crm/pam/pam_accelerate.h | 89 +++++++++---------- 1 file changed, 41 insertions(+), 48 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_accelerate.h b/components/eam/src/physics/crm/pam/pam_accelerate.h index a7b179a85c9d..5b277fdca966 100644 --- a/components/eam/src/physics/crm/pam/pam_accelerate.h +++ b/components/eam/src/physics/crm/pam/pam_accelerate.h @@ -17,16 +17,22 @@ inline void pam_accelerate_init( pam::PamCoupler &coupler ) { using yakl::c::parallel_for; using yakl::c::SimpleBounds; using yakl::atomicAdd; + //------------------------------------------------------------------------------------------------ + // The ability to accelerate dry density is implemented, but tests suggested + // that this was contributing to instabilities in PAM. Combined with the + // uncertainty of how to handle dry density in the anelastic case, we will + // disable dry density mean-state acceleration by default. + coupler.set_option("crm_accel_rd",false); + //------------------------------------------------------------------------------------------------ auto &dm_device = coupler.get_data_manager_device_readwrite(); auto nens = coupler.get_option("ncrms"); auto nz = coupler.get_option("crm_nz"); auto ny = coupler.get_option("crm_ny"); auto nx = coupler.get_option("crm_nx"); - auto crm_accel_uv = coupler.get_option("crm_accel_uv"); //------------------------------------------------------------------------------------------------ dm_device.register_and_allocate("accel_save_t", "saved temperature for MSA", {nz,nens}, {"z","nens"} ); - // dm_device.register_and_allocate("accel_save_r", "saved dry density for MSA", {nz,nens}, {"z","nens"} ); dm_device.register_and_allocate("accel_save_q", "saved total water for MSA", {nz,nens}, {"z","nens"} ); + dm_device.register_and_allocate("accel_save_r", "saved dry density for MSA", {nz,nens}, {"z","nens"} ); dm_device.register_and_allocate("accel_save_u", "saved uvel for MSA", {nz,nens}, {"z","nens"} ); dm_device.register_and_allocate("accel_save_v", "saved vvel for MSA", {nz,nens}, {"z","nens"} ); //------------------------------------------------------------------------------------------------ @@ -42,40 +48,37 @@ inline void pam_accelerate_diagnose( pam::PamCoupler &coupler ) { auto nz = coupler.get_option("crm_nz"); auto ny = coupler.get_option("crm_ny"); auto nx = coupler.get_option("crm_nx"); + auto crm_accel_rd = coupler.get_option("crm_accel_rd"); auto crm_accel_uv = coupler.get_option("crm_accel_uv"); //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); - // auto rhod = dm_device.get("density_dry"); auto rhov = dm_device.get("water_vapor"); auto rhol = dm_device.get("cloud_water"); auto rhoi = dm_device.get("ice" ); + auto rhod = dm_device.get("density_dry"); auto uvel = dm_device.get("uvel" ); auto vvel = dm_device.get("vvel" ); auto accel_save_t = dm_device.get("accel_save_t"); - // auto accel_save_r = dm_device.get("accel_save_r"); auto accel_save_q = dm_device.get("accel_save_q"); + auto accel_save_r = dm_device.get("accel_save_r"); auto accel_save_u = dm_device.get("accel_save_u"); auto accel_save_v = dm_device.get("accel_save_v"); //------------------------------------------------------------------------------------------------ // compute horizontal means needed later for mean-state acceleration parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { accel_save_t(k,n) = 0.0; - // accel_save_r(k,n) = 0.0; accel_save_q(k,n) = 0.0; - if (crm_accel_uv) { - accel_save_u(k,n) = 0.0; - accel_save_v(k,n) = 0.0; - } + if (crm_accel_rd) { accel_save_r(k,n) = 0.0; } + if (crm_accel_uv) { accel_save_u(k,n) = 0.0; } + if (crm_accel_uv) { accel_save_v(k,n) = 0.0; } }); real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { yakl::atomicAdd( accel_save_t(k,n), temp(k,j,i,n) * r_nx_ny ); - // yakl::atomicAdd( accel_save_r(k,n), rhod(k,j,i,n) * r_nx_ny ); yakl::atomicAdd( accel_save_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny ); - if (crm_accel_uv) { - yakl::atomicAdd( accel_save_u(k,n), uvel(k,j,i,n) * r_nx_ny ); - yakl::atomicAdd( accel_save_v(k,n), vvel(k,j,i,n) * r_nx_ny ); - } + if (crm_accel_rd) { yakl::atomicAdd( accel_save_r(k,n), rhod(k,j,i,n) * r_nx_ny ); } + if (crm_accel_uv) { yakl::atomicAdd( accel_save_u(k,n), uvel(k,j,i,n) * r_nx_ny ); } + if (crm_accel_uv) { yakl::atomicAdd( accel_save_v(k,n), vvel(k,j,i,n) * r_nx_ny ); } }); //------------------------------------------------------------------------------------------------ } @@ -86,35 +89,35 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { using yakl::c::SimpleBounds; using yakl::atomicAdd; using yakl::ScalarLiveOut; - auto &dm_device = coupler.get_data_manager_device_readwrite(); - auto nens = coupler.get_option("ncrms"); - auto nz = coupler.get_option("crm_nz"); - auto ny = coupler.get_option("crm_ny"); - auto nx = coupler.get_option("crm_nx"); + auto &dm_device = coupler.get_data_manager_device_readwrite(); + auto nens = coupler.get_option("ncrms"); + auto nz = coupler.get_option("crm_nz"); + auto ny = coupler.get_option("crm_ny"); + auto nx = coupler.get_option("crm_nx"); + auto crm_accel_rd = coupler.get_option("crm_accel_rd"); + auto crm_accel_uv = coupler.get_option("crm_accel_uv"); + real crm_accel_factor = coupler.get_option("crm_accel_factor"); //------------------------------------------------------------------------------------------------ auto temp = dm_device.get("temp" ); - // auto rhod = dm_device.get("density_dry"); auto rhov = dm_device.get("water_vapor"); auto rhol = dm_device.get("cloud_water"); auto rhoi = dm_device.get("ice" ); + auto rhod = dm_device.get("density_dry"); auto uvel = dm_device.get("uvel" ); auto vvel = dm_device.get("vvel" ); auto accel_save_t = dm_device.get("accel_save_t"); - // auto accel_save_r = dm_device.get("accel_save_r"); auto accel_save_q = dm_device.get("accel_save_q"); + auto accel_save_r = dm_device.get("accel_save_r"); auto accel_save_u = dm_device.get("accel_save_u"); auto accel_save_v = dm_device.get("accel_save_v"); //------------------------------------------------------------------------------------------------ - bool crm_accel_uv = coupler.get_option("crm_accel_uv"); - real crm_accel_factor = coupler.get_option("crm_accel_factor"); - //------------------------------------------------------------------------------------------------ real2d hmean_t ("hmean_t", nz,nens); - // real2d hmean_r ("hmean_r", nz,nens); + real2d hmean_r ("hmean_r", nz,nens); real2d hmean_q ("hmean_q", nz,nens); real2d hmean_u ("hmean_u", nz,nens); real2d hmean_v ("hmean_v", nz,nens); real2d ttend_acc("ttend_acc", nz,nens); - // real2d rtend_acc("rtend_acc", nz,nens); + real2d rtend_acc("rtend_acc", nz,nens); real2d qtend_acc("qtend_acc", nz,nens); real2d utend_acc("utend_acc", nz,nens); real2d vtend_acc("vtend_acc", nz,nens); @@ -127,22 +130,18 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { // Compute the horizontal mean for each variable parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { hmean_t(k,n) = 0.0; - // hmean_r(k,n) = 0.0; hmean_q(k,n) = 0.0; - if (crm_accel_uv) { - hmean_u(k,n) = 0.0; - hmean_v(k,n) = 0.0; - } + if (crm_accel_rd) { hmean_r(k,n) = 0.0; } + if (crm_accel_uv) { hmean_u(k,n) = 0.0; } + if (crm_accel_uv) { hmean_v(k,n) = 0.0; } }); real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { yakl::atomicAdd( hmean_t(k,n), temp(k,j,i,n) * r_nx_ny ); - // yakl::atomicAdd( hmean_r(k,n), rhod(k,j,i,n) * r_nx_ny ); yakl::atomicAdd( hmean_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny ); - if (crm_accel_uv) { - yakl::atomicAdd( hmean_u(k,n), uvel(k,j,i,n) * r_nx_ny ); - yakl::atomicAdd( hmean_v(k,n), vvel(k,j,i,n) * r_nx_ny ); - } + if (crm_accel_rd) { yakl::atomicAdd( hmean_r(k,n), rhod(k,j,i,n) * r_nx_ny ); } + if (crm_accel_uv) { yakl::atomicAdd( hmean_u(k,n), uvel(k,j,i,n) * r_nx_ny ); } + if (crm_accel_uv) { yakl::atomicAdd( hmean_v(k,n), vvel(k,j,i,n) * r_nx_ny ); } }); //------------------------------------------------------------------------------------------------ // Compute the accelerated tendencies @@ -150,12 +149,10 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { ScalarLiveOut ceaseflag_liveout(false); parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { ttend_acc(k,n) = hmean_t(k,n) - accel_save_t(k,n); - // rtend_acc(k,n) = hmean_r(k,n) - accel_save_r(k,n); qtend_acc(k,n) = hmean_q(k,n) - accel_save_q(k,n); - if (crm_accel_uv) { - utend_acc(k,n) = hmean_u(k,n) - accel_save_u(k,n); - vtend_acc(k,n) = hmean_v(k,n) - accel_save_v(k,n); - } + if (crm_accel_rd) { rtend_acc(k,n) = hmean_r(k,n) - accel_save_r(k,n); } + if (crm_accel_uv) { utend_acc(k,n) = hmean_u(k,n) - accel_save_u(k,n); } + if (crm_accel_uv) { vtend_acc(k,n) = hmean_v(k,n) - accel_save_v(k,n); } if (abs(ttend_acc(k,n)) > dtemp_max) { ceaseflag_liveout = true; } @@ -191,14 +188,11 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { //------------------------------------------------------------------------------------------------ // Apply the accelerated tendencies parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { - temp(k,j,i,n) = temp(k,j,i,n) + crm_accel_factor * ttend_acc(k,n); - // rhod(k,j,i,n) = rhod(k,j,i,n) + crm_accel_factor * rtend_acc(k,n); rhov(k,j,i,n) = rhov(k,j,i,n) + crm_accel_factor * qtend_acc(k,n); - if (crm_accel_uv) { - uvel(k,j,i,n) = uvel(k,j,i,n) + crm_accel_factor * utend_acc(k,n); - vvel(k,j,i,n) = vvel(k,j,i,n) + crm_accel_factor * vtend_acc(k,n); - } + if (crm_accel_rd) { rhod(k,j,i,n) = rhod(k,j,i,n) + crm_accel_factor * rtend_acc(k,n); } + if (crm_accel_uv) { uvel(k,j,i,n) = uvel(k,j,i,n) + crm_accel_factor * utend_acc(k,n); } + if (crm_accel_uv) { vvel(k,j,i,n) = vvel(k,j,i,n) + crm_accel_factor * vtend_acc(k,n); } // apply limiter on temperature temp(k,j,i,n) = std::max( temp_min, temp(k,j,i,n) ); }); @@ -231,4 +225,3 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) { }); //------------------------------------------------------------------------------------------------ } - From 8e65b050739206ac2490374af86e0721562194a8 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 18 Apr 2024 16:24:38 -0400 Subject: [PATCH 34/45] remove MMF_PAM_dyn_per_phys namelist option --- components/eam/bld/build-namelist | 3 --- .../eam/bld/namelist_files/namelist_defaults_eam.xml | 1 - components/eam/bld/namelist_files/namelist_definition.xml | 6 ------ components/eam/src/physics/cam/phys_control.F90 | 8 ++------ components/eam/src/physics/crm/crm_physics.F90 | 7 ------- components/eam/src/physics/crm/pam/pam_driver.cpp | 1 + 6 files changed, 3 insertions(+), 23 deletions(-) diff --git a/components/eam/bld/build-namelist b/components/eam/bld/build-namelist index ac9001c3f7b1..8dc532b6a3df 100755 --- a/components/eam/bld/build-namelist +++ b/components/eam/bld/build-namelist @@ -3941,9 +3941,6 @@ if ($cfg->get('use_MMF')) { # MMF CRM domain orientation add_default($nl, 'MMF_orientation_angle'); - - # PAM dynamics sub-stepping parameter - add_default($nl, 'MMF_PAM_dyn_per_phys'); } add_default($nl, 'do_aerocom_ind3'); diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index d7e9afecf20e..fd7a91f19a1b 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -871,7 +871,6 @@ 0 .false. .true. - 1 90.0 0.0 diff --git a/components/eam/bld/namelist_files/namelist_definition.xml b/components/eam/bld/namelist_files/namelist_definition.xml index cc59697d5bbf..19e356ca6ccb 100644 --- a/components/eam/bld/namelist_files/namelist_definition.xml +++ b/components/eam/bld/namelist_files/namelist_definition.xml @@ -4824,12 +4824,6 @@ Defaults: 3D CRM: 0 - -Number of PAM dycor calls per physics time step. -Default: 2 - - Wavenumber cutoff for filtered MMF variance transport. diff --git a/components/eam/src/physics/cam/phys_control.F90 b/components/eam/src/physics/cam/phys_control.F90 index b87ea442c4e7..b7c9b37fa817 100644 --- a/components/eam/src/physics/cam/phys_control.F90 +++ b/components/eam/src/physics/cam/phys_control.F90 @@ -72,7 +72,6 @@ module phys_control logical :: use_crm_accel = .false. ! true => use MMF CRM mean-state acceleration (MSA) real(r8) :: crm_accel_factor = 2.D0 ! CRM acceleration factor logical :: crm_accel_uv = .true. ! true => apply MMF CRM MSA to momentum fields -integer :: MMF_PAM_dyn_per_phys = 2 ! PAM CRM dynamics steps per CRM physics steps logical :: use_subcol_microp = .false. ! if .true. then use sub-columns in microphysics @@ -229,7 +228,7 @@ subroutine phys_ctl_readnl(nlfile) eddy_scheme, microp_scheme, macrop_scheme, radiation_scheme, srf_flux_avg, & MMF_microphysics_scheme, MMF_orientation_angle, use_MMF, use_ECPP, & use_MMF_VT, MMF_VT_wn_max, use_MMF_ESMT, & - use_crm_accel, crm_accel_factor, crm_accel_uv, MMF_PAM_dyn_per_phys, & + use_crm_accel, crm_accel_factor, crm_accel_uv, & use_subcol_microp, atm_dep_flux, history_amwg, history_verbose, history_vdiag, & get_presc_aero_data,history_aerosol, history_aero_optics, & is_output_interactive_volc, & @@ -314,7 +313,6 @@ subroutine phys_ctl_readnl(nlfile) call mpibcast(use_crm_accel, 1 , mpilog, 0, mpicom) call mpibcast(crm_accel_factor, 1 , mpir8, 0, mpicom) call mpibcast(crm_accel_uv, 1 , mpilog, 0, mpicom) - call mpibcast(MMF_PAM_dyn_per_phys, 1 , mpiint, 0, mpicom) call mpibcast(use_subcol_microp, 1 , mpilog, 0, mpicom) call mpibcast(atm_dep_flux, 1 , mpilog, 0, mpicom) call mpibcast(history_amwg, 1 , mpilog, 0, mpicom) @@ -598,7 +596,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, & prog_modal_aero_out, macrop_scheme_out, ideal_phys_option_out, & use_MMF_out, use_ECPP_out, MMF_microphysics_scheme_out, & MMF_orientation_angle_out, use_MMF_VT_out, MMF_VT_wn_max_out, use_MMF_ESMT_out, & - use_crm_accel_out, crm_accel_factor_out, crm_accel_uv_out, MMF_PAM_dyn_per_phys_out, & + use_crm_accel_out, crm_accel_factor_out, crm_accel_uv_out, & do_clubb_sgs_out, do_shoc_sgs_out, do_tms_out, state_debug_checks_out, & linearize_pbl_winds_out, & do_aerocom_ind3_out, & @@ -641,7 +639,6 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, & logical, intent(out), optional :: use_crm_accel_out real(r8), intent(out), optional :: crm_accel_factor_out logical, intent(out), optional :: crm_accel_uv_out - integer, intent(out), optional :: MMF_PAM_dyn_per_phys_out logical, intent(out), optional :: use_subcol_microp_out logical, intent(out), optional :: atm_dep_flux_out logical, intent(out), optional :: history_amwg_out @@ -745,7 +742,6 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, & if ( present(use_crm_accel_out ) ) use_crm_accel_out = use_crm_accel if ( present(crm_accel_factor_out ) ) crm_accel_factor_out = crm_accel_factor if ( present(crm_accel_uv_out ) ) crm_accel_uv_out = crm_accel_uv - if ( present(MMF_PAM_dyn_per_phys_out) ) MMF_PAM_dyn_per_phys_out = MMF_PAM_dyn_per_phys if ( present(use_subcol_microp_out ) ) use_subcol_microp_out = use_subcol_microp if ( present(macrop_scheme_out ) ) macrop_scheme_out = macrop_scheme diff --git a/components/eam/src/physics/crm/crm_physics.F90 b/components/eam/src/physics/crm/crm_physics.F90 index 51628d5ea74e..39ed6a0876c7 100644 --- a/components/eam/src/physics/crm/crm_physics.F90 +++ b/components/eam/src/physics/crm/crm_physics.F90 @@ -620,8 +620,6 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, logical :: use_MMF_VT_tmp ! flag for MMF variance transport (for Fortran CRM) logical :: use_MMF_ESMT_tmp ! flag for MMF scalar momentum transport (for Fortran CRM) integer :: MMF_VT_wn_max ! wavenumber cutoff for filtered variance transport - - integer :: MMF_PAM_dyn_per_phys ! PAM CRM dynamics steps per CRM physics steps real(r8) :: tmp_e_sat ! temporary saturation vapor pressure real(r8) :: tmp_q_sat ! temporary saturation specific humidity @@ -734,9 +732,6 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, use_crm_accel = use_crm_accel_tmp crm_accel_uv = crm_accel_uv_tmp - ! PAM namelist options - call phys_getopts(MMF_PAM_dyn_per_phys_out = MMF_PAM_dyn_per_phys) - nstep = get_nstep() itim = pbuf_old_tim_idx() ! "Old" pbuf time index (what does all this mean?) @@ -1466,8 +1461,6 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, call pam_set_option('enable_physics_tend_stats', .false. ) - call pam_set_option('crm_dyn_per_phys', MMF_PAM_dyn_per_phys ) - call pam_set_option('is_first_step', (nstep<=1) ) call pam_set_option('is_restart', (nsrest>0) ) call pam_set_option('am_i_root', masterproc ) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index a05b5cf8b037..f8f575ec8a1d 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -122,6 +122,7 @@ extern "C" void pam_driver() { coupler.set_option("crm_acceleration_ceaseflag",false); //------------------------------------------------------------------------------------------------ // coupler options for SPAM dycor + coupler.set_option("crm_dyn_per_phys",1); coupler.set_option("spam_couple_wind_exact_inverse",true); coupler.set_option("spam_clip_negative_densities",true); coupler.set_option("spam_clip_vertical_velocities",true); From 743b8ff2bdcc8a50672e5135ea3fa5c2e6a4ea4d Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 18 Apr 2024 16:26:14 -0400 Subject: [PATCH 35/45] fix white space --- components/eam/bld/namelist_files/namelist_defaults_eam.xml | 1 + 1 file changed, 1 insertion(+) diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index fd7a91f19a1b..56c816b1753c 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -871,6 +871,7 @@ 0 .false. .true. + 90.0 0.0 From 56c11f70c5137f887f204d00cd6bd7de9ed9e29a Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 18 Apr 2024 16:37:56 -0400 Subject: [PATCH 36/45] pam_debug.h clean-up --- .../eam/src/physics/crm/pam/pam_debug.h | 235 ++++++++++-------- 1 file changed, 129 insertions(+), 106 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_debug.h b/components/eam/src/physics/crm/pam/pam_debug.h index 43d2b5bf4f5d..7e54f437f691 100644 --- a/components/eam/src/physics/crm/pam/pam_debug.h +++ b/components/eam/src/physics/crm/pam/pam_debug.h @@ -91,121 +91,144 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { auto input_phis = dm_host.get("input_phis").createDeviceCopy(); real grav = 9.80616; //------------------------------------------------------------------------------------------------ + // logical switches to set which checks are active + bool chk_nan = true; // check for NaN values in select fields + bool chk_neg = true; // check for negative values in select fields + bool chk_low_t = false; // check for low temperatures + bool chk_drop_t = true; // check for large drops in temperature + bool chk_wvel = false; // check for large vertical velocity + //------------------------------------------------------------------------------------------------ // Check for invalid values parallel_for("", SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) { auto phis = input_phis(iens)/grav; + //-------------------------------------------------------------------------- // Check for NaNs - const auto is_nan_t_atm = isnan( temp(k,j,i,iens) ); - const auto is_nan_d_atm = isnan( rhod(k,j,i,iens) ); - const auto is_nan_q_atm = isnan( rhov(k,j,i,iens) ); - if ( is_nan_t_atm || is_nan_q_atm || is_nan_d_atm ) { - auto phis = input_phis(iens)/grav; - printf("PAM-DEBUG nan-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", - nstep,id,k,i,iens,lat(iens),lon(iens),phis, - temp(k,j,i,iens), - rhod(k,j,i,iens), - rhov(k,j,i,iens), - rhoc(k,j,i,iens), - rhoi(k,j,i,iens), - uvel(k,j,i,iens), - wvel(k,j,i,iens), - debug_save_temp(k,j,i,iens), - debug_save_rhod(k,j,i,iens), - debug_save_rhov(k,j,i,iens), - debug_save_rhoc(k,j,i,iens), - debug_save_rhoi(k,j,i,iens), - debug_save_uvel(k,j,i,iens), - debug_save_wvel(k,j,i,iens) - ); + if (chk_nan) { + const auto is_nan_t_atm = isnan( temp(k,j,i,iens) ); + const auto is_nan_d_atm = isnan( rhod(k,j,i,iens) ); + const auto is_nan_q_atm = isnan( rhov(k,j,i,iens) ); + if ( is_nan_t_atm || is_nan_q_atm || is_nan_d_atm ) { + auto phis = input_phis(iens)/grav; + printf("PAM-DEBUG nan-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + nstep,id,k,i,iens,lat(iens),lon(iens),phis, + temp(k,j,i,iens), + rhod(k,j,i,iens), + rhov(k,j,i,iens), + rhoc(k,j,i,iens), + rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), + debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), + debug_save_rhov(k,j,i,iens), + debug_save_rhoc(k,j,i,iens), + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) + ); + } } + //-------------------------------------------------------------------------- // Check for negative values - const auto is_neg_t_atm = temp(k,j,i,iens)<0; - const auto is_neg_d_atm = rhod(k,j,i,iens)<0; - const auto is_neg_q_atm = rhov(k,j,i,iens)<0; - if ( is_neg_t_atm || is_neg_q_atm || is_neg_d_atm ) { - auto phis = input_phis(iens)/grav; - printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", - nstep,id,k,i,iens,lat(iens),lon(iens),phis, - temp(k,j,i,iens), - rhod(k,j,i,iens), - rhov(k,j,i,iens), - rhoc(k,j,i,iens), - rhoi(k,j,i,iens), - uvel(k,j,i,iens), - wvel(k,j,i,iens), - debug_save_temp(k,j,i,iens), - debug_save_rhod(k,j,i,iens), - debug_save_rhov(k,j,i,iens), - debug_save_rhoc(k,j,i,iens), - debug_save_rhoi(k,j,i,iens), - debug_save_uvel(k,j,i,iens), - debug_save_wvel(k,j,i,iens) - ); + if (chk_neg) { + const auto is_neg_t_atm = temp(k,j,i,iens)<0; + const auto is_neg_d_atm = rhod(k,j,i,iens)<0; + const auto is_neg_q_atm = rhov(k,j,i,iens)<0; + if ( is_neg_t_atm || is_neg_q_atm || is_neg_d_atm ) { + auto phis = input_phis(iens)/grav; + printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + nstep,id,k,i,iens,lat(iens),lon(iens),phis, + temp(k,j,i,iens), + rhod(k,j,i,iens), + rhov(k,j,i,iens), + rhoc(k,j,i,iens), + rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), + debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), + debug_save_rhov(k,j,i,iens), + debug_save_rhoc(k,j,i,iens), + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) + ); + } + } + //-------------------------------------------------------------------------- + // Check for low temperature + if (chk_low_t) { + const auto is_low_t = temp(k,j,i,iens)<100; + if ( is_low_t ) { + printf("PAM-DEBUG low-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + nstep,id,k,i,iens,lat(iens),lon(iens),phis, + temp(k,j,i,iens), + rhod(k,j,i,iens), + rhov(k,j,i,iens), + rhoc(k,j,i,iens), + rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), + debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), + debug_save_rhov(k,j,i,iens), + debug_save_rhoc(k,j,i,iens), + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) + ); + } } - // // Check for low temperature - // const auto is_low_t = temp(k,j,i,iens)<100; - // if ( is_low_t ) { - // printf("PAM-DEBUG low-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", - // nstep,id,k,i,iens,lat(iens),lon(iens),phis, - // temp(k,j,i,iens), - // rhod(k,j,i,iens), - // rhov(k,j,i,iens), - // rhoc(k,j,i,iens), - // rhoi(k,j,i,iens), - // uvel(k,j,i,iens), - // wvel(k,j,i,iens), - // debug_save_temp(k,j,i,iens), - // debug_save_rhod(k,j,i,iens), - // debug_save_rhov(k,j,i,iens), - // debug_save_rhoc(k,j,i,iens), - // debug_save_rhoi(k,j,i,iens), - // debug_save_uvel(k,j,i,iens), - // debug_save_wvel(k,j,i,iens) - // ); - // } + //-------------------------------------------------------------------------- // Check for large temperature drops - const auto is_drop_t = (temp(k,j,i,iens)-debug_save_temp(k,j,i,iens))<-50; - if ( is_drop_t ) { - printf("PAM-DEBUG drop-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", - nstep,id,k,i,iens,lat(iens),lon(iens),phis, - temp(k,j,i,iens), - rhod(k,j,i,iens), - rhov(k,j,i,iens), - rhoc(k,j,i,iens), - rhoi(k,j,i,iens), - uvel(k,j,i,iens), - wvel(k,j,i,iens), - debug_save_temp(k,j,i,iens), - debug_save_rhod(k,j,i,iens), - debug_save_rhov(k,j,i,iens), - debug_save_rhoc(k,j,i,iens), - debug_save_rhoi(k,j,i,iens), - debug_save_uvel(k,j,i,iens), - debug_save_wvel(k,j,i,iens) - ); + if (chk_drop_t) { + const auto is_drop_t = (temp(k,j,i,iens)-debug_save_temp(k,j,i,iens))<-50; + if ( is_drop_t ) { + printf("PAM-DEBUG drop-T - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + nstep,id,k,i,iens,lat(iens),lon(iens),phis, + temp(k,j,i,iens), + rhod(k,j,i,iens), + rhov(k,j,i,iens), + rhoc(k,j,i,iens), + rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), + debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), + debug_save_rhov(k,j,i,iens), + debug_save_rhoc(k,j,i,iens), + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) + ); + } + } + //-------------------------------------------------------------------------- + // Check for large vertical velocity + if (chk_wvel) { + const auto is_large_pos_w = wvel(k,j,i,iens)> 40; + const auto is_large_neg_w = wvel(k,j,i,iens)<-40; + if ( is_large_pos_w || is_large_neg_w ) { + printf("PAM-DEBUG large-W - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", + nstep,id,k,i,iens,lat(iens),lon(iens),phis, + temp(k,j,i,iens), + rhod(k,j,i,iens), + rhov(k,j,i,iens), + rhoc(k,j,i,iens), + rhoi(k,j,i,iens), + uvel(k,j,i,iens), + wvel(k,j,i,iens), + debug_save_temp(k,j,i,iens), + debug_save_rhod(k,j,i,iens), + debug_save_rhov(k,j,i,iens), + debug_save_rhoc(k,j,i,iens), + debug_save_rhoi(k,j,i,iens), + debug_save_uvel(k,j,i,iens), + debug_save_wvel(k,j,i,iens) + ); + } } - // // Check for large vertical velocity - // const auto is_large_pos_w = wvel(k,j,i,iens)> 40; - // const auto is_large_neg_w = wvel(k,j,i,iens)<-40; - // if ( is_large_pos_w || is_large_neg_w ) { - // printf("PAM-DEBUG large-W - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", - // nstep,id,k,i,iens,lat(iens),lon(iens),phis, - // temp(k,j,i,iens), - // rhod(k,j,i,iens), - // rhov(k,j,i,iens), - // rhoc(k,j,i,iens), - // rhoi(k,j,i,iens), - // uvel(k,j,i,iens), - // wvel(k,j,i,iens), - // debug_save_temp(k,j,i,iens), - // debug_save_rhod(k,j,i,iens), - // debug_save_rhov(k,j,i,iens), - // debug_save_rhoc(k,j,i,iens), - // debug_save_rhoi(k,j,i,iens), - // debug_save_uvel(k,j,i,iens), - // debug_save_wvel(k,j,i,iens) - // ); - // } + //-------------------------------------------------------------------------- // update saved previous values debug_save_temp(k,j,i,iens) = temp(k,j,i,iens); debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens); From 8eae43e1a3df9f9ddde9a52d37ecad06ec77b199 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 25 Apr 2024 15:13:25 -0400 Subject: [PATCH 37/45] update PAM submodule --- components/eam/src/physics/crm/pam/external | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/external b/components/eam/src/physics/crm/pam/external index 8f9538db7138..5931666c95a3 160000 --- a/components/eam/src/physics/crm/pam/external +++ b/components/eam/src/physics/crm/pam/external @@ -1 +1 @@ -Subproject commit 8f9538db7138baf26fa2b0f163fae3050d3fe29d +Subproject commit 5931666c95a3b75ae63fbe5bff6e743e74cb2c0a From e6e4fb7bb4d1bccfd65e096c31edc3a4b28500b1 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 6 Jun 2024 15:39:18 -0400 Subject: [PATCH 38/45] remove line --- components/eam/src/physics/crm/pam/pam_state.h | 1 - 1 file changed, 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/pam_state.h b/components/eam/src/physics/crm/pam/pam_state.h index 270fef7f6c23..980fb7f4eb0a 100644 --- a/components/eam/src/physics/crm/pam/pam_state.h +++ b/components/eam/src/physics/crm/pam/pam_state.h @@ -156,7 +156,6 @@ inline void pam_state_set_reference_state( pam::PamCoupler &coupler ) { parallel_for(SimpleBounds<2>(nz+1,nens), YAKL_LAMBDA (int k, int iens) { hmean_pint(k,iens) = 0; if (k < nz) { - hmean_pmid (k,iens) = 0; hmean_rho_d(k,iens) = 0; hmean_rho_v(k,iens) = 0; hmean_rho_c(k,iens) = 0; From 30faaba50681d0821a80ccb9b5592c9d81b6074e Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 6 Jun 2024 15:39:41 -0400 Subject: [PATCH 39/45] make variables const --- .../eam/src/physics/crm/pam/pam_variance_transport.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_variance_transport.h b/components/eam/src/physics/crm/pam/pam_variance_transport.h index 1eb7779f3469..70b985cc2b7b 100644 --- a/components/eam/src/physics/crm/pam/pam_variance_transport.h +++ b/components/eam/src/physics/crm/pam/pam_variance_transport.h @@ -36,11 +36,11 @@ inline void pam_variance_transport_diagnose( pam::PamCoupler &coupler ) { auto ny = coupler.get_option("crm_ny"); auto nx = coupler.get_option("crm_nx"); //------------------------------------------------------------------------------------------------ - auto temp = dm_device.get("temp" ); - auto rhov = dm_device.get("water_vapor"); - auto rhoc = dm_device.get("cloud_water"); - auto rhoi = dm_device.get("ice" ); - auto uvel = dm_device.get("uvel" ); + auto temp = dm_device.get("temp" ); + auto rhov = dm_device.get("water_vapor"); + auto rhoc = dm_device.get("cloud_water"); + auto rhoi = dm_device.get("ice" ); + auto uvel = dm_device.get("uvel" ); auto vt_temp = dm_device.get("vt_temp" ); auto vt_rhov = dm_device.get("vt_rhov" ); auto vt_uvel = dm_device.get("vt_uvel" ); From f9e6cc20fbe98699adcbae004932871bf198c943 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 6 Jun 2024 15:40:13 -0400 Subject: [PATCH 40/45] switch to std::abs --- components/eam/src/physics/crm/pam/pam_driver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index f8f575ec8a1d..bed8980ea67d 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -61,8 +61,8 @@ inline int pam_driver_set_subcycle_timestep( pam::PamCoupler &coupler, real crm_ }); // calculate max U and W parallel_for( SimpleBounds<4>(crm_nz,crm_ny,crm_nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) { - yakl::atomicMax(uvel_max(k,n), sqrt(uvel(k,j,i,n)*uvel(k,j,i,n)) ); - yakl::atomicMax(wvel_max(k,n), fabs(wvel(k,j,i,n)) ); + yakl::atomicMax(uvel_max(k,n), std::abs(uvel(k,j,i,n)) ); + yakl::atomicMax(wvel_max(k,n), std::abs(wvel(k,j,i,n)) ); }); // find max CFL between horizontal and vertical CFL values parallel_for( SimpleBounds<2>(crm_nz,nens) , YAKL_LAMBDA (int k, int n) { From 2590f226f43c316ecb224649c00e597a7dce8b43 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Thu, 6 Jun 2024 15:42:58 -0400 Subject: [PATCH 41/45] remove pmax --- components/eam/src/physics/crm/pam/pam_driver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index bed8980ea67d..5fd1c5071d16 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -73,7 +73,7 @@ inline int pam_driver_set_subcycle_timestep( pam::PamCoupler &coupler, real crm_ cfl_max(k,n) = max(cfl_u,cfl_w); }); // calculate final CFL across ensemble - real cfl_loc = pmax(cfl_max.data()); + real cfl_loc = yakl::intrinsics::maxval(cfl_max); cfl = max(cfl,cfl_loc); // update number of subcycles and time step num_subcycle = max(num_subcycle,max(1,static_cast(ceil(cfl/0.7)))); From 56e9d2efd44924b09ed99eca0331d97999ab3326 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Wed, 10 Jul 2024 13:58:10 -0500 Subject: [PATCH 42/45] PAM fixes following SCREAM merge into E3SM --- components/eam/src/physics/crm/pam/CMakeLists.txt | 6 ++++-- components/eam/src/physics/crm/pam/pam_driver.cpp | 6 +++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/components/eam/src/physics/crm/pam/CMakeLists.txt b/components/eam/src/physics/crm/pam/CMakeLists.txt index 95075debd4fd..4a8f88f669b0 100644 --- a/components/eam/src/physics/crm/pam/CMakeLists.txt +++ b/components/eam/src/physics/crm/pam/CMakeLists.txt @@ -11,6 +11,9 @@ set(PAM_DRIVER_SRC add_library(pam_driver ${PAM_DRIVER_SRC}) +set(SCREAM_HOME ${CMAKE_CURRENT_SOURCE_DIR}/../../../../../..) +add_library(eamxx_physics INTERFACE ${SCREAM_HOME}/components/eamxx/src/physics/) + if (${USE_CUDA}) # PAM will be CUDA-linked with device symbols resolved at library creation set_target_properties(pam_driver @@ -37,7 +40,6 @@ set(PAM_MICRO p3) set(PAM_SGS shoc) set(PAM_RAD forced) set(PAM_SCREAM_USE_CXX True) -set(SCREAM_HOME ${CMAKE_CURRENT_SOURCE_DIR}/../../../../../..) # removing this once the cime cmake_macros aren't using obsolete # Kokkos settings like KOKKOS_OPTIONS. @@ -124,7 +126,7 @@ target_compile_options(pam_driver PUBLIC ) if (PAM_SCREAM_USE_CXX) - target_link_libraries(pam_driver pam_core physics dynamics pam_scream_cxx_interfaces ekat p3 shoc physics_share scream_share) + target_link_libraries(pam_driver pam_core physics dynamics pam_scream_cxx_interfaces ekat p3 shoc physics_share scream_share eamxx_physics) else() target_link_libraries(pam_driver pam_core physics dynamics ) endif() diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index 5fd1c5071d16..c7038343f08a 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -70,13 +70,13 @@ inline int pam_driver_set_subcycle_timestep( pam::PamCoupler &coupler, real crm_ real crm_dz = input_zint(k_gcm,n) - input_zint(k_gcm+1,n); real cfl_u = uvel_max(k,n)*crm_dt_fixed/crm_dx; real cfl_w = wvel_max(k,n)*crm_dt_fixed/crm_dz; - cfl_max(k,n) = max(cfl_u,cfl_w); + cfl_max(k,n) = std::max(cfl_u,cfl_w); }); // calculate final CFL across ensemble real cfl_loc = yakl::intrinsics::maxval(cfl_max); - cfl = max(cfl,cfl_loc); + cfl = std::max(cfl,cfl_loc); // update number of subcycles and time step - num_subcycle = max(num_subcycle,max(1,static_cast(ceil(cfl/0.7)))); + num_subcycle = std::max(num_subcycle,std::max(1,static_cast(ceil(cfl/0.7)))); real crm_dt_subcycle = crm_dt_fixed / num_subcycle; coupler.set_option("crm_dt",crm_dt_subcycle); // check for excessive subcylcing - don't exit, just print From 48fb3c0faee1391a078c76dcd1d80b2cf2a5de38 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Wed, 31 Jul 2024 17:17:03 -0500 Subject: [PATCH 43/45] update PAM submodule --- components/eam/src/physics/crm/pam/external | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/external b/components/eam/src/physics/crm/pam/external index 5931666c95a3..3ea20ad38f28 160000 --- a/components/eam/src/physics/crm/pam/external +++ b/components/eam/src/physics/crm/pam/external @@ -1 +1 @@ -Subproject commit 5931666c95a3b75ae63fbe5bff6e743e74cb2c0a +Subproject commit 3ea20ad38f286730973429e0d491420a6f599f11 From 7d69007b0ec4a15872584ebd6faa38642cfe4729 Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Wed, 31 Jul 2024 17:17:35 -0500 Subject: [PATCH 44/45] add negative mass checks for pam_debug.h --- components/eam/src/physics/crm/pam/pam_debug.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/components/eam/src/physics/crm/pam/pam_debug.h b/components/eam/src/physics/crm/pam/pam_debug.h index 7e54f437f691..a168710afc6d 100644 --- a/components/eam/src/physics/crm/pam/pam_debug.h +++ b/components/eam/src/physics/crm/pam/pam_debug.h @@ -134,7 +134,9 @@ void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) { const auto is_neg_t_atm = temp(k,j,i,iens)<0; const auto is_neg_d_atm = rhod(k,j,i,iens)<0; const auto is_neg_q_atm = rhov(k,j,i,iens)<0; - if ( is_neg_t_atm || is_neg_q_atm || is_neg_d_atm ) { + const auto is_neg_c_atm = rhoc(k,j,i,iens)<0; + const auto is_neg_i_atm = rhoi(k,j,i,iens)<0; + if ( is_neg_t_atm || is_neg_q_atm || is_neg_d_atm || is_neg_c_atm || is_neg_i_atm ) { auto phis = input_phis(iens)/grav; printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g u:%8.2g w:%8.2g \n", nstep,id,k,i,iens,lat(iens),lon(iens),phis, From a39296ce127815b05ef80253d1a986e75fc75dae Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Tue, 20 Aug 2024 16:46:20 -0400 Subject: [PATCH 45/45] re-enable MMF2 test --- cime_config/tests.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cime_config/tests.py b/cime_config/tests.py index 059c280e2bb4..0248077934d6 100644 --- a/cime_config/tests.py +++ b/cime_config/tests.py @@ -370,7 +370,7 @@ "ERP_Ln9.ne4pg2_oQU480.WCYCL20TRNS-MMF1.allactive-mmf_fixed_subcycle", "ERS_Ln9.ne4pg2_ne4pg2.FRCE-MMF1.eam-cosp_nhtfrq9", "SMS_Ln5.ne4_ne4.FSCM-ARM97-MMF1", - # "SMS_Ln3.ne4pg2_ne4pg2.F2010-MMF2", - temporarily disabled due to ongoing dev issues + "SMS_Ln3.ne4pg2_oQU480.F2010-MMF2", ) },