ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Static Public Attributes | List of all members
MYNNPBLH Struct Reference

Diagnostic utility for the planetary boundary layer height. More...

#include <ERF_PBLHeight.H>

Collaboration diagram for MYNNPBLH:
Collaboration graph
[legend]

Public Member Functions

AMREX_GPU_HOST AMREX_FORCE_INLINE void compute_pblh (const amrex::Geometry &geom, const amrex::MultiFab *z_phys_cc, amrex::MultiFab *pblh, const amrex::MultiFab &cons, const amrex::iMultiFab *lmask, const MoistureComponentIndices &moisture_indices) const
 

Static Public Attributes

static constexpr amrex::Real thetamin_height = amrex::Real(200.0)
 Height below which minimum theta-v is determined [m]. More...
 
static constexpr amrex::Real theta_incr_land = amrex::Real(1.25)
 Theta increase determining the capping inversion height over land [K]. More...
 
static constexpr amrex::Real theta_incr_water = one
 Theta increase determining the capping inversion height over water [K]. More...
 
static constexpr amrex::Real sbl_lim = amrex::Real(200.0)
 Upper limit of the stable boundary layer height [m]. More...
 
static constexpr amrex::Real sbl_damp = amrex::Real(400.0)
 Transition length for blending [m]. More...
 

Detailed Description

Diagnostic utility for the planetary boundary layer height.

Member Function Documentation

◆ compute_pblh()

