Skip to content

Commit

Permalink
MetGrid balance HSE with moisture. (#1424)
Browse files Browse the repository at this point in the history 
* MetGrid balance HSE with moisture.

* Minor changes.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Feb 5, 2024
1 parent b4c5d04 commit f9bc46d
Show file tree
Hide file tree
Showing 6 changed files with 120 additions and 50 deletions.
4 changes: 1 addition & 3 deletions Exec/DevTests/MetGrid/inputs_flowmas_metgrid
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,6 @@ zhi.type = "SlipWall"
# TIME STEP CONTROL
erf.fixed_dt = 1.0 # fixed time step depending on grid resolution

erf.use_terrain = true
erf.terrain_smoothing = 2

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = -1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
Expand Down Expand Up @@ -54,6 +51,7 @@ erf.moisture_model = "Kessler"

erf.use_terrain = true
erf.terrain_smoothing = 0
erf.buoyancy_type = 1
erf.terrain_z_levels = 0 53.25373067 116.3299934 190.8583428 278.6665766 381.7196783 502.1312033 642.3454199 \
805.2193582 992.889535 1207.388303 1450.761988 1724.488474 2029.43012 2365.936136 \
2733.669456 3131.532666 3557.60783 4009.253851 4483.532727 4978.420557 5494.911024 \
Expand Down
63 changes: 63 additions & 0 deletions Exec/DevTests/MetGrid/inputs_flowmas_metgrid_NoZlev
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 1000

amrex.fpe_trap_invalid = 1
amrex.fpe_trap_zero = 1
amrex.fpe_trap_overflow = 1

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 447000 387000 18000
amr.n_cell = 149 129 40

geometry.is_periodic = 0 0 0

xlo.type = "Outflow"
xhi.type = "Outflow"
ylo.type = "Outflow"
yhi.type = "Outflow"
zlo.type = "NoSlipWall"
zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.fixed_dt = 1.0 # fixed time step depending on grid resolution

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = -1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
amr.v = 1 # verbosity in Amr.cpp

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed

# CHECKPOINT FILES
erf.check_file = chk # root name of checkpoint file
erf.check_int = 25 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 25 # number of timesteps between plotfiles
erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse dens_hse pert_pres pert_dens rhoQ1 rhoQ2 rhoQ3 qv qc

# SOLVER CHOICE
erf.alpha_T = 1.0
erf.alpha_C = 1.0
erf.use_gravity = true

erf.molec_diff_type = "None"
erf.les_type = "Smagorinsky"
erf.Cs = 0.1

erf.moisture_model = "Kessler"

erf.use_terrain = true
erf.terrain_smoothing = 2
erf.buoyancy_type = 1

# INITIALIZATION WITH ATM DATA
erf.metgrid_bdy_width = 7
erf.metgrid_bdy_set_width = 1
erf.init_type = "metgrid"
erf.nc_init_file_0 = "met_em.d01.2015-04-01_00_00_00.nc" "met_em.d01.2015-04-01_08_00_00.nc"

#There will be no OpenMP tiling
fabarray.mfiter_tile_size = 1024 1024 1024
9 changes: 5 additions & 4 deletions Exec/DevTests/MetGrid/inputs_flowmas_wrfinput
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,6 @@ zhi.type = "SlipWall"
# TIME STEP CONTROL
erf.fixed_dt = 1.0 # fixed time step depending on grid resolution

erf.use_terrain = true
erf.terrain_smoothing = 2

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = -1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
Expand All @@ -38,7 +35,7 @@ erf.check_int = 25 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 1 # number of timesteps between plotfiles
erf.plot_int_1 = 25 # number of timesteps between plotfiles
erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse dens_hse pert_pres pert_dens rhoQ1 rhoQ2 rhoQ3 qv qc

# SOLVER CHOICE
Expand All @@ -52,6 +49,10 @@ erf.Cs = 0.1

erf.moisture_model = "Kessler"

erf.use_terrain = true
erf.terrain_smoothing = 0
erf.buoyancy_type = 1

# INITIALIZATION WITH ATM DATA
erf.wrfbdy_width = 5
erf.wrfbdy_set_width = 1
Expand Down
5 changes: 3 additions & 2 deletions Source/IO/ReadFromMetgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,21 +29,22 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
std::vector<int> attr;
ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0];

/* UNCOMMENT FOR FLOWMAS
/* UNCOMMENT FOR FLOWMAS */
flag_msfu = 0; //ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0];
flag_msfv = 0; //ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0];
flag_msfm = 0; //ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0];
flag_hgt = 1; //ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0];
flag_sst = 0; //ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0];
flag_lmask = 0; //ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0];
*/

/*
ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0];
ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0];
ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0];
ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0];
ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0];
ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0];
*/

ncf.get_attr("WEST-EAST_GRID_DIMENSION", attr); NC_nx = attr[0];
ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0];
Expand Down
57 changes: 30 additions & 27 deletions Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,7 +281,7 @@ ERF::init_from_metgrid (int lev)
FArrayBox &xvel_fab = lev_new[Vars::xvel][mfi];
FArrayBox &yvel_fab = lev_new[Vars::yvel][mfi];
FArrayBox &zvel_fab = lev_new[Vars::zvel][mfi];
FArrayBox &z_phys_nd_fab = (*z_phys)[mfi];
FArrayBox &z_phys_cc_fab = (*z_phys_cc[lev])[mfi];

const Array4<const int>& mask_c_arr = mask_c->const_array(mfi);
const Array4<const int>& mask_u_arr = mask_u->const_array(mfi);
Expand All @@ -296,7 +296,7 @@ ERF::init_from_metgrid (int lev)
init_state_from_metgrid(use_moisture, l_rdOcp,
tbxc, tbxu, tbxv,
cons_fab, xvel_fab, yvel_fab, zvel_fab,
z_phys_nd_fab,
z_phys_cc_fab,
NC_hgt_fab, NC_ght_fab, NC_xvel_fab,
NC_yvel_fab, NC_temp_fab, NC_rhum_fab,
NC_pres_fab, theta_fab, mxrat_fab,
Expand Down Expand Up @@ -823,7 +823,7 @@ init_base_state_from_metgrid (const bool use_moisture,
FArrayBox& r_hse_fab,
FArrayBox& p_hse_fab,
FArrayBox& pi_hse_fab,
FArrayBox& z_phys_nd_fab,
FArrayBox& z_phys_cc_fab,
const Vector<FArrayBox>& NC_ght_fab,
const Vector<FArrayBox>& NC_psfc_fab,
Vector<Vector<FArrayBox>>& fabs_for_bcs,
Expand Down Expand Up @@ -876,7 +876,7 @@ init_base_state_from_metgrid (const bool use_moisture,
valid_bx2d.setRange(2,0);
auto const orig_psfc = NC_psfc_fab[0].const_array();
auto new_data = state_fab.array();
auto const new_z = z_phys_nd_fab.const_array();
auto const new_z = z_phys_cc_fab.const_array();

ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
Expand All @@ -887,29 +887,32 @@ init_base_state_from_metgrid (const bool use_moisture,
}
z_vec[kmax+1] = new_z(i,j,kmax+1);

calc_rho_p(kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Thetad_vec,Thetam_vec,
Q_vec,
z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec);
calc_rho_p(kmax, flag_psfc_vec[0], orig_psfc(i,j,0),
Thetad_vec, Thetam_vec, Q_vec, z_vec,
Rhod_vec, Rhom_vec, Pd_vec, Pm_vec);

for (int k=0; k<=kmax; k++) {
p_hse_arr(i,j,k) = Pd_vec[k];
r_hse_arr(i,j,k) = Rhod_vec[k];
p_hse_arr(i,j,k) = Pm_vec[k];
r_hse_arr(i,j,k) = Rhom_vec[k];
}
});

ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Multiply by Rho to get conserved vars
Real Qv = 0.0;
new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k);
new_data(i,j,k,RhoScalar_comp) = 0.0;
// RhoTheta and RhoQ1 or RhoQv currently hold Theta and Q1 or Qv. Multiply by Rho.
Real RhoTheta = r_hse_arr(i,j,k)*new_data(i,j,k,RhoTheta_comp);
new_data(i,j,k,RhoTheta_comp) = RhoTheta;
new_data(i,j,k,RhoTheta_comp) *= r_hse_arr(i,j,k);
if (use_moisture){
Real RhoQ = r_hse_arr(i,j,k)*new_data(i,j,k,RhoQ_comp);
new_data(i,j,k,RhoQ_comp) = RhoQ;
Qv = new_data(i,j,k,RhoQ_comp);
new_data(i,j,k,RhoQ_comp) *= r_hse_arr(i,j,k);
}
new_data(i,j,k,RhoScalar_comp) = 0.0;

// r_hse needs to include the moisture (account for that here)
r_hse_arr(i,j,k) *= (1.0 + Qv);

pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
//pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp);
});

// FOEXTRAP hse arrays
Expand Down Expand Up @@ -975,11 +978,11 @@ init_base_state_from_metgrid (const bool use_moisture,
Box valid_bx2d = valid_bx;
valid_bx2d.setRange(2,0);
auto const orig_psfc = NC_psfc_fab[it].const_array();
auto const new_z = z_phys_nd_fab.const_array();
auto const new_z = z_phys_cc_fab.const_array();
auto r_arr = fabs_for_bcs[it][MetGridBdyVars::R].array();
auto Theta_arr = fabs_for_bcs[it][MetGridBdyVars::T].array();
auto Q_arr = (use_moisture ) ? fabs_for_bcs[it][MetGridBdyVars::QV].array() : Array4<Real>{};
auto r_hse_arr = fabs_for_bcs[it][MetGridBdyVars::R].array();
auto p_hse_arr = p_hse_bcs_fab.array();
auto p_hse_arr = p_hse_bcs_fab.array();

ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
Expand All @@ -988,18 +991,18 @@ init_base_state_from_metgrid (const bool use_moisture,
Thetad_vec[k] = Theta_arr(i,j,k);
Q_vec[k] = (use_moisture) ? Q_arr(i,j,k) : 0.0;
}
z_vec[kmax+1] = new_z(i,j,kmax+1);
z_vec[kmax+1] = new_z(i,j,kmax+1);

calc_rho_p(kmax,flag_psfc_vec[it],orig_psfc(i,j,0),Thetad_vec,Thetam_vec,
Q_vec,
z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec);
calc_rho_p(kmax, flag_psfc_vec[it], orig_psfc(i,j,0),
Thetad_vec, Thetam_vec, Q_vec, z_vec,
Rhod_vec, Rhom_vec, Pd_vec, Pm_vec);

for (int k=0; k<=kmax; k++) {
p_hse_arr(i,j,k) = Pd_vec[k];
p_hse_arr(i,j,k) = Pm_vec[k];
if (mask_c_arr(i,j,k)) {
r_hse_arr(i,j,k) = Rhod_vec[k];
if (use_moisture) Q_arr(i,j,k) = Rhod_vec[k]*Q_vec[k];
Theta_arr(i,j,k) = Rhod_vec[k]*Thetad_vec[k];
r_arr(i,j,k) = Rhom_vec[k];
if (use_moisture) Q_arr(i,j,k) = Rhom_vec[k]*Q_vec[k];
Theta_arr(i,j,k) = Rhom_vec[k]*Thetad_vec[k];
}
} // k
});
Expand Down
32 changes: 18 additions & 14 deletions Source/Initialization/Metgrid_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ init_base_state_from_metgrid (const bool use_moisture,
amrex::FArrayBox& r_hse_fab,
amrex::FArrayBox& p_hse_fab,
amrex::FArrayBox& pi_hse_fab,
amrex::FArrayBox& z_phys_nd_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,
Expand All @@ -113,46 +113,50 @@ calc_rho_p (const int& kmax,
const int maxiter = 10;

// Calculate or use moist pressure at the surface.
amrex::Real Psurf;
if (flag_psfc == 1) {
Pm_vec[0] = psfc;
Psurf = psfc;
} else {
amrex::Real t_0 = 290.0; // WRF's model_config_rec%base_temp
amrex::Real a = 50.0; // WRF's model_config_rec%base_lapse
Pm_vec[0] = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*CONST_GRAV*z_vec[0]/(a*R_d)), 0.5));
Psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*CONST_GRAV*z_vec[0]/(a*R_d)), 0.5));
}

// calculate virtual potential temperature at the surface.
// Iterations for the first CC point that is 1/2 dz off the surface
{
amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0];
amrex::Real half_dz = z_vec[0];
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[0];
Thetam_vec[0] = Thetad_vec[0]*qvf;
for (int it=0; it<maxiter; it++) {
Pm_vec[0] = Psurf-half_dz*(Rhom_vec[0])*(1.0+Q_vec[0])*CONST_GRAV;
if (Pm_vec[0] < 0.0) Pm_vec[0] = 0.0;
Rhom_vec[0] = (p_0/(R_d*Thetam_vec[0]))*std::pow(Pm_vec[0]/p_0, iGamma);
} // it
}

// calculate moist density at the surface.
Rhom_vec[0] = 1.0/(R_d/p_0*Thetam_vec[0]*std::pow(Pm_vec[0]/p_0, -iGamma));

// integrate from the surface to the top boundary.
// Integrate from the first CC point to the top boundary.
for (int k=1; k<=kmax; k++) {
amrex::Real dz = z_vec[k]-z_vec[k-1];
amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[k];
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[k];
Thetam_vec[k] = Thetad_vec[k]*qvf;
Rhom_vec[k] = Rhom_vec[k-1]; // an initial guess.
for (int it=0; it<maxiter; it++) {
Pm_vec[k] = Pm_vec[k-1]-0.5*dz*(Rhom_vec[k]+Rhom_vec[k-1])*CONST_GRAV;
Pm_vec[k] = Pm_vec[k-1]-0.5*dz*(Rhom_vec[k]+Rhom_vec[k-1])*(1.0+Q_vec[k])*CONST_GRAV;
if (Pm_vec[k] < 0.0) Pm_vec[k] = 0.0;
Rhom_vec[k] = 1.0/(R_d/p_0*Thetam_vec[k]*std::pow(Pm_vec[k]/p_0, -iGamma));
Rhom_vec[k] = (p_0/(R_d*Thetam_vec[k]))*std::pow(Pm_vec[k]/p_0, iGamma);
} // it
} // k

// integrate from the top back down to get dry pressure and density.
Pd_vec[kmax] = Pm_vec[kmax];
Rhod_vec[kmax] = 1.0/(R_d/p_0*Thetam_vec[kmax]*std::pow(Pd_vec[kmax]/p_0, -iGamma));
Rhod_vec[kmax] = (p_0/(R_d*Thetam_vec[kmax]))*std::pow(Pd_vec[kmax]/p_0, iGamma);
for (int k=kmax-1; k>=0; k--) {
amrex::Real dz = z_vec[k+1]-z_vec[k];
Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess.
for (int it=0; it < maxiter; it++) {
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV;
if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0;
Rhod_vec[k] = 1.0/(R_d/p_0*Thetam_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma));
Rhod_vec[k] = (p_0/(R_d*Thetam_vec[k]))*std::pow(Pd_vec[k]/p_0, iGamma);
} // it
} // k
}
Expand Down

0 comments on commit f9bc46d

Please sign in to comment.