Skip to content

Commit

Permalink
Get lat/lon and bounds check with met and wrf init. Need to verify on…
Browse files Browse the repository at this point in the history 
… GPU. (#1497)
  • Loading branch information
@AMLattanzi
AMLattanzi authored Mar 15, 2024
1 parent b7c8ccc commit 24ddedf
Show file tree
Hide file tree
Showing 7 changed files with 234 additions and 54 deletions.
2 changes: 2 additions & 0 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -493,6 +493,8 @@ private:
amrex::Vector<amrex::Vector<amrex::FArrayBox>> bdy_data_yhi;

amrex::Real bdy_time_interval;
amrex::Real Latitude;
amrex::Real Longitude;
#endif // ERF_USE_NETCDF

// Struct for working with the sounding data we take as an input
Expand Down
16 changes: 15 additions & 1 deletion Source/IO/NCWpsFile.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,6 +223,10 @@ void ReadNetCDFFile (const std::string& fname, amrex::Vector<std::string> names,
template<class FAB,typename DType>
void
fill_fab_from_arrays (int iv,
amrex::Real& Latitude,
amrex::Real& Longitude,
std::string& Lat_var_name,
std::string& Lon_var_name,
amrex::Vector<NDArray<float>>& nc_arrays,
const std::string& var_name,
NC_Data_Dims_Type& NC_dim_type,
Expand Down Expand Up @@ -280,6 +284,9 @@ fill_fab_from_arrays (int iv,
int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff;
fab_arr(i,j,k,0) = static_cast<DType>(*(nc_arrays[iv].get_data()+n));
}

if (var_name == Lat_var_name) Latitude = fab_arr(0,0,0);
if (var_name == Lon_var_name) Longitude = fab_arr(0,0,0);
}

/**
Expand All @@ -293,6 +300,10 @@ fill_fab_from_arrays (int iv,
template<class FAB,typename DType>
void
BuildFABsFromNetCDFFile (const amrex::Box& domain,
amrex::Real& Latitude,
amrex::Real& Longitude,
std::string& Lat_var_name,
std::string& Lon_var_name,
const std::string &fname,
amrex::Vector<std::string> nc_var_names,
amrex::Vector<enum NC_Data_Dims_Type> NC_dim_types,
Expand All @@ -311,7 +322,10 @@ BuildFABsFromNetCDFFile (const amrex::Box& domain,
{
FAB tmp;
if (amrex::ParallelDescriptor::IOProcessor()) {
fill_fab_from_arrays<FAB,DType>(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp);
fill_fab_from_arrays<FAB,DType>(iv, Latitude, Longitude,
Lat_var_name, Lon_var_name,
nc_arrays, nc_var_names[iv],
NC_dim_types[iv], tmp);
}

int ncomp = tmp.nComp();
Expand Down
40 changes: 36 additions & 4 deletions Source/IO/ReadFromMetgrid.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include <NCWpsFile.H>
#include <AMReX_FArrayBox.H>
#include <AMReX_IArrayBox.H>
#include <Metgrid_utils.H>

using namespace amrex;

Expand All @@ -19,9 +20,13 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
FArrayBox& NC_hgt_fab, FArrayBox& NC_psfc_fab,
FArrayBox& NC_msfu_fab, FArrayBox& NC_msfv_fab,
FArrayBox& NC_msfm_fab, FArrayBox& NC_sst_fab,
IArrayBox& NC_lmask_iab)
FArrayBox& NC_LAT_fab, FArrayBox& NC_LON_fab,
IArrayBox& NC_lmask_iab,
Real& Latitude,
Real& Longitude,
Geometry& geom)
{
Print() << "Loading initial data from NetCDF file at level " << lev << std::endl;
Print() << "Loading header data from NetCDF file at level " << lev << std::endl;

if (ParallelDescriptor::IOProcessor()) {
auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
Expand Down Expand Up @@ -60,6 +65,23 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
ncf.get_attr("DY", attr); NC_dy = attr[0];
}
ncf.close();

// Verify the inputs geometry matches what the NETCDF file has
Real tol = 1.0e-3;
Real Len_x = NC_dx * Real(NC_nx-1);
Real Len_y = NC_dy * Real(NC_ny-1);
if (std::fabs(Len_x - (geom.ProbHi(0) - geom.ProbLo(0))) > tol) {
Print() << "X problem extent " << (geom.ProbHi(0) - geom.ProbLo(0)) << " does not match NETCDF file "
<< Len_x << "!\n";
Print() << "dx: " << NC_dx << ' ' << "Nx: " << NC_nx-1 << "\n";
Abort("Domain specification error");
}
if (std::fabs(Len_y - (geom.ProbHi(1) - geom.ProbLo(1))) > tol) {
Print() << "Y problem extent " << (geom.ProbHi(1) - geom.ProbLo(1)) << " does not match NETCDF file "
<< Len_y << "!\n";
Print() << "dy: " << NC_dy << ' ' << "Ny: " << NC_ny-1 << "\n";
Abort("Domain specification error");
}
}
int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank
ParallelDescriptor::Bcast(&flag_psfc, 1, ioproc);
Expand All @@ -75,6 +97,8 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
ParallelDescriptor::Bcast(&NC_dx, 1, ioproc);
ParallelDescriptor::Bcast(&NC_dy, 1, ioproc);

Print() << "Loading initial data from NetCDF file at level " << lev << std::endl;

Vector<FArrayBox*> NC_fabs;
Vector<IArrayBox*> NC_iabs;
Vector<std::string> NC_fnames;
Expand All @@ -88,6 +112,8 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
NC_fabs.push_back(&NC_rhum_fab); NC_fnames.push_back("RH"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);
NC_fabs.push_back(&NC_pres_fab); NC_fnames.push_back("PRES"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);
NC_fabs.push_back(&NC_ght_fab); NC_fnames.push_back("GHT"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);
NC_fabs.push_back(&NC_LAT_fab); NC_fnames.push_back("XLAT_V"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE);
NC_fabs.push_back(&NC_LON_fab); NC_fnames.push_back("XLONG_U"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE);

if (flag_psfc) { NC_fabs.push_back(&NC_psfc_fab); NC_fnames.push_back("PSFC"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); }
if (flag_msfu) { NC_fabs.push_back(&NC_msfu_fab); NC_fnames.push_back("MAPFAC_U"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); }
Expand All @@ -99,12 +125,18 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
if (flag_lmask) { NC_iabs.push_back(&NC_lmask_iab); NC_inames.push_back("LANDMASK"); NC_idim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); }

// Read the netcdf file and fill these FABs
std::string Lat_var_name = "XLAT_V";
std::string Lon_var_name = "XLONG_U";
Print() << "Building initial FABS from file " << fname << std::endl;
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, fname, NC_fnames, NC_fdim_types, NC_fabs);
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, Latitude, Longitude,
Lat_var_name, Lon_var_name,
fname, NC_fnames, NC_fdim_types, NC_fabs);

// Read the netcdf file and fill these IABs
Print() << "Building initial IABS from file " << fname << std::endl;
BuildFABsFromNetCDFFile<IArrayBox,int>(domain, fname, NC_inames, NC_idim_types, NC_iabs);
BuildFABsFromNetCDFFile<IArrayBox,int>(domain, Latitude, Longitude,
Lat_var_name, Lon_var_name,
fname, NC_inames, NC_idim_types, NC_iabs);

// TODO: FIND OUT IF WE NEED TO DIVIDE VELS BY MAPFAC
//
Expand Down
105 changes: 81 additions & 24 deletions Source/IO/ReadFromWRFInput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,43 +25,100 @@ read_from_wrfinput (int lev,
FArrayBox& NC_PHB_fab,
FArrayBox& NC_ALB_fab,
FArrayBox& NC_PB_fab,
MoistureType moisture_type)
FArrayBox& NC_LAT_fab,
FArrayBox& NC_LON_fab,
MoistureType moisture_type,
Real& Latitude,
Real& Longitude,
Geometry& geom)
{
Print() << "Loading header data from NetCDF file at level " << lev << std::endl;

if (ParallelDescriptor::IOProcessor()) {
int NC_nx, NC_ny;
Real NC_dx, NC_dy;
Real NC_epochTime;
std::string NC_dateTime;
auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
{ // Global Attributes (int)
std::vector<int> attr;
ncf.get_attr("WEST-EAST_GRID_DIMENSION", attr); NC_nx = attr[0];
ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0];
}
{ // Global Attributes (Real)
std::vector<Real> attr;
ncf.get_attr("DX", attr); NC_dx = attr[0];
ncf.get_attr("DY", attr); NC_dy = attr[0];
}
{ // Global Attributes (string)
NC_dateTime = ncf.get_attr("SIMULATION_START_DATE")+"UTC";
const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S%Z";
NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat);
}
ncf.close();

amrex::ignore_unused(NC_dateTime); amrex::ignore_unused(NC_epochTime);

// Verify the inputs geometry matches what the NETCDF file has
Real tol = 1.0e-3;
Real Len_x = NC_dx * Real(NC_nx-1);
Real Len_y = NC_dy * Real(NC_ny-1);
if (std::fabs(Len_x - (geom.ProbHi(0) - geom.ProbLo(0))) > tol) {
Print() << "X problem extent " << (geom.ProbHi(0) - geom.ProbLo(0)) << " does not match NETCDF file "
<< Len_x << "!\n";
Print() << "dx: " << NC_dx << ' ' << "Nx: " << NC_nx-1 << "\n";
Abort("Domain specification error");
}
if (std::fabs(Len_y - (geom.ProbHi(1) - geom.ProbLo(1))) > tol) {
Print() << "Y problem extent " << (geom.ProbHi(1) - geom.ProbLo(1)) << " does not match NETCDF file "
<< Len_y << "!\n";
Print() << "dy: " << NC_dy << ' ' << "Ny: " << NC_ny-1 << "\n";
Abort("Domain specification error");
}
}

Print() << "Loading initial data from NetCDF file at level " << lev << std::endl;

Vector<FArrayBox*> NC_fabs;
Vector<std::string> NC_names;
Vector<enum NC_Data_Dims_Type> NC_dim_types;

NC_fabs.push_back(&NC_xvel_fab); NC_names.push_back("U"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 0
NC_fabs.push_back(&NC_yvel_fab); NC_names.push_back("V"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 1
NC_fabs.push_back(&NC_zvel_fab); NC_names.push_back("W"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 2
NC_fabs.push_back(&NC_rho_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 3
NC_fabs.push_back(&NC_rhop_fab), NC_names.push_back("AL"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 4
NC_fabs.push_back(&NC_rhotheta_fab); NC_names.push_back("T"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 5
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_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
NC_fabs.push_back(&NC_xvel_fab); NC_names.push_back("U"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 0
NC_fabs.push_back(&NC_yvel_fab); NC_names.push_back("V"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 1
NC_fabs.push_back(&NC_zvel_fab); NC_names.push_back("W"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 2
NC_fabs.push_back(&NC_rho_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 3
NC_fabs.push_back(&NC_rhop_fab), NC_names.push_back("AL"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 4
NC_fabs.push_back(&NC_rhotheta_fab); NC_names.push_back("T"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 5
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_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
NC_fabs.push_back(&NC_LAT_fab); NC_names.push_back("XLAT_V"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 20
NC_fabs.push_back(&NC_LON_fab); NC_names.push_back("XLONG_U"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 21

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); // 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
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); // 22
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); // 23
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); // 24
}

// Read the netcdf file and fill these FABs
std::string Lat_var_name = "XLAT_V";
std::string Lon_var_name = "XLONG_U";
Print() << "Building initial FABS from file " << fname << std::endl;
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, fname, NC_names, NC_dim_types, NC_fabs);
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, Latitude, Longitude,
Lat_var_name, Lon_var_name,
fname, NC_names, NC_dim_types, NC_fabs);


//
// Convert the velocities using the map factors
Expand Down
5 changes: 4 additions & 1 deletion Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ ERF::init_from_metgrid (int lev)
Vector<FArrayBox> NC_MSFV_fab; NC_MSFV_fab.resize(ntimes);
Vector<FArrayBox> NC_MSFM_fab; NC_MSFM_fab.resize(ntimes);
Vector<FArrayBox> NC_sst_fab; NC_sst_fab.resize (ntimes);
Vector<FArrayBox> NC_LAT_fab; NC_LAT_fab.resize (ntimes);
Vector<FArrayBox> NC_LON_fab; NC_LON_fab.resize (ntimes);

// *** IArrayBox's at this level for holding mask data
Vector<IArrayBox> NC_lmask_iab; NC_lmask_iab.resize(ntimes);
Expand Down Expand Up @@ -87,7 +89,8 @@ ERF::init_from_metgrid (int lev)
NC_temp_fab[it], NC_rhum_fab[it], NC_pres_fab[it],
NC_ght_fab[it], NC_hgt_fab[it], NC_psfc_fab[it],
NC_MSFU_fab[it], NC_MSFV_fab[it], NC_MSFM_fab[it],
NC_sst_fab[it], NC_lmask_iab[it]);
NC_sst_fab[it], NC_LAT_fab[it], NC_LON_fab[it],
NC_lmask_iab[it], Latitude, Longitude, geom[lev]);
} // it

// Verify that files in nc_init_file[lev] are ordered from earliest to latest.
Expand Down
Loading

0 comments on commit 24ddedf

Please sign in to comment.