ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ComputeTurbulentViscosity.cpp File Reference
#include <ABLMost.H>
#include <EddyViscosity.H>
#include <Diffusion.H>
#include <TileNoZ.H>
#include <TerrainMetrics.H>
Include dependency graph for ComputeTurbulentViscosity.cpp:

Functions

void ComputeTurbulentViscosityPBL (const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< ABLMost > &most, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd)
 
void ComputeTurbulentViscosityLES (const MultiFab &Tau11, const MultiFab &Tau22, const MultiFab &Tau33, const MultiFab &Tau12, const MultiFab &Tau13, const MultiFab &Tau23, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, const MultiFab &mapfac_u, const MultiFab &mapfac_v, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< ABLMost > &most, const bool &exp_most)
 
void ComputeTurbulentViscosity (const MultiFab &xvel, const MultiFab &yvel, const MultiFab &Tau11, const MultiFab &Tau22, const MultiFab &Tau33, const MultiFab &Tau12, const MultiFab &Tau13, const MultiFab &Tau23, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, const MultiFab &mapfac_u, const MultiFab &mapfac_v, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< ABLMost > &most, const bool &exp_most, int level, const BCRec *bc_ptr, bool vert_only)
 

Function Documentation

◆ ComputeTurbulentViscosity()

void ComputeTurbulentViscosity ( const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  Tau11,
const MultiFab &  Tau22,
const MultiFab &  Tau33,
const MultiFab &  Tau12,
const MultiFab &  Tau13,
const MultiFab &  Tau23,
const MultiFab &  cons_in,
MultiFab &  eddyViscosity,
MultiFab &  Hfx1,
MultiFab &  Hfx2,
MultiFab &  Hfx3,
MultiFab &  Diss,
const Geometry &  geom,
const MultiFab &  mapfac_u,
const MultiFab &  mapfac_v,
const std::unique_ptr< MultiFab > &  z_phys_nd,
const TurbChoice &  turbChoice,
const Real  const_grav,
std::unique_ptr< ABLMost > &  most,
const bool &  exp_most,
int  level,
const BCRec *  bc_ptr,
bool  vert_only 
)

Wrapper to compute turbulent viscosity with LES or PBL.

Parameters
[in]xvelvelocity in x-dir
[in]yvelvelocity in y-dir
[in]Tau1111 strain
[in]Tau2222 strain
[in]Tau3333 strain
[in]Tau1212 strain
[in]Tau1313 strain
[in]Tau2323 strain
[in]cons_incell center conserved quantities
[out]eddyViscosityturbulent viscosity
[in]Hfx1heat flux in x-dir
[in]Hfx2heat flux in y-dir
[in]Hfx3heat flux in z-dir
[in]Dissdissipation of turbulent kinetic energy
[in]geomproblem geometry
[in]mapfac_umap factor at x-face
[in]mapfac_vmap factor at y-face
[in]turbChoicecontainer with turbulence parameters
[in]mostpointer to Monin-Obukhov class if instantiated
[in]vert_onlyflag for vertical components of eddyViscosity
421 {
422  BL_PROFILE_VAR("ComputeTurbulentViscosity()",ComputeTurbulentViscosity);
423  //
424  // In LES mode, the turbulent viscosity is isotropic, so the LES model sets both horizontal and vertical viscosities
425  // In PBL mode, the primary purpose of the PBL model is to control vertical transport, so the PBL model sets the vertical viscosity.
426  // Optionally, the PBL model can be run in conjunction with an LES model that sets the horizontal viscosity
427  // (this isn’t truly LES, but the model form is the same as Smagorinsky).
428  //
429  // ComputeTurbulentViscosityLES populates the LES viscosity for both horizontal and vertical components.
430  // ComputeTurbulentViscosityPBL computes the PBL viscosity just for the vertical component.
431  //
432 
433  if (most) {
434  bool l_use_turb = ( turbChoice.les_type == LESType::Smagorinsky ||
435  turbChoice.les_type == LESType::Deardorff ||
436  turbChoice.pbl_type == PBLType::MYNN25 ||
437  turbChoice.pbl_type == PBLType::YSU );
438  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l_use_turb,
439  "An LES or PBL model must be utilized with MOST boundaries to compute the turbulent viscosity");
440  } else {
441  AMREX_ALWAYS_ASSERT(!vert_only);
442  }
443 
444  if (turbChoice.les_type != LESType::None) {
445  ComputeTurbulentViscosityLES(Tau11, Tau22, Tau33,
446  Tau12, Tau13, Tau23,
447  cons_in, eddyViscosity,
448  Hfx1, Hfx2, Hfx3, Diss,
449  geom, mapfac_u, mapfac_v,
450  z_phys_nd, turbChoice, const_grav,
451  most, exp_most);
452  }
453 
454  if (turbChoice.pbl_type != PBLType::None) {
455  // NOTE: state_new is passed in for Cons_old (due to ptr swap in advance)
456  ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity,
457  geom, turbChoice, most, level, bc_ptr, vert_only, z_phys_nd);
458  }
459 }
ComputeTurbulentViscosityLES
void ComputeTurbulentViscosityLES(const MultiFab &Tau11, const MultiFab &Tau22, const MultiFab &Tau33, const MultiFab &Tau12, const MultiFab &Tau13, const MultiFab &Tau23, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, const MultiFab &mapfac_u, const MultiFab &mapfac_v, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< ABLMost > &most, const bool &exp_most)
Definition: ComputeTurbulentViscosity.cpp:44
ComputeTurbulentViscosity
void ComputeTurbulentViscosity(const MultiFab &xvel, const MultiFab &yvel, const MultiFab &Tau11, const MultiFab &Tau22, const MultiFab &Tau33, const MultiFab &Tau12, const MultiFab &Tau13, const MultiFab &Tau23, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, const MultiFab &mapfac_u, const MultiFab &mapfac_v, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< ABLMost > &most, const bool &exp_most, int level, const BCRec *bc_ptr, bool vert_only)
Definition: ComputeTurbulentViscosity.cpp:406
ComputeTurbulentViscosityPBL
void ComputeTurbulentViscosityPBL(const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< ABLMost > &most, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd)
Definition: PBLModels.cpp:22
Vars::xvel
@ xvel
Definition: IndexDefines.H:100
Vars::yvel
@ yvel
Definition: IndexDefines.H:101
TurbChoice::pbl_type
PBLType pbl_type
Definition: TurbStruct.H:279
TurbChoice::les_type
LESType les_type
Definition: TurbStruct.H:256

Referenced by ERF::advance_dycore().

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

◆ ComputeTurbulentViscosityLES()

void ComputeTurbulentViscosityLES ( const MultiFab &  Tau11,
const MultiFab &  Tau22,
const MultiFab &  Tau33,
const MultiFab &  Tau12,
const MultiFab &  Tau13,
const MultiFab &  Tau23,
const MultiFab &  cons_in,
MultiFab &  eddyViscosity,
MultiFab &  Hfx1,
MultiFab &  Hfx2,
MultiFab &  Hfx3,
MultiFab &  Diss,
const Geometry &  geom,
const MultiFab &  mapfac_u,
const MultiFab &  mapfac_v,
const std::unique_ptr< MultiFab > &  z_phys_nd,
const TurbChoice &  turbChoice,
const Real  const_grav,
std::unique_ptr< ABLMost > &  most,
const bool &  exp_most 
)

Function for computing the turbulent viscosity with LES.

Parameters
[in]Tau1111 strain
[in]Tau2222 strain
[in]Tau3333 strain
[in]Tau1212 strain
[in]Tau1313 strain
[in]Tau2323 strain
[in]cons_incell center conserved quantities
[out]eddyViscosityturbulent viscosity
[in]Hfx1heat flux in x-dir
[in]Hfx2heat flux in y-dir
[in]Hfx3heat flux in z-dir
[in]Dissdissipation of turbulent kinetic energy
[in]geomproblem geometry
[in]mapfac_umap factor at x-face
[in]mapfac_vmap factor at y-face
[in]turbChoicecontainer with turbulence parameters
53 {
54  const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
55  const Box& domain = geom.Domain();
56  const int& klo = domain.smallEnd(2);
57  const bool use_most = (most != nullptr);
58  const bool use_terrain = (z_phys_nd != nullptr);
59 
60  // SMAGORINSKY: Fill Kturb for momentum in horizontal and vertical
61  //***********************************************************************************
62  if (turbChoice.les_type == LESType::Smagorinsky)
63  {
64  Real Cs = turbChoice.Cs;
65 
66 #ifdef _OPENMP
67 #pragma omp parallel if (Gpu::notInLaunchRegion())
68 #endif
69  for (MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
70  {
71  // NOTE: This gets us the lateral ghost cells for lev>0; which
72  // have been filled from FP Two Levels.
73  Box bxcc = mfi.growntilebox() & domain;
74 
75  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
76  const Array4<Real const > &cell_data = cons_in.array(mfi);
77 
78  Array4<Real const> tau11 = Tau11.array(mfi);
79  Array4<Real const> tau22 = Tau22.array(mfi);
80  Array4<Real const> tau33 = Tau33.array(mfi);
81  Array4<Real const> tau12 = Tau12.array(mfi);
82  Array4<Real const> tau13 = Tau13.array(mfi);
83  Array4<Real const> tau23 = Tau23.array(mfi);
84 
85  Array4<Real const> mf_u = mapfac_u.array(mfi);
86  Array4<Real const> mf_v = mapfac_v.array(mfi);
87 
88  Array4<Real const> z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4<Real const>{};
89 
90  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
91  {
92  Real SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,klo,use_most,exp_most);
93  Real dxInv = cellSizeInv[0];
94  Real dyInv = cellSizeInv[1];
95  Real dzInv = cellSizeInv[2];
96  if (use_terrain) {
97  // the terrain grid is only deformed in z for now
98  dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
99  }
100  Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
101  Real DeltaMsf = std::pow(cellVolMsf,1.0/3.0);
102  Real CsDeltaSqrMsf = Cs*Cs*DeltaMsf*DeltaMsf;
103 
104  mu_turb(i, j, k, EddyDiff::Mom_h) = CsDeltaSqrMsf * cell_data(i, j, k, Rho_comp) * std::sqrt(2.0*SmnSmn);
105  mu_turb(i, j, k, EddyDiff::Mom_v) = mu_turb(i, j, k, EddyDiff::Mom_h);
106  });
107  }
108  }
109  // DEARDORFF: Fill Kturb for momentum in horizontal and vertical
110  //***********************************************************************************
111  else if (turbChoice.les_type == LESType::Deardorff)
112  {
113  const Real l_C_k = turbChoice.Ck;
114  const Real l_C_e = turbChoice.Ce;
115  const Real l_C_e_wall = turbChoice.Ce_wall;
116  const Real Ce_lcoeff = amrex::max(0.0, l_C_e - 1.9*l_C_k);
117  const Real l_abs_g = const_grav;
118  const Real l_inv_theta0 = 1.0 / turbChoice.theta_ref;
119 
120 #ifdef _OPENMP
121 #pragma omp parallel if (Gpu::notInLaunchRegion())
122 #endif
123  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
124  {
125  Box bxcc = mfi.tilebox();
126 
127  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
128  const Array4<Real>& hfx_x = Hfx1.array(mfi);
129  const Array4<Real>& hfx_y = Hfx2.array(mfi);
130  const Array4<Real>& hfx_z = Hfx3.array(mfi);
131  const Array4<Real>& diss = Diss.array(mfi);
132 
133  const Array4<Real const > &cell_data = cons_in.array(mfi);
134 
135  Array4<Real const> mf_u = mapfac_u.array(mfi);
136  Array4<Real const> mf_v = mapfac_v.array(mfi);
137 
138  Array4<Real const> z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4<Real const>{};
139 
140  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
141  {
142  Real dxInv = cellSizeInv[0];
143  Real dyInv = cellSizeInv[1];
144  Real dzInv = cellSizeInv[2];
145  if (use_terrain) {
146  // the terrain grid is only deformed in z for now
147  dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
148  }
149  Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
150  Real DeltaMsf = std::pow(cellVolMsf,1.0/3.0);
151 
152  // Calculate stratification-dependent mixing length (Deardorff 1980)
153  Real eps = std::numeric_limits<Real>::epsilon();
154  Real dtheta_dz;
155  if (use_most && k==klo) {
156  if (exp_most) {
157  dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
158  - cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) )*dzInv;
159  } else {
160  dtheta_dz = 0.5 * (-3 * cell_data(i,j,k ,RhoTheta_comp)
161  / cell_data(i,j,k ,Rho_comp)
162  + 4 * cell_data(i,j,k+1,RhoTheta_comp)
163  / cell_data(i,j,k+1,Rho_comp)
164  - cell_data(i,j,k+2,RhoTheta_comp)
165  / cell_data(i,j,k+2,Rho_comp) ) * dzInv;
166  }
167  } else {
168  dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
169  - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv;
170  }
171  Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp);
172  Real strat = l_abs_g * dtheta_dz * l_inv_theta0; // stratification
173  Real length;
174  if (strat <= eps) {
175  length = DeltaMsf;
176  } else {
177  length = 0.76 * std::sqrt(E / strat);
178  // mixing length should be _reduced_ for stable stratification
179  length = amrex::min(length, DeltaMsf);
180  // following WRF, make sure the mixing length isn't too small
181  length = amrex::max(length, 0.001 * DeltaMsf);
182  }
183 
184  // Calculate eddy diffusivities
185  // K = rho * C_k * l * KE^(1/2)
186  mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E);
187  mu_turb(i,j,k,EddyDiff::Mom_v) = mu_turb(i,j,k,EddyDiff::Mom_h);
188  // KH = (1 + 2*l/delta) * mu_turb
189  mu_turb(i,j,k,EddyDiff::Theta_v) = (1.+2.*length/DeltaMsf) * mu_turb(i,j,k,EddyDiff::Mom_v);
190 
191  // Calculate SFS quantities
192  // - dissipation
193  Real Ce;
194  if ((l_C_e_wall > 0) && (k==0))
195  Ce = l_C_e_wall;
196  else
197  Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf;
198  diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length;
199  // - heat flux
200  // (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will
201  // be overwritten when BCs are applied)
202  hfx_x(i,j,k) = 0.0;
203  hfx_y(i,j,k) = 0.0;
204  hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K]
205  });
206  }
207  }
208 
209  // Extrapolate Kturb in x/y, fill remaining elements (relevent to lev==0)
210  //***********************************************************************************
211  int ngc(1);
212  Real inv_Pr_t = turbChoice.Pr_t_inv;
213  Real inv_Sc_t = turbChoice.Sc_t_inv;
214  Real inv_sigma_k = 1.0 / turbChoice.sigma_k;
215  // EddyDiff mapping : Theta_h KE_h QKE_h Scalar_h Q_h
216  Vector<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
217  Gpu::AsyncVector<Real> d_Factors; d_Factors.resize(Factors.size());
218  Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin());
219  Real* fac_ptr = d_Factors.data();
220 
221  bool use_KE = (turbChoice.les_type == LESType::Deardorff);
222  bool use_QKE = (turbChoice.use_QKE && turbChoice.diffuse_QKE_3D);
223 
224 #ifdef _OPENMP
225 #pragma omp parallel if (Gpu::notInLaunchRegion())
226 #endif
227  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
228  {
229  Box bxcc = mfi.tilebox();
230  Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1);
231  Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); planey.grow(0,1);
232  int i_lo = bxcc.smallEnd(0); int i_hi = bxcc.bigEnd(0);
233  int j_lo = bxcc.smallEnd(1); int j_hi = bxcc.bigEnd(1);
234  bxcc.growLo(0,ngc); bxcc.growHi(0,ngc);
235  bxcc.growLo(1,ngc); bxcc.growHi(1,ngc);
236 
237  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
238 
239  // Extrapolate outside the domain in lateral directions (planex owns corner cells)
240  if (i_lo == domain.smallEnd(0)) {
241  ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
242  {
243  int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1));
244  mu_turb(i_lo-i, j, k, EddyDiff::Mom_h) = mu_turb(i_lo, lj, k, EddyDiff::Mom_h);
245  mu_turb(i_lo-i, j, k, EddyDiff::Mom_v) = mu_turb(i_lo, lj, k, EddyDiff::Mom_v);
246  });
247  }
248  if (i_hi == domain.bigEnd(0)) {
249  ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
250  {
251  int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1));
252  mu_turb(i_hi+i, j, k, EddyDiff::Mom_h) = mu_turb(i_hi, lj, k, EddyDiff::Mom_h);
253  mu_turb(i_hi+i, j, k, EddyDiff::Mom_v) = mu_turb(i_hi, lj, k, EddyDiff::Mom_v);
254  });
255  }
256  if (j_lo == domain.smallEnd(1)) {
257  ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
258  {
259  int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0));
260  mu_turb(i, j_lo-j, k, EddyDiff::Mom_h) = mu_turb(li, j_lo, k, EddyDiff::Mom_h);
261  mu_turb(i, j_lo-j, k, EddyDiff::Mom_v) = mu_turb(li, j_lo, k, EddyDiff::Mom_v);
262  });
263  }
264  if (j_hi == domain.bigEnd(1)) {
265  ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
266  {
267  int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0));
268  mu_turb(i, j_hi+j, k, EddyDiff::Mom_h) = mu_turb(li, j_hi, k, EddyDiff::Mom_h);
269  mu_turb(i, j_hi+j, k, EddyDiff::Mom_v) = mu_turb(li, j_hi, k, EddyDiff::Mom_v);
270  });
271  }
272 
273  // refactor the code to eliminate the need for ifdef's
274  for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) {
275  int offset = (EddyDiff::NumDiffs-1)/2;
276  switch (n)
277  {
278  case EddyDiff::QKE_h:
279  // Populate element other than mom_h/v on the whole grid
280  if(use_QKE) {
281  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
282  {
283  int indx = n;
284  int indx_v = indx + offset;
285  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
286  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
287  });
288  }
289  break;
290  case EddyDiff::KE_h:
291  if (use_KE) {
292  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
293  {
294  int indx = n;
295  int indx_v = indx + offset;
296  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
297  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
298  });
299  }
300  break;
301  default:
302  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
303  {
304  int indx = n;
305  int indx_v = indx + offset;
306  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
307  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
308  });
309  break;
310  }
311  }
312  }
313 
314  // Fill interior ghost cells and any ghost cells outside a periodic domain
315  //***********************************************************************************
316  eddyViscosity.FillBoundary(geom.periodicity());
317 
318 
319  // Extrapolate top & bottom
320  //***********************************************************************************
321 #ifdef _OPENMP
322 #pragma omp parallel if (Gpu::notInLaunchRegion())
323 #endif
324  for ( MFIter mfi(eddyViscosity,TileNoZ()); mfi.isValid(); ++mfi)
325  {
326  Box bxcc = mfi.tilebox();
327  Box planez = bxcc; planez.setSmall(2, 1); planez.setBig(2, ngc);
328  int k_lo = bxcc.smallEnd(2); int k_hi = bxcc.bigEnd(2);
329  planez.growLo(0,ngc); planez.growHi(0,ngc);
330  planez.growLo(1,ngc); planez.growHi(1,ngc);
331 
332  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
333 
334  // refactor the code to eliminate the need for ifdef's
335  for (auto n = 0; n < (EddyDiff::NumDiffs-1)/2; ++n) {
336  int offset = (EddyDiff::NumDiffs-1)/2;
337  switch (n)
338  {
339  case EddyDiff::QKE_h:
340  // Extrap all components at top & bottom
341  if(use_QKE) {
342  ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
343  {
344  int indx = n;
345  int indx_v = indx + offset;
346  mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx );
347  mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx );
348  mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v);
349  mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v);
350  });
351  }
352  break;
353  case EddyDiff::KE_h:
354  if (use_KE) {
355  ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
356  {
357  int indx = n;
358  int indx_v = indx + offset;
359  mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx );
360  mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx );
361  mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v);
362  mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v);
363  });
364  }
365  break;
366  default:
367  ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
368  {
369  int indx = n ;
370  int indx_v = indx + offset;
371  mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx );
372  mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx );
373  mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v);
374  mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v);
375  });
376  break;
377  }
378  }
379  }
380 }
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
ComputeSmnSmn
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeSmnSmn(int &i, int &j, int &k, const amrex::Array4< amrex::Real const > &tau11, const amrex::Array4< amrex::Real const > &tau22, const amrex::Array4< amrex::Real const > &tau33, const amrex::Array4< amrex::Real const > &tau12, const amrex::Array4< amrex::Real const > &tau13, const amrex::Array4< amrex::Real const > &tau23, const int &klo, const bool &use_most, const bool &exp_most)
Definition: EddyViscosity.H:30
Rho_comp
#define Rho_comp
Definition: IndexDefines.H:12
RhoTheta_comp
#define RhoTheta_comp
Definition: IndexDefines.H:13
RhoKE_comp
#define RhoKE_comp
Definition: IndexDefines.H:14
Compute_h_zeta_AtCellCenter
AMREX_FORCE_INLINE AMREX_GPU_DEVICE 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: TerrainMetrics.H:31
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: TileNoZ.H:11
EddyDiff::Theta_v
@ Theta_v
Definition: IndexDefines.H:128
EddyDiff::Mom_h
@ Mom_h
Definition: IndexDefines.H:121
EddyDiff::Mom_v
@ Mom_v
Definition: IndexDefines.H:127
EddyDiff::NumDiffs
@ NumDiffs
Definition: IndexDefines.H:134
EddyDiff::KE_h
@ KE_h
Definition: IndexDefines.H:123
EddyDiff::QKE_h
@ QKE_h
Definition: IndexDefines.H:124
TurbChoice::sigma_k
amrex::Real sigma_k
Definition: TurbStruct.H:274
TurbChoice::Sc_t_inv
amrex::Real Sc_t_inv
Definition: TurbStruct.H:266
TurbChoice::Ck
amrex::Real Ck
Definition: TurbStruct.H:272
TurbChoice::use_QKE
bool use_QKE
Definition: TurbStruct.H:299
TurbChoice::Cs
amrex::Real Cs
Definition: TurbStruct.H:258
TurbChoice::Pr_t_inv
amrex::Real Pr_t_inv
Definition: TurbStruct.H:263
TurbChoice::Ce_wall
amrex::Real Ce_wall
Definition: TurbStruct.H:270
TurbChoice::Ce
amrex::Real Ce
Definition: TurbStruct.H:269
TurbChoice::diffuse_QKE_3D
bool diffuse_QKE_3D
Definition: TurbStruct.H:300
TurbChoice::theta_ref
amrex::Real theta_ref
Definition: TurbStruct.H:276

Referenced by ComputeTurbulentViscosity().

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

◆ ComputeTurbulentViscosityPBL()

void ComputeTurbulentViscosityPBL ( const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  cons_in,
MultiFab &  eddyViscosity,
const Geometry &  geom,
const TurbChoice &  turbChoice,
std::unique_ptr< ABLMost > &  most,
int  level,
const BCRec *  bc_ptr,
bool  ,
const std::unique_ptr< MultiFab > &  z_phys_nd 
)

Function to compute turbulent viscosity with PBL.

Parameters
[in]xvelvelocity in x-dir
[in]yvelvelocity in y-dir
[in]cons_incell center conserved quantities
[out]eddyViscosityholds turbulent viscosity
[in]geomproblem geometry
[in]turbChoicecontainer with turbulence parameters
[in]mostpointer to Monin-Obukhov class if instantiated
33 {
34  const bool use_terrain = (z_phys_nd != nullptr);
35 
36  // MYNN Level 2.5 PBL Model
37  if (turbChoice.pbl_type == PBLType::MYNN25) {
38 
39  const Real A1 = turbChoice.pbl_mynn_A1;
40  const Real A2 = turbChoice.pbl_mynn_A2;
41  const Real B1 = turbChoice.pbl_mynn_B1;
42  const Real B2 = turbChoice.pbl_mynn_B2;
43  const Real C1 = turbChoice.pbl_mynn_C1;
44  const Real C2 = turbChoice.pbl_mynn_C2;
45  const Real C3 = turbChoice.pbl_mynn_C3;
46  //const Real C4 = turbChoice.pbl_mynn_C4;
47  const Real C5 = turbChoice.pbl_mynn_C5;
48 
49  // Dirichlet flags to switch derivative stencil
50  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
51  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) );
52  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
53  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) );
54  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
55  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) );
56 
57  // Epsilon
58  Real eps = std::numeric_limits<Real>::epsilon();
59 
60 #ifdef _OPENMP
61 #pragma omp parallel if (Gpu::notInLaunchRegion())
62 #endif
63  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
64 
65  const Box &bx = mfi.growntilebox(1);
66  const Array4<Real const> &cell_data = cons_in.array(mfi);
67  const Array4<Real > &K_turb = eddyViscosity.array(mfi);
68  const Array4<Real const> &uvel = xvel.array(mfi);
69  const Array4<Real const> &vvel = yvel.array(mfi);
70 
71  // Compute some quantities that are constant in each column
72  // Sbox is shrunk to only include the interior of the domain in the vertical direction to compute integrals
73  // Box includes one ghost cell in each direction
74  const Box &dbx = geom.Domain();
75  Box sbx(bx.smallEnd(), bx.bigEnd());
76  sbx.grow(2,-1);
77  AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
78 
79  const GeometryData gdata = geom.data();
80 
81  const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
82  FArrayBox qintegral(xybx,2);
83  qintegral.setVal<RunOn::Device>(0.0);
84  FArrayBox qturb(bx,1); FArrayBox qturb_old(bx,1);
85  const Array4<Real> qint = qintegral.array();
86  const Array4<Real> qvel = qturb.array();
87  const Array4<Real> qvel_old = qturb_old.array();
88 
89  // vertical integrals to compute lengthscale
90  if (use_terrain) {
91  const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
92  const auto invCellSize = geom.InvCellSizeArray();
93  ParallelFor(Gpu::KernelInfo().setReduction(true), bx,
94  [=] AMREX_GPU_DEVICE (int i, int j, int k, Gpu::Handler const& handler) noexcept
95  {
96  qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
97  qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps);
98  AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value");
99  AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value");
100 
101  if (sbx.contains(i,j,k)) {
102  const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr);
103  const Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd_arr);
104  Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz, handler);
105  Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k)*dz, handler);
106  }
107  });
108  } else {
109  ParallelFor(Gpu::KernelInfo().setReduction(true), bx,
110  [=] AMREX_GPU_DEVICE (int i, int j, int k, Gpu::Handler const& handler) noexcept
111  {
112  qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
113  qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps);
114  AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value");
115  AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value");
116 
117  // Not multiplying by dz: its constant and would fall out when we divide qint0/qint1 anyway
118  if (sbx.contains(i,j,k)) {
119  const Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
120  Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k), handler);
121  Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k), handler);
122  }
123  });
124  }
125 
126  Real dz_inv = geom.InvCellSize(2);
127  const auto& dxInv = geom.InvCellSizeArray();
128  int izmin = geom.Domain().smallEnd(2);
129  int izmax = geom.Domain().bigEnd(2);
130 
131  // Spatially varying MOST
132  Real d_kappa = KAPPA;
133  Real d_gravity = CONST_GRAV;
134 
135  const auto& t_mean_mf = most->get_mac_avg(level,2); // TODO: IS THIS ACTUALLY RHOTHETA
136  const auto& u_star_mf = most->get_u_star(level); // Use desired level
137  const auto& t_star_mf = most->get_t_star(level); // Use desired level
138 
139  const auto& tm_arr = t_mean_mf->array(mfi); // TODO: IS THIS ACTUALLY RHOTHETA
140  const auto& u_star_arr = u_star_mf->array(mfi);
141  const auto& t_star_arr = t_star_mf->array(mfi);
142 
143  const Array4<Real const> z_nd_arr = use_terrain ? z_phys_nd->array(mfi) : Array4<Real>{};
144 
145  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
146  {
147  // Compute some partial derivatives that we will need (second order)
148  // U and V derivatives are interpolated to account for staggered grid
149  const Real met_h_zeta = use_terrain ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : 1.0;
150  Real dthetadz, dudz, dvdz;
151  ComputeVerticalDerivativesPBL(i, j, k,
152  uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta,
153  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
154  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
155  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
156  dthetadz, dudz, dvdz);
157 
158  // Spatially varying MOST
159  Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
160  Real theta0 = tm_arr(i,j,0); // TODO: IS THIS ACTUALLY RHOTHETA
161  Real l_obukhov;
162  if (std::abs(surface_heat_flux) > eps) {
163  l_obukhov = ( theta0 * u_star_arr(i,j,0) * u_star_arr(i,j,0) ) /
164  ( d_kappa * d_gravity * t_star_arr(i,j,0) );
165  } else {
166  l_obukhov = std::numeric_limits<Real>::max();
167  }
168 
169  // First Length Scale
170  AMREX_ASSERT(l_obukhov != 0);
171  int lk = amrex::max(k,0);
172  const Real zval = use_terrain ? Compute_Zrel_AtCellCenter(i,j,lk,z_nd_arr) : gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
173  const Real zeta = zval/l_obukhov;
174  Real l_S;
175  if (zeta >= 1.0) {
176  l_S = KAPPA*zval/3.7;
177  } else if (zeta >= 0) {
178  l_S = KAPPA*zval/(1+2.7*zeta);
179  } else {
180  l_S = KAPPA*zval*std::pow(1.0 - 100.0 * zeta, 0.2);
181  }
182 
183  // Second Length Scale
184  Real l_T;
185  if (qint(i,j,0,1) > 0.0) {
186  l_T = 0.23*qint(i,j,0,0)/qint(i,j,0,1);
187  } else {
188  l_T = std::numeric_limits<Real>::max();
189  }
190 
191  // Third Length Scale
192  Real l_B;
193  if (dthetadz > 0) {
194  Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetadz);
195  if (zeta < 0) {
196  Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T;
197  qc = std::pow(qc,1.0/3.0);
198  l_B = (1.0 + 5.0*std::sqrt(qc/(N_brunt_vaisala * l_T))) * qvel(i,j,k)/N_brunt_vaisala;
199  } else {
200  l_B = qvel(i,j,k) / N_brunt_vaisala;
201  }
202  } else {
203  l_B = std::numeric_limits<Real>::max();
204  }
205 
206  // Overall Length Scale
207  Real l_comb = 1.0 / (1.0/l_S + 1.0/l_T + 1.0/l_B);
208 
209  // NOTE: Level 2 limiting from balance of production and dissipation.
210  // K_turb has a setval of 0.0 when the MF is created (NOT EACH STEP).
211  // We do this inline to avoid storing qe^2 at each cell.
212  Real l_comb_old = K_turb(i,j,k,EddyDiff::PBL_lengthscale);
213  Real shearProd = dudz*dudz + dvdz*dvdz;
214  Real buoyProd = -(CONST_GRAV/theta0) * dthetadz;
215  Real lSM = K_turb(i,j,k,EddyDiff::Mom_v) / (qvel_old(i,j,k) + eps);
216  Real lSH = K_turb(i,j,k,EddyDiff::Theta_v) / (qvel_old(i,j,k) + eps);
217  Real qe2 = B1 * l_comb_old * ( lSM * shearProd + lSH * buoyProd );
218  Real qe = (qe2 < 0.0) ? 0.0 : std::sqrt(qe2);
219  Real one_m_alpha = (qvel(i,j,k) > qe) ? 1.0 : qvel(i,j,k) / (qe + eps);
220  Real one_m_alpha2 = one_m_alpha * one_m_alpha;
221 
222  // Compute non-dimensional parameters
223  Real l2_over_q2 = l_comb*l_comb/(qvel(i,j,k)*qvel(i,j,k));
224  Real GM = l2_over_q2 * shearProd;
225  Real GH = l2_over_q2 * buoyProd;
226  Real E1 = 1.0 + one_m_alpha2 * ( 6.0*A1*A1*GM - 9.0*A1*A2*(1.0-C2)*GH );
227  Real E2 = one_m_alpha2 * ( -3.0*A1*(4.0*A1 + 3.0*A2*(1.0-C5))*(1.0-C2)*GH );
228  Real E3 = one_m_alpha2 * ( 6.0*A1*A2*GM );
229  Real E4 = 1.0 + one_m_alpha2 * ( -12.0*A1*A2*(1.0-C2)*GH - 3.0*A2*B2*(1.0-C3)*GH );
230  Real R1 = one_m_alpha * ( A1*(1.0-3.0*C1) );
231  Real R2 = one_m_alpha * A2;
232 
233  Real SM = (R2*E2 - R1*E4)/(E2*E3 - E1*E4);
234  Real SH = (R1*E3 - R2*E1)/(E2*E3 - E1*E4);
235  Real SQ = 3.0 * SM; // Nakanishi & Niino 2009
236 
237  // Finally, compute the eddy viscosity/diffusivities
238  const Real rho = cell_data(i,j,k,Rho_comp);
239  K_turb(i,j,k,EddyDiff::Mom_v) = rho * l_comb * qvel(i,j,k) * SM * 0.5; // 0.5 for mu_turb
240  K_turb(i,j,k,EddyDiff::Theta_v) = rho * l_comb * qvel(i,j,k) * SH;
241  K_turb(i,j,k,EddyDiff::QKE_v) = rho * l_comb * qvel(i,j,k) * SQ;
242 
243  K_turb(i,j,k,EddyDiff::PBL_lengthscale) = l_comb;
244  // TODO: How should this be done for other components (scalars, moisture)
245  });
246  }
247  } else if (turbChoice.pbl_type == PBLType::YSU) {
248  Error("YSU Model not implemented yet");
249 
250  /*
251  YSU PBL initially introduced by S.-Y. Hong, Y. Noh, and J. Dudhia, MWR, 2006 [HND06]
252 
253  Further Modifications from S.-Y. Hong, Q. J. R. Meteorol. Soc., 2010 [H10]
254 
255  Implementation follows WRF as of early 2024 with some simplifications
256  */
257 
258 #ifdef _OPENMP
259 #pragma omp parallel if (Gpu::notInLaunchRegion())
260 #endif
261  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
262 
263  // Pull out the box we're working on, make sure it covers full domain in z-direction
264  const Box &bx = mfi.growntilebox(1);
265  const Box &dbx = geom.Domain();
266  Box sbx(bx.smallEnd(), bx.bigEnd());
267  sbx.grow(2,-1);
268  AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
269 
270  // Get some data in arrays
271  const auto& cell_data = cons_in.const_array(mfi);
272  //const auto& K_turb = eddyViscosity.array(mfi);
273  const auto& uvel = xvel.const_array(mfi);
274  const auto& vvel = yvel.const_array(mfi);
275 
276  const auto& z0_arr = most->get_z0(level)->const_array();
277  const auto& ws10av_arr = most->get_mac_avg(level,4)->const_array(mfi);
278  const auto& t10av_arr = most->get_mac_avg(level,2)->const_array(mfi);
279  //const auto& u_star_arr = most->get_u_star(level)->const_array(mfi);
280  //const auto& t_star_arr = most->get_t_star(level)->const_array(mfi);
281  const auto& t_surf_arr = most->get_t_surf(level)->const_array(mfi);
282  //const auto& l_obuk_arr = most->get_olen(level)->const_array(mfi);
283  const auto& z_nd_arr = z_phys_nd->array(mfi);
284  const Real most_zref = most->get_zref();
285 
286  // Require that MOST zref is 10 m so we get the wind speed at 10 m from most
287  if (most_zref != 10.0) {
288  amrex::Abort("MOST Zref must be 10m for YSU PBL scheme");
289  }
290 
291  // create flattened boxes to store PBL height
292  const GeometryData gdata = geom.data();
293  const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
294  FArrayBox pbl_height(xybx,1);
295  const auto& pblh_arr = pbl_height.array();
296 
297  // -- Diagnose PBL height - starting out assuming non-moist
298  // loop is only over i,j in order to find height at each x,y
299  const Real f0 = turbChoice.pbl_ysu_coriolis_freq;
300  const bool over_land = turbChoice.pbl_ysu_over_land; // TODO: make this local and consistent
301  const Real land_Ribcr = turbChoice.pbl_ysu_land_Ribcr;
302  const Real unst_Ribcr = turbChoice.pbl_ysu_unst_Ribcr;
303  ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
304  {
305  // Reconstruct a surface bulk Richardson number from the surface layer model
306  // In WRF, this value is supplied to YSU by the MM5 surface layer model
307  const Real t_surf = t_surf_arr(i,j,0);
308  const Real t_layer = t10av_arr(i,j,0);
309  const Real ws_layer = ws10av_arr(i,j,0);
310  const Real Rib_layer = CONST_GRAV * most_zref / (ws_layer*ws_layer) * (t_layer - t_surf)/(t_layer);
311 
312  // For now, we only support stable boundary layers
313  if (Rib_layer < unst_Ribcr) {
314  Abort("For now, YSU PBL only supports stable conditions");
315  }
316 
317  // TODO: unstable BLs
318 
319  // PBL Height: Stable Conditions
320  Real Rib_cr;
321  if (over_land) {
322  Rib_cr = land_Ribcr;
323  } else { // over water
324  // Velocity at z=10 m comes from MOST -> currently the average using whatever averaging MOST uses.
325  // TODO: Revisit this calculation with local ws10?
326  const Real z0 = z0_arr(i,j,0);
327  const Real Rossby = ws_layer/(f0*z0);
328  Rib_cr = min(0.16*std::pow(1.0e-7*Rossby,-0.18),0.3); // Note: upper bound in WRF code, but not H10 paper
329  }
330 
331  bool above_critical = false;
332  int kpbl = 0;
333  Real Rib_up = Rib_layer, Rib_dn;
334  const Real base_theta = cell_data(i,j,0,RhoTheta_comp) / cell_data(i,j,0,Rho_comp);
335  while (!above_critical and bx.contains(i,j,kpbl+1)) {
336  kpbl += 1;
337  const Real zval = use_terrain ? Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr) : gdata.ProbLo(2) + (kpbl + 0.5)*gdata.CellSize(2);
338  const Real ws2_level = (uvel(i,j,kpbl)+uvel(i+1,j,kpbl))*(uvel(i,j,kpbl)+uvel(i+1,j,kpbl)) + (vvel(i,j,kpbl)+vvel(i,j+1,kpbl))*(uvel(i,j,kpbl)+uvel(i,j+1,kpbl));
339  const Real theta = cell_data(i,j,kpbl,RhoTheta_comp) / cell_data(i,j,kpbl,Rho_comp);
340  Rib_dn = Rib_up;
341  Rib_up = (theta-base_theta)/base_theta * CONST_GRAV * zval / ws2_level;
342  above_critical = Rib_up >= Rib_cr;
343  }
344 
345  Real interp_fact;
346  if (Rib_dn >= Rib_cr) {
347  interp_fact = 0.0;
348  } else if (Rib_up <= Rib_cr)
349  interp_fact = 1.0;
350  else {
351  interp_fact = (Rib_cr - Rib_dn) / (Rib_up - Rib_dn);
352  }
353 
354  const Real zval_up = use_terrain ? Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr) : gdata.ProbLo(2) + (kpbl + 0.5)*gdata.CellSize(2);
355  const Real zval_dn = use_terrain ? Compute_Zrel_AtCellCenter(i,j,kpbl-1,z_nd_arr) : gdata.ProbLo(2) + (kpbl-1 + 0.5)*gdata.CellSize(2);
356  pblh_arr(i,j,0) = zval_dn + interp_fact*(zval_up-zval_dn);
357  if (pblh_arr(i,j,0) < 0.5*(zval_up+zval_dn) ) {
358  kpbl = 0;
359  }
360  });
361 
362  // -- Compute nonlocal/countergradient mixing parameters
363 
364  // -- Compute entrainment parameters
365 
366  // -- Compute diffusion coefficients within PBL
367 
368  // -- Compute coefficients in free stream above PBL
369 
370  }
371 
372  }
373 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
RhoQKE_comp
#define RhoQKE_comp
Definition: IndexDefines.H:15
ComputeVerticalDerivativesPBL
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 amrex::Real &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)
Definition: PBLModels.H:14
Compute_Zrel_AtCellCenter
AMREX_GPU_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: TerrainMetrics.H:344
BCVars::yvel_bc
@ yvel_bc
Definition: IndexDefines.H:56
BCVars::cons_bc
@ cons_bc
Definition: IndexDefines.H:47
BCVars::xvel_bc
@ xvel_bc
Definition: IndexDefines.H:55
ERFBCType::ext_dir
@ ext_dir
Definition: IndexDefines.H:152
EddyDiff::PBL_lengthscale
@ PBL_lengthscale
Definition: IndexDefines.H:133
EddyDiff::QKE_v
@ QKE_v
Definition: IndexDefines.H:130
LsmVar_MM5::theta
@ theta
Definition: MM5.H:20
MicVar_Kess::rho
@ rho
Definition: Kessler.H:30
TurbChoice::pbl_mynn_A2
amrex::Real pbl_mynn_A2
Definition: TurbStruct.H:282
TurbChoice::pbl_mynn_C2
amrex::Real pbl_mynn_C2
Definition: TurbStruct.H:286
TurbChoice::pbl_mynn_A1
amrex::Real pbl_mynn_A1
Definition: TurbStruct.H:281
TurbChoice::pbl_mynn_C3
amrex::Real pbl_mynn_C3
Definition: TurbStruct.H:287
TurbChoice::pbl_mynn_B2
amrex::Real pbl_mynn_B2
Definition: TurbStruct.H:284
TurbChoice::pbl_ysu_over_land
bool pbl_ysu_over_land
Definition: TurbStruct.H:294
TurbChoice::pbl_mynn_C1
amrex::Real pbl_mynn_C1
Definition: TurbStruct.H:285
TurbChoice::pbl_mynn_C5
amrex::Real pbl_mynn_C5
Definition: TurbStruct.H:289
TurbChoice::pbl_ysu_land_Ribcr
amrex::Real pbl_ysu_land_Ribcr
Definition: TurbStruct.H:295
TurbChoice::pbl_ysu_coriolis_freq
amrex::Real pbl_ysu_coriolis_freq
Definition: TurbStruct.H:293
TurbChoice::pbl_mynn_B1
amrex::Real pbl_mynn_B1
Definition: TurbStruct.H:283
TurbChoice::pbl_ysu_unst_Ribcr
amrex::Real pbl_ysu_unst_Ribcr
Definition: TurbStruct.H:296

Referenced by ComputeTurbulentViscosity().

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