ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_ComputeDiffusivityMYNN25.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 dependency graph for ERF_ComputeDiffusivityMYNN25.cpp:

Macros

#define EXTRA_MYNN25_CHECKS   0
 

Functions

void ComputeDiffusivityMYNN25 (const MultiFab &xvel, const MultiFab &yvel, const 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 MoistureComponentIndices &moisture_indices)
 

Macro Definition Documentation

◆ EXTRA_MYNN25_CHECKS

#define EXTRA_MYNN25_CHECKS   0

Function Documentation

◆ ComputeDiffusivityMYNN25()

void ComputeDiffusivityMYNN25 ( const MultiFab &  xvel,
const MultiFab &  yvel,
const 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 MoistureComponentIndices moisture_indices 
)
27 {
28  auto mynn = turbChoice.pbl_mynn;
29  auto level2 = turbChoice.pbl_mynn_level2;
30 
31  Real Lt_alpha = (mynn.config == MYNNConfigType::CHEN2021) ? 0.1 : 0.23;
32 
33  // Dirichlet flags to switch derivative stencil
34  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
35  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) );
36  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
37  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) );
38  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
39  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) );
40 
41  // Epsilon
43 
44 #ifdef _OPENMP
45 #pragma omp parallel if (Gpu::notInLaunchRegion())
46 #endif
47  for ( MFIter mfi(eddyViscosity,false); mfi.isValid(); ++mfi) {
48 
49  const Box &bx = mfi.growntilebox(1);
50  const Array4<Real const>& cell_data = cons_in.array(mfi);
51  const Array4<Real >& K_turb = eddyViscosity.array(mfi);
52  const Array4<Real const>& uvel = xvel.array(mfi);
53  const Array4<Real const>& vvel = yvel.array(mfi);
54 
55  // Compute some quantities that are constant in each column
56  // Sbox is shrunk to only include the interior of the domain in the vertical direction to compute integrals
57  // Box includes one ghost cell in each direction
58  const Box &dbx = geom.Domain();
59  Box sbx(bx.smallEnd(), bx.bigEnd());
60  sbx.grow(2,-1);
61  AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
62 
63  const GeometryData gdata = geom.data();
64 
65  const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
66  FArrayBox qintegral(xybx,2);
67  qintegral.setVal<RunOn::Device>(0.0);
68  FArrayBox qturb(bx,1); FArrayBox qturb_old(bx,1);
69  const Array4<Real> qint = qintegral.array();
70  const Array4<Real> qvel = qturb.array();
71 
72  // vertical integrals to compute lengthscale
73  if (use_terrain_fitted_coords) {
74  const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
75  const auto invCellSize = geom.InvCellSizeArray();
76  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
77  {
78  // q^2 / 2 is the TKE
79  qvel(i,j,k) = std::sqrt(2.0 * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp));
80  AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "KE must have a positive value");
81 
82  Real fac = (sbx.contains(i,j,k)) ? 1.0 : 0.0;
83  const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr);
84  const Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd_arr);
85  Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz*fac);
86  Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*dz*fac);
87  });
88  } else {
89  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
90  {
91  // q^2 / 2 is the TKE
92  qvel(i,j,k) = std::sqrt(2.0 * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp));
93  AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "KE must have a positive value");
94 
95  // Not multiplying by dz: it's constant and would fall out when we divide qint0/qint1 anyway
96 
97  Real fac = (sbx.contains(i,j,k)) ? 1.0 : 0.0;
98  const Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
99  Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*fac);
100  Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*fac);
101  });
102  }
103 
104  Real dz_inv = geom.InvCellSize(2);
105  const auto& dxInv = geom.InvCellSizeArray();
106  int izmin = geom.Domain().smallEnd(2);
107  int izmax = geom.Domain().bigEnd(2);
108 
109  // Spatially varying MOST
110  Real d_kappa = KAPPA;
111  Real d_gravity = CONST_GRAV;
112 
113  const auto& t_mean_mf = SurfLayer->get_mac_avg(level,4); // theta_v
114  const auto& q_mean_mf = SurfLayer->get_mac_avg(level,3); // q_v
115  const auto& u_star_mf = SurfLayer->get_u_star(level);
116  const auto& t_star_mf = SurfLayer->get_t_star(level);
117  const auto& q_star_mf = SurfLayer->get_q_star(level);
118 
119  const auto& tm_arr = t_mean_mf->const_array(mfi);
120  const auto& qm_arr = q_mean_mf->const_array(mfi);
121  const auto& u_star_arr = u_star_mf->const_array(mfi);
122  const auto& t_star_arr = t_star_mf->const_array(mfi);
123  const auto& q_star_arr = (use_moisture) ? q_star_mf->const_array(mfi) : Array4<Real>{};
124 
125  const Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
126 
127  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
128  {
129  // Compute some partial derivatives that we will need (second order)
130  // U and V derivatives are interpolated to account for staggered grid
131  const Real met_h_zeta = use_terrain_fitted_coords ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : 1.0;
132 
133  Real dthetavdz, dudz, dvdz;
135  uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta,
136  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
137  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
138  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
139  dthetavdz, dudz, dvdz,
140  moisture_indices);
141 
142  // Spatially varying MOST
143  Real theta0 = tm_arr(i,j,0);
144  Real qv0 = qm_arr(i,j,0);
145  Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
146  Real surface_latent_heat{0};
147  if (use_moisture) {
148  // Compute buoyancy flux (Stull Eqn. 4.4.5d)
149  surface_latent_heat = -u_star_arr(i,j,0) * q_star_arr(i,j,0);
150  surface_heat_flux *= (1.0 + 0.61*qv0);
151  surface_heat_flux += 0.61 * theta0 * surface_latent_heat;
152  }
153 
154  Real l_obukhov;
155  if (std::abs(surface_heat_flux) > eps) {
156  l_obukhov = -( theta0 * u_star_arr(i,j,0)*u_star_arr(i,j,0)*u_star_arr(i,j,0) )
157  / ( d_kappa * d_gravity * surface_heat_flux );
158  } else {
159  l_obukhov = std::numeric_limits<Real>::max();
160  }
161 
162  // Surface-layer length scale (NN09, Eqn. 53)
163  AMREX_ASSERT(l_obukhov != 0);
164  int lk = amrex::max(k,0);
165  const Real zval = use_terrain_fitted_coords ? Compute_Zrel_AtCellCenter(i,j,lk,z_nd_arr)
166  : gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
167  const Real zeta = zval/l_obukhov;
168  Real l_S;
169  if (zeta >= 1.0) {
170  l_S = KAPPA*zval/3.7;
171  } else if (zeta >= 0) {
172  l_S = KAPPA*zval/(1+2.7*zeta);
173  } else {
174  l_S = KAPPA*zval*std::pow(1.0 - 100.0 * zeta, 0.2);
175  }
176 
177  // ABL-depth length scale (NN09, Eqn. 54)
178  Real l_T;
179  if (qint(i,j,0,1) > 0.0) {
180  l_T = Lt_alpha*qint(i,j,0,0)/qint(i,j,0,1);
181  } else {
182  l_T = std::numeric_limits<Real>::max();
183  }
184 
185  // Buoyancy length scale (NN09, Eqn. 55)
186  Real l_B;
187  if (dthetavdz > 0) {
188  Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetavdz);
189  if (zeta < 0) {
190  Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T; // velocity scale
191  qc = std::pow(qc,1.0/3.0);
192  l_B = (1.0 + 5.0*std::sqrt(qc/(N_brunt_vaisala * l_T))) * qvel(i,j,k)/N_brunt_vaisala;
193  } else {
194  l_B = qvel(i,j,k) / N_brunt_vaisala;
195  }
196  } else {
197  l_B = std::numeric_limits<Real>::max();
198  }
199 
200  // Master length scale
201  Real Lm;
202  if (mynn.config == MYNNConfigType::CHEN2021) {
203  Lm = std::pow(1.0/(l_S*l_S) + 1.0/(l_T*l_T) + 1.0/(l_B*l_B), -0.5);
204  } else {
205  // NN09, Eqn 52
206  Lm = 1.0 / (1.0/l_S + 1.0/l_T + 1.0/l_B);
207  }
208 
209  // Calculate nondimensional production terms
210  Real shearProd = dudz*dudz + dvdz*dvdz;
211  Real buoyProd = -(CONST_GRAV/theta0) * dthetavdz;
212  Real L2_over_q2 = Lm*Lm/(qvel(i,j,k)*qvel(i,j,k));
213  Real GM = L2_over_q2 * shearProd;
214  Real GH = L2_over_q2 * buoyProd;
215 
216  // Equilibrium (Level-2) q calculation follows NN09, Appendix A
217  Real Rf = level2.calc_Rf(GM, GH);
218  Real SM2 = level2.calc_SM(Rf);
219  Real qe2 = mynn.B1 * Lm*Lm * SM2 * (1.0-Rf) * shearProd;
220  Real qe = (qe2 < 0.0) ? 0.0 : std::sqrt(qe2);
221 
222  // Level 2 limiting (Helfand and Labraga 1988)
223  Real alphac = (qvel(i,j,k) >= qe) ? 1.0 : qvel(i,j,k) / (qe + eps);
224 #ifdef EXTRA_MYNN25_CHECKS
225  Real Ri = -GH/(GM+level2.eps);
226  if (alphac < 1 && (Ri > 1 || Ri < -1)) {
227  Warning("Level 2 limiting being applied with Ri out of expected range");
228  //AllPrint() << "alphac"<<IntVect(i,j,k)<<"= " << alphac
229  // << " Ri,SM2,SH2= " << Ri << " " << SM2 << " " << level2.calc_SH(Rf)
230  // << std::endl;
231  }
232 #endif
233 
234  // Level 2.5 stability functions
235  Real SM, SH, SQ;
236  mynn.calc_stability_funcs(SM,SH,SQ,GM,GH,alphac);
237 
238  // Clip SM, SH following WRF
239  SM = amrex::min(amrex::max(SM, mynn.SMmin), mynn.SMmax);
240  SH = amrex::min(amrex::max(SH, mynn.SHmin), mynn.SHmax);
241 #ifdef EXTRA_MYNN25_CHECKS
242  if (SM == mynn.SMmin) {
243  Warning("SM clipped at min val");
244  } else if (SM == mynn.SMmax) {
245  Warning("SM clipped at max val");
246  }
247  if (SH == mynn.SHmin) {
248  Warning("SH clipped at min val");
249  } else if (SH == mynn.SHmax) {
250  Warning("SH clipped at max val");
251  }
252 #endif
253 
254  // Finally, compute the eddy viscosity/diffusivities
255  const Real rho = cell_data(i,j,k,Rho_comp);
256  K_turb(i,j,k,EddyDiff::Mom_v) = rho * Lm * qvel(i,j,k) * SM;
257  K_turb(i,j,k,EddyDiff::Theta_v) = rho * Lm * qvel(i,j,k) * SH;
258  K_turb(i,j,k,EddyDiff::KE_v) = rho * Lm * qvel(i,j,k) * SQ;
259 
260  // TODO: implement partial-condensation scheme?
261  // Currently, implementation matches NN09 without rain (i.e.,
262  // the liquid water potential temperature is equal to the
263  // potential temperature.
264 
265  // NN09 gives the total water content flux; this assumes that
266  // all the species have the same eddy diffusivity
267  if (mynn.diffuse_moistvars) {
268  K_turb(i,j,k,EddyDiff::Q_v) = rho * Lm * qvel(i,j,k) * SH;
269  }
270 
271  K_turb(i,j,k,EddyDiff::Turb_lengthscale) = Lm;
272  });
273  }
274 }
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
TurbChoice turbChoice
Definition: ERF_DiffSetup.H:2
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
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, const MoistureComponentIndices &moisture_indices)
Definition: ERF_PBLModels.H:147
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: ERF_TerrainMetrics.H:47
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: ERF_TerrainMetrics.H:390
@ yvel_bc
Definition: ERF_IndexDefines.H:88
@ cons_bc
Definition: ERF_IndexDefines.H:76
@ xvel_bc
Definition: ERF_IndexDefines.H:87
@ ext_dir
Definition: ERF_IndexDefines.H:209
@ Theta_v
Definition: ERF_IndexDefines.H:176
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:180
@ Q_v
Definition: ERF_IndexDefines.H:179
@ Mom_v
Definition: ERF_IndexDefines.H:175
@ KE_v
Definition: ERF_IndexDefines.H:177
@ rho
Definition: ERF_Kessler.H:22
@ qc
Definition: ERF_SatAdj.H:36
@ xvel
Definition: ERF_IndexDefines.H:141
@ yvel
Definition: ERF_IndexDefines.H:142
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
MYNNLevel25 pbl_mynn
Definition: ERF_TurbStruct.H:372
MYNNLevel2 pbl_mynn_level2
Definition: ERF_TurbStruct.H:373

Referenced by ComputeTurbulentViscosity().

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