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& z_phys_cc_fab,
108 const amrex::Vector<amrex::FArrayBox>& NC_psfc_fab);
124 amrex::Real Px = 0.0;
125 for (
int i=0; i <= order; i++) {
128 for (
int k=0; k <= order; k++) {
129 if (k == i)
continue;
144 const bool& exp_interp,
150 amrex::Real* orig_x_z,
151 amrex::Real* orig_x_p,
153 amrex::Real* new_x_z,
154 amrex::Real* new_x_p,
158 amrex::Real CRC_const1 = 11880.516;
159 amrex::Real CRC_const2 = 0.1902632;
160 amrex::Real CRC_const3 = 0.0065;
162 amrex::ignore_unused(i,j);
163 #ifndef AMREX_USE_GPU
167 for (
int new_k=0; new_k < new_n; new_k++) {
168 #ifndef AMREX_USE_GPU
169 if (debug) amrex::Print() <<
"new_k=" << new_k;
173 bool extrapolating =
true;
175 for (
int ko=0; ko < orig_n-1; ko++) {
176 amrex::Real a = new_x_z[new_k]-orig_x_z[ko];
177 amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1];
181 extrapolating =
false;
187 if (var_type ==
'T') {
190 amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[0];
191 amrex::Real avg_of_extrap_p = 0.5*(new_x_p[new_k]+orig_x_p[0]);
192 amrex::Real temp_extrap_starting_point = orig_y[0]*std::pow(orig_x_p[0]/100000.0,
R_d/
Cp_d);
193 amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/100.0, CRC_const2-1.0);
194 amrex::Real dZ = dZdP*(depth_of_extrap_in_p/100.0);
195 amrex::Real dT = dZ*CRC_const3;
196 new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(100000.0/new_x_p[new_k],
R_d/
Cp_d);
200 new_y[new_k] = orig_y[0];
206 if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_n-1)) {
208 int ksta = kl-(((order+1)/2)-1);
209 int kend = ksta+order;
210 #ifndef AMREX_USE_GPU
211 int ksize = kend-ksta;
212 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;
215 amrex::GpuArray<amrex::Real,9> orig_x_sub;
216 amrex::Real* orig_x_sub_p = orig_x_sub.data();
218 new_x = new_x_p[new_k];
219 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
221 new_x = new_x_z[new_k];
222 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
224 amrex::GpuArray<amrex::Real,9> orig_y_sub;
225 amrex::Real* orig_y_sub_p = orig_y_sub.data();
226 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
227 #ifndef AMREX_USE_GPU
229 amrex::Print() <<
" orig_x_sub_p = [";
230 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
231 amrex::Print() <<
"]" << std::endl;
232 amrex::Print() <<
" orig_y_sub_p = [";
233 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
234 amrex::Print() <<
"]" << std::endl;
237 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
243 #ifndef AMREX_USE_GPU
244 int ksize = kend-ksta+1;
245 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;
248 amrex::GpuArray<amrex::Real,2> orig_x_sub;
249 amrex::GpuArray<amrex::Real,2> orig_y_sub;
250 amrex::Real* orig_x_sub_p = orig_x_sub.data();
251 amrex::Real* orig_y_sub_p = orig_y_sub.data();
253 new_x = new_x_p[new_k];
254 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
256 new_x = new_x_z[new_k];
257 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
259 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
260 #ifndef AMREX_USE_GPU
262 amrex::Print() <<
" orig_x_sub_p = [";
263 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
264 amrex::Print() <<
"]" << std::endl;
265 amrex::Print() <<
" orig_y_sub_p = [";
266 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
267 amrex::Print() <<
"]" << std::endl;
272 }
else if (order%2 == 0) {
273 if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_n-1)) {
276 amrex::Real new_y_l, new_y_r;
278 int ksta = kl-(order/2-1);
279 int kend = ksta+order;
280 #ifndef AMREX_USE_GPU
281 int ksize = kend-ksta;
282 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;
285 amrex::GpuArray<amrex::Real,10> orig_x_sub;
286 amrex::GpuArray<amrex::Real,10> orig_y_sub;
287 amrex::Real* orig_x_sub_p = orig_x_sub.data();
288 amrex::Real* orig_y_sub_p = orig_y_sub.data();
290 new_x = new_x_p[new_k];
291 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
293 new_x = new_x_z[new_k];
294 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
296 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
297 #ifndef AMREX_USE_GPU
299 amrex::Print() <<
" orig_x_sub_p = [";
300 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
301 amrex::Print() <<
"]" << std::endl;
302 amrex::Print() <<
" orig_y_sub_p = [";
303 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
304 amrex::Print() <<
"]" << std::endl;
310 int ksta = kl-order/2;
311 int kend = ksta+order;
312 #ifndef AMREX_USE_GPU
313 int ksize = kend-ksta;
314 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;
317 amrex::GpuArray<amrex::Real,10> orig_x_sub;
318 amrex::GpuArray<amrex::Real,10> orig_y_sub;
319 amrex::Real* orig_x_sub_p = orig_x_sub.data();
320 amrex::Real* orig_y_sub_p = orig_y_sub.data();
322 new_x = new_x_p[new_k];
323 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
325 new_x = new_x_z[new_k];
326 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
328 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
329 #ifndef AMREX_USE_GPU
331 amrex::Print() <<
" orig_x_sub_p = [";
332 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
333 amrex::Print() <<
"]" << std::endl;
334 amrex::Print() <<
" orig_y_sub_p = [";
335 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
336 amrex::Print() <<
"]" << std::endl;
341 new_y[new_k] = 0.5*(new_y_l+new_y_r);
342 }
else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_n-1)) {
343 int ksta = kl-(order/2-1);
344 int kend = ksta+order;
345 #ifndef AMREX_USE_GPU
346 int ksize = kend-ksta;
347 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;
350 amrex::GpuArray<amrex::Real,10> orig_x_sub;
351 amrex::GpuArray<amrex::Real,10> orig_y_sub;
352 amrex::Real* orig_x_sub_p = orig_x_sub.data();
353 amrex::Real* orig_y_sub_p = orig_y_sub.data();
355 new_x = new_x_p[new_k];
356 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
358 new_x = new_x_z[new_k];
359 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
361 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
362 #ifndef AMREX_USE_GPU
364 amrex::Print() <<
" orig_x_sub_p = [";
365 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
366 amrex::Print() <<
"]" << std::endl;
367 amrex::Print() <<
" orig_y_sub_p = [";
368 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
369 amrex::Print() <<
"]" << std::endl;
372 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
373 }
else if ((kl-order/2 >= 0) && (kr+order/2-1 <= orig_n-1)) {
374 int ksta = kl-order/2;
375 int kend = ksta+order;
376 #ifndef AMREX_USE_GPU
377 int ksize = kend-ksta;
378 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;
381 amrex::GpuArray<amrex::Real,10> orig_x_sub;
382 amrex::GpuArray<amrex::Real,10> orig_y_sub;
383 amrex::Real* orig_x_sub_p = orig_x_sub.data();
384 amrex::Real* orig_y_sub_p = orig_y_sub.data();
386 new_x = new_x_p[new_k];
387 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
389 new_x = new_x_z[new_k];
390 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
392 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
393 #ifndef AMREX_USE_GPU
395 amrex::Print() <<
" orig_x_sub_p = [";
396 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
397 amrex::Print() <<
"]" << std::endl;
398 amrex::Print() <<
" orig_y_sub_p = [";
399 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
400 amrex::Print() <<
"]" << std::endl;
403 lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
408 #ifndef AMREX_USE_GPU
409 int ksize = kend-ksta+1;
410 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;
413 amrex::GpuArray<amrex::Real,2> orig_x_sub;
414 amrex::GpuArray<amrex::Real,2> orig_y_sub;
415 amrex::Real* orig_x_sub_p = orig_x_sub.data();
416 amrex::Real* orig_y_sub_p = orig_y_sub.data();
418 new_x = new_x_p[new_k];
419 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
421 new_x = new_x_z[new_k];
422 for (
int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
424 for (
int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
425 #ifndef AMREX_USE_GPU
427 amrex::Print() <<
" orig_x_sub_p = [";
428 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_x_sub_p[k];
429 amrex::Print() <<
"]" << std::endl;
430 amrex::Print() <<
" orig_y_sub_p = [";
431 for (
int k=0; k < ksize; k++) amrex::Print() <<
" " << orig_y_sub_p[k];
432 amrex::Print() <<
"]" << std::endl;
438 #ifndef AMREX_USE_GPU
439 if (debug) amrex::Print() <<
" new_y[" << new_k <<
"]=" << new_y[new_k] << std::endl;
457 const bool& metgrid_use_sfc,
458 const bool& exp_interp,
459 const bool& metgrid_retain_sfc,
460 const amrex::Real& metgrid_proximity,
461 const int& metgrid_order,
462 const int& metgrid_force_sfc_k,
469 const amrex::Array4<amrex::Real const>& orig_z_full,
470 const amrex::Array4<amrex::Real const>& orig_data,
471 const amrex::Array4<amrex::Real const>& new_z_full,
472 const amrex::Array4<amrex::Real>& new_data_full,
473 const bool& update_bc_data,
474 const amrex::Array4<amrex::Real>& bc_data_xlo,
475 const amrex::Array4<amrex::Real>& bc_data_xhi,
476 const amrex::Array4<amrex::Real>& bc_data_ylo,
477 const amrex::Array4<amrex::Real>& bc_data_yhi,
478 const amrex::Box& bx_xlo,
479 const amrex::Box& bx_xhi,
480 const amrex::Box& bx_ylo,
481 const amrex::Box& bx_yhi,
482 const amrex::Array4<const int>& mask)
488 int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
489 int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
490 int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
491 int kmax_new = amrex::ubound(amrex::Box(new_z_full)).z;
493 AMREX_ASSERT(kmax_orig < 256);
494 AMREX_ASSERT(kmax_new < 256);
496 amrex::GpuArray<amrex::Real,256> new_z;
497 amrex::GpuArray<amrex::Real,256> new_p;
498 amrex::GpuArray<amrex::Real,256> new_data;
499 amrex::Real* new_z_p = new_z.data();
500 amrex::Real* new_p_p = new_p.data();
501 amrex::Real* new_data_p = new_data.data();
502 for (
int k=0; k < kmax_new; k++) {
504 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));
505 }
else if (stag ==
'Y') {
506 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));
507 }
else if (stag ==
'M') {
508 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 )+
509 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));
514 amrex::GpuArray<amrex::Real,256> orig_z;
515 amrex::Real* orig_z_p = orig_z.data();
516 for (
int k=0; k < kmax_orig; k++) {
518 orig_z_p[k] = orig_z_full(i,j,k);
519 }
else if (stag ==
'X') {
521 orig_z_p[k] = orig_z_full(i,j,k);
522 }
else if (i == imax_orig) {
523 orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
525 orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
527 }
else if (stag ==
'Y') {
529 orig_z_p[k] = orig_z_full(i,j,k);
530 }
else if (j == jmax_orig) {
531 orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
533 orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
539 bool flip_data_required =
false;
540 if (orig_z[1] > orig_z[kmax_orig-1]) flip_data_required =
true;
541 if (flip_data_required) amrex::Abort(
"metgrid initialization flip_data_required. Not yet implemented.");
547 for (
int k=1; k < kmax_orig; k++) {
548 if (orig_z_p[k] > orig_z_p[0]) {
555 amrex::GpuArray<amrex::Real,256> ordered_z;
556 amrex::GpuArray<amrex::Real,256> ordered_data;
557 amrex::Real* ordered_z_p = ordered_z.data();
558 amrex::Real* ordered_data_p = ordered_data.data();
559 if (k_above_sfc > 1) {
565 for (
int k=1; k < k_above_sfc; k++) {
566 ordered_z_p[count] = orig_z_p[k];
567 ordered_data_p[count] = orig_data(i,j,k);
581 if (Pl-Pu < metgrid_proximity) {
586 ordered_z_p[count] = orig_z_p[0];
587 ordered_data_p[count] = orig_data(i,j,0);
592 int knext = k_above_sfc;
598 if (metgrid_force_sfc_k > 0) {
599 for (
int k=k_above_sfc; k < kmax_orig; k++) {
600 if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
611 if (Pl-Pu < metgrid_proximity) {
616 for (
int k=knext; k < kmax_orig; k++) {
617 ordered_z_p[count] = orig_z_p[k];
618 ordered_data_p[count] = orig_data(i,j,k);
627 ordered_z_p[0] = orig_z[0];
628 ordered_data_p[0] = orig_data(i,j,0);
634 if (metgrid_force_sfc_k > 0) {
635 for (
int k=knext; k < kmax_orig; k++) {
636 if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k]) {
645 for (
int k=knext; k < kmax_orig; k++) {
649 if (Pl-Pu < metgrid_proximity) {
652 ordered_z_p[count] = orig_z_p[k];
653 ordered_data_p[count] = orig_data(i,j,k);
659 int ksta(0), kend(0);
660 if (metgrid_use_below_sfc && metgrid_use_sfc) {
664 }
else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
667 for (
int k=0; k < kmax_orig; k++) {
668 if (ordered_z_p[k] == orig_z_p[0]) {
673 for (
int k=ksfc; k < kmax_orig-1; k++) {
674 ordered_z_p[k] = ordered_z_p[k+1];
675 ordered_data_p[k] = ordered_data_p[k+1];
679 }
else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
682 for (
int k=0; k < kend_order; k++) {
683 if (ordered_z_p[k] == orig_z_p[0]) {
692 amrex::Abort(
"metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
708 amrex::GpuArray<amrex::Real,256> ordered_p;
709 amrex::Real* ordered_p_p = ordered_p.data();
710 for (
int k=0; k < kend_order; k++) {
716 amrex::GpuArray<amrex::Real,256> final_z;
717 amrex::GpuArray<amrex::Real,256> final_p;
718 amrex::GpuArray<amrex::Real,256> final_data;
719 amrex::Real* final_z_p = final_z.data();
720 amrex::Real* final_p_p = final_p.data();
721 amrex::Real* final_data_p = final_data.data();
722 final_z_p[0] = ordered_z[ksta];
723 final_p_p[0] = ordered_p[ksta];
724 final_data_p[0] = ordered_data[ksta];
725 for (
int k=ksta+1; k <= kend; k++) {
726 if ((final_p_p[new_n]-ordered_p_p[k]) < metgrid_proximity) {
730 final_z_p[new_n] = ordered_z_p[k];
731 final_p_p[new_n] = ordered_p_p[k];
732 final_data_p[new_n] = ordered_data_p[k];
754 if (metgrid_retain_sfc) new_data[0] = ordered_data[0];
757 for (
int k=0; k < kmax_new; k++) {
758 if (mask(i,j,k) && update_bc_data && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = new_data[k];
759 if (mask(i,j,k) && update_bc_data && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = new_data[k];
760 if (mask(i,j,k) && update_bc_data && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = new_data[k];
761 if (mask(i,j,k) && update_bc_data && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = new_data[k];
762 if (itime == 0) new_data_full(i,j,k,src_comp) = new_data[k];
775 const amrex::Array4<amrex::Real const>& orig_z,
776 const amrex::Array4<amrex::Real const>& orig_data,
777 const amrex::Array4<amrex::Real const>& new_z)
780 int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
781 int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
782 int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
786 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));
788 else if (stag ==
'Y') {
789 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));
791 else if (stag ==
'M') {
792 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 )+
793 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));
799 amrex::Real dzlow = 1.0e12;
800 amrex::Real dzhi0 = -1.0e12;
801 for (
int kk = 0; kk < kmax_orig; kk++) {
802 amrex::Real orig_z_stag = 0.0;
804 orig_z_stag = orig_z(i,j,kk);
808 orig_z_stag = orig_z(i,j,kk);
810 else if (i == imax_orig) {
811 orig_z_stag = orig_z(imax_orig-1,j,kk);
814 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
817 else if (stag ==
'Y') {
819 orig_z_stag = orig_z(i,j,kk);
821 else if (j == jmax_orig) {
822 orig_z_stag = orig_z(i,jmax_orig-1,kk);
825 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
829 amrex::Real dz =
z - orig_z_stag;
830 if ((dz < 0.0) && (dz > dzhi0)) {
835 if ((dz >= 0.0) && (dz < dzlow)) {
845 amrex::Real dzhi1 = -1.0e12;
846 for (
int kk = 0; kk < kmax_orig; kk++) {
847 amrex::Real orig_z_stag = 0.0;
849 orig_z_stag = orig_z(i,j,kk);
851 else if (stag ==
'X') {
853 orig_z_stag = orig_z(i,j,kk);
855 else if (i == imax_orig) {
856 orig_z_stag = orig_z(imax_orig-1,j,kk);
859 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
862 else if (stag ==
'Y') {
864 orig_z_stag = orig_z(i,j,kk);
866 else if (j == jmax_orig) {
867 orig_z_stag = orig_z(i,jmax_orig-1,kk);
870 orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
873 amrex::Real dz =
z - orig_z_stag;
874 if ((dz < 0.0) && (dz > dzhi1) && (kk != khi0)) {
880 amrex::Real y0 = orig_data(i,j,khi0,src_comp);
881 amrex::Real y1 = orig_data(i,j,khi1,src_comp);
882 return ( y0-(y1-y0)/(z1-z0)*(z0-
z) );
885 }
else if (khi0 == -1) {
889 z0 = orig_z(i,j,khi0);
891 else if (stag ==
'X') {
893 z0 = orig_z(i,j,khi0);
895 else if (i == imax_orig) {
896 z0 = orig_z(imax_orig-1,j,khi0);
899 z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
902 else if (stag ==
'Y') {
904 z0 = orig_z(i,j,khi0);
906 else if (j == jmax_orig) {
907 z0 = orig_z(i,jmax_orig-1,khi0);
910 z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
913 amrex::Real y0 = orig_data(i,j,khi0,src_comp);
914 amrex::Real y1 = orig_data(i,j,khi1,src_comp);
915 return ( y0+(y1-y0)/(z1-z0)*(
z-z0) );
918 amrex::Real y0 = orig_data(i,j,klow,src_comp);
919 amrex::Real y1 = orig_data(i,j,khi0,src_comp);
920 return ( y0+(y1-y0)/(z1-z0)*(
z-z0) );
931 const amrex::Array4<amrex::Real const>& rhum,
932 const amrex::Array4<amrex::Real const>& temp,
933 const amrex::Array4<amrex::Real const>&
pres,
934 const amrex::Array4<amrex::Real>& mxrat)
936 amrex::Real qv_max_p_safe = 10000.0;
937 amrex::Real qv_max_flag = 1.0e-5;
938 amrex::Real qv_max_value = 3.0e-6;
939 amrex::Real qv_min_p_safe = 110000.0;
940 amrex::Real qv_min_flag = 1.0e-6;
941 amrex::Real qv_min_value = 1.0e-6;
942 amrex::Real eps = 0.622;
943 amrex::Real svp1 = 0.6112;
944 amrex::Real svp2 = 17.67;
945 amrex::Real svp3 = 29.65;
946 amrex::Real svpt0 = 273.15;
948 if (temp(i,j,k) != 0.0) {
949 amrex::Real es=0.01*rhum(i,j,k)*svp1*10.0*exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3));
950 if (es >=
pres(i,j,k)/100.0) {
952 mxrat(i,j,k) = std::pow(10.0, -6);
955 mxrat(i,j,k) = amrex::max(eps*es/(
pres(i,j,k)/100.0-es), 1.0e-6);
961 mxrat(i,j,k) = 1.0e-6;
969 if (
pres(i,j,k) < qv_max_p_safe) {
970 if (mxrat(i,j,k) > qv_max_flag) {
971 mxrat(i,j,k) = qv_max_value;
974 if (
pres(i,j,k) < qv_min_p_safe) {
975 if (mxrat(i,j,k) < qv_min_flag) {
976 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:456
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 &z_phys_cc_fab, const amrex::Vector< amrex::FArrayBox > &NC_psfc_fab)
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:113
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:770
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:447
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:928
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:143
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:33