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  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
615  int kmax_new = kmax + 1;
616 
617  AMREX_ASSERT(kmax_orig < 256);
618  AMREX_ASSERT(kmax_new < 256);
619 
620  amrex::GpuArray<amrex::Real,256> new_z;
621  amrex::GpuArray<amrex::Real,256> new_p;
622  amrex::GpuArray<amrex::Real,256> new_data;
623  amrex::Real* new_z_p = new_z.data();
624  amrex::Real* new_p_p = new_p.data();
625  amrex::Real* new_data_p = new_data.data();
626  for (int k=0; k < kmax_new; k++) {
627  if (stag == 'X') {
628  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));
629  } else if (stag == 'Y') {
630  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));
631  } else if (stag == 'M') {
632  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 )+
633  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));
634  }
635  calc_p_isothermal(new_z_p[k], new_p_p[k]);
636  }
637 
638  amrex::GpuArray<amrex::Real,256> orig_z;
639  amrex::Real* orig_z_p = orig_z.data();
640  for (int k=0; k < kmax_orig; k++) {
641  if (stag == 'M') {
642  orig_z_p[k] = orig_z_full(i,j,k);
643  } else if (stag == 'X') {
644  if (i <= imin_orig) {
645  orig_z_p[k] = orig_z_full(i,j,k);
646  } else if (i >= imax_orig) {
647  orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
648  } else {
649  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
650  }
651  } else if (stag == 'Y') {
652  if (j <= jmin_orig) {
653  orig_z_p[k] = orig_z_full(i,j,k);
654  } else if (j >= jmax_orig) {
655  orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
656  } else {
657  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
658  }
659  }
660  }
661 
662  // Check if the data is top-down instead of bottom-up.
663  bool flip_data_required = false;
664  if (orig_z[1] > orig_z[kmax_orig-1]) flip_data_required = true;
665  if (flip_data_required) amrex::Abort("metgrid initialization flip_data_required. Not yet implemented.");
666 
667  // Search for the first level above the surface in the origin data.
668  // This is needed since the origin model topography will be
669  // different than the topography processed by WPS.
670  int k_above_sfc = 0;
671  for (int k=1; k < kmax_orig; k++) {
672  if (orig_z_p[k] > orig_z_p[0]) {
673  k_above_sfc = k;
674  break;
675  }
676  }
677 
678  int kend_order;
679  amrex::GpuArray<amrex::Real,256> ordered_z;
680  amrex::GpuArray<amrex::Real,256> ordered_data;
681  amrex::Real* ordered_z_p = ordered_z.data();
682  amrex::Real* ordered_data_p = ordered_data.data();
683  if (k_above_sfc > 1) {
684  // The levels are not monotonically increasing in height, so
685  // we sort and then make "artistic" quality control choices.
686  int count = 0;
687 
688  // Insert levels that are below the surface.
689  for (int k=1; k < k_above_sfc; k++) {
690  ordered_z_p[count] = orig_z_p[k];
691  ordered_data_p[count] = orig_data(i,j,k,src_comp);
692  count++;
693  }
694 
695  // Check if the level that is nearest to and below the surface
696  // is "too close". If so, we'll ignore the upper level and keep
697  // the lower. Origin data is likely to be on pressure levels
698  // with higher spatial resolution near-surface, which supports
699  // the choice of eliminating levels that are "too close" in
700  // pressure-space. For simplicity, calculate delta P assuming a
701  // baroclinic atmosphere.
702  amrex::Real Pu, Pl;
703  calc_p_isothermal(orig_z_p[0], Pu);
704  calc_p_isothermal(ordered_z_p[count-1], Pl);
705  if (Pl-Pu < metgrid_proximity) {
706  count--;
707  }
708 
709  // Insert the surface level.
710  ordered_z_p[count] = orig_z_p[0];
711  ordered_data_p[count] = orig_data(i,j,0,src_comp);
712  count++;
713 
714  // Quoting WRF's comments, the next level to use is at,
715  // "... ta da, the first level above the surface. I know, wow."
716  int knext = k_above_sfc;
717  // Conditionally more strongly use the surface data by removing
718  // levels between the surface and the height corresponding to a
719  // set number of ERF levels from the surface. This forces the
720  // interpolator to use the surface data up through a number of
721  // ERF levels from the surface.
722  if (metgrid_force_sfc_k > 0) {
723  for (int k=k_above_sfc; k < kmax_orig; k++) {
724  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
725  knext = k;
726  break;
727  }
728  }
729  }
730 
731  // Check if the level that is nearest to and above the surface
732  // is "too close". If so, we'll ignore that level.
733  calc_p_isothermal(orig_z_p[knext], Pu);
734  calc_p_isothermal(ordered_z_p[count-1], Pl);
735  if (Pl-Pu < metgrid_proximity) {
736  knext++;
737  }
738 
739  // Insert levels that are above the surface.
740  for (int k=knext; k < kmax_orig; k++) {
741  ordered_z_p[count] = orig_z_p[k];
742  ordered_data_p[count] = orig_data(i,j,k,src_comp);
743  count++;
744  }
745 
746  kend_order = count;
747  } else {
748  // The surface is the lowest level in the origin data.
749 
750  // Insert the surface.
751  ordered_z_p[0] = orig_z[0];
752  ordered_data_p[0] = orig_data(i,j,0,src_comp);
753 
754  // Similar to above, conditionally more strongly use the
755  // surface data.
756  int count = 1;
757  int knext = count;
758  if (metgrid_force_sfc_k > 0) {
759  for (int k=knext; k < kmax_orig; k++) {
760  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k]) {
761  knext = k;
762  break;
763  }
764  }
765  }
766 
767  // Insert the remaining levels, again ignoring levels that are
768  // "too close" to the prior valid level.
769  for (int k=knext; k < kmax_orig; k++) {
770  amrex::Real Pu, Pl;
771  calc_p_isothermal(orig_z_p[k], Pu);
772  calc_p_isothermal(ordered_z_p[count-1], Pl);
773  if (Pl-Pu < metgrid_proximity) {
774  continue;
775  }
776  ordered_z_p[count] = orig_z_p[k];
777  ordered_data_p[count] = orig_data(i,j,k,src_comp);
778  count++;
779  }
780  kend_order = count;
781  }
782 
783  int ksta(0), kend(0);
784  if (metgrid_use_below_sfc && metgrid_use_sfc) {
785  // Use all levels.
786  ksta = 0;
787  kend = kend_order-1;
788  } else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
789  // Use all levels except for the surface.
790  int ksfc = 0;
791  for (int k=0; k < kmax_orig; k++) {
792  if (ordered_z_p[k] == orig_z_p[0]) {
793  ksfc = k;
794  break;
795  }
796  }
797  for (int k=ksfc; k < kmax_orig-1; k++) {
798  ordered_z_p[k] = ordered_z_p[k+1];
799  ordered_data_p[k] = ordered_data_p[k+1];
800  }
801  ksta = 0;
802  kend = kend_order-2;
803  } else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
804  // Use all levels above and including the surface.
805  int ksfc = 0;
806  for (int k=0; k < kend_order; k++) {
807  if (ordered_z_p[k] == orig_z_p[0]) {
808  ksfc = k;
809  break;
810  }
811  }
812  ksta = ksfc;
813  kend = kend_order-1;
814  } else {
815  // We shouldn't be in here!
816  amrex::Abort("metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
817  }
818 
819  // Insert the level of maximum winds.
820 // amrex::Real maxw_above_this_level = amrex::Real(30000.0);
821 // amrex::Real maxw_horiz_pres_diff = amrex::Real(5000.0);
822 // if ((flag_maxw == 1) && (use_maxw)) {
823 // amrex::Abort("metgrid initialization, use_maxw not yet implemented");
824 // }
825 
826  // Insert the level of the tropopause.
827 // amrex::Real trop_horiz_pres_diff = amrex::Real(5000.0);
828 // if ((flag_trop == 1) && (use_trop)) {
829 // amrex::Abort("metgrid initialization, use_trop not yet implemented");
830 // }
831 
832  amrex::GpuArray<amrex::Real,256> ordered_p;
833  amrex::Real* ordered_p_p = ordered_p.data();
834  for (int k=0; k < kend_order; k++) {
835  calc_p_isothermal(ordered_z_p[k], ordered_p_p[k]);
836  }
837 
838  // kfinal is the index of the topmost level retained in the final_*
839  // arrays, so those arrays hold kfinal+1 points.
840  int kfinal = 0;
841  amrex::GpuArray<amrex::Real,256> final_z;
842  amrex::GpuArray<amrex::Real,256> final_p;
843  amrex::GpuArray<amrex::Real,256> final_data;
844  amrex::Real* final_z_p = final_z.data();
845  amrex::Real* final_p_p = final_p.data();
846  amrex::Real* final_data_p = final_data.data();
847  final_z_p[0] = ordered_z[ksta];
848  final_p_p[0] = ordered_p[ksta];
849  final_data_p[0] = ordered_data[ksta];
850  for (int k=ksta+1; k <= kend; k++) {
851  // Skip levels that are "too close" to the prior retained level.
852  if ((final_p_p[kfinal]-ordered_p_p[k]) < metgrid_proximity) continue;
853  kfinal++;
854  final_z_p[kfinal] = ordered_z_p[k];
855  final_p_p[kfinal] = ordered_p_p[k];
856  final_data_p[kfinal] = ordered_data_p[k];
857  }
858 
859  // Call the interpolator.
860  lagrange_setup(var_type,
861  exp_interp,
862  kfinal+1,
863  kmax_new,
864  metgrid_order,
865  i,
866  j,
867  final_z_p,
868  final_p_p,
869  final_data_p,
870  new_z_p,
871  new_p_p,
872  new_data_p);
873 
874  // Optionally replace the lowest level of data with the surface
875  // field from the origin data.
876  if (metgrid_retain_sfc) new_data[0] = ordered_data[0];
877 
878  // Save the interpolated data.
879  for (int k=0; k < kmax_new; k++) {
880  new_data_full(i,j,k,dst_comp) = new_data[k];
881  }
882 
883 }
884 
885 /**
886  * Linearly interpolate one metgrid column value onto an ERF point.
887  */
888 AMREX_FORCE_INLINE
889 AMREX_GPU_DEVICE
890 amrex::Real
891 interpolate_column_metgrid_linear (const int& i,
892  const int& j,
893  const int& k,
894  char stag,
895  int src_comp,
896  const amrex::Array4<amrex::Real const>& orig_z,
897  const amrex::Array4<amrex::Real const>& orig_data,
898  const amrex::Array4<amrex::Real const>& new_z)
899 {
900  // This subroutine is a bit ham-handed and can be cleaned up later.
901  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
902  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
903  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
904 
905  amrex::Real z;
906  if (stag == 'X') {
907  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));
908  }
909  else if (stag == 'Y') {
910  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));
911  }
912  else if (stag == 'M') {
913  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 )+
914  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));
915  }
916 
917  // Initialized to NaN so that any path which fails to set them fails loudly
918  // (with amrex.fpe_trap_invalid=1) rather than returning plausible garbage.
919  amrex::Real z0 = std::numeric_limits<amrex::Real>::quiet_NaN();
920  amrex::Real z1 = std::numeric_limits<amrex::Real>::quiet_NaN();
921  int klow = -1;
922  int khi0 = -1;
923  amrex::Real dzlow = amrex::Real(1.0e12);
924  amrex::Real dzhi0 = -amrex::Real(1.0e12);
925  for (int kk = 0; kk < kmax_orig; kk++) {
926  amrex::Real orig_z_stag = zero;
927  if (stag == 'M') {
928  orig_z_stag = orig_z(i,j,kk);
929  }
930  if (stag == 'X') {
931  if (i == 0) {
932  orig_z_stag = orig_z(i,j,kk);
933  }
934  else if (i == imax_orig) {
935  orig_z_stag = orig_z(imax_orig-1,j,kk);
936  }
937  else {
938  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
939  }
940  }
941  else if (stag == 'Y') {
942  if (j == 0) {
943  orig_z_stag = orig_z(i,j,kk);
944  }
945  else if (j == jmax_orig) {
946  orig_z_stag = orig_z(i,jmax_orig-1,kk);
947  }
948  else {
949  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
950  }
951  }
952 
953  amrex::Real dz = z - orig_z_stag;
954  if ((dz < zero) && (dz > dzhi0)) {
955  dzhi0 = dz;
956  khi0 = kk;
957  z1 = orig_z_stag;
958  }
959  if ((dz >= zero) && (dz < dzlow)) {
960  dzlow = dz;
961  klow = kk;
962  z0 = orig_z_stag;
963  }
964  } // kk
965 
966  // extrapolate below the bottom surface
967  if (klow == -1) {
968  // klow was never found, so z0 was never set above. The two levels used
969  // here are khi0 and khi1, both above z, so z0 is the height at khi0 --
970  // which is exactly what the search loop left in z1 before we clobber it.
971  z0 = z1;
972  int khi1 = -1;
973  amrex::Real dzhi1 = -amrex::Real(1.0e12);
974  for (int kk = 0; kk < kmax_orig; kk++) {
975  amrex::Real orig_z_stag = zero;
976  if (stag == 'M') {
977  orig_z_stag = orig_z(i,j,kk);
978  }
979  else if (stag == 'X') {
980  if (i == 0) {
981  orig_z_stag = orig_z(i,j,kk);
982  }
983  else if (i == imax_orig) {
984  orig_z_stag = orig_z(imax_orig-1,j,kk);
985  }
986  else {
987  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
988  }
989  }
990  else if (stag == 'Y') {
991  if (j == 0) {
992  orig_z_stag = orig_z(i,j,kk);
993  }
994  else if (j == jmax_orig) {
995  orig_z_stag = orig_z(i,jmax_orig-1,kk);
996  }
997  else {
998  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
999  }
1000  }
1001  amrex::Real dz = z - orig_z_stag;
1002  if ((dz < zero) && (dz > dzhi1) && (kk != khi0)) {
1003  dzhi1 = dz;
1004  khi1 = kk;
1005  z1 = orig_z_stag;
1006  }
1007  }
1008  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
1009  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
1010  return ( y0-(y1-y0)/(z1-z0)*(z0-z) );
1011 
1012  // Extrapolate above the top surface
1013  } else if (khi0 == -1) {
1014  khi0 = klow - 1;
1015  int khi1 = klow;
1016  // khi0 was never found, so z1 was never set above. The upper level used
1017  // here is khi1 = klow, whose height is currently held in z0 -- capture
1018  // it before z0 is recomputed at the new (lower) khi0.
1019  z1 = z0;
1020  if (stag == 'M') {
1021  z0 = orig_z(i,j,khi0);
1022  }
1023  else if (stag == 'X') {
1024  if (i == 0) {
1025  z0 = orig_z(i,j,khi0);
1026  }
1027  else if (i == imax_orig) {
1028  z0 = orig_z(imax_orig-1,j,khi0);
1029  }
1030  else {
1031  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
1032  }
1033  }
1034  else if (stag == 'Y') {
1035  if (j == 0) {
1036  z0 = orig_z(i,j,khi0);
1037  }
1038  else if (j == jmax_orig) {
1039  z0 = orig_z(i,jmax_orig-1,khi0);
1040  }
1041  else {
1042  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
1043  }
1044  }
1045  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
1046  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
1047  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
1048  } else {
1049  // interpolate
1050  amrex::Real y0 = orig_data(i,j,klow,src_comp);
1051  amrex::Real y1 = orig_data(i,j,khi0,src_comp);
1052  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
1053 
1054  }
1055 }
1056 
1057 /**
1058  * Convert relative humidity to vapor mixing ratio using WRF-compatible bounds.
1059  */
1060 AMREX_FORCE_INLINE
1061 AMREX_GPU_DEVICE
1062 void
1063 rh_to_mxrat (int i,
1064  int j,
1065  int k,
1066  const amrex::Array4<amrex::Real const>& rhum,
1067  const amrex::Array4<amrex::Real const>& temp,
1068  const amrex::Array4<amrex::Real const>& pres,
1069  int src_indx,
1070  const amrex::Array4<amrex::Real>& mxrat)
1071 {
1072  amrex::Real qv_max_p_safe = amrex::Real(10000.0); // WRF default value
1073  amrex::Real qv_max_flag = amrex::Real(1.0e-5); // WRF default value
1074  amrex::Real qv_max_value = amrex::Real(3.0e-6); // WRF default value
1075  amrex::Real qv_min_p_safe = amrex::Real(110000.0); // WRF default value
1076  amrex::Real qv_min_flag = amrex::Real(1.0e-6); // WRF default value
1077  amrex::Real qv_min_value = amrex::Real(1.0e-6); // WRF default value
1078  amrex::Real eps = RdoRv;
1079  amrex::Real svp1 = amrex::Real(0.6112);
1080  amrex::Real svp2 = amrex::Real(17.67);
1081  amrex::Real svp3 = amrex::Real(29.65);
1082  amrex::Real svpt0 = amrex::Real(273.15);
1083  // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior)
1084  if (temp(i,j,k) != zero) {
1085  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));
1086  if (es >= pres(i,j,k)/amrex::Real(100.0)) {
1087  // vapor pressure exceeds total pressure
1088  mxrat(i,j,k,src_indx) = amrex::Math::powi<-6>(amrex::Real(10.0));
1089  }
1090  else {
1091  mxrat(i,j,k,src_indx) = amrex::max(eps*es/(pres(i,j,k)/amrex::Real(100.0)-es), amrex::Real(1.0e-6));
1092  }
1093  }
1094  else {
1095  // I don't know why there's a fringe case handled in WRF where T is absolute zero...
1096  // Let's just deal with it here in case we also end up needing it.
1097  mxrat(i,j,k,src_indx) = amrex::Real(1.0e-6);
1098  }
1099  // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1.
1100  // For pressures above a defined level, reasonable Qv values should be
1101  // a certain value or smaller. If they are larger than this, the input data
1102  // probably had "missing" RH, and we filled in some values. This is an
1103  // attempt to catch those. Also, set the minimum value for the entire
1104  // domain that is above the selected pressure level.
1105  if (pres(i,j,k) < qv_max_p_safe) {
1106  if (mxrat(i,j,k,src_indx) > qv_max_flag) {
1107  mxrat(i,j,k,src_indx) = qv_max_value;
1108  }
1109  }
1110  if (pres(i,j,k) < qv_min_p_safe) {
1111  if (mxrat(i,j,k,src_indx) < qv_min_flag) {
1112  mxrat(i,j,k,src_indx) = qv_min_value;
1113  }
1114  }
1115 }
1116 #endif
ERF_BaseStateParams.H
ERF_Constants.H
p_0
constexpr amrex::Real p_0
Definition: ERF_Constants.H:53
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:56
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:34
RdoCp
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:41
RdoRv
constexpr amrex::Real RdoRv
Definition: ERF_Constants.H:44
Coord::y
@ y
Coord::x
@ x
Coord::z
@ z
ERF_EOS.H
ERF_HSEUtils.H
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
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:891
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:1063
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:35
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
ERF_ProbCommon.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
ERF_SurfaceTemperature.H
ERF_Utils.H
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:19
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:28
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
WSM6Ind::p
@ p
Definition: ERF_WSM6.H:280
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:73
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