AMREX_GPU_HOST AMREX_FORCE_INLINE void MYNNPBLH::compute_pblh ( const amrex::Geometry &  geom,
const amrex::MultiFab *  z_phys_cc,
amrex::MultiFab *  pblh,
const amrex::MultiFab &  cons,
const amrex::iMultiFab *  lmask,
const MoistureComponentIndices &  moisture_indices 
) const
inline
32  {
33 #if 0
34  // NOTE: Cannot use ReduceToPlane because it clips the box to the
35  // validbox only, i.e., lateral ghost cells aren't updated
36  int dir = 2; // z
37  auto const& cons_arrs = cons.const_arrays();
38  auto thetav_min = amrex::ReduceToPlane<amrex::ReduceOpMin,amrex::Real>(dir, bxlow, cons,
39  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> amrex::Real
40  {
41  return GetThetav(i,j,k,cons_arrs[box_no],moisture_indices);
42  });
43 #endif
44 
45  // Create 2D multifabs like pblh
46  auto const& ba = pblh->boxArray();
47  auto const& dm = pblh->DistributionMap();
48  auto const& ng = pblh->nGrowVect();
49 
50  amrex::MultiFab min_thetav(ba,dm,1,ng);
51  min_thetav.setVal(bogus_large_value);
52 
53  amrex::MultiFab pblh_tke(ba,dm,1,ng);
54  pblh_tke.setVal(zero);
55 
56  pblh->setVal(zero);
57 
58  // The scan below ends at the top cell of each box and reads the cell above it, so
59  // the fields it reads must carry a ghost cell in z. On grids that all start at the
60  // ground that is the state's own halo; on grids that do not, it is the halo
61  // SurfaceLayer::compute_pblh copies onto the columns it hands over.
62  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(cons.nGrowVect()[2] >= 1,
63  "MYNNPBLH::compute_pblh: the state must carry a ghost cell in z; the scan reads "
64  "one cell above the top of every box");
65  AMREX_ALWAYS_ASSERT_WITH_MESSAGE((z_phys_cc == nullptr) || (z_phys_cc->nGrowVect()[2] >= 1),
66  "MYNNPBLH::compute_pblh: z_phys_cc must carry a ghost cell in z; the scan reads "
67  "one cell above the top of every box");
68 
69  // Now, loop over columns...
70  for (amrex::MFIter mfi(cons,TileNoZ()); mfi.isValid(); ++mfi)
71  {
72  const amrex::Box& domain = geom.Domain();
73 
74  // NOTE: this MFIter is over cons, which carries more ghost cells than the
75  // 2D fields we write here (pblh, pblh_tke, min_thetav all have the
76  // ghost cells of pblh). Grow by *their* ghost vector, not by the
77  // default, or we write past the end of their allocations.
78  //
79  // The scan covers the cells this box holds and no more. A box that stops
80  // below the top of the domain (a refined level of partial height) has no
81  // data above its top ghost cell, so the scan must stop with it; the
82  // k+1 reads below then end in that ghost cell. The caller hands over
83  // boxes that start at the ground (SurfaceLayer::compute_pblh builds them
84  // when the grids are stacked in z).
85  const amrex::Box& vbx = mfi.validbox();
86  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vbx.smallEnd(2) == domain.smallEnd(2),
87  "MYNNPBLH::compute_pblh: every box must start at the bottom of the domain");
88  amrex::Box gtbx = mfi.growntilebox(ng);
89  gtbx.setSmall(2,vbx.smallEnd(2)); // don't loop over ghost cells
90  gtbx.setBig(2,vbx.bigEnd(2)); // in z
91 
92  // Columns are parallelized over; the k loops below are serialized within
93  // each thread. Everything in this routine is a per-column scan -- the
94  // minimum theta_v, and the lowest k at which each PBLH criterion is met --
95  // so parallelizing over k would both race on the 2D output and lose the
96  // "lowest k wins" semantics.
97  amrex::Box gtbx2d = gtbx; gtbx2d.setRange(2,0);
98  const int klo = gtbx.smallEnd(2);
99  const int khi = gtbx.bigEnd(2);
100 
101  auto min_thv_arr = min_thetav.array(mfi);
102  auto pblh_arr = pblh->array(mfi);
103  auto pblh_tke_arr = pblh_tke.array(mfi);
104 
105  const auto cons_arr = cons.const_array(mfi);
106  const auto lmask_arr = (lmask) ? lmask->const_array(mfi) : amrex::Array4<int> {};
107 
108  // -----------------------------------------------------
109  // WITH terrain/grid stretching
110  // -----------------------------------------------------
111  if (z_phys_cc)
112  {
113  const auto zphys_arr = z_phys_cc->const_array(mfi);
114 
115  // Need to sort out ghost cell differences (z_phys_cc has ng=1)
116  int imin = lbound(zphys_arr).x;
117  int jmin = lbound(zphys_arr).y;
118  int imax = ubound(zphys_arr).x;
119  int jmax = ubound(zphys_arr).y;
120 
121  // Heights here are above the ground. The ground is the bottom face of the
122  // column, where z_phys_cc in the ghost cell below mirrors the bottom cell.
123 
124  // Find minimum thetav in the surface layer (this updates
125  // ghost cells, too)
126  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
127  {
128  int ii = amrex::max(amrex::min(i,imax),imin);
129  int jj = amrex::max(amrex::min(j,jmax),jmin);
130  const amrex::Real z_sfc = myhalf * (zphys_arr(ii,jj,klo) + zphys_arr(ii,jj,klo-1));
131 
132  amrex::Real min_thv = min_thv_arr(i,j,0);
133  for (int k(klo); k <= khi; ++k) {
134  if (zphys_arr(ii,jj,k) - z_sfc < thetamin_height) {
135  amrex::Real thv = GetThetav(i, j, k, cons_arr, moisture_indices);
136  min_thv = amrex::min(min_thv, thv);
137  }
138  }
139  min_thv_arr(i,j,0) = min_thv;
140  });
141 
142  // This depends on TileNoZ and k increasing monotonically
143  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
144  {
145  int ii = amrex::max(amrex::min(i,imax),imin);
146  int jj = amrex::max(amrex::min(j,jmax),jmin);
147  const amrex::Real z_sfc = myhalf * (zphys_arr(ii,jj,klo) + zphys_arr(ii,jj,klo-1));
148 
149  const amrex::Real min_thv = min_thv_arr(i,j,0);
150  const int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
151 
152  amrex::Real zi = pblh_arr(i,j,0);
153  amrex::Real zi_tke = pblh_tke_arr(i,j,0);
154 
155  for (int k(klo); k <= khi; ++k)
156  {
157  if (zi == 0)
158  {
159  //
160  // Find PBL height based on thetav increase (best for CBLs)
161  //
162  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
163  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
164 
165  if (is_land && (thv1 >= min_thv + theta_incr_land)
166  && (thv < min_thv + theta_incr_land))
167  {
168  // Interpolate to get lowest height where theta = min_theta + theta_incr
169  zi = zphys_arr(ii,jj,k) - z_sfc
170  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
171  * (min_thv + theta_incr_land - thv);
172  }
173  else if (!is_land && (thv1 >= min_thv + theta_incr_water)
174  && (thv < min_thv + theta_incr_water))
175  {
176  // Interpolate to get lowest height where theta = min_theta + theta_incr
177  zi = zphys_arr(ii,jj,k) - z_sfc
178  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
179  * (min_thv + theta_incr_water - thv);
180  }
181  }
182  if (zi_tke == 0)
183  {
184  //
185  // Find PBL height based on TKE (for SBLs only)
186  //
187  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
188  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
189  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
190  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
191  amrex::Real TKEeps = amrex::Real(0.05) * maxtke;
192  TKEeps = amrex::max(TKEeps, amrex::Real(0.02)); // min val from WRF
193 
194  if ((tke1 <= TKEeps) && (tke > TKEeps))
195  {
196  // Interpolate to get lowest height where TKE -> 0
197  zi_tke = zphys_arr(ii,jj,k) - z_sfc
198  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(tke1-tke)
199  * (TKEeps - tke);
200  }
201  }
202  // Both criteria are satisfied; the remaining levels are no-ops
203  if ((zi != 0) && (zi_tke != 0)) { break; }
204  }
205 
206  pblh_arr(i,j,0) = zi;
207  pblh_tke_arr(i,j,0) = zi_tke;
208  });
209  }
210  else
211  // -----------------------------------------------------
212  // NO terrain
213  // -----------------------------------------------------
214  {
215  const amrex::Real dz_no_terrain = geom.CellSize(2);
216 
217  // Find minimum thetav in the surface layer (this updates
218  // ghost cells, too)
219  // - box size is known a priori
220  int kmax = static_cast<int>(thetamin_height / dz_no_terrain);
221  AMREX_ASSERT(kmax > 0);
222  const int khi_low = amrex::min(kmax, khi);
223 
224  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
225  {
226  amrex::Real min_thv = min_thv_arr(i,j,0);
227  for (int k(klo); k <= khi_low; ++k) {
228  amrex::Real thv = GetThetav(i, j, k, cons_arr, moisture_indices);
229  min_thv = amrex::min(min_thv, thv);
230  }
231  min_thv_arr(i,j,0) = min_thv;
232  });
233 
234  // This depends on TileNoZ and k increasing monotonically
235  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
236  {
237  const amrex::Real min_thv = min_thv_arr(i,j,0);
238  const int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
239 
240  amrex::Real zi = pblh_arr(i,j,0);
241  amrex::Real zi_tke = pblh_tke_arr(i,j,0);
242 
243  for (int k(klo); k <= khi; ++k)
244  {
245  if (zi == 0)
246  {
247  //
248  // Find PBL height based on thetav increase (best for CBLs)
249  //
250  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
251  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
252 
253  if (is_land && (thv1 >= min_thv + theta_incr_land)
254  && (thv < min_thv + theta_incr_land))
255  {
256  // Interpolate to get lowest height where theta = min_theta + theta_incr
257  zi = (k+myhalf)*dz_no_terrain
258  + dz_no_terrain/(thv1-thv)
259  * (min_thv + theta_incr_land - thv);
260  }
261  else if (!is_land && (thv1 >= min_thv + theta_incr_water)
262  && (thv < min_thv + theta_incr_water))
263  {
264  // Interpolate to get lowest height where theta = min_theta + theta_incr
265  zi = (k+myhalf)*dz_no_terrain
266  + dz_no_terrain/(thv1-thv)
267  * (min_thv + theta_incr_water - thv);
268  }
269  }
270  if (zi_tke == 0)
271  {
272  //
273  // Find PBL height based on TKE (for SBLs only)
274  //
275  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
276  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
277  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
278  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
279  amrex::Real TKEeps = amrex::Real(0.05) * maxtke;
280  TKEeps = amrex::max(TKEeps, amrex::Real(0.02)); // min val from WRF
281 
282  if ((tke1 <= TKEeps) && (tke > TKEeps))
283  {
284  // Interpolate to get lowest height where TKE -> 0
285  zi_tke = (k+myhalf)*dz_no_terrain
286  + dz_no_terrain/(tke1-tke) * (TKEeps - tke);
287  }
288  }
289  // Both criteria are satisfied; the remaining levels are no-ops
290  if ((zi != 0) && (zi_tke != 0)) { break; }
291  }
292 
293  pblh_arr(i,j,0) = zi;
294  pblh_tke_arr(i,j,0) = zi_tke;
295  });
296  }
297  }// MFIter
298 
299  //
300  // Calculate hybrid PBL height
301  //
302  for (amrex::MFIter mfi(*pblh); mfi.isValid(); ++mfi)
303  {
304  const auto cons_arr = cons.const_array(mfi);
305  auto pblh_tke_arr = pblh_tke.array(mfi);
306  auto pblh_arr = pblh->array(mfi);
307 
308  amrex::Box gtbx = mfi.growntilebox();
309  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
310  {
311  //
312  // Clip PBLH_TKE to more realistic values
313  //
314  // Note from WRF MYNN-EDMF: TKE-based PBLH can be very large in cells
315  // with convective precipitation (> 8km!), so an artificial limit is
316  // imposed to not let PBLH_TKE exceed the theta_v-based PBL height
317  // +/- 350 m. This has no impact on 98-99% of the domain, but is the
318  // simplest patch that adequately addresses these extremely large
319  // PBLHs.
320  amrex::Real zi = pblh_arr(i,j,0);
321  pblh_tke_arr(i,j,0) = amrex::max(
322  amrex::min(pblh_tke_arr(i,j,0), zi+amrex::Real(350.)),
323  amrex::max(zi-amrex::Real(350.), amrex::Real(10.)));
324 
325  //
326  // Finally, blend between the two PBLH estimates
327  //
328  amrex::Real maxqke = two * cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
329  if (maxqke > amrex::Real(0.05)) {
330  amrex::Real wt = myhalf*std::tanh((zi - sbl_lim)/sbl_damp) + myhalf;
331  pblh_arr(i,j,0) = (one-wt)*pblh_tke_arr(i,j,0) + wt*pblh_arr(i,j,0);
332  }
333  });
334  }//MFIter
335  }
bogus_large_value
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:19
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:53
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:55
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(axis >=0 &&axis< AMREX_SPACEDIM, "prob.axis must be 0, 1 or 2")
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:25
klo
const int klo
Definition: ERF_InitCustomPert_KE.H:9
GetThetav
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetav(const int &i, const int &j, const int &k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:93
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);})
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:51
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:254
WDM6SedNodeScratch::zi
@ zi
Definition: ERF_AdvanceWDM6.cpp:277
MYNNPBLH::sbl_lim
static constexpr amrex::Real sbl_lim
Upper limit of the stable boundary layer height [m].
Definition: ERF_PBLHeight.H:352
MYNNPBLH::sbl_damp
static constexpr amrex::Real sbl_damp
Transition length for blending [m].
Definition: ERF_PBLHeight.H:353
MYNNPBLH::theta_incr_water
static constexpr amrex::Real theta_incr_water
Theta increase determining the capping inversion height over water [K].
Definition: ERF_PBLHeight.H:351
MYNNPBLH::thetamin_height
static constexpr amrex::Real thetamin_height
Height below which minimum theta-v is determined [m].
Definition: ERF_PBLHeight.H:349
MYNNPBLH::theta_incr_land
static constexpr amrex::Real theta_incr_land
Theta increase determining the capping inversion height over land [K].
Definition: ERF_PBLHeight.H:350
Here is the call graph for this function:

Member Data Documentation

◆ sbl_damp

MYNNPBLH::sbl_damp = amrex::Real(400.0)
staticconstexpr

Transition length for blending [m].

Referenced by compute_pblh().

◆ sbl_lim

MYNNPBLH::sbl_lim = amrex::Real(200.0)
staticconstexpr

Upper limit of the stable boundary layer height [m].

Referenced by compute_pblh().

◆ theta_incr_land

MYNNPBLH::theta_incr_land = amrex::Real(1.25)
staticconstexpr

Theta increase determining the capping inversion height over land [K].

Referenced by compute_pblh().

◆ theta_incr_water

MYNNPBLH::theta_incr_water = one
staticconstexpr

Theta increase determining the capping inversion height over water [K].

Referenced by compute_pblh().

◆ thetamin_height

MYNNPBLH::thetamin_height = amrex::Real(200.0)
staticconstexpr

Height below which minimum theta-v is determined [m].

Referenced by compute_pblh().

The documentation for this struct was generated from the following file: