From 13f3fb5efa9b0d1f53e4cc0a4c0b468a0749a0fb Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Thu, 15 Aug 2024 13:56:19 -0600 Subject: [PATCH] Compute Obukhov length based on virtual temperature/flux --- Source/BoundaryConditions/ABLMost.H | 2 +- Source/BoundaryConditions/ABLMost.cpp | 23 +++--- Source/BoundaryConditions/MOSTAverage.cpp | 3 + Source/BoundaryConditions/MOSTStress.H | 92 ++++++++++++++++------- 4 files changed, 82 insertions(+), 38 deletions(-) diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 3be1aa3c0..267068684 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -217,7 +217,7 @@ public: t_star[lev]->setVal(1.E34); q_star[lev] = std::make_unique(ba2d,dm,ncomp,ng); - q_star[lev]->setVal(1.E34); + q_star[lev]->setVal(0.0); // default to dry olen[lev] = std::make_unique(ba2d,dm,ncomp,ng); olen[lev]->setVal(1.E34); diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index 3425a7b4d..3af23510d 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -138,8 +138,10 @@ ABLMost::compute_fluxes (const int& lev, bool is_land) { // Pointers to the computed averages - const auto *const tm_ptr = m_ma.get_average(lev,2); - const auto *const umm_ptr = m_ma.get_average(lev,5); + const auto *const tm_ptr = m_ma.get_average(lev,2); // potential temperature + const auto *const qvm_ptr = m_ma.get_average(lev,3); // water vapor mixing ratio + const auto *const tvm_ptr = m_ma.get_average(lev,4); // virtual potential temperature + const auto *const umm_ptr = m_ma.get_average(lev,5); // horizontal velocity magnitude for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi) { @@ -147,10 +149,13 @@ ABLMost::compute_fluxes (const int& lev, auto u_star_arr = u_star[lev]->array(mfi); auto t_star_arr = t_star[lev]->array(mfi); + auto q_star_arr = q_star[lev]->array(mfi); auto t_surf_arr = t_surf[lev]->array(mfi); auto olen_arr = olen[lev]->array(mfi); const auto tm_arr = tm_ptr->array(mfi); + const auto tvm_arr = tvm_ptr->array(mfi); + const auto qvm_arr = qvm_ptr->array(mfi); const auto umm_arr = umm_ptr->array(mfi); const auto z0_arr = z_0[lev].array(); @@ -164,13 +169,13 @@ ABLMost::compute_fluxes (const int& lev, ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - if (is_land && lmask_arr(i,j,k) == 1) { - most_flux.iterate_flux(i, j, k, max_iters, z0_arr, umm_arr, tm_arr, - u_star_arr, t_star_arr, t_surf_arr, olen_arr, - Hwave_arr, Lwave_arr, eta_arr); - } else if (!is_land && lmask_arr(i,j,k) == 0) { - most_flux.iterate_flux(i, j, k, max_iters, z0_arr, umm_arr, tm_arr, - u_star_arr, t_star_arr, t_surf_arr, olen_arr, + if (( is_land && lmask_arr(i,j,k) == 1) || + (!is_land && lmask_arr(i,j,k) == 0)) + { + most_flux.iterate_flux(i, j, k, max_iters, + z0_arr, umm_arr, tm_arr, tvm_arr, qvm_arr, + u_star_arr, t_star_arr, q_star_arr, // to be updated + t_surf_arr, olen_arr, // to be updated Hwave_arr, Lwave_arr, eta_arr); } }); diff --git a/Source/BoundaryConditions/MOSTAverage.cpp b/Source/BoundaryConditions/MOSTAverage.cpp index f3af8519a..c5d57f2a2 100644 --- a/Source/BoundaryConditions/MOSTAverage.cpp +++ b/Source/BoundaryConditions/MOSTAverage.cpp @@ -124,6 +124,9 @@ MOSTAverage::MOSTAverage (Vector geom, m_averages[lev][iavg]->setVal(1.E34); } + // Default to dry + m_averages[lev][3]->setVal(0.0); + if (m_rotate) { m_rot_fields[lev][2] = std::make_unique(ba,dm,ncomp,ng); m_rot_fields[lev][3] = std::make_unique(ba,dm,ncomp,ng); diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index ebdf08ca0..2fbede49f 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -85,8 +85,11 @@ struct adiabatic const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& /*tm_arr*/, + const amrex::Array4& /*tvm_arr*/, + const amrex::Array4& /*qvm_arr*/, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& /*q_star_arr*/, const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -128,8 +131,11 @@ struct adiabatic_charnock const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& /*tm_arr*/, + const amrex::Array4& /*tvm_arr*/, + const amrex::Array4& /*qvm_arr*/, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& /*q_star_arr*/, const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -184,8 +190,11 @@ struct adiabatic_mod_charnock const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& /*tm_arr*/, + const amrex::Array4& /*tvm_arr*/, + const amrex::Array4& /*qvm_arr*/, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& /*q_star_arr*/, const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -237,8 +246,11 @@ struct adiabatic_wave_coupled const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& /*tm_arr*/, + const amrex::Array4& /*tvm_arr*/, + const amrex::Array4& /*qvm_arr*/, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& /*q_star_arr*/, const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& olen_arr, const amrex::Array4& Hwave_arr, @@ -299,8 +311,11 @@ struct surface_flux const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& t_surf_arr, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -309,6 +324,7 @@ struct surface_flux { int iter = 0; amrex::Real ustar = 0.0; + amrex::Real tflux = 0.0; amrex::Real zeta = 0.0; amrex::Real psi_m = 0.0; amrex::Real psi_h = 0.0; @@ -316,9 +332,8 @@ struct surface_flux u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k)); do { ustar = u_star_arr(i,j,k); - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * mdata.surf_temp_flux); + tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k); + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta); @@ -363,8 +378,11 @@ struct surface_flux_charnock const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& t_surf_arr, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -373,6 +391,7 @@ struct surface_flux_charnock { int iter = 0; amrex::Real ustar = 0.0; + amrex::Real tflux = 0.0; amrex::Real z0 = 0.0; amrex::Real zeta = 0.0; amrex::Real psi_m = 0.0; @@ -382,9 +401,8 @@ struct surface_flux_charnock do { ustar = u_star_arr(i,j,k); z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar; - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * mdata.surf_temp_flux); + tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k); + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta); @@ -431,8 +449,11 @@ struct surface_flux_mod_charnock const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& t_surf_arr, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -441,6 +462,7 @@ struct surface_flux_mod_charnock { int iter = 0; amrex::Real ustar = 0.0; + amrex::Real tflux = 0.0; amrex::Real z0 = 0.0; amrex::Real zeta = 0.0; amrex::Real psi_m = 0.0; @@ -450,9 +472,8 @@ struct surface_flux_mod_charnock do { ustar = u_star_arr(i,j,k); z0 = std::exp( (2.7*ustar - 1.8/mdata.Cnk_b) / (ustar + 0.17/mdata.Cnk_b) ); - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * mdata.surf_temp_flux); + tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k); + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta); @@ -496,8 +517,11 @@ struct surface_flux_wave_coupled const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& t_surf_arr, const amrex::Array4& olen_arr, const amrex::Array4& Hwave_arr, @@ -506,6 +530,7 @@ struct surface_flux_wave_coupled { int iter = 0; amrex::Real ustar = 0.0; + amrex::Real tflux = 0.0; amrex::Real z0 = 0.0; amrex::Real zeta = 0.0; amrex::Real psi_m = 0.0; @@ -521,9 +546,8 @@ struct surface_flux_wave_coupled ustar = u_star_arr(i,j,k); z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * mdata.surf_temp_flux); + tflux = mdata.surf_temp_flux*(1 + 0.61*qvm_arr(i,j,k)) - 0.61*tm_arr(i,j,k)*ustar*q_star_arr(i,j,k); + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta); @@ -570,8 +594,11 @@ struct surface_temp const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& t_surf_arr, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -589,10 +616,10 @@ struct surface_temp do { ustar = u_star_arr(i,j,k); tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa / - (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h); - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * tflux); + (std::log(mdata.zref / z0_arr(i,j,k)) - psi_h); // + tflux *= (1 + 0.61*qvm_arr(i,j,k)); + tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta); @@ -636,9 +663,12 @@ struct surface_temp_charnock const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, const amrex::Array4& t_surf_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, const amrex::Array4& /*Lwave_arr*/, @@ -657,10 +687,10 @@ struct surface_temp_charnock ustar = u_star_arr(i,j,k); z0 = (mdata.Cnk_a / mdata.gravity) * ustar * ustar; tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa / - (std::log(mdata.zref / z0) - psi_h); - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * tflux); + (std::log(mdata.zref / z0) - psi_h); // + tflux *= (1 + 0.61*qvm_arr(i,j,k)); + tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta); @@ -706,8 +736,11 @@ struct surface_temp_mod_charnock const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& t_surf_arr, const amrex::Array4& olen_arr, const amrex::Array4& /*Hwave_arr*/, @@ -727,10 +760,10 @@ struct surface_temp_mod_charnock ustar = u_star_arr(i,j,k); z0 = std::exp( (2.7*ustar - 1.8/mdata.Cnk_b) / (ustar + 0.17/mdata.Cnk_b) ); tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa / - (std::log(mdata.zref / z0) - psi_h); - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * tflux); + (std::log(mdata.zref / z0) - psi_h); // + tflux *= (1 + 0.61*qvm_arr(i,j,k)); + tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta); @@ -773,8 +806,11 @@ struct surface_temp_wave_coupled const amrex::Array4& z0_arr, const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, + const amrex::Array4& tvm_arr, + const amrex::Array4& qvm_arr, const amrex::Array4& u_star_arr, const amrex::Array4& t_star_arr, + const amrex::Array4& q_star_arr, const amrex::Array4& t_surf_arr, const amrex::Array4& olen_arr, const amrex::Array4& Hwave_arr, @@ -800,10 +836,10 @@ struct surface_temp_wave_coupled z0 = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 ) + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max ); tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa / - (std::log(mdata.zref / z0) - psi_h); - // TODO: calculate virtual temperature flux if we have q* - Olen = -ustar * ustar * ustar * tm_arr(i,j,k) / - (mdata.kappa * mdata.gravity * tflux); + (std::log(mdata.zref / z0) - psi_h); // + tflux *= (1 + 0.61*qvm_arr(i,j,k)); + tflux += 0.61*tm_arr(i,j,k) * -ustar*q_star_arr(i,j,k); // ~= + Olen = -ustar * ustar * ustar * tvm_arr(i,j,k) / (mdata.kappa * mdata.gravity * tflux); zeta = mdata.zref / Olen; psi_m = sfuns.calc_psi_m(zeta); psi_h = sfuns.calc_psi_h(zeta);