From ace1d2acfd6378bb80a2988631175587f4015ca2 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Wed, 14 Aug 2024 23:24:51 -0600 Subject: [PATCH] Make the QKE source term update moisture-aware --- Source/Diffusion/Diffusion.H | 6 ++++-- Source/Diffusion/DiffusionSrcForState_N.cpp | 5 ++++- Source/Diffusion/DiffusionSrcForState_T.cpp | 8 ++++++-- Source/Diffusion/PBLModels.H | 6 ++++-- Source/Diffusion/PBLModels.cpp | 5 +++++ Source/TimeIntegration/ERF_slow_rhs_post.cpp | 6 ++++-- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 4 ++-- 7 files changed, 29 insertions(+), 11 deletions(-) diff --git a/Source/Diffusion/Diffusion.H b/Source/Diffusion/Diffusion.H index e62eca609..4a0aec188 100644 --- a/Source/Diffusion/Diffusion.H +++ b/Source/Diffusion/Diffusion.H @@ -66,7 +66,8 @@ void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& tm_arr, const amrex::GpuArray grav_gpu, const amrex::BCRec* bc_ptr, - const bool use_most); + const bool use_most, + const bool use_moisture); void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, int start_comp, int num_comp, @@ -99,7 +100,8 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& tm_arr, const amrex::GpuArray grav_gpu, const amrex::BCRec* bc_ptr, - const bool use_most); + const bool use_most, + const bool use_moisture); diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index deaa6c243..669dbbce8 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -35,6 +35,7 @@ using namespace amrex; * @param[in] grav_gpu gravity vector * @param[in] bc_ptr container with boundary conditions * @param[in] use_most whether we have turned on MOST BCs + * @param[in] use_moisture whether we account for moisture in the QKE update */ void DiffusionSrcForState_N (const Box& bx, const Box& domain, @@ -63,7 +64,8 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, const Array4& tm_arr, const GpuArray grav_gpu, const BCRec* bc_ptr, - const bool use_most) + const bool use_most, + const bool use_moisture) { BL_PROFILE_VAR("DiffusionSrcForState_N()",DiffusionSrcForState_N); @@ -535,6 +537,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim, mu_turb,cellSizeInv,domain, pbl_mynn_B1_l,tm_arr(i,j,0), + use_moisture, c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi); diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 5ae702293..59b0deda1 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -38,6 +38,7 @@ using namespace amrex; * @param[in] grav_gpu gravity vector * @param[in] bc_ptr container with boundary conditions * @param[in] use_most whether we have turned on MOST BCs + * @param[in] use_moisture whether we account for moisture in the QKE update */ void DiffusionSrcForState_T (const Box& bx, const Box& domain, @@ -71,7 +72,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, const Array4& tm_arr, const GpuArray grav_gpu, const BCRec* bc_ptr, - const bool use_most) + const bool use_most, + const bool use_moisture) { BL_PROFILE_VAR("DiffusionSrcForState_T()",DiffusionSrcForState_T); @@ -656,7 +658,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, { const Real met_h_zeta = detJ(i,j,k); cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim, - mu_turb,cellSizeInv,domain,pbl_mynn_B1_l,tm_arr(i,j,0), + mu_turb,cellSizeInv,domain, + pbl_mynn_B1_l,tm_arr(i,j,0), + use_moisture, c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, diff --git a/Source/Diffusion/PBLModels.H b/Source/Diffusion/PBLModels.H index cadefe26d..31090f174 100644 --- a/Source/Diffusion/PBLModels.H +++ b/Source/Diffusion/PBLModels.H @@ -103,6 +103,7 @@ ComputeQKESourceTerms (int i, int j, int k, const amrex::Box& domain, amrex::Real pbl_mynn_B1_l, const amrex::Real theta_mean, + bool use_moisture, bool c_ext_dir_on_zlo, bool c_ext_dir_on_zhi, bool u_ext_dir_on_zlo, @@ -124,7 +125,8 @@ ComputeQKESourceTerms (int i, int j, int k, c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, - dthetadz, dudz, dvdz); + dthetadz, dudz, dvdz, + use_moisture); // Note: Transport terms due to turbulence and pressure are included when // DiffusionSrcForState_* is called. @@ -132,7 +134,7 @@ ComputeQKESourceTerms (int i, int j, int k, // Shear Production (Note: We store mu_turb, which is 0.5*K_turb) source_term += 4.0*K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz); - // Buoyancy (TODO: implement partial-condensation scheme) + // Buoyancy Production source_term -= 2.0*(CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz; // Dissipation diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 610d20979..8bd9b9796 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -267,6 +267,11 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, K_turb(i,j,k,EddyDiff::Theta_v) = rho * Lturb * qvel(i,j,k) * SH; K_turb(i,j,k,EddyDiff::QKE_v) = rho * Lturb * qvel(i,j,k) * SQ; + // TODO: implement partial-condensation scheme? + // Currently, implementation matches NN09 without rain (i.e., + // the liquid water potential temperature is equal to the + // potential temperature. + // NN09 gives the total water content flux; this assumes that // all the species have the same eddy diffusivity if (update_moist_eddydiff) { diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 71aa0e8e3..39e4912e9 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -394,14 +394,16 @@ void erf_slow_rhs_post (int level, int finest_level, z_nd, ax_arr, ay_arr, az_arr, detJ_arr, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, q1fx_z, q2fx_z, diss, - mu_turb, dc, tc, tm_arr, grav_gpu, bc_ptr_d, use_most); + mu_turb, dc, tc, + tm_arr, grav_gpu, bc_ptr_d, use_most, l_use_moisture); } else { DiffusionSrcForState_N(tbx, domain, start_comp, num_comp, exp_most, u, v, new_cons, cur_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, q1fx_z, q2fx_z, diss, - mu_turb, dc, tc, tm_arr, grav_gpu, bc_ptr_d, use_most); + mu_turb, dc, tc, + tm_arr, grav_gpu, bc_ptr_d, use_most, l_use_moisture); } } // use_diff } // valid slow var diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 0fa8bc791..54889c805 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -481,7 +481,7 @@ if (cell_data(i,j,k,RhoTheta_comp) < 0.) printf("BAD THETA AT %d %d %d %e %e \n" z_nd, ax_arr, ay_arr, az_arr, detJ_arr, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, dc, tc, - tm_arr, grav_gpu, bc_ptr_d, l_use_most); + tm_arr, grav_gpu, bc_ptr_d, l_use_most, l_use_moisture); } else { DiffusionSrcForState_N(bx, domain, n_start, n_comp, l_exp_most, u, v, cell_data, cell_prim, cell_rhs, @@ -489,7 +489,7 @@ if (cell_data(i,j,k,RhoTheta_comp) < 0.) printf("BAD THETA AT %d %d %d %e %e \n" dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, dc, tc, - tm_arr, grav_gpu, bc_ptr_d, l_use_most); + tm_arr, grav_gpu, bc_ptr_d, l_use_most, l_use_moisture); } }