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,
25 amrex::FArrayBox& NC_xvel_fab,
26 amrex::FArrayBox& NC_yvel_fab,
27 amrex::FArrayBox& NC_temp_fab,
28 amrex::FArrayBox& NC_rhum_fab,
29 amrex::FArrayBox& NC_pres_fab,
30 amrex::FArrayBox& NC_ght_fab,
31 amrex::FArrayBox& NC_hgt_fab,
32 amrex::FArrayBox& NC_psfc_fab,
33 amrex::FArrayBox& NC_msfu_fab,
34 amrex::FArrayBox& NC_msfv_fab,
35 amrex::FArrayBox& NC_msfm_fab,
36 amrex::FArrayBox& NC_sst_fab,
37 amrex::FArrayBox& NC_tsk_fab,
38 amrex::FArrayBox& NC_LAT_fab,
39 amrex::FArrayBox& NC_LON_fab,
40 amrex::IArrayBox& NC_lmask_iab,
41 amrex::Geometry& geom);
45 const amrex::Vector<amrex::FArrayBox>& NC_hgt_fab);
49 const bool interp_theta,
50 const bool metgrid_debug_quiescent,
51 const bool metgrid_debug_isothermal,
52 const bool metgrid_debug_dry,
53 const bool metgrid_basic_linear,
54 const bool metgrid_use_below_sfc,
55 const bool metgrid_use_sfc,
56 const bool metgrid_retain_sfc,
57 const amrex::Real metgrid_proximity,
58 const int metgrid_order,
59 const int metgrid_metgrid_force_sfc_k,
60 const amrex::Real l_rdOcp,
64 amrex::FArrayBox& state_fab,
65 amrex::FArrayBox& x_vel_fab,
66 amrex::FArrayBox& y_vel_fab,
67 amrex::FArrayBox& z_vel_fab,
68 amrex::FArrayBox& z_phys_nd_fab,
69 const amrex::Vector<amrex::FArrayBox>& NC_hgt_fab,
70 const amrex::Vector<amrex::FArrayBox>& NC_ght_fab,
71 const amrex::Vector<amrex::FArrayBox>& NC_xvel_fab,
72 const amrex::Vector<amrex::FArrayBox>& NC_yvel_fab,
73 const amrex::Vector<amrex::FArrayBox>& NC_temp_fab,
74 const amrex::Vector<amrex::FArrayBox>& NC_rhum_fab,
75 const amrex::Vector<amrex::FArrayBox>& NC_pres_fab,
76 amrex::FArrayBox& p_interp_fab,
77 amrex::FArrayBox& t_interp_fab,
78 amrex::FArrayBox& theta_fab,
79 amrex::FArrayBox& mxrat_fab,
80 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xlo,
81 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xhi,
82 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_ylo,
83 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_yhi,
84 const amrex::Array4<const int>& mask_c_arr,
85 const amrex::Array4<const int>& mask_u_arr,
86 const amrex::Array4<const int>& mask_v_arr);
90 amrex::FArrayBox& msfu_fab,
91 amrex::FArrayBox& msfv_fab,
92 amrex::FArrayBox& msfm_fab,
94 const amrex::Vector<amrex::FArrayBox>& NC_MSFU_fab,
95 const amrex::Vector<amrex::FArrayBox>& NC_MSFV_fab,
96 const amrex::Vector<amrex::FArrayBox>& NC_MSFM_fab);
100 const bool metgrid_debug_psfc,
101 const amrex::Real l_rdOcp,
102 const amrex::Box& valid_bx,
103 const amrex::Vector<int>& flag_psfc,
104 amrex::FArrayBox& state,
105 amrex::FArrayBox& r_hse_fab,
106 amrex::FArrayBox& p_hse_fab,
107 amrex::FArrayBox& pi_hse_fab,
108 amrex::FArrayBox& th_hse_fab,
109 amrex::FArrayBox& qv_hse_fab,
110 amrex::FArrayBox& z_phys_cc_fab,
111 const amrex::Vector<amrex::FArrayBox>& NC_psfc_fab);
127 amrex::Real Px = 0.0;
128 for (
int i=0; i <= order; i++) {
131 for (
int k=0; k <= order; k++) {
132 if (k == i)
continue;
147 const bool& exp_interp,
153 amrex::Real* orig_x_z,
154 amrex::Real* orig_x_p,
156 amrex::Real* new_x_z,
157 amrex::Real* new_x_p,
161 amrex::Real CRC_const1 = 11880.516;
162 amrex::Real CRC_const2 = 0.1902632;
163 amrex::Real CRC_const3 = 0.0065;
165 amrex::ignore_unused(i,j);
166 #ifndef AMREX_USE_GPU
170 for (
int new_k=0; new_k < new_n; new_k++) {
171 #ifndef AMREX_USE_GPU
172 if (debug) amrex::Print() <<
"new_k=" << new_k;
176 bool extrapolating =
true;
178 for (
int ko=0; ko < orig_n-1; ko++) {
179 amrex::Real a = new_x_z[new_k]-orig_x_z[ko];
180 amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1];
184 extrapolating =
false;
190 if (var_type ==
'T') {
193 amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[0];
194 amrex::Real avg_of_extrap_p = 0.5*(new_x_p[new_k]+orig_x_p[0]);
195 amrex::Real temp_extrap_starting_point = orig_y[0]*std::pow(orig_x_p[0]/100000.0,
R_d/
Cp_d);
196 amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/100.0, CRC_const2-1.0);
197 amrex::Real dZ = dZdP*(depth_of_extrap_in_p/100.0);
198 amrex::Real dT = dZ*CRC_const3;
199 new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(100000.0/new_x_p[new_k],
R_d/
Cp_d);
203 new_y[new_k] = orig_y[0];
209 if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_n-1)) {
211 int ksta = kl-(((order+1)/2)-1);
212 int kend = ksta+order;
213 #ifndef AMREX_USE_GPU
214 int ksize = kend-ksta;
215 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;
218 amrex::GpuArray<amrex::Real,9> orig_x_sub;
219 amrex::Real* orig_x_sub_p = orig_x_sub.data();
221 new_x = new_x_p[new_k];
222 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
224 new_x = new_x_z[new_k];
225 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
227 amrex::GpuArray<amrex::Real,9> orig_y_sub;
228 amrex::Real* orig_y_sub_p = orig_y_sub.data();
229 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
230 #ifndef AMREX_USE_GPU
232 amrex::Print() <<
" orig_x_sub_p = [";
233 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
234 amrex::Print() <<
"]" << std::endl;
235 amrex::Print() <<
" orig_y_sub_p = [";
236 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
237 amrex::Print() <<
"]" << std::endl;
240 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
246 #ifndef AMREX_USE_GPU
247 int ksize = kend-ksta+1;
248 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;
251 amrex::GpuArray<amrex::Real,2> orig_x_sub;
252 amrex::GpuArray<amrex::Real,2> orig_y_sub;
253 amrex::Real* orig_x_sub_p = orig_x_sub.data();
254 amrex::Real* orig_y_sub_p = orig_y_sub.data();
256 new_x = new_x_p[new_k];
257 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
259 new_x = new_x_z[new_k];
260 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
262 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
263 #ifndef AMREX_USE_GPU
265 amrex::Print() <<
" orig_x_sub_p = [";
266 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
267 amrex::Print() <<
"]" << std::endl;
268 amrex::Print() <<
" orig_y_sub_p = [";
269 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
270 amrex::Print() <<
"]" << std::endl;
275 }
else if (order%2 == 0) {
276 if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_n-1)) {
279 amrex::Real new_y_l, new_y_r;
281 int ksta = kl-(order/2-1);
282 int kend = ksta+order;
283 #ifndef AMREX_USE_GPU
284 int ksize = kend-ksta;
285 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;
288 amrex::GpuArray<amrex::Real,10> orig_x_sub;
289 amrex::GpuArray<amrex::Real,10> orig_y_sub;
290 amrex::Real* orig_x_sub_p = orig_x_sub.data();
291 amrex::Real* orig_y_sub_p = orig_y_sub.data();
293 new_x = new_x_p[new_k];
294 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
296 new_x = new_x_z[new_k];
297 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
299 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
300 #ifndef AMREX_USE_GPU
302 amrex::Print() <<
" orig_x_sub_p = [";
303 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
304 amrex::Print() <<
"]" << std::endl;
305 amrex::Print() <<
" orig_y_sub_p = [";
306 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
307 amrex::Print() <<
"]" << std::endl;
313 int ksta = kl-order/2;
314 int kend = ksta+order;
315 #ifndef AMREX_USE_GPU
316 int ksize = kend-ksta;
317 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;
320 amrex::GpuArray<amrex::Real,10> orig_x_sub;
321 amrex::GpuArray<amrex::Real,10> orig_y_sub;
322 amrex::Real* orig_x_sub_p = orig_x_sub.data();
323 amrex::Real* orig_y_sub_p = orig_y_sub.data();
325 new_x = new_x_p[new_k];
326 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
328 new_x = new_x_z[new_k];
329 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
331 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
332 #ifndef AMREX_USE_GPU
334 amrex::Print() <<
" orig_x_sub_p = [";
335 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
336 amrex::Print() <<
"]" << std::endl;
337 amrex::Print() <<
" orig_y_sub_p = [";
338 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
339 amrex::Print() <<
"]" << std::endl;
344 new_y[new_k] = 0.5*(new_y_l+new_y_r);
345 }
else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_n-1)) {
346 int ksta = kl-(order/2-1);
347 int kend = ksta+order;
348 #ifndef AMREX_USE_GPU
349 int ksize = kend-ksta;
350 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;
353 amrex::GpuArray<amrex::Real,10> orig_x_sub;
354 amrex::GpuArray<amrex::Real,10> orig_y_sub;
355 amrex::Real* orig_x_sub_p = orig_x_sub.data();
356 amrex::Real* orig_y_sub_p = orig_y_sub.data();
358 new_x = new_x_p[new_k];
359 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
361 new_x = new_x_z[new_k];
362 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
364 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
365 #ifndef AMREX_USE_GPU
367 amrex::Print() <<
" orig_x_sub_p = [";
368 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
369 amrex::Print() <<
"]" << std::endl;
370 amrex::Print() <<
" orig_y_sub_p = [";
371 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
372 amrex::Print() <<
"]" << std::endl;
375 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
376 }
else if ((kl-order/2 >= 0) && (kr+order/2-1 <= orig_n-1)) {
377 int ksta = kl-order/2;
378 int kend = ksta+order;
379 #ifndef AMREX_USE_GPU
380 int ksize = kend-ksta;
381 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;
384 amrex::GpuArray<amrex::Real,10> orig_x_sub;
385 amrex::GpuArray<amrex::Real,10> orig_y_sub;
386 amrex::Real* orig_x_sub_p = orig_x_sub.data();
387 amrex::Real* orig_y_sub_p = orig_y_sub.data();
389 new_x = new_x_p[new_k];
390 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
392 new_x = new_x_z[new_k];
393 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
395 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
396 #ifndef AMREX_USE_GPU
398 amrex::Print() <<
" orig_x_sub_p = [";
399 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
400 amrex::Print() <<
"]" << std::endl;
401 amrex::Print() <<
" orig_y_sub_p = [";
402 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
403 amrex::Print() <<
"]" << std::endl;
406 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
411 #ifndef AMREX_USE_GPU
412 int ksize = kend-ksta+1;
413 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;
416 amrex::GpuArray<amrex::Real,2> orig_x_sub;
417 amrex::GpuArray<amrex::Real,2> orig_y_sub;
418 amrex::Real* orig_x_sub_p = orig_x_sub.data();
419 amrex::Real* orig_y_sub_p = orig_y_sub.data();
421 new_x = new_x_p[new_k];
422 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
424 new_x = new_x_z[new_k];
425 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
427 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
428 #ifndef AMREX_USE_GPU
430 amrex::Print() <<
" orig_x_sub_p = [";
431 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
432 amrex::Print() <<
"]" << std::endl;
433 amrex::Print() <<
" orig_y_sub_p = [";
434 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
435 amrex::Print() <<
"]" << std::endl;
441 #ifndef AMREX_USE_GPU
442 if (debug) amrex::Print() <<
" new_y[" << new_k <<
"]=" << new_y[new_k] << std::endl;
460 const bool& metgrid_use_sfc,
461 const bool& exp_interp,
462 const bool& metgrid_retain_sfc,
463 const amrex::Real& metgrid_proximity,
464 const int& metgrid_order,
465 const int& metgrid_force_sfc_k,
472 const amrex::Array4<amrex::Real const>& orig_z_full,
473 const amrex::Array4<amrex::Real const>& orig_data,
474 const amrex::Array4<amrex::Real const>& new_z_full,
475 const amrex::Array4<amrex::Real>& new_data_full,
476 const bool& update_bc_data,
477 const amrex::Array4<amrex::Real>& bc_data_xlo,
478 const amrex::Array4<amrex::Real>& bc_data_xhi,
479 const amrex::Array4<amrex::Real>& bc_data_ylo,
480 const amrex::Array4<amrex::Real>& bc_data_yhi,
481 const amrex::Box& bx_xlo,
482 const amrex::Box& bx_xhi,
483 const amrex::Box& bx_ylo,
484 const amrex::Box& bx_yhi,
485 const amrex::Array4<const int>& mask)
491 int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
492 int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
493 int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
494 int kmax_new = amrex::ubound(amrex::Box(new_z_full)).z;
496 AMREX_ASSERT(kmax_orig < 256);
497 AMREX_ASSERT(kmax_new < 256);
499 amrex::GpuArray<amrex::Real,256> new_z;
500 amrex::GpuArray<amrex::Real,256> new_p;
501 amrex::GpuArray<amrex::Real,256> new_data;
502 amrex::Real* new_z_p = new_z.data();
503 amrex::Real* new_p_p = new_p.data();
504 amrex::Real* new_data_p = new_data.data();
505 for (
int k=0; k < kmax_new; k++) {
507 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));
508 }
else if (stag ==
'Y') {
509 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));
510 }
else if (stag ==
'M') {
511 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 )+
512 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));
517 amrex::GpuArray<amrex::Real,256> orig_z;
518 amrex::Real* orig_z_p = orig_z.data();
519 for (
int k=0; k < kmax_orig; k++) {
521 orig_z_p[k] = orig_z_full(i,j,k);
522 }
else if (stag ==
'X') {
524 orig_z_p[k] = orig_z_full(i,j,k);
525 }
else if (i == imax_orig) {
526 orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
528 orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
530 }
else if (stag ==
'Y') {
532 orig_z_p[k] = orig_z_full(i,j,k);
533 }
else if (j == jmax_orig) {
534 orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
536 orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
542 bool flip_data_required =
false;
543 if (orig_z[1] > orig_z[kmax_orig-1]) flip_data_required =
true;
544 if (flip_data_required) amrex::Abort(
"metgrid initialization flip_data_required. Not yet implemented.");
550 for (
int k=1; k < kmax_orig; k++) {
551 if (orig_z_p[k] > orig_z_p[0]) {
558 amrex::GpuArray<amrex::Real,256> ordered_z;
559 amrex::GpuArray<amrex::Real,256> ordered_data;
560 amrex::Real* ordered_z_p = ordered_z.data();
561 amrex::Real* ordered_data_p = ordered_data.data();
562 if (k_above_sfc > 1) {
568 for (
int k=1; k < k_above_sfc; k++) {
569 ordered_z_p[count] = orig_z_p[k];
570 ordered_data_p[count] = orig_data(i,j,k);
584 if (Pl-Pu < metgrid_proximity) {
589 ordered_z_p[count] = orig_z_p[0];
590 ordered_data_p[count] = orig_data(i,j,0);
595 int knext = k_above_sfc;
601 if (metgrid_force_sfc_k > 0) {
602 for (
int k=k_above_sfc; k < kmax_orig; k++) {
603 if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
614 if (Pl-Pu < metgrid_proximity) {
619 for (
int k=knext; k < kmax_orig; k++) {
620 ordered_z_p[count] = orig_z_p[k];
621 ordered_data_p[count] = orig_data(i,j,k);
630 ordered_z_p[0] = orig_z[0];
631 ordered_data_p[0] = orig_data(i,j,0);
637 if (metgrid_force_sfc_k > 0) {
638 for (
int k=knext; k < kmax_orig; k++) {
639 if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k]) {
648 for (
int k=knext; k < kmax_orig; k++) {
652 if (Pl-Pu < metgrid_proximity) {
655 ordered_z_p[count] = orig_z_p[k];
656 ordered_data_p[count] = orig_data(i,j,k);
662 int ksta(0), kend(0);
663 if (metgrid_use_below_sfc && metgrid_use_sfc) {
667 }
else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
670 for (
int k=0; k < kmax_orig; k++) {
671 if (ordered_z_p[k] == orig_z_p[0]) {
676 for (
int k=ksfc; k < kmax_orig-1; k++) {
677 ordered_z_p[k] = ordered_z_p[k+1];
678 ordered_data_p[k] = ordered_data_p[k+1];
682 }
else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
685 for (
int k=0; k < kend_order; k++) {
686 if (ordered_z_p[k] == orig_z_p[0]) {
695 amrex::Abort(
"metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
711 amrex::GpuArray<amrex::Real,256> ordered_p;
712 amrex::Real* ordered_p_p = ordered_p.data();
713 for (
int k=0; k < kend_order; k++) {
719 amrex::GpuArray<amrex::Real,256> final_z;
720 amrex::GpuArray<amrex::Real,256> final_p;
721 amrex::GpuArray<amrex::Real,256> final_data;
722 amrex::Real* final_z_p = final_z.data();
723 amrex::Real* final_p_p = final_p.data();
724 amrex::Real* final_data_p = final_data.data();
725 final_z_p[0] = ordered_z[ksta];
726 final_p_p[0] = ordered_p[ksta];
727 final_data_p[0] = ordered_data[ksta];
728 for (
int k=ksta+1; k <= kend; k++) {
729 if ((final_p_p[new_n]-ordered_p_p[k]) < metgrid_proximity) {
733 final_z_p[new_n] = ordered_z_p[k];
734 final_p_p[new_n] = ordered_p_p[k];
735 final_data_p[new_n] = ordered_data_p[k];
757 if (metgrid_retain_sfc) new_data[0] = ordered_data[0];
760 for (
int k=0; k < kmax_new; k++) {
761 if (mask(i,j,k) && update_bc_data && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = new_data[k];
762 if (mask(i,j,k) && update_bc_data && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = new_data[k];
763 if (mask(i,j,k) && update_bc_data && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = new_data[k];
764 if (mask(i,j,k) && update_bc_data && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = new_data[k];
765 if (itime == 0) new_data_full(i,j,k,src_comp) = new_data[k];
778 const amrex::Array4<amrex::Real const>& orig_z,
779 const amrex::Array4<amrex::Real const>& orig_data,
780 const amrex::Array4<amrex::Real const>& new_z)
783 int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
784 int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
785 int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
789 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));
791 else if (stag ==
'Y') {
792 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));
794 else if (stag ==
'M') {
795 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 )+
796 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));
802 amrex::Real dzlow = 1.0e12;
803 amrex::Real dzhi0 = -1.0e12;
804 for (
int kk = 0; kk < kmax_orig; kk++) {
805 amrex::Real orig_z_stag = 0.0;
807 orig_z_stag = orig_z(i,j,kk);
811 orig_z_stag = orig_z(i,j,kk);
813 else if (i == imax_orig) {
814 orig_z_stag = orig_z(imax_orig-1,j,kk);
817 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
820 else if (stag ==
'Y') {
822 orig_z_stag = orig_z(i,j,kk);
824 else if (j == jmax_orig) {
825 orig_z_stag = orig_z(i,jmax_orig-1,kk);
828 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
832 amrex::Real dz =
z - orig_z_stag;
833 if ((dz < 0.0) && (dz > dzhi0)) {
838 if ((dz >= 0.0) && (dz < dzlow)) {
848 amrex::Real dzhi1 = -1.0e12;
849 for (
int kk = 0; kk < kmax_orig; kk++) {
850 amrex::Real orig_z_stag = 0.0;
852 orig_z_stag = orig_z(i,j,kk);
854 else if (stag ==
'X') {
856 orig_z_stag = orig_z(i,j,kk);
858 else if (i == imax_orig) {
859 orig_z_stag = orig_z(imax_orig-1,j,kk);
862 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
865 else if (stag ==
'Y') {
867 orig_z_stag = orig_z(i,j,kk);
869 else if (j == jmax_orig) {
870 orig_z_stag = orig_z(i,jmax_orig-1,kk);
873 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
876 amrex::Real dz =
z - orig_z_stag;
877 if ((dz < 0.0) && (dz > dzhi1) && (kk != khi0)) {
883 amrex::Real y0 = orig_data(i,j,khi0,src_comp);
884 amrex::Real y1 = orig_data(i,j,khi1,src_comp);
885 return ( y0-(y1-y0)/(z1-z0)*(z0-
z) );
888 }
else if (khi0 == -1) {
892 z0 = orig_z(i,j,khi0);
894 else if (stag ==
'X') {
896 z0 = orig_z(i,j,khi0);
898 else if (i == imax_orig) {
899 z0 = orig_z(imax_orig-1,j,khi0);
902 z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
905 else if (stag ==
'Y') {
907 z0 = orig_z(i,j,khi0);
909 else if (j == jmax_orig) {
910 z0 = orig_z(i,jmax_orig-1,khi0);
913 z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
916 amrex::Real y0 = orig_data(i,j,khi0,src_comp);
917 amrex::Real y1 = orig_data(i,j,khi1,src_comp);
918 return ( y0+(y1-y0)/(z1-z0)*(
z-z0) );
921 amrex::Real y0 = orig_data(i,j,klow,src_comp);
922 amrex::Real y1 = orig_data(i,j,khi0,src_comp);
923 return ( y0+(y1-y0)/(z1-z0)*(
z-z0) );
934 const amrex::Array4<amrex::Real const>& rhum,
935 const amrex::Array4<amrex::Real const>& temp,
936 const amrex::Array4<amrex::Real const>&
pres,
937 const amrex::Array4<amrex::Real>& mxrat)
939 amrex::Real qv_max_p_safe = 10000.0;
940 amrex::Real qv_max_flag = 1.0e-5;
941 amrex::Real qv_max_value = 3.0e-6;
942 amrex::Real qv_min_p_safe = 110000.0;
943 amrex::Real qv_min_flag = 1.0e-6;
944 amrex::Real qv_min_value = 1.0e-6;
945 amrex::Real eps = 0.622;
946 amrex::Real
svp1 = 0.6112;
947 amrex::Real
svp2 = 17.67;
948 amrex::Real
svp3 = 29.65;
949 amrex::Real
svpt0 = 273.15;
951 if (temp(i,j,k) != 0.0) {
952 amrex::Real es=0.01*rhum(i,j,k)*
svp1*10.0*exp(
svp2*(temp(i,j,k)-
svpt0)/(temp(i,j,k)-
svp3));
953 if (es >=
pres(i,j,k)/100.0) {
955 mxrat(i,j,k) = std::pow(10.0, -6);
958 mxrat(i,j,k) = amrex::max(eps*es/(
pres(i,j,k)/100.0-es), 1.0e-6);
964 mxrat(i,j,k) = 1.0e-6;
972 if (
pres(i,j,k) < qv_max_p_safe) {
973 if (mxrat(i,j,k) > qv_max_flag) {
974 mxrat(i,j,k) = qv_max_value;
977 if (
pres(i,j,k) < qv_min_p_safe) {
978 if (mxrat(i,j,k) < qv_min_flag) {
979 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:459
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:116
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:773
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:450
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:931
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:146
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 init_terrain_from_metgrid(amrex::FArrayBox &z_phys_nd_fab, const amrex::Vector< amrex::FArrayBox > &NC_hgt_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_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)
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
real(c_double), parameter svp1
Definition: ERF_module_model_constants.F90:78
real(c_double), parameter svp3
Definition: ERF_module_model_constants.F90:80
real(c_double), parameter svp2
Definition: ERF_module_model_constants.F90:79
real(c_double), parameter svpt0
Definition: ERF_module_model_constants.F90:81