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  // Step 1: Compute the height of the PBL without thermal excess
48  // h = Rib_cf * theta_va * | U(h) |^2 / (g * (theta_va - theta_surf))
49  // create flattened boxes to store PBL height
50  const GeometryData gdata = geom.data();
51  const Box xybx = PerpendicularBox<ZDir>(gbx, IntVect{0, 0, 0});
52  FArrayBox pbl_height_predictor(xybx, 1, The_Async_Arena());
53  FArrayBox pbl_height_corrector(xybx, 1, The_Async_Arena());
54  IArrayBox pbl_index(xybx, 1, The_Async_Arena());
55  const auto& pblh_arr = pbl_height_predictor.array();
56  const auto& pblh_corr_arr = pbl_height_corrector.array();
57  const auto& pbli_arr = pbl_index.array();
58 
59  // Get some data in arrays
60  const auto& cell_data = cons_in.const_array(mfi);
61  const auto& uvel = xvel.const_array(mfi);
62  const auto& vvel = yvel.const_array(mfi);
63 
64  const Real Ribcr = turbChoice.pbl_mrf_Ribcr;
65  //const Real f0 = turbChoice.pbl_mrf_coriolis_freq;
66  const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
67  const auto& t_star_arr = SurfLayer->get_t_star(level)->const_array(mfi);
68  const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
69  const auto& t10av_arr = SurfLayer->get_mac_avg(level, 2)->const_array(mfi);
70  const auto& t_surf_arr = SurfLayer->get_t_surf(level)->const_array(mfi);
71  const Array4<Real const> z_nd_arr = z_phys_nd->array(mfi);
72 
73  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
74  {
75  const Real t_surf = t_surf_arr(i, j, 0);
76  const Real t_layer = t10av_arr(i, j, 0);
77 
78  Real zval;
79  int kpbl = klo;
80  bool above_critical = false;
81  while (!above_critical && ((kpbl + 1) <= khi)) {
82  kpbl += 1;
83 
84  // height above ground level
85  zval = (use_terrain_fitted_coords)
86  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
87  : (kpbl + 0.5) * gdata.CellSize(2);
88 
89  const Real theta = cell_data(i, j, kpbl, RhoTheta_comp) /
90  cell_data(i, j, kpbl, Rho_comp);
91  const Real ws2 = 0.25 * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
92  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
93  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
94  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
95  const Real Rib = CONST_GRAV * zval * (theta - t_surf) / (ws2 * t_layer);
96  above_critical = (Rib >= Ribcr);
97  }
98 
99  // Initial PBL Height
100  // Avoiding detailed interpolation here and just using map-nearest
101  // neighbor Empirical expression for PBLH is given by h = c u* / f
102  // Garratt (1994) and Tennekes (1982)
103  //const Real c_pblh = (l_obuk_arr(i, j, 0) > 0) ? 0.16 : 0.60;
104  //const Real pblh_emp = c_pblh * u_star_arr(i, j, 0) / f0;
105  const Real pblh_emp = gdata.ProbLo(2) + 0.5 * gdata.CellSize(2);
106  pblh_arr(i, j, 0) = (above_critical) ? zval : pblh_emp;
107  pbli_arr(i, j, 0) = (above_critical) ? kpbl : klo+1; // k < kpbl is considered the PBL
108  });
109 
110  // Corrector PBL height for thermal excess
111  const Real const_b = turbChoice.pbl_mrf_const_b;
112  const Real sf = turbChoice.pbl_mrf_sf;
113  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
114  {
115  const Real t_layer = t10av_arr(i, j, 0);
116  const Real phiM = (l_obuk_arr(i, j, 0) > 0)
117  ? (1 + 5 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0))
118  : std::pow(
119  (1 - 8 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0)),
120  -1.0 / 3.0);
121  const Real wstar = u_star_arr(i, j, 0) / phiM;
122  const Real t_excess = -const_b * u_star_arr(i, j, 0) * t_star_arr(i, j, 0) / wstar;
123  const Real t_surf = t_layer + std::max(std::min(t_excess, 3.0), 0.0);
124 
125  int kpbl = klo;
126  Real zval = 10;
127  bool above_critical = false;
128  while (!above_critical && ((kpbl + 1) <= khi)) {
129  zval = (use_terrain_fitted_coords)
130  ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr)
131  : gdata.ProbLo(2) + (kpbl + 0.5) * gdata.CellSize(2);
132  kpbl += 1;
133  const Real theta = cell_data(i, j, kpbl, RhoTheta_comp) /
134  cell_data(i, j, kpbl, Rho_comp);
135  const Real ws2 = 0.25 * ( (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) *
136  (uvel(i, j, kpbl) + uvel(i + 1, j, kpbl)) +
137  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) *
138  (vvel(i, j, kpbl) + vvel(i, j + 1, kpbl)) );
139  const Real Rib = CONST_GRAV * zval * (theta - t_surf) / (ws2 * t_layer);
140  above_critical = (Rib >= Ribcr);
141  }
142  //const Real c_pblh = (l_obuk_arr(i, j, 0) > 0) ? 0.16 : 0.60;
143  //const Real pblh_emp = c_pblh * u_star_arr(i, j, 0) / f0;
144  const Real pblh_emp = gdata.ProbLo(2) + 0.5 * gdata.CellSize(2);
145  pblh_corr_arr(i, j, 0) = (above_critical) ? zval : pblh_emp;
146  pbli_arr(i, j, 0) = (above_critical) ? kpbl : klo+1; // k < kpbl is considered the PBL
147  });
148  /*
149  amrex::Print() << "PBL height computed for MRF scheme at level "
150  << pblh_arr(2, 2, 0) << " " << pblh_corr_arr(2, 2, 0)
151  << std::endl;
152  amrex::Print() << "PBL Temp:" << t_surf_arr(2, 2, 0) << " "
153  << t10av_arr(2, 2, 0) << std::endl;
154  */
155 
156  // -- Compute diffusion coefficients --
157 
158  const Array4<Real>& K_turb = eddyViscosity.array(mfi);
159 
160  // Dirichlet flags to switch derivative stencil
161  bool c_ext_dir_on_zlo = ((bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir));
162  bool c_ext_dir_on_zhi = ((bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir));
163  bool u_ext_dir_on_zlo = ((bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir));
164  bool u_ext_dir_on_zhi = ((bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir));
165  bool v_ext_dir_on_zlo = ((bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir));
166  bool v_ext_dir_on_zhi = ((bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir));
167 
168  const auto& dxInv = geom.InvCellSizeArray();
169  const Real dz_inv = geom.InvCellSize(2);
170  const int izmin = geom.Domain().smallEnd(2);
171  const int izmax = geom.Domain().bigEnd(2);
172 
173  ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
174  {
175  const Real zval = (use_terrain_fitted_coords)
176  ? Compute_Zrel_AtCellCenter(i, j, k, z_nd_arr)
177  : gdata.ProbLo(2) + (k + 0.5) * gdata.CellSize(2);
178  const Real rho = cell_data(i, j, k, Rho_comp);
179  const Real met_h_zeta = (use_terrain_fitted_coords)
180  ? Compute_h_zeta_AtCellCenter(i, j, k, dxInv, z_nd_arr) : 1.0;
181  const Real dz_terrain = met_h_zeta / dz_inv;
182  if (k < pbli_arr(i, j, 0)) {
183  const Real phiM = (l_obuk_arr(i, j, 0) > 0)
184  ? (1 + 5 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0))
185  : std::pow(
186  (1 - 8 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0)),
187  -1.0 / 3.0);
188  const Real phit = (l_obuk_arr(i, j, 0) > 0)
189  ? (1 + 5 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0))
190  : std::pow(
191  (1 - 16 * sf * pblh_arr(i, j, 0) / l_obuk_arr(i, j, 0)),
192  -1.0 / 2.0);
193  const Real Prt = phit / phiM + const_b * KAPPA * sf;
194  const Real wstar = u_star_arr(i, j, 0) / phiM;
195  K_turb(i, j, k, EddyDiff::Mom_v) = rho * wstar * KAPPA * zval
196  * (1 - zval / pblh_corr_arr(i, j, 0))
197  * (1 - zval / pblh_corr_arr(i, j, 0));
198  K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Prt;
199  } else {
200  const Real lambda = 150.0;
201  const Real lscale = (KAPPA * zval * lambda) / (KAPPA * zval + lambda);
202  Real dthetadz, dudz, dvdz;
203  ComputeVerticalDerivativesPBL(i, j, k, uvel, vvel, cell_data, izmin, izmax, 1.0 / dz_terrain,
204  c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo,
205  u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, dthetadz,
206  dudz, dvdz, moisture_indices);
207  const Real wind_shear = dudz * dudz + dvdz * dvdz + 1.0e-9;
208  const Real theta = cell_data(i, j, k, RhoTheta_comp) / cell_data(i, j, k, Rho_comp);
209  const Real grad_Ri = std::max(CONST_GRAV / theta * dthetadz / wind_shear, -100.0);
210  /*
211  const Real Pr = 1.5 + 3.08 * grad_Ri;
212  const Real fm =
213  (grad_Ri > 0)
214  ? (std::exp(-8.5 * grad_Ri) + (0.15 / (grad_Ri + 3.0)) * Pr)
215  : std::pow((1 - 12 * grad_Ri), -1.0 / 3.0);
216  const Real ft =
217  (grad_Ri > 0)
218  ? (std::exp(-8.5 * grad_Ri) + (0.15 / (grad_Ri + 3.0)))
219  : std::pow((1 - 16 * grad_Ri), -1.0 / 2.0);
220  */
221  // Using YSU model instead of MRF model
222  const Real Pr = 1.0 + 2.1 * grad_Ri;
223  const Real fm = (grad_Ri > 0)
224  ? 1.0 / ((1.0 + 5.0 * grad_Ri) * (1.0 + 5.0 * grad_Ri))
225  : 1 - 8 * grad_Ri / (1 + 1.746 * std::sqrt(-grad_Ri));
226  const Real ft = (grad_Ri > 0)
227  ? 1.0 / ((1.0 + 5.0 * grad_Ri) * (1.0 + 5.0 * grad_Ri))
228  : 1 - 8 * grad_Ri / (1 + 1.286 * std::sqrt(-grad_Ri));
229  const Real rl2wsp = rho * lscale * lscale * std::sqrt(wind_shear);
230  K_turb(i, j, k, EddyDiff::Mom_v) = rl2wsp * fm * Pr;
231  K_turb(i, j, k, EddyDiff::Theta_v) = rl2wsp * ft;
232  }
233 
234  // limit both diffusion coefficients - from WRF, not documented in
235  // papers
236  constexpr Real ckz = 0.001;
237  constexpr Real Kmax = 1000.0;
238  const Real rhoKmin = ckz * dz_terrain * rho;
239  const Real rhoKmax = rho * Kmax;
240  K_turb(i, j, k, EddyDiff::Mom_v) = std::max(
241  std::min(K_turb(i, j, k, EddyDiff::Mom_v), rhoKmax), rhoKmin);
242  K_turb(i, j, k, EddyDiff::Theta_v) = std::max(
243  std::min(K_turb(i, j, k, EddyDiff::Theta_v), rhoKmax), rhoKmin);
244  K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v);
245  K_turb(i, j, k, EddyDiff::Turb_lengthscale) = pblh_arr(i, j, 0);
246  });
247 
248  // FOEXTRAP top and bottom ghost cells
249  ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int ) noexcept
250  {
251  K_turb(i, j, klo-1, EddyDiff::Mom_v ) = K_turb(i, j, klo, EddyDiff::Mom_v );
252  K_turb(i, j, klo-1, EddyDiff::Theta_v) = K_turb(i, j, klo, EddyDiff::Theta_v);
253  K_turb(i, j, klo-1, EddyDiff::Q_v ) = K_turb(i, j, klo, EddyDiff::Q_v );
254  K_turb(i, j, khi+1, EddyDiff::Mom_v ) = K_turb(i, j, khi, EddyDiff::Mom_v );
255  K_turb(i, j, khi+1, EddyDiff::Theta_v) = K_turb(i, j, khi, EddyDiff::Theta_v);
256  K_turb(i, j, khi+1, EddyDiff::Q_v ) = K_turb(i, j, khi, EddyDiff::Q_v );
257  });
258  }// mfi
259 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
turbChoice
TurbChoice turbChoice
Definition: ERF_DiffSetup.H:2
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
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:47
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:381
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:421
TurbChoice::pbl_mrf_Ribcr
amrex::Real pbl_mrf_Ribcr
Definition: ERF_TurbStruct.H:419
TurbChoice::pbl_mrf_sf
amrex::Real pbl_mrf_sf
Definition: ERF_TurbStruct.H:423

Referenced by ComputeTurbulentViscosity().

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