ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_MetgridUtils.H File Reference
#include <ERF.H>
#include <ERF_EOS.H>
#include <ERF_Utils.H>
#include <ERF_ProbCommon.H>
#include <ERF_HSEUtils.H>
Include dependency graph for ERF_MetgridUtils.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

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_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)
 
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)
 
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)
 
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)
 
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)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal (const amrex::Real &z, amrex::Real &p)
 
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)
 
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)
 
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)
 

Function Documentation

◆ calc_p_isothermal()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal ( const amrex::Real &  z,
amrex::Real &  p 
)
450 {
451  p = p_0*exp(-CONST_GRAV*z/(290.0*R_d));
452 }
p_0
constexpr amrex::Real p_0
Definition: ERF_Constants.H:18
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10

Referenced by interpolate_column_metgrid().

Here is the caller graph for this function:

◆ 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 &  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 
)

◆ 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,
const amrex::Vector< amrex::FArrayBox > &  NC_MSFU_fab,
const amrex::Vector< amrex::FArrayBox > &  NC_MSFV_fab,
const amrex::Vector< amrex::FArrayBox > &  NC_MSFM_fab 
)

◆ init_state_from_metgrid()

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 
)

◆ init_terrain_from_metgrid()

void init_terrain_from_metgrid ( amrex::FArrayBox &  z_phys_nd_fab,
const amrex::Vector< amrex::FArrayBox > &  NC_hgt_fab 
)

