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  // Find minimum thetav in the surface layer (this updates
122  // ghost cells, too)
123  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
124  {
125  int ii = amrex::max(amrex::min(i,imax),imin);
126  int jj = amrex::max(amrex::min(j,jmax),jmin);
127 
128  amrex::Real min_thv = min_thv_arr(i,j,0);
129  for (int k(klo); k <= khi; ++k) {
130  if (zphys_arr(ii,jj,k) < thetamin_height) {
131  amrex::Real thv = GetThetav(i, j, k, cons_arr, moisture_indices);
132  min_thv = amrex::min(min_thv, thv);
133  }
134  }
135  min_thv_arr(i,j,0) = min_thv;
136  });
137 
138  // This depends on TileNoZ and k increasing monotonically
139  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
140  {
141  int ii = amrex::max(amrex::min(i,imax),imin);
142  int jj = amrex::max(amrex::min(j,jmax),jmin);
143 
144  const amrex::Real min_thv = min_thv_arr(i,j,0);
145  const int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
146 
147  amrex::Real zi = pblh_arr(i,j,0);
148  amrex::Real zi_tke = pblh_tke_arr(i,j,0);
149 
150  for (int k(klo); k <= khi; ++k)
151  {
152  if (zi == 0)
153  {
154  //
155  // Find PBL height based on thetav increase (best for CBLs)
156  //
157  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
158  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
159 
160  if (is_land && (thv1 >= min_thv + theta_incr_land)
161  && (thv < min_thv + theta_incr_land))
162  {
163  // Interpolate to get lowest height where theta = min_theta + theta_incr
164  zi = zphys_arr(ii,jj,k)
165  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
166  * (min_thv + theta_incr_land - thv);
167  }
168  else if (!is_land && (thv1 >= min_thv + theta_incr_water)
169  && (thv < min_thv + theta_incr_water))
170  {
171  // Interpolate to get lowest height where theta = min_theta + theta_incr
172  zi = zphys_arr(ii,jj,k)
173  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
174  * (min_thv + theta_incr_water - thv);
175  }
176  }
177  if (zi_tke == 0)
178  {
179  //
180  // Find PBL height based on TKE (for SBLs only)
181  //
182  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
183  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
184  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
185  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
186  amrex::Real TKEeps = amrex::Real(0.05) * maxtke;
187  TKEeps = amrex::max(TKEeps, amrex::Real(0.02)); // min val from WRF
188 
189  if ((tke1 <= TKEeps) && (tke > TKEeps))
190  {
191  // Interpolate to get lowest height where TKE -> 0
192  zi_tke = zphys_arr(ii,jj,k)
193  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(tke1-tke)
194  * (TKEeps - tke);
195  }
196  }
197  // Both criteria are satisfied; the remaining levels are no-ops
198  if ((zi != 0) && (zi_tke != 0)) { break; }
199  }
200 
201  pblh_arr(i,j,0) = zi;
202  pblh_tke_arr(i,j,0) = zi_tke;
203  });
204  }
205  else
206  // -----------------------------------------------------
207  // NO terrain
208  // -----------------------------------------------------
209  {
210  const amrex::Real dz_no_terrain = geom.CellSize(2);
211 
212  // Find minimum thetav in the surface layer (this updates
213  // ghost cells, too)
214  // - box size is known a priori
215  int kmax = static_cast<int>(thetamin_height / dz_no_terrain);
216  AMREX_ASSERT(kmax > 0);
217  const int khi_low = amrex::min(kmax, khi);
218 
219  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
220  {
221  amrex::Real min_thv = min_thv_arr(i,j,0);
222  for (int k(klo); k <= khi_low; ++k) {
223  amrex::Real thv = GetThetav(i, j, k, cons_arr, moisture_indices);
224  min_thv = amrex::min(min_thv, thv);
225  }
226  min_thv_arr(i,j,0) = min_thv;
227  });
228 
229  // This depends on TileNoZ and k increasing monotonically
230  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
231  {
232  const amrex::Real min_thv = min_thv_arr(i,j,0);
233  const int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
234 
235  amrex::Real zi = pblh_arr(i,j,0);
236  amrex::Real zi_tke = pblh_tke_arr(i,j,0);
237 
238  for (int k(klo); k <= khi; ++k)
239  {
240  if (zi == 0)
241  {
242  //
243  // Find PBL height based on thetav increase (best for CBLs)
244  //
245  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
246  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
247 
248  if (is_land && (thv1 >= min_thv + theta_incr_land)
249  && (thv < min_thv + theta_incr_land))
250  {
251  // Interpolate to get lowest height where theta = min_theta + theta_incr
252  zi = (k+myhalf)*dz_no_terrain
253  + dz_no_terrain/(thv1-thv)
254  * (min_thv + theta_incr_land - thv);
255  }
256  else if (!is_land && (thv1 >= min_thv + theta_incr_water)
257  && (thv < min_thv + theta_incr_water))
258  {
259  // Interpolate to get lowest height where theta = min_theta + theta_incr
260  zi = (k+myhalf)*dz_no_terrain
261  + dz_no_terrain/(thv1-thv)
262  * (min_thv + theta_incr_water - thv);
263  }
264  }
265  if (zi_tke == 0)
266  {
267  //
268  // Find PBL height based on TKE (for SBLs only)
269  //
270  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
271  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
272  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
273  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
274  amrex::Real TKEeps = amrex::Real(0.05) * maxtke;
275  TKEeps = amrex::max(TKEeps, amrex::Real(0.02)); // min val from WRF
276 
277  if ((tke1 <= TKEeps) && (tke > TKEeps))
278  {
279  // Interpolate to get lowest height where TKE -> 0
280  zi_tke = (k+myhalf)*dz_no_terrain
281  + dz_no_terrain/(tke1-tke) * (TKEeps - tke);
282  }
283  }
284  // Both criteria are satisfied; the remaining levels are no-ops
285  if ((zi != 0) && (zi_tke != 0)) { break; }
286  }
287 
288  pblh_arr(i,j,0) = zi;
289  pblh_tke_arr(i,j,0) = zi_tke;
290  });
291  }
292  }// MFIter
293 
294  //
295  // Calculate hybrid PBL height
296  //
297  for (amrex::MFIter mfi(*pblh); mfi.isValid(); ++mfi)
298  {
299  const auto cons_arr = cons.const_array(mfi);
300  auto pblh_tke_arr = pblh_tke.array(mfi);
301  auto pblh_arr = pblh->array(mfi);
302 
303  amrex::Box gtbx = mfi.growntilebox();
304  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
305  {
306  //
307  // Clip PBLH_TKE to more realistic values
308  //
309  // Note from WRF MYNN-EDMF: TKE-based PBLH can be very large in cells
310  // with convective precipitation (> 8km!), so an artificial limit is
311  // imposed to not let PBLH_TKE exceed the theta_v-based PBL height
312  // +/- 350 m. This has no impact on 98-99% of the domain, but is the
313  // simplest patch that adequately addresses these extremely large
314  // PBLHs.
315  amrex::Real zi = pblh_arr(i,j,0);
316  pblh_tke_arr(i,j,0) = amrex::max(
317  amrex::min(pblh_tke_arr(i,j,0), zi+amrex::Real(350.)),
318  amrex::max(zi-amrex::Real(350.), amrex::Real(10.)));
319 
320  //
321  // Finally, blend between the two PBLH estimates
322  //
323  amrex::Real maxqke = two * cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
324  if (maxqke > amrex::Real(0.05)) {
325  amrex::Real wt = myhalf*std::tanh((zi - sbl_lim)/sbl_damp) + myhalf;
326  pblh_arr(i,j,0) = (one-wt)*pblh_tke_arr(i,j,0) + wt*pblh_arr(i,j,0);
327  }
328  });
329  }//MFIter
330  }
bogus_large_value
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:17
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:41
klo
const int klo
Definition: ERF_InitCustomPert_ABL.H:75
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
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:74
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:31
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
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:50
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
WDM6SedNodeScratch::zi
@ zi
Definition: ERF_AdvanceWDM6.cpp:276
MYNNPBLH::sbl_lim
static constexpr amrex::Real sbl_lim
Upper limit of the stable boundary layer height [m].
Definition: ERF_PBLHeight.H:347
MYNNPBLH::sbl_damp
static constexpr amrex::Real sbl_damp
Transition length for blending [m].
Definition: ERF_PBLHeight.H:348
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:346
MYNNPBLH::thetamin_height
static constexpr amrex::Real thetamin_height
Height below which minimum theta-v is determined [m].
Definition: ERF_PBLHeight.H:344
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:345
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: