1 #ifndef ERF_METGRIDUTIL_H_
2 #define ERF_METGRIDUTIL_H_
12 const amrex::Box& domain,
13 const std::string& fname,
14 std::string& NC_dateTime,
15 amrex::Real& NC_epochTime,
24 amrex::FArrayBox& NC_xvel_fab,
25 amrex::FArrayBox& NC_yvel_fab,
26 amrex::FArrayBox& NC_temp_fab,
27 amrex::FArrayBox& NC_rhum_fab,
28 amrex::FArrayBox& NC_pres_fab,
29 amrex::FArrayBox& NC_ght_fab,
30 amrex::FArrayBox& NC_hgt_fab,
31 amrex::FArrayBox& NC_psfc_fab,
32 amrex::FArrayBox& NC_msfu_fab,
33 amrex::FArrayBox& NC_msfv_fab,
34 amrex::FArrayBox& NC_msfm_fab,
35 amrex::FArrayBox& NC_sst_fab,
36 amrex::FArrayBox& NC_LAT_fab,
37 amrex::FArrayBox& NC_LON_fab,
38 amrex::IArrayBox& NC_lmask_iab,
39 amrex::Geometry& geom);
43 const amrex::Vector<amrex::FArrayBox>& NC_hgt_fab);
47 const bool interp_theta,
48 const bool metgrid_debug_quiescent,
49 const bool metgrid_debug_isothermal,
50 const bool metgrid_debug_dry,
51 const bool metgrid_basic_linear,
52 const bool metgrid_use_below_sfc,
53 const bool metgrid_use_sfc,
54 const bool metgrid_retain_sfc,
55 const amrex::Real metgrid_proximity,
56 const int metgrid_order,
57 const int metgrid_metgrid_force_sfc_k,
58 const amrex::Real l_rdOcp,
62 amrex::FArrayBox& state_fab,
63 amrex::FArrayBox& x_vel_fab,
64 amrex::FArrayBox& y_vel_fab,
65 amrex::FArrayBox& z_vel_fab,
66 amrex::FArrayBox& z_phys_nd_fab,
67 const amrex::Vector<amrex::FArrayBox>& NC_hgt_fab,
68 const amrex::Vector<amrex::FArrayBox>& NC_ght_fab,
69 const amrex::Vector<amrex::FArrayBox>& NC_xvel_fab,
70 const amrex::Vector<amrex::FArrayBox>& NC_yvel_fab,
71 const amrex::Vector<amrex::FArrayBox>& NC_temp_fab,
72 const amrex::Vector<amrex::FArrayBox>& NC_rhum_fab,
73 const amrex::Vector<amrex::FArrayBox>& NC_pres_fab,
74 amrex::FArrayBox& p_interp_fab,
75 amrex::FArrayBox& t_interp_fab,
76 amrex::FArrayBox& theta_fab,
77 amrex::FArrayBox& mxrat_fab,
78 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xlo,
79 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xhi,
80 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_ylo,
81 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_yhi,
82 const amrex::Array4<const int>& mask_c_arr,
83 const amrex::Array4<const int>& mask_u_arr,
84 const amrex::Array4<const int>& mask_v_arr);
88 amrex::FArrayBox& msfu_fab,
89 amrex::FArrayBox& msfv_fab,
90 amrex::FArrayBox& msfm_fab,
92 const amrex::Vector<amrex::FArrayBox>& NC_MSFU_fab,
93 const amrex::Vector<amrex::FArrayBox>& NC_MSFV_fab,
94 const amrex::Vector<amrex::FArrayBox>& NC_MSFM_fab);
98 const bool metgrid_debug_psfc,
99 const amrex::Real l_rdOcp,
100 const amrex::Box& valid_bx,
101 const amrex::Vector<int>& flag_psfc,
102 amrex::FArrayBox& state,
103 amrex::FArrayBox& r_hse_fab,
104 amrex::FArrayBox& p_hse_fab,
105 amrex::FArrayBox& pi_hse_fab,
106 amrex::FArrayBox& th_hse_fab,
107 amrex::FArrayBox& qv_hse_fab,
108 amrex::FArrayBox& z_phys_cc_fab,
109 const amrex::Vector<amrex::FArrayBox>& NC_psfc_fab);
125 amrex::Real Px = 0.0;
126 for (
int i=0; i <= order; i++) {
129 for (
int k=0; k <= order; k++) {
130 if (k == i)
continue;
145 const bool& exp_interp,
151 amrex::Real* orig_x_z,
152 amrex::Real* orig_x_p,
154 amrex::Real* new_x_z,
155 amrex::Real* new_x_p,
159 amrex::Real CRC_const1 = 11880.516;
160 amrex::Real CRC_const2 = 0.1902632;
161 amrex::Real CRC_const3 = 0.0065;
163 amrex::ignore_unused(i,j);
164 #ifndef AMREX_USE_GPU
168 for (
int new_k=0; new_k < new_n; new_k++) {
169 #ifndef AMREX_USE_GPU
170 if (debug) amrex::Print() <<
"new_k=" << new_k;
174 bool extrapolating =
true;
176 for (
int ko=0; ko < orig_n-1; ko++) {
177 amrex::Real a = new_x_z[new_k]-orig_x_z[ko];
178 amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1];
182 extrapolating =
false;
188 if (var_type ==
'T') {
191 amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[0];
192 amrex::Real avg_of_extrap_p = 0.5*(new_x_p[new_k]+orig_x_p[0]);
193 amrex::Real temp_extrap_starting_point = orig_y[0]*std::pow(orig_x_p[0]/100000.0,
R_d/
Cp_d);
194 amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/100.0, CRC_const2-1.0);
195 amrex::Real dZ = dZdP*(depth_of_extrap_in_p/100.0);
196 amrex::Real dT = dZ*CRC_const3;
197 new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(100000.0/new_x_p[new_k],
R_d/
Cp_d);
201 new_y[new_k] = orig_y[0];
207 if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_n-1)) {
209 int ksta = kl-(((order+1)/2)-1);
210 int kend = ksta+order;
211 #ifndef AMREX_USE_GPU
212 int ksize = kend-ksta;
213 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;
216 amrex::GpuArray<amrex::Real,9> orig_x_sub;
217 amrex::Real* orig_x_sub_p = orig_x_sub.data();
219 new_x = new_x_p[new_k];
220 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
222 new_x = new_x_z[new_k];
223 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
225 amrex::GpuArray<amrex::Real,9> orig_y_sub;
226 amrex::Real* orig_y_sub_p = orig_y_sub.data();
227 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
228 #ifndef AMREX_USE_GPU
230 amrex::Print() <<
" orig_x_sub_p = [";
231 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
232 amrex::Print() <<
"]" << std::endl;
233 amrex::Print() <<
" orig_y_sub_p = [";
234 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
235 amrex::Print() <<
"]" << std::endl;
238 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
244 #ifndef AMREX_USE_GPU
245 int ksize = kend-ksta+1;
246 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;
249 amrex::GpuArray<amrex::Real,2> orig_x_sub;
250 amrex::GpuArray<amrex::Real,2> orig_y_sub;
251 amrex::Real* orig_x_sub_p = orig_x_sub.data();
252 amrex::Real* orig_y_sub_p = orig_y_sub.data();
254 new_x = new_x_p[new_k];
255 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
257 new_x = new_x_z[new_k];
258 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
260 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
261 #ifndef AMREX_USE_GPU
263 amrex::Print() <<
" orig_x_sub_p = [";
264 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
265 amrex::Print() <<
"]" << std::endl;
266 amrex::Print() <<
" orig_y_sub_p = [";
267 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
268 amrex::Print() <<
"]" << std::endl;
273 }
else if (order%2 == 0) {
274 if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_n-1)) {
277 amrex::Real new_y_l, new_y_r;
279 int ksta = kl-(order/2-1);
280 int kend = ksta+order;
281 #ifndef AMREX_USE_GPU
282 int ksize = kend-ksta;
283 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;
286 amrex::GpuArray<amrex::Real,10> orig_x_sub;
287 amrex::GpuArray<amrex::Real,10> orig_y_sub;
288 amrex::Real* orig_x_sub_p = orig_x_sub.data();
289 amrex::Real* orig_y_sub_p = orig_y_sub.data();
291 new_x = new_x_p[new_k];
292 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
294 new_x = new_x_z[new_k];
295 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
297 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
298 #ifndef AMREX_USE_GPU
300 amrex::Print() <<
" orig_x_sub_p = [";
301 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
302 amrex::Print() <<
"]" << std::endl;
303 amrex::Print() <<
" orig_y_sub_p = [";
304 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
305 amrex::Print() <<
"]" << std::endl;
311 int ksta = kl-order/2;
312 int kend = ksta+order;
313 #ifndef AMREX_USE_GPU
314 int ksize = kend-ksta;
315 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;
318 amrex::GpuArray<amrex::Real,10> orig_x_sub;
319 amrex::GpuArray<amrex::Real,10> orig_y_sub;
320 amrex::Real* orig_x_sub_p = orig_x_sub.data();
321 amrex::Real* orig_y_sub_p = orig_y_sub.data();
323 new_x = new_x_p[new_k];
324 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
326 new_x = new_x_z[new_k];
327 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
329 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
330 #ifndef AMREX_USE_GPU
332 amrex::Print() <<
" orig_x_sub_p = [";
333 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
334 amrex::Print() <<
"]" << std::endl;
335 amrex::Print() <<
" orig_y_sub_p = [";
336 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
337 amrex::Print() <<
"]" << std::endl;
342 new_y[new_k] = 0.5*(new_y_l+new_y_r);
343 }
else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_n-1)) {
344 int ksta = kl-(order/2-1);
345 int kend = ksta+order;
346 #ifndef AMREX_USE_GPU
347 int ksize = kend-ksta;
348 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;
351 amrex::GpuArray<amrex::Real,10> orig_x_sub;
352 amrex::GpuArray<amrex::Real,10> orig_y_sub;
353 amrex::Real* orig_x_sub_p = orig_x_sub.data();
354 amrex::Real* orig_y_sub_p = orig_y_sub.data();
356 new_x = new_x_p[new_k];
357 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
359 new_x = new_x_z[new_k];
360 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
362 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
363 #ifndef AMREX_USE_GPU
365 amrex::Print() <<
" orig_x_sub_p = [";
366 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
367 amrex::Print() <<
"]" << std::endl;
368 amrex::Print() <<
" orig_y_sub_p = [";
369 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
370 amrex::Print() <<
"]" << std::endl;
373 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
374 }
else if ((kl-order/2 >= 0) && (kr+order/2-1 <= orig_n-1)) {
375 int ksta = kl-order/2;
376 int kend = ksta+order;
377 #ifndef AMREX_USE_GPU
378 int ksize = kend-ksta;
379 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;
382 amrex::GpuArray<amrex::Real,10> orig_x_sub;
383 amrex::GpuArray<amrex::Real,10> orig_y_sub;
384 amrex::Real* orig_x_sub_p = orig_x_sub.data();
385 amrex::Real* orig_y_sub_p = orig_y_sub.data();
387 new_x = new_x_p[new_k];
388 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
390 new_x = new_x_z[new_k];
391 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
393 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
394 #ifndef AMREX_USE_GPU
396 amrex::Print() <<
" orig_x_sub_p = [";
397 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
398 amrex::Print() <<
"]" << std::endl;
399 amrex::Print() <<
" orig_y_sub_p = [";
400 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
401 amrex::Print() <<
"]" << std::endl;
404 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
409 #ifndef AMREX_USE_GPU
410 int ksize = kend-ksta+1;
411 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;
414 amrex::GpuArray<amrex::Real,2> orig_x_sub;
415 amrex::GpuArray<amrex::Real,2> orig_y_sub;
416 amrex::Real* orig_x_sub_p = orig_x_sub.data();
417 amrex::Real* orig_y_sub_p = orig_y_sub.data();
419 new_x = new_x_p[new_k];
420 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
422 new_x = new_x_z[new_k];
423 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
425 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
426 #ifndef AMREX_USE_GPU
428 amrex::Print() <<
" orig_x_sub_p = [";
429 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
430 amrex::Print() <<
"]" << std::endl;
431 amrex::Print() <<
" orig_y_sub_p = [";
432 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
433 amrex::Print() <<
"]" << std::endl;
439 #ifndef AMREX_USE_GPU
440 if (debug) amrex::Print() <<
" new_y[" << new_k <<
"]=" << new_y[new_k] << std::endl;
458 const bool& metgrid_use_sfc,
459 const bool& exp_interp,
460 const bool& metgrid_retain_sfc,
461 const amrex::Real& metgrid_proximity,
462 const int& metgrid_order,
463 const int& metgrid_force_sfc_k,
470 const amrex::Array4<amrex::Real const>& orig_z_full,
471 const amrex::Array4<amrex::Real const>& orig_data,
472 const amrex::Array4<amrex::Real const>& new_z_full,
473 const amrex::Array4<amrex::Real>& new_data_full,
474 const bool& update_bc_data,
475 const amrex::Array4<amrex::Real>& bc_data_xlo,
476 const amrex::Array4<amrex::Real>& bc_data_xhi,
477 const amrex::Array4<amrex::Real>& bc_data_ylo,
478 const amrex::Array4<amrex::Real>& bc_data_yhi,
479 const amrex::Box& bx_xlo,
480 const amrex::Box& bx_xhi,
481 const amrex::Box& bx_ylo,
482 const amrex::Box& bx_yhi,
483 const amrex::Array4<const int>& mask)
489 int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
490 int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
491 int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
492 int kmax_new = amrex::ubound(amrex::Box(new_z_full)).z;
494 AMREX_ASSERT(kmax_orig < 256);
495 AMREX_ASSERT(kmax_new < 256);
497 amrex::GpuArray<amrex::Real,256> new_z;
498 amrex::GpuArray<amrex::Real,256> new_p;
499 amrex::GpuArray<amrex::Real,256> new_data;
500 amrex::Real* new_z_p = new_z.data();
501 amrex::Real* new_p_p = new_p.data();
502 amrex::Real* new_data_p = new_data.data();
503 for (
int k=0; k < kmax_new; k++) {
505 new_z_p[k] = 0.25*(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));
506 }
else if (stag ==
'Y') {
507 new_z_p[k] = 0.25*(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));
508 }
else if (stag ==
'M') {
509 new_z_p[k] = 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 )+
510 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));
515 amrex::GpuArray<amrex::Real,256> orig_z;
516 amrex::Real* orig_z_p = orig_z.data();
517 for (
int k=0; k < kmax_orig; k++) {
519 orig_z_p[k] = orig_z_full(i,j,k);
520 }
else if (stag ==
'X') {
522 orig_z_p[k] = orig_z_full(i,j,k);
523 }
else if (i == imax_orig) {
524 orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
526 orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
528 }
else if (stag ==
'Y') {
530 orig_z_p[k] = orig_z_full(i,j,k);
531 }
else if (j == jmax_orig) {
532 orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
534 orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
540 bool flip_data_required =
false;
541 if (orig_z[1] > orig_z[kmax_orig-1]) flip_data_required =
true;
542 if (flip_data_required) amrex::Abort(
"metgrid initialization flip_data_required. Not yet implemented.");
548 for (
int k=1; k < kmax_orig; k++) {
549 if (orig_z_p[k] > orig_z_p[0]) {
556 amrex::GpuArray<amrex::Real,256> ordered_z;
557 amrex::GpuArray<amrex::Real,256> ordered_data;
558 amrex::Real* ordered_z_p = ordered_z.data();
559 amrex::Real* ordered_data_p = ordered_data.data();
560 if (k_above_sfc > 1) {
566 for (
int k=1; k < k_above_sfc; k++) {
567 ordered_z_p[count] = orig_z_p[k];
568 ordered_data_p[count] = orig_data(i,j,k);
582 if (Pl-Pu < metgrid_proximity) {
587 ordered_z_p[count] = orig_z_p[0];
588 ordered_data_p[count] = orig_data(i,j,0);
593 int knext = k_above_sfc;
599 if (metgrid_force_sfc_k > 0) {
600 for (
int k=k_above_sfc; k < kmax_orig; k++) {
601 if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
612 if (Pl-Pu < metgrid_proximity) {
617 for (
int k=knext; k < kmax_orig; k++) {
618 ordered_z_p[count] = orig_z_p[k];
619 ordered_data_p[count] = orig_data(i,j,k);
628 ordered_z_p[0] = orig_z[0];
629 ordered_data_p[0] = orig_data(i,j,0);
635 if (metgrid_force_sfc_k > 0) {
636 for (
int k=knext; k < kmax_orig; k++) {
637 if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k]) {
646 for (
int k=knext; k < kmax_orig; k++) {
650 if (Pl-Pu < metgrid_proximity) {
653 ordered_z_p[count] = orig_z_p[k];
654 ordered_data_p[count] = orig_data(i,j,k);
660 int ksta(0), kend(0);
661 if (metgrid_use_below_sfc && metgrid_use_sfc) {
665 }
else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
668 for (
int k=0; k < kmax_orig; k++) {
669 if (ordered_z_p[k] == orig_z_p[0]) {
674 for (
int k=ksfc; k < kmax_orig-1; k++) {
675 ordered_z_p[k] = ordered_z_p[k+1];
676 ordered_data_p[k] = ordered_data_p[k+1];
680 }
else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
683 for (
int k=0; k < kend_order; k++) {
684 if (ordered_z_p[k] == orig_z_p[0]) {
693 amrex::Abort(
"metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
709 amrex::GpuArray<amrex::Real,256> ordered_p;
710 amrex::Real* ordered_p_p = ordered_p.data();
711 for (
int k=0; k < kend_order; k++) {
717 amrex::GpuArray<amrex::Real,256> final_z;
718 amrex::GpuArray<amrex::Real,256> final_p;
719 amrex::GpuArray<amrex::Real,256> final_data;
720 amrex::Real* final_z_p = final_z.data();
721 amrex::Real* final_p_p = final_p.data();
722 amrex::Real* final_data_p = final_data.data();
723 final_z_p[0] = ordered_z[ksta];
724 final_p_p[0] = ordered_p[ksta];
725 final_data_p[0] = ordered_data[ksta];
726 for (
int k=ksta+1; k <= kend; k++) {
727 if ((final_p_p[new_n]-ordered_p_p[k]) < metgrid_proximity) {
731 final_z_p[new_n] = ordered_z_p[k];
732 final_p_p[new_n] = ordered_p_p[k];
733 final_data_p[new_n] = ordered_data_p[k];
755 if (metgrid_retain_sfc) new_data[0] = ordered_data[0];
758 for (
int k=0; k < kmax_new; k++) {
759 if (mask(i,j,k) && update_bc_data && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = new_data[k];
760 if (mask(i,j,k) && update_bc_data && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = new_data[k];
761 if (mask(i,j,k) && update_bc_data && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = new_data[k];
762 if (mask(i,j,k) && update_bc_data && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = new_data[k];
763 if (itime == 0) new_data_full(i,j,k,src_comp) = new_data[k];
776 const amrex::Array4<amrex::Real const>& orig_z,
777 const amrex::Array4<amrex::Real const>& orig_data,
778 const amrex::Array4<amrex::Real const>& new_z)
781 int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
782 int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
783 int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
787 z = 0.25*(new_z(i,j,k)+new_z(i,j+1,k)+new_z(i,j,k+1)+new_z(i,j+1,k+1));
789 else if (stag ==
'Y') {
790 z = 0.25*(new_z(i,j,k)+new_z(i+1,j,k)+new_z(i,j,k+1)+new_z(i+1,j,k+1));
792 else if (stag ==
'M') {
793 z = 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 )+
794 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));
800 amrex::Real dzlow = 1.0e12;
801 amrex::Real dzhi0 = -1.0e12;
802 for (
int kk = 0; kk < kmax_orig; kk++) {
803 amrex::Real orig_z_stag = 0.0;
805 orig_z_stag = orig_z(i,j,kk);
809 orig_z_stag = orig_z(i,j,kk);
811 else if (i == imax_orig) {
812 orig_z_stag = orig_z(imax_orig-1,j,kk);
815 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
818 else if (stag ==
'Y') {
820 orig_z_stag = orig_z(i,j,kk);
822 else if (j == jmax_orig) {
823 orig_z_stag = orig_z(i,jmax_orig-1,kk);
826 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
830 amrex::Real dz =
z - orig_z_stag;
831 if ((dz < 0.0) && (dz > dzhi0)) {
836 if ((dz >= 0.0) && (dz < dzlow)) {
846 amrex::Real dzhi1 = -1.0e12;
847 for (
int kk = 0; kk < kmax_orig; kk++) {
848 amrex::Real orig_z_stag = 0.0;
850 orig_z_stag = orig_z(i,j,kk);
852 else if (stag ==
'X') {
854 orig_z_stag = orig_z(i,j,kk);
856 else if (i == imax_orig) {
857 orig_z_stag = orig_z(imax_orig-1,j,kk);
860 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
863 else if (stag ==
'Y') {
865 orig_z_stag = orig_z(i,j,kk);
867 else if (j == jmax_orig) {
868 orig_z_stag = orig_z(i,jmax_orig-1,kk);
871 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
874 amrex::Real dz =
z - orig_z_stag;
875 if ((dz < 0.0) && (dz > dzhi1) && (kk != khi0)) {
881 amrex::Real y0 = orig_data(i,j,khi0,src_comp);
882 amrex::Real y1 = orig_data(i,j,khi1,src_comp);
883 return ( y0-(y1-y0)/(z1-z0)*(z0-
z) );
886 }
else if (khi0 == -1) {
890 z0 = orig_z(i,j,khi0);
892 else if (stag ==
'X') {
894 z0 = orig_z(i,j,khi0);
896 else if (i == imax_orig) {
897 z0 = orig_z(imax_orig-1,j,khi0);
900 z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
903 else if (stag ==
'Y') {
905 z0 = orig_z(i,j,khi0);
907 else if (j == jmax_orig) {
908 z0 = orig_z(i,jmax_orig-1,khi0);
911 z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
914 amrex::Real y0 = orig_data(i,j,khi0,src_comp);
915 amrex::Real y1 = orig_data(i,j,khi1,src_comp);
916 return ( y0+(y1-y0)/(z1-z0)*(
z-z0) );
919 amrex::Real y0 = orig_data(i,j,klow,src_comp);
920 amrex::Real y1 = orig_data(i,j,khi0,src_comp);
921 return ( y0+(y1-y0)/(z1-z0)*(
z-z0) );
932 const amrex::Array4<amrex::Real const>& rhum,
933 const amrex::Array4<amrex::Real const>& temp,
934 const amrex::Array4<amrex::Real const>&
pres,
935 const amrex::Array4<amrex::Real>& mxrat)
937 amrex::Real qv_max_p_safe = 10000.0;
938 amrex::Real qv_max_flag = 1.0e-5;
939 amrex::Real qv_max_value = 3.0e-6;
940 amrex::Real qv_min_p_safe = 110000.0;
941 amrex::Real qv_min_flag = 1.0e-6;
942 amrex::Real qv_min_value = 1.0e-6;
943 amrex::Real eps = 0.622;
944 amrex::Real svp1 = 0.6112;
945 amrex::Real svp2 = 17.67;
946 amrex::Real svp3 = 29.65;
947 amrex::Real svpt0 = 273.15;
949 if (temp(i,j,k) != 0.0) {
950 amrex::Real es=0.01*rhum(i,j,k)*svp1*10.0*exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3));
951 if (es >=
pres(i,j,k)/100.0) {
953 mxrat(i,j,k) = std::pow(10.0, -6);
956 mxrat(i,j,k) = amrex::max(eps*es/(
pres(i,j,k)/100.0-es), 1.0e-6);
962 mxrat(i,j,k) = 1.0e-6;
970 if (
pres(i,j,k) < qv_max_p_safe) {
971 if (mxrat(i,j,k) > qv_max_flag) {
972 mxrat(i,j,k) = qv_max_value;
975 if (
pres(i,j,k) < qv_min_p_safe) {
976 if (mxrat(i,j,k) < qv_min_flag) {
977 mxrat(i,j,k) = qv_min_value;
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
constexpr amrex::Real p_0
Definition: ERF_Constants.H:18
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10
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 &src_comp, const int &itime, 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, const bool &update_bc_data, const amrex::Array4< amrex::Real > &bc_data_xlo, const amrex::Array4< amrex::Real > &bc_data_xhi, const amrex::Array4< amrex::Real > &bc_data_ylo, const amrex::Array4< amrex::Real > &bc_data_yhi, const amrex::Box &bx_xlo, const amrex::Box &bx_xhi, const amrex::Box &bx_ylo, const amrex::Box &bx_yhi, const amrex::Array4< const int > &mask)
Definition: ERF_MetgridUtils.H:457
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:114
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:771
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:448
void init_state_from_metgrid(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::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::Vector< amrex::FArrayBox > &NC_hgt_fab, const amrex::Vector< amrex::FArrayBox > &NC_ght_fab, const amrex::Vector< amrex::FArrayBox > &NC_xvel_fab, const amrex::Vector< amrex::FArrayBox > &NC_yvel_fab, const amrex::Vector< amrex::FArrayBox > &NC_temp_fab, const amrex::Vector< amrex::FArrayBox > &NC_rhum_fab, const amrex::Vector< amrex::FArrayBox > &NC_pres_fab, amrex::FArrayBox &p_interp_fab, amrex::FArrayBox &t_interp_fab, amrex::FArrayBox &theta_fab, amrex::FArrayBox &mxrat_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)
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, const amrex::Array4< amrex::Real > &mxrat)
Definition: ERF_MetgridUtils.H:929
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_setup(char var_type, const bool &exp_interp, const int &orig_n, const int &new_n, 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:144
void init_base_state_from_metgrid(const bool use_moisture, const bool metgrid_debug_psfc, const amrex::Real l_rdOcp, const amrex::Box &valid_bx, const amrex::Vector< int > &flag_psfc, amrex::FArrayBox &state, 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_cc_fab, const amrex::Vector< amrex::FArrayBox > &NC_psfc_fab)
void read_from_metgrid(int lev, const amrex::Box &domain, const std::string &fname, std::string &NC_dateTime, amrex::Real &NC_epochTime, int &flag_psfc, int &flag_msf, int &flag_sst, 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_LAT_fab, amrex::FArrayBox &NC_LON_fab, amrex::IArrayBox &NC_lmask_iab, amrex::Geometry &geom)
void init_terrain_from_metgrid(amrex::FArrayBox &z_phys_nd_fab, const amrex::Vector< amrex::FArrayBox > &NC_hgt_fab)
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, const amrex::Vector< amrex::FArrayBox > &NC_MSFU_fab, const amrex::Vector< amrex::FArrayBox > &NC_MSFV_fab, const amrex::Vector< amrex::FArrayBox > &NC_MSFM_fab)
@ pres
Definition: ERF_Kessler.H:25