◆ 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 &  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 
)
484 {
485  // Here we closely follow WRF's vert_interp from
486  // dyn_em/module_initialize_real.F, although changes have been
487  // made to accommodate interpolation relative to height instead of
488  // pressure.
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;
493 
494  AMREX_ASSERT(kmax_orig < 256);
495  AMREX_ASSERT(kmax_new < 256);
496 
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++) {
504  if (stag == 'X') {
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));
511  }
512  calc_p_isothermal(new_z_p[k], new_p_p[k]);
513  }
514 
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++) {
518  if (stag == 'M') {
519  orig_z_p[k] = orig_z_full(i,j,k);
520  } else if (stag == 'X') {
521  if (i == 0) {
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);
525  } else {
526  orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
527  }
528  } else if (stag == 'Y') {
529  if (j == 0) {
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);
533  } else {
534  orig_z_p[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
535  }
536  }
537  }
538 
539  // Check if the data is top-down instead of bottom-up.
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.");
543 
544  // Search for the first level above the surface in the origin data.
545  // This is needed since the origin model topography will be
546  // different than the topography processed by WPS.
547  int k_above_sfc = 0;
548  for (int k=1; k < kmax_orig; k++) {
549  if (orig_z_p[k] > orig_z_p[0]) {
550  k_above_sfc = k;
551  break;
552  }
553  }
554 
555  int kend_order;
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) {
561  // The levels are not monotonically increasing in height, so
562  // we sort and then make "artistic" quality control choices.
563  int count = 0;
564 
565  // Insert levels that are below the surface.
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);
569  count++;
570  }
571 
572  // Check if the level that is nearest to and below the surface
573  // is "too close". If so, we'll ignore the upper level and keep
574  // the lower. Origin data is likely to be on pressure levels
575  // with higher spatial resolution near-surface, which supports
576  // the choice of eliminating levels that are "too close" in
577  // pressure-space. For simplicity, calculate delta P assuming a
578  // baroclinic atmosphere.
579  amrex::Real Pu, Pl;
580  calc_p_isothermal(orig_z_p[0], Pu);
581  calc_p_isothermal(ordered_z_p[count-1], Pl);
582  if (Pl-Pu < metgrid_proximity) {
583  count--;
584  }
585 
586  // Insert the surface level.
587  ordered_z_p[count] = orig_z_p[0];
588  ordered_data_p[count] = orig_data(i,j,0);
589  count++;
590 
591  // Quoting WRF's comments, the next level to use is at,
592  // "... ta da, the first level above the surface. I know, wow."
593  int knext = k_above_sfc;
594  // Conditionally more strongly use the surface data by removing
595  // levels between the surface and the height corresponding to a
596  // set number of ERF levels from the surface. This forces the
597  // interpolator to use the surface data up through a number of
598  // ERF levels from the surface.
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]) {
602  knext = k;
603  break;
604  }
605  }
606  }
607 
608  // Check if the level that is nearest to and above the surface
609  // is "too close". If so, we'll ignore that level.
610  calc_p_isothermal(orig_z_p[knext], Pu);
611  calc_p_isothermal(ordered_z_p[count-1], Pl);
612  if (Pl-Pu < metgrid_proximity) {
613  knext++;
614  }
615 
616  // Insert levels that are above the surface.
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);
620  count++;
621  }
622 
623  kend_order = count;
624  } else {
625  // The surface is the lowest level in the origin data.
626 
627  // Insert the surface.
628  ordered_z_p[0] = orig_z[0];
629  ordered_data_p[0] = orig_data(i,j,0);
630 
631  // Similar to above, conditionally more strongly use the
632  // surface data.
633  int count = 1;
634  int knext = count;
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]) {
638  knext = k;
639  break;
640  }
641  }
642  }
643 
644  // Insert the remaining levels, again ignoring levels that are
645  // "too close" to the prior valid level.
646  for (int k=knext; k < kmax_orig; k++) {
647  amrex::Real Pu, Pl;
648  calc_p_isothermal(orig_z_p[k], Pu);
649  calc_p_isothermal(ordered_z_p[count-1], Pl);
650  if (Pl-Pu < metgrid_proximity) {
651  continue;
652  }
653  ordered_z_p[count] = orig_z_p[k];
654  ordered_data_p[count] = orig_data(i,j,k);
655  count++;
656  }
657  kend_order = count;
658  }
659 
660  int ksta(0), kend(0);
661  if (metgrid_use_below_sfc && metgrid_use_sfc) {
662  // Use all levels.
663  ksta = 0;
664  kend = kend_order-1;
665  } else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
666  // Use all levels except for the surface.
667  int ksfc = 0;
668  for (int k=0; k < kmax_orig; k++) {
669  if (ordered_z_p[k] == orig_z_p[0]) {
670  ksfc = k;
671  break;
672  }
673  }
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];
677  }
678  ksta = 0;
679  kend = kend_order-2;
680  } else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
681  // Use all levels above and including the surface.
682  int ksfc = 0;
683  for (int k=0; k < kend_order; k++) {
684  if (ordered_z_p[k] == orig_z_p[0]) {
685  ksfc = k;
686  break;
687  }
688  }
689  ksta = ksfc;
690  kend = kend_order-1;
691  } else {
692  // We shouldn't be in here!
693  amrex::Abort("metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
694  }
695 
696  // Insert the level of maximum winds.
697 // amrex::Real maxw_above_this_level = 30000.0;
698 // amrex::Real maxw_horiz_pres_diff = 5000.0;
699 // if ((flag_maxw == 1) && (use_maxw)) {
700 // amrex::Abort("metgrid initialization, use_maxw not yet implemented");
701 // }
702 
703  // Insert the level of the tropopause.
704 // amrex::Real trop_horiz_pres_diff = 5000.0;
705 // if ((flag_trop == 1) && (use_trop)) {
706 // amrex::Abort("metgrid initialization, use_trop not yet implemented");
707 // }
708 
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++) {
712  calc_p_isothermal(ordered_z_p[k], ordered_p_p[k]);
713  }
714 
715  int new_n = 0;
716  int zap_final = 0;
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) {
728  zap_final++;
729  } else {
730  new_n++;
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];
734  }
735  }
736  kend -= zap_final;
737 
738  // Call the interpolator.
739  lagrange_setup(var_type,
740  exp_interp,
741  kend-ksta,
742  kmax_new,
743  metgrid_order,
744  i,
745  j,
746  final_z_p,
747  final_p_p,
748  final_data_p,
749  new_z_p,
750  new_p_p,
751  new_data_p);
752 
753  // Optionally replace the lowest level of data with the surface
754  // field from the origin data.
755  if (metgrid_retain_sfc) new_data[0] = ordered_data[0];
756 
757  // Save the interpolated data.
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];
764  }
765 
766 }
calc_p_isothermal
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:448
lagrange_setup
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
Here is the call graph for this function:

◆ 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 
)
779 {
780  // This subroutine is a bit ham-handed and can be cleaned up later.
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;
784 
785  amrex::Real z;
786  if (stag == 'X') {
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));
788  }
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));
791  }
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));
795  }
796 
797  amrex::Real z0, z1;
798  int klow = -1;
799  int khi0 = -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;
804  if (stag == 'M') {
805  orig_z_stag = orig_z(i,j,kk);
806  }
807  if (stag == 'X') {
808  if (i == 0) {
809  orig_z_stag = orig_z(i,j,kk);
810  }
811  else if (i == imax_orig) {
812  orig_z_stag = orig_z(imax_orig-1,j,kk);
813  }
814  else {
815  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
816  }
817  }
818  else if (stag == 'Y') {
819  if (j == 0) {
820  orig_z_stag = orig_z(i,j,kk);
821  }
822  else if (j == jmax_orig) {
823  orig_z_stag = orig_z(i,jmax_orig-1,kk);
824  }
825  else {
826  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
827  }
828  }
829 
830  amrex::Real dz = z - orig_z_stag;
831  if ((dz < 0.0) && (dz > dzhi0)) {
832  dzhi0 = dz;
833  khi0 = kk;
834  z1 = orig_z_stag;
835  }
836  if ((dz >= 0.0) && (dz < dzlow)) {
837  dzlow = dz;
838  klow = kk;
839  z0 = orig_z_stag;
840  }
841  } // kk
842 
843  // extrapolate below the bottom surface
844  if (klow == -1) {
845  int khi1 = -1;
846  amrex::Real dzhi1 = -1.0e12;
847  for (int kk = 0; kk < kmax_orig; kk++) {
848  amrex::Real orig_z_stag = 0.0;
849  if (stag == 'M') {
850  orig_z_stag = orig_z(i,j,kk);
851  }
852  else if (stag == 'X') {
853  if (i == 0) {
854  orig_z_stag = orig_z(i,j,kk);
855  }
856  else if (i == imax_orig) {
857  orig_z_stag = orig_z(imax_orig-1,j,kk);
858  }
859  else {
860  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
861  }
862  }
863  else if (stag == 'Y') {
864  if (j == 0) {
865  orig_z_stag = orig_z(i,j,kk);
866  }
867  else if (j == jmax_orig) {
868  orig_z_stag = orig_z(i,jmax_orig-1,kk);
869  }
870  else {
871  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
872  }
873  }
874  amrex::Real dz = z - orig_z_stag;
875  if ((dz < 0.0) && (dz > dzhi1) && (kk != khi0)) {
876  dzhi1 = dz;
877  khi1 = kk;
878  z1 = orig_z_stag;
879  }
880  }
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) );
884 
885  // Extrapolate above the top surface
886  } else if (khi0 == -1) {
887  khi0 = klow - 1;
888  int khi1 = klow;
889  if (stag == 'M') {
890  z0 = orig_z(i,j,khi0);
891  }
892  else if (stag == 'X') {
893  if (i == 0) {
894  z0 = orig_z(i,j,khi0);
895  }
896  else if (i == imax_orig) {
897  z0 = orig_z(imax_orig-1,j,khi0);
898  }
899  else {
900  z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
901  }
902  }
903  else if (stag == 'Y') {
904  if (j == 0) {
905  z0 = orig_z(i,j,khi0);
906  }
907  else if (j == jmax_orig) {
908  z0 = orig_z(i,jmax_orig-1,khi0);
909  }
910  else {
911  z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
912  }
913  }
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) );
917  } else {
918  // interpolate
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) );
922 
923  }
924 }
Coord::z
@ z

◆ 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 
)
119 {
120  // Interpolation using Lagrange polynomials.
121  // P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
122  // where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
123  // ---------------------------------------------
124  // (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
125  amrex::Real Px = 0.0;
126  for (int i=0; i <= order; i++) {
127  amrex::Real n = 1.0;
128  amrex::Real d = 1.0;
129  for (int k=0; k <= order; k++) {
130  if (k == i) continue;
131  n *= new_x-x[k];
132  d *= x[i]-x[k];
133  }
134  if (d != 0.0) {
135  Px += y[i]*n/d;
136  }
137  }
138  new_y = Px;
139 }
Coord::y
@ y
Coord::x
@ x

Referenced by lagrange_setup().

Here is the caller graph for this function:

