ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_MetgridUtils.H
Go to the documentation of this file.
1 /**
2  * \file ERF_MetgridUtils.H
3  */
4 #ifndef ERF_METGRIDUTIL_H_
5 #define ERF_METGRIDUTIL_H_
6 
7 #include <AMReX_Math.H>
8 
9 #include <cmath>
10 #include "ERF_Constants.H"
11 #include <limits>
12 
13 #include <ERF.H>
14 #include <ERF_EOS.H>
15 #include <ERF_SurfaceTemperature.H>
16 #include <ERF_Utils.H>
17 #include <ERF_ProbCommon.H>
18 #include <ERF_HSEUtils.H>
19 #include <ERF_BaseStateParams.H>
20 
21 /**
22  * Return the physical surface pressure used by the Metgrid initialization.
23  *
24  * A debug pressure deliberately overrides file data. Otherwise a PSFC value
25  * from the current forcing file is preferred; when PSFC is absent, invert the
26  * standard-atmosphere temperature profile at the physical terrain height.
27  */
28 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
29 bool
30 metgrid_surface_pressure (const bool debug_psfc,
31  const int flag_psfc,
32  const amrex::Real file_psfc,
33  const amrex::Real z_sfc,
34  const amrex::Real P00,
35  const amrex::Real T00,
36  const amrex::Real TLP,
37  amrex::Real& pressure)
38 {
39  pressure = std::numeric_limits<amrex::Real>::quiet_NaN();
40  if (debug_psfc) {
41  pressure = amrex::Real(100000.0);
42  return true;
43  }
44  if (flag_psfc == 1) {
45  if (!amrex::Math::isfinite(file_psfc) || file_psfc <= amrex::Real(0.0)) {
46  return false;
47  }
48  pressure = file_psfc;
49  return true;
50  }
51 
52  if (!amrex::Math::isfinite(z_sfc) || !amrex::Math::isfinite(P00) || P00 <= amrex::Real(0.0) ||
53  !amrex::Math::isfinite(T00) || T00 <= amrex::Real(0.0) ||
54  !amrex::Math::isfinite(TLP) || TLP <= amrex::Real(0.0)) {
55  return false;
56  }
57  const amrex::Real ToA = T00 / TLP;
58  const amrex::Real radicand =
59  ToA*ToA - amrex::Real(2.0) * CONST_GRAV * z_sfc / (TLP * R_d);
60  if (!amrex::Math::isfinite(radicand) || radicand < amrex::Real(0.0)) {
61  return false;
62  }
63  const amrex::Real exponent = -ToA + std::sqrt(radicand);
64  const amrex::Real result = P00 * std::exp(exponent);
65  if (!amrex::Math::isfinite(result) || result <= amrex::Real(0.0)) {
66  return false;
67  }
68  pressure = result;
69  return true;
70 }
71 
72 /**
73  * Convert an absolute Metgrid surface temperature to the canonical
74  * SurfaceLayer potential-temperature representation using the selected
75  * physical surface pressure.
76  */
77 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
78 bool
79 metgrid_surface_theta (const amrex::Real temperature,
80  const bool debug_psfc,
81  const int flag_psfc,
82  const amrex::Real file_psfc,
83  const amrex::Real z_sfc,
84  const amrex::Real P00,
85  const amrex::Real T00,
86  const amrex::Real TLP,
87  const amrex::Real rdOcp,
88  amrex::Real& theta)
89 {
90  amrex::Real surface_pressure;
91  if (!metgrid_surface_pressure(debug_psfc, flag_psfc, file_psfc, z_sfc,
92  P00, T00, TLP, surface_pressure)) {
93  return false;
94  }
95  return erf_surface_temperature::temperature_to_theta(
96  temperature, surface_pressure, rdOcp, theta);
97 }
98 
99 /**
100  * Read metgrid NetCDF data and populate ERF initialization FABs.
101  */
102 void
103 read_from_metgrid (int lev,
104  int itime,
105  const amrex::Box& domain,
106  const std::string& fname,
107  std::string& NC_dateTime,
108  double& NC_epochTime,
109  int& flag_psfc,
110  int& flag_msf,
111  int& flag_sst,
112  int& flag_tsk,
113  int& flag_lmask,
114  int& NC_nx,
115  int& NC_ny,
116  amrex::Real& NC_dx,
117  amrex::Real& NC_dy,
118  amrex::FArrayBox& NC_xvel_fab,
119  amrex::FArrayBox& NC_yvel_fab,
120  amrex::FArrayBox& NC_temp_fab,
121  amrex::FArrayBox& NC_rhum_fab,
122  amrex::FArrayBox& NC_pres_fab,
123  amrex::FArrayBox& NC_ght_fab,
124  amrex::FArrayBox& NC_hgt_fab,
125  amrex::FArrayBox& NC_psfc_fab,
126  amrex::FArrayBox& NC_msfu_fab,
127  amrex::FArrayBox& NC_msfv_fab,
128  amrex::FArrayBox& NC_msfm_fab,
129  amrex::FArrayBox& NC_sst_fab,
130  amrex::FArrayBox& NC_tsk_fab,
131  amrex::FArrayBox& NC_LAT_fab,
132  amrex::FArrayBox& NC_LON_fab,
133  amrex::IArrayBox& NC_lmask_iab,
134  amrex::Geometry& geom);
135 
136 /**
137  * Initialize nodal terrain heights from metgrid surface-height data.
138  */
139 void
140 init_terrain_from_metgrid (const bool& avg_grid_faces_to_nodes,
141  amrex::Geometry& geom,
142  amrex::MultiFab* z_phys_nd,
143  amrex::FArrayBox& NC_hgt_fab);
144 
145 /**
146  * Initialize ERF state, velocity, and boundary data from metgrid fields.
147  */
148 void
149 init_state_from_metgrid (const int lev,
150  const int itime,
151  const bool use_moisture,
152  const bool interp_theta,
153  const bool metgrid_debug_quiescent,
154  const bool metgrid_debug_isothermal,
155  const bool metgrid_debug_dry,
156  const bool metgrid_basic_linear,
157  const bool metgrid_use_below_sfc,
158  const bool metgrid_use_sfc,
159  const bool metgrid_retain_sfc,
160  const amrex::Real metgrid_proximity,
161  const int metgrid_order,
162  const int metgrid_metgrid_force_sfc_k,
163  const amrex::Real l_rdOcp,
164  amrex::Box& tbxc,
165  amrex::Box& tbxu,
166  amrex::Box& tbxv,
167  amrex::Box& tbxw,
168  amrex::FArrayBox& state_fab,
169  amrex::FArrayBox& x_vel_fab,
170  amrex::FArrayBox& y_vel_fab,
171  amrex::FArrayBox& z_vel_fab,
172  amrex::FArrayBox& z_phys_nd_fab,
173  const amrex::FArrayBox& NC_ght_fab,
174  const amrex::FArrayBox& NC_xvel_fab,
175  const amrex::FArrayBox& NC_yvel_fab,
176  const amrex::FArrayBox& NC_temp_fab,
177  const amrex::FArrayBox& NC_rhum_fab,
178  const amrex::FArrayBox& NC_pres_fab,
179  amrex::FArrayBox& tmp_src_fab,
180  amrex::FArrayBox& tmp_dst_fab,
181  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xlo,
182  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xhi,
183  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_ylo,
184  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_yhi,
185  const amrex::Array4<const int>& mask_c_arr,
186  const amrex::Array4<const int>& mask_u_arr,
187  const amrex::Array4<const int>& mask_v_arr);
188 
189 /**
190  * Initialize map-scale factors from metgrid fields or defaults.
191  */
192 void
193 init_msfs_from_metgrid (const bool metgrid_debug_msf,
194  amrex::FArrayBox& msfu_fab,
195  amrex::FArrayBox& msfv_fab,
196  amrex::FArrayBox& msfm_fab,
197  const int& flag_msf,
198  amrex::FArrayBox& NC_MSFU_fab,
199  amrex::FArrayBox& NC_MSFV_fab,
200  amrex::FArrayBox& NC_MSFM_fab);
201 
202 
203 /**
204  * Initialize hydrostatic base-state quantities from metgrid data.
205  */
206 void
207 init_base_state_from_metgrid (const bool use_moisture,
208  const bool metgrid_debug_psfc,
209  const amrex::Real l_rdOcp,
210  const amrex::Box& domain,
211  const amrex::Box& valid_bx,
212  const int& flag_psfc,
213  amrex::FArrayBox& state_fab,
214  amrex::FArrayBox& r_hse_fab,
215  amrex::FArrayBox& p_hse_fab,
216  amrex::FArrayBox& pi_hse_fab,
217  amrex::FArrayBox& th_hse_fab,
218  amrex::FArrayBox& qv_hse_fab,
219  amrex::FArrayBox& z_phys_nd_fab,
220  amrex::FArrayBox& z_phys_cc_fab,
221  const amrex::FArrayBox& NC_psfc_fab,
222  int* base_state_column_failed,
223  const BaseStateParams& bsp);
224 
225 /**
226  * Interpolate a value with a Lagrange polynomial stencil.
227  */
228 AMREX_FORCE_INLINE
229 AMREX_GPU_DEVICE
230 void
231 lagrange_interp (const int& order,
232  amrex::Real* x,
233  amrex::Real* y,
234  amrex::Real& new_x,
235  amrex::Real& new_y)
236 {
237  // Interpolation using Lagrange polynomials.
238  // P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
239  // where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
240  // ---------------------------------------------
241  // (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
242  amrex::Real Px = zero;
243  for (int i=0; i <= order; i++) {
244  amrex::Real n = one;
245  amrex::Real d = one;
246  for (int k=0; k <= order; k++) {
247  if (k == i) continue;
248  n *= new_x-x[k];
249  d *= x[i]-x[k];
250  }
251  if (d != zero) {
252  Px += y[i]*n/d;
253  }
254  }
255  new_y = Px;
256 }
257 
258 /**
259  * Prepare interpolation stencils and apply Lagrange interpolation or extrapolation.
260  *
261  * Both orig_npts and new_npts are point counts, not last indices; valid indices
262  * into the orig_* arrays are 0 through orig_npts-1, and into the new_* arrays
263  * are 0 through new_npts-1.
264  */
265 AMREX_FORCE_INLINE
266 AMREX_GPU_DEVICE
267 void
268 lagrange_setup (char var_type,
269  const bool& exp_interp,
270  const int& orig_npts,
271  const int& new_npts,
272  const int& order,
273  const int& i,
274  const int& j,
275  amrex::Real* orig_x_z,
276  amrex::Real* orig_x_p,
277  amrex::Real* orig_y,
278  amrex::Real* new_x_z,
279  amrex::Real* new_x_p,
280  amrex::Real* new_y)
281 {
282 
283  amrex::Real CRC_const1 = amrex::Real(11880.516); // m
284  amrex::Real CRC_const2 = amrex::Real(0.1902632);
285  amrex::Real CRC_const3 = amrex::Real(0.0065); // K km-1
286 
287  amrex::ignore_unused(i,j);
288 #ifndef AMREX_USE_GPU
289  bool debug = false;
290 #endif
291 
292  for (int new_k=0; new_k < new_npts; new_k++) {
293 #ifndef AMREX_USE_GPU
294  if (debug) amrex::Print() << "new_k=" << new_k;
295 #endif
296  // Determine if interpolating or extrapolating.
297  // If interpolating, find bounding x values and store the indices.
298  bool extrapolating = true;
299  int kl, kr;
300  for (int ko=0; ko < orig_npts-1; ko++) {
301  amrex::Real a = new_x_z[new_k]-orig_x_z[ko];
302  amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1];
303  if (a*b <= zero) {
304  kl = ko;
305  kr = ko+1;
306  extrapolating = false;
307  break;
308  }
309  }
310 
311  if (extrapolating) {
312  if (var_type == 'T') {
313  // Assume a standard atmosphere -amrex::Real(6.5) K km-1 lapse rate.
314  // Comparable to the WRF default, t_extrap_type=two
315  amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[0];
316  amrex::Real avg_of_extrap_p = myhalf*(new_x_p[new_k]+orig_x_p[0]);
317  amrex::Real temp_extrap_starting_point = orig_y[0]*std::pow(orig_x_p[0]/amrex::Real(100000.0), RdoCp);
318  amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/amrex::Real(100.0), CRC_const2-one);
319  amrex::Real dZ = dZdP*(depth_of_extrap_in_p/amrex::Real(100.0));
320  amrex::Real dT = dZ*CRC_const3;
321  new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(amrex::Real(100000.0)/new_x_p[new_k], RdoCp);
322  } else {
323  // Use a constant value below ground.
324  // Comparable to the WRF default, extrap_type=two
325  new_y[new_k] = orig_y[0];
326  }
327  continue;
328  }
329 
330  if (order%2 != 0) {
331  if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_npts-1)) {
332  // Odd order interpolation.
333  int ksta = kl-(((order+1)/2)-1);
334  int kend = ksta+order;
335 #ifndef AMREX_USE_GPU
336  int ksize = kend-ksta;
337  if (debug) amrex::Print() << " (1a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
338 #endif
339  amrex::Real new_x;
340  amrex::GpuArray<amrex::Real,9> orig_x_sub;
341  amrex::Real* orig_x_sub_p = orig_x_sub.data();
342  if (exp_interp) {
343  new_x = new_x_p[new_k];
344  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
345  } else {
346  new_x = new_x_z[new_k];
347  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
348  }
349  amrex::GpuArray<amrex::Real,9> orig_y_sub;
350  amrex::Real* orig_y_sub_p = orig_y_sub.data();
351  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
352 #ifndef AMREX_USE_GPU
353  if (debug) {
354  amrex::Print() << " orig_x_sub_p = [";
355  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
356  amrex::Print() << "]" << std::endl;
357  amrex::Print() << " orig_y_sub_p = [";
358  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
359  amrex::Print() << "]" << std::endl;
360  }
361 #endif
362  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
363  } else {
364  // Interpolation stencil is too big due to proximity to the surface or model top.
365  // We have bounding points so resort to a simple linear interpolation.
366  int ksta = kl;
367  int kend = kr;
368 #ifndef AMREX_USE_GPU
369  int ksize = kend-ksta+1;
370  if (debug) amrex::Print() << " (1b) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
371 #endif
372  amrex::Real new_x;
373  amrex::GpuArray<amrex::Real,2> orig_x_sub;
374  amrex::GpuArray<amrex::Real,2> orig_y_sub;
375  amrex::Real* orig_x_sub_p = orig_x_sub.data();
376  amrex::Real* orig_y_sub_p = orig_y_sub.data();
377  if (exp_interp) {
378  new_x = new_x_p[new_k];
379  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
380  } else {
381  new_x = new_x_z[new_k];
382  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
383  }
384  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
385 #ifndef AMREX_USE_GPU
386  if (debug) {
387  amrex::Print() << " orig_x_sub_p = [";
388  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
389  amrex::Print() << "]" << std::endl;
390  amrex::Print() << " orig_y_sub_p = [";
391  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
392  amrex::Print() << "]" << std::endl;
393  }
394 #endif
395  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
396  }
397  } else if (order%2 == 0) {
398  if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_npts-1)) {
399  // Even order interpolation.
400  // Interpolate twice. Offset the range by 1 the 2nd time. Average the results.
401  amrex::Real new_y_l, new_y_r;
402  {
403  int ksta = kl-(order/2-1);
404  int kend = ksta+order;
405 #ifndef AMREX_USE_GPU
406  int ksize = kend-ksta;
407  if (debug) amrex::Print() << " (2a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
408 #endif
409  amrex::Real new_x;
410  amrex::GpuArray<amrex::Real,10> orig_x_sub;
411  amrex::GpuArray<amrex::Real,10> orig_y_sub;
412  amrex::Real* orig_x_sub_p = orig_x_sub.data();
413  amrex::Real* orig_y_sub_p = orig_y_sub.data();
414  if (exp_interp) {
415  new_x = new_x_p[new_k];
416  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
417  } else {
418  new_x = new_x_z[new_k];
419  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
420  }
421  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
422 #ifndef AMREX_USE_GPU
423  if (debug) {
424  amrex::Print() << " orig_x_sub_p = [";
425  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
426  amrex::Print() << "]" << std::endl;
427  amrex::Print() << " orig_y_sub_p = [";
428  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
429  amrex::Print() << "]" << std::endl;
430  }
431 #endif
432  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_l);
433  }
434  {
435  int ksta = kl-order/2;
436  int kend = ksta+order;
437 #ifndef AMREX_USE_GPU
438  int ksize = kend-ksta;
439  if (debug) amrex::Print() << "new_k=" << new_k << " (2b) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
440 #endif
441  amrex::Real new_x;
442  amrex::GpuArray<amrex::Real,10> orig_x_sub;
443  amrex::GpuArray<amrex::Real,10> orig_y_sub;
444  amrex::Real* orig_x_sub_p = orig_x_sub.data();
445  amrex::Real* orig_y_sub_p = orig_y_sub.data();
446  if (exp_interp) {
447  new_x = new_x_p[new_k];
448  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
449  } else {
450  new_x = new_x_z[new_k];
451  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
452  }
453  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
454 #ifndef AMREX_USE_GPU
455  if (debug) {
456  amrex::Print() << " orig_x_sub_p = [";
457  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
458  amrex::Print() << "]" << std::endl;
459  amrex::Print() << " orig_y_sub_p = [";
460  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
461  amrex::Print() << "]" << std::endl;
462  }
463 #endif
464  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_r);
465  }
466  new_y[new_k] = myhalf*(new_y_l+new_y_r);
467  } else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_npts-1)) {
468  int ksta = kl-(order/2-1);
469  int kend = ksta+order;
470 #ifndef AMREX_USE_GPU
471  int ksize = kend-ksta;
472  if (debug) amrex::Print() << " (3) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
473 #endif
474  amrex::Real new_x;
475  amrex::GpuArray<amrex::Real,10> orig_x_sub;
476  amrex::GpuArray<amrex::Real,10> orig_y_sub;
477  amrex::Real* orig_x_sub_p = orig_x_sub.data();
478  amrex::Real* orig_y_sub_p = orig_y_sub.data();
479  if (exp_interp) {
480  new_x = new_x_p[new_k];
481  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
482  } else {
483  new_x = new_x_z[new_k];
484  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
485  }
486  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
487 #ifndef AMREX_USE_GPU
488  if (debug) {
489  amrex::Print() << " orig_x_sub_p = [";
490  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
491  amrex::Print() << "]" << std::endl;
492  amrex::Print() << " orig_y_sub_p = [";
493  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
494  amrex::Print() << "]" << std::endl;
495  }
496 #endif
497  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
498  } else if ((kl-order/2 >= 0) && (kr+order/2-1 <= orig_npts-1)) {
499  int ksta = kl-order/2;
500  int kend = ksta+order;
501 #ifndef AMREX_USE_GPU
502  int ksize = kend-ksta;
503  if (debug) amrex::Print() << " (4) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
504 #endif
505  amrex::Real new_x;
506  amrex::GpuArray<amrex::Real,10> orig_x_sub;
507  amrex::GpuArray<amrex::Real,10> orig_y_sub;
508  amrex::Real* orig_x_sub_p = orig_x_sub.data();
509  amrex::Real* orig_y_sub_p = orig_y_sub.data();
510  if (exp_interp) {
511  new_x = new_x_p[new_k];
512  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
513  } else {
514  new_x = new_x_z[new_k];
515  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
516  }
517  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
518 #ifndef AMREX_USE_GPU
519  if (debug) {
520  amrex::Print() << " orig_x_sub_p = [";
521  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
522  amrex::Print() << "]" << std::endl;
523  amrex::Print() << " orig_y_sub_p = [";
524  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
525  amrex::Print() << "]" << std::endl;
526  }
527 #endif
528  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
529  } else {
530  // Linear interpolation.
531  int ksta = kl;
532  int kend = kr;
533 #ifndef AMREX_USE_GPU
534  int ksize = kend-ksta+1;
535  if (debug) amrex::Print() << " (5) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
536 #endif
537  amrex::Real new_x;
538  amrex::GpuArray<amrex::Real,2> orig_x_sub;
539  amrex::GpuArray<amrex::Real,2> orig_y_sub;
540  amrex::Real* orig_x_sub_p = orig_x_sub.data();
541  amrex::Real* orig_y_sub_p = orig_y_sub.data();
542  if (exp_interp) {
543  new_x = new_x_p[new_k];
544  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
545  } else {
546  new_x = new_x_z[new_k];
547  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
548  }
549  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
550 #ifndef AMREX_USE_GPU
551  if (debug) {
552  amrex::Print() << " orig_x_sub_p = [";
553  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
554  amrex::Print() << "]" << std::endl;
555  amrex::Print() << " orig_y_sub_p = [";
556  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
557  amrex::Print() << "]" << std::endl;
558  }
559 #endif
560  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
561  }
562  }
563 #ifndef AMREX_USE_GPU
564  if (debug) amrex::Print() << " new_y[" << new_k << "]=" << new_y[new_k] << std::endl;
565 #endif
566  }
567 }
568 
569 /**
570  * Compute isothermal pressure at height z.
571  */
572 AMREX_FORCE_INLINE
573 AMREX_GPU_DEVICE
574 void
575 calc_p_isothermal (const amrex::Real& z,
576  amrex::Real& p)
577 {
578  p = p_0*std::exp(-CONST_GRAV*z/(amrex::Real(290.0)*R_d));
579 }
580 
581 /**
582  * Interpolate one metgrid column onto ERF vertical levels with quality control.
583  */
584 AMREX_FORCE_INLINE
585 AMREX_GPU_DEVICE
586 void
587 interpolate_column_metgrid (const bool& metgrid_use_below_sfc,
588  const bool& metgrid_use_sfc,
589  const bool& exp_interp,
590  const bool& metgrid_retain_sfc,
591  const amrex::Real& metgrid_proximity,
592  const int& metgrid_order,
593  const int& metgrid_force_sfc_k,
594  const int& i,
595  const int& j,
596  const int& kmax,
597  const int& src_comp,
598  const int& dst_comp,
599  char var_type,
600  char stag,
601  const amrex::Array4<amrex::Real const>& orig_z_full,
602  const amrex::Array4<amrex::Real const>& orig_data,
603  const amrex::Array4<amrex::Real const>& new_z_full,
604  const amrex::Array4<amrex::Real>& new_data_full)
605 {
606  // Here we closely follow WRF's vert_interp from
607  // dyn_em/module_initialize_real.F, although changes have been
608  // made to accommodate interpolation relative to height instead of
609  // pressure.
610  int imin_orig = amrex::lbound(amrex::Box(orig_data)).x;
611  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
612  int jmin_orig = amrex::lbound(amrex::Box(orig_data)).y;
613  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
614  // The number of origin levels, 0 .. ubound: the loops below run k < kmax_orig, so taking
615  // ubound itself here would leave the top origin level out
616  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z + 1;
617  int kmax_new = kmax + 1;
618 
619  AMREX_ASSERT(kmax_orig <= 256);
620  AMREX_ASSERT(kmax_new < 256);
621 
622  amrex::GpuArray<amrex::Real,256> new_z;
623  amrex::GpuArray<amrex::Real,256> new_p;
624  amrex::GpuArray<amrex::Real,256> new_data;
625  amrex::Real* new_z_p = new_z.data();
626  amrex::Real* new_p_p = new_p.data();
627  amrex::Real* new_data_p = new_data.data();
628  for (int k=0; k < kmax_new; k++) {
629  if (stag == 'X') {
630  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i,j+1,k)+new_z_full(i,j,k+1)+new_z_full(i,j+1,k+1));
631  } else if (stag == 'Y') {
632  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i+1,j,k)+new_z_full(i,j,k+1)+new_z_full(i+1,j,k+1));
633  } else if (stag == 'M') {
634  new_z_p[k] = amrex::Real(0.125)*(new_z_full(i,j,k )+new_z_full(i,j+1,k )+new_z_full(i+1,j,k )+new_z_full(i+1,j+1,k )+
635  new_z_full(i,j,k+1)+new_z_full(i,j+1,k+1)+new_z_full(i+1,j,k+1)+new_z_full(i+1,j+1,k+1));
636  }
637  calc_p_isothermal(new_z_p[k], new_p_p[k]);
638  }
639 
640  amrex::GpuArray<amrex::Real,256> orig_z;
641  amrex::Real* orig_z_p = orig_z.data();
642  for (int k=0; k < kmax_orig; k++) {
643  if (stag == 'M') {
644  orig_z_p[k] = orig_z_full(i,j,k);
645  } else if (stag == 'X') {
646  if (i <= imin_orig) {
647  orig_z_p[k] = orig_z_full(i,j,k);
648  } else if (i >= imax_orig) {
649  orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
650  } else {
651  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
652  }
653  } else if (stag == 'Y') {
654  if (j <= jmin_orig) {
655  orig_z_p[k] = orig_z_full(i,j,k);
656  } else if (j >= jmax_orig) {
657  orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
658  } else {
659  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
660  }
661  }
662  }
663 
664  // Check if the data is top-down instead of bottom-up.
665  bool flip_data_required = false;
666  if (orig_z_p[1] > orig_z_p[kmax_orig-1]) flip_data_required = true;
667  if (flip_data_required) amrex::Abort("metgrid initialization flip_data_required. Not yet implemented.");
668 
669  // Search for the first level above the surface in the origin data.
670  // This is needed since the origin model topography will be
671  // different than the topography processed by WPS.
672  int k_above_sfc = 0;
673  for (int k=1; k < kmax_orig; k++) {
674  if (orig_z_p[k] > orig_z_p[0]) {
675  k_above_sfc = k;
676  break;
677  }
678  }
679 
680  int kend_order;
681  amrex::GpuArray<amrex::Real,256> ordered_z;
682  amrex::GpuArray<amrex::Real,256> ordered_data;
683  amrex::Real* ordered_z_p = ordered_z.data();
684  amrex::Real* ordered_data_p = ordered_data.data();
685  if (k_above_sfc > 1) {
686  // The levels are not monotonically increasing in height, so
687  // we sort and then make "artistic" quality control choices.
688  int count = 0;
689 
690  // Insert levels that are below the surface.
691  for (int k=1; k < k_above_sfc; k++) {
692  ordered_z_p[count] = orig_z_p[k];
693  ordered_data_p[count] = orig_data(i,j,k,src_comp);
694  count++;
695  }
696 
697  // Check if the level that is nearest to and below the surface
698  // is "too close". If so, we'll ignore the upper level and keep
699  // the lower. Origin data is likely to be on pressure levels
700  // with higher spatial resolution near-surface, which supports
701  // the choice of eliminating levels that are "too close" in
702  // pressure-space. For simplicity, calculate delta P assuming a
703  // baroclinic atmosphere.
704  amrex::Real Pu, Pl;
705  calc_p_isothermal(orig_z_p[0], Pu);
706  calc_p_isothermal(ordered_z_p[count-1], Pl);
707  if (Pl-Pu < metgrid_proximity) {
708  count--;
709  }
710 
711  // Insert the surface level.
712  ordered_z_p[count] = orig_z_p[0];
713  ordered_data_p[count] = orig_data(i,j,0,src_comp);
714  count++;
715 
716  // Quoting WRF's comments, the next level to use is at,
717  // "... ta da, the first level above the surface. I know, wow."
718  int knext = k_above_sfc;
719  // Conditionally more strongly use the surface data by removing
720  // levels between the surface and the height corresponding to a
721  // set number of ERF levels from the surface. This forces the
722  // interpolator to use the surface data up through a number of
723  // ERF levels from the surface.
724  if (metgrid_force_sfc_k > 0) {
725  for (int k=k_above_sfc; k < kmax_orig; k++) {
726  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
727  knext = k;
728  break;
729  }
730  }
731  }
732 
733  // Check if the level that is nearest to and above the surface
734  // is "too close". If so, we'll ignore that level.
735  calc_p_isothermal(orig_z_p[knext], Pu);
736  calc_p_isothermal(ordered_z_p[count-1], Pl);
737  if (Pl-Pu < metgrid_proximity) {
738  knext++;
739  }
740 
741  // Insert levels that are above the surface.
742  for (int k=knext; k < kmax_orig; k++) {
743  ordered_z_p[count] = orig_z_p[k];
744  ordered_data_p[count] = orig_data(i,j,k,src_comp);
745  count++;
746  }
747 
748  kend_order = count;
749  } else {
750  // The surface is the lowest level in the origin data.
751 
752  // Insert the surface.
753  ordered_z_p[0] = orig_z_p[0];
754  ordered_data_p[0] = orig_data(i,j,0,src_comp);
755 
756  // Similar to above, conditionally more strongly use the
757  // surface data.
758  int count = 1;
759  int knext = count;
760  // The same ERF level as in the branch above: the surface is used through
761  // metgrid_force_sfc_k ERF levels, k = 0 .. metgrid_force_sfc_k-1
762  if (metgrid_force_sfc_k > 0) {
763  for (int k=knext; k < kmax_orig; k++) {
764  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
765  knext = k;
766  break;
767  }
768  }
769  }
770 
771  // Insert the remaining levels, again ignoring levels that are
772  // "too close" to the prior valid level.
773  for (int k=knext; k < kmax_orig; k++) {
774  amrex::Real Pu, Pl;
775  calc_p_isothermal(orig_z_p[k], Pu);
776  calc_p_isothermal(ordered_z_p[count-1], Pl);
777  if (Pl-Pu < metgrid_proximity) {
778  continue;
779  }
780  ordered_z_p[count] = orig_z_p[k];
781  ordered_data_p[count] = orig_data(i,j,k,src_comp);
782  count++;
783  }
784  kend_order = count;
785  }
786 
787  int ksta(0), kend(0);
788  if (metgrid_use_below_sfc && metgrid_use_sfc) {
789  // Use all levels.
790  ksta = 0;
791  kend = kend_order-1;
792  } else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
793  // Use all levels except for the surface.
794  int ksfc = 0;
795  for (int k=0; k < kmax_orig; k++) {
796  if (ordered_z_p[k] == orig_z_p[0]) {
797  ksfc = k;
798  break;
799  }
800  }
801  for (int k=ksfc; k < kmax_orig-1; k++) {
802  ordered_z_p[k] = ordered_z_p[k+1];
803  ordered_data_p[k] = ordered_data_p[k+1];
804  }
805  ksta = 0;
806  kend = kend_order-2;
807  } else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
808  // Use all levels above and including the surface.
809  int ksfc = 0;
810  for (int k=0; k < kend_order; k++) {
811  if (ordered_z_p[k] == orig_z_p[0]) {
812  ksfc = k;
813  break;
814  }
815  }
816  ksta = ksfc;
817  kend = kend_order-1;
818  } else {
819  // We shouldn't be in here!
820  amrex::Abort("metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
821  }
822 
823  // Insert the level of maximum winds.
824 // amrex::Real maxw_above_this_level = amrex::Real(30000.0);
825 // amrex::Real maxw_horiz_pres_diff = amrex::Real(5000.0);
826 // if ((flag_maxw == 1) && (use_maxw)) {
827 // amrex::Abort("metgrid initialization, use_maxw not yet implemented");
828 // }
829 
830  // Insert the level of the tropopause.
831 // amrex::Real trop_horiz_pres_diff = amrex::Real(5000.0);
832 // if ((flag_trop == 1) && (use_trop)) {
833 // amrex::Abort("metgrid initialization, use_trop not yet implemented");
834 // }
835 
836  amrex::GpuArray<amrex::Real,256> ordered_p;
837  amrex::Real* ordered_p_p = ordered_p.data();
838  for (int k=0; k < kend_order; k++) {
839  calc_p_isothermal(ordered_z_p[k], ordered_p_p[k]);
840  }
841 
842  // kfinal is the index of the topmost level retained in the final_*
843  // arrays, so those arrays hold kfinal+1 points.
844  int kfinal = 0;
845  amrex::GpuArray<amrex::Real,256> final_z;
846  amrex::GpuArray<amrex::Real,256> final_p;
847  amrex::GpuArray<amrex::Real,256> final_data;
848  amrex::Real* final_z_p = final_z.data();
849  amrex::Real* final_p_p = final_p.data();
850  amrex::Real* final_data_p = final_data.data();
851  final_z_p[0] = ordered_z[ksta];
852  final_p_p[0] = ordered_p[ksta];
853  final_data_p[0] = ordered_data[ksta];
854  for (int k=ksta+1; k <= kend; k++) {
855  // Skip levels that are "too close" to the prior retained level.
856  if ((final_p_p[kfinal]-ordered_p_p[k]) < metgrid_proximity) continue;
857  kfinal++;
858  final_z_p[kfinal] = ordered_z_p[k];
859  final_p_p[kfinal] = ordered_p_p[k];
860  final_data_p[kfinal] = ordered_data_p[k];
861  }
862 
863  // Call the interpolator.
864  lagrange_setup(var_type,
865  exp_interp,
866  kfinal+1,
867  kmax_new,
868  metgrid_order,
869  i,
870  j,
871  final_z_p,
872  final_p_p,
873  final_data_p,
874  new_z_p,
875  new_p_p,
876  new_data_p);
877 
878  // Optionally replace the lowest level of data with the surface
879  // field from the origin data. That is origin level 0: ordered_data[0] is
880  // the surface only when no origin level lies below it.
881  if (metgrid_retain_sfc) new_data[0] = orig_data(i,j,0,src_comp);
882 
883  // Save the interpolated data.
884  for (int k=0; k < kmax_new; k++) {
885  new_data_full(i,j,k,dst_comp) = new_data[k];
886  }
887 
888 }
889 
890 /**
891  * Linearly interpolate one metgrid column value onto an ERF point.
892  */
893 AMREX_FORCE_INLINE
894 AMREX_GPU_DEVICE
895 amrex::Real
896 interpolate_column_metgrid_linear (const int& i,
897  const int& j,
898  const int& k,
899  char stag,
900  int src_comp,
901  const amrex::Array4<amrex::Real const>& orig_z,
902  const amrex::Array4<amrex::Real const>& orig_data,
903  const amrex::Array4<amrex::Real const>& new_z)
904 {
905  // This subroutine is a bit ham-handed and can be cleaned up later.
906  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
907  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
908  // The number of origin levels, 0 .. ubound, so that the top one is searched too
909  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z + 1;
910 
911  amrex::Real z;
912  if (stag == 'X') {
913  z = fourth*(new_z(i,j,k)+new_z(i,j+1,k)+new_z(i,j,k+1)+new_z(i,j+1,k+1));
914  }
915  else if (stag == 'Y') {
916  z = fourth*(new_z(i,j,k)+new_z(i+1,j,k)+new_z(i,j,k+1)+new_z(i+1,j,k+1));
917  }
918  else if (stag == 'M') {
919  z = amrex::Real(0.125)*(new_z(i,j,k )+new_z(i,j+1,k )+new_z(i+1,j,k )+new_z(i+1,j+1,k )+
920  new_z(i,j,k+1)+new_z(i,j+1,k+1)+new_z(i+1,j,k+1)+new_z(i+1,j+1,k+1));
921  }
922 
923  // Initialized to NaN so that any path which fails to set them fails loudly
924  // (with amrex.fpe_trap_invalid=1) rather than returning plausible garbage.
925  amrex::Real z0 = std::numeric_limits<amrex::Real>::quiet_NaN();
926  amrex::Real z1 = std::numeric_limits<amrex::Real>::quiet_NaN();
927  int klow = -1;
928  int khi0 = -1;
929  amrex::Real dzlow = amrex::Real(1.0e12);
930  amrex::Real dzhi0 = -amrex::Real(1.0e12);
931  for (int kk = 0; kk < kmax_orig; kk++) {
932  amrex::Real orig_z_stag = zero;
933  if (stag == 'M') {
934  orig_z_stag = orig_z(i,j,kk);
935  }
936  if (stag == 'X') {
937  if (i == 0) {
938  orig_z_stag = orig_z(i,j,kk);
939  }
940  else if (i == imax_orig) {
941  orig_z_stag = orig_z(imax_orig-1,j,kk);
942  }
943  else {
944  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
945  }
946  }
947  else if (stag == 'Y') {
948  if (j == 0) {
949  orig_z_stag = orig_z(i,j,kk);
950  }
951  else if (j == jmax_orig) {
952  orig_z_stag = orig_z(i,jmax_orig-1,kk);
953  }
954  else {
955  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
956  }
957  }
958 
959  amrex::Real dz = z - orig_z_stag;
960  if ((dz < zero) && (dz > dzhi0)) {
961  dzhi0 = dz;
962  khi0 = kk;
963  z1 = orig_z_stag;
964  }
965  if ((dz >= zero) && (dz < dzlow)) {
966  dzlow = dz;
967  klow = kk;
968  z0 = orig_z_stag;
969  }
970  } // kk
971 
972  // extrapolate below the bottom surface
973  if (klow == -1) {
974  // klow was never found, so z0 was never set above. The two levels used
975  // here are khi0 and khi1, both above z, so z0 is the height at khi0 --
976  // which is exactly what the search loop left in z1 before we clobber it.
977  z0 = z1;
978  int khi1 = -1;
979  amrex::Real dzhi1 = -amrex::Real(1.0e12);
980  for (int kk = 0; kk < kmax_orig; kk++) {
981  amrex::Real orig_z_stag = zero;
982  if (stag == 'M') {
983  orig_z_stag = orig_z(i,j,kk);
984  }
985  else if (stag == 'X') {
986  if (i == 0) {
987  orig_z_stag = orig_z(i,j,kk);
988  }
989  else if (i == imax_orig) {
990  orig_z_stag = orig_z(imax_orig-1,j,kk);
991  }
992  else {
993  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
994  }
995  }
996  else if (stag == 'Y') {
997  if (j == 0) {
998  orig_z_stag = orig_z(i,j,kk);
999  }
1000  else if (j == jmax_orig) {
1001  orig_z_stag = orig_z(i,jmax_orig-1,kk);
1002  }
1003  else {
1004  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
1005  }
1006  }
1007  amrex::Real dz = z - orig_z_stag;
1008  if ((dz < zero) && (dz > dzhi1) && (kk != khi0)) {
1009  dzhi1 = dz;
1010  khi1 = kk;
1011  z1 = orig_z_stag;
1012  }
1013  }
1014  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
1015  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
1016  return ( y0-(y1-y0)/(z1-z0)*(z0-z) );
1017 
1018  // Extrapolate above the top surface
1019  } else if (khi0 == -1) {
1020  khi0 = klow - 1;
1021  int khi1 = klow;
1022  // khi0 was never found, so z1 was never set above. The upper level used
1023  // here is khi1 = klow, whose height is currently held in z0 -- capture
1024  // it before z0 is recomputed at the new (lower) khi0.
1025  z1 = z0;
1026  if (stag == 'M') {
1027  z0 = orig_z(i,j,khi0);
1028  }
1029  else if (stag == 'X') {
1030  if (i == 0) {
1031  z0 = orig_z(i,j,khi0);
1032  }
1033  else if (i == imax_orig) {
1034  z0 = orig_z(imax_orig-1,j,khi0);
1035  }
1036  else {
1037  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
1038  }
1039  }
1040  else if (stag == 'Y') {
1041  if (j == 0) {
1042  z0 = orig_z(i,j,khi0);
1043  }
1044  else if (j == jmax_orig) {
1045  z0 = orig_z(i,jmax_orig-1,khi0);
1046  }
1047  else {
1048  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
1049  }
1050  }
1051  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
1052  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
1053  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
1054  } else {
1055  // interpolate
1056  amrex::Real y0 = orig_data(i,j,klow,src_comp);
1057  amrex::Real y1 = orig_data(i,j,khi0,src_comp);
1058  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
1059 
1060  }
1061 }
1062 
1063 /**
1064  * Convert relative humidity to vapor mixing ratio using WRF-compatible bounds.
1065  */
1066 AMREX_FORCE_INLINE
1067 AMREX_GPU_DEVICE
1068 void
1069 rh_to_mxrat (int i,
1070  int j,
1071  int k,
1072  const amrex::Array4<amrex::Real const>& rhum,
1073  const amrex::Array4<amrex::Real const>& temp,
1074  const amrex::Array4<amrex::Real const>& pres,
1075  int src_indx,
1076  const amrex::Array4<amrex::Real>& mxrat)
1077 {
1078  amrex::Real qv_max_p_safe = amrex::Real(10000.0); // WRF default value
1079  amrex::Real qv_max_flag = amrex::Real(1.0e-5); // WRF default value
1080  amrex::Real qv_max_value = amrex::Real(3.0e-6); // WRF default value
1081  amrex::Real qv_min_p_safe = amrex::Real(110000.0); // WRF default value
1082  amrex::Real qv_min_flag = amrex::Real(1.0e-6); // WRF default value
1083  amrex::Real qv_min_value = amrex::Real(1.0e-6); // WRF default value
1084  amrex::Real eps = RdoRv;
1085  amrex::Real svp1 = amrex::Real(0.6112);
1086  amrex::Real svp2 = amrex::Real(17.67);
1087  amrex::Real svp3 = amrex::Real(29.65);
1088  amrex::Real svpt0 = amrex::Real(273.15);
1089  // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior)
1090  if (temp(i,j,k) != zero) {
1091  amrex::Real es=amrex::Real(0.01)*rhum(i,j,k)*svp1*amrex::Real(10.0)*std::exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3));
1092  if (es >= pres(i,j,k)/amrex::Real(100.0)) {
1093  // vapor pressure exceeds total pressure
1094  mxrat(i,j,k,src_indx) = amrex::Math::powi<-6>(amrex::Real(10.0));
1095  }
1096  else {
1097  mxrat(i,j,k,src_indx) = amrex::max(eps*es/(pres(i,j,k)/amrex::Real(100.0)-es), amrex::Real(1.0e-6));
1098  }
1099  }
1100  else {
1101  // I don't know why there's a fringe case handled in WRF where T is absolute zero...
1102  // Let's just deal with it here in case we also end up needing it.
1103  mxrat(i,j,k,src_indx) = amrex::Real(1.0e-6);
1104  }
1105  // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1.
1106  // For pressures above a defined level, reasonable Qv values should be
1107  // a certain value or smaller. If they are larger than this, the input data
1108  // probably had "missing" RH, and we filled in some values. This is an
1109  // attempt to catch those. Also, set the minimum value for the entire
1110  // domain that is above the selected pressure level.
1111  if (pres(i,j,k) < qv_max_p_safe) {
1112  if (mxrat(i,j,k,src_indx) > qv_max_flag) {
1113  mxrat(i,j,k,src_indx) = qv_max_value;
1114  }
1115  }
1116  if (pres(i,j,k) < qv_min_p_safe) {
1117  if (mxrat(i,j,k,src_indx) < qv_min_flag) {
1118  mxrat(i,j,k,src_indx) = qv_min_value;
1119  }
1120  }
1121 }
1122 #endif
ERF_BaseStateParams.H
ERF_Constants.H
p_0
constexpr amrex::Real p_0
Definition: ERF_Constants.H:69
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:72
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:49
RdoCp
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:57
RdoRv
constexpr amrex::Real RdoRv
Definition: ERF_Constants.H:60
Coord::y
@ y
Coord::x
@ x
Coord::z
@ z
ERF_EOS.H
ERF_HSEUtils.H
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:12
use_moisture
const bool use_moisture
Definition: ERF_InitCustomPert_ABL.H:71
rdOcp
const Real rdOcp
Definition: ERF_InitCustomPert_ABL.H:72
init_terrain_from_metgrid
void init_terrain_from_metgrid(const bool &avg_grid_faces_to_nodes, amrex::Geometry &geom, amrex::MultiFab *z_phys_nd, amrex::FArrayBox &NC_hgt_fab)
lagrange_interp
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_interp(const int &order, amrex::Real *x, amrex::Real *y, amrex::Real &new_x, amrex::Real &new_y)
Definition: ERF_MetgridUtils.H:231
interpolate_column_metgrid_linear
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real interpolate_column_metgrid_linear(const int &i, const int &j, const int &k, char stag, int src_comp, const amrex::Array4< amrex::Real const > &orig_z, const amrex::Array4< amrex::Real const > &orig_data, const amrex::Array4< amrex::Real const > &new_z)
Definition: ERF_MetgridUtils.H:896
interpolate_column_metgrid
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void interpolate_column_metgrid(const bool &metgrid_use_below_sfc, const bool &metgrid_use_sfc, const bool &exp_interp, const bool &metgrid_retain_sfc, const amrex::Real &metgrid_proximity, const int &metgrid_order, const int &metgrid_force_sfc_k, const int &i, const int &j, const int &kmax, const int &src_comp, const int &dst_comp, char var_type, char stag, const amrex::Array4< amrex::Real const > &orig_z_full, const amrex::Array4< amrex::Real const > &orig_data, const amrex::Array4< amrex::Real const > &new_z_full, const amrex::Array4< amrex::Real > &new_data_full)
Definition: ERF_MetgridUtils.H:587
read_from_metgrid
void read_from_metgrid(int lev, int itime, const amrex::Box &domain, const std::string &fname, std::string &NC_dateTime, double &NC_epochTime, int &flag_psfc, int &flag_msf, int &flag_sst, int &flag_tsk, int &flag_lmask, int &NC_nx, int &NC_ny, amrex::Real &NC_dx, amrex::Real &NC_dy, amrex::FArrayBox &NC_xvel_fab, amrex::FArrayBox &NC_yvel_fab, amrex::FArrayBox &NC_temp_fab, amrex::FArrayBox &NC_rhum_fab, amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &NC_ght_fab, amrex::FArrayBox &NC_hgt_fab, amrex::FArrayBox &NC_psfc_fab, amrex::FArrayBox &NC_msfu_fab, amrex::FArrayBox &NC_msfv_fab, amrex::FArrayBox &NC_msfm_fab, amrex::FArrayBox &NC_sst_fab, amrex::FArrayBox &NC_tsk_fab, amrex::FArrayBox &NC_LAT_fab, amrex::FArrayBox &NC_LON_fab, amrex::IArrayBox &NC_lmask_iab, amrex::Geometry &geom)
metgrid_surface_pressure
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool metgrid_surface_pressure(const bool debug_psfc, const int flag_psfc, const amrex::Real file_psfc, const amrex::Real z_sfc, const amrex::Real P00, const amrex::Real T00, const amrex::Real TLP, amrex::Real &pressure)
Definition: ERF_MetgridUtils.H:30
metgrid_surface_theta
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool metgrid_surface_theta(const amrex::Real temperature, const bool debug_psfc, const int flag_psfc, const amrex::Real file_psfc, const amrex::Real z_sfc, const amrex::Real P00, const amrex::Real T00, const amrex::Real TLP, const amrex::Real rdOcp, amrex::Real &theta)
Definition: ERF_MetgridUtils.H:79
calc_p_isothermal
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:575
init_msfs_from_metgrid
void init_msfs_from_metgrid(const bool metgrid_debug_msf, amrex::FArrayBox &msfu_fab, amrex::FArrayBox &msfv_fab, amrex::FArrayBox &msfm_fab, const int &flag_msf, amrex::FArrayBox &NC_MSFU_fab, amrex::FArrayBox &NC_MSFV_fab, amrex::FArrayBox &NC_MSFM_fab)
lagrange_setup
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_setup(char var_type, const bool &exp_interp, const int &orig_npts, const int &new_npts, const int &order, const int &i, const int &j, amrex::Real *orig_x_z, amrex::Real *orig_x_p, amrex::Real *orig_y, amrex::Real *new_x_z, amrex::Real *new_x_p, amrex::Real *new_y)
Definition: ERF_MetgridUtils.H:268
init_base_state_from_metgrid
void init_base_state_from_metgrid(const bool use_moisture, const bool metgrid_debug_psfc, const amrex::Real l_rdOcp, const amrex::Box &domain, const amrex::Box &valid_bx, const int &flag_psfc, amrex::FArrayBox &state_fab, amrex::FArrayBox &r_hse_fab, amrex::FArrayBox &p_hse_fab, amrex::FArrayBox &pi_hse_fab, amrex::FArrayBox &th_hse_fab, amrex::FArrayBox &qv_hse_fab, amrex::FArrayBox &z_phys_nd_fab, amrex::FArrayBox &z_phys_cc_fab, const amrex::FArrayBox &NC_psfc_fab, int *base_state_column_failed, const BaseStateParams &bsp)
init_state_from_metgrid
void init_state_from_metgrid(const int lev, const int itime, const bool use_moisture, const bool interp_theta, const bool metgrid_debug_quiescent, const bool metgrid_debug_isothermal, const bool metgrid_debug_dry, const bool metgrid_basic_linear, const bool metgrid_use_below_sfc, const bool metgrid_use_sfc, const bool metgrid_retain_sfc, const amrex::Real metgrid_proximity, const int metgrid_order, const int metgrid_metgrid_force_sfc_k, const amrex::Real l_rdOcp, amrex::Box &tbxc, amrex::Box &tbxu, amrex::Box &tbxv, amrex::Box &tbxw, amrex::FArrayBox &state_fab, amrex::FArrayBox &x_vel_fab, amrex::FArrayBox &y_vel_fab, amrex::FArrayBox &z_vel_fab, amrex::FArrayBox &z_phys_nd_fab, const amrex::FArrayBox &NC_ght_fab, const amrex::FArrayBox &NC_xvel_fab, const amrex::FArrayBox &NC_yvel_fab, const amrex::FArrayBox &NC_temp_fab, const amrex::FArrayBox &NC_rhum_fab, const amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &tmp_src_fab, amrex::FArrayBox &tmp_dst_fab, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xlo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xhi, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_ylo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_yhi, const amrex::Array4< const int > &mask_c_arr, const amrex::Array4< const int > &mask_u_arr, const amrex::Array4< const int > &mask_v_arr)
rh_to_mxrat
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void rh_to_mxrat(int i, int j, int k, const amrex::Array4< amrex::Real const > &rhum, const amrex::Array4< amrex::Real const > &temp, const amrex::Array4< amrex::Real const > &pres, int src_indx, const amrex::Array4< amrex::Real > &mxrat)
Definition: ERF_MetgridUtils.H:1069
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
ERF_ProbCommon.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
ERF_SurfaceTemperature.H
ERF_Utils.H
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:28
MicVar_Kess::theta
@ theta
Definition: ERF_Kessler.H:26
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:273
WSM6Ind::p
@ p
Definition: ERF_WSM6.H:301
erf_surface_temperature::temperature_to_theta
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool temperature_to_theta(const amrex::Real temperature, const amrex::Real pressure, const amrex::Real rdOcp, amrex::Real &theta)
Definition: ERF_SurfaceTemperature.H:79
module_model_constants::svp1
real(c_double), parameter svp1
Definition: ERF_module_model_constants.F90:78
module_model_constants::svp3
real(c_double), parameter svp3
Definition: ERF_module_model_constants.F90:80
module_model_constants::svp2
real(c_double), parameter svp2
Definition: ERF_module_model_constants.F90:79
module_model_constants::svpt0
real(c_double), parameter svpt0
Definition: ERF_module_model_constants.F90:81
BaseStateParams
Definition: ERF_BaseStateParams.H:35