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", ) }, diff --git a/components/eam/cime_config/config_component.xml b/components/eam/cime_config/config_component.xml index 1a961b84affe..e2ef2ef7173e 100755 --- a/components/eam/cime_config/config_component.xml +++ b/components/eam/cime_config/config_component.xml @@ -71,14 +71,15 @@ -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 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_dx 2000 -crm_nx 64 -crm_ny 1 -crm_nx_rad 4 -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 - -crm_nx_rad 4 -crm_ny_rad 1 -rad rrtmgp -rrtmgpxx - -use_MMF_VT + -rad rrtmgp -rrtmgpxx + -use_MMF_VT -use_MMF_ESMT -aquaplanet -aquaplanet -rce diff --git a/components/eam/src/physics/crm/crm_physics.F90 b/components/eam/src/physics/crm/crm_physics.F90 index 61ef03f8ca65..39ed6a0876c7 100644 --- a/components/eam/src/physics/crm/crm_physics.F90 +++ b/components/eam/src/physics/crm/crm_physics.F90 @@ -1379,8 +1379,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, '' ) 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/external b/components/eam/src/physics/crm/pam/external index 87731d56aeee..3ea20ad38f28 160000 --- a/components/eam/src/physics/crm/pam/external +++ b/components/eam/src/physics/crm/pam/external @@ -1 +1 @@ -Subproject commit 87731d56aeee4bfae4750e17732828ba186367a7 +Subproject commit 3ea20ad38f286730973429e0d491420a6f599f11 diff --git a/components/eam/src/physics/crm/pam/pam_accelerate.h b/components/eam/src/physics/crm/pam/pam_accelerate.h index a48ca628f691..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,28 +89,28 @@ 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_q ("hmean_q", 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 ) { }); //------------------------------------------------------------------------------------------------ } - diff --git a/components/eam/src/physics/crm/pam/pam_debug.h b/components/eam/src/physics/crm/pam/pam_debug.h index 113944365d7c..a168710afc6d 100644 --- a/components/eam/src/physics/crm/pam/pam_debug.h +++ b/components/eam/src/physics/crm/pam/pam_debug.h @@ -31,17 +31,23 @@ inline void pam_debug_init( pam::PamCoupler &coupler ) { 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 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) { @@ -50,6 +56,8 @@ inline void pam_debug_init( pam::PamCoupler &coupler ) { 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); }); //------------------------------------------------------------------------------------------------ } @@ -59,81 +67,185 @@ 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 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"); - 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; + //------------------------------------------------------------------------------------------------ + // 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_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", - 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), - 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) - ); + 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_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", - 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), - 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) - ); + 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; + 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, + 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 large temperature drops + 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) + ); + } + } + //-------------------------------------------------------------------------- // 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"); @@ -195,18 +307,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; diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index bfba64696163..c7038343f08a 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,12 +13,11 @@ #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" @@ -27,13 +25,78 @@ #include "pam_debug.h" 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), 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) { + 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) = std::max(cfl_u,cfl_w); + }); + // calculate final CFL across ensemble + real cfl_loc = yakl::intrinsics::maxval(cfl_max); + cfl = std::max(cfl,cfl_loc); + // update number of subcycles and time step + 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 + 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"); @@ -52,16 +117,19 @@ 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); //------------------------------------------------------------------------------------------------ + // 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); + 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",30.0); + //------------------------------------------------------------------------------------------------ // Allocate the coupler state and retrieve host/device data managers coupler.allocate_coupler_state( crm_nz , crm_ny , crm_nx , nens ); @@ -89,8 +157,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) { @@ -150,12 +220,15 @@ 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 //------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------ // 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) { @@ -164,50 +237,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); } - 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); } - - // 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, 3, 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); - - // 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, 5, 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); } + // 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")) { @@ -219,7 +300,6 @@ extern "C" void pam_driver() { pam_radiation_timestep_aggregation(coupler); pam_statistics_timestep_aggregation(coupler); - etime_crm += crm_dt; nstep += 1; } //------------------------------------------------------------------------------------------------ diff --git a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h index 83dc83252c68..b632b0e49093 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" ); //------------------------------------------------------------------------------------------------ - real constexpr hd_timescale = 10.0; // damping time scale [sec] + #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] + #endif //------------------------------------------------------------------------------------------------ real4d hd_temp("hd_temp",nz,ny,nx,nens); real4d hd_rhod("hd_rhod",nz,ny,nx,nens); diff --git a/components/eam/src/physics/crm/pam/pam_state.h b/components/eam/src/physics/crm/pam/pam_state.h index e8f6e4ab668a..980fb7f4eb0a 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 @@ -144,10 +153,9 @@ 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; hmean_rho_d(k,iens) = 0; hmean_rho_v(k,iens) = 0; hmean_rho_c(k,iens) = 0; @@ -157,33 +165,37 @@ 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) { - 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); - 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); - 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 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); + }); + + // check that interface pressure is reasonable + 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 -- " + "pint:%12.4f pmid:%12.4f \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); @@ -258,9 +270,17 @@ 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; + 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,12 +289,31 @@ 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); + // 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); + // }); + //------------------------------------------------------------------------------------------------ } 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..70b985cc2b7b 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; @@ -34,9 +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 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" ); @@ -62,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); }); //------------------------------------------------------------------------------------------------ @@ -83,6 +89,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; @@ -98,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" ); @@ -114,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 ; }); //------------------------------------------------------------------------------------------------ } @@ -131,13 +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. - real constexpr pert_scale_min = 1.0 - 0.005; - real constexpr pert_scale_max = 1.0 + 0.005; + 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" ); @@ -157,19 +164,22 @@ 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) { // 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 + 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_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 ); @@ -208,9 +218,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" ); @@ -229,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