ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_MetgridUtils.H File Reference
#include <AMReX_Math.H>
#include <cmath>
#include "ERF_Constants.H"
#include <limits>
#include <ERF.H>
#include <ERF_EOS.H>
#include <ERF_SurfaceTemperature.H>
#include <ERF_Utils.H>
#include <ERF_ProbCommon.H>
#include <ERF_HSEUtils.H>
#include <ERF_BaseStateParams.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

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool metgrid_surface_pressure (const bool debug_psfc, const int flag_psfc, const amrex::Real file_psfc, const amrex::Real z_sfc, const amrex::Real P00, const amrex::Real T00, const amrex::Real TLP, amrex::Real &pressure)
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool metgrid_surface_theta (const amrex::Real temperature, const bool debug_psfc, const int flag_psfc, const amrex::Real file_psfc, const amrex::Real z_sfc, const amrex::Real P00, const amrex::Real T00, const amrex::Real TLP, const amrex::Real rdOcp, amrex::Real &theta)
 
void read_from_metgrid (int lev, int itime, const amrex::Box &domain, const std::string &fname, std::string &NC_dateTime, double &NC_epochTime, int &flag_psfc, int &flag_msf, int &flag_sst, int &flag_tsk, int &flag_lmask, int &NC_nx, int &NC_ny, amrex::Real &NC_dx, amrex::Real &NC_dy, amrex::FArrayBox &NC_xvel_fab, amrex::FArrayBox &NC_yvel_fab, amrex::FArrayBox &NC_temp_fab, amrex::FArrayBox &NC_rhum_fab, amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &NC_ght_fab, amrex::FArrayBox &NC_hgt_fab, amrex::FArrayBox &NC_psfc_fab, amrex::FArrayBox &NC_msfu_fab, amrex::FArrayBox &NC_msfv_fab, amrex::FArrayBox &NC_msfm_fab, amrex::FArrayBox &NC_sst_fab, amrex::FArrayBox &NC_tsk_fab, amrex::FArrayBox &NC_LAT_fab, amrex::FArrayBox &NC_LON_fab, amrex::IArrayBox &NC_lmask_iab, amrex::Geometry &geom)
 
void init_terrain_from_metgrid (const bool &avg_grid_faces_to_nodes, amrex::Geometry &geom, amrex::MultiFab *z_phys_nd, amrex::FArrayBox &NC_hgt_fab)
 
void init_state_from_metgrid (const int lev, const int itime, const bool use_moisture, const bool interp_theta, const bool metgrid_debug_quiescent, const bool metgrid_debug_isothermal, const bool metgrid_debug_dry, const bool metgrid_basic_linear, const bool metgrid_use_below_sfc, const bool metgrid_use_sfc, const bool metgrid_retain_sfc, const amrex::Real metgrid_proximity, const int metgrid_order, const int metgrid_metgrid_force_sfc_k, const amrex::Real l_rdOcp, amrex::Box &tbxc, amrex::Box &tbxu, amrex::Box &tbxv, amrex::Box &tbxw, amrex::FArrayBox &state_fab, amrex::FArrayBox &x_vel_fab, amrex::FArrayBox &y_vel_fab, amrex::FArrayBox &z_vel_fab, amrex::FArrayBox &z_phys_nd_fab, const amrex::FArrayBox &NC_ght_fab, const amrex::FArrayBox &NC_xvel_fab, const amrex::FArrayBox &NC_yvel_fab, const amrex::FArrayBox &NC_temp_fab, const amrex::FArrayBox &NC_rhum_fab, const amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &tmp_src_fab, amrex::FArrayBox &tmp_dst_fab, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xlo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xhi, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_ylo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_yhi, const amrex::Array4< const int > &mask_c_arr, const amrex::Array4< const int > &mask_u_arr, const amrex::Array4< const int > &mask_v_arr)
 
void init_msfs_from_metgrid (const bool metgrid_debug_msf, amrex::FArrayBox &msfu_fab, amrex::FArrayBox &msfv_fab, amrex::FArrayBox &msfm_fab, const int &flag_msf, amrex::FArrayBox &NC_MSFU_fab, amrex::FArrayBox &NC_MSFV_fab, amrex::FArrayBox &NC_MSFM_fab)
 
void init_base_state_from_metgrid (const bool use_moisture, const bool metgrid_debug_psfc, const amrex::Real l_rdOcp, const amrex::Box &domain, const amrex::Box &valid_bx, const int &flag_psfc, amrex::FArrayBox &state_fab, amrex::FArrayBox &r_hse_fab, amrex::FArrayBox &p_hse_fab, amrex::FArrayBox &pi_hse_fab, amrex::FArrayBox &th_hse_fab, amrex::FArrayBox &qv_hse_fab, amrex::FArrayBox &z_phys_nd_fab, amrex::FArrayBox &z_phys_cc_fab, const amrex::FArrayBox &NC_psfc_fab, int *base_state_column_failed, const BaseStateParams &bsp)
 
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_npts, const int &new_npts, const int &order, const int &i, const int &j, amrex::Real *orig_x_z, amrex::Real *orig_x_p, amrex::Real *orig_y, amrex::Real *new_x_z, amrex::Real *new_x_p, amrex::Real *new_y)
 
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 &kmax, const int &src_comp, const int &dst_comp, char var_type, char stag, const amrex::Array4< amrex::Real const > &orig_z_full, const amrex::Array4< amrex::Real const > &orig_data, const amrex::Array4< amrex::Real const > &new_z_full, const amrex::Array4< amrex::Real > &new_data_full)
 
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, int src_indx, 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 
)

Compute isothermal pressure at height z.

577 {
578  p = p_0*std::exp(-CONST_GRAV*z/(amrex::Real(290.0)*R_d));
579 }
p_0
constexpr amrex::Real p_0
Definition: ERF_Constants.H:69
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:72
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:49
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
WSM6Ind::p
@ p
Definition: ERF_WSM6.H:301

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 &  domain,
const amrex::Box &  valid_bx,
const int &  flag_psfc,
amrex::FArrayBox &  state_fab,
amrex::FArrayBox &  r_hse_fab,
amrex::FArrayBox &  p_hse_fab,
amrex::FArrayBox &  pi_hse_fab,
amrex::FArrayBox &  th_hse_fab,
amrex::FArrayBox &  qv_hse_fab,
amrex::FArrayBox &  z_phys_nd_fab,
amrex::FArrayBox &  z_phys_cc_fab,
const amrex::FArrayBox &  NC_psfc_fab,
int *  base_state_column_failed,
const BaseStateParams &  bsp 
)

Initialize hydrostatic base-state quantities from metgrid data.

◆ init_msfs_from_metgrid()

void init_msfs_from_metgrid ( const bool  metgrid_debug_msf,
amrex::FArrayBox &  msfu_fab,
amrex::FArrayBox &  msfv_fab,
amrex::FArrayBox &  msfm_fab,
const int &  flag_msf,
amrex::FArrayBox &  NC_MSFU_fab,
amrex::FArrayBox &  NC_MSFV_fab,
amrex::FArrayBox &  NC_MSFM_fab 
)

Initialize map-scale factors from metgrid fields or defaults.

◆ init_state_from_metgrid()

void init_state_from_metgrid ( const int  lev,
const int  itime,
const bool  use_moisture,
const bool  interp_theta,
const bool  metgrid_debug_quiescent,
const bool  metgrid_debug_isothermal,
const bool  metgrid_debug_dry,
const bool  metgrid_basic_linear,
const bool  metgrid_use_below_sfc,
const bool  metgrid_use_sfc,
const bool  metgrid_retain_sfc,
const amrex::Real  metgrid_proximity,
const int  metgrid_order,
const int  metgrid_metgrid_force_sfc_k,
const amrex::Real  l_rdOcp,
amrex::Box &  tbxc,
amrex::Box &  tbxu,
amrex::Box &  tbxv,
amrex::Box &  tbxw,
amrex::FArrayBox &  state_fab,
amrex::FArrayBox &  x_vel_fab,
amrex::FArrayBox &  y_vel_fab,
amrex::FArrayBox &  z_vel_fab,
amrex::FArrayBox &  z_phys_nd_fab,
const amrex::FArrayBox &  NC_ght_fab,
const amrex::FArrayBox &  NC_xvel_fab,
const amrex::FArrayBox &  NC_yvel_fab,
const amrex::FArrayBox &  NC_temp_fab,
const amrex::FArrayBox &  NC_rhum_fab,
const amrex::FArrayBox &  NC_pres_fab,
amrex::FArrayBox &  tmp_src_fab,
amrex::FArrayBox &  tmp_dst_fab,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_xlo,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_xhi,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_ylo,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_yhi,
const amrex::Array4< const int > &  mask_c_arr,
const amrex::Array4< const int > &  mask_u_arr,
const amrex::Array4< const int > &  mask_v_arr 
)

Initialize ERF state, velocity, and boundary data from metgrid fields.

◆ init_terrain_from_metgrid()

void init_terrain_from_metgrid ( const bool &  avg_grid_faces_to_nodes,
amrex::Geometry &  geom,
amrex::MultiFab *  z_phys_nd,
amrex::FArrayBox &  NC_hgt_fab 
)

Initialize nodal terrain heights from metgrid surface-height data.

◆ interpolate_column_metgrid()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void interpolate_column_metgrid ( const bool &  metgrid_use_below_sfc,
const bool &  metgrid_use_sfc,
const bool &  exp_interp,
const bool &  metgrid_retain_sfc,
const amrex::Real &  metgrid_proximity,
const int &  metgrid_order,
const int &  metgrid_force_sfc_k,
const int &  i,
const int &  j,
const int &  kmax,
const int &  src_comp,
const int &  dst_comp,
char  var_type,
char  stag,
const amrex::Array4< amrex::Real const > &  orig_z_full,
const amrex::Array4< amrex::Real const > &  orig_data,
const amrex::Array4< amrex::Real const > &  new_z_full,
const amrex::Array4< amrex::Real > &  new_data_full 
)

Interpolate one metgrid column onto ERF vertical levels with quality control.

605 {
606  // Here we closely follow WRF's vert_interp from
607  // dyn_em/module_initialize_real.F, although changes have been
608  // made to accommodate interpolation relative to height instead of
609  // pressure.
610  int imin_orig = amrex::lbound(amrex::Box(orig_data)).x;
611  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
612  int jmin_orig = amrex::lbound(amrex::Box(orig_data)).y;
613  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
614  // The number of origin levels, 0 .. ubound: the loops below run k < kmax_orig, so taking
615  // ubound itself here would leave the top origin level out
616  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z + 1;
617  int kmax_new = kmax + 1;
618 
619  AMREX_ASSERT(kmax_orig <= 256);
620  AMREX_ASSERT(kmax_new < 256);
621 
622  amrex::GpuArray<amrex::Real,256> new_z;
623  amrex::GpuArray<amrex::Real,256> new_p;
624  amrex::GpuArray<amrex::Real,256> new_data;
625  amrex::Real* new_z_p = new_z.data();
626  amrex::Real* new_p_p = new_p.data();
627  amrex::Real* new_data_p = new_data.data();
628  for (int k=0; k < kmax_new; k++) {
629  if (stag == 'X') {
630  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i,j+1,k)+new_z_full(i,j,k+1)+new_z_full(i,j+1,k+1));
631  } else if (stag == 'Y') {
632  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i+1,j,k)+new_z_full(i,j,k+1)+new_z_full(i+1,j,k+1));
633  } else if (stag == 'M') {
634  new_z_p[k] = amrex::Real(0.125)*(new_z_full(i,j,k )+new_z_full(i,j+1,k )+new_z_full(i+1,j,k )+new_z_full(i+1,j+1,k )+
635  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));
636  }
637  calc_p_isothermal(new_z_p[k], new_p_p[k]);
638  }
639 
640  amrex::GpuArray<amrex::Real,256> orig_z;
641  amrex::Real* orig_z_p = orig_z.data();
642  for (int k=0; k < kmax_orig; k++) {
643  if (stag == 'M') {
644  orig_z_p[k] = orig_z_full(i,j,k);
645  } else if (stag == 'X') {
646  if (i <= imin_orig) {
647  orig_z_p[k] = orig_z_full(i,j,k);
648  } else if (i >= imax_orig) {
649  orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
650  } else {
651  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
652  }
653  } else if (stag == 'Y') {
654  if (j <= jmin_orig) {
655  orig_z_p[k] = orig_z_full(i,j,k);
656  } else if (j >= jmax_orig) {
657  orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
658  } else {
659  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
660  }
661  }
662  }
663 
664  // Check if the data is top-down instead of bottom-up.
665  bool flip_data_required = false;
666  if (orig_z_p[1] > orig_z_p[kmax_orig-1]) flip_data_required = true;
667  if (flip_data_required) amrex::Abort("metgrid initialization flip_data_required. Not yet implemented.");
668 
669  // Search for the first level above the surface in the origin data.
670  // This is needed since the origin model topography will be
671  // different than the topography processed by WPS.
672  int k_above_sfc = 0;
673  for (int k=1; k < kmax_orig; k++) {
674  if (orig_z_p[k] > orig_z_p[0]) {
675  k_above_sfc = k;
676  break;
677  }
678  }
679 
680  int kend_order;
681  amrex::GpuArray<amrex::Real,256> ordered_z;
682  amrex::GpuArray<amrex::Real,256> ordered_data;
683  amrex::Real* ordered_z_p = ordered_z.data();
684  amrex::Real* ordered_data_p = ordered_data.data();
685  if (k_above_sfc > 1) {
686  // The levels are not monotonically increasing in height, so
687  // we sort and then make "artistic" quality control choices.
688  int count = 0;
689 
690  // Insert levels that are below the surface.
691  for (int k=1; k < k_above_sfc; k++) {
692  ordered_z_p[count] = orig_z_p[k];
693  ordered_data_p[count] = orig_data(i,j,k,src_comp);
694  count++;
695  }
696 
697  // Check if the level that is nearest to and below the surface
698  // is "too close". If so, we'll ignore the upper level and keep
699  // the lower. Origin data is likely to be on pressure levels
700  // with higher spatial resolution near-surface, which supports
701  // the choice of eliminating levels that are "too close" in
702  // pressure-space. For simplicity, calculate delta P assuming a
703  // baroclinic atmosphere.
704  amrex::Real Pu, Pl;
705  calc_p_isothermal(orig_z_p[0], Pu);
706  calc_p_isothermal(ordered_z_p[count-1], Pl);
707  if (Pl-Pu < metgrid_proximity) {
708  count--;
709  }
710 
711  // Insert the surface level.
712  ordered_z_p[count] = orig_z_p[0];
713  ordered_data_p[count] = orig_data(i,j,0,src_comp);
714  count++;
715 
716  // Quoting WRF's comments, the next level to use is at,
717  // "... ta da, the first level above the surface. I know, wow."
718  int knext = k_above_sfc;
719  // Conditionally more strongly use the surface data by removing
720  // levels between the surface and the height corresponding to a
721  // set number of ERF levels from the surface. This forces the
722  // interpolator to use the surface data up through a number of
723  // ERF levels from the surface.
724  if (metgrid_force_sfc_k > 0) {
725  for (int k=k_above_sfc; k < kmax_orig; k++) {
726  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
727  knext = k;
728  break;
729  }
730  }
731  }
732 
733  // Check if the level that is nearest to and above the surface
734  // is "too close". If so, we'll ignore that level.
735  calc_p_isothermal(orig_z_p[knext], Pu);
736  calc_p_isothermal(ordered_z_p[count-1], Pl);
737  if (Pl-Pu < metgrid_proximity) {
738  knext++;
739  }
740 
741  // Insert levels that are above the surface.
742  for (int k=knext; k < kmax_orig; k++) {
743  ordered_z_p[count] = orig_z_p[k];
744  ordered_data_p[count] = orig_data(i,j,k,src_comp);
745  count++;
746  }
747 
748  kend_order = count;
749  } else {
750  // The surface is the lowest level in the origin data.
751 
752  // Insert the surface.
753  ordered_z_p[0] = orig_z_p[0];
754  ordered_data_p[0] = orig_data(i,j,0,src_comp);
755 
756  // Similar to above, conditionally more strongly use the
757  // surface data.
758  int count = 1;
759  int knext = count;
760  // The same ERF level as in the branch above: the surface is used through
761  // metgrid_force_sfc_k ERF levels, k = 0 .. metgrid_force_sfc_k-1
762  if (metgrid_force_sfc_k > 0) {
763  for (int k=knext; k < kmax_orig; k++) {
764  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
765  knext = k;
766  break;
767  }
768  }
769  }
770 
771  // Insert the remaining levels, again ignoring levels that are
772  // "too close" to the prior valid level.
773  for (int k=knext; k < kmax_orig; k++) {
774  amrex::Real Pu, Pl;
775  calc_p_isothermal(orig_z_p[k], Pu);
776  calc_p_isothermal(ordered_z_p[count-1], Pl);
777  if (Pl-Pu < metgrid_proximity) {
778  continue;
779  }
780  ordered_z_p[count] = orig_z_p[k];
781  ordered_data_p[count] = orig_data(i,j,k,src_comp);
782  count++;
783  }
784  kend_order = count;
785  }
786 
787  int ksta(0), kend(0);
788  if (metgrid_use_below_sfc && metgrid_use_sfc) {
789  // Use all levels.
790  ksta = 0;
791  kend = kend_order-1;
792  } else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
793  // Use all levels except for the surface.
794  int ksfc = 0;
795  for (int k=0; k < kmax_orig; k++) {
796  if (ordered_z_p[k] == orig_z_p[0]) {
797  ksfc = k;
798  break;
799  }
800  }
801  for (int k=ksfc; k < kmax_orig-1; k++) {
802  ordered_z_p[k] = ordered_z_p[k+1];
803  ordered_data_p[k] = ordered_data_p[k+1];
804  }
805  ksta = 0;
806  kend = kend_order-2;
807  } else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
808  // Use all levels above and including the surface.
809  int ksfc = 0;
810  for (int k=0; k < kend_order; k++) {
811  if (ordered_z_p[k] == orig_z_p[0]) {
812  ksfc = k;
813  break;
814  }
815  }
816  ksta = ksfc;
817  kend = kend_order-1;
818  } else {
819  // We shouldn't be in here!
820  amrex::Abort("metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
821  }
822 
823  // Insert the level of maximum winds.
824 // amrex::Real maxw_above_this_level = amrex::Real(30000.0);
825 // amrex::Real maxw_horiz_pres_diff = amrex::Real(5000.0);
826 // if ((flag_maxw == 1) && (use_maxw)) {
827 // amrex::Abort("metgrid initialization, use_maxw not yet implemented");
828 // }
829 
830  // Insert the level of the tropopause.
831 // amrex::Real trop_horiz_pres_diff = amrex::Real(5000.0);
832 // if ((flag_trop == 1) && (use_trop)) {
833 // amrex::Abort("metgrid initialization, use_trop not yet implemented");
834 // }
835 
836  amrex::GpuArray<amrex::Real,256> ordered_p;
837  amrex::Real* ordered_p_p = ordered_p.data();
838  for (int k=0; k < kend_order; k++) {
839  calc_p_isothermal(ordered_z_p[k], ordered_p_p[k]);
840  }
841 
842  // kfinal is the index of the topmost level retained in the final_*
843  // arrays, so those arrays hold kfinal+1 points.
844  int kfinal = 0;
845  amrex::GpuArray<amrex::Real,256> final_z;
846  amrex::GpuArray<amrex::Real,256> final_p;
847  amrex::GpuArray<amrex::Real,256> final_data;
848  amrex::Real* final_z_p = final_z.data();
849  amrex::Real* final_p_p = final_p.data();
850  amrex::Real* final_data_p = final_data.data();
851  final_z_p[0] = ordered_z[ksta];
852  final_p_p[0] = ordered_p[ksta];
853  final_data_p[0] = ordered_data[ksta];
854  for (int k=ksta+1; k <= kend; k++) {
855  // Skip levels that are "too close" to the prior retained level.
856  if ((final_p_p[kfinal]-ordered_p_p[k]) < metgrid_proximity) continue;
857  kfinal++;
858  final_z_p[kfinal] = ordered_z_p[k];
859  final_p_p[kfinal] = ordered_p_p[k];
860  final_data_p[kfinal] = ordered_data_p[k];
861  }
862 
863  // Call the interpolator.
864  lagrange_setup(var_type,
865  exp_interp,
866  kfinal+1,
867  kmax_new,
868  metgrid_order,
869  i,
870  j,
871  final_z_p,
872  final_p_p,
873  final_data_p,
874  new_z_p,
875  new_p_p,
876  new_data_p);
877 
878  // Optionally replace the lowest level of data with the surface
879  // field from the origin data. That is origin level 0: ordered_data[0] is
880  // the surface only when no origin level lies below it.
881  if (metgrid_retain_sfc) new_data[0] = orig_data(i,j,0,src_comp);
882 
883  // Save the interpolated data.
884  for (int k=0; k < kmax_new; k++) {
885  new_data_full(i,j,k,dst_comp) = new_data[k];
886  }
887 
888 }
calc_p_isothermal
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:575
lagrange_setup
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_setup(char var_type, const bool &exp_interp, const int &orig_npts, const int &new_npts, const int &order, const int &i, const int &j, amrex::Real *orig_x_z, amrex::Real *orig_x_p, amrex::Real *orig_y, amrex::Real *new_x_z, amrex::Real *new_x_p, amrex::Real *new_y)
Definition: ERF_MetgridUtils.H:268
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
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 
)

Linearly interpolate one metgrid column value onto an ERF point.

904 {
905  // This subroutine is a bit ham-handed and can be cleaned up later.
906  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
907  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
908  // The number of origin levels, 0 .. ubound, so that the top one is searched too
909  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z + 1;
910 
911  amrex::Real z;
912  if (stag == 'X') {
913  z = fourth*(new_z(i,j,k)+new_z(i,j+1,k)+new_z(i,j,k+1)+new_z(i,j+1,k+1));
914  }
915  else if (stag == 'Y') {
916  z = fourth*(new_z(i,j,k)+new_z(i+1,j,k)+new_z(i,j,k+1)+new_z(i+1,j,k+1));
917  }
918  else if (stag == 'M') {
919  z = amrex::Real(0.125)*(new_z(i,j,k )+new_z(i,j+1,k )+new_z(i+1,j,k )+new_z(i+1,j+1,k )+
920  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));
921  }
922 
923  // Initialized to NaN so that any path which fails to set them fails loudly
924  // (with amrex.fpe_trap_invalid=1) rather than returning plausible garbage.
925  amrex::Real z0 = std::numeric_limits<amrex::Real>::quiet_NaN();
926  amrex::Real z1 = std::numeric_limits<amrex::Real>::quiet_NaN();
927  int klow = -1;
928  int khi0 = -1;
929  amrex::Real dzlow = amrex::Real(1.0e12);
930  amrex::Real dzhi0 = -amrex::Real(1.0e12);
931  for (int kk = 0; kk < kmax_orig; kk++) {
932  amrex::Real orig_z_stag = zero;
933  if (stag == 'M') {
934  orig_z_stag = orig_z(i,j,kk);
935  }
936  if (stag == 'X') {
937  if (i == 0) {
938  orig_z_stag = orig_z(i,j,kk);
939  }
940  else if (i == imax_orig) {
941  orig_z_stag = orig_z(imax_orig-1,j,kk);
942  }
943  else {
944  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
945  }
946  }
947  else if (stag == 'Y') {
948  if (j == 0) {
949  orig_z_stag = orig_z(i,j,kk);
950  }
951  else if (j == jmax_orig) {
952  orig_z_stag = orig_z(i,jmax_orig-1,kk);
953  }
954  else {
955  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
956  }
957  }
958 
959  amrex::Real dz = z - orig_z_stag;
960  if ((dz < zero) && (dz > dzhi0)) {
961  dzhi0 = dz;
962  khi0 = kk;
963  z1 = orig_z_stag;
964  }
965  if ((dz >= zero) && (dz < dzlow)) {
966  dzlow = dz;
967  klow = kk;
968  z0 = orig_z_stag;
969  }
970  } // kk
971 
972  // extrapolate below the bottom surface
973  if (klow == -1) {
974  // klow was never found, so z0 was never set above. The two levels used
975  // here are khi0 and khi1, both above z, so z0 is the height at khi0 --
976  // which is exactly what the search loop left in z1 before we clobber it.
977  z0 = z1;
978  int khi1 = -1;
979  amrex::Real dzhi1 = -amrex::Real(1.0e12);
980  for (int kk = 0; kk < kmax_orig; kk++) {
981  amrex::Real orig_z_stag = zero;
982  if (stag == 'M') {
983  orig_z_stag = orig_z(i,j,kk);
984  }
985  else if (stag == 'X') {
986  if (i == 0) {
987  orig_z_stag = orig_z(i,j,kk);
988  }
989  else if (i == imax_orig) {
990  orig_z_stag = orig_z(imax_orig-1,j,kk);
991  }
992  else {
993  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
994  }
995  }
996  else if (stag == 'Y') {
997  if (j == 0) {
998  orig_z_stag = orig_z(i,j,kk);
999  }
1000  else if (j == jmax_orig) {
1001  orig_z_stag = orig_z(i,jmax_orig-1,kk);
1002  }
1003  else {
1004  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
1005  }
1006  }
1007  amrex::Real dz = z - orig_z_stag;
1008  if ((dz < zero) && (dz > dzhi1) && (kk != khi0)) {
1009  dzhi1 = dz;
1010  khi1 = kk;
1011  z1 = orig_z_stag;
1012  }
1013  }
1014  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
1015  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
1016  return ( y0-(y1-y0)/(z1-z0)*(z0-z) );
1017 
1018  // Extrapolate above the top surface
1019  } else if (khi0 == -1) {
1020  khi0 = klow - 1;
1021  int khi1 = klow;
1022  // khi0 was never found, so z1 was never set above. The upper level used
1023  // here is khi1 = klow, whose height is currently held in z0 -- capture
1024  // it before z0 is recomputed at the new (lower) khi0.
1025  z1 = z0;
1026  if (stag == 'M') {
1027  z0 = orig_z(i,j,khi0);
1028  }
1029  else if (stag == 'X') {
1030  if (i == 0) {
1031  z0 = orig_z(i,j,khi0);
1032  }
1033  else if (i == imax_orig) {
1034  z0 = orig_z(imax_orig-1,j,khi0);
1035  }
1036  else {
1037  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
1038  }
1039  }
1040  else if (stag == 'Y') {
1041  if (j == 0) {
1042  z0 = orig_z(i,j,khi0);
1043  }
1044  else if (j == jmax_orig) {
1045  z0 = orig_z(i,jmax_orig-1,khi0);
1046  }
1047  else {
1048  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
1049  }
1050  }
1051  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
1052  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
1053  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
1054  } else {
1055  // interpolate
1056  amrex::Real y0 = orig_data(i,j,klow,src_comp);
1057  amrex::Real y1 = orig_data(i,j,khi0,src_comp);
1058  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
1059 
1060  }
1061 }
Coord::z
@ z
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:12
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:273

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

Interpolate a value with a Lagrange polynomial stencil.

236 {
237  // Interpolation using Lagrange polynomials.
238  // P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
239  // where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
240  // ---------------------------------------------
241  // (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
242  amrex::Real Px = zero;
243  for (int i=0; i <= order; i++) {
244  amrex::Real n = one;
245  amrex::Real d = one;
246  for (int k=0; k <= order; k++) {
247  if (k == i) continue;
248  n *= new_x-x[k];
249  d *= x[i]-x[k];
250  }
251  if (d != zero) {
252  Px += y[i]*n/d;
253  }
254  }
255  new_y = Px;
256 }
Coord::y
@ y
Coord::x
@ x
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37

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_npts,
const int &  new_npts,
const int &  order,
const int &  i,
const int &  j,
amrex::Real *  orig_x_z,
amrex::Real *  orig_x_p,
amrex::Real *  orig_y,
amrex::Real *  new_x_z,
amrex::Real *  new_x_p,
amrex::Real *  new_y 
)

Prepare interpolation stencils and apply Lagrange interpolation or extrapolation.

Both orig_npts and new_npts are point counts, not last indices; valid indices into the orig_* arrays are 0 through orig_npts-1, and into the new_* arrays are 0 through new_npts-1.

281 {
282 
283  amrex::Real CRC_const1 = amrex::Real(11880.516); // m
284  amrex::Real CRC_const2 = amrex::Real(0.1902632);
285  amrex::Real CRC_const3 = amrex::Real(0.0065); // K km-1
286 
287  amrex::ignore_unused(i,j);
288 #ifndef AMREX_USE_GPU
289  bool debug = false;
290 #endif
291 
292  for (int new_k=0; new_k < new_npts; new_k++) {
293 #ifndef AMREX_USE_GPU
294  if (debug) amrex::Print() << "new_k=" << new_k;
295 #endif
296  // Determine if interpolating or extrapolating.
297  // If interpolating, find bounding x values and store the indices.
298  bool extrapolating = true;
299  int kl, kr;
300  for (int ko=0; ko < orig_npts-1; ko++) {
301  amrex::Real a = new_x_z[new_k]-orig_x_z[ko];
302  amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1];
303  if (a*b <= zero) {
304  kl = ko;
305  kr = ko+1;
306  extrapolating = false;
307  break;
308  }
309  }
310 
311  if (extrapolating) {
312  if (var_type == 'T') {
313  // Assume a standard atmosphere -amrex::Real(6.5) K km-1 lapse rate.
314  // Comparable to the WRF default, t_extrap_type=two
315  amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[0];
316  amrex::Real avg_of_extrap_p = myhalf*(new_x_p[new_k]+orig_x_p[0]);
317  amrex::Real temp_extrap_starting_point = orig_y[0]*std::pow(orig_x_p[0]/amrex::Real(100000.0), RdoCp);
318  amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/amrex::Real(100.0), CRC_const2-one);
319  amrex::Real dZ = dZdP*(depth_of_extrap_in_p/amrex::Real(100.0));
320  amrex::Real dT = dZ*CRC_const3;
321  new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(amrex::Real(100000.0)/new_x_p[new_k], RdoCp);
322  } else {
323  // Use a constant value below ground.
324  // Comparable to the WRF default, extrap_type=two
325  new_y[new_k] = orig_y[0];
326  }
327  continue;
328  }
329 
330  if (order%2 != 0) {
331  if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_npts-1)) {
332  // Odd order interpolation.
333  int ksta = kl-(((order+1)/2)-1);
334  int kend = ksta+order;
335 #ifndef AMREX_USE_GPU
336  int ksize = kend-ksta;
337  if (debug) amrex::Print() << " (1a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
338 #endif
339  amrex::Real new_x;
340  amrex::GpuArray<amrex::Real,9> orig_x_sub;
341  amrex::Real* orig_x_sub_p = orig_x_sub.data();
342  if (exp_interp) {
343  new_x = new_x_p[new_k];
344  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
345  } else {
346  new_x = new_x_z[new_k];
347  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
348  }
349  amrex::GpuArray<amrex::Real,9> orig_y_sub;
350  amrex::Real* orig_y_sub_p = orig_y_sub.data();
351  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
352 #ifndef AMREX_USE_GPU
353  if (debug) {
354  amrex::Print() << " orig_x_sub_p = [";
355  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
356  amrex::Print() << "]" << std::endl;
357  amrex::Print() << " orig_y_sub_p = [";
358  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
359  amrex::Print() << "]" << std::endl;
360  }
361 #endif
362  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
363  } else {
364  // Interpolation stencil is too big due to proximity to the surface or model top.
365  // We have bounding points so resort to a simple linear interpolation.
366  int ksta = kl;
367  int kend = kr;
368 #ifndef AMREX_USE_GPU
369  int ksize = kend-ksta+1;
370  if (debug) amrex::Print() << " (1b) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
371 #endif
372  amrex::Real new_x;
373  amrex::GpuArray<amrex::Real,2> orig_x_sub;
374  amrex::GpuArray<amrex::Real,2> orig_y_sub;
375  amrex::Real* orig_x_sub_p = orig_x_sub.data();
376  amrex::Real* orig_y_sub_p = orig_y_sub.data();
377  if (exp_interp) {
378  new_x = new_x_p[new_k];
379  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
380  } else {
381  new_x = new_x_z[new_k];
382  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
383  }
384  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
385 #ifndef AMREX_USE_GPU
386  if (debug) {
387  amrex::Print() << " orig_x_sub_p = [";
388  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
389  amrex::Print() << "]" << std::endl;
390  amrex::Print() << " orig_y_sub_p = [";
391  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
392  amrex::Print() << "]" << std::endl;
393  }
394 #endif
395  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
396  }
397  } else if (order%2 == 0) {
398  if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_npts-1)) {
399  // Even order interpolation.
400  // Interpolate twice. Offset the range by 1 the 2nd time. Average the results.
401  amrex::Real new_y_l, new_y_r;
402  {
403  int ksta = kl-(order/2-1);
404  int kend = ksta+order;
405 #ifndef AMREX_USE_GPU
406  int ksize = kend-ksta;
407  if (debug) amrex::Print() << " (2a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
408 #endif
409  amrex::Real new_x;
410  amrex::GpuArray<amrex::Real,10> orig_x_sub;
411  amrex::GpuArray<amrex::Real,10> orig_y_sub;
412  amrex::Real* orig_x_sub_p = orig_x_sub.data();
413  amrex::Real* orig_y_sub_p = orig_y_sub.data();
414  if (exp_interp) {
415  new_x = new_x_p[new_k];
416  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
417  } else {
418  new_x = new_x_z[new_k];
419  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
420  }
421  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
422 #ifndef AMREX_USE_GPU
423  if (debug) {
424  amrex::Print() << " orig_x_sub_p = [";
425  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
426  amrex::Print() << "]" << std::endl;
427  amrex::Print() << " orig_y_sub_p = [";
428  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
429  amrex::Print() << "]" << std::endl;
430  }
431 #endif
432  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_l);
433  }
434  {
435  int ksta = kl-order/2;
436  int kend = ksta+order;
437 #ifndef AMREX_USE_GPU
438  int ksize = kend-ksta;
439  if (debug) amrex::Print() << "new_k=" << new_k << " (2b) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
440 #endif
441  amrex::Real new_x;
442  amrex::GpuArray<amrex::Real,10> orig_x_sub;
443  amrex::GpuArray<amrex::Real,10> orig_y_sub;
444  amrex::Real* orig_x_sub_p = orig_x_sub.data();
445  amrex::Real* orig_y_sub_p = orig_y_sub.data();
446  if (exp_interp) {
447  new_x = new_x_p[new_k];
448  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
449  } else {
450  new_x = new_x_z[new_k];
451  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
452  }
453  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
454 #ifndef AMREX_USE_GPU
455  if (debug) {
456  amrex::Print() << " orig_x_sub_p = [";
457  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
458  amrex::Print() << "]" << std::endl;
459  amrex::Print() << " orig_y_sub_p = [";
460  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
461  amrex::Print() << "]" << std::endl;
462  }
463 #endif
464  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_r);
465  }
466  new_y[new_k] = myhalf*(new_y_l+new_y_r);
467  } else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_npts-1)) {
468  int ksta = kl-(order/2-1);
469  int kend = ksta+order;
470 #ifndef AMREX_USE_GPU
471  int ksize = kend-ksta;
472  if (debug) amrex::Print() << " (3) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
473 #endif
474  amrex::Real new_x;
475  amrex::GpuArray<amrex::Real,10> orig_x_sub;
476  amrex::GpuArray<amrex::Real,10> orig_y_sub;
477  amrex::Real* orig_x_sub_p = orig_x_sub.data();
478  amrex::Real* orig_y_sub_p = orig_y_sub.data();
479  if (exp_interp) {
480  new_x = new_x_p[new_k];
481  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
482  } else {
483  new_x = new_x_z[new_k];
484  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
485  }
486  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
487 #ifndef AMREX_USE_GPU
488  if (debug) {
489  amrex::Print() << " orig_x_sub_p = [";
490  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
491  amrex::Print() << "]" << std::endl;
492  amrex::Print() << " orig_y_sub_p = [";
493  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
494  amrex::Print() << "]" << std::endl;
495  }
496 #endif
497  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
498  } else if ((kl-order/2 >= 0) && (kr+order/2-1 <= orig_npts-1)) {
499  int ksta = kl-order/2;
500  int kend = ksta+order;
501 #ifndef AMREX_USE_GPU
502  int ksize = kend-ksta;
503  if (debug) amrex::Print() << " (4) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
504 #endif
505  amrex::Real new_x;
506  amrex::GpuArray<amrex::Real,10> orig_x_sub;
507  amrex::GpuArray<amrex::Real,10> orig_y_sub;
508  amrex::Real* orig_x_sub_p = orig_x_sub.data();
509  amrex::Real* orig_y_sub_p = orig_y_sub.data();
510  if (exp_interp) {
511  new_x = new_x_p[new_k];
512  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
513  } else {
514  new_x = new_x_z[new_k];
515  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
516  }
517  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
518 #ifndef AMREX_USE_GPU
519  if (debug) {
520  amrex::Print() << " orig_x_sub_p = [";
521  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
522  amrex::Print() << "]" << std::endl;
523  amrex::Print() << " orig_y_sub_p = [";
524  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
525  amrex::Print() << "]" << std::endl;
526  }
527 #endif
528  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
529  } else {
530  // Linear interpolation.
531  int ksta = kl;
532  int kend = kr;
533 #ifndef AMREX_USE_GPU
534  int ksize = kend-ksta+1;
535  if (debug) amrex::Print() << " (5) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
536 #endif
537  amrex::Real new_x;
538  amrex::GpuArray<amrex::Real,2> orig_x_sub;
539  amrex::GpuArray<amrex::Real,2> orig_y_sub;
540  amrex::Real* orig_x_sub_p = orig_x_sub.data();
541  amrex::Real* orig_y_sub_p = orig_y_sub.data();
542  if (exp_interp) {
543  new_x = new_x_p[new_k];
544  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
545  } else {
546  new_x = new_x_z[new_k];
547  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
548  }
549  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
550 #ifndef AMREX_USE_GPU
551  if (debug) {
552  amrex::Print() << " orig_x_sub_p = [";
553  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
554  amrex::Print() << "]" << std::endl;
555  amrex::Print() << " orig_y_sub_p = [";
556  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
557  amrex::Print() << "]" << std::endl;
558  }
559 #endif
560  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
561  }
562  }
563 #ifndef AMREX_USE_GPU
564  if (debug) amrex::Print() << " new_y[" << new_k << "]=" << new_y[new_k] << std::endl;
565 #endif
566  }
567 }
RdoCp
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:57
lagrange_interp
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_interp(const int &order, amrex::Real *x, amrex::Real *y, amrex::Real &new_x, amrex::Real &new_y)
Definition: ERF_MetgridUtils.H:231

Referenced by interpolate_column_metgrid().

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

◆ metgrid_surface_pressure()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool metgrid_surface_pressure ( const bool  debug_psfc,
const int  flag_psfc,
const amrex::Real  file_psfc,
const amrex::Real  z_sfc,
const amrex::Real  P00,
const amrex::Real  T00,
const amrex::Real  TLP,
amrex::Real &  pressure 
)

Return the physical surface pressure used by the Metgrid initialization.

A debug pressure deliberately overrides file data. Otherwise a PSFC value from the current forcing file is preferred; when PSFC is absent, invert the standard-atmosphere temperature profile at the physical terrain height.

38 {
39  pressure = std::numeric_limits<amrex::Real>::quiet_NaN();
40  if (debug_psfc) {
41  pressure = amrex::Real(100000.0);
42  return true;
43  }
44  if (flag_psfc == 1) {
45  if (!amrex::Math::isfinite(file_psfc) || file_psfc <= amrex::Real(0.0)) {
46  return false;
47  }
48  pressure = file_psfc;
49  return true;
50  }
51 
52  if (!amrex::Math::isfinite(z_sfc) || !amrex::Math::isfinite(P00) || P00 <= amrex::Real(0.0) ||
53  !amrex::Math::isfinite(T00) || T00 <= amrex::Real(0.0) ||
54  !amrex::Math::isfinite(TLP) || TLP <= amrex::Real(0.0)) {
55  return false;
56  }
57  const amrex::Real ToA = T00 / TLP;
58  const amrex::Real radicand =
59  ToA*ToA - amrex::Real(2.0) * CONST_GRAV * z_sfc / (TLP * R_d);
60  if (!amrex::Math::isfinite(radicand) || radicand < amrex::Real(0.0)) {
61  return false;
62  }
63  const amrex::Real exponent = -ToA + std::sqrt(radicand);
64  const amrex::Real result = P00 * std::exp(exponent);
65  if (!amrex::Math::isfinite(result) || result <= amrex::Real(0.0)) {
66  return false;
67  }
68  pressure = result;
69  return true;
70 }

Referenced by metgrid_surface_theta().

Here is the caller graph for this function:

◆ metgrid_surface_theta()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool metgrid_surface_theta ( const amrex::Real  temperature,
const bool  debug_psfc,
const int  flag_psfc,
const amrex::Real  file_psfc,
const amrex::Real  z_sfc,
const amrex::Real  P00,
const amrex::Real  T00,
const amrex::Real  TLP,
const amrex::Real  rdOcp,
amrex::Real &  theta 
)

Convert an absolute Metgrid surface temperature to the canonical SurfaceLayer potential-temperature representation using the selected physical surface pressure.

89 {
90  amrex::Real surface_pressure;
91  if (!metgrid_surface_pressure(debug_psfc, flag_psfc, file_psfc, z_sfc,
92  P00, T00, TLP, surface_pressure)) {
93  return false;
94  }
95  return erf_surface_temperature::temperature_to_theta(
96  temperature, surface_pressure, rdOcp, theta);
97 }
rdOcp
const Real rdOcp
Definition: ERF_InitCustomPert_ABL.H:72
metgrid_surface_pressure
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool metgrid_surface_pressure(const bool debug_psfc, const int flag_psfc, const amrex::Real file_psfc, const amrex::Real z_sfc, const amrex::Real P00, const amrex::Real T00, const amrex::Real TLP, amrex::Real &pressure)
Definition: ERF_MetgridUtils.H:30
MicVar_Kess::theta
@ theta
Definition: ERF_Kessler.H:26
erf_surface_temperature::temperature_to_theta
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool temperature_to_theta(const amrex::Real temperature, const amrex::Real pressure, const amrex::Real rdOcp, amrex::Real &theta)
Definition: ERF_SurfaceTemperature.H:79
Here is the call graph for this function:

◆ read_from_metgrid()

void read_from_metgrid ( int  lev,
int  itime,
const amrex::Box &  domain,
const std::string &  fname,
std::string &  NC_dateTime,
double &  NC_epochTime,
int &  flag_psfc,
int &  flag_msf,
int &  flag_sst,
int &  flag_tsk,
int &  flag_lmask,
int &  NC_nx,
int &  NC_ny,
amrex::Real &  NC_dx,
amrex::Real &  NC_dy,
amrex::FArrayBox &  NC_xvel_fab,
amrex::FArrayBox &  NC_yvel_fab,
amrex::FArrayBox &  NC_temp_fab,
amrex::FArrayBox &  NC_rhum_fab,
amrex::FArrayBox &  NC_pres_fab,
amrex::FArrayBox &  NC_ght_fab,
amrex::FArrayBox &  NC_hgt_fab,
amrex::FArrayBox &  NC_psfc_fab,
amrex::FArrayBox &  NC_msfu_fab,
amrex::FArrayBox &  NC_msfv_fab,
amrex::FArrayBox &  NC_msfm_fab,
amrex::FArrayBox &  NC_sst_fab,
amrex::FArrayBox &  NC_tsk_fab,
amrex::FArrayBox &  NC_LAT_fab,
amrex::FArrayBox &  NC_LON_fab,
amrex::IArrayBox &  NC_lmask_iab,
amrex::Geometry &  geom 
)

Read metgrid NetCDF data and populate ERF initialization FABs.

◆ rh_to_mxrat()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void rh_to_mxrat ( int  i,
int  j,
int  k,
const amrex::Array4< amrex::Real const > &  rhum,
const amrex::Array4< amrex::Real const > &  temp,
const amrex::Array4< amrex::Real const > &  pres,
int  src_indx,
const amrex::Array4< amrex::Real > &  mxrat 
)

Convert relative humidity to vapor mixing ratio using WRF-compatible bounds.

1077 {
1078  amrex::Real qv_max_p_safe = amrex::Real(10000.0); // WRF default value
1079  amrex::Real qv_max_flag = amrex::Real(1.0e-5); // WRF default value
1080  amrex::Real qv_max_value = amrex::Real(3.0e-6); // WRF default value
1081  amrex::Real qv_min_p_safe = amrex::Real(110000.0); // WRF default value
1082  amrex::Real qv_min_flag = amrex::Real(1.0e-6); // WRF default value
1083  amrex::Real qv_min_value = amrex::Real(1.0e-6); // WRF default value
1084  amrex::Real eps = RdoRv;
1085  amrex::Real svp1 = amrex::Real(0.6112);
1086  amrex::Real svp2 = amrex::Real(17.67);
1087  amrex::Real svp3 = amrex::Real(29.65);
1088  amrex::Real svpt0 = amrex::Real(273.15);
1089  // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior)
1090  if (temp(i,j,k) != zero) {
1091  amrex::Real es=amrex::Real(0.01)*rhum(i,j,k)*svp1*amrex::Real(10.0)*std::exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3));
1092  if (es >= pres(i,j,k)/amrex::Real(100.0)) {
1093  // vapor pressure exceeds total pressure
1094  mxrat(i,j,k,src_indx) = amrex::Math::powi<-6>(amrex::Real(10.0));
1095  }
1096  else {
1097  mxrat(i,j,k,src_indx) = amrex::max(eps*es/(pres(i,j,k)/amrex::Real(100.0)-es), amrex::Real(1.0e-6));
1098  }
1099  }
1100  else {
1101  // I don't know why there's a fringe case handled in WRF where T is absolute zero...
1102  // Let's just deal with it here in case we also end up needing it.
1103  mxrat(i,j,k,src_indx) = amrex::Real(1.0e-6);
1104  }
1105  // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1.
1106  // For pressures above a defined level, reasonable Qv values should be
1107  // a certain value or smaller. If they are larger than this, the input data
1108  // probably had "missing" RH, and we filled in some values. This is an
1109  // attempt to catch those. Also, set the minimum value for the entire
1110  // domain that is above the selected pressure level.
1111  if (pres(i,j,k) < qv_max_p_safe) {
1112  if (mxrat(i,j,k,src_indx) > qv_max_flag) {
1113  mxrat(i,j,k,src_indx) = qv_max_value;
1114  }
1115  }
1116  if (pres(i,j,k) < qv_min_p_safe) {
1117  if (mxrat(i,j,k,src_indx) < qv_min_flag) {
1118  mxrat(i,j,k,src_indx) = qv_min_value;
1119  }
1120  }
1121 }
RdoRv
constexpr amrex::Real RdoRv
Definition: ERF_Constants.H:60
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:28
module_model_constants::svp1
real(c_double), parameter svp1
Definition: ERF_module_model_constants.F90:78
module_model_constants::svp3
real(c_double), parameter svp3
Definition: ERF_module_model_constants.F90:80
module_model_constants::svp2
real(c_double), parameter svp2
Definition: ERF_module_model_constants.F90:79
module_model_constants::svpt0
real(c_double), parameter svpt0
Definition: ERF_module_model_constants.F90:81