From 74e265bd91b6ef802420e9ba2b246bd0b653f29c Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 16 Jan 2024 18:00:33 -0700 Subject: [PATCH 1/3] make PBL lengthscale a plotable variable --- Source/ERF.H | 4 +++- Source/IO/Plotfile.cpp | 4 ++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/Source/ERF.H b/Source/ERF.H index 4d2c65827..d3be2249a 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -656,7 +656,9 @@ private: // eddy viscosity "Kmv","Kmh", // eddy diffusivity of heat - "Khv","Khh" + "Khv","Khh", + // mynn pbl lengthscale + "Lpbl" ,"qt", "qp", "qv", "qc", "qi", "qrain", "qsnow", "qgraup" #ifdef ERF_COMPUTE_ERROR ,"xvel_err", "yvel_err", "zvel_err", "pp_err" diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index a688fa5d2..1f75e5742 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -631,6 +631,10 @@ ERF::WritePlotFile (int which, Vector plot_var_names) MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Theta_h,mf_comp,1,0); mf_comp ++; } + if (containerHasElement(plot_var_names, "Lpbl")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::PBL_lengthscale,mf_comp,1,0); + mf_comp ++; + } // NOTE: Protect against accessing non-existent data if (use_moisture) { From 556105663bb5c10b684155615d36f0ba784bc67b Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 16 Jan 2024 18:01:39 -0700 Subject: [PATCH 2/3] fix zeta computation for terrain in MYNN --- Source/Diffusion/PBLModels.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 53b7d3a21..81740311c 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -168,7 +168,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, // First Length Scale AMREX_ASSERT(l_obukhov != 0); int lk = amrex::max(k,0); - const Real zval = gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2); + const Real zval = use_terrain ? Compute_Zrel_AtCellCenter(i,j,lk,z_nd_arr) : gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2); const Real zeta = zval/l_obukhov; Real l_S; if (zeta >= 1.0) { From b0b72654ea5b04a900def9e43318b8ad1ac7cad6 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 16 Jan 2024 18:13:56 -0700 Subject: [PATCH 3/3] terrain initialization for GABLS1 --- Exec/RegTests/GABLS1/prob.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp index 5aed34a33..436da69ce 100644 --- a/Exec/RegTests/GABLS1/prob.cpp +++ b/Exec/RegTests/GABLS1/prob.cpp @@ -55,7 +55,7 @@ Problem::init_custom_pert( amrex::Array4 const& r_hse, amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, - amrex::Array4 const& /*z_cc*/, + amrex::Array4 const& z_cc, amrex::GeometryData const& geomdata, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, @@ -71,21 +71,35 @@ Problem::init_custom_pert( amrex::Print() << "Initializing QKE" << std::endl; } + const bool use_terrain = sc.use_terrain; + if (parms.pert_ref_height > 0) { if (parms.pert_deltaU != 0.0) { amrex::Print() << "Adding divergence-free x-velocity perturbations" << std::endl; + if (use_terrain) { + amrex::Abort("u perturbations not supported for terrain"); + } } if (parms.pert_deltaV != 0.0) { amrex::Print() << "Adding divergence-free y-velocity perturbations" << std::endl; + if (use_terrain) { + amrex::Abort("v perturbations not supported for terrain"); + } } if (parms.T_0_Pert_Mag != 0.0) { amrex::Print() << "Adding random theta perturbations" << std::endl; } if (parms.U_0_Pert_Mag != 0.0) { amrex::Print() << "Adding random x-velocity perturbations" << std::endl; + if (use_terrain) { + amrex::Abort("u perturbations not supported for terrain"); + } } if (parms.V_0_Pert_Mag != 0.0) { amrex::Print() << "Adding random y-velocity perturbations" << std::endl; + if (use_terrain) { + amrex::Abort("v perturbations not supported for terrain"); + } } } @@ -96,7 +110,7 @@ Problem::init_custom_pert( const Real* dx = geomdata.CellSize(); const Real x = prob_lo[0] + (i + 0.5) * dx[0]; const Real y = prob_lo[1] + (j + 0.5) * dx[1]; - const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + const Real z = use_terrain ? z_cc(i,j,k) : prob_lo[2] + (k + 0.5) * dx[2]; // Define a point (xc,yc,zc) at the center of the domain const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]); @@ -196,4 +210,3 @@ Problem::init_custom_pert( } }); } -