◆ lagrange_setup()

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 
)
157 {
158 
159  amrex::Real CRC_const1 = 11880.516; // m
160  amrex::Real CRC_const2 = 0.1902632;
161  amrex::Real CRC_const3 = 0.0065; // K km-1
162 
163  amrex::ignore_unused(i,j);
164 #ifndef AMREX_USE_GPU
165  bool debug = false;
166 #endif
167 
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;
171 #endif
172  // Determine if interpolating or extrapolating.
173  // If interpolating, find bounding x values and store the indices.
174  bool extrapolating = true;
175  int kl, kr;
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];
179  if (a*b <= 0.0) {
180  kl = ko;
181  kr = ko+1;
182  extrapolating = false;
183  break;
184  }
185  }
186 
187  if (extrapolating) {
188  if (var_type == 'T') {
189  // Assume a standard atmosphere -6.5 K km-1 lapse rate.
190  // Comparable to the WRF default, t_extrap_type=2.
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);
198  } else {
199  // Use a constant value below ground.
200  // Comparable to the WRF default, extrap_type=2.
201  new_y[new_k] = orig_y[0];
202  }
203  continue;
204  }
205 
206  if (order%2 != 0) {
207  if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_n-1)) {
208  // Odd order interpolation.
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;
214 #endif
215  amrex::Real new_x;
216  amrex::GpuArray<amrex::Real,9> orig_x_sub;
217  amrex::Real* orig_x_sub_p = orig_x_sub.data();
218  if (exp_interp) {
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]; }
221  } else {
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]; }
224  }
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
229  if (debug) {
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;
236  }
237 #endif
238  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
239  } else {
240  // Interpolation stencil is too big due to proximity to the surface or model top.
241  // We have bounding points so resort to a simple linear interpolation.
242  int ksta = kl;
243  int kend = kr;
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;
247 #endif
248  amrex::Real new_x;
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();
253  if (exp_interp) {
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]; }
256  } else {
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]; }
259  }
260  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
261 #ifndef AMREX_USE_GPU
262  if (debug) {
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;
269  }
270 #endif
271  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
272  }
273  } else if (order%2 == 0) {
274  if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_n-1)) {
275  // Even order interpolation.
276  // Interpolate twice. Offset the range by 1 the 2nd time. Average the results.
277  amrex::Real new_y_l, new_y_r;
278  {
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;
284 #endif
285  amrex::Real new_x;
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();
290  if (exp_interp) {
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]; }
293  } else {
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]; }
296  }
297  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
298 #ifndef AMREX_USE_GPU
299  if (debug) {
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;
306  }
307 #endif
308  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_l);
309  }
310  {
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;
316 #endif
317  amrex::Real new_x;
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();
322  if (exp_interp) {
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]; }
325  } else {
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]; }
328  }
329  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
330 #ifndef AMREX_USE_GPU
331  if (debug) {
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;
338  }
339 #endif
340  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_r);
341  }
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;
349 #endif
350  amrex::Real new_x;
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();
355  if (exp_interp) {
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]; }
358  } else {
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]; }
361  }
362  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
363 #ifndef AMREX_USE_GPU
364  if (debug) {
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;
371  }
372 #endif
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;
380 #endif
381  amrex::Real new_x;
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();
386  if (exp_interp) {
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]; }
389  } else {
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]; }
392  }
393  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
394 #ifndef AMREX_USE_GPU
395  if (debug) {
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;
402  }
403 #endif
404  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
405  } else {
406  // Linear interpolation.
407  int ksta = kl;
408  int kend = kr;
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;
412 #endif
413  amrex::Real new_x;
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();
418  if (exp_interp) {
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]; }
421  } else {
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]; }
424  }
425  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
426 #ifndef AMREX_USE_GPU
427  if (debug) {
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;
434  }
435 #endif
436  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
437  }
438  }
439 #ifndef AMREX_USE_GPU
440  if (debug) amrex::Print() << " new_y[" << new_k << "]=" << new_y[new_k] << std::endl;
441 #endif
442  }
443 }
Cp_d
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
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:114

Referenced by interpolate_column_metgrid().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ read_from_metgrid()

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 
)

◆ 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,
const amrex::Array4< amrex::Real > &  mxrat 
)
936 {
937  amrex::Real qv_max_p_safe = 10000.0; // WRF default value
938  amrex::Real qv_max_flag = 1.0e-5; // WRF default value
939  amrex::Real qv_max_value = 3.0e-6; // WRF default value
940  amrex::Real qv_min_p_safe = 110000.0; // WRF default value
941  amrex::Real qv_min_flag = 1.0e-6; // WRF default value
942  amrex::Real qv_min_value = 1.0e-6; // WRF default value
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;
948  // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior)
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) {
952  // vapor pressure exceeds total pressure
953  mxrat(i,j,k) = std::pow(10.0, -6);
954  }
955  else {
956  mxrat(i,j,k) = amrex::max(eps*es/(pres(i,j,k)/100.0-es), 1.0e-6);
957  }
958  }
959  else {
960  // I don't know why there's a fringe case handled in WRF where T is absolute zero...
961  // Let's just deal with it here in case we also end up needing it.
962  mxrat(i,j,k) = 1.0e-6;
963  }
964  // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1.
965  // For pressures above a defined level, reasonable Qv values should be
966  // a certain value or smaller. If they are larger than this, the input data
967  // probably had "missing" RH, and we filled in some values. This is an
968  // attempt to catch those. Also, set the minimum value for the entire
969  // domain that is above the selected pressure level.
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;
973  }
974  }
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;
978  }
979  }
980 }
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:25