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

Referenced by ComputeTurbulentViscosity().

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