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