ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeDiffusivityYSU.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_ComputeDiffusivityYSU.cpp:

Functions

void ComputeDiffusivityYSU (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 std::unique_ptr< MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
 

Function Documentation

◆ ComputeDiffusivityYSU()

void ComputeDiffusivityYSU ( 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 std::unique_ptr< MultiFab > &  z_phys_cc,
const MoistureComponentIndices &  moisture_indices 
)

Compute eddy diffusivity using the YSU PBL scheme.

Parameters
[in]xvelx-velocity field
[in]yvely-velocity field
[in]cons_inConservative state
[out]eddyViscosityComputed eddy viscosity coefficients
[in]geomGeometry used for spacing and domain info
[in]turbChoiceTurbulence model options
[in]SurfLayerSurface layer model data
[in]use_terrain_fitted_coordsFlag to use terrain-fitted coordinates
[in]levelAMR level
[in]bc_ptrBoundary condition records
[in]z_phys_ndNodal physical heights
[in]z_phys_ccCell-centered physical heights
[in]moisture_indicesIndices for moisture components
44 {
45  /*
46  YSU PBL initially introduced by S.-Y. Hong, Y. Noh, and J. Dudhia, MWR, 2006 [HND06]
47 
48  Further Modifications from S.-Y. Hong, Q. J. R. Meteorol. Soc., 2010 [H10]
49 
50  Implementation follows WRF as of early 2024 with some simplifications
51  */
52 
53  const Real most_zref = SurfLayer->get_zref(level);
54 
55  // Require that MOST zref is 10 m so we get the wind speed at 10 m from most
56  bool invalid_zref = false;
57  if (use_terrain_fitted_coords) {
58  invalid_zref = most_zref != Real(10.0);
59  } else {
60  // zref gets reset to nearest cell center, so assert that zref is in the same cell as the 10m point
61  Real dz = geom.CellSize(2);
62  invalid_zref = int((most_zref - myhalf*dz)/dz) != int((Real(10.0) - myhalf*dz)/dz);
63  }
64  if (invalid_zref) {
65  Print() << "most_zref = " << most_zref << std::endl;
66  Abort("MOST Zref must be 10m for YSU PBL scheme");
67  }
68 
69 #ifdef _OPENMP
70 #pragma omp parallel if (Gpu::notInLaunchRegion())
71 #endif
72  // TileNoZ, not TilingIfNotGPU: every iterate must span the full column (asserted
73  // below), and the CPU default tile size would split any box of 16+ cells in z.
74  for ( MFIter mfi(eddyViscosity,TileNoZ()); mfi.isValid(); ++mfi) {
75 
76  // Pull out the box we're working on, make sure it covers full domain in z-direction
77  const Box &bx = mfi.growntilebox(1);
78  const Box &dbx = geom.Domain();
79  Box sbx(bx.smallEnd(), bx.bigEnd());
80  sbx.grow(2,-1);
81  AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
82 
83  // Get some data in arrays
84  const auto& cell_data = cons_in.const_array(mfi);
85  const auto& uvel = xvel.const_array(mfi);
86  const auto& vvel = yvel.const_array(mfi);
87 
88  const auto& z0_arr = SurfLayer->get_z0(level)->const_array(mfi);
89  const auto& ws10av_arr = SurfLayer->get_mac_avg(level,6)->const_array(mfi);
90  const auto& t10av_arr = SurfLayer->get_mac_avg(level,3)->const_array(mfi);
91  const auto& t_surf_arr = SurfLayer->get_t_surf(level)->const_array(mfi);
92  const auto& over_land_arr = (SurfLayer->get_lmask(level)) ? SurfLayer->get_lmask(level)->const_array(mfi) :
93  Array4<int> {};
94  const Array4<Real const> z_nd_arr = z_phys_nd->array(mfi);
95  const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
96 
97  // create flattened boxes to store PBL height
98  const GeometryData gdata = geom.data();
99  const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
100  FArrayBox pbl_height(xybx,1,The_Async_Arena());
101  IArrayBox pbl_index(xybx,1,The_Async_Arena());
102  const auto& pblh_arr = pbl_height.array();
103  const auto& pbli_arr = pbl_index.array();
104 
105  // -- Diagnose PBL height - starting out assuming non-moist --
106  // loop is only over i,j in order to find height at each x,y
107  const Real f0 = turbChoice.pbl_ysu_coriolis_freq;
108  const bool force_over_water = turbChoice.pbl_ysu_force_over_water;
109  const Real land_Ribcr = turbChoice.pbl_ysu_land_Ribcr;
110  const Real unst_Ribcr = turbChoice.pbl_ysu_unst_Ribcr;
111  ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
112  {
113  // Reconstruct a surface bulk Richardson number from the surface layer model
114  // In WRF, this value is supplied to YSU by the surface layer model
115  const Real t_surf = t_surf_arr(i,j,0);
116  const Real t_layer = t10av_arr(i,j,0);
117  const Real ws_layer = ws10av_arr(i,j,0);
118  const Real Rib_layer = CONST_GRAV * most_zref / (ws_layer*ws_layer) * (t_layer - t_surf)/(t_layer);
119 
120  // For now, we only support stable boundary layers
121  if (Rib_layer < unst_Ribcr) {
122  Abort("For now, YSU PBL only supports stable conditions");
123  }
124 
125  // TODO: unstable BLs
126 
127  // PBL Height: Stable Conditions
128  Real Rib_cr;
129  bool over_land = (over_land_arr) ? over_land_arr(i,j,0) : 1;
130  if (over_land && !force_over_water) {
131  Rib_cr = land_Ribcr;
132  } else { // over water
133  // Velocity at z=10 m comes from MOST -> currently the average using whatever averaging MOST uses.
134  // TODO: Revisit this calculation with local ws10?
135  const Real z0 = z0_arr(i,j,0);
136  const Real Rossby = ws_layer/(f0*z0);
137  Rib_cr = min(Real(0.16)*std::pow(Real(1.0e-7)*Rossby,-Real(0.18)),Real(0.3)); // Note: upper bound in WRF code, but not H10 paper
138  }
139 
140  bool above_critical = false;
141  int kpbl = 0;
142  Real Rib_up = Rib_layer, Rib_dn;
143  const Real base_theta = cell_data(i,j,0,RhoTheta_comp) / cell_data(i,j,0,Rho_comp);
144  while (!above_critical and bx.contains(i,j,kpbl+1)) {
145  kpbl += 1;
146  const Real zval = use_terrain_fitted_coords ?
147  Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr) : gdata.ProbLo(2) + (kpbl + myhalf)*gdata.CellSize(2);
148  const Real ws2_level = fourth*( (uvel(i,j,kpbl)+uvel(i+1,j ,kpbl))*(uvel(i,j,kpbl)+uvel(i+1,j ,kpbl))
149  + (vvel(i,j,kpbl)+vvel(i ,j+1,kpbl))*(vvel(i,j,kpbl)+vvel(i ,j+1,kpbl)) );
150  const Real theta = cell_data(i,j,kpbl,RhoTheta_comp) / cell_data(i,j,kpbl,Rho_comp);
151  Rib_dn = Rib_up;
152  Rib_up = (theta-base_theta)/base_theta * CONST_GRAV * zval / ws2_level;
153  above_critical = Rib_up >= Rib_cr;
154  }
155 
156  Real interp_fact;
157  if (Rib_dn >= Rib_cr) {
158  interp_fact = zero;
159  } else if (Rib_up <= Rib_cr)
160  interp_fact = one;
161  else {
162  interp_fact = (Rib_cr - Rib_dn) / (Rib_up - Rib_dn);
163  }
164 
165  const Real zval_up = use_terrain_fitted_coords ?
166  Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr) : gdata.ProbLo(2) + (kpbl + myhalf)*gdata.CellSize(2);
167  const Real zval_dn = use_terrain_fitted_coords ?
168  Compute_Zrel_AtCellCenter(i,j,kpbl-1,z_nd_arr) : gdata.ProbLo(2) + (kpbl-1 + myhalf)*gdata.CellSize(2);
169  pblh_arr(i,j,0) = zval_dn + interp_fact*(zval_up-zval_dn);
170 
171  const Real zval_0 = use_terrain_fitted_coords ?
172  Compute_Zrel_AtCellCenter(i,j,0,z_nd_arr) : gdata.ProbLo(2) + (myhalf)*gdata.CellSize(2);
173  const Real zval_1 = use_terrain_fitted_coords ?
174  Compute_Zrel_AtCellCenter(i,j,1,z_nd_arr) : gdata.ProbLo(2) + (Real(1.5))*gdata.CellSize(2);
175  if (pblh_arr(i,j,0) < myhalf*(zval_0+zval_1) ) {
176  kpbl = 0;
177  }
178  pbli_arr(i,j,0) = kpbl;
179  });
180 
181  // -- Compute nonlocal/countergradient mixing parameters --
182  // Not included for stable so nothing to do until unstable treatment is added
183 
184  // -- Compute entrainment parameters --
185  // 0 for stable so nothing to do?
186 
187  // -- Compute diffusion coefficients --
188 
189  const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
190  const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
191  const Array4<Real > &K_turb = eddyViscosity.array(mfi);
192 
193  // Dirichlet flags to switch derivative stencil
194  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
195  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].hi(2) == ERFBCType::ext_dir) );
196  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
197  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) );
198  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
199  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) );
200 
201  const auto& dxInv = geom.InvCellSizeArray();
202  const Real dz_inv = geom.InvCellSize(2);
203  const int izmin = geom.Domain().smallEnd(2);
204  const int izmax = geom.Domain().bigEnd(2);
205 
206  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
207  {
208  const Real zval = use_terrain_fitted_coords ?
209  Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr) : gdata.ProbLo(2) + (k + myhalf)*gdata.CellSize(2);
210  const Real rho = cell_data(i,j,k,Rho_comp);
211  const Real met_h_zeta = use_terrain_fitted_coords ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : one;
212  const Real dz_terrain = met_h_zeta/dz_inv;
213  if (k < pbli_arr(i,j,0)) {
214  // -- Compute diffusion coefficients within PBL
215  constexpr Real zfacmin = Real(1e-8); // value from WRF
216  constexpr Real phifac = Real(8.0); // value from H10 and WRF
217  constexpr Real wstar3 = Real(0.); // only nonzero for unstable
218  constexpr Real pfac = Real(2.); // profile exponent
219  const Real zfac = std::min(std::max(amrex::Real(1) - zval / pblh_arr(i,j,0), zfacmin ), amrex::Real(1));
220  // Not including YSU top down PBL term (not in H10, added to WRF later)
221  const Real ust3 = u_star_arr(i,j,0) * u_star_arr(i,j,0) * u_star_arr(i,j,0);
222  Real wscalek = ust3 + phifac * KAPPA * wstar3 * (amrex::Real(1) - zfac);
223  wscalek = std::pow(wscalek, amrex::Real(1.0/3.0));
224  // stable only
225  const Real phi_term = amrex::Real(1) + amrex::Real(5) * zval / l_obuk_arr(i,j,0); // phi_term appears in WRF but not papers
226  wscalek = std::max(u_star_arr(i,j,0) / phi_term, Real(0.001)); // Real(0.001) limit appears in WRF but not papers
227  K_turb(i,j,k,EddyDiff::Mom_v) = rho * wscalek * KAPPA * zval * std::pow(zfac, pfac);
228  K_turb(i,j,k,EddyDiff::Theta_v) = K_turb(i,j,k,EddyDiff::Mom_v);
229  } else {
230  // -- Compute coefficients in free stream above PBL
231  constexpr Real lam0 = Real(30.0);
232  constexpr Real min_richardson = -Real(100.0);
233  constexpr Real prandtl_max = Real(4.0);
234  Real dthetadz, dudz, dvdz;
235  ComputeVerticalDerivativesPBL(i, j, k,
236  uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
237  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
238  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
239  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
240  dthetadz, dudz, dvdz, moisture_indices);
241  const Real shear_squared = dudz*dudz + dvdz*dvdz + Real(1.0e-9); // Real(1.0e-9) from WRF to avoid divide by zero
242  const Real theta = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);
243  Real richardson = CONST_GRAV / theta * dthetadz / shear_squared;
244  const Real lambdadz = std::min(std::max(Real(0.1)*dz_terrain , lam0), Real(300.0)); // in WRF, H10 paper just says use lam0
245  const Real lengthscale = lambdadz * KAPPA * zval / (lambdadz + KAPPA * zval);
246  const Real turbfact = lengthscale * lengthscale * std::sqrt(shear_squared);
247 
248  if (richardson < 0) {
249  richardson = max(richardson, min_richardson);
250  Real sqrt_richardson = std::sqrt(-richardson);
251  K_turb(i,j,k,EddyDiff::Mom_v) = rho * turbfact * (one - Real(8.0) * richardson / (one + Real(1.746) * sqrt_richardson));
252  K_turb(i,j,k,EddyDiff::Theta_v) = rho * turbfact * (one - Real(8.0) * richardson / (one + Real(1.286) * sqrt_richardson));
253  } else {
254  const Real oneplus5ri = one + Real(5.0) * richardson;
255  K_turb(i,j,k,EddyDiff::Theta_v) = rho * turbfact / (oneplus5ri * oneplus5ri);
256  const Real prandtl = std::min(one+Real(2.1)*richardson, prandtl_max); // limit from WRF
257  K_turb(i,j,k,EddyDiff::Mom_v) = K_turb(i,j,k,EddyDiff::Theta_v) * prandtl;
258  }
259  }
260 
261  // limit both diffusion coefficients - from WRF, not documented in papers
262  constexpr Real ckz = Real(0.001);
263  constexpr Real Kmax = Real(1000.0);
264  const Real rhoKmin = ckz * dz_terrain * rho;
265  const Real rhoKmax = rho * Kmax;
266  K_turb(i,j,k,EddyDiff::Mom_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Mom_v) ,rhoKmax), rhoKmin);
267  K_turb(i,j,k,EddyDiff::Theta_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Theta_v) ,rhoKmax), rhoKmin);
268  K_turb(i,j,k,EddyDiff::Turb_lengthscale) = pblh_arr(i,j,0);
269  });
270 
271  // HACK set bottom ghost cell to 1st cell
272  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
273  {
274  if (k==-1) {
275  K_turb(i,j,k,EddyDiff::Mom_v) = K_turb(i,j,0,EddyDiff::Mom_v);
276  K_turb(i,j,k,EddyDiff::Theta_v) = K_turb(i,j,0,EddyDiff::Theta_v);
277  }
278  });
279  } // mfi
280 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:55
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:56
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:35
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
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 PBLDerivativeDzInv &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:281
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Compute_Zrel_AtCellCenter
AMREX_GPU_HOST_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:751
Compute_h_zeta_AtCellCenter
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:190
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
AL01::richardson
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real richardson(amrex::Real length, amrex::Real N2, amrex::Real tke, amrex::Real Cmu0_pow3)
Definition: ERF_RANSClosure.H:110
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:106
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:89
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:105
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:297
EddyDiff::Theta_v
@ Theta_v
Definition: ERF_IndexDefines.H:250
EddyDiff::Turb_lengthscale
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:254
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:249
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:19
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
module_model_constants::prandtl
real(c_double), parameter prandtl
Definition: ERF_module_model_constants.F90:88
PBLDerivativeDzInv_T
Functor for inverse vertical spacings for terrain-following grids using cell-center heights.
Definition: ERF_PBLModels.H:461
TurbChoice::pbl_ysu_force_over_water
bool pbl_ysu_force_over_water
Whether YSU is forced to use over-water behavior.
Definition: ERF_TurbStruct.H:859
TurbChoice::pbl_ysu_land_Ribcr
amrex::Real pbl_ysu_land_Ribcr
Critical bulk Richardson number over land for stable YSU conditions.
Definition: ERF_TurbStruct.H:862
TurbChoice::pbl_ysu_unst_Ribcr
amrex::Real pbl_ysu_unst_Ribcr
Critical bulk Richardson number for unstable YSU conditions.
Definition: ERF_TurbStruct.H:864
TurbChoice::pbl_ysu_coriolis_freq
amrex::Real pbl_ysu_coriolis_freq
Coriolis frequency used by YSU-family PBL schemes.
Definition: ERF_TurbStruct.H:853

Referenced by ComputeTurbulentViscosity().

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