diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index b04b2c756..9e18901e1 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -136,6 +136,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/TimeIntegration/ERF_TimeStep.cpp ${SRC_DIR}/TimeIntegration/ERF_advance_dycore.cpp ${SRC_DIR}/TimeIntegration/ERF_advance_microphysics.cpp + ${SRC_DIR}/TimeIntegration/ERF_advance_lsm.cpp ${SRC_DIR}/TimeIntegration/ERF_advance_radiation.cpp ${SRC_DIR}/TimeIntegration/ERF_make_buoyancy.cpp ${SRC_DIR}/TimeIntegration/ERF_make_fast_coeffs.cpp @@ -164,6 +165,8 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Microphysics/FastEddy/FastEddy.cpp ${SRC_DIR}/Microphysics/FastEddy/Diagnose_FE.cpp ${SRC_DIR}/Microphysics/FastEddy/Update_FE.cpp + ${SRC_DIR}/LandSurfaceModel/SLM/SLM.cpp + ${SRC_DIR}/LandSurfaceModel/MM5/MM5.cpp ) if(NOT "${erf_exe_name}" STREQUAL "erf_unit_tests") @@ -208,8 +211,12 @@ function(build_erf_lib erf_lib_name) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler) - target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/FastEddy) - + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/FastEddy) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/LandSurfaceModel) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/LandSurfaceModel/Null) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/LandSurfaceModel/SLM) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/LandSurfaceModel/MM5) + if(ERF_ENABLE_RRTMGP) target_link_libraries(${erf_lib_name} PUBLIC yakl) target_link_libraries(${erf_lib_name} PUBLIC rrtmgp) diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index 13da633d3..8b0fca2e8 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -27,5 +27,6 @@ else () add_subdirectory(DevTests/ParticlesOverWoA) add_subdirectory(DevTests/MiguelDev) add_subdirectory(DevTests/MetGrid) + add_subdirectory(DevTests/LandSurfaceModel) add_subdirectory(DevTests/TemperatureSource) endif() diff --git a/Exec/DevTests/LandSurfaceModel/CMakeLists.txt b/Exec/DevTests/LandSurfaceModel/CMakeLists.txt new file mode 100644 index 000000000..ca744373e --- /dev/null +++ b/Exec/DevTests/LandSurfaceModel/CMakeLists.txt @@ -0,0 +1,15 @@ +set(erf_exe_name LandSurfaceModel) + +add_executable(${erf_exe_name} "") +target_sources(${erf_exe_name} + PRIVATE + prob.cpp +) + +target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + +include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake) +build_erf_exe(${erf_exe_name}) + +#find_package(AMReX) +#target_link_libraries( ${_target} AMReX::amrex) diff --git a/Exec/DevTests/LandSurfaceModel/GNUmakefile b/Exec/DevTests/LandSurfaceModel/GNUmakefile new file mode 100644 index 000000000..76fc1049d --- /dev/null +++ b/Exec/DevTests/LandSurfaceModel/GNUmakefile @@ -0,0 +1,35 @@ +# AMReX +COMP = gnu +PRECISION = DOUBLE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = FALSE +COMM_PROFILE = FALSE +TRACE_PROFILE = FALSE +MEM_PROFILE = FALSE +USE_GPROF = FALSE + +# Performance +USE_MPI = TRUE +USE_OMP = FALSE + +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +# Debugging +DEBUG = FALSE + +TEST = TRUE +USE_ASSERTION = TRUE + +USE_NETCDF = TRUE + +# GNU Make +Bpack := ./Make.package +Blocs := . + +ERF_HOME := ../../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/DevTests/LandSurfaceModel +include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/DevTests/LandSurfaceModel/Make.package b/Exec/DevTests/LandSurfaceModel/Make.package new file mode 100644 index 000000000..aeacb72f0 --- /dev/null +++ b/Exec/DevTests/LandSurfaceModel/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += prob.H +CEXE_sources += prob.cpp diff --git a/Exec/DevTests/LandSurfaceModel/README b/Exec/DevTests/LandSurfaceModel/README new file mode 100644 index 000000000..499fab443 --- /dev/null +++ b/Exec/DevTests/LandSurfaceModel/README @@ -0,0 +1 @@ +This problem setup is to test the ability of ERF to handle a land surface model. diff --git a/Exec/DevTests/LandSurfaceModel/inputs b/Exec/DevTests/LandSurfaceModel/inputs new file mode 100644 index 000000000..f138b541d --- /dev/null +++ b/Exec/DevTests/LandSurfaceModel/inputs @@ -0,0 +1,74 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 100 + +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 2430000 2673000 12000 +amr.n_cell = 90 99 50 + +geometry.is_periodic = 0 0 0 + +xlo.type = "Outflow" +xhi.type = "Outflow" + +ylo.type = "Outflow" +yhi.type = "Outflow" + +zlo.type = "Most" +erf.most.z0 = 0.1 +erf.most.zref = 120.0 +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.01 +erf.fixed_mri_dt_ratio = 4 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = -1 # timesteps between computing mass +#erf.data_log = my_data my_1d_data +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints +erf.restart_type = "native" +# 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 rhotheta rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens + +erf.plot_lsm = true + +# SOLVER CHOICE +erf.alpha_T = 1.0 +erf.alpha_C = 1.0 +erf.use_gravity = true + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 +#erf.pbl_type = "MYNN2.5" +#erf.QKE_0 = 0.5 + +erf.use_terrain = true +erf.terrain_smoothing = 2 + +erf.moisture_model = "NullMoist" + +erf.land_surface_model = "MM5" + +# INITIALIZATION WITH METGRID DATA +erf.metgrid_bdy_width = 5 +erf.metgrid_bdy_set_width = 1 +erf.init_type = "metgrid" +erf.nc_init_file_0 = "met_em.d01.2016-10-06_00_00_00.nc" "met_em.d01.2016-10-06_06_00_00.nc" + +# NO CACHE TILING +fabarray.mfiter_tile_size = 1024 1024 1024 diff --git a/Exec/DevTests/LandSurfaceModel/prob.H b/Exec/DevTests/LandSurfaceModel/prob.H new file mode 100644 index 000000000..92ead6942 --- /dev/null +++ b/Exec/DevTests/LandSurfaceModel/prob.H @@ -0,0 +1,24 @@ +#ifndef _PROB_H_ +#define _PROB_H_ + +#include + +#include "AMReX_REAL.H" + +#include "prob_common.H" + +struct ProbParm : ProbParmDefaults { +}; // namespace ProbParm +class Problem : public ProblemBase +{ +public: + Problem(const amrex::Real* /*problo*/, const amrex::Real* /*probhi*/); + +protected: + std::string name() override { return "LandSurfaceModel"; } + +private: + ProbParm parms; +}; + +#endif diff --git a/Exec/DevTests/LandSurfaceModel/prob.cpp b/Exec/DevTests/LandSurfaceModel/prob.cpp new file mode 100644 index 000000000..c01e06833 --- /dev/null +++ b/Exec/DevTests/LandSurfaceModel/prob.cpp @@ -0,0 +1,12 @@ +#include "prob.H" + +using namespace amrex; + +std::unique_ptr +amrex_probinit(const amrex_real* problo, const amrex_real* probhi) +{ + return std::make_unique(problo, probhi); +} + +Problem::Problem(const amrex::Real* /*problo*/, const amrex::Real* /*probhi*/) +{} diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 787166d3c..6b2a09dd1 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -127,6 +127,26 @@ include $(ERF_MOISTURE_KESSLER_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) +ERF_LSM_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel +include $(ERF_LSM_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_LSM_DIR) +INCLUDE_LOCATIONS += $(ERF_LSM_DIR) + +ERF_LSM_NULL_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel/Null +include $(ERF_LSM_NULL_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_LSM_NULL_DIR) +INCLUDE_LOCATIONS += $(ERF_LSM_NULL_DIR) + +ERF_LSM_SLM_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel/SLM +include $(ERF_LSM_SLM_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_LSM_SLM_DIR) +INCLUDE_LOCATIONS += $(ERF_LSM_SLM_DIR) + +ERF_LSM_MM5_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel/MM5 +include $(ERF_LSM_MM5_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_LSM_MM5_DIR) +INCLUDE_LOCATIONS += $(ERF_LSM_MM5_DIR) + ifeq ($(COMPUTE_ERROR), TRUE) DEFINES += -DERF_COMPUTE_ERROR endif diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 2b805a060..c0cdeb802 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -35,6 +35,8 @@ public: amrex::Vector>& z_phys_nd, amrex::Vector>>& sst_lev, amrex::Vector>>& lmask_lev, + amrex::Vector> lsm_data, + amrex::Vector> lsm_flux, amrex::Real start_bdy_time = 0.0, amrex::Real bdy_time_interval = 0.0) : m_start_bdy_time(start_bdy_time), @@ -113,6 +115,19 @@ public: } } + // Get pointers to LSM data and Fluxes + m_lsm_data_lev.resize(nlevs); + m_lsm_flux_lev.resize(nlevs); + for (int lev(0); lev(ba2d,dm,ncomp,ng); - if (m_sst_lev[lev][0]) { // Valid SST data at t==0 + // TODO: Do we want lsm_data to have theta at 0 index always? + // Do we want an external enum struct for indexing? + if (m_sst_lev[lev][0] || m_lsm_data_lev[lev][0]) { + // Valid SST or LSM data; t_surf set before computing fluxes (avoids extended lambda capture) theta_type = ThetaCalcType::SURFACE_TEMPERATURE; - amrex::MultiFab::Copy(*(t_surf[lev]), *(m_sst_lev[lev][0]), 0, 0, 1, ng); - } else if (erf_st) { // Constant temp + } else if (erf_st) { t_surf[lev]->setVal(surf_temp); } else { t_surf[lev]->setVal(0.0); @@ -187,8 +204,11 @@ public: const FluxCalc& flux_comp); void - time_interp_tsurf(const int& lev, - const amrex::Real& time); + time_interp_sst (const int& lev, + const amrex::Real& time); + + void + get_lsm_tsurf (const int& lev); void update_surf_temp (const amrex::Real& time) @@ -265,6 +285,8 @@ private: amrex::Vector> m_sst_lev; amrex::Vector> m_lmask_lev; + amrex::Vector> m_lsm_data_lev; + amrex::Vector> m_lsm_flux_lev; }; #endif /* ABLMOST_H */ diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index cedb9202e..2bb8778aa 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -15,7 +15,14 @@ ABLMost::update_fluxes (const int& lev, int max_iters) { // Update SST data if we have a valid pointer - if (m_sst_lev[lev][0]) time_interp_tsurf(lev, time); + if (m_sst_lev[lev][0]) time_interp_sst(lev, time); + + // TODO: we want 0 index to always be theta? + // Update land surface temp if we have a valid pointer + if (m_lsm_data_lev[lev][0]) get_lsm_tsurf(lev); + + // Fill interior ghost cells + t_surf[lev]->FillBoundary(m_geom[lev].periodicity()); // Compute plane averages for all vars (regardless of flux type) m_ma.compute_averages(lev); @@ -145,6 +152,10 @@ ABLMost::compute_most_bcs (const int& lev, const int icomp = 0; for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi) { + // TODO: No LSM lateral ghost cells, should this change? + // Valid CC box + Box vbx = mfi.validbox(); vbx.makeSlab(2,zlo-1); + // Get field arrays const auto cons_arr = mfs[Vars::cons]->array(mfi); const auto velx_arr = mfs[Vars::xvel]->array(mfi); @@ -168,6 +179,12 @@ ABLMost::compute_most_bcs (const int& lev, const auto t_star_arr = t_star[lev]->array(mfi); const auto t_surf_arr = t_surf[lev]->array(mfi); + // Get LSM fluxes + auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) : + Array4 {}; + auto lsm_flux_arr = (m_lsm_flux_lev[lev][0]) ? m_lsm_flux_lev[lev][0]->array(mfi) : + Array4 {}; + for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx) { const Box& bx = (*mfs[var_idx])[mfi].box(); @@ -181,10 +198,15 @@ ABLMost::compute_most_bcs (const int& lev, ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain; - flux_comp.compute_t_flux(i, j, k, n, icomp, dz, - cons_arr, velx_arr, vely_arr, eta_arr, - umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr, - dest_arr); + Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, + cons_arr, velx_arr, vely_arr, eta_arr, + umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr, + dest_arr); + + int is_land = (lmask_arr) ? lmask_arr(i,j,zlo) : 1; + if (is_land && vbx.contains(i,j,k)) { + lsm_flux_arr(i,j,zlo) = Tflux; + } }); } else if (var_idx == Vars::xvel || var_idx == Vars::xmom) { @@ -220,8 +242,8 @@ ABLMost::compute_most_bcs (const int& lev, } void -ABLMost::time_interp_tsurf(const int& lev, - const Real& time) +ABLMost::time_interp_sst (const int& lev, + const Real& time) { // Time interpolation Real dT = m_bdy_time_interval; @@ -240,11 +262,49 @@ ABLMost::time_interp_tsurf(const int& lev, auto t_surf_arr = t_surf[lev]->array(mfi); const auto sst_hi_arr = m_sst_lev[lev][n_time+1]->const_array(mfi); const auto sst_lo_arr = m_sst_lev[lev][n_time ]->const_array(mfi); + auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) : + Array4 {}; + + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1; + if (!is_land) { + t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k) + + alpha * sst_hi_arr(i,j,k); + } + }); + } +} + +void +ABLMost::get_lsm_tsurf (const int& lev) +{ + for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi) + { + Box gtbx = mfi.growntilebox(); + + // TODO: LSM does not carry lateral ghost cells. + // This copies the valid box into the ghost cells. + // Fillboundary is called after this to pick up the + // interior ghost and periodic directions. Is there + // a better approach? + Box vbx = mfi.validbox(); + int i_lo = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0); + int j_lo = vbx.smallEnd(1); int j_hi = vbx.bigEnd(1); + + auto t_surf_arr = t_surf[lev]->array(mfi); + auto lmask_arr = (m_lmask_lev[lev][0]) ? m_lmask_lev[lev][0]->array(mfi) : + Array4 {}; + const auto lsm_arr = m_lsm_data_lev[lev][0]->const_array(mfi); ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k) - + alpha * sst_hi_arr(i,j,k); + int is_land = (lmask_arr) ? lmask_arr(i,j,k) : 1; + if (is_land) { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + t_surf_arr(i,j,k) = lsm_arr(li,lj,k); + } }); } } diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index c654546f0..59a3ef4e1 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -593,7 +593,7 @@ struct moeng_flux AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void + amrex::Real compute_t_flux (const int& i, const int& j, const int& k, @@ -655,6 +655,8 @@ struct moeng_flux amrex::Real deltaz = dz * (zlo - k); dest_arr(i,j,k,icomp+n) = rho*(theta - moflux*rho/eta*deltaz); + + return moflux; } AMREX_GPU_DEVICE @@ -802,7 +804,7 @@ struct donelan_flux AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void + amrex::Real compute_t_flux (const int& i, const int& j, const int& k, @@ -846,6 +848,8 @@ struct donelan_flux amrex::Real deltaz = dz * (zlo - k); dest_arr(i,j,k,icomp+n) = rho*(theta - moflux*rho/eta*deltaz); + + return moflux; } AMREX_GPU_DEVICE diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index abf5e6b6a..bc21a6244 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -32,6 +32,10 @@ enum struct MoistureType { Kessler, SAM, FastEddy, None }; +enum struct LandSurfaceType { + SLM, MM5, None +}; + enum struct Coord { x, y, z }; @@ -93,7 +97,7 @@ struct SolverChoice { amrex::Abort("buoyancy_type must be 1, 2, 3 or 4"); } - /* + // Set a different default for moist vs dry if (moisture_type != MoistureType::None) { if (moisture_type == MoistureType::FastEddy) { @@ -102,7 +106,17 @@ struct SolverChoice { buoyancy_type = 2; // uses Tprime } } - */ + + // What type of land surface model to use + static std::string lsm_model_string = "None"; + pp.query("land_surface_model", lsm_model_string); + if (lsm_model_string == "SLM") { + lsm_type = LandSurfaceType::SLM; + } else if (lsm_model_string == "MM5") { + lsm_type = LandSurfaceType::MM5; + } else { + lsm_type = LandSurfaceType::None; + } // Is the terrain static or moving? static std::string terrain_type_string = "Static"; @@ -372,6 +386,7 @@ struct SolverChoice { CouplingType coupling_type; TerrainType terrain_type; MoistureType moisture_type; + LandSurfaceType lsm_type; ABLDriverType abl_driver_type; amrex::GpuArray abl_pressure_grad; diff --git a/Source/ERF.H b/Source/ERF.H index 5152d2f75..c7cf09297 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -45,6 +45,7 @@ #endif #include "Microphysics.H" +#include "LandSurface.H" #ifdef ERF_USE_RRTMGP #include "Radiation.H" @@ -213,9 +214,14 @@ public: amrex::Real dt, amrex::Real time, amrex::InterpFaceRegister* ifr); - void advance_microphysics (int lev, amrex::MultiFab& cons_in, + void advance_microphysics (int lev, + amrex::MultiFab& cons_in, const amrex::Real& dt_advance); + void advance_lsm (int lev, + amrex::MultiFab& /*cons_in*/, + const amrex::Real& dt_advance); + #if defined(ERF_USE_RRTMGP) void advance_radiation (int lev, amrex::MultiFab& cons_in, @@ -527,6 +533,10 @@ private: Microphysics micro; amrex::Vector> qmoist; // (lev,ncomp) This has up to 6 components: qv, qc, qi, qr, qs, qg + LandSurface lsm; + amrex::Vector> lsm_data; // (lev,ncomp) Components: theta, q1, q2 + amrex::Vector> lsm_flux; // (lev,ncomp) Components: theta, q1, q2 + #if defined(ERF_USE_RRTMGP) Radiation rad; amrex::Vector qheating_rates; // radiation heating rate source terms @@ -634,6 +644,7 @@ private: std::string plot_file_2 {"plt_2_"}; int plot_int_1 = -1; int plot_int_2 = -1; + bool plot_lsm = false; // other sampling output control int profile_int = -1; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index b8dfd971a..40a6895ee 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -122,6 +122,11 @@ ERF::ERF () qheating_rates.resize(nlevs_max); #endif + // NOTE: size lsm before readparams (chooses the model at all levels) + lsm.ReSize(nlevs_max); + lsm_data.resize(nlevs_max); + lsm_flux.resize(nlevs_max); + ReadParameters(); const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1); const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2); @@ -617,7 +622,7 @@ ERF::InitData () if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { m_most = std::make_unique(geom, vars_old, Theta_prim, z_phys_nd, - sst_lev, lmask_lev + sst_lev, lmask_lev, lsm_data, lsm_flux #ifdef ERF_USE_NETCDF ,start_bdy_time, bdy_time_interval #endif @@ -1103,6 +1108,8 @@ ERF::ReadParameters () pp.query("profile_int", profile_int); + pp.query("plot_lsm", plot_lsm); + pp.query("output_1d_column", output_1d_column); pp.query("column_per", column_per); pp.query("column_interval", column_interval); @@ -1174,6 +1181,21 @@ ERF::ReadParameters () amrex::Abort("Dont know this moisture_type!") ; } + // What type of land surface model to use + // NOTE: Must be checked after init_params + if (solverChoice.lsm_type == LandSurfaceType::SLM) { + lsm.SetModel(); + amrex::Print() << "SLM land surface model!\n"; + } else if (solverChoice.lsm_type == LandSurfaceType::MM5) { + lsm.SetModel(); + amrex::Print() << "MM5 land surface model!\n"; + } else if (solverChoice.lsm_type == LandSurfaceType::None) { + lsm.SetModel(); + amrex::Print() << "Null land surface model!\n"; + } else { + amrex::Abort("Dont know this moisture_type!") ; + } + if (verbose > 0) { solverChoice.display(); } @@ -1514,6 +1536,11 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, qheating_rates.resize(nlevs_max); #endif + // NOTE: size micro before readparams (chooses the model at all levels) + lsm.ReSize(nlevs_max); + lsm_data.resize(nlevs_max); + lsm_flux.resize(nlevs_max); + ReadParameters(); const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1); const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 9ae9d1f1a..a7d50bf98 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -73,6 +73,22 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar); } + //******************************************************************************************** + // Land Surface Model + // ******************************************************************************************* + int lsm_size = lsm.Get_Data_Size(); + lsm_data[lev].resize(lsm_size); + lsm_flux[lev].resize(lsm_size); + lsm.Define(lev, solverChoice); + if (solverChoice.lsm_type != LandSurfaceType::None) + { + lsm.Init(lev, vars_new[lev][Vars::cons], Geom(lev), 0.0); // dummy dt value + } + for (int mvar(0); mvarnGrowVect(); - MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); - MultiFab::Copy(z_height,*z_phys_nd[lev],0,0,1,ng); - VisMF::Write(z_height, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Z_Phys_nd")); - } + if (solverChoice.use_terrain) { + // Note that we also write the ghost cells of z_phys_nd + ng = z_phys_nd[lev]->nGrowVect(); + MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); + MultiFab::Copy(z_height,*z_phys_nd[lev],0,0,1,ng); + VisMF::Write(z_height, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Z_Phys_nd")); + } - // Note that we also write the ghost cells of the mapfactors (2D) - BoxList bl2d = grids[lev].boxList(); - for (auto& b : bl2d) { - b.setRange(2,0); - } - BoxArray ba2d(std::move(bl2d)); - - ng = mapfac_m[lev]->nGrowVect(); - MultiFab mf_m(ba2d,dmap[lev],1,ng); - MultiFab::Copy(mf_m,*mapfac_m[lev],0,0,1,ng); - VisMF::Write(mf_m, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_m")); - - ng = mapfac_u[lev]->nGrowVect(); - MultiFab mf_u(convert(ba2d,IntVect(1,0,0)),dmap[lev],1,ng); - MultiFab::Copy(mf_u,*mapfac_u[lev],0,0,1,ng); - VisMF::Write(mf_u, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_u")); - - ng = mapfac_v[lev]->nGrowVect(); - MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ng); - MultiFab::Copy(mf_v,*mapfac_v[lev],0,0,1,ng); - VisMF::Write(mf_v, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_v")); - } + if (solverChoice.lsm_type != LandSurfaceType::None) { + for (int mvar(0); mvarboxArray(); + DistributionMapping dm = lsm_data[lev][mvar]->DistributionMap(); + ng = lsm_data[lev][mvar]->nGrowVect(); + int nvar = lsm_data[lev][mvar]->nComp(); + MultiFab lsm_vars(ba,dm,nvar,ng); + MultiFab::Copy(lsm_vars,*(lsm_data[lev][mvar]),0,0,nvar,ng); + VisMF::Write(lsm_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "LsmVars")); + } + } + + // Note that we also write the ghost cells of the mapfactors (2D) + BoxList bl2d = grids[lev].boxList(); + for (auto& b : bl2d) { + b.setRange(2,0); + } + BoxArray ba2d(std::move(bl2d)); + + ng = mapfac_m[lev]->nGrowVect(); + MultiFab mf_m(ba2d,dmap[lev],1,ng); + MultiFab::Copy(mf_m,*mapfac_m[lev],0,0,1,ng); + VisMF::Write(mf_m, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_m")); + + ng = mapfac_u[lev]->nGrowVect(); + MultiFab mf_u(convert(ba2d,IntVect(1,0,0)),dmap[lev],1,ng); + MultiFab::Copy(mf_u,*mapfac_u[lev],0,0,1,ng); + VisMF::Write(mf_u, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_u")); + + ng = mapfac_v[lev]->nGrowVect(); + MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ng); + MultiFab::Copy(mf_v,*mapfac_v[lev],0,0,1,ng); + VisMF::Write(mf_v, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_v")); + } #ifdef ERF_USE_PARTICLES particleData.Checkpoint(checkpointname); @@ -213,9 +225,7 @@ ERF::WriteCheckpointFile () const void ERF::ReadCheckpointFile () { - amrex::Print() << "Restart from checkpoint " << restart_chkfile << "\n"; - - + amrex::Print() << "Restart from native checkpoint " << restart_chkfile << "\n"; // Header std::string File(restart_chkfile + "/Header"); @@ -343,6 +353,18 @@ ERF::ReadCheckpointFile () update_terrain_arrays(lev, t_new[lev]); } + if (solverChoice.lsm_type != LandSurfaceType::None) { + for (int mvar(0); mvarboxArray(); + DistributionMapping dm = lsm_data[lev][mvar]->DistributionMap(); + ng = lsm_data[lev][mvar]->nGrowVect(); + int nvar = lsm_data[lev][mvar]->nComp(); + MultiFab lsm_vars(ba,dm,nvar,ng); + VisMF::Read(lsm_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "LsmVars")); + MultiFab::Copy(*(lsm_data[lev][mvar]),lsm_vars,0,0,nvar,ng); + } + } + // Note that we read the ghost cells of the mapfactors BoxList bl2d = grids[lev].boxList(); for (auto& b : bl2d) { diff --git a/Source/IO/NCCheckpoint.cpp b/Source/IO/NCCheckpoint.cpp index 24df09c6c..e489a44a8 100644 --- a/Source/IO/NCCheckpoint.cpp +++ b/Source/IO/NCCheckpoint.cpp @@ -125,7 +125,7 @@ ERF::WriteNCCheckpointFile () const void ERF::ReadNCCheckpointFile () { - amrex::Print() << "Restart from checkpoint " << restart_chkfile << "\n"; + amrex::Print() << "Restart from NetCDF checkpoint " << restart_chkfile << "\n"; // Header std::string HeaderFileName(restart_chkfile + "/Header.nc"); diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 1f75e5742..77040ae25 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -895,10 +895,15 @@ ERF::WritePlotFile (int which, Vector plot_var_names) else if (which == 2) plotfilename = Concatenate(plot_file_2, istep[0], 5); + // LSM writes it's own data + if (which==1 && plot_lsm) { + lsm.Plot_Lsm_Data(t_new[0], istep, refRatio()); + } + if (finest_level == 0) { if (plotfile_type == "amrex") { - amrex::Print() << "Writing plotfile " << plotfilename << "\n"; + amrex::Print() << "Writing native plotfile " << plotfilename << "\n"; if (solverChoice.use_terrain) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, GetVecOfConstPtrs(mf), diff --git a/Source/LandSurfaceModel/LandSurface.H b/Source/LandSurfaceModel/LandSurface.H new file mode 100644 index 000000000..37f56d919 --- /dev/null +++ b/Source/LandSurfaceModel/LandSurface.H @@ -0,0 +1,128 @@ +#ifndef LANDSURFACE_H +#define LANDSURFACE_H + +#include + +#include +#include +#include + +class LandSurface { + +public: + + LandSurface () { } + + ~LandSurface () = default; + + void + ReSize (const int& nlev) { m_lsm_model.resize(nlev); } + + template + void + SetModel () + { + for (int lev(0); lev(); + } + } + + void + Define (const int& lev, + SolverChoice& sc) + { + m_lsm_model[lev]->Define(sc); + } + + void + Init (const int& lev, + const amrex::MultiFab& cons_in, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) + { + m_lsm_model[lev]->Init(cons_in, geom, dt_advance); + } + + void + Advance (const int& lev, const amrex::Real& dt_advance) + { + m_lsm_model[lev]->Advance(dt_advance); + } + + void + Update_Micro_Vars_Lev (const int& lev, amrex::MultiFab& cons_in) + { + m_lsm_model[lev]->Update_Micro_Vars(cons_in); + } + + void + Update_State_Vars_Lev (const int& lev, amrex::MultiFab& cons_in) + { + m_lsm_model[lev]->Update_State_Vars(cons_in); + } + + amrex::MultiFab* + Get_Data_Ptr (const int& lev, const int& varIdx) { return m_lsm_model[lev]->Lsm_Data_Ptr(varIdx); } + + amrex::MultiFab* + Get_Flux_Ptr (const int& lev, const int& varIdx) { return m_lsm_model[lev]->Lsm_Flux_Ptr(varIdx); } + + amrex::Geometry + Get_Lsm_Geom (const int& lev ) { return m_lsm_model[lev]->Lsm_Geom(); } + + int + Get_Data_Size () { return m_lsm_model[0]->Lsm_Data_Size(); } + + std::string + Get_VarName (const int& varIdx) { return m_lsm_model[0]->Lsm_VarName(varIdx); } + + void + Plot_Lsm_Data (amrex::Real time, + const amrex::Vector &level_steps, + const amrex::Vector &ref_ratio) + { + int nlev = m_lsm_model.size(); + int nvar = this->Get_Data_Size(); + amrex::Vector varnames; + + // Only write if we have valid pointers + if (this->Get_Data_Ptr(0,0)) { + varnames.resize(nvar); + m_lsm_geom_lev.resize(nlev); + m_lsm_data_lev.resize(nlev); + std::string plotfilename = amrex::Concatenate(plot_file_lsm, level_steps[0], 5); + for (int lev(0); levGet_Lsm_Geom(lev); + + amrex::MultiFab* mf_lsm = this->Get_Data_Ptr(lev,0); + amrex::IntVect ng(0,0,0); + amrex::BoxArray ba = mf_lsm->boxArray(); + amrex::DistributionMapping dm = mf_lsm->DistributionMap(); + m_lsm_data_lev[lev].define(ba, dm, nvar, ng); + for (int n(0); nGet_Data_Ptr(lev,n); + amrex::MultiFab::Copy(m_lsm_data_lev[lev],*(mf_lsm),0,n,1,0); + if (lev==0) varnames[n] = this->Get_VarName(n); + } + } + WriteMultiLevelPlotfile (plotfilename, nlev, GetVecOfConstPtrs(m_lsm_data_lev), + varnames, m_lsm_geom_lev, time, level_steps, ref_ratio); + m_lsm_geom_lev.clear(); + m_lsm_data_lev.clear(); + } + } + +private: + // lsm model at each level + amrex::Vector> m_lsm_model; + + // plotfile prefix + std::string plot_file_lsm {"plt_lsm_"}; + + // Vector of geometry + amrex::Vector m_lsm_geom_lev; + + // Vector of data pointers + amrex::Vector m_lsm_data_lev; +}; +#endif diff --git a/Source/LandSurfaceModel/MM5/MM5.H b/Source/LandSurfaceModel/MM5/MM5.H new file mode 100644 index 000000000..325a2f091 --- /dev/null +++ b/Source/LandSurfaceModel/MM5/MM5.H @@ -0,0 +1,163 @@ +#ifndef MM5_H +#define MM5_H + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace LsmVar_MM5 { + enum { + // independent variables + theta = 0, + NumVars + }; +} + +/* See Chen & Dudhia (2001) https://doi.org/10.1175/1520-0493(2001)129<0569:CAALSH>2.0.CO;2 */ +class MM5 : public NullSurf { + + using FabPtr = std::shared_ptr; + +public: + // Constructor + MM5 () {} + + // Destructor + virtual ~MM5 () = default; + + // Set thermo and grid properties + void + Define (SolverChoice& /*sc*/) override + { + // NOTE: We should parse things from sc here, + // but they are hard coded because this + // is a demonstration for now. + } + + // Initialize data structures + void + Init (const amrex::MultiFab& cons_in, + const amrex::Geometry& geom, + const amrex::Real& dt) override; + + // Wrapper to do all the updating + void + Advance (const amrex::Real& dt) override + { + m_dt = dt; + this->ComputeFluxes(); + this->AdvanceMM5(); + this->ComputeTsurf(); + } + + // Compute surface temperature + void + ComputeTsurf (); + + // Compute diffusive fluxes + void + ComputeFluxes (); + + // Advance the lsm state vars + void + AdvanceMM5 (); + + // Get state vars from lsm class + amrex::MultiFab* + Lsm_Data_Ptr (const int& varIdx) override + { + int lsmIdx = LsmVarMap[varIdx]; + AMREX_ALWAYS_ASSERT(lsmIdx < MM5::m_lsm_size && lsmIdx>=0); + return lsm_fab_vars[lsmIdx].get(); + } + + // Get flux vars from lsm class + amrex::MultiFab* + Lsm_Flux_Ptr (const int& varIdx) override + { + int lsmIdx = LsmVarMap[varIdx]; + AMREX_ALWAYS_ASSERT(lsmIdx < MM5::m_lsm_size && lsmIdx>=0); + return lsm_fab_flux[lsmIdx].get(); + } + + // Get lsm geometry + amrex::Geometry + Lsm_Geom ( ) override { return m_lsm_geom; } + + // Get number of vars lsm class contains + int + Lsm_Data_Size () override { return MM5::m_lsm_size; } + + // Get variable names + std::string + Lsm_VarName (const int& varIdx) override + { + int lsmIdx = LsmVarMap[varIdx]; + AMREX_ALWAYS_ASSERT(lsmIdx < MM5::m_lsm_size && lsmIdx>=0); + return LsmVarName[lsmIdx]; + } + +private: + // number of lsm variables (theta) + int m_lsm_size = 1; + + // LsmVar map (state indices -> LsmVar enum) + amrex::Vector LsmVarMap; + + // Lsm varnames + amrex::Vector LsmVarName; + + // geometry for atmosphere + amrex::Geometry m_geom; + + // geometry for lsm + amrex::Geometry m_lsm_geom; + + // timestep + amrex::Real m_dt; + + // domain klo-1 or lsm khi + int khi_lsm; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // independent variables + amrex::Array lsm_fab_vars; + + // flux array for conjugate transfer + amrex::Array lsm_fab_flux; + + // Vars that should be parsed + // ========================== + // Number of grid points in z + int m_nz_lsm = 30; + + // Size of grid spacing in z + amrex::Real m_dz_lsm = 0.1; // 3 [m] below surface + + // Specific heat + amrex::Real m_cp_soil = 1.26e6; // [J/m^3 K] + + // Conductivity + amrex::Real m_k_soil = 0.2; // dry soil [W/m K] + + // Theta Dirichlet value at lowest point below surface + amrex::Real m_theta_dir = 300.0; // [K] + + // Thermal diffusivity + amrex::Real m_d_soil = m_k_soil / m_cp_soil; // [m^2/s] +}; +#endif diff --git a/Source/LandSurfaceModel/MM5/MM5.cpp b/Source/LandSurfaceModel/MM5/MM5.cpp new file mode 100644 index 000000000..b548b1d54 --- /dev/null +++ b/Source/LandSurfaceModel/MM5/MM5.cpp @@ -0,0 +1,124 @@ +#include + +using namespace amrex; + +/* Initialize lsm data structures */ +void +MM5::Init (const MultiFab& cons_in, + const Geometry& geom, + const Real& dt) +{ + m_dt = dt; + m_geom = geom; + + Box domain = geom.Domain(); + khi_lsm = domain.smallEnd(2) - 1; + + LsmVarMap.resize(m_lsm_size); + LsmVarMap = {LsmVar_MM5::theta}; + + LsmVarName.resize(m_lsm_size); + LsmVarName = {"theta"}; + + // NOTE: All boxes in ba extend from zlo to zhi, so this transform is valid. + // If that were to change, the dm and new ba are no longer valid and + // direct copying between lsm data/flux vars cannot be done in a parfor. + + // Set box array for lsm data + IntVect ng(0,0,1); + BoxArray ba = cons_in.boxArray(); + DistributionMapping dm = cons_in.DistributionMap(); + BoxList bl_lsm = ba.boxList(); + for (auto& b : bl_lsm) { + b.setBig(2,khi_lsm); // First point below the surface + b.setSmall(2,khi_lsm - m_nz_lsm + 1); // Last point below the surface + } + BoxArray ba_lsm(std::move(bl_lsm)); + + // Set up lsm geometry + const RealBox& dom_rb = m_geom.ProbDomain(); + const Real* dom_dx = m_geom.CellSize(); + RealBox lsm_rb = dom_rb; + Real lsm_dx[AMREX_SPACEDIM] = {AMREX_D_DECL(dom_dx[0],dom_dx[1],m_dz_lsm)}; + Real lsm_z_hi = dom_rb.lo(2); + Real lsm_z_lo = lsm_z_hi - Real(m_nz_lsm)*lsm_dx[2]; + lsm_rb.setHi(2,lsm_z_hi); lsm_rb.setLo(2,lsm_z_lo); + m_lsm_geom.define( ba_lsm.minimalBox(), lsm_rb, m_geom.Coord(), m_geom.isPeriodic() ); + + // Create the data and fluxes + for (auto ivar = 0; ivar < LsmVar_MM5::NumVars; ++ivar) { + // State vars are CC + Real theta_0 = m_theta_dir; + lsm_fab_vars[ivar] = std::make_shared(ba_lsm, dm, 1, ng); + lsm_fab_vars[ivar]->setVal(theta_0); + + // Fluxes are nodal in z + lsm_fab_flux[ivar] = std::make_shared(convert(ba_lsm, IntVect(0,0,1)), dm, 1, IntVect(0,0,0)); + lsm_fab_flux[ivar]->setVal(0.); + } +} + +/* Extrapolate surface temperature and store in ghost cell */ +void +MM5::ComputeTsurf () +{ + // Expose for GPU copy + int khi = khi_lsm; + + for ( MFIter mfi(*(lsm_fab_vars[LsmVar_MM5::theta])); mfi.isValid(); ++mfi) { + auto box2d = mfi.tilebox(); box2d.makeSlab(2,khi); + + auto theta_array = lsm_fab_vars[LsmVar_MM5::theta]->array(mfi); + + ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) + { + theta_array(i,j,khi+1) = 1.5*theta_array(i,j,khi) - 0.5*theta_array(i,j,khi-1); + }); + } +} + +/* Compute the diffusive fluxes */ +void +MM5::ComputeFluxes () +{ + // Expose for GPU copy + int khi = khi_lsm; + Real Dsoil = m_d_soil; + Real dzInv = m_lsm_geom.InvCellSize(2); + + for ( MFIter mfi(*(lsm_fab_flux[LsmVar_MM5::theta])); mfi.isValid(); ++mfi) { + auto box3d = mfi.tilebox(); + + // Do not overwrite the flux at the top (comes from MOST BC) + if (box3d.bigEnd(2) == khi+1) box3d.setBig(2,khi); + + auto theta_array = lsm_fab_vars[LsmVar_MM5::theta]->array(mfi); + auto theta_flux = lsm_fab_flux[LsmVar_MM5::theta]->array(mfi); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + theta_flux(i,j,k) = Dsoil * ( theta_array(i,j,k) - theta_array(i,j,k-1) ) * dzInv; + }); + } +} + +/* Advance the solution with a simple explicit update (should use tridiagonal solve) */ +void +MM5::AdvanceMM5 () +{ + // Expose for GPU copy + Real dt = m_dt; + Real dzInv = m_lsm_geom.InvCellSize(2); + + for ( MFIter mfi(*(lsm_fab_vars[LsmVar_MM5::theta])); mfi.isValid(); ++mfi) { + auto box3d = mfi.tilebox(); + + auto theta_array = lsm_fab_vars[LsmVar_MM5::theta]->array(mfi); + auto theta_flux = lsm_fab_flux[LsmVar_MM5::theta]->array(mfi); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + theta_array(i,j,k) += dt * ( theta_flux(i,j,k+1) - theta_flux(i,j,k) ) * dzInv; + }); + } +} diff --git a/Source/LandSurfaceModel/MM5/Make.package b/Source/LandSurfaceModel/MM5/Make.package new file mode 100644 index 000000000..c844c2c58 --- /dev/null +++ b/Source/LandSurfaceModel/MM5/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += MM5.cpp +CEXE_headers += MM5.H diff --git a/Source/LandSurfaceModel/Make.package b/Source/LandSurfaceModel/Make.package new file mode 100644 index 000000000..4bdd1f1c0 --- /dev/null +++ b/Source/LandSurfaceModel/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += LandSurface.H + diff --git a/Source/LandSurfaceModel/Null/Make.package b/Source/LandSurfaceModel/Null/Make.package new file mode 100644 index 000000000..c8876d592 --- /dev/null +++ b/Source/LandSurfaceModel/Null/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += NullSurf.H + diff --git a/Source/LandSurfaceModel/Null/NullSurf.H b/Source/LandSurfaceModel/Null/NullSurf.H new file mode 100644 index 000000000..ccdc485dc --- /dev/null +++ b/Source/LandSurfaceModel/Null/NullSurf.H @@ -0,0 +1,72 @@ +#ifndef NULLSURF_H +#define NULLSURF_H + +#include +#include +#include + +class NullSurf { + + public: + NullSurf () {} + + virtual ~NullSurf () = default; + + virtual + void + Define (SolverChoice& /*sc*/) { } + + virtual + void Init (const amrex::MultiFab& /*cons_in*/, + const amrex::Geometry& /*geom*/, + const amrex::Real& /*dt_advance*/) { } + + virtual + void + Advance (const amrex::Real& /*dt_advance*/) { } + + virtual + void + Update_Micro_Vars (amrex::MultiFab& /*cons_in*/) { } + + virtual + void + Update_State_Vars (amrex::MultiFab& /*cons_in*/) { } + + virtual + void + Copy_State_to_Micro (const amrex::MultiFab& /*cons_in*/) { } + + virtual + void + Copy_Micro_to_State (amrex::MultiFab& /*cons_in*/) { } + + virtual + amrex::MultiFab* + Lsm_Data_Ptr (const int& /*varIdx*/ ) { return nullptr; } + + virtual + amrex::MultiFab* + Lsm_Flux_Ptr (const int& /*varIdx*/ ) { return nullptr; } + + virtual + amrex::Geometry + Lsm_Geom ( ) { return m_lsm_geom; } + + virtual + int + Lsm_Data_Size () { return NullSurf::m_lsm_size; } + + virtual + std::string + Lsm_VarName (const int& /*varIdx*/) + { + return varname; + } + + private: + int m_lsm_size = 1; + amrex::Geometry m_lsm_geom; + std::string varname = " "; +}; +#endif diff --git a/Source/LandSurfaceModel/SLM/Make.package b/Source/LandSurfaceModel/SLM/Make.package new file mode 100644 index 000000000..cc472527e --- /dev/null +++ b/Source/LandSurfaceModel/SLM/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += SLM.cpp +CEXE_headers += SLM.H diff --git a/Source/LandSurfaceModel/SLM/SLM.H b/Source/LandSurfaceModel/SLM/SLM.H new file mode 100644 index 000000000..22dcf592e --- /dev/null +++ b/Source/LandSurfaceModel/SLM/SLM.H @@ -0,0 +1,163 @@ +#ifndef SLM_H +#define SLM_H + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace LsmVar_SLM { + enum { + // independent variables + theta = 0, + NumVars + }; +} + +/* See Chen & Dudhia (2001) https://doi.org/10.1175/1520-0493(2001)129<0569:CAALSH>2.0.CO;2 */ +class SLM : public NullSurf { + + using FabPtr = std::shared_ptr; + +public: + // Constructor + SLM () {} + + // Destructor + virtual ~SLM () = default; + + // Set thermo and grid properties + void + Define (SolverChoice& /*sc*/) override + { + // NOTE: We should parse things from sc here, + // but they are hard coded because this + // is a demonstration for now. + } + + // Initialize data structures + void + Init (const amrex::MultiFab& cons_in, + const amrex::Geometry& geom, + const amrex::Real& dt) override; + + // Wrapper to do all the updating + void + Advance (const amrex::Real& dt) override + { + m_dt = dt; + this->ComputeFluxes(); + this->AdvanceSLM(); + this->ComputeTsurf(); + } + + // Compute surface temperature + void + ComputeTsurf (); + + // Compute diffusive fluxes + void + ComputeFluxes (); + + // Advance the lsm state vars + void + AdvanceSLM (); + + // Get state vars from lsm class + amrex::MultiFab* + Lsm_Data_Ptr (const int& varIdx) override + { + int lsmIdx = LsmVarMap[varIdx]; + AMREX_ALWAYS_ASSERT(lsmIdx < SLM::m_lsm_size && lsmIdx>=0); + return lsm_fab_vars[lsmIdx].get(); + } + + // Get flux vars from lsm class + amrex::MultiFab* + Lsm_Flux_Ptr (const int& varIdx) override + { + int lsmIdx = LsmVarMap[varIdx]; + AMREX_ALWAYS_ASSERT(lsmIdx < SLM::m_lsm_size && lsmIdx>=0); + return lsm_fab_flux[lsmIdx].get(); + } + + // Get lsm geometry + amrex::Geometry + Lsm_Geom ( ) override { return m_lsm_geom; } + + // Get number of vars lsm class contains + int + Lsm_Data_Size () override { return SLM::m_lsm_size; } + + // Get variable names + std::string + Lsm_VarName (const int& varIdx) override + { + int lsmIdx = LsmVarMap[varIdx]; + AMREX_ALWAYS_ASSERT(lsmIdx < SLM::m_lsm_size && lsmIdx>=0); + return LsmVarName[lsmIdx]; + } + +private: + // number of lsm variables (theta) + int m_lsm_size = 1; + + // LsmVar map (state indices -> LsmVar enum) + amrex::Vector LsmVarMap; + + // Lsm varnames + amrex::Vector LsmVarName; + + // geometry for atmosphere + amrex::Geometry m_geom; + + // geometry for lsm + amrex::Geometry m_lsm_geom; + + // timestep + amrex::Real m_dt; + + // domain klo-1 or lsm khi + int khi_lsm; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // independent variables + amrex::Array lsm_fab_vars; + + // flux array for conjugate transfer + amrex::Array lsm_fab_flux; + + // Vars that should be parsed + // ========================== + // Number of grid points in z + int m_nz_lsm = 30; + + // Size of grid spacing in z + amrex::Real m_dz_lsm = 0.1; // 3 [m] below surface + + // Specific heat + amrex::Real m_cp_soil = 1.26e6; // [J/m^3 K] + + // Conductivity + amrex::Real m_k_soil = 0.2; // dry soil [W/m K] + + // Theta Dirichlet value at lowest point below surface + amrex::Real m_theta_dir = 400.0; // [K] + + // Thermal diffusivity + amrex::Real m_d_soil = m_k_soil / m_cp_soil; // [m^2/s] +}; +#endif diff --git a/Source/LandSurfaceModel/SLM/SLM.cpp b/Source/LandSurfaceModel/SLM/SLM.cpp new file mode 100644 index 000000000..d9e7905dd --- /dev/null +++ b/Source/LandSurfaceModel/SLM/SLM.cpp @@ -0,0 +1,124 @@ +#include + +using namespace amrex; + +/* Initialize lsm data structures */ +void +SLM::Init (const MultiFab& cons_in, + const Geometry& geom, + const Real& dt) +{ + m_dt = dt; + m_geom = geom; + + Box domain = geom.Domain(); + khi_lsm = domain.smallEnd(2) - 1; + + LsmVarMap.resize(m_lsm_size); + LsmVarMap = {LsmVar_SLM::theta}; + + LsmVarName.resize(m_lsm_size); + LsmVarName = {"theta"}; + + // NOTE: All boxes in ba extend from zlo to zhi, so this transform is valid. + // If that were to change, the dm and new ba are no longer valid and + // direct copying between lsm data/flux vars cannot be done in a parfor. + + // Set box array for lsm data + IntVect ng(0,0,1); + BoxArray ba = cons_in.boxArray(); + DistributionMapping dm = cons_in.DistributionMap(); + BoxList bl_lsm = ba.boxList(); + for (auto& b : bl_lsm) { + b.setBig(2,khi_lsm); // First point below the surface + b.setSmall(2,khi_lsm - m_nz_lsm + 1); // Last point below the surface + } + BoxArray ba_lsm(std::move(bl_lsm)); + + // Set up lsm geometry + const RealBox& dom_rb = m_geom.ProbDomain(); + const Real* dom_dx = m_geom.CellSize(); + RealBox lsm_rb = dom_rb; + Real lsm_dx[AMREX_SPACEDIM] = {AMREX_D_DECL(dom_dx[0],dom_dx[1],m_dz_lsm)}; + Real lsm_z_hi = dom_rb.lo(2); + Real lsm_z_lo = lsm_z_hi - Real(m_nz_lsm)*lsm_dx[2]; + lsm_rb.setHi(2,lsm_z_hi); lsm_rb.setLo(2,lsm_z_lo); + m_lsm_geom.define( ba_lsm.minimalBox(), lsm_rb, m_geom.Coord(), m_geom.isPeriodic() ); + + // Create the data and fluxes + for (auto ivar = 0; ivar < LsmVar_SLM::NumVars; ++ivar) { + // State vars are CC + Real theta_0 = m_theta_dir; + lsm_fab_vars[ivar] = std::make_shared(ba_lsm, dm, 1, ng); + lsm_fab_vars[ivar]->setVal(theta_0); + + // Fluxes are nodal in z + lsm_fab_flux[ivar] = std::make_shared(convert(ba_lsm, IntVect(0,0,1)), dm, 1, IntVect(0,0,0)); + lsm_fab_flux[ivar]->setVal(0.); + } +} + +/* Extrapolate surface temperature and store in ghost cell */ +void +SLM::ComputeTsurf () +{ + // Expose for GPU copy + int khi = khi_lsm; + + for ( MFIter mfi(*(lsm_fab_vars[LsmVar_SLM::theta])); mfi.isValid(); ++mfi) { + auto box2d = mfi.tilebox(); box2d.makeSlab(2,khi); + + auto theta_array = lsm_fab_vars[LsmVar_SLM::theta]->array(mfi); + + ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) + { + theta_array(i,j,khi+1) = 1.5*theta_array(i,j,khi) - 0.5*theta_array(i,j,khi-1); + }); + } +} + +/* Compute the diffusive fluxes */ +void +SLM::ComputeFluxes () +{ + // Expose for GPU copy + int khi = khi_lsm; + Real Dsoil = m_d_soil; + Real dzInv = m_lsm_geom.InvCellSize(2); + + for ( MFIter mfi(*(lsm_fab_flux[LsmVar_SLM::theta])); mfi.isValid(); ++mfi) { + auto box3d = mfi.tilebox(); + + // Do not overwrite the flux at the top (comes from MOST BC) + if (box3d.bigEnd(2) == khi+1) box3d.setBig(2,khi); + + auto theta_array = lsm_fab_vars[LsmVar_SLM::theta]->array(mfi); + auto theta_flux = lsm_fab_flux[LsmVar_SLM::theta]->array(mfi); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + theta_flux(i,j,k) = Dsoil * ( theta_array(i,j,k) - theta_array(i,j,k-1) ) * dzInv; + }); + } +} + +/* Advance the solution with a simple explicit update (should use tridiagonal solve) */ +void +SLM::AdvanceSLM () +{ + // Expose for GPU copy + Real dt = m_dt; + Real dzInv = m_lsm_geom.InvCellSize(2); + + for ( MFIter mfi(*(lsm_fab_vars[LsmVar_SLM::theta])); mfi.isValid(); ++mfi) { + auto box3d = mfi.tilebox(); + + auto theta_array = lsm_fab_vars[LsmVar_SLM::theta]->array(mfi); + auto theta_flux = lsm_fab_flux[LsmVar_SLM::theta]->array(mfi); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + theta_array(i,j,k) += dt * ( theta_flux(i,j,k+1) - theta_flux(i,j,k) ) * dzInv; + }); + } +} diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 17391a0f5..92e48c347 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -128,6 +128,9 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ // Update the microphysics (moisture) advance_microphysics(lev, S_new, dt_lev); + // Update the land surface model + advance_lsm(lev, S_new, dt_lev); + #if defined(ERF_USE_RRTMGP) // Update the radiation advance_radiation(lev, S_new, dt_lev); diff --git a/Source/TimeIntegration/ERF_advance_lsm.cpp b/Source/TimeIntegration/ERF_advance_lsm.cpp new file mode 100644 index 000000000..fd33cefc8 --- /dev/null +++ b/Source/TimeIntegration/ERF_advance_lsm.cpp @@ -0,0 +1,12 @@ +#include + +using namespace amrex; + +void ERF::advance_lsm (int lev, + MultiFab& /*cons*/, + const Real& dt_advance) +{ + if (solverChoice.lsm_type != LandSurfaceType::None) { + lsm.Advance(lev, dt_advance); + } +} diff --git a/Source/TimeIntegration/Make.package b/Source/TimeIntegration/Make.package index 162bff16b..fb58562b6 100644 --- a/Source/TimeIntegration/Make.package +++ b/Source/TimeIntegration/Make.package @@ -3,6 +3,7 @@ CEXE_sources += ERF_Advance.cpp CEXE_sources += ERF_TimeStep.cpp CEXE_sources += ERF_advance_dycore.cpp CEXE_sources += ERF_advance_microphysics.cpp +CEXE_sources += ERF_advance_lsm.cpp CEXE_sources += ERF_make_buoyancy.cpp CEXE_sources += ERF_make_fast_coeffs.cpp CEXE_sources += ERF_slow_rhs_pre.cpp