ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Public Member Functions | Static Public Attributes | List of all members
MYNNPBLH Struct Reference

#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 = 200.0
 
static constexpr amrex::Real theta_incr_land = 1.25
 
static constexpr amrex::Real theta_incr_water = 1.0
 
static constexpr amrex::Real sbl_lim = 200.0
 
static constexpr amrex::Real sbl_damp = 400.0
 

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
28  {
29 #if 0
30  // NOTE: Cannot use ReduceToPlane because it clips the box to the
31  // validbox only, i.e., lateral ghost cells aren't updated
32  int dir = 2; // z
33  auto const& cons_arrs = cons.const_arrays();
34  auto thetav_min = amrex::ReduceToPlane<amrex::ReduceOpMin,amrex::Real>(dir, bxlow, cons,
35  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> amrex::Real
36  {
37  return GetThetav(i,j,k,cons_arrs[box_no],moisture_indices);
38  });
39 #endif
40 
41  // Create 2D multifabs like pblh
42  auto const& ba = pblh->boxArray();
43  auto const& dm = pblh->DistributionMap();
44  auto const& ng = pblh->nGrowVect();
45 
46  amrex::MultiFab min_thetav(ba,dm,1,ng);
47  min_thetav.setVal(1.E34);
48 
49  amrex::MultiFab pblh_tke(ba,dm,1,ng);
50  pblh_tke.setVal(0);
51 
52  pblh->setVal(0);
53 
54  // Now, loop over columns...
55  for (amrex::MFIter mfi(cons,TileNoZ()); mfi.isValid(); ++mfi)
56  {
57  const amrex::Box& domain = geom.Domain();
58  amrex::Box gtbx = mfi.growntilebox();
59  gtbx.setSmall(2,domain.smallEnd(2)); // don't loop over ghost cells
60  gtbx.setBig(2,domain.bigEnd(2)); // in z
61 
62  auto min_thv_arr = min_thetav.array(mfi);
63  auto pblh_arr = pblh->array(mfi);
64  auto pblh_tke_arr = pblh_tke.array(mfi);
65 
66  const auto cons_arr = cons.const_array(mfi);
67  const auto lmask_arr = (lmask) ? lmask->const_array(mfi) : amrex::Array4<int> {};
68 
69  // -----------------------------------------------------
70  // WITH terrain/grid stretching
71  // -----------------------------------------------------
72  if (z_phys_cc)
73  {
74  const auto zphys_arr = z_phys_cc->const_array(mfi);
75 
76  // Need to sort out ghost cell differences (z_phys_cc has ng=1)
77  int imin = lbound(zphys_arr).x;
78  int jmin = lbound(zphys_arr).y;
79  int imax = ubound(zphys_arr).x;
80  int jmax = ubound(zphys_arr).y;
81 
82  // Find minimum thetav in the surface layer (this updates
83  // ghost cells, too)
84  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
85  {
86  int ii = amrex::max(amrex::min(i,imax),imin);
87  int jj = amrex::max(amrex::min(j,jmax),jmin);
88 
89  if (zphys_arr(ii,jj,k) < thetamin_height) {
90  amrex::Real thv = GetThetav(i, j, k, cons_arr, moisture_indices);
91  if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
92  }
93  });
94 
95  // This depends on TileNoZ and k increasing monotonically
96  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
97  {
98  int ii = amrex::max(amrex::min(i,imax),imin);
99  int jj = amrex::max(amrex::min(j,jmax),jmin);
100 
101  if (pblh_arr(i,j,0) == 0)
102  {
103  //
104  // Find PBL height based on thetav increase (best for CBLs)
105  //
106  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
107  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
108 
109  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
110  if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)
111  && (thv < min_thv_arr(i,j,0) + theta_incr_land))
112  {
113  // Interpolate to get lowest height where theta = min_theta + theta_incr
114  pblh_arr(i,j,0) = zphys_arr(ii,jj,k)
115  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
116  * (min_thv_arr(i,j,0) + theta_incr_land - thv);
117  }
118  else if (!is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_water)
119  && (thv < min_thv_arr(i,j,0) + theta_incr_water))
120  {
121  // Interpolate to get lowest height where theta = min_theta + theta_incr
122  pblh_arr(i,j,0) = zphys_arr(ii,jj,k)
123  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
124  * (min_thv_arr(i,j,0) + theta_incr_water - thv);
125  }
126  }
127  if (pblh_tke_arr(i,j,0) == 0)
128  {
129  //
130  // Find PBL height based on TKE (for SBLs only)
131  //
132  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
133  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
134  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
135  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
136  amrex::Real TKEeps = 0.05 * maxtke;
137  TKEeps = amrex::max(TKEeps, 0.02); // min val from WRF
138 
139  if ((tke1 <= TKEeps) && (tke > TKEeps))
140  {
141  // Interpolate to get lowest height where TKE -> 0
142  pblh_tke_arr(i,j,0) = zphys_arr(ii,jj,k)
143  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(tke1-tke)
144  * (TKEeps - tke);
145  }
146  }
147  });
148  }
149  else
150  // -----------------------------------------------------
151  // NO terrain
152  // -----------------------------------------------------
153  {
154  const amrex::Real dz_no_terrain = geom.CellSize(2);
155 
156  // Find minimum thetav in the surface layer (this updates
157  // ghost cells, too)
158  // - box size is known a priori
159  int kmax = static_cast<int>(thetamin_height / dz_no_terrain);
160  AMREX_ASSERT(kmax > 0);
161  amrex::Box gtbxlow = gtbx;
162  gtbxlow.setBig(2,kmax);
163 
164  ParallelFor(gtbxlow, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
165  {
166  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
167  if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
168  });
169 
170  // This depends on TileNoZ and k increasing monotonically
171  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
172  {
173  if (pblh_arr(i,j,0) == 0)
174  {
175  //
176  // Find PBL height based on thetav increase (best for CBLs)
177  //
178  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
179  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
180 
181  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
182  if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)
183  && (thv < min_thv_arr(i,j,0) + theta_incr_land))
184  {
185  // Interpolate to get lowest height where theta = min_theta + theta_incr
186  pblh_arr(i,j,0) = (k+0.5)*dz_no_terrain
187  + dz_no_terrain/(thv1-thv)
188  * (min_thv_arr(i,j,0) + theta_incr_land - thv);
189  }
190  else if (!is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_water)
191  && (thv < min_thv_arr(i,j,0) + theta_incr_water))
192  {
193  // Interpolate to get lowest height where theta = min_theta + theta_incr
194  pblh_arr(i,j,0) = (k+0.5)*dz_no_terrain
195  + dz_no_terrain/(thv1-thv)
196  * (min_thv_arr(i,j,0) + theta_incr_water - thv);
197  }
198  }
199  if (pblh_tke_arr(i,j,0) == 0)
200  {
201  //
202  // Find PBL height based on TKE (for SBLs only)
203  //
204  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
205  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
206  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
207  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
208  amrex::Real TKEeps = 0.05 * maxtke;
209  TKEeps = amrex::max(TKEeps, 0.02); // min val from WRF
210 
211  if ((tke1 <= TKEeps) && (tke > TKEeps))
212  {
213  // Interpolate to get lowest height where TKE -> 0
214  pblh_tke_arr(i,j,0) = (k+0.5)*dz_no_terrain
215  + dz_no_terrain/(tke1-tke) * (TKEeps - tke);
216  }
217  }
218  });
219  }
220  }// MFIter
221 
222  //
223  // Calculate hybrid PBL height
224  //
225  for (amrex::MFIter mfi(*pblh); mfi.isValid(); ++mfi)
226  {
227  const auto cons_arr = cons.const_array(mfi);
228  auto pblh_tke_arr = pblh_tke.array(mfi);
229  auto pblh_arr = pblh->array(mfi);
230 
231  amrex::Box gtbx = mfi.growntilebox();
232  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
233  {
234  //
235  // Clip PBLH_TKE to more realistic values
236  //
237  // Note from WRF MYNN-EDMF: TKE-based PBLH can be very large in cells
238  // with convective precipiation (> 8km!), so an artificial limit is
239  // imposed to not let PBLH_TKE exceed the theta_v-based PBL height
240  // +/- 350 m. This has no impact on 98-99% of the domain, but is the
241  // simplest patch that adequately addresses these extremely large
242  // PBLHs.
243  amrex::Real zi = pblh_arr(i,j,0);
244  pblh_tke_arr(i,j,0) = amrex::max(
245  amrex::min(pblh_tke_arr(i,j,0), zi+350.),
246  amrex::max(zi-350., 10.));
247 
248  //
249  // Finally, blend between the two PBLH estimates
250  //
251  amrex::Real maxqke = 2.0 * cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
252  if (maxqke > 0.05) {
253  amrex::Real wt = 0.5*std::tanh((zi - sbl_lim)/sbl_damp) + 0.5;
254  pblh_arr(i,j,0) = (1.-wt)*pblh_tke_arr(i,j,0) + wt*pblh_arr(i,j,0);
255  }
256  });
257  }//MFIter
258  }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
GetThetav
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetav(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:72
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:48
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:140
MYNNPBLH::sbl_lim
static constexpr amrex::Real sbl_lim
Definition: ERF_PBLHeight.H:263
MYNNPBLH::theta_incr_land
static constexpr amrex::Real theta_incr_land
Definition: ERF_PBLHeight.H:261
MYNNPBLH::thetamin_height
static constexpr amrex::Real thetamin_height
Definition: ERF_PBLHeight.H:260
MYNNPBLH::theta_incr_water
static constexpr amrex::Real theta_incr_water
Definition: ERF_PBLHeight.H:262
MYNNPBLH::sbl_damp
static constexpr amrex::Real sbl_damp
Definition: ERF_PBLHeight.H:264
Here is the call graph for this function:

Member Data Documentation

◆ sbl_damp

constexpr amrex::Real MYNNPBLH::sbl_damp = 400.0
staticconstexpr

Referenced by compute_pblh().

◆ sbl_lim

constexpr amrex::Real MYNNPBLH::sbl_lim = 200.0
staticconstexpr

Referenced by compute_pblh().

◆ theta_incr_land

constexpr amrex::Real MYNNPBLH::theta_incr_land = 1.25
staticconstexpr

Referenced by compute_pblh().

◆ theta_incr_water

constexpr amrex::Real MYNNPBLH::theta_incr_water = 1.0
staticconstexpr

Referenced by compute_pblh().

◆ thetamin_height

constexpr amrex::Real MYNNPBLH::thetamin_height = 200.0
staticconstexpr

Referenced by compute_pblh().

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