Skip to content

Commit

Permalink
Balance pressure.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi committed Feb 6, 2024
1 parent e336d57 commit 76dd3f2
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 53 deletions.
32 changes: 18 additions & 14 deletions Source/IO/ReadFromWRFInput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,11 @@ read_from_wrfinput (int lev,
FArrayBox& NC_QVAPOR_fab,
FArrayBox& NC_QCLOUD_fab,
FArrayBox& NC_QRAIN_fab,
FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab,
FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab,
FArrayBox& NC_PH_fab,
FArrayBox& NC_P_fab,
FArrayBox& NC_PHB_fab,
FArrayBox& NC_ALB_fab,
FArrayBox& NC_PB_fab,
MoistureType moisture_type)
{
amrex::Print() << "Loading initial data from NetCDF file at level " << lev << std::endl;
Expand All @@ -39,20 +42,21 @@ read_from_wrfinput (int lev,
NC_fabs.push_back(&NC_PH_fab); NC_names.push_back("PH"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 6
NC_fabs.push_back(&NC_PHB_fab); NC_names.push_back("PHB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 7
NC_fabs.push_back(&NC_PB_fab); NC_names.push_back("PB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 9
NC_fabs.push_back(&NC_ALB_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 0
NC_fabs.push_back(&NC_MUB_fab); NC_names.push_back("MUB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 10
NC_fabs.push_back(&NC_MSFU_fab); NC_names.push_back("MAPFAC_UY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 11
NC_fabs.push_back(&NC_MSFV_fab); NC_names.push_back("MAPFAC_VY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 12
NC_fabs.push_back(&NC_MSFM_fab); NC_names.push_back("MAPFAC_MY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 13
NC_fabs.push_back(&NC_SST_fab); NC_names.push_back("SST"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 14
NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 15
NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 16
NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17
NC_fabs.push_back(&NC_P_fab); NC_names.push_back("P"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 10
NC_fabs.push_back(&NC_ALB_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 11
NC_fabs.push_back(&NC_MUB_fab); NC_names.push_back("MUB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 12
NC_fabs.push_back(&NC_MSFU_fab); NC_names.push_back("MAPFAC_UY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 13
NC_fabs.push_back(&NC_MSFV_fab); NC_names.push_back("MAPFAC_VY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 14
NC_fabs.push_back(&NC_MSFM_fab); NC_names.push_back("MAPFAC_MY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 15
NC_fabs.push_back(&NC_SST_fab); NC_names.push_back("SST"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 16
NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17
NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 18
NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 19

if (moisture_type != MoistureType::None) {
NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 18
NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 19
NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20
NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20
NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 21
NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 22
}

// Read the netcdf file and fill these FABs
Expand Down
119 changes: 80 additions & 39 deletions Source/Initialization/ERF_init_from_wrfinput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,11 @@ read_from_wrfinput (int lev, const Box& domain, const std::string& fname,
FArrayBox& NC_QVAPOR_fab,
FArrayBox& NC_QCLOUD_fab,
FArrayBox& NC_QRAIN_fab,
FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab,
FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab,
FArrayBox& NC_PH_fab,
FArrayBox& NC_P_fab,
FArrayBox& NC_PHB_fab,
FArrayBox& NC_ALB_fab,
FArrayBox& NC_PB_fab,
MoistureType moisture_type);

Real
Expand All @@ -53,8 +56,10 @@ convert_wrfbdy_data (int which, const Box& domain,
const FArrayBox& NC_rhoth_fab);

void
init_state_from_wrfinput (int lev, FArrayBox& state_fab,
FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
init_state_from_wrfinput (int lev,
FArrayBox& cons_fab,
FArrayBox& x_vel_fab,
FArrayBox& y_vel_fab,
FArrayBox& z_vel_fab,
const Vector<FArrayBox>& NC_QVAPOR_fab,
const Vector<FArrayBox>& NC_QCLOUD_fab,
Expand All @@ -80,11 +85,18 @@ init_terrain_from_wrfinput (int lev, const Real& z_top,
const Vector<FArrayBox>& NC_PHB_fab);

void
init_base_state_from_wrfinput (int lev, const Box& bx, Real l_rdOcp,
FArrayBox& p_hse, FArrayBox& pi_hse,
init_base_state_from_wrfinput (int lev,
const Box& bx,
Real l_rdOcp,
MoistureType moisture_type,
const int& n_qstate,
FArrayBox& cons_fab,
FArrayBox& p_hse,
FArrayBox& pi_hse,
FArrayBox& r_hse,
const Vector<FArrayBox>& NC_ALB_fab,
const Vector<FArrayBox>& NC_PB_fab);
const Vector<FArrayBox>& NC_PB_fab,
const Vector<FArrayBox>& NC_P_fab);

/**
* ERF function that initializes data from a WRF dataset
Expand All @@ -110,6 +122,7 @@ ERF::init_from_wrfinput (int lev)
Vector<FArrayBox> NC_C2H_fab ; NC_C2H_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_RDNW_fab ; NC_RDNW_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_PH_fab ; NC_PH_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_P_fab ; NC_P_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_PHB_fab ; NC_PHB_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_ALB_fab ; NC_ALB_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_PB_fab ; NC_PB_fab.resize(num_boxes_at_level[lev]);
Expand All @@ -124,12 +137,13 @@ ERF::init_from_wrfinput (int lev)
for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
{
read_from_wrfinput(lev, boxes_at_level[lev][idx], nc_init_file[lev][idx],
NC_xvel_fab[idx], NC_yvel_fab[idx], NC_zvel_fab[idx], NC_rho_fab[idx],
NC_rhop_fab[idx], NC_rhoth_fab[idx], NC_MUB_fab[idx],
NC_MSFU_fab[idx], NC_MSFV_fab[idx], NC_MSFM_fab[idx],
NC_SST_fab[idx], NC_C1H_fab[idx], NC_C2H_fab[idx], NC_RDNW_fab[idx],
NC_xvel_fab[idx] , NC_yvel_fab[idx] , NC_zvel_fab[idx] , NC_rho_fab[idx],
NC_rhop_fab[idx] , NC_rhoth_fab[idx] , NC_MUB_fab[idx] ,
NC_MSFU_fab[idx] , NC_MSFV_fab[idx] , NC_MSFM_fab[idx] ,
NC_SST_fab[idx] , NC_C1H_fab[idx] , NC_C2H_fab[idx] , NC_RDNW_fab[idx],
NC_QVAPOR_fab[idx], NC_QCLOUD_fab[idx], NC_QRAIN_fab[idx],
NC_PH_fab[idx],NC_PHB_fab[idx],NC_ALB_fab[idx],NC_PB_fab[idx],
NC_PH_fab[idx] , NC_P_fab[idx] , NC_PHB_fab[idx] ,
NC_ALB_fab[idx] , NC_PB_fab[idx] ,
solverChoice.moisture_type);
}

Expand Down Expand Up @@ -195,14 +209,15 @@ ERF::init_from_wrfinput (int lev)
if (init_type == "real") {
for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
FArrayBox& cons_fab = lev_new[Vars::cons][mfi];
FArrayBox& p_hse_fab = p_hse[mfi];
FArrayBox& pi_hse_fab = pi_hse[mfi];
FArrayBox& r_hse_fab = r_hse[mfi];

const Box valid_bx = mfi.validbox();
init_base_state_from_wrfinput(lev, valid_bx, l_rdOcp,
p_hse_fab, pi_hse_fab, r_hse_fab,
NC_ALB_fab, NC_PB_fab);
init_base_state_from_wrfinput(lev, valid_bx, l_rdOcp, solverChoice.moisture_type, n_qstate,
cons_fab, p_hse_fab, pi_hse_fab, r_hse_fab,
NC_ALB_fab, NC_PB_fab, NC_P_fab);
}
}

Expand All @@ -218,21 +233,21 @@ ERF::init_from_wrfinput (int lev)
<< " and relaxation width: " << wrfbdy_width - wrfbdy_set_width << std::endl;

convert_wrfbdy_data(0,domain,bdy_data_xlo,
NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0],
NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0],
NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]);
NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] ,
NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0],
NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0]);
convert_wrfbdy_data(1,domain,bdy_data_xhi,
NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0],
NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0],
NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]);
NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] ,
NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0],
NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0]);
convert_wrfbdy_data(2,domain,bdy_data_ylo,
NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0],
NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0],
NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]);
NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] ,
NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0],
NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0]);
convert_wrfbdy_data(3,domain,bdy_data_yhi,
NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0],
NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0],
NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]);
NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] ,
NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0],
NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0]);
}

// Start at the earliest time (read_from_wrfbdy)
Expand Down Expand Up @@ -360,10 +375,18 @@ init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab,
* @param NC_PB_fab Vector of FArrayBox objects containing WRF data specifying pressure
*/
void
init_base_state_from_wrfinput (int lev, const Box& valid_bx, const Real l_rdOcp,
FArrayBox& p_hse, FArrayBox& pi_hse, FArrayBox& r_hse,
init_base_state_from_wrfinput (int lev,
const Box& valid_bx,
const Real l_rdOcp,
MoistureType moisture_type,
const int& n_qstate,
FArrayBox& cons_fab,
FArrayBox& p_hse,
FArrayBox& pi_hse,
FArrayBox& r_hse,
const Vector<FArrayBox>& NC_ALB_fab,
const Vector<FArrayBox>& NC_PB_fab)
const Vector<FArrayBox>& NC_PB_fab,
const Vector<FArrayBox>& NC_P_fab)
{
int nboxes = NC_ALB_fab.size();
for (int idx = 0; idx < nboxes; idx++)
Expand All @@ -372,16 +395,34 @@ init_base_state_from_wrfinput (int lev, const Box& valid_bx, const Real l_rdOcp,
// FArrayBox to FArrayBox copy does "copy on intersection"
// This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
//
const Array4<Real >& p_hse_arr = p_hse.array();
const Array4<Real >& pi_hse_arr = pi_hse.array(); const Array4<Real >& r_hse_arr = r_hse.array();
const Array4<Real const>& alpha_arr = NC_ALB_fab[idx].const_array();
const Array4<Real const>& nc_pb_arr = NC_PB_fab[idx].const_array();

ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
p_hse_arr(i,j,k) = nc_pb_arr(i,j,k);
const Array4<Real >& cons_arr = cons_fab.array();
const Array4<Real >& p_hse_arr = p_hse.array();
const Array4<Real >& pi_hse_arr = pi_hse.array();
const Array4<Real >& r_hse_arr = r_hse.array();
const Array4<Real const>& alpha_arr = NC_ALB_fab[idx].const_array();
const Array4<Real const>& nc_pb_arr = NC_PB_fab[idx].const_array();
const Array4<Real const>& nc_p_arr = NC_P_fab[idx].const_array();

ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Base plus perturbational pressure
Real Ptot = nc_pb_arr(i,j,k) + nc_p_arr(i,j,k);

// Compute pressure from EOS
Real Qv = (moisture_type != MoistureType::None) ?
cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp) : 0.0;
Real RT = cons_arr(i,j,k,RhoTheta_comp);
Real P_eos = getPgivenRTh(RT, Qv);
Real DelP = std::fabs(Ptot - P_eos);
AMREX_ASSERT_WITH_MESSAGE((DelP < 1.0), "Initial state is inconsistent with EOS!");

// Compute rhse
Real Rhse_Sum = cons_arr(i,j,k,Rho_comp);
for (int q_offset(0); q_offset<n_qstate; ++q_offset) Rhse_Sum += cons_arr(i,j,k,RhoQ1_comp+q_offset);

p_hse_arr(i,j,k) = Ptot;
pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
r_hse_arr(i,j,k) = 1.0 / alpha_arr(i,j,k);

r_hse_arr(i,j,k) = Rhse_Sum;
});
} // idx
}
Expand Down

0 comments on commit 76dd3f2

Please sign in to comment.