diff --git a/Source/IO/ReadFromWRFInput.cpp b/Source/IO/ReadFromWRFInput.cpp index 50724fb46..f1a6aa63a 100644 --- a/Source/IO/ReadFromWRFInput.cpp +++ b/Source/IO/ReadFromWRFInput.cpp @@ -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; @@ -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 diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index 0da111d0d..7af6cca4f 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -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 @@ -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& NC_QVAPOR_fab, const Vector& NC_QCLOUD_fab, @@ -80,11 +85,18 @@ init_terrain_from_wrfinput (int lev, const Real& z_top, const Vector& 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& NC_ALB_fab, - const Vector& NC_PB_fab); + const Vector& NC_PB_fab, + const Vector& NC_P_fab); /** * ERF function that initializes data from a WRF dataset @@ -110,6 +122,7 @@ ERF::init_from_wrfinput (int lev) Vector NC_C2H_fab ; NC_C2H_fab.resize(num_boxes_at_level[lev]); Vector NC_RDNW_fab ; NC_RDNW_fab.resize(num_boxes_at_level[lev]); Vector NC_PH_fab ; NC_PH_fab.resize(num_boxes_at_level[lev]); + Vector NC_P_fab ; NC_P_fab.resize(num_boxes_at_level[lev]); Vector NC_PHB_fab ; NC_PHB_fab.resize(num_boxes_at_level[lev]); Vector NC_ALB_fab ; NC_ALB_fab.resize(num_boxes_at_level[lev]); Vector NC_PB_fab ; NC_PB_fab.resize(num_boxes_at_level[lev]); @@ -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); } @@ -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); } } @@ -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) @@ -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& NC_ALB_fab, - const Vector& NC_PB_fab) + const Vector& NC_PB_fab, + const Vector& NC_P_fab) { int nboxes = NC_ALB_fab.size(); for (int idx = 0; idx < nboxes; idx++) @@ -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& p_hse_arr = p_hse.array(); - const Array4& pi_hse_arr = pi_hse.array(); const Array4& r_hse_arr = r_hse.array(); - const Array4& alpha_arr = NC_ALB_fab[idx].const_array(); - const Array4& 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& cons_arr = cons_fab.array(); + const Array4& p_hse_arr = p_hse.array(); + const Array4& pi_hse_arr = pi_hse.array(); + const Array4& r_hse_arr = r_hse.array(); + const Array4& alpha_arr = NC_ALB_fab[idx].const_array(); + const Array4& nc_pb_arr = NC_PB_fab[idx].const_array(); + const Array4& 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