diff --git a/components/eam/src/physics/crm/crm_physics.F90 b/components/eam/src/physics/crm/crm_physics.F90 index cc65d559bb88..61ef03f8ca65 100644 --- a/components/eam/src/physics/crm/crm_physics.F90 +++ b/components/eam/src/physics/crm/crm_physics.F90 @@ -1302,14 +1302,23 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, call pam_mirror_array_readonly( 'input_tau00', crm_input%tau00 ) ! call pam_mirror_array_readonly( 'input_ul_esmt', crm_input%ul_esmt ) ! call pam_mirror_array_readonly( 'input_vl_esmt', crm_input%vl_esmt ) - ! call pam_mirror_array_readonly( 'input_t_vt', crm_input%t_vt ) - ! call pam_mirror_array_readonly( 'input_q_vt', crm_input%q_vt ) - ! call pam_mirror_array_readonly( 'input_u_vt', crm_input%u_vt ) - call pam_mirror_array_readonly( 'input_nccn_prescribed',crm_input%nccn_prescribed ) call pam_mirror_array_readonly( 'input_nc_nuceat_tend', crm_input%nc_nuceat_tend ) call pam_mirror_array_readonly( 'input_ni_activated', crm_input%ni_activated ) + ! Variance transport inputs and outputs + if (use_MMF_VT) then + call pam_mirror_array_readonly( 'input_vt_t', crm_input%t_vt ) + call pam_mirror_array_readonly( 'input_vt_q', crm_input%q_vt ) + call pam_mirror_array_readonly( 'input_vt_u', crm_input%u_vt ) + call pam_mirror_array_readwrite( 'output_t_vt_tend', crm_output%t_vt_tend ) + call pam_mirror_array_readwrite( 'output_q_vt_tend', crm_output%q_vt_tend ) + call pam_mirror_array_readwrite( 'output_u_vt_tend', crm_output%u_vt_tend ) + call pam_mirror_array_readwrite( 'output_t_vt_ls', crm_output%t_vt_ls ) + call pam_mirror_array_readwrite( 'output_q_vt_ls', crm_output%q_vt_ls ) + call pam_mirror_array_readwrite( 'output_u_vt_ls', crm_output%u_vt_ls ) + end if + call pam_mirror_array_readwrite( 'state_u_wind', crm_state%u_wind ) call pam_mirror_array_readwrite( 'state_v_wind', crm_state%v_wind ) call pam_mirror_array_readwrite( 'state_w_wind', crm_state%w_wind ) @@ -1389,12 +1398,6 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, call pam_mirror_array_readwrite( 'output_qltend', crm_output%qltend, '' ) call pam_mirror_array_readwrite( 'output_qcltend', crm_output%qcltend, '' ) call pam_mirror_array_readwrite( 'output_qiltend', crm_output%qiltend, '' ) - ! call pam_mirror_array_readwrite( 'output_t_vt_tend', crm_output%t_vt_tend, '' ) - ! call pam_mirror_array_readwrite( 'output_q_vt_tend', crm_output%q_vt_tend, '' ) - ! call pam_mirror_array_readwrite( 'output_u_vt_tend', crm_output%u_vt_tend, '' ) - ! call pam_mirror_array_readwrite( 'output_t_vt_ls', crm_output%t_vt_ls, '' ) - ! call pam_mirror_array_readwrite( 'output_q_vt_ls', crm_output%q_vt_ls, '' ) - ! call pam_mirror_array_readwrite( 'output_u_vt_ls', crm_output%u_vt_ls, '' ) call pam_mirror_array_readwrite( 'output_ultend', crm_output%ultend, '' ) call pam_mirror_array_readwrite( 'output_vltend', crm_output%vltend, '' ) ! call pam_mirror_array_readwrite( 'output_tk', crm_output%tk, '' ) @@ -1450,7 +1453,6 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out, call pam_register_dimension('gcm_lev',pver) call pam_set_option('use_MMF_VT', use_MMF_VT_tmp ) - call pam_set_option('MMF_VT_wn_max', MMF_VT_wn_max ) call pam_set_option('use_MMF_ESMT', use_MMF_ESMT_tmp ) call pam_set_option('use_crm_accel', use_crm_accel_tmp ) diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp index f131bbd44f09..bfba64696163 100644 --- a/components/eam/src/physics/crm/pam/pam_driver.cpp +++ b/components/eam/src/physics/crm/pam/pam_driver.cpp @@ -13,6 +13,7 @@ #include "pam_statistics.h" #include "pam_output.h" #include "pam_accelerate.h" +#include "pam_variance_transport.h" #include "sponge_layer.h" #include "surface_friction.h" #include "scream_cxx_interface_finalize.h" @@ -46,6 +47,7 @@ extern "C" void pam_driver() { 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"); + bool use_MMF_VT = coupler.get_option("use_MMF_VT"); bool enable_physics_tend_stats = coupler.get_option("enable_physics_tend_stats"); //------------------------------------------------------------------------------------------------ // set various coupler options @@ -111,6 +113,11 @@ extern "C" void pam_driver() { // initialize variables for CRM mean-state acceleration if (use_crm_accel) { pam_accelerate_init(coupler); } + if (use_MMF_VT) { + pam_variance_transport_init(coupler); + pam_variance_transport_compute_forcing(coupler); + } + // initilize surface "psuedo-friction" (psuedo => doesn't match "real" GCM friction) auto input_tau = dm_host.get("input_tau00").createDeviceCopy(); auto input_bflx = dm_host.get("input_bflxls").createDeviceCopy(); @@ -166,7 +173,8 @@ extern "C" void pam_driver() { if (enable_check_state) { pam_debug_check_state(coupler, 1, nstep); } - // run a PAM time step + // 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); } @@ -237,6 +245,11 @@ extern "C" void pam_driver() { pam_statistics_compute_means(coupler); pam_statistics_copy_to_host(coupler); + if (use_MMF_VT) { + pam_variance_transport_compute_feedback(coupler); + pam_variance_transport_copy_to_host(coupler); + } + //------------------------------------------------------------------------------------------------ // Finalize and clean up micro .finalize(coupler); diff --git a/components/eam/src/physics/crm/pam/pam_variance_transport.h b/components/eam/src/physics/crm/pam/pam_variance_transport.h new file mode 100644 index 000000000000..b7067631203f --- /dev/null +++ b/components/eam/src/physics/crm/pam/pam_variance_transport.h @@ -0,0 +1,267 @@ +#pragma once + +#include "pam_coupler.h" + +inline void pam_variance_transport_init( 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"); + auto ny = coupler.get_option("crm_ny"); + auto nx = coupler.get_option("crm_nx"); + //------------------------------------------------------------------------------------------------ + dm_device.register_and_allocate("vt_temp", "temperature variance", {nz,nens}, {"z","nens"} ); + dm_device.register_and_allocate("vt_rhov", "water vapor variance", {nz,nens}, {"z","nens"} ); + dm_device.register_and_allocate("vt_uvel", "u momentum variance", {nz,nens}, {"z","nens"} ); + dm_device.register_and_allocate("vt_temp_pert", "temperature perturbation from horz mean", {nz,ny,nx,nens}, {"z","y","x","nens"} ); + dm_device.register_and_allocate("vt_rhov_pert", "water vapor perturbation from horz mean", {nz,ny,nx,nens}, {"z","y","x","nens"} ); + dm_device.register_and_allocate("vt_uvel_pert", "u momentum perturbation from horz mean", {nz,ny,nx,nens}, {"z","y","x","nens"} ); + dm_device.register_and_allocate("vt_temp_forcing_tend", "temperature variance forcing tendency", {nz,nens}, {"z","nens"} ); + dm_device.register_and_allocate("vt_rhov_forcing_tend", "water vapor variance forcing tendency", {nz,nens}, {"z","nens"} ); + dm_device.register_and_allocate("vt_uvel_forcing_tend", "u momentum variance forcing tendency", {nz,nens}, {"z","nens"} ); + //------------------------------------------------------------------------------------------------ +} + +inline void pam_variance_transport_diagnose( 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"); + 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 vt_temp = dm_device.get("vt_temp" ); + auto vt_rhov = dm_device.get("vt_rhov" ); + auto vt_uvel = dm_device.get("vt_uvel" ); + auto vt_temp_pert = dm_device.get("vt_temp_pert" ); + auto vt_rhov_pert = dm_device.get("vt_rhov_pert" ); + auto vt_uvel_pert = dm_device.get("vt_uvel_pert" ); + //------------------------------------------------------------------------------------------------ + // local variables + real2d temp_mean("temp_mean", nz, nens); + real2d rhov_mean("rhov_mean", nz, nens); + real2d uvel_mean("uvel_mean", nz, nens); + //------------------------------------------------------------------------------------------------ + // initialize variance and horizontal mean + parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) { + temp_mean(k,n) = 0.0; + rhov_mean(k,n) = 0.0; + uvel_mean(k,n) = 0.0; + vt_temp(k,n) = 0.0; + vt_rhov(k,n) = 0.0; + vt_uvel(k,n) = 0.0; + }); + //------------------------------------------------------------------------------------------------ + // 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) { + 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( 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) { + 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_uvel_pert(k,j,i,n) = uvel(k,j,i,n) - uvel_mean(k,n); + }); + //------------------------------------------------------------------------------------------------ + // calculate variance + parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int n) { + yakl::atomicAdd( vt_temp(k,n), ( vt_temp_pert(k,j,i,n) * vt_temp_pert(k,j,i,n) ) * r_nx_ny ); + yakl::atomicAdd( vt_rhov(k,n), ( vt_rhov_pert(k,j,i,n) * vt_rhov_pert(k,j,i,n) ) * r_nx_ny ); + yakl::atomicAdd( vt_uvel(k,n), ( vt_uvel_pert(k,j,i,n) * vt_uvel_pert(k,j,i,n) ) * r_nx_ny ); + }); + //------------------------------------------------------------------------------------------------ +} + +inline void pam_variance_transport_compute_forcing( pam::PamCoupler &coupler ) { + 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 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 gcm_nlev = coupler.get_option("gcm_nlev"); + auto gcm_dt = coupler.get_option("gcm_dt"); + //------------------------------------------------------------------------------------------------ + // 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" ); + auto vt_temp_forcing_tend = dm_device.get("vt_temp_forcing_tend"); + auto vt_rhov_forcing_tend = dm_device.get("vt_rhov_forcing_tend"); + auto vt_uvel_forcing_tend = dm_device.get("vt_uvel_forcing_tend"); + auto gcm_vt_temp = dm_host.get("input_vt_t").createDeviceCopy(); + auto gcm_vt_rhov = dm_host.get("input_vt_q").createDeviceCopy(); + auto gcm_vt_uvel = dm_host.get("input_vt_u").createDeviceCopy(); + //------------------------------------------------------------------------------------------------ + // 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 ; + }); + //------------------------------------------------------------------------------------------------ +} + +inline void pam_variance_transport_apply_forcing( pam::PamCoupler &coupler ) { + using yakl::c::parallel_for; + using yakl::c::SimpleBounds; + 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_dt = coupler.get_option("crm_dt"); + //------------------------------------------------------------------------------------------------ + // 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; + //------------------------------------------------------------------------------------------------ + 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" ); + auto vt_temp_pert = dm_device.get("vt_temp_pert" ); + auto vt_rhov_pert = dm_device.get("vt_rhov_pert" ); + auto vt_uvel_pert = dm_device.get("vt_uvel_pert" ); + auto vt_temp_forcing_tend = dm_device.get("vt_temp_forcing_tend"); + auto vt_rhov_forcing_tend = dm_device.get("vt_rhov_forcing_tend"); + auto vt_uvel_forcing_tend = dm_device.get("vt_uvel_forcing_tend"); + //------------------------------------------------------------------------------------------------ + // local variables + real2d temp_pert_scale("temp_pert_scale", nz, nens); + 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) ); } + // 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 ); + uvel_pert_scale(k,n) = std::max( uvel_pert_scale(k,n), pert_scale_min ); + // enforce maximum scaling + temp_pert_scale(k,n) = std::min( temp_pert_scale(k,n), pert_scale_max ); + rhov_pert_scale(k,n) = std::min( rhov_pert_scale(k,n), pert_scale_max ); + uvel_pert_scale(k,n) = std::min( uvel_pert_scale(k,n), pert_scale_max ); + }); + //------------------------------------------------------------------------------------------------ + // apply variance forcing tendency + parallel_for(SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int n) { + real temp_tend_loc = ( temp_pert_scale(k,n) * vt_temp_pert(k,j,i,n) - vt_temp_pert(k,j,i,n) ); + real rhov_tend_loc = ( rhov_pert_scale(k,n) * vt_rhov_pert(k,j,i,n) - vt_rhov_pert(k,j,i,n) ); + real uvel_tend_loc = ( uvel_pert_scale(k,n) * vt_uvel_pert(k,j,i,n) - vt_uvel_pert(k,j,i,n) ); + temp(k,j,i,n) = temp(k,j,i,n) + temp_tend_loc; + rhov(k,j,i,n) = rhov(k,j,i,n) + rhov_tend_loc; + uvel(k,j,i,n) = uvel(k,j,i,n) + uvel_tend_loc; + }); + //------------------------------------------------------------------------------------------------ +} + +inline void pam_variance_transport_compute_feedback( 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 &dm_host = coupler.get_data_manager_host_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 gcm_nlev = coupler.get_option("gcm_nlev"); + auto gcm_dt = coupler.get_option("gcm_dt"); + //------------------------------------------------------------------------------------------------ + // 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" ); + auto gcm_vt_temp = dm_host.get("input_vt_t").createDeviceCopy(); + auto gcm_vt_rhov = dm_host.get("input_vt_q").createDeviceCopy(); + auto gcm_vt_uvel = dm_host.get("input_vt_u").createDeviceCopy(); + //------------------------------------------------------------------------------------------------ + dm_device.register_and_allocate("vt_temp_feedback_tend", "feedback tend of temp variance", {gcm_nlev,nens},{"gcm_lev","nens"}); + dm_device.register_and_allocate("vt_rhov_feedback_tend", "feedback tend of rhov variance", {gcm_nlev,nens},{"gcm_lev","nens"}); + dm_device.register_and_allocate("vt_uvel_feedback_tend", "feedback tend of uvel variance", {gcm_nlev,nens},{"gcm_lev","nens"}); + auto vt_temp_feedback_tend = dm_device.get("vt_temp_feedback_tend" ); + auto vt_rhov_feedback_tend = dm_device.get("vt_rhov_feedback_tend" ); + auto vt_uvel_feedback_tend = dm_device.get("vt_uvel_feedback_tend" ); + //------------------------------------------------------------------------------------------------ + // 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; + if (k_crm("vt_temp_feedback_tend" ); + auto vt_rhov_feedback_tend = dm_device.get("vt_rhov_feedback_tend" ); + auto vt_uvel_feedback_tend = dm_device.get("vt_uvel_feedback_tend" ); + //------------------------------------------------------------------------------------------------ + auto output_vt_temp_tend_host = dm_host.get("output_t_vt_tend" ); + auto output_vt_rhov_tend_host = dm_host.get("output_q_vt_tend" ); + auto output_vt_uvel_tend_host = dm_host.get("output_u_vt_tend" ); + //------------------------------------------------------------------------------------------------ + // Copy the data to host + 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(); + //------------------------------------------------------------------------------------------------ +} +