Compute turbulent diffusivities using the Mellor-Yamada-Janjic (MYJ) closure model.
70 Real EPSQ1 = std::sqrt(EPSQ2);
81 Real A1 =
Real(0.659888514560862645);
82 Real A2x =
Real(0.6574209922667784586);
83 Real B1 =
Real(11.87799326209552761);
84 Real B2 =
Real(7.226971804046074028);
85 Real C1 =
Real(0.000830955950095854396);
92 Real ANMH = -
Real(9.)*A1*A2x*A2x*BTG*BTG;
100 Real BSHH =
Real(9.)*A1*A2x*A2x*BTG;
106 Real AEQH =
Real(9.)*A1*A2x*A2x*B1*BTG*BTG
109 +
Real(18.)*A1*A1*A2x*(B2-
three*A2x)*BTG;
111 Real REQU = -AEQH/AEQM;
113 Real EPSGM = REQU*EPSGH;
115 Real UBRYL = (
Real(18.)*REQU*A1*A1*A2x*B2*C1*BTG +
Real(9.)*A1*A2x*A2x*B2*BTG*BTG)
120 Real AUBH =
Real(27.)*A1*A2x*A2x*B2*BTG*BTG-ADNH*UBRY3;
121 Real AUBM =
Real(54.)*A1*A1*A2x*B2*C1*BTG-ADNM*UBRY3;
123 Real BUBM =
Real(18.)*A1*A1*C1-BDNM*UBRY3;
128 #pragma omp parallel if (Gpu::notInLaunchRegion())
130 for (MFIter mfi(eddyViscosity,
false); mfi.isValid(); ++mfi) {
132 const Box& bx = mfi.validbox();
133 const Array4<Real >& cell_data = cons_in.array(mfi);
134 const Array4<Real >& K_turb = eddyViscosity.array(mfi);
135 const Array4<Real const>& uvel =
xvel.array(mfi);
136 const Array4<Real const>& vvel =
yvel.array(mfi);
139 const Box& dbx = geom.Domain();
143 int klo = bx.smallEnd(2);
144 int khi = bx.bigEnd(2);
145 Box planexy = makeSlab(bx,2,
klo);
148 const GeometryData gdata = geom.data();
151 const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
152 FArrayBox qturb(bx,1,The_Async_Arena());
153 FArrayBox qintegral(xybx,2,The_Async_Arena());
154 IArrayBox pbl_k(xybx,1,The_Async_Arena());
155 qintegral.setVal<RunOn::Device>(0);
156 pbl_k.setVal<RunOn::Device>(
khi);
157 const Array4<Real> qint = qintegral.array();
158 const Array4<Real> qvel = qturb.array();
159 const Array4<int> k_arr = pbl_k.array();
162 const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
164 const auto&
dxInv = geom.InvCellSizeArray();
165 int izmin = geom.Domain().smallEnd(2);
166 int izmax = geom.Domain().bigEnd(2);
171 const bool surface_layer_on_zlo = (SurfLayer !=
nullptr);
172 const auto u_star_arr = (surface_layer_on_zlo) ?
173 SurfLayer->get_u_star(level)->const_array(mfi) : Array4<const Real>{};
174 const auto t_star_arr = (surface_layer_on_zlo) ?
175 SurfLayer->get_t_star(level)->const_array(mfi) : Array4<const Real>{};
176 const auto q_star_arr = (surface_layer_on_zlo &&
use_moisture) ?
177 SurfLayer->get_q_star(level)->const_array(mfi) : Array4<const Real>{};
178 const auto l_obuk_arr = (surface_layer_on_zlo) ?
179 SurfLayer->get_olen(level)->const_array(mfi) : Array4<const Real>{};
182 if (use_terrain_fitted_coords) {
183 ParallelFor(planexy, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ) noexcept
186 for (
int k(
klo); k<=
khi; ++k) {
189 qvel(i,j,k) = std::sqrt(q2);
191 k_arr(i,j,0) = std::min(k,k_arr(i,j,0));
196 for (
int k(
klo); k<=k_arr(i,j,0); ++k) {
199 Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*
dz);
200 Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*
dz);
204 ParallelFor(planexy, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ) noexcept
207 for (
int k(
klo); k<=
khi; ++k) {
210 qvel(i,j,k) = std::sqrt(q2);
212 k_arr(i,j,0) = std::min(k,k_arr(i,j,0));
217 for (
int k(
klo); k<=k_arr(i,j,0); ++k) {
219 const Real Zval = (k +
myhalf)*gdata.CellSize(2);
220 Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k));
221 Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k));
227 ParallelFor(planexy, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ) noexcept
230 int kpbl = k_arr(i,j,0);
233 Real l0 = std::max(std::min(ALPHA*qint(i,j,0,0)/qint(i,j,0,1),EL0MAX),EL0MIN);
236 for (
int k(
klo); k<=
khi; ++k) {
239 Real dthetavdz, dudz, dvdz;
241 uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
242 c_ext_dir_on_zlo, c_ext_dir_on_zhi,
243 u_ext_dir_on_zlo, u_ext_dir_on_zhi,
244 v_ext_dir_on_zlo, v_ext_dir_on_zhi,
245 dthetavdz, dudz, dvdz,
246 moisture_indices, pbl_strat);
250 if (surface_layer_on_zlo && k == izmin) {
251 const Real zval = use_terrain_fitted_coords ?
253 (k +
myhalf)*gdata.CellSize(2);
256 const Real qv = (moisture_indices.
qv >= 0) ?
257 cell_data(i,j,k,moisture_indices.
qv) /
rho :
zero;
259 const Real l_obuk = std::copysign(std::max(std::fabs(l_obuk_arr(i,j,0)),l_min),l_obuk_arr(i,j,0));
261 sl.
u_star = u_star_arr(i,j,0);
263 (q_star_arr) ? q_star_arr(i,j,0) :
zero,
266 sl.
zeta = zval / l_obuk;
271 Real GML = std::max(dudz*dudz + dvdz*dvdz, EPSGM);
274 Real GHL = dthetavdz;
275 if (std::fabs(GHL)<=EPSGH) { GHL=EPSGH; }
280 if (GML/GHL <= REQU) {
283 Real AUBR = (AUBM*GML+AUBH*GHL)*GHL;
284 Real BUBR = BUBM*GML+BUBH*GHL;
287 ELM = std::max(std::sqrt(ELOQ2X*qvel(i,j,k)*qvel(i,j,k)),EPSL);
290 Real ADEN = (ADNM*GML+ADNH*GHL)*GHL;
291 Real BDEN = BDNM*GML+BDNH*GHL;
293 Real ELOQ2X =
one/(QOL2UN+EPSRU);
294 ELM = std::max(std::sqrt(ELOQ2X*qvel(i,j,k)*qvel(i,j,k)),EPSL);
300 L = std::min((met_h_zeta/
dxInv[2])*ELFC, ELM);
303 : (k +
myhalf)*gdata.CellSize(2);
304 L = std::min(l0*d_kappa*zval / (d_kappa*zval + l0), ELM);
308 Real AEQU = (AEQM*GML+AEQH*GHL)*GHL;
309 Real BEQU = BEQM*GML+BEQH*GHL;
313 if ( ((GML+GHL*GHL)<=EPSTRB) ||
314 ((GHL>=EPSGH) && ((GML/GHL)<=REQU)) ||
319 Real ANUM=(ANMM*GML+ANMH*GHL)*GHL;
320 Real BNUM= BNMM*GML+BNMH*GHL;
322 Real ADEN=(ADNM*GML+ADNH*GHL)*GHL;
323 Real BDEN= BDNM*GML+BDNH*GHL;
326 Real ARHS=-(ANUM*BDEN-BNUM*ADEN)*
two;
330 Real DLOQ1=L/qvel(i,j,k);
333 Real ELOQ11=std::sqrt(ELOQ21);
334 Real ELOQ31=ELOQ21*ELOQ11;
335 Real ELOQ41=ELOQ21*ELOQ21;
336 Real ELOQ51=ELOQ21*ELOQ31;
338 Real RDEN1=
one/(ADEN*ELOQ41+BDEN*ELOQ21+CDEN);
340 Real RHSP1=(ARHS*ELOQ51+BRHS*ELOQ31+CRHS*ELOQ11)*RDEN1*RDEN1;
342 Real DTTURBL =
static_cast<Real>(dt);
343 Real ELOQ12 = amrex::max(ELOQ11+(DLOQ1-ELOQ11)*std::exp(RHSP1*DTTURBL),EPS1);
345 Real ELOQ22=ELOQ12*ELOQ12;
346 Real ELOQ32=ELOQ22*ELOQ12;
347 Real ELOQ42=ELOQ22*ELOQ22;
348 Real ELOQ52=ELOQ22*ELOQ32;
350 Real RDEN2=
one/(ADEN*ELOQ42+BDEN*ELOQ22+CDEN);
351 Real RHS2 =-(ANUM*ELOQ42+BNUM*ELOQ22)*RDEN2+RB1;
352 Real RHSP2= (ARHS*ELOQ52+BRHS*ELOQ32+CRHS*ELOQ12)*RDEN2*RDEN2;
353 Real RHST2=RHS2/RHSP2;
355 Real ELOQ13 = amrex::max(ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)*std::exp(RHSP2*DTTURBL),EPS1);
359 qvel(i,j,k) = std::max(L/ELOQN,EPSQ1);
360 if (qvel(i,j,k)==EPSQ1) {L = EPSL; }
377 cell_data(i,j,k,
RhoKE_comp) =
myhalf*cell_data(i,j,k,
Rho_comp)*qvel(i,j,k)*qvel(i,j,k);
380 Real ELOQ2 = L*L/(qvel(i,j,k)*qvel(i,j,k));
381 Real ELOQ4 = ELOQ2*ELOQ2;
384 Real ADEN=(ADNM*GML+ADNH*GHL)*GHL;
385 Real BDEN= BDNM*GML+BDNH*GHL;
392 Real BESH=BSHM*GML+BSHH*GHL;
395 Real RDEN=
one/(ADEN*ELOQ4+BDEN*ELOQ2+CDEN);
398 Real SM=(BESM*ELOQ2+CESM)*RDEN;
399 Real SH=(BESH*ELOQ2+CESH)*RDEN;
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:67
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:68
#define Rho_comp
Definition: ERF_IndexDefines.H:42
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:43
#define RhoKE_comp
Definition: ERF_IndexDefines.H:44
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
const int klo
Definition: ERF_InitCustomPert_ABL.H:75
const bool use_moisture
Definition: ERF_InitCustomPert_ABL.H:71
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:39
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void ComputeVerticalDerivativesPBL(int i, int j, int k, const amrex::Array4< const amrex::Real > &uvel, const amrex::Array4< const amrex::Real > &vvel, const amrex::Array4< const amrex::Real > &cell_data, const int izmin, const int izmax, const PBLDerivativeDzInv &dz_inv, const bool c_ext_dir_on_zlo, const bool c_ext_dir_on_zhi, const bool u_ext_dir_on_zlo, const bool u_ext_dir_on_zhi, const bool v_ext_dir_on_zlo, const bool v_ext_dir_on_zhi, amrex::Real &dthetadz, amrex::Real &dudz, amrex::Real &dvdz, const MoistureComponentIndices &moisture_indices, const StratType strat=StratType::thetav)
Definition: ERF_PBLModels.H:304
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void ApplySurfaceLayerGradientsPBL(const PBLSurfaceLayerGradient &sl, amrex::Real &dthetadz, amrex::Real &dudz, amrex::Real &dvdz)
Definition: ERF_PBLModels.H:423
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeVirtualTStarPBL(amrex::Real t_star, amrex::Real q_star, amrex::Real theta, amrex::Real qv, bool use_moisture)
Definition: ERF_PBLModels.H:386
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Zrel_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:758
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:197
@ yvel_bc
Definition: ERF_IndexDefines.H:109
@ cons_bc
Definition: ERF_IndexDefines.H:92
@ xvel_bc
Definition: ERF_IndexDefines.H:108
@ ext_dir
Definition: ERF_IndexDefines.H:300
@ Theta_v
Definition: ERF_IndexDefines.H:253
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:257
@ Q_v
Definition: ERF_IndexDefines.H:256
@ Mom_v
Definition: ERF_IndexDefines.H:252
@ KE_v
Definition: ERF_IndexDefines.H:254
@ rho
Definition: ERF_Kessler.H:25
@ theta
Definition: ERF_Kessler.H:26
@ qv
Definition: ERF_Kessler.H:31
@ xvel
Definition: ERF_IndexDefines.H:218
@ yvel
Definition: ERF_IndexDefines.H:219
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
int qv
water vapor
Definition: ERF_DataStruct.H:249
Functor for inverse vertical spacings for terrain-following grids using cell-center heights.
Definition: ERF_PBLModels.H:663
Definition: ERF_PBLModels.H:366
amrex::Real tstar_v
virtual temperature scale
Definition: ERF_PBLModels.H:368
amrex::Real zeta
stability coordinate z/L
Definition: ERF_PBLModels.H:370
amrex::Real u_star
friction velocity
Definition: ERF_PBLModels.H:367
amrex::Real zval
height of the cell center above the surface
Definition: ERF_PBLModels.H:369
StratType pbl_strat_type
Definition: ERF_TurbStruct.H:819