From c87d8216c4f10e59c39d29c5627216e5667791ac Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 12 Oct 2023 13:45:38 -0700 Subject: [PATCH] WIP debugging initialization and forcing from met_em files. --- Exec/DevTests/MetGrid/GNUmakefile | 6 +- Exec/DevTests/MetGrid/inputs | 2 +- .../BoundaryConditions_metgrid.cpp | 27 +- .../Initialization/ERF_init_from_metgrid.cpp | 251 +++++++++--------- Source/TimeIntegration/ERF_fast_rhs_T.cpp | 2 +- 5 files changed, 152 insertions(+), 136 deletions(-) diff --git a/Exec/DevTests/MetGrid/GNUmakefile b/Exec/DevTests/MetGrid/GNUmakefile index a3a528fe1..3f42df89f 100644 --- a/Exec/DevTests/MetGrid/GNUmakefile +++ b/Exec/DevTests/MetGrid/GNUmakefile @@ -27,8 +27,10 @@ USE_ASSERTION = TRUE USE_NETCDF = TRUE -#USE_MOISTURE = FALSE -USE_MOISTURE = TRUE +USE_MOISTURE = FALSE +#USE_MOISTURE = TRUE + +USE_WARM_NO_PRECIP = TRUE # GNU Make Bpack := ./Make.package diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index afdd66ac7..dd885fe18 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -40,7 +40,7 @@ erf.check_int = 100 # 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_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres KE QKE +erf.plot_vars_1 = RhoQv density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres KE QKE # SOLVER CHOICE erf.alpha_T = 1.0 diff --git a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp index 3adfb40bc..65800f3c3 100644 --- a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp @@ -26,12 +26,19 @@ ERF::fill_from_metgrid (const Vector& mfs, const Real time) amrex::Real oma = 1.0 - alpha; // Flags for read vars and index mapping -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + // See IndexDefines.H +#if defined(ERF_USE_MOISTURE) + // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQt RhoQp NumVars] + Vector cons_read = {1, 1, 0, 0, 0, 1, 0}; + Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0}; +#elif defined(ERF_USE_WARM_NO_PRECIP) + // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQv RhoQc NumVars] Vector cons_read = {1, 1, 0, 0, 0, 1, 0}; Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0}; # else - Vector cons_read = {1, 1, 0, 0, 0}; // R, RT, RKE, RQKE, RS - Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0}; // R, RT, RKE, RQKE, RS + // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar NumVars] + Vector cons_read = {1, 1, 0, 0, 0}; + Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0}; #endif Vector> is_read; @@ -75,6 +82,20 @@ ERF::fill_from_metgrid (const Vector& mfs, const Real time) int ivar = ind_map[var_idx][comp_idx]; IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; + if (ivar == MetGridBdyVars::U) { + amrex::Print() << "fill_from_metgrid U var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; + } else if (ivar == MetGridBdyVars::V) { + amrex::Print() << "fill_from_metgrid V var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; + } else if (ivar == MetGridBdyVars::R) { + amrex::Print() << "fill_from_metgrid R var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; + } else if (ivar == MetGridBdyVars::T) { + amrex::Print() << "fill_from_metgrid T var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; + } else if (ivar == MetGridBdyVars::QV) { + amrex::Print() << "fill_from_metgrid QV var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; + } else { + amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl; + } + // We have data at fixed time intervals we will call dT // Then to interpolate, given time, we can define n = (time/dT) // and alpha = (time - n*dT) / dT, then we define the data at time diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 44fd06ad4..a76806707 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -51,7 +51,7 @@ init_state_from_metgrid(const Real l_rdOcp, const Vector& NC_rhum_fab, Vector& theta_fab, Vector& mxrat_fab, - Vector>& fabs_for_bcs); + amrex::Vector>& fabs_for_bcs); void init_msfs_from_metgrid(FArrayBox& msfu_fab, @@ -75,7 +75,7 @@ init_base_state_from_metgrid(const Real l_rdOcp, FArrayBox& z_phys_nd_fab, const Vector& NC_ght_fab, const Vector& NC_psfc_fab, - Vector>& fabs_for_bcs); + amrex::Vector>& fabs_for_bcs); void rh_to_mxrat(int i, int j, int k, @@ -258,7 +258,8 @@ ERF::init_from_metgrid(int lev) #else int MetGridBdyEnd = MetGridBdyVars::NumTypes-1; #endif - Vector> fabs_for_bcs; + //amrex::Vector > fabs_for_bcs; + amrex::Vector> fabs_for_bcs; fabs_for_bcs.resize(ntimes); for (int it(0); it < ntimes; it++) { fabs_for_bcs[it].resize(MetGridBdyEnd); @@ -330,7 +331,6 @@ ERF::init_from_metgrid(int lev) fabs_for_bcs); } // mf - // Set up boxes for lateral boundary arrays. bdy_data_xlo.resize(ntimes); bdy_data_xhi.resize(ntimes); @@ -433,129 +433,113 @@ ERF::init_from_metgrid(int lev) // we only need the whole domain processed at initialization and the lateral boundaries // at subsequent times. We can optimize this later if needed. For now, we need to fill // the lateral boundary arrays using the info set aside earlier. - if (amrex::ParallelDescriptor::IOProcessor()) { - - for (int it(0); it < ntimes; it++) { + for (int it(0); it < ntimes; it++) { - for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { - - auto xlo_arr = bdy_data_xlo[it][ivar].array(); - auto xhi_arr = bdy_data_xhi[it][ivar].array(); - auto ylo_arr = bdy_data_ylo[it][ivar].array(); - auto yhi_arr = bdy_data_yhi[it][ivar].array(); - const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); - - if (ivar == MetGridBdyVars::U) { - // xvel at west boundary - amrex::ParallelFor(xlo_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // xvel at east boundary - amrex::ParallelFor(xhi_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // xvel at south boundary - amrex::ParallelFor(ylo_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // xvel at north boundary - amrex::ParallelFor(yhi_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - } // MetGridBdyVars::U - - if (ivar == MetGridBdyVars::V) { - // yvel at west boundary - amrex::ParallelFor(xlo_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // yvel at east boundary - amrex::ParallelFor(xhi_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // yvel at south boundary - amrex::ParallelFor(ylo_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // yvel at north boundary - amrex::ParallelFor(yhi_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - } // MetGridBdyVars::V - - if (ivar == MetGridBdyVars::R) { - // density at west boundary - amrex::ParallelFor(xlo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // density at east boundary - amrex::ParallelFor(xhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // density at south boundary - amrex::ParallelFor(ylo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // density at north boundary - amrex::ParallelFor(yhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - } // MetGridBdyVars::R - - if (ivar == MetGridBdyVars::T) { - // potential temperature at west boundary - amrex::ParallelFor(xlo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // potential temperature at east boundary - amrex::ParallelFor(xhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // potential temperature at south boundary - amrex::ParallelFor(ylo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // potential temperature at north boundary - amrex::ParallelFor(yhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - } // MetGridBdyVars::T - - if (ivar == MetGridBdyVars::QV) { - // mixing ratio at west boundary - amrex::ParallelFor(xlo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // mixing ratio at east boundary - amrex::ParallelFor(xhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // mixing ratio at south boundary - amrex::ParallelFor(ylo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - // mixing ratio at north boundary - amrex::ParallelFor(yhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - } // MetGridBdyVars::QV - - } // ivar + for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { - } // it + auto xlo_arr = bdy_data_xlo[it][ivar].array(); + auto xhi_arr = bdy_data_xhi[it][ivar].array(); + auto ylo_arr = bdy_data_ylo[it][ivar].array(); + auto yhi_arr = bdy_data_yhi[it][ivar].array(); + const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); - } // ioproc + if (ivar == MetGridBdyVars::U) { + // xvel at west boundary + amrex::ParallelFor(xlo_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // xvel at east boundary + amrex::ParallelFor(xhi_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // xvel at south boundary + amrex::ParallelFor(ylo_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // xvel at north boundary + amrex::ParallelFor(yhi_plane_x_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + } // MetGridBdyVars::U + + if (ivar == MetGridBdyVars::V) { + // yvel at west boundary + amrex::ParallelFor(xlo_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // yvel at east boundary + amrex::ParallelFor(xhi_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // yvel at south boundary + amrex::ParallelFor(ylo_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // yvel at north boundary + amrex::ParallelFor(yhi_plane_y_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + } // MetGridBdyVars::V + + if (ivar == MetGridBdyVars::R) { + // density at west boundary + amrex::ParallelFor(xlo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // density at east boundary + amrex::ParallelFor(xhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // density at south boundary + amrex::ParallelFor(ylo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // density at north boundary + amrex::ParallelFor(yhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + } // MetGridBdyVars::R + + if (ivar == MetGridBdyVars::T) { + // potential temperature at west boundary + amrex::ParallelFor(xlo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // potential temperature at east boundary + amrex::ParallelFor(xhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // potential temperature at south boundary + amrex::ParallelFor(ylo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // potential temperature at north boundary + amrex::ParallelFor(yhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + } // MetGridBdyVars::T + + if (ivar == MetGridBdyVars::QV) { + // mixing ratio at west boundary + amrex::ParallelFor(xlo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // mixing ratio at east boundary + amrex::ParallelFor(xhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // mixing ratio at south boundary + amrex::ParallelFor(ylo_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + // mixing ratio at north boundary + amrex::ParallelFor(yhi_plane_no_stag, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); + }); + } // MetGridBdyVars::QV - // Broadcast boundary data from the IO processor. - amrex::ParallelDescriptor::Barrier(); - int ioproc = ParallelDescriptor::IOProcessorNumber(); - for (int it(0); it < ntimes; it++) { - for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { - ParallelDescriptor::Bcast(bdy_data_xlo[it][ivar].dataPtr(),bdy_data_xlo[it][ivar].box().numPts(),ioproc); - ParallelDescriptor::Bcast(bdy_data_xhi[it][ivar].dataPtr(),bdy_data_xhi[it][ivar].box().numPts(),ioproc); - ParallelDescriptor::Bcast(bdy_data_ylo[it][ivar].dataPtr(),bdy_data_ylo[it][ivar].box().numPts(),ioproc); - ParallelDescriptor::Bcast(bdy_data_yhi[it][ivar].dataPtr(),bdy_data_yhi[it][ivar].box().numPts(),ioproc); } // ivar + } // it } @@ -649,7 +633,7 @@ init_state_from_metgrid(const Real l_rdOcp, const Vector& NC_pres_fab, Vector& theta_fab, Vector& mxrat_fab, - Vector>& fabs_for_bcs) + amrex::Vector>& fabs_for_bcs) { int ntimes = NC_hgt_fab.size(); for (int it = 0; it < ntimes; it++) @@ -844,6 +828,14 @@ init_state_from_metgrid(const Real l_rdOcp, #endif } + // TODO: TEMPORARY CODE TO RUN QUIESCENT, REMOVE WHEN NOT NEEDED. +// if (it == 0) { +// x_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. +// y_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. +// } +// fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. +// fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. + } // it } @@ -976,8 +968,8 @@ calc_rho_p(const int& kmax, 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; - if (Pm_vec[k] <= 0.0) Pm_vec[k] = 0.0; - Rhom_vec[k] = 1.0/(R_d/p_0*theta_v[k]*qvf*std::pow(Pm_vec[k]/p_0, -iGamma)); + if (Pm_vec[k] < 0.0) Pm_vec[k] = 0.0; + Rhom_vec[k] = 1.0/(R_d/p_0*theta_v[k]*std::pow(Pm_vec[k]/p_0, -iGamma)); } // it } // k @@ -989,7 +981,7 @@ calc_rho_p(const int& kmax, 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; + if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; Rhod_vec[k] = 1.0/(R_d/p_0*Theta_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); } // it } // k @@ -1086,7 +1078,7 @@ init_base_state_from_metgrid(const Real l_rdOcp, FArrayBox& z_phys_nd_fab, const Vector& NC_ght_fab, const Vector& NC_psfc_fab, - Vector>& fabs_for_bcs) + amrex::Vector>& fabs_for_bcs) { #if defined(ERF_USE_MOISTURE) int RhoQ_comp = RhoQt_comp; @@ -1145,8 +1137,8 @@ init_base_state_from_metgrid(const Real l_rdOcp, Real RhoQ = r_hse_arr(i,j,k)*new_data(i,j,k,RhoQ_comp); new_data(i,j,k,RhoQ_comp) = RhoQ; #endif - 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); +// 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); }); } @@ -1198,6 +1190,7 @@ init_base_state_from_metgrid(const Real l_rdOcp, z_vec,Rhod_vec,Pd_vec); for (int k=0; k < kmax; k++) { r_hse_arr(i,j,k) = Rhod_vec[k]; + p_hse_arr(i,j,k) = Pd_vec[k]; #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) Q_arr(i,j,k) = Rhod_vec[k]*Q_vec[k]; #endif diff --git a/Source/TimeIntegration/ERF_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index 4a1b61822..06fd50df7 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -200,7 +200,7 @@ void erf_fast_rhs_T (int step, int /*level*/, const Array4 & stage_xmom = S_stage_data[IntVar::xmom].const_array(mfi); const Array4 & stage_ymom = S_stage_data[IntVar::ymom].const_array(mfi); -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_MOISTURE) +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) const Array4 & prim = S_stage_prim.const_array(mfi); #endif