ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeDiffusivityMRF.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 "ERF_TileNoZ.H"
Include dependency graph for ERF_ComputeDiffusivityMRF.cpp:

Functions

void ComputeDiffusivityMRF (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, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
 

Function Documentation

◆ ComputeDiffusivityMRF()

void ComputeDiffusivityMRF ( 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  ,
int  level,
const BCRec *  bc_ptr,
bool  ,
const std::unique_ptr< MultiFab > &  z_phys_nd,
const MoistureComponentIndices &  moisture_indices 
)
26 {
27  /*
28  Implementation of the older MRF Scheme based on Hong and Pan (1996)
29  " Nonlocal Boundary Layer Vertical Diffusion in a Medium-Range Forecast
30  Model"
31  */
32 
33  // Domain extent in z-dir
34  int klo = geom.Domain().smallEnd(2);
35  int khi = geom.Domain().bigEnd(2);
36 
37 #ifdef _OPENMP
38 #pragma omp parallel if (Gpu::notInLaunchRegion())
39 #endif
40  for (MFIter mfi(eddyViscosity, TileNoZ()); mfi.isValid(); ++mfi) {
41 
42  // Box operated on must span fill domain in z-dir
43  const Box& gbx = mfi.growntilebox(IntVect(1,1,0));
44  AMREX_ALWAYS_ASSERT( gbx.smallEnd(2) == klo &&
45  gbx.bigEnd(2) == khi );
46 
47  //
48  // Step 1: Compute the height of the PBL without thermal excess
49  // From Hong et al. 2006, Eqns. 1 & 2:
50  //
51  // h = Rib_cf * theta_va * | U(h) |^2 / (g * (theta_v(h) - theta_s))
52  //
53  // where the surface virtual potential temperature is
54  //
55  // theta_s = theta_va + theta_T
56  //
57  // and here the thermal excess theta_T = 0.
58  //
59 
60  // create flattened boxes to store PBL height
61  const GeometryData gdata = geom.data();
62  const Box xybx = PerpendicularBox<ZDir>(gbx, IntVect{0, 0, 0});
63  FArrayBox pbl_height_predictor(xybx, 1, The_Async_Arena());
64  FArrayBox pbl_height_corrector(xybx, 1, The_Async_Arena());
65  IArrayBox pbl_index(xybx, 1, The_Async_Arena());
66  const auto& pblh_arr = pbl_height_predictor.array();
67  const auto& pblh_corr_arr = pbl_height_corrector.array();
68  const auto& pbli_arr = pbl_index.array();
69 
70  // Get some data in arrays
71  const auto& cell_data = cons_in.const_array(mfi);
72  const auto& uvel = xvel.const_array(mfi);
73  const auto& vvel = yvel.const_array(mfi);
74 
75  const Real Ribcr = turbChoice.pbl_mrf_Ribcr;
76  //const Real f0 = turbChoice.pbl_mrf_coriolis_freq;
77  const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
78  const auto& t_star_arr = SurfLayer->get_t_star(level)->const_array(mfi);
79  const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
80  const auto& t10av_arr = SurfLayer->get_mac_avg(level, 2)->const_array(mfi);
81  //const auto& t_surf_arr = SurfLayer->get_t_surf(level)->const_array(mfi);
82  const Array4<Real const> z_nd_arr = z_phys_nd->array(mfi);
83 
84  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
85  {
86  // note: if using a fixed surf_heat_flux (default), t_surf is not updated
87  //const Real t_surf = t_surf_arr(i, j, 0);
88  const Real t_layer = t10av_arr(i, j, 0);
89 
90  Real zval;
91  int kpbl = klo;
92  bool above_critical = false;
93  while (!above_critical && ((kpbl + 1) <= khi)) {
94  kpbl += 1;
95 
96  // height above ground level
97  zval = (use_terrain_fitted_coords)
98  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
99  : (kpbl + 0.5) * gdata.CellSize(2);
100 
101  const Real theta = cell_data(i, j, kpbl, RhoTheta_comp) /
102  cell_data(i, j, kpbl, Rho_comp);
103  const Real ws2 = 0.25 * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
104  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
105  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
106  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
107  const Real Rib = CONST_GRAV * zval * (theta - t_layer) / (ws2 * t_layer);
108  above_critical = (Rib >= Ribcr);
109  }
110 
111  // Empirical expression for PBLH is given by h = c u* / f
112  // Garratt (1994) and Tennekes (1982)
113  // Also, c.f. Zilitinkevitch et al 2012 referenced in Pedersen et al. 2014.
114  //const Real c_pblh = (l_obuk_arr(i, j, 0) > 0) ? 0.16 : 0.60;
115  //const Real pblh_emp = c_pblh * u_star_arr(i, j, 0) / f0;
116 
117  // Fallback to first cell
118  const Real pblh_emp = (use_terrain_fitted_coords)
119  ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr)
120  : gdata.ProbLo(2) + 0.5 * gdata.CellSize(2);
121 
122  // Initial PBL Height
123  // Avoiding detailed interpolation here
124  pblh_arr(i, j, 0) = (above_critical) ? zval : pblh_emp;
125  pbli_arr(i, j, 0) = (above_critical) ? kpbl : klo+1; // k < kpbl is considered the PBL
126  });
127 
128  //
129  // Step 2: Corrector PBL height for thermal excess
130  // where
131  //
132  // theta_T = b * surf_temp_flux / w_star
133  //
134  const Real const_b = turbChoice.pbl_mrf_const_b;
135  const Real sf = turbChoice.pbl_mrf_sf;
136  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
137  {
138  const Real t_layer = t10av_arr(i, j, 0);
139  const Real phiM = (l_obuk_arr(i, j, 0) > 0)
140  ? (1 + 5 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0))
141  : std::pow(
142  (1 - 8 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0)),
143  -1.0 / 3.0);
144  const Real wstar = u_star_arr(i, j, 0) / phiM;
145  const Real t_excess = -const_b * u_star_arr(i, j, 0) * t_star_arr(i, j, 0) / wstar;
146  const Real t_surf = t_layer + std::max(std::min(t_excess, 3.0), 0.0);
147 
148  int kpbl = klo;
149  Real zval0, zval, Rib0, Rib;
150  {
151  zval = (use_terrain_fitted_coords)
152  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
153  : gdata.ProbLo(2) + (kpbl + 0.5) * gdata.CellSize(2);
154  const Real theta = cell_data(i, j, kpbl, RhoTheta_comp) /
155  cell_data(i, j, kpbl, Rho_comp);
156  const Real ws2 = 0.25 * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
157  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
158  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
159  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
160  Rib = CONST_GRAV * zval * (theta - t_surf) / (ws2 * t_layer);
161  }
162 
163  bool above_critical = false;
164  while (!above_critical && ((kpbl + 1) <= khi)) {
165  zval0 = zval;
166  Rib0 = Rib;
167  kpbl += 1;
168 
169  zval = (use_terrain_fitted_coords)
170  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
171  : gdata.ProbLo(2) + (kpbl + 0.5) * gdata.CellSize(2);
172  const Real theta = cell_data(i, j, kpbl, RhoTheta_comp) /
173  cell_data(i, j, kpbl, Rho_comp);
174  const Real ws2 = 0.25 * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
175  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
176  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
177  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
178  Rib = CONST_GRAV * zval * (theta - t_surf) / (ws2 * t_layer);
179  above_critical = (Rib >= Ribcr);
180  }
181 
182  if (above_critical) {
183  // Interpolate to height at which Rib == Ribcr
184  Real pblh_interp = zval0 + (zval - zval0) / (Rib - Rib0) * (Ribcr - Rib0);
185  pblh_corr_arr(i, j, 0) = pblh_interp;
186  pbli_arr(i, j, 0) = kpbl; // k < kpbl is considered the PBL
187  } else {
188  // Empirical expression for PBLH is given by h = c u* / f
189  // Garratt (1994) and Tennekes (1982)
190  // Also, c.f. Zilitinkevitch et al 2012 referenced in Pedersen et al. 2014.
191  //const Real c_pblh = (l_obuk_arr(i, j, 0) > 0) ? 0.16 : 0.60;
192  //const Real pblh_emp = c_pblh * u_star_arr(i, j, 0) / f0;
193 
194  // Fallback to first cell
195  pblh_corr_arr(i, j, 0) = (use_terrain_fitted_coords)
196  ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr)
197  : gdata.ProbLo(2) + 0.5 * gdata.CellSize(2);
198  pbli_arr(i, j, 0) = klo + 1;
199  }
200 
201  });
202  /*
203  amrex::Print() << "PBL height computed for MRF scheme at level "
204  << pblh_arr(2, 2, 0) << " " << pblh_corr_arr(2, 2, 0)
205  << std::endl;
206  amrex::Print() << "PBL Temp:" << t_surf_arr(2, 2, 0) << " "
207  << t10av_arr(2, 2, 0) << std::endl;
208  */
209 
210  // -- Compute diffusion coefficients --
211 
212  const Array4<Real>& K_turb = eddyViscosity.array(mfi);
213 
214  // Dirichlet flags to switch derivative stencil
215  bool c_ext_dir_on_zlo = ((bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir));
216  bool c_ext_dir_on_zhi = ((bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir));
217  bool u_ext_dir_on_zlo = ((bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir));
218  bool u_ext_dir_on_zhi = ((bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir));
219  bool v_ext_dir_on_zlo = ((bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir));
220  bool v_ext_dir_on_zhi = ((bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir));
221 
222  const auto& dxInv = geom.InvCellSizeArray();
223  const Real dz_inv = geom.InvCellSize(2);
224  const int izmin = geom.Domain().smallEnd(2);
225  const int izmax = geom.Domain().bigEnd(2);
226 
227  ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
228  {
229  const Real zval = (use_terrain_fitted_coords)
230  ? Compute_Zrel_AtCellCenter(i, j, k, z_nd_arr)
231  : gdata.ProbLo(2) + (k + 0.5) * gdata.CellSize(2);
232  const Real rho = cell_data(i, j, k, Rho_comp);
233  const Real met_h_zeta = (use_terrain_fitted_coords)
234  ? Compute_h_zeta_AtCellCenter(i, j, k, dxInv, z_nd_arr) : 1.0;
235  const Real dz_terrain = met_h_zeta / dz_inv;
236  if (k < pbli_arr(i, j, 0)) {
237  const Real phiM = (l_obuk_arr(i, j, 0) > 0)
238  ? (1 + 5 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0))
239  : std::pow(
240  (1 - 8 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0)),
241  -1.0 / 3.0);
242  const Real phit = (l_obuk_arr(i, j, 0) > 0)
243  ? (1 + 5 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0))
244  : std::pow(
245  (1 - 16 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0)),
246  -1.0 / 2.0);
247  const Real Prt = phit / phiM + const_b * KAPPA * sf;
248  const Real wstar = u_star_arr(i, j, 0) / phiM;
249  K_turb(i, j, k, EddyDiff::Mom_v) = rho * wstar * KAPPA * zval
250  * (1 - zval / pblh_corr_arr(i, j, 0))
251  * (1 - zval / pblh_corr_arr(i, j, 0));
252  K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Prt;
253  } else {
254  const Real lambda = 150.0;
255  const Real lscale = (KAPPA * zval * lambda) / (KAPPA * zval + lambda);
256  Real dthetadz, dudz, dvdz;
257  ComputeVerticalDerivativesPBL(i, j, k, uvel, vvel, cell_data, izmin, izmax, 1.0 / dz_terrain,
258  c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo,
259  u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, dthetadz,
260  dudz, dvdz, moisture_indices);
261  const Real wind_shear = dudz * dudz + dvdz * dvdz + 1.0e-9;
262  const Real theta = cell_data(i, j, k, RhoTheta_comp) / cell_data(i, j, k, Rho_comp);
263  Real grad_Ri = CONST_GRAV / theta * dthetadz / wind_shear; // clear sky -- TODO: reduce stability in cloudy air
264  grad_Ri = std::max(grad_Ri, -100.0); // Hong et al. 2006, MWR, Appendix A
265  /*
266  const Real Pr = 1.5 + 3.08 * grad_Ri;
267  const Real fm =
268  (grad_Ri > 0)
269  ? (std::exp(-8.5 * grad_Ri) + (0.15 / (grad_Ri + 3.0)) * Pr)
270  : std::pow((1 - 12 * grad_Ri), -1.0 / 3.0);
271  const Real ft =
272  (grad_Ri > 0)
273  ? (std::exp(-8.5 * grad_Ri) + (0.15 / (grad_Ri + 3.0)))
274  : std::pow((1 - 16 * grad_Ri), -1.0 / 2.0);
275  */
276  // Using YSU model instead of MRF model
277  Real Pr = 1.0 + 2.1 * grad_Ri; // Hong et al. 2006, MWR, Eqn. A19
278  const Real fm = (grad_Ri > 0)
279  ? 1.0 / ((1.0 + 5.0 * grad_Ri) * (1.0 + 5.0 * grad_Ri))
280  : 1 - 8 * grad_Ri / (1 + 1.746 * std::sqrt(-grad_Ri)); // Hong et al. 2006, MWR, Eqn. A20b
281  const Real ft = (grad_Ri > 0)
282  ? 1.0 / ((1.0 + 5.0 * grad_Ri) * (1.0 + 5.0 * grad_Ri))
283  : 1 - 8 * grad_Ri / (1 + 1.286 * std::sqrt(-grad_Ri)); // Hong et al. 2006, MWR, Eqn. A20a
284  const Real rl2wsp = rho * lscale * lscale * std::sqrt(wind_shear);
285 
286  Pr = std::max(0.25, std::min(Pr, 4.0)); // Hong et al. 2006, MWR, Appendix A
287 
288  K_turb(i, j, k, EddyDiff::Mom_v) = rl2wsp * fm * Pr;
289  K_turb(i, j, k, EddyDiff::Theta_v) = rl2wsp * ft;
290  }
291 
292  // limit both diffusion coefficients
293 #if 0
294  // Hong et al. 2006, MWR, Appendix A
295  constexpr Real ckz = 0.001;
296  constexpr Real Kmax = 1000.0;
297  const Real rhoKmin = ckz * dz_terrain * rho;
298  const Real rhoKmax = rho * Kmax;
299 #endif
300  // Hong & Pan 1996, MWR
301  constexpr Real Kmin = 0.1;
302  constexpr Real Kmax = 300.0;
303  const Real rhoKmin = rho * Kmin;
304  const Real rhoKmax = rho * Kmax;
305 
306  K_turb(i, j, k, EddyDiff::Mom_v) = std::max(
307  std::min(K_turb(i, j, k, EddyDiff::Mom_v), rhoKmax), rhoKmin);
308  K_turb(i, j, k, EddyDiff::Theta_v) = std::max(
309  std::min(K_turb(i, j, k, EddyDiff::Theta_v), rhoKmax), rhoKmin);
310  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
311  K_turb(i, j, k, EddyDiff::Turb_lengthscale) = pblh_arr(i, j, 0);
312  });
313 
314  // FOEXTRAP top and bottom ghost cells
315  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
316  {
317  K_turb(i, j, klo-1, EddyDiff::Mom_v ) = K_turb(i, j, klo, EddyDiff::Mom_v );
318  K_turb(i, j, klo-1, EddyDiff::Theta_v) = K_turb(i, j, klo, EddyDiff::Theta_v);
319  K_turb(i, j, klo-1, EddyDiff::Q_v ) = K_turb(i, j, klo, EddyDiff::Q_v );
320  K_turb(i, j, khi+1, EddyDiff::Mom_v ) = K_turb(i, j, khi, EddyDiff::Mom_v );
321  K_turb(i, j, khi+1, EddyDiff::Theta_v) = K_turb(i, j, khi, EddyDiff::Theta_v);
322  K_turb(i, j, khi+1, EddyDiff::Q_v ) = K_turb(i, j, khi, EddyDiff::Q_v );
323  });
324  }// mfi
325 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
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, const MoistureComponentIndices &moisture_indices)
Definition: ERF_PBLModels.H:181
turbChoice
TurbChoice turbChoice
Definition: ERF_SetupVertDiff.H:5
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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: ERF_TerrainMetrics.H:53
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: ERF_TerrainMetrics.H:387
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:88
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:76
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:87
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:209
EddyDiff::Theta_v
@ Theta_v
Definition: ERF_IndexDefines.H:176
EddyDiff::Turb_lengthscale
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:180
EddyDiff::Q_v
@ Q_v
Definition: ERF_IndexDefines.H:179
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:175
LsmVar_MM5::theta
@ theta
Definition: ERF_MM5.H:20
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:22
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142
TurbChoice::pbl_mrf_const_b
amrex::Real pbl_mrf_const_b
Definition: ERF_TurbStruct.H:428
TurbChoice::pbl_mrf_Ribcr
amrex::Real pbl_mrf_Ribcr
Definition: ERF_TurbStruct.H:426
TurbChoice::pbl_mrf_sf
amrex::Real pbl_mrf_sf
Definition: ERF_TurbStruct.H:430

Referenced by ComputeTurbulentViscosity().

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