ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_ComputeDiffusivityMYJ.cpp File Reference
#include "ERF_SurfaceLayer.H"
#include "ERF_DirectionSelector.H"
#include "ERF_Diffusion.H"
#include "ERF_Constants.H"
#include "ERF_TurbStruct.H"
#include "ERF_PBLModels.H"
#include <cmath>
Include dependency graph for ERF_ComputeDiffusivityMYJ.cpp:

Functions

void ComputeDiffusivityMYJ (double dt, const MultiFab &xvel, const MultiFab &yvel, MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const std::unique_ptr< MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
 

Function Documentation

◆ ComputeDiffusivityMYJ()

void ComputeDiffusivityMYJ ( double  dt,
const MultiFab &  xvel,
const MultiFab &  yvel,
MultiFab &  cons_in,
MultiFab &  eddyViscosity,
const Geometry &  geom,
const TurbChoice &  turbChoice,
std::unique_ptr< SurfaceLayer > &  SurfLayer,
bool  use_terrain_fitted_coords,
bool  use_moisture,
int  level,
const BCRec *  bc_ptr,
bool  ,
const std::unique_ptr< MultiFab > &  z_phys_nd,
const std::unique_ptr< MultiFab > &  z_phys_cc,
const MoistureComponentIndices &  moisture_indices 
)

Compute turbulent diffusivities using the Mellor-Yamada-Janjic (MYJ) closure model.

Parameters
[in]dtTime step.
[in]xvelX-velocity MultiFab.
[in]yvelY-velocity MultiFab.
[in,out]cons_inConservative variables.
[out]eddyViscosityTurbulent diffusivity MultiFab.
[in]geomGrid geometry.
[in]turbChoiceTurbulence choice parameters.
[in]SurfLayerPointer to surface layer data.
[in]use_terrain_fitted_coordsWhether to use terrain-fitted coordinates.
[in]use_moistureWhether moisture variables are used.
[in]levelLevel index.
[in]bc_ptrBoundary condition pointers.
[in]vert_onlyFlag for vertical-only computation.
[in]z_phys_ndNode-centered physical height field.
[in]z_phys_ccCell-centered physical height field.
[in]moisture_indicesIndices for moisture components.
49 {
50  const StratType pbl_strat = turbChoice.pbl_strat_type;
51  // Dirichlet flags to switch derivative stencil
52  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
53  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].hi(2) == ERFBCType::ext_dir) );
54  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
55  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) );
56  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
57  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) );
58 
59  // Expose constants
60  Real d_kappa = KAPPA;
61 
62  // Closure coefficients (from Janjic (2002), NCEP Office Note 437)
63  Real EPS1 = Real(1.0e-12);
64  Real EPS2 = zero;
65  Real EPSRS = Real(1.0e-7);
66  Real EPSRU = Real(1.0e-7);
67  Real EPSTRB = Real(1.0e-24);
68  Real EPSL = Real(0.32);
69  Real EPSQ2 = Real(0.2);
70  Real EPSQ1 = std::sqrt(EPSQ2);
71  Real ESQHF = Real(2.5);
72  Real FH = Real(1.01);
73 
74  Real G = CONST_GRAV;
75  Real ALPHA = Real(0.3);
76  Real BETA = one/Real(273.);
77  Real EL0MAX = Real(1000.);
78  Real EL0MIN = one;
79  Real ELFC = Real(0.23)*myhalf;
80 
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);
86 
87  Real BTG = BETA*G;
88  Real RB1 = one/B1;
89 
90  Real ADNH = Real(9.)*A1*A2x*A2x*(Real(12.)*A1+three*B2)*BTG*BTG;
91  Real ADNM = Real(18.)*A1*A1*A2x*(B2-three*A2x)*BTG;
92  Real ANMH = -Real(9.)*A1*A2x*A2x*BTG*BTG;
93  Real ANMM = -three*A1*A2x*(three*A2x+three*B2*C1+Real(18.)*A1*C1-B2)*BTG;
94  Real BDNH = three*A2x*(Real(7.)*A1+B2)*BTG;
95  Real BDNM = Real(6.)*A1*A1;
96  Real BEQH = A2x*B1*BTG+three*A2x*(Real(7.)*A1+B2)*BTG;
97  Real BEQM = -A1*B1*(one-three*C1)+Real(6.)*A1*A1;
98  Real BNMH = -A2x*BTG;
99  Real BNMM = A1*(one-three*C1);
100  Real BSHH = Real(9.)*A1*A2x*A2x*BTG;
101  Real BSHM = Real(18.)*A1*A1*A2x*C1;
102  Real BSMH = -three*A1*A2x*(three*A2x+three*B2*C1+Real(12.)*A1*C1-B2)*BTG;
103  Real CESH = A2x;
104  Real CESM = A1*(one-three*C1);
105 
106  Real AEQH = Real(9.)*A1*A2x*A2x*B1*BTG*BTG
107  + Real(9.)*A1*A2x*A2x*(Real(12.)*A1+three*B2)*BTG*BTG;
108  Real AEQM = three*A1*A2x*B1*(three*A2x+three*B2*C1+Real(18.)*A1*C1-B2)*BTG
109  + Real(18.)*A1*A1*A2x*(B2-three*A2x)*BTG;
110 
111  Real REQU = -AEQH/AEQM;
112  Real EPSGH = Real(1.E-9);
113  Real EPSGM = REQU*EPSGH;
114 
115  Real UBRYL = (Real(18.)*REQU*A1*A1*A2x*B2*C1*BTG + Real(9.)*A1*A2x*A2x*B2*BTG*BTG)
116  / (REQU*ADNM+ADNH);
117  Real UBRY = (one+EPSRS)*UBRYL;
118  Real UBRY3 = three*UBRY;
119 
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;
122  Real BUBH = (Real(9.)*A1*A2x+three*A2x*B2)*BTG-BDNH*UBRY3;
123  Real BUBM = Real(18.)*A1*A1*C1-BDNM*UBRY3;
124  Real CUBR = one-UBRY3;
125  Real RCUBR = one/CUBR;
126 
127 #ifdef _OPENMP
128 #pragma omp parallel if (Gpu::notInLaunchRegion())
129 #endif
130  for (MFIter mfi(eddyViscosity,false); mfi.isValid(); ++mfi) {
131 
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);
137 
138  // Ensure the box spans the vertical domain
139  const Box& dbx = geom.Domain();
140  AMREX_ALWAYS_ASSERT(bx.smallEnd(2) == dbx.smallEnd(2) && bx.bigEnd(2) == dbx.bigEnd(2));
141 
142  // Create a plane box
143  int klo = bx.smallEnd(2);
144  int khi = bx.bigEnd(2);
145  Box planexy = makeSlab(bx,2,klo);
146 
147  // Expose for GPU capture
148  const GeometryData gdata = geom.data();
149 
150  // Allocate space for integrals
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();
160 
161  // Terrain and gradient calcs
162  const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
163  const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
164  const auto& dxInv = geom.InvCellSizeArray();
165  int izmin = geom.Domain().smallEnd(2);
166  int izmax = geom.Domain().bigEnd(2);
167 
168  // Surface-layer scales, used to rebuild the gradients in the first cell
169  // above the surface (ERF #4037). MYJ can run without a SurfaceLayer, so
170  // these stay empty when there is none and the resolved stencil is kept.
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>{};
180 
181  // Vertical integrals to compute l0
182  if (use_terrain_fitted_coords) {
183  ParallelFor(planexy, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) noexcept
184  {
185  // Locate PBL k index and set qvel
186  for (int k(klo); k<=khi; ++k) {
187  Real q2 = two * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp);
188  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(q2 > zero, "KE must have a positive value");
189  qvel(i,j,k) = std::sqrt(q2);
190  if (q2<=EPSQ2*FH) {
191  k_arr(i,j,0) = std::min(k,k_arr(i,j,0));
192  }
193  }
194 
195  // Perform integral over PBL height
196  for (int k(klo); k<=k_arr(i,j,0); ++k) {
197  const Real dz = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr);
198  const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr);
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);
201  }
202  });
203  } else {
204  ParallelFor(planexy, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) noexcept
205  {
206  // Locate PBL k index and set qvel
207  for (int k(klo); k<=khi; ++k) {
208  Real q2 = two * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp);
209  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(q2 > zero, "KE must have a positive value");
210  qvel(i,j,k) = std::sqrt(q2);
211  if (q2<=EPSQ2*FH) {
212  k_arr(i,j,0) = std::min(k,k_arr(i,j,0));
213  }
214  }
215 
216  // Perform integral over PBL height
217  for (int k(klo); k<=k_arr(i,j,0); ++k) {
218  // Not multiplying by dz: it's constant and would fall out when we divide qint0/qint1 anyway
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));
222  }
223  });
224  }
225 
226  // Main work to fill diffusivities
227  ParallelFor(planexy, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) noexcept
228  {
229  // Get the PBL k index
230  int kpbl = k_arr(i,j,0);
231 
232  // Compute the integral length scale
233  Real l0 = std::max(std::min(ALPHA*qint(i,j,0,0)/qint(i,j,0,1),EL0MAX),EL0MIN);
234 
235  // Compute diffusivities in each column
236  for (int k(klo); k<=khi; ++k) {
237  // Gradients for shear and buoy production
238  const Real met_h_zeta = use_terrain_fitted_coords ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : one;
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);
247 
248  // Replace the resolved gradients in the first cell with the MOST
249  // profile gradients; see ApplySurfaceLayerGradientsPBL (ERF #4037)
250  if (surface_layer_on_zlo && k == izmin) {
251  const Real zval = use_terrain_fitted_coords ?
252  Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr) :
253  (k + myhalf)*gdata.CellSize(2);
254  const Real rho = cell_data(i,j,k,Rho_comp);
255  const Real theta = cell_data(i,j,k,RhoTheta_comp) / rho;
256  const Real qv = (moisture_indices.qv >= 0) ?
257  cell_data(i,j,k,moisture_indices.qv) / rho : zero;
258  const Real l_min = Real(1.0e-10); // same floor on |L| as MRF and YSUNew
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);
262  sl.tstar_v = ComputeVirtualTStarPBL(t_star_arr(i,j,0),
263  (q_star_arr) ? q_star_arr(i,j,0) : zero,
264  theta, qv, use_moisture);
265  sl.zval = zval;
266  sl.zeta = zval / l_obuk;
267  ApplySurfaceLayerGradientsPBL(sl, dthetavdz, dudz, dvdz);
268  }
269 
270  // Calculate dimensional production terms
271  Real GML = std::max(dudz*dudz + dvdz*dvdz, EPSGM);
272  // NOTE: model uses BTG = beta*g in coeffs above
273  // NOTE: sign convention follows code but theory differs
274  Real GHL = dthetavdz;
275  if (std::fabs(GHL)<=EPSGH) { GHL=EPSGH; }
276 
277  // Find the maximum mixing length
278  Real ELM;
279  if (GHL >= EPSGH) {
280  if (GML/GHL <= REQU) {
281  ELM = EPSL;
282  } else {
283  Real AUBR = (AUBM*GML+AUBH*GHL)*GHL;
284  Real BUBR = BUBM*GML+BUBH*GHL;
285  Real QOL2ST = (-myhalf*BUBR+std::sqrt(BUBR*BUBR*fourth-AUBR*CUBR))*RCUBR;
286  Real ELOQ2X = one/QOL2ST;
287  ELM = std::max(std::sqrt(ELOQ2X*qvel(i,j,k)*qvel(i,j,k)),EPSL);
288  }
289  } else {
290  Real ADEN = (ADNM*GML+ADNH*GHL)*GHL;
291  Real BDEN = BDNM*GML+BDNH*GHL;
292  Real QOL2UN = -myhalf*BDEN+std::sqrt(BDEN*BDEN*fourth-ADEN);
293  Real ELOQ2X = one/(QOL2UN+EPSRU);
294  ELM = std::max(std::sqrt(ELOQ2X*qvel(i,j,k)*qvel(i,j,k)),EPSL);
295  }
296 
297  // Compute master length scale
298  Real L;
299  if (k>kpbl) {
300  L = std::min((met_h_zeta/dxInv[2])*ELFC, ELM);
301  } else {
302  const Real zval = use_terrain_fitted_coords ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr)
303  : (k + myhalf)*gdata.CellSize(2);
304  L = std::min(l0*d_kappa*zval / (d_kappa*zval + l0), ELM);
305  }
306 
307  // Update qvel from production and dissipation
308  Real AEQU = (AEQM*GML+AEQH*GHL)*GHL;
309  Real BEQU = BEQM*GML+BEQH*GHL;
310 
311  Real EQOL2 = -myhalf*BEQU+std::sqrt(BEQU*BEQU*fourth-AEQU);
312 
313  if ( ((GML+GHL*GHL)<=EPSTRB) ||
314  ((GHL>=EPSGH) && ((GML/GHL)<=REQU)) ||
315  (EQOL2<=EPS2) ) {
316  L = EPSL;
317  qvel(i,j,k) = EPSQ1;
318  } else {
319  Real ANUM=(ANMM*GML+ANMH*GHL)*GHL;
320  Real BNUM= BNMM*GML+BNMH*GHL;
321 
322  Real ADEN=(ADNM*GML+ADNH*GHL)*GHL;
323  Real BDEN= BDNM*GML+BDNH*GHL;
324  Real CDEN= one;
325 
326  Real ARHS=-(ANUM*BDEN-BNUM*ADEN)*two;
327  Real BRHS=- ANUM*Real(4.);
328  Real CRHS=- BNUM*two;
329 
330  Real DLOQ1=L/qvel(i,j,k);
331 
332  Real ELOQ21=one/EQOL2;
333  Real ELOQ11=std::sqrt(ELOQ21);
334  Real ELOQ31=ELOQ21*ELOQ11;
335  Real ELOQ41=ELOQ21*ELOQ21;
336  Real ELOQ51=ELOQ21*ELOQ31;
337 
338  Real RDEN1=one/(ADEN*ELOQ41+BDEN*ELOQ21+CDEN);
339 
340  Real RHSP1=(ARHS*ELOQ51+BRHS*ELOQ31+CRHS*ELOQ11)*RDEN1*RDEN1;
341 
342  Real DTTURBL = static_cast<Real>(dt);
343  Real ELOQ12 = amrex::max(ELOQ11+(DLOQ1-ELOQ11)*std::exp(RHSP1*DTTURBL),EPS1);
344 
345  Real ELOQ22=ELOQ12*ELOQ12;
346  Real ELOQ32=ELOQ22*ELOQ12;
347  Real ELOQ42=ELOQ22*ELOQ22;
348  Real ELOQ52=ELOQ22*ELOQ32;
349 
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;
354 
355  Real ELOQ13 = amrex::max(ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)*std::exp(RHSP2*DTTURBL),EPS1);
356 
357  Real ELOQN=ELOQ13;
358  if (ELOQN>EPS1) {
359  qvel(i,j,k) = std::max(L/ELOQN,EPSQ1);
360  if (qvel(i,j,k)==EPSQ1) {L = EPSL; }
361  } else {
362  L = EPSL;
363  qvel(i,j,k) = EPSQ1;
364  }
365  }
366 
367  // NOTE: MYJ in WRF places the TKE at the w-faces and not the CC.
368  // Therefore, WRF sets the bottom face in accordance with
369  // the surface layer friction velocity (below). We do not
370  // do that since TKE lives at CC here.
371  //if (k==klo) {
372  // Real q2 = std::pow(B1,(two/three))*ustar_arr(i,j,k)*ustar_arr(i,j,k);
373  // Real q = std::max(std::sqrt(q2),EPSQ1);
374  // qvel(i,j,k) = myhalf * (q + qvel(i,j,k));
375  //}
376 
377  cell_data(i,j,k,RhoKE_comp) = myhalf*cell_data(i,j,k,Rho_comp)*qvel(i,j,k)*qvel(i,j,k);
378 
379  // L^n/Q^n
380  Real ELOQ2 = L*L/(qvel(i,j,k)*qvel(i,j,k));
381  Real ELOQ4 = ELOQ2*ELOQ2;
382 
383  // COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
384  Real ADEN=(ADNM*GML+ADNH*GHL)*GHL;
385  Real BDEN= BDNM*GML+BDNH*GHL;
386  Real CDEN= one;
387 
388  // COEFFICIENTS FOR THE SM DETERMINANT
389  Real BESM=BSMH*GHL;
390 
391  // COEFFICIENTS FOR THE SH DETERMINANT
392  Real BESH=BSHM*GML+BSHH*GHL;
393 
394  // 1/DENOMINATOR
395  Real RDEN=one/(ADEN*ELOQ4+BDEN*ELOQ2+CDEN);
396 
397  // SM, SH, SQ
398  Real SM=(BESM*ELOQ2+CESM)*RDEN;
399  Real SH=(BESH*ELOQ2+CESH)*RDEN;
400  Real SQ=ESQHF*SH;
401 
402  // Finally, compute the eddy viscosity/diffusivities
403  const Real rho = cell_data(i,j,k,Rho_comp);
404  K_turb(i,j,k,EddyDiff::Mom_v ) = rho * L * qvel(i,j,k) * SM;
405  K_turb(i,j,k,EddyDiff::Theta_v) = rho * L * qvel(i,j,k) * SH;
406  K_turb(i,j,k,EddyDiff::KE_v ) = rho * L * qvel(i,j,k) * SQ;
407  K_turb(i,j,k,EddyDiff::Q_v ) = rho * L * qvel(i,j,k) * SH;
408  K_turb(i,j,k,EddyDiff::Turb_lengthscale) = L;
409 
410  // NOTE: Ghost cells are handled at end of ERF_ComputeTurbulentViscosity.cpp
411 
412  } // for k
413  }); // ParFor
414  } // mfi
415 }
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

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:
Here is the caller graph for this function: