ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_ComputeDiffusivityMYNN25.cpp File Reference
#include "ERF_ABLMost.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:

Functions

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

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< ABLMost > &  most,
bool  use_moisture,
int  level,
const BCRec *  bc_ptr,
bool  ,
const std::unique_ptr< MultiFab > &  z_phys_nd,
const int  RhoQv_comp,
const int  RhoQc_comp,
const int  RhoQr_comp 
)
26 {
27  const bool use_terrain = (z_phys_nd != nullptr);
28  const bool use_most = (most != nullptr);
29 
30  auto mynn = turbChoice.pbl_mynn;
31  auto level2 = turbChoice.pbl_mynn_level2;
32 
33  Real Lt_alpha = (mynn.config == MYNNConfigType::CHEN2021) ? 0.1 : 0.23;
34 
35  // Dirichlet flags to switch derivative stencil
36  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
37  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) );
38  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
39  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) );
40  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
41  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) );
42 
43  // Epsilon
44  Real eps = std::numeric_limits<Real>::epsilon();
45 
46 #ifdef _OPENMP
47 #pragma omp parallel if (Gpu::notInLaunchRegion())
48 #endif
49  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
50 
51  const Box &bx = mfi.growntilebox(1);
52  const Array4<Real const>& cell_data = cons_in.array(mfi);
53  const Array4<Real >& K_turb = eddyViscosity.array(mfi);
54  const Array4<Real const>& uvel = xvel.array(mfi);
55  const Array4<Real const>& vvel = yvel.array(mfi);
56 
57  // Compute some quantities that are constant in each column
58  // Sbox is shrunk to only include the interior of the domain in the vertical direction to compute integrals
59  // Box includes one ghost cell in each direction
60  const Box &dbx = geom.Domain();
61  Box sbx(bx.smallEnd(), bx.bigEnd());
62  sbx.grow(2,-1);
63  AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
64 
65  const GeometryData gdata = geom.data();
66 
67  const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
68  FArrayBox qintegral(xybx,2);
69  qintegral.setVal<RunOn::Device>(0.0);
70  FArrayBox qturb(bx,1); FArrayBox qturb_old(bx,1);
71  const Array4<Real> qint = qintegral.array();
72  const Array4<Real> qvel = qturb.array();
73 
74  // vertical integrals to compute lengthscale
75  if (use_terrain) {
76  const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
77  const auto invCellSize = geom.InvCellSizeArray();
78  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
79  {
80  qvel(i,j,k) = std::sqrt(2.0 * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp));
81  AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "KE must have a positive value");
82 
83  Real fac = (sbx.contains(i,j,k)) ? 1.0 : 0.0;
84  const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr);
85  const Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd_arr);
86  Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz*fac);
87  Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*dz*fac);
88  });
89  } else {
90  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
91  {
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: its 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 = most->get_mac_avg(level,4); // theta_v
114  const auto& q_mean_mf = most->get_mac_avg(level,3); // q_v
115  const auto& u_star_mf = most->get_u_star(level);
116  const auto& t_star_mf = most->get_t_star(level);
117  const auto& q_star_mf = most->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 = use_terrain ? z_phys_nd->const_array(mfi) : Array4<Real>{};
126 
127  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
128  {
129  // NOTE: With MOST, the ghost cells are filled AFTER k_turb is computed
130  // so that the non-explicit pathway works. Therefore, at this
131  // point we do NOT have valid ghost cells from MOST. We need to
132  // pass the MOST flag to use one-sided diffs here.
133 
134  // Compute some partial derivatives that we will need (second order)
135  // U and V derivatives are interpolated to account for staggered grid
136  const Real met_h_zeta = use_terrain ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : 1.0;
137  Real dthetadz, dudz, dvdz;
139  uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta,
140  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
141  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
142  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
143  dthetadz, dudz, dvdz,
144  RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most);
145 
146  // Spatially varying MOST
147  Real theta0 = tm_arr(i,j,0);
148  Real qv0 = qm_arr(i,j,0);
149  Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
150  Real surface_latent_heat{0};
151  if (use_moisture) {
152  // Compute buoyancy flux (Stull Eqn. 4.4.5d)
153  surface_latent_heat = -u_star_arr(i,j,0) * q_star_arr(i,j,0);
154  surface_heat_flux *= (1.0 + 0.61*qv0);
155  surface_heat_flux += 0.61 * theta0 * surface_latent_heat;
156  }
157 
158  Real l_obukhov;
159  if (std::abs(surface_heat_flux) > eps) {
160  l_obukhov = -( theta0 * u_star_arr(i,j,0)*u_star_arr(i,j,0)*u_star_arr(i,j,0) )
161  / ( d_kappa * d_gravity * surface_heat_flux );
162  } else {
163  l_obukhov = std::numeric_limits<Real>::max();
164  }
165 
166  // Surface-layer length scale (NN09, Eqn. 53)
167  AMREX_ASSERT(l_obukhov != 0);
168  int lk = amrex::max(k,0);
169  const Real zval = use_terrain ? Compute_Zrel_AtCellCenter(i,j,lk,z_nd_arr)
170  : gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
171  const Real zeta = zval/l_obukhov;
172  Real l_S;
173  if (zeta >= 1.0) {
174  l_S = KAPPA*zval/3.7;
175  } else if (zeta >= 0) {
176  l_S = KAPPA*zval/(1+2.7*zeta);
177  } else {
178  l_S = KAPPA*zval*std::pow(1.0 - 100.0 * zeta, 0.2);
179  }
180 
181  // ABL-depth length scale (NN09, Eqn. 54)
182  Real l_T;
183  if (qint(i,j,0,1) > 0.0) {
184  l_T = Lt_alpha*qint(i,j,0,0)/qint(i,j,0,1);
185  } else {
186  l_T = std::numeric_limits<Real>::max();
187  }
188 
189  // Buoyancy length scale (NN09, Eqn. 55)
190  Real l_B;
191  if (dthetadz > 0) {
192  Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetadz);
193  if (zeta < 0) {
194  Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T; // velocity scale
195  qc = std::pow(qc,1.0/3.0);
196  l_B = (1.0 + 5.0*std::sqrt(qc/(N_brunt_vaisala * l_T))) * qvel(i,j,k)/N_brunt_vaisala;
197  } else {
198  l_B = qvel(i,j,k) / N_brunt_vaisala;
199  }
200  } else {
201  l_B = std::numeric_limits<Real>::max();
202  }
203 
204  // Master length scale
205  Real Lm;
206  if (mynn.config == MYNNConfigType::CHEN2021) {
207  Lm = std::pow(1.0/(l_S*l_S) + 1.0/(l_T*l_T) + 1.0/(l_B*l_B), -0.5);
208  } else {
209  // NN09, Eqn 52
210  Lm = 1.0 / (1.0/l_S + 1.0/l_T + 1.0/l_B);
211  }
212 
213  // Calculate nondimensional production terms
214  Real shearProd = dudz*dudz + dvdz*dvdz;
215  Real buoyProd = -(CONST_GRAV/theta0) * dthetadz;
216  Real L2_over_q2 = Lm*Lm/(qvel(i,j,k)*qvel(i,j,k));
217  Real GM = L2_over_q2 * shearProd;
218  Real GH = L2_over_q2 * buoyProd;
219 
220  // Equilibrium (Level-2) q calculation follows NN09, Appendix 2
221  Real Rf = level2.calc_Rf(GM, GH);
222  Real SM2 = level2.calc_SM(Rf);
223  Real qe2 = mynn.B1*Lm*Lm*SM2*(1.0-Rf)*shearProd;
224  Real qe = (qe2 < 0.0) ? 0.0 : std::sqrt(qe2);
225 
226  // Level 2 limiting (Helfand and Labraga 1988)
227  Real alphac = (qvel(i,j,k) > qe) ? 1.0 : qvel(i,j,k) / (qe + eps);
228 
229  // Level 2.5 stability functions
230  Real SM, SH, SQ;
231  mynn.calc_stability_funcs(SM,SH,SQ,GM,GH,alphac);
232 
233  // Clip SM, SH following WRF
234  SM = amrex::min(amrex::max(SM,mynn.SMmin), mynn.SMmax);
235  SH = amrex::min(amrex::max(SH,mynn.SHmin), mynn.SHmax);
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 * Lm * qvel(i,j,k) * SM;
240  K_turb(i,j,k,EddyDiff::Theta_v) = rho * Lm * qvel(i,j,k) * SH;
241  K_turb(i,j,k,EddyDiff::KE_v) = rho * Lm * qvel(i,j,k) * SQ;
242 
243  // TODO: implement partial-condensation scheme?
244  // Currently, implementation matches NN09 without rain (i.e.,
245  // the liquid water potential temperature is equal to the
246  // potential temperature.
247 
248  // NN09 gives the total water content flux; this assumes that
249  // all the species have the same eddy diffusivity
250  if (mynn.diffuse_moistvars) {
251  K_turb(i,j,k,EddyDiff::Q_v) = rho * Lm * qvel(i,j,k) * SH;
252  }
253 
254  K_turb(i,j,k,EddyDiff::Turb_lengthscale) = Lm;
255  });
256  }
257 }
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
#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 int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp, const bool use_most=true)
Definition: ERF_PBLModels.H:82
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:39
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:352
@ yvel_bc
Definition: ERF_IndexDefines.H:89
@ cons_bc
Definition: ERF_IndexDefines.H:76
@ xvel_bc
Definition: ERF_IndexDefines.H:88
@ ext_dir
Definition: ERF_IndexDefines.H:180
@ Theta_v
Definition: ERF_IndexDefines.H:157
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:161
@ Q_v
Definition: ERF_IndexDefines.H:160
@ Mom_v
Definition: ERF_IndexDefines.H:156
@ KE_v
Definition: ERF_IndexDefines.H:158
@ rho
Definition: ERF_Kessler.H:30
@ qc
Definition: ERF_SatAdj.H:36
@ xvel
Definition: ERF_IndexDefines.H:130
@ yvel
Definition: ERF_IndexDefines.H:131
MYNNLevel25 pbl_mynn
Definition: ERF_TurbStruct.H:203
MYNNLevel2 pbl_mynn_level2
Definition: ERF_TurbStruct.H:204

Referenced by ComputeTurbulentViscosity().

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