ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_Metgrid_utils.H File Reference
#include <ERF.H>
#include <ERF_EOS.H>
#include <ERF_Utils.H>
#include <ERF_prob_common.H>
#include <ERF_HSE_utils.H>
Include dependency graph for ERF_Metgrid_utils.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_msfu, int &flag_msfv, int &flag_msfm, int &flag_hgt, 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::Real &Latitude, amrex::Real &Longitude, 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 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_zvel_fab, const amrex::Vector< amrex::FArrayBox > &NC_temp_fab, const amrex::Vector< amrex::FArrayBox > &NC_rhum_fab, amrex::Vector< amrex::FArrayBox > &theta_fab, amrex::Vector< amrex::FArrayBox > &mxrat_fab, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs, 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 (amrex::FArrayBox &msfu_fab, amrex::FArrayBox &msfv_fab, amrex::FArrayBox &msfm_fab, const int &flag_msfu, const int &flag_msfv, const int &flag_msfm, 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 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_ght_fab, const amrex::Vector< amrex::FArrayBox > &NC_psfc_fab, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real interpolate_column_metgrid (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

◆ init_base_state_from_metgrid()

void init_base_state_from_metgrid ( const bool  use_moisture,
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_ght_fab,
const amrex::Vector< amrex::FArrayBox > &  NC_psfc_fab,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs 
)

◆ init_msfs_from_metgrid()

void init_msfs_from_metgrid ( amrex::FArrayBox &  msfu_fab,
amrex::FArrayBox &  msfv_fab,
amrex::FArrayBox &  msfm_fab,
const int &  flag_msfu,
const int &  flag_msfv,
const int &  flag_msfm,
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 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_zvel_fab,
const amrex::Vector< amrex::FArrayBox > &  NC_temp_fab,
const amrex::Vector< amrex::FArrayBox > &  NC_rhum_fab,
amrex::Vector< amrex::FArrayBox > &  theta_fab,
amrex::Vector< amrex::FArrayBox > &  mxrat_fab,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs,
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 amrex::Real interpolate_column_metgrid ( 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 
)
114 {
115  // This subroutine is a bit ham-handed and can be cleaned up later.
116  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
117  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
118  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
119 
120  amrex::Real z;
121  if (stag == 'X') {
122  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));
123  }
124  else if (stag == 'Y') {
125  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));
126  }
127  else if (stag == 'M') {
128  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 )+
129  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));
130  }
131 
132  amrex::Real z0, z1;
133  int klow = -1;
134  int khi0 = -1;
135  amrex::Real dzlow = 1.0e12;
136  amrex::Real dzhi0 = -1.0e12;
137  for (int kk = 0; kk < kmax_orig; kk++) {
138  amrex::Real orig_z_stag = 0.0;
139  if (stag == 'M') {
140  orig_z_stag = orig_z(i,j,kk);
141  }
142  if (stag == 'X') {
143  if (i == 0) {
144  orig_z_stag = orig_z(i,j,kk);
145  }
146  else if (i == imax_orig) {
147  orig_z_stag = orig_z(imax_orig-1,j,kk);
148  }
149  else {
150  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
151  }
152  }
153  else if (stag == 'Y') {
154  if (j == 0) {
155  orig_z_stag = orig_z(i,j,kk);
156  }
157  else if (j == jmax_orig) {
158  orig_z_stag = orig_z(i,jmax_orig-1,kk);
159  }
160  else {
161  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
162  }
163  }
164 
165  amrex::Real dz = z - orig_z_stag;
166  if ((dz < 0.0) && (dz > dzhi0)) {
167  dzhi0 = dz;
168  khi0 = kk;
169  z1 = orig_z_stag;
170  }
171  if ((dz >= 0.0) && (dz < dzlow)) {
172  dzlow = dz;
173  klow = kk;
174  z0 = orig_z_stag;
175  }
176  } // kk
177 
178  // extrapolate below the bottom surface
179  if (klow == -1) {
180  int khi1 = -1;
181  amrex::Real dzhi1 = -1.0e12;
182  for (int kk = 0; kk < kmax_orig; kk++) {
183  amrex::Real orig_z_stag = 0.0;
184  if (stag == 'M') {
185  orig_z_stag = orig_z(i,j,kk);
186  }
187  else if (stag == 'X') {
188  if (i == 0) {
189  orig_z_stag = orig_z(i,j,kk);
190  }
191  else if (i == imax_orig) {
192  orig_z_stag = orig_z(imax_orig-1,j,kk);
193  }
194  else {
195  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
196  }
197  }
198  else if (stag == 'Y') {
199  if (j == 0) {
200  orig_z_stag = orig_z(i,j,kk);
201  }
202  else if (j == jmax_orig) {
203  orig_z_stag = orig_z(i,jmax_orig-1,kk);
204  }
205  else {
206  orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
207  }
208  }
209  amrex::Real dz = z - orig_z_stag;
210  if ((dz < 0.0) && (dz > dzhi1) && (kk != khi0)) {
211  dzhi1 = dz;
212  khi1 = kk;
213  z1 = orig_z_stag;
214  }
215  }
216  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
217  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
218  return ( y0-(y1-y0)/(z1-z0)*(z0-z) );
219 
220  // Extrapolate above the top surface
221  } else if (khi0 == -1) {
222  khi0 = klow - 1;
223  int khi1 = klow;
224  if (stag == 'M') {
225  z0 = orig_z(i,j,khi0);
226  }
227  else if (stag == 'X') {
228  if (i == 0) {
229  z0 = orig_z(i,j,khi0);
230  }
231  else if (i == imax_orig) {
232  z0 = orig_z(imax_orig-1,j,khi0);
233  }
234  else {
235  z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
236  }
237  }
238  else if (stag == 'Y') {
239  if (j == 0) {
240  z0 = orig_z(i,j,khi0);
241  }
242  else if (j == jmax_orig) {
243  z0 = orig_z(i,jmax_orig-1,khi0);
244  }
245  else {
246  z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
247  }
248  }
249  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
250  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
251  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
252  } else {
253  // interpolate
254  amrex::Real y0 = orig_data(i,j,klow,src_comp);
255  amrex::Real y1 = orig_data(i,j,khi0,src_comp);
256  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
257 
258  }
259 }

◆ 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_msfu,
int &  flag_msfv,
int &  flag_msfm,
int &  flag_hgt,
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::Real &  Latitude,
amrex::Real &  Longitude,
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 
)
271 {
272  amrex::Real qv_max_p_safe = 10000.0; // WRF default value
273  amrex::Real qv_max_flag = 1.0e-5; // WRF default value
274  amrex::Real qv_max_value = 3.0e-6; // WRF default value
275  amrex::Real qv_min_p_safe = 110000.0; // WRF default value
276  amrex::Real qv_min_flag = 1.0e-6; // WRF default value
277  amrex::Real qv_min_value = 1.0e-6; // WRF default value
278  amrex::Real eps = 0.622;
279  amrex::Real svp1 = 0.6112;
280  amrex::Real svp2 = 17.67;
281  amrex::Real svp3 = 29.65;
282  amrex::Real svpt0 = 273.15;
283  // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior)
284  if (temp(i,j,k) != 0.0) {
285  amrex::Real es=0.01*rhum(i,j,k)*svp1*10.0*exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3));
286  if (es >= pres(i,j,k)/100.0) {
287  // vapor pressure exceeds total pressure
288  mxrat(i,j,k) = std::pow(10.0, -6);
289  }
290  else {
291  mxrat(i,j,k) = amrex::max(eps*es/(pres(i,j,k)/100.0-es), 1.0e-6);
292  }
293  }
294  else {
295  // I don't know why there's a fringe case handled in WRF where T is absolute zero...
296  // Let's just deal with it here in case we also end up needing it.
297  mxrat(i,j,k) = 1.0e-6;
298  }
299  // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1.
300  // For pressures above a defined level, reasonable Qv values should be
301  // a certain value or smaller. If they are larger than this, the input data
302  // probably had "missing" RH, and we filled in some values. This is an
303  // attempt to catch those. Also, set the minimum value for the entire
304  // domain that is above the selected pressure level.
305  if (pres(i,j,k) < qv_max_p_safe) {
306  if (mxrat(i,j,k) > qv_max_flag) {
307  mxrat(i,j,k) = qv_max_value;
308  }
309  }
310  if (pres(i,j,k) < qv_min_p_safe) {
311  if (mxrat(i,j,k) < qv_min_flag) {
312  mxrat(i,j,k) = qv_min_value;
313  }
314  }
315 }
@ pres
Definition: ERF_Kessler.H:33