From ed07455d14617ce2fcc1d9360e14bd341ff2cf99 Mon Sep 17 00:00:00 2001 From: Hannah Klion Date: Fri, 22 Nov 2024 17:11:49 -0800 Subject: [PATCH] add attributes to netcdf file, as well as s_w --- Source/IO/NCPlotFile.cpp | 273 +++++++++++++++--- Source/Initialization/REMORA_init.cpp | 22 -- .../Initialization/REMORA_make_new_level.cpp | 8 +- Source/REMORA.H | 3 + Source/Utils/DepthStretchTransform.H | 5 + 5 files changed, 250 insertions(+), 61 deletions(-) diff --git a/Source/IO/NCPlotFile.cpp b/Source/IO/NCPlotFile.cpp index 7732a65..ff1cf15 100644 --- a/Source/IO/NCPlotFile.cpp +++ b/Source/IO/NCPlotFile.cpp @@ -93,8 +93,6 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head // Number of cells in this "domain" at this level std::vector n_cells; - int nblocks = grids[lev].size(); - // We only do single-level writes when using NetCDF format int flev = 1; //max_level; @@ -119,26 +117,19 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head n_cells.push_back(nz); int num_pts = nx * ny * nz; - int num_u_pts = (nx + 1) * (ny + 2) * nz; - int num_v_pts = (nx + 2) * (ny + 1) * nz; const std::string nt_name = "ocean_time"; const std::string ndim_name = "num_geo_dimensions"; - const std::string np_name = "num_points_per_block"; - const std::string np_u_name = "num_points_per_block_u"; - const std::string np_v_name = "num_points_per_block_v"; - - const std::string nb_name = "num_blocks"; const std::string flev_name = "FINEST_LEVEL"; const std::string nx_name = "NX"; const std::string ny_name = "NY"; const std::string nz_name = "NZ"; - const std::string nx_s_name = "xi_rho"; - const std::string ny_s_name = "eta_rho"; - const std::string nz_s_name = "s_rho"; + const std::string nx_r_name = "xi_rho"; + const std::string ny_r_name = "eta_rho"; + const std::string nz_r_name = "s_rho"; const std::string nx_u_name = "xi_u"; const std::string ny_u_name = "eta_u"; @@ -148,6 +139,7 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head const std::string nx_p_name = "xi_psi"; const std::string ny_p_name = "eta_psi"; + const std::string nz_w_name = "s_w"; const Real fill_value = 1.0e37_rt; @@ -157,9 +149,9 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head ncf.def_dim(nt_name, nt); ncf.def_dim(ndim_name, AMREX_SPACEDIM); - ncf.def_dim(nx_s_name, nx + 2); - ncf.def_dim(ny_s_name, ny + 2); - ncf.def_dim(nz_s_name, nz); + ncf.def_dim(nx_r_name, nx + 2); + ncf.def_dim(ny_r_name, ny + 2); + ncf.def_dim(nz_r_name, nz); ncf.def_dim(nx_u_name, nx + 1); ncf.def_dim(ny_u_name, ny + 2); @@ -170,11 +162,8 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head ncf.def_dim(nx_p_name, nx + 1); ncf.def_dim(ny_p_name, ny + 1); - ncf.def_dim(np_name, num_pts); - ncf.def_dim(np_u_name, num_u_pts); - ncf.def_dim(np_v_name, num_v_pts); + ncf.def_dim(nz_w_name, nz + 1); - ncf.def_dim(nb_name, nblocks); ncf.def_dim(flev_name, flev); ncf.def_dim(nx_name, n_cells[0]); @@ -182,46 +171,213 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head ncf.def_dim(nz_name, n_cells[2]); ncf.def_var("probLo", ncutils::NCDType::Real, { ndim_name }); + ncf.var("probLo").put_attr("long_name","Low side of problem domain in internal AMReX grid"); + ncf.var("probLo").put_attr("units","meter"); ncf.def_var("probHi", ncutils::NCDType::Real, { ndim_name }); + ncf.var("probHi").put_attr("long_name","High side of problem domain in internal AMReX grid"); + ncf.var("probHi").put_attr("units","meter"); ncf.def_var("Geom.smallend", NC_INT, { flev_name, ndim_name }); + ncf.var("Geom.smallend").put_attr("long_name","Low side of problem domain in index space"); ncf.def_var("Geom.bigend", NC_INT, { flev_name, ndim_name }); + ncf.var("Geom.bigend").put_attr("long_name","High side of problem domain in index space"); ncf.def_var("CellSize", ncutils::NCDType::Real, { flev_name, ndim_name }); + ncf.var("CellSize").put_attr("long_name","Cell size on internal AMReX grid"); + ncf.var("CellSize").put_attr("units","meter"); + + ncf.def_var("theta_s",ncutils::NCDType::Real,{}); + ncf.var("theta_s").put_attr("long_name","S-coordinate surface control parameter"); + ncf.def_var("theta_b",ncutils::NCDType::Real,{}); + ncf.var("theta_b").put_attr("long_name","S-coordinate bottom control parameter"); + ncf.def_var("hc",ncutils::NCDType::Real,{}); + ncf.var("hc").put_attr("long_name","S-coordinate parameter, critical depth"); + ncf.var("hc").put_attr("units","meter"); + + ncf.def_var("grid",NC_INT, {}); + ncf.var("grid").put_attr("cf_role","grid_topology"); + ncf.var("grid").put_attr("topology_dimension","2"); + ncf.var("grid").put_attr("node_dimensions", "xi_psi eta_psi"); + ncf.var("grid").put_attr("face_dimensions", "xi_rho: xi_psi (padding: both) eta_rho: eta_psi (padding: both)"); + ncf.var("grid").put_attr("edge1_dimensions", "xi_u: xi_psi eta_u: eta_psi (padding: both)"); + ncf.var("grid").put_attr("edge2_dimensions", "xi_v: xi_psi (padding: both) eta_v: eta_psi"); + ncf.var("grid").put_attr("node_coordinates", "x_psi y_psi"); + ncf.var("grid").put_attr("face_coordinates", "x_rho y_rho"); + ncf.var("grid").put_attr("edge1_coordinates", "x_u y_u"); + ncf.var("grid").put_attr("edge2_coordinates", "x_v y_v"); + ncf.var("grid").put_attr("vertical_dimensions", "s_rho: s_w (padding: none)"); + + ncf.def_var("s_rho",ncutils::NCDType::Real, {nz_r_name}); + ncf.var("s_rho").put_attr("long_name","S-coordinate at RHO-points"); + ncf.var("s_rho").put_attr("field","s_rho, scalar"); + + ncf.def_var("s_w",ncutils::NCDType::Real, {nz_w_name}); + ncf.var("s_w").put_attr("long_name","S-coordinate at W-points"); + ncf.var("s_w").put_attr("field","s_w, scalar"); + + ncf.def_var("pm",ncutils::NCDType::Real, {ny_r_name, nx_r_name}); + ncf.var("pm").put_attr("long_name","curvilinear coordinate metric in XI"); + ncf.var("pm").put_attr("units","meter-1"); + ncf.var("pm").put_attr("grid","grid"); + ncf.var("pm").put_attr("location","face"); + ncf.var("pm").put_attr("coordinates","x_rho y_rho"); + ncf.var("pm").put_attr("field","pm, scalar"); + + ncf.def_var("pn",ncutils::NCDType::Real, {ny_r_name, nx_r_name}); + ncf.var("pn").put_attr("long_name","curvilinear coordinate metric in ETA"); + ncf.var("pn").put_attr("units","meter-1"); + ncf.var("pn").put_attr("grid","grid"); + ncf.var("pn").put_attr("location","face"); + ncf.var("pn").put_attr("coordinates","x_rho y_rho"); + ncf.var("pn").put_attr("field","pn, scalar"); + + ncf.def_var("f",ncutils::NCDType::Real, {ny_r_name, nx_r_name}); + ncf.var("f").put_attr("long_name","Coriolis parameter at RHO-points"); + ncf.var("f").put_attr("units","second-1"); + ncf.var("f").put_attr("grid","grid"); + ncf.var("f").put_attr("location","face"); + ncf.var("f").put_attr("coordinates","x_rho y_rho"); + ncf.var("f").put_attr("field","coriolis, scalar"); + + ncf.def_var("x_rho",ncutils::NCDType::Real, {ny_r_name, nx_r_name}); + ncf.var("x_rho").put_attr("long_name","x-locations of RHO-points"); + ncf.var("x_rho").put_attr("units","meter"); + ncf.var("x_rho").put_attr("field","x_rho, scalar"); + + ncf.def_var("y_rho",ncutils::NCDType::Real, {ny_r_name, nx_r_name}); + ncf.var("y_rho").put_attr("long_name","y-locations of RHO-points"); + ncf.var("y_rho").put_attr("units","meter"); + ncf.var("y_rho").put_attr("field","y_rho, scalar"); - ncf.def_var("pm",ncutils::NCDType::Real, {ny_s_name, nx_s_name}); - ncf.def_var("pn",ncutils::NCDType::Real, {ny_s_name, nx_s_name}); - ncf.def_var("x_rho",ncutils::NCDType::Real, {ny_s_name, nx_s_name}); - ncf.def_var("y_rho",ncutils::NCDType::Real, {ny_s_name, nx_s_name}); ncf.def_var("x_u",ncutils::NCDType::Real, {ny_u_name, nx_u_name}); + ncf.var("x_u").put_attr("long_name","x-locations of U-points"); + ncf.var("x_u").put_attr("units","meter"); + ncf.var("x_u").put_attr("field","x_u, scalar"); + ncf.def_var("y_u",ncutils::NCDType::Real, {ny_u_name, nx_u_name}); + ncf.var("y_u").put_attr("long_name","y-locations of U-points"); + ncf.var("y_u").put_attr("units","meter"); + ncf.var("y_u").put_attr("field","y_u, scalar"); + ncf.def_var("x_v",ncutils::NCDType::Real, {ny_v_name, nx_v_name}); + ncf.var("x_v").put_attr("long_name","x-locations of V-points"); + ncf.var("x_v").put_attr("units","meter"); + ncf.var("x_v").put_attr("field","x_v, scalar"); + ncf.def_var("y_v",ncutils::NCDType::Real, {ny_v_name, nx_v_name}); + ncf.var("y_v").put_attr("long_name","y-locations of V-points"); + ncf.var("y_v").put_attr("units","meter"); + ncf.var("y_v").put_attr("field","y_v, scalar"); + ncf.def_var("x_psi",ncutils::NCDType::Real, {ny_p_name, nx_p_name}); + ncf.var("x_psi").put_attr("long_name","x-locations of PSI-points"); + ncf.var("x_psi").put_attr("units","meter"); + ncf.var("x_psi").put_attr("field","x_psi, scalar"); + ncf.def_var("y_psi",ncutils::NCDType::Real, {ny_p_name, nx_p_name}); + ncf.var("y_psi").put_attr("long_name","y-locations of PSI-points"); + ncf.var("y_psi").put_attr("units","meter"); + ncf.var("y_psi").put_attr("field","y_psi, scalar"); - ncf.def_var("x_grid", ncutils::NCDType::Real, { np_name }); - ncf.def_var("y_grid", ncutils::NCDType::Real, { np_name }); - ncf.def_var("z_grid", ncutils::NCDType::Real, { np_name }); ncf.def_var("ocean_time", ncutils::NCDType::Real, { nt_name }); + ncf.var("ocean_time").put_attr("long_name","time since initialization"); + ncf.var("ocean_time").put_attr("units","seconds since 0001-01-01 00:00:00"); + ncf.var("ocean_time").put_attr("field","time, scalar, series"); + + ncf.def_var("h", ncutils::NCDType::Real, { ny_r_name, nx_r_name }); + ncf.var("h").put_attr("long_name","bathymetry at RHO-points"); + ncf.var("h").put_attr("units","meter"); + ncf.var("h").put_attr("grid","grid"); + ncf.var("h").put_attr("location","face"); + ncf.var("h").put_attr("coordinates","x_rho y_rho"); + ncf.var("h").put_attr("field","bath, scalar"); + + ncf.def_var_fill("zeta", ncutils::NCDType::Real, { nt_name, ny_r_name, nx_r_name }, &fill_value); + ncf.var("zeta").put_attr("long_name","free-surface"); + ncf.var("zeta").put_attr("units","meter"); + ncf.var("zeta").put_attr("time","ocean_time"); + ncf.var("zeta").put_attr("grid","grid"); + ncf.var("zeta").put_attr("location","face"); + ncf.var("zeta").put_attr("coordinates","x_rho y_rho ocean_time"); + ncf.var("zeta").put_attr("field","free-surface, scalar, series"); + + ncf.def_var_fill("temp", ncutils::NCDType::Real, { nt_name, nz_r_name, ny_r_name, nx_r_name }, &fill_value); + ncf.var("temp").put_attr("long_name","potential temperature"); + ncf.var("temp").put_attr("units","Celsius"); + ncf.var("temp").put_attr("time","ocean_time"); + ncf.var("temp").put_attr("grid","grid"); + ncf.var("temp").put_attr("location","face"); + ncf.var("temp").put_attr("coordinates","x_rho y_rho s_rho ocean_time"); + ncf.var("temp").put_attr("field","temperature, scalar, series"); + + ncf.def_var_fill("salt", ncutils::NCDType::Real, { nt_name, nz_r_name, ny_r_name, nx_r_name }, &fill_value); + ncf.var("salt").put_attr("long_name","salinity"); + ncf.var("salt").put_attr("time","ocean_time"); + ncf.var("salt").put_attr("grid","grid"); + ncf.var("salt").put_attr("location","face"); + ncf.var("salt").put_attr("coordinates","x_rho y_rho s_rho ocean_time"); + ncf.var("salt").put_attr("field","salinity, scalar, series"); + + ncf.def_var_fill("u", ncutils::NCDType::Real, { nt_name, nz_r_name, ny_u_name, nx_u_name }, &fill_value); + ncf.var("u").put_attr("long_name","u-momentum component"); + ncf.var("u").put_attr("units","meter second-1"); + ncf.var("u").put_attr("time","ocean_time"); + ncf.var("u").put_attr("grid","grid"); + ncf.var("u").put_attr("location","edge1"); + ncf.var("u").put_attr("coordinates","x_u y_u s_rho ocean_time"); + ncf.var("u").put_attr("field","u-velocity, scalar, series"); + + ncf.def_var_fill("v", ncutils::NCDType::Real, { nt_name, nz_r_name, ny_v_name, nx_v_name }, &fill_value); + ncf.var("v").put_attr("long_name","v-momentum component"); + ncf.var("v").put_attr("units","meter second-1"); + ncf.var("v").put_attr("time","ocean_time"); + ncf.var("v").put_attr("grid","grid"); + ncf.var("v").put_attr("location","edge2"); + ncf.var("v").put_attr("coordinates","x_v y_v s_rho ocean_time"); + ncf.var("v").put_attr("field","v-velocity, scalar, series"); - ncf.def_var("h", ncutils::NCDType::Real, { ny_s_name, nx_s_name }); - ncf.def_var_fill("zeta", ncutils::NCDType::Real, { nt_name, ny_s_name, nx_s_name }, &fill_value); - ncf.def_var_fill("temp", ncutils::NCDType::Real, { nt_name, nz_s_name, ny_s_name, nx_s_name }, &fill_value); - ncf.def_var_fill("salt", ncutils::NCDType::Real, { nt_name, nz_s_name, ny_s_name, nx_s_name }, &fill_value); - ncf.def_var_fill("u", ncutils::NCDType::Real, { nt_name, nz_s_name, ny_u_name, nx_u_name }, &fill_value); - ncf.def_var_fill("v", ncutils::NCDType::Real, { nt_name, nz_s_name, ny_v_name, nx_v_name }, &fill_value); ncf.def_var_fill("ubar", ncutils::NCDType::Real, { nt_name, ny_u_name, nx_u_name }, &fill_value); + ncf.var("ubar").put_attr("long_name","vertically integrated u-momentum component"); + ncf.var("ubar").put_attr("units","meter second-1"); + ncf.var("ubar").put_attr("time","ocean_time"); + ncf.var("ubar").put_attr("grid","grid"); + ncf.var("ubar").put_attr("location","edge1"); + ncf.var("ubar").put_attr("coordinates","x_u y_u ocean_time"); + ncf.var("ubar").put_attr("field","ubar-velocity, scalar, series"); + ncf.def_var_fill("vbar", ncutils::NCDType::Real, { nt_name, ny_v_name, nx_v_name }, &fill_value); + ncf.var("vbar").put_attr("long_name","vertically integrated v-momentum component"); + ncf.var("vbar").put_attr("units","meter second-1"); + ncf.var("vbar").put_attr("time","ocean_time"); + ncf.var("vbar").put_attr("grid","grid"); + ncf.var("vbar").put_attr("location","edge2"); + ncf.var("vbar").put_attr("coordinates","x_v y_v ocean_time"); + ncf.var("vbar").put_attr("field","vbar-velocity, scalar, series"); + ncf.def_var("sustr", ncutils::NCDType::Real, { nt_name, ny_u_name, nx_u_name }); + ncf.var("sustr").put_attr("long_name","surface u-momentum stress"); + ncf.var("sustr").put_attr("units","newton meter-2"); + ncf.var("sustr").put_attr("time","ocean_time"); + ncf.var("sustr").put_attr("grid","grid"); + ncf.var("sustr").put_attr("location","edge1"); + ncf.var("sustr").put_attr("coordinates","x_u y_u ocean_time"); + ncf.var("sustr").put_attr("field","surface u-momentum stress, scalar, series"); + ncf.def_var("svstr", ncutils::NCDType::Real, { nt_name, ny_v_name, nx_v_name }); + ncf.var("svstr").put_attr("long_name","surface v-momentum stress"); + ncf.var("svstr").put_attr("units","newton meter-2"); + ncf.var("svstr").put_attr("time","ocean_time"); + ncf.var("svstr").put_attr("grid","grid"); + ncf.var("svstr").put_attr("location","edge2"); + ncf.var("svstr").put_attr("coordinates","x_v y_v ocean_time"); + ncf.var("svstr").put_attr("field","surface v-momentum stress, scalar, series"); Real time = 0.; // Right now this is hard-wired to {temp, salt, u, v} int n_data_items = 4; - ncf.put_attr("number_variables", std::vector { n_data_items }); +// ncf.put_attr("number_variables", std::vector { n_data_items }); ncf.put_attr("space_dimension", std::vector { AMREX_SPACEDIM }); - ncf.put_attr("current_time", std::vector { time }); +// ncf.put_attr("current_time", std::vector { time }); ncf.put_attr("start_time", std::vector { start_bdy_time }); ncf.put_attr("CurrentLevel", std::vector { flev }); ncf.put_attr("DefaultGeometry", std::vector { amrex::DefaultGeometry().Coord() }); @@ -291,6 +447,12 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head nc_CellSize.put(CellSize.data(), { static_cast(i - lev), 0 }, { 1, AMREX_SPACEDIM }); } + Real hc = solverChoice.tcline; + Real theta_s = solverChoice.theta_s; + Real theta_b = solverChoice.theta_b; + ncf.var("hc").put(&hc); + ncf.var("theta_s").put(&theta_s); + ncf.var("theta_b").put(&theta_b); } ncmpi_end_indep_data(ncf.ncid); @@ -337,6 +499,10 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head Box tmp_bx_2d(tmp_bx); tmp_bx_2d.makeSlab(2, 0); + Box tmp_bx_1d(tmp_bx); + tmp_bx_1d.makeSlab(0, 0); + tmp_bx_1d.makeSlab(1, 0); + // // These are the dimensions of the data we write for only this box // @@ -350,6 +516,33 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head long long local_start_z = static_cast(tmp_bx.smallEnd()[2]); if (write_header) { + // Only write out s_rho and s_w at x=0,y=0 to avoid NaNs + if (bx.contains(IntVect(0,0,0))) + { + { + FArrayBox tmp_srho; + tmp_srho.resize(tmp_bx_1d, 1, amrex::The_Pinned_Arena()); + + tmp_srho.template copy((*vec_s_r[lev])[mfi.index()], 0, 0, 1); + Gpu::streamSynchronize(); + + auto nc_plot_var = ncf.var("s_rho"); + //nc_plot_var.par_access(NC_INDEPENDENT); + nc_plot_var.put(tmp_srho.dataPtr(), { local_start_z }, { local_nz }); + } + { + FArrayBox tmp_sw; + tmp_sw.resize(convert(tmp_bx_1d,IntVect(0,0,1)), 1, amrex::The_Pinned_Arena()); + + tmp_sw.template copy((*vec_s_w[lev])[mfi.index()], 0, 0, 1); + Gpu::streamSynchronize(); + + auto nc_plot_var = ncf.var("s_w"); + //nc_plot_var.par_access(NC_INDEPENDENT); + nc_plot_var.put(tmp_sw.dataPtr(), { local_start_z }, { local_nz + 1}); + } + } + { FArrayBox tmp_bathy; tmp_bathy.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena()); @@ -387,6 +580,18 @@ void REMORA::WriteNCPlotFile_which(int lev, int which_subdomain, bool write_head nc_plot_var.put(tmp_pn.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx }); } + { + FArrayBox tmp_f; + tmp_f.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena()); + + tmp_f.template copy((*vec_fcor[lev])[mfi.index()], 0, 0, 1); + Gpu::streamSynchronize(); + + auto nc_plot_var = ncf.var("f"); + //nc_plot_var.par_access(NC_INDEPENDENT); + nc_plot_var.put(tmp_f.dataPtr(), { local_start_y, local_start_x }, { local_ny, local_nx }); + } + { FArrayBox tmp_xr; tmp_xr.resize(tmp_bx_2d, 1, amrex::The_Pinned_Arena()); diff --git a/Source/Initialization/REMORA_init.cpp b/Source/Initialization/REMORA_init.cpp index 44ea786..cd6cfc6 100644 --- a/Source/Initialization/REMORA_init.cpp +++ b/Source/Initialization/REMORA_init.cpp @@ -98,30 +98,8 @@ REMORA::set_zeta_average (int lev) void REMORA::set_2darrays (int lev) { - std::unique_ptr& mf_x_r = vec_x_r[lev]; - std::unique_ptr& mf_y_r = vec_y_r[lev]; auto N = Geom(lev).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ - for ( MFIter mfi(*(mf_x_r), TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - - Array4 const& x_r = (mf_x_r)->array(mfi); - Array4 const& y_r = (mf_y_r)->array(mfi); - const Box& bx = mfi.growntilebox(); - const auto & geomdata = Geom(lev).data(); - Gpu::synchronize(); - ParallelFor(amrex::makeSlab(bx,2,0), - [=] AMREX_GPU_DEVICE (int i, int j, int ) - { - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); - - x_r(i,j,0) = prob_lo[0] + (i + 0.5_rt) * dx[0]; - y_r(i,j,0) = prob_lo[1] + (j + 0.5_rt) * dx[1]; - - }); - } - vec_ubar[lev]->setVal(0.0_rt); vec_vbar[lev]->setVal(0.0_rt); diff --git a/Source/Initialization/REMORA_make_new_level.cpp b/Source/Initialization/REMORA_make_new_level.cpp index d7c26e4..1ca3c13 100644 --- a/Source/Initialization/REMORA_make_new_level.cpp +++ b/Source/Initialization/REMORA_make_new_level.cpp @@ -305,10 +305,9 @@ void REMORA::resize_stuff(int lev) vec_hOfTheConfusingName.resize(lev+1); vec_Zt_avg1.resize(lev+1); vec_s_r.resize(lev+1); + vec_s_w.resize(lev+1); vec_z_w.resize(lev+1); vec_z_r.resize(lev+1); - vec_y_r.resize(lev+1); - vec_x_r.resize(lev+1); vec_Hz.resize(lev+1); vec_Huon.resize(lev+1); vec_Hvom.resize(lev+1); @@ -427,11 +426,10 @@ void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& vec_z_phys_nd[lev].reset (new MultiFab(ba_nd,dm,1,IntVect(NGROW,NGROW,1))); // z at psi points (nodes) MIGHT NEED NGROW+1 - vec_x_r[lev].reset (new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+1,0))); // x at r points (cell center) - vec_y_r[lev].reset (new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+1,0))); // y at r points (cell center) - vec_s_r[lev].reset (new MultiFab(ba1d,dm,1,IntVect( 0, 0,0))); // scaled vertical coordinate [0,1] , transforms to z + vec_s_w[lev].reset (new MultiFab(convert(ba1d,IntVect(0,0,1)),dm,1,IntVect( 0, 0,0))); // scaled vertical coordinate at w-points [0,1] , transforms to z + vec_z_w[lev].reset (new MultiFab(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW+1,NGROW+1,0))); // z at w points (cell faces) vec_z_r[lev].reset (new MultiFab(ba,dm,1,IntVect(NGROW+1,NGROW+1,0))); // z at r points (cell center) vec_Hz[lev].reset (new MultiFab(ba,dm,1,IntVect(NGROW+1,NGROW+1,NGROW+1))); // like in ROMS, thickness of cell in z diff --git a/Source/REMORA.H b/Source/REMORA.H index c35e02d..265fc86 100644 --- a/Source/REMORA.H +++ b/Source/REMORA.H @@ -246,6 +246,9 @@ public: /** Scaled vertical coordinate (range [0,1]) that transforms to z, defined at rho points (cell centers) */ amrex::Vector> vec_s_r; + /** Scaled vertical coordinate (range [0,1]) that transforms to z, defined at w-points (cell faces) */ + amrex::Vector> vec_s_w; + /** z coordinates at psi points (cell nodes) **/ amrex::Vector> vec_z_phys_nd; diff --git a/Source/Utils/DepthStretchTransform.H b/Source/Utils/DepthStretchTransform.H index d59fc76..948e3c7 100644 --- a/Source/Utils/DepthStretchTransform.H +++ b/Source/Utils/DepthStretchTransform.H @@ -14,6 +14,7 @@ REMORA::stretch_transform (int lev) std::unique_ptr& mf_z_w = vec_z_w[lev]; std::unique_ptr& mf_z_r = vec_z_r[lev]; std::unique_ptr& mf_s_r = vec_s_r[lev]; + std::unique_ptr& mf_s_w = vec_s_w[lev]; std::unique_ptr& mf_Hz = vec_Hz[lev]; std::unique_ptr& mf_h = vec_hOfTheConfusingName[lev]; std::unique_ptr& mf_Zt_avg1 = vec_Zt_avg1[lev]; @@ -25,6 +26,7 @@ REMORA::stretch_transform (int lev) Array4 const& z_w = (mf_z_w)->array(mfi); Array4 const& z_r = (mf_z_r)->array(mfi); Array4 const& s_r = (mf_s_r)->array(mfi); + Array4 const& s_w = (mf_s_w)->array(mfi); Array4 const& Hz = (mf_Hz)->array(mfi); Array4 const& h = (mf_h)->array(mfi); Array4 const& Zt_avg1 = (mf_Zt_avg1)->array(mfi); @@ -139,6 +141,9 @@ REMORA::stretch_transform (int lev) if (i==0&&j==0&&k=0) { s_r(0,0,k) = sc_r; } + if (i==0&&j==0&&k<=N&&k>=0) { + s_w(0,0,k) = sc_w; + } cff_r=hc*sc_r; cff1_r=Cs_r;