From e248a47774d55e5cdbe72d460e74e9c3ceb759f6 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Thu, 16 Nov 2023 13:07:05 -0800 Subject: [PATCH] Does not compile with templated Microphysics.H but does compile with the case select version (commented out). --- CMake/BuildERFExe.cmake | 16 +- Exec/Make.ERF | 10 + .../WPS_Test/inputs_real_ChisholmView | 2 + Source/ERF.H | 3 +- Source/ERF.cpp | 11 + Source/Microphysics/Make.package | 7 - Source/Microphysics/Microphysics.H | 305 ++++++++--------- Source/Microphysics/Null/Make.package | 2 + Source/Microphysics/Null/Make.package~ | 2 + Source/Microphysics/Null/NullMoist.H | 35 ++ Source/Microphysics/Null/NullMoist.H~ | 25 ++ Source/Microphysics/SAM/Cloud.cpp | 134 ++++++++ .../{Cloud.cpp => SAM/Cloud.cpp~} | 0 Source/Microphysics/SAM/Diagnose.cpp | 59 ++++ .../{Diagnose.cpp => SAM/Diagnose.cpp~} | 0 Source/Microphysics/SAM/IceFall.cpp | 171 ++++++++++ .../{IceFall.cpp => SAM/IceFall.cpp~} | 0 Source/Microphysics/SAM/Init.cpp | 245 ++++++++++++++ .../Microphysics/{Init.cpp => SAM/Init.cpp~} | 0 Source/Microphysics/SAM/Make.package | 9 + Source/Microphysics/SAM/Precip.cpp | 148 +++++++++ .../{Precip.cpp => SAM/Precip.cpp~} | 0 Source/Microphysics/SAM/PrecipFall.cpp | 307 ++++++++++++++++++ .../{PrecipFall.cpp => SAM/PrecipFall.cpp~} | 0 Source/Microphysics/SAM/SAM.H | 177 ++++++++++ Source/Microphysics/SAM/SAM.H~ | 163 ++++++++++ Source/Microphysics/SAM/Update.cpp | 60 ++++ .../{Update.cpp => SAM/Update.cpp~} | 0 .../ERF_advance_microphysics.cpp | 8 +- 29 files changed, 1729 insertions(+), 170 deletions(-) create mode 100644 Source/Microphysics/Null/Make.package create mode 100644 Source/Microphysics/Null/Make.package~ create mode 100644 Source/Microphysics/Null/NullMoist.H create mode 100644 Source/Microphysics/Null/NullMoist.H~ create mode 100644 Source/Microphysics/SAM/Cloud.cpp rename Source/Microphysics/{Cloud.cpp => SAM/Cloud.cpp~} (100%) create mode 100644 Source/Microphysics/SAM/Diagnose.cpp rename Source/Microphysics/{Diagnose.cpp => SAM/Diagnose.cpp~} (100%) create mode 100644 Source/Microphysics/SAM/IceFall.cpp rename Source/Microphysics/{IceFall.cpp => SAM/IceFall.cpp~} (100%) create mode 100644 Source/Microphysics/SAM/Init.cpp rename Source/Microphysics/{Init.cpp => SAM/Init.cpp~} (100%) create mode 100644 Source/Microphysics/SAM/Make.package create mode 100644 Source/Microphysics/SAM/Precip.cpp rename Source/Microphysics/{Precip.cpp => SAM/Precip.cpp~} (100%) create mode 100644 Source/Microphysics/SAM/PrecipFall.cpp rename Source/Microphysics/{PrecipFall.cpp => SAM/PrecipFall.cpp~} (100%) create mode 100644 Source/Microphysics/SAM/SAM.H create mode 100644 Source/Microphysics/SAM/SAM.H~ create mode 100644 Source/Microphysics/SAM/Update.cpp rename Source/Microphysics/{Update.cpp => SAM/Update.cpp~} (100%) diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index b9219ad4b..01c3e196e 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -58,13 +58,13 @@ function(build_erf_lib erf_lib_name) if(ERF_ENABLE_MOISTURE) target_sources(${erf_lib_name} PRIVATE - ${SRC_DIR}/Microphysics/Init.cpp - ${SRC_DIR}/Microphysics/Cloud.cpp - ${SRC_DIR}/Microphysics/IceFall.cpp - ${SRC_DIR}/Microphysics/Precip.cpp - ${SRC_DIR}/Microphysics/PrecipFall.cpp - ${SRC_DIR}/Microphysics/Diagnose.cpp - ${SRC_DIR}/Microphysics/Update.cpp) + ${SRC_DIR}/Microphysics/SAM/Init.cpp + ${SRC_DIR}/Microphysics/SAM/Cloud.cpp + ${SRC_DIR}/Microphysics/SAM/IceFall.cpp + ${SRC_DIR}/Microphysics/SAM/Precip.cpp + ${SRC_DIR}/Microphysics/SAM/PrecipFall.cpp + ${SRC_DIR}/Microphysics/SAM/Diagnose.cpp + ${SRC_DIR}/Microphysics/SAM/Update.cpp) target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MOISTURE) endif() @@ -184,6 +184,8 @@ function(build_erf_lib erf_lib_name) if(ERF_ENABLE_MOISTURE) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM) endif() if(ERF_ENABLE_RRTMGP) diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 34d4b8371..0505ccf8d 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -107,6 +107,16 @@ ifeq ($(USE_MOISTURE), TRUE) include $(ERF_MOISTURE_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_DIR) INCLUDE_LOCATIONS += $(ERF_MOISTURE_DIR) + + ERF_MOISTURE_NULL_DIR = $(ERF_SOURCE_DIR)/Microphysics/Null + include $(ERF_MOISTURE_NULL_DIR)/Make.package + VPATH_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) + INCLUDE_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) + + ERF_MOISTURE_SAM_DIR = $(ERF_SOURCE_DIR)/Microphysics/SAM + include $(ERF_MOISTURE_SAM_DIR)/Make.package + VPATH_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) + INCLUDE_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) endif ifeq ($(USE_WARM_NO_PRECIP), TRUE) diff --git a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView index 8c3c14032..d077f9e88 100644 --- a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView +++ b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView @@ -53,6 +53,8 @@ erf.use_gravity = true erf.use_terrain = true erf.use_terrain = false +erf.moisture_model = "SAM" + erf.les_type = "None" erf.molec_diff_type = "Constant" erf.dynamicViscosity = 5.0 diff --git a/Source/ERF.H b/Source/ERF.H index 455eb8a0a..bf2ebfbe0 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -543,7 +543,8 @@ private: amrex::Vector rW_new; #if defined(ERF_USE_MOISTURE) - Microphysics micro; + std::string moisture_model = "NullMoist"; + Microphysics micro; amrex::Vector qmoist; // This has 6 components: qv, qc, qi, qr, qs, qg #endif diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 95bc5b29d..a4d41f030 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -967,6 +967,17 @@ ERF::ReadParameters () pp.query("use_tracer_particles", use_tracer_particles); #endif +#ifdef ERF_USE_MOISTURE + // What type of moisture model to use + pp.query("moisture_model", moisture_model); + if (moisture_model == "SAM") { + SAM T; + micro.SetModel(T); + } + + //micro.SetModel(moisture_model); +#endif + // If this is set, it must be even if (fixed_mri_dt_ratio > 0 && (fixed_mri_dt_ratio%2 != 0) ) { diff --git a/Source/Microphysics/Make.package b/Source/Microphysics/Make.package index 02f599048..36572d9e5 100644 --- a/Source/Microphysics/Make.package +++ b/Source/Microphysics/Make.package @@ -1,9 +1,2 @@ -CEXE_sources += Cloud.cpp -CEXE_sources += Init.cpp -CEXE_sources += Diagnose.cpp -CEXE_sources += Update.cpp -CEXE_sources += IceFall.cpp -CEXE_sources += Precip.cpp -CEXE_sources += PrecipFall.cpp CEXE_headers += Microphysics.H diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index 09e538320..419c19c31 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -1,163 +1,172 @@ -/* - * Implementation 1-moment microphysics model - * NOTE: this model is based on the SAM code, and the Klemp's paper - * 1): Joseph, Klemp, the simulation of three-dimensional convective storm dynamics, - * Journal of the atmospheric sciences, vol35, p1070 - * 2): Marat Khairoutdinov and David Randall, cloud resolving modeling of the ARM summer 1997 IOP: - * model formulation, results, unvertainties, and sensitivities, Journal of the atmospheric sciences, vol60, p607 - */ #ifndef MICROPHYSICS_H #define MICROPHYSICS_H -#include -#include -#include +#include +#include + +template +class Microphysics : public MoistModel { + +public: + + Microphysics () { } + + ~Microphysics () = default; + + template + void + SetModel (NewMoistModel& T ) + { + m_moist_model.reset(new NewMoistModel()); + m_moist_model = std::make_unique(); + } + + void + define (SolverChoice& sc) + { + m_moist_model->define(sc); + } + + void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) + { + m_moist_model->Init(cons_in, qmoist, grids, geom, dt_advance); + } + + void + Advance ( ) + { + m_moist_model->Advance(); + } + + void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist) + { + m_moist_model->Update(cons_in, qmoist); + } + + +private: + std::unique_ptr m_moist_model; +}; +#endif + -#include -#include -#include -#include +/* +#ifndef MICROPHYSICS_H +#define MICROPHYSICS_H -#include "ERF_Constants.H" -#include "Microphysics_Utils.H" -#include "IndexDefines.H" -#include "DataStruct.H" +#include +#include -namespace MicVar { +// Available models +namespace MicMod { enum { - // independent variables - qt = 0, - qp, - theta, // liquid/ice water potential temperature - tabs, // temperature - rho, // density - pres, // pressure - // derived variables - qr, // rain water - qv, // water vapor - qn, // cloud condensate (liquid+ice), initial to zero - qci, // cloud ice - qcl, // cloud water - qpl, // precip rain - qpi, // precip ice - // temporary variable - omega, + Null = 0, + SAM, + Kessler, NumVars }; } -// -// use MultiFab for 3D data, but table for 1D data -// class Microphysics { - using FabPtr = std::shared_ptr; - - public: - // constructor - Microphysics () {} - - // Set up for first time - void define (SolverChoice& sc) - { - docloud = sc.do_cloud; - doprecip = sc.do_precip; - m_fac_cond = lcond / sc.c_p; - m_fac_fus = lfus / sc.c_p; - m_fac_sub = lsub / sc.c_p; - m_gOcp = CONST_GRAV / sc.c_p; - m_axis = sc.ave_plane; - } - - // destructor - ~Microphysics () = default; - - // cloud physics - void Cloud (); - - // ice physics - void IceFall (); - - // precip - void Precip (); - - // precip fall - void PrecipFall (int hydro_type); - - // micro interface for precip fall - void MicroPrecipFall (); - - // init - void Init (const amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist, - const amrex::BoxArray& grids, - const amrex::Geometry& geom, - const amrex::Real& dt_advance); - - // update ERF variables - void Update (amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist); - - // diagnose - void Diagnose (); - - // process microphysics - void Proc (); - - private: - // geometry - amrex::Geometry m_geom; - // valid boxes on which to evolve the solution - amrex::BoxArray m_gtoe; - - // timestep - amrex::Real dt; - - // number of vertical levels - int nlev, zlo, zhi; - - // plane average axis - int m_axis; - - // model options - bool docloud, doprecip; - - // constants - amrex::Real m_fac_cond; - amrex::Real m_fac_fus; - amrex::Real m_fac_sub; - amrex::Real m_gOcp; - - // microphysics parameters/coefficients - amrex::TableData accrrc; - amrex::TableData accrsi; - amrex::TableData accrsc; - amrex::TableData coefice; - amrex::TableData evaps1; - amrex::TableData evaps2; - amrex::TableData accrgi; - amrex::TableData accrgc; - amrex::TableData evapg1; - amrex::TableData evapg2; - amrex::TableData evapr1; - amrex::TableData evapr2; - - // vertical plane average data - amrex::TableData rho1d; - amrex::TableData pres1d; - amrex::TableData tabs1d; - amrex::TableData qt1d; - amrex::TableData qv1d; - amrex::TableData qn1d; - - // independent variables - amrex::Array mic_fab_vars; - - amrex::TableData gamaz; - amrex::TableData zmid; // mid value of vertical coordinate in physical domain - - // data (output) - amrex::TableData qifall; - amrex::TableData tlatqi; +public: + + Microphysics ( ) { } + + ~Microphysics ( ) = default; + + void + SetModel (std::string moisture_model) + { + if ( moisture_model == "SAM" ) { + m_MicMod = MicMod::SAM; + m_sam_model = std::make_unique(); + } else if ( moisture_model == "SAM" ) { + // TODO: Code block + }else { + amrex::Print() << "WARNING: Compiled with moisture support but no valid moisture model found\n"; + m_MicMod = MicMod::Null; + m_null_model = std::make_unique(); + } + } + + void + define (SolverChoice& sc) + { + switch(m_MicMod) { + case MicMod::SAM: + m_sam_model->define(sc); + break; + case MicMod::Kessler: + // TODO: code block + break; + default: + m_null_model->define(sc); + } + } + + void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) + { + switch(m_MicMod) { + case MicMod::SAM: + m_sam_model->Init(cons_in, qmoist, grids, geom, dt_advance); + break; + case MicMod::Kessler: + // TODO: code block + break; + default: + m_null_model->Init(cons_in, qmoist, grids, geom, dt_advance); + } + } + + void + Advance ( ) + { + switch(m_MicMod) { + case MicMod::SAM: + m_sam_model->Advance(); + break; + case MicMod::Kessler: + // TODO: code block + break; + default: + m_null_model->Advance(); + } + } + + void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist) + { + switch(m_MicMod) { + case MicMod::SAM: + m_sam_model->Update(cons_in, qmoist); + break; + case MicMod::Kessler: + // TODO: code block + break; + default: + m_null_model->Update(cons_in, qmoist); + } + } + +private: + int m_MicMod; + std::unique_ptr m_null_model = nullptr; + std::unique_ptr m_sam_model = nullptr; + //std::unique_ptr m_kessler_model = nullptr; }; #endif +*/ diff --git a/Source/Microphysics/Null/Make.package b/Source/Microphysics/Null/Make.package new file mode 100644 index 000000000..d528e5844 --- /dev/null +++ b/Source/Microphysics/Null/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += NullMoist.H + diff --git a/Source/Microphysics/Null/Make.package~ b/Source/Microphysics/Null/Make.package~ new file mode 100644 index 000000000..36572d9e5 --- /dev/null +++ b/Source/Microphysics/Null/Make.package~ @@ -0,0 +1,2 @@ +CEXE_headers += Microphysics.H + diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H new file mode 100644 index 000000000..6361f00e7 --- /dev/null +++ b/Source/Microphysics/Null/NullMoist.H @@ -0,0 +1,35 @@ +#ifndef NULLMOIST_H +#define NULLMOIST_H + +#include +#include +#include + +class NullMoist { + + public: + NullMoist () {} + + virtual ~NullMoist () = default; + + virtual void + define (SolverChoice& /*sc*/) { }; + + virtual + void Init (const amrex::MultiFab& /*cons_in*/, + amrex::MultiFab& /*qmoist*/, + const amrex::BoxArray& /*grids*/, + const amrex::Geometry& /*geom*/, + const amrex::Real& /*dt_advance*/) { }; + + virtual void + Advance ( ) { }; + + virtual void + Update (amrex::MultiFab& /*cons_in*/, + amrex::MultiFab& /*qmoist*/) { }; + + private: + +}; +#endif diff --git a/Source/Microphysics/Null/NullMoist.H~ b/Source/Microphysics/Null/NullMoist.H~ new file mode 100644 index 000000000..4d1a3109c --- /dev/null +++ b/Source/Microphysics/Null/NullMoist.H~ @@ -0,0 +1,25 @@ +#ifndef NULLMOIST_H +#define NULLMOIST_H + +class NullMoist { + + public: + NullMoist () {} + + ~NullMoist () = default; + + void define (SolverChoice& /*sc*/) { }; + + void Init (const amrex::MultiFab& /*cons_in*/, + amrex::MultiFab& /*qmoist*/, + const amrex::BoxArray& /*grids*/, + const amrex::Geometry& /*geom*/, + const amrex::Real& /*dt_advance*/) { }; + + void Update (amrex::MultiFab& /*cons_in*/, + amrex::MultiFab& /*qmoist*/) { }; + + private: + +}; +#endif diff --git a/Source/Microphysics/SAM/Cloud.cpp b/Source/Microphysics/SAM/Cloud.cpp new file mode 100644 index 000000000..6ef1cc15b --- /dev/null +++ b/Source/Microphysics/SAM/Cloud.cpp @@ -0,0 +1,134 @@ + +#include "Microphysics.H" +#include "IndexDefines.H" +#include "TileNoZ.H" + +using namespace amrex; + +/** + * Compute Cloud-related Microphysics quantities. + */ +void SAM::Cloud () { + + constexpr Real an = 1.0/(tbgmax-tbgmin); + constexpr Real bn = tbgmin*an; + constexpr Real ap = 1.0/(tprmax-tprmin); + constexpr Real bp = tprmin*ap; + + auto pres1d_t = pres1d.table(); + + auto qt = mic_fab_vars[MicVar::qt]; + auto qp = mic_fab_vars[MicVar::qp]; + auto qn = mic_fab_vars[MicVar::qn]; + auto rho = mic_fab_vars[MicVar::rho]; + auto tabs = mic_fab_vars[MicVar::tabs]; + + Real fac_cond = m_fac_cond; + Real fac_sub = m_fac_sub; + Real fac_fus = m_fac_fus; + + for ( MFIter mfi(*tabs, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qt_array = qt->array(mfi); + auto qp_array = qp->array(mfi); + auto qn_array = qn->array(mfi); + auto tabs_array = tabs->array(mfi); + + const auto& box3d = mfi.tilebox() & m_gtoe[mfi.index()]; + + ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + qt_array(i,j,k) = std::max(0.0,qt_array(i,j,k)); + // Initial guess for temperature assuming no cloud water/ice: + Real tabs1 = tabs_array(i,j,k); + + Real qsatt; + Real om; + Real qsatt1; + Real qsatt2; + + // Warm cloud: + if(tabs1 > tbgmax) { + tabs1 = tabs_array(i,j,k)+fac_cond*qp_array(i,j,k); + erf_qsatw(tabs1, pres1d_t(k), qsatt); + } + // Ice cloud: + else if(tabs1 <= tbgmin) { + tabs1 = tabs_array(i,j,k)+fac_sub*qp_array(i,j,k); + erf_qsati(tabs1, pres1d_t(k), qsatt); + } + // Mixed-phase cloud: + else { + om = an*tabs1-bn; + erf_qsatw(tabs1, pres1d_t(k), qsatt1); + erf_qsati(tabs1, pres1d_t(k), qsatt2); + qsatt = om*qsatt1 + (1.-om)*qsatt2; + } +//if(i==2 && j==2) +// printf("cloud: %d, %d, %d, %13.6f, %13.6f, %13.6f, %13.6f, %13.6f, %13.6f\n", +// i,j,k,tabs1,pres1d_t(k),qt_array(i,j,k),qsatt,qn_array(i,j,k),qp_array(i,j,k)); + + int niter; + Real dtabs, lstarn, dlstarn, omp, lstarp, dlstarp, fff, dfff, dqsat; + // Test if condensation is possible: + if(qt_array(i,j,k) > qsatt) { + niter = 0; + dtabs = 1; + do { + if(tabs1 >= tbgmax) { + om=1.0; + lstarn = fac_cond; + dlstarn = 0.0; + erf_qsatw(tabs1, pres1d_t(k), qsatt); + erf_dtqsatw(tabs1, pres1d_t(k), dqsat); + } + else if(tabs1 <= tbgmin) { + om = 0.0; + lstarn = fac_sub; + dlstarn = 0.0; + erf_qsati(tabs1, pres1d_t(k), qsatt); + erf_dtqsati(tabs1, pres1d_t(k), dqsat); + } + else { + om=an*tabs1-bn; + lstarn = fac_cond+(1.0-om)*fac_fus; + dlstarn = an*fac_fus; + erf_qsatw(tabs1, pres1d_t(k), qsatt1); + erf_qsati(tabs1, pres1d_t(k), qsatt2); + + qsatt = om*qsatt1+(1.-om)*qsatt2; + erf_dtqsatw(tabs1, pres1d_t(k), qsatt1); + erf_dtqsati(tabs1, pres1d_t(k), qsatt2); + dqsat = om*qsatt1+(1.-om)*qsatt2; + } + + if(tabs1 >= tprmax) { + omp = 1.0; + lstarp = fac_cond; + dlstarp = 0.0; + } + else if(tabs1 <= tprmin) { + omp = 0.0; + lstarp = fac_sub; + dlstarp = 0.0; + } + else { + omp=ap*tabs1-bp; + lstarp = fac_cond+(1.0-omp)*fac_fus; + dlstarp = ap*fac_fus; + } + fff = tabs_array(i,j,k)-tabs1+lstarn*(qt_array(i,j,k)-qsatt)+lstarp*qp_array(i,j,k); + dfff = dlstarn*(qt_array(i,j,k)-qsatt)+dlstarp*qp_array(i,j,k)-lstarn*dqsat-1.0; + dtabs = -fff/dfff; + niter = niter+1; + tabs1 = tabs1+dtabs; + } while(std::abs(dtabs) > 0.01 && niter < 10); + qsatt = qsatt + dqsat*dtabs; + qn_array(i,j,k) = std::max(0.0, qt_array(i,j,k)-qsatt); + } + else { + qn_array(i,j,k) = 0.0; + } + tabs_array(i,j,k) = tabs1; + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); // just in case + }); + } +} diff --git a/Source/Microphysics/Cloud.cpp b/Source/Microphysics/SAM/Cloud.cpp~ similarity index 100% rename from Source/Microphysics/Cloud.cpp rename to Source/Microphysics/SAM/Cloud.cpp~ diff --git a/Source/Microphysics/SAM/Diagnose.cpp b/Source/Microphysics/SAM/Diagnose.cpp new file mode 100644 index 000000000..784f21a8e --- /dev/null +++ b/Source/Microphysics/SAM/Diagnose.cpp @@ -0,0 +1,59 @@ +#include "Microphysics.H" + +/** + * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid + * from the existing Microphysics variables. + */ +void SAM::Diagnose () { + + auto qt = mic_fab_vars[MicVar::qt]; + auto qp = mic_fab_vars[MicVar::qp]; + auto qv = mic_fab_vars[MicVar::qv]; + auto qn = mic_fab_vars[MicVar::qn]; + auto qcl = mic_fab_vars[MicVar::qcl]; + auto qci = mic_fab_vars[MicVar::qci]; + auto qpl = mic_fab_vars[MicVar::qpl]; + auto qpi = mic_fab_vars[MicVar::qpi]; + auto tabs = mic_fab_vars[MicVar::tabs]; + + for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qt_array = qt->array(mfi); + auto qp_array = qp->array(mfi); + auto qv_array = qv->array(mfi); + auto qn_array = qn->array(mfi); + auto qcl_array = qcl->array(mfi); + auto qci_array = qci->array(mfi); + auto qpl_array = qpl->array(mfi); + auto qpi_array = qpi->array(mfi); + auto tabs_array = tabs->array(mfi); + + const auto& box3d = mfi.tilebox(); + + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + qv_array(i,j,k) = qt_array(i,j,k) - qn_array(i,j,k); + amrex::Real omn = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); + qcl_array(i,j,k) = qn_array(i,j,k)*omn; + qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn); + amrex::Real omp = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + qpl_array(i,j,k) = qp_array(i,j,k)*omp; + qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp); + }); + } +} + +/** + * Wrapper for the PrecipFall, Cloud, Precipitation, and Diagnostics routines. + */ +void SAM::Proc () { + + MicroPrecipFall(); + + if (docloud) { + Cloud(); + if (doprecip) { + Precip(); + } + } + + Diagnose(); +} diff --git a/Source/Microphysics/Diagnose.cpp b/Source/Microphysics/SAM/Diagnose.cpp~ similarity index 100% rename from Source/Microphysics/Diagnose.cpp rename to Source/Microphysics/SAM/Diagnose.cpp~ diff --git a/Source/Microphysics/SAM/IceFall.cpp b/Source/Microphysics/SAM/IceFall.cpp new file mode 100644 index 000000000..e3635e284 --- /dev/null +++ b/Source/Microphysics/SAM/IceFall.cpp @@ -0,0 +1,171 @@ + +#include +#include "Microphysics.H" +#include "TileNoZ.H" + +using namespace amrex; + +/** + * Computes contributions to Microphysics and thermodynamic variables from falling cloud ice in each column. + */ +void SAM::IceFall () { + + Real dz = m_geom.CellSize(2); + Real dtn = dt; + int nz = nlev; + + int kmax, kmin; + auto qcl = mic_fab_vars[MicVar::qcl]; + auto qci = mic_fab_vars[MicVar::qci]; + auto qt = mic_fab_vars[MicVar::qt]; + auto tabs = mic_fab_vars[MicVar::tabs]; + auto theta = mic_fab_vars[MicVar::theta]; + + MultiFab fz; + fz.define(qcl->boxArray(),qcl->DistributionMap(), 1, qcl->nGrowVect()); + fz.setVal(0.); + + auto qifall_t = qifall.table(); + auto tlatqi_t = tlatqi.table(); + auto rho1d_t = rho1d.table(); + + kmin = nz; + kmax = -1; + + { // calculate maximum and minium ice fall vertical region + auto const& qcl_arrays = qcl->const_arrays(); + auto const& qci_arrays = qci->const_arrays(); + auto const& tabs_arrays = tabs->const_arrays(); + + GpuTuple k_max_min = ParReduce(TypeList{}, + TypeList{}, + *qcl, IntVect::TheZeroVector(), + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept + -> GpuTuple + { + bool is_positive = qcl_arrays[box_no](i,j,k)+qci_arrays[box_no](i,j,k) > 0.0; + bool smaller_than_zero = tabs_arrays[box_no](i,j,k) < 273.15; + int mkmin = is_positive && smaller_than_zero ? k : nz-1; + int mkmax = is_positive && smaller_than_zero ? k : 0; + return {mkmax, mkmin}; + }); + kmax = amrex::get<0>(k_max_min); + kmin = amrex::get<1>(k_max_min); + } + +//std::cout << "ice_fall: " << kmin << "; " << kmax << std::endl; + + // for (int k=0; karray(mfi); + auto qci_array = qci->array(mfi); + auto qt_array = qt->array(mfi); + //auto tabs_array = tabs->array(mfi); + auto theta_array = theta->array(mfi); + auto fz_array = fz.array(mfi); + + const auto& box3d = mfi.tilebox(); + + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (k >= std::max(0,kmin-1) && k <= kmax ) { + // Set up indices for x-y planes above and below current plane. + int kc = std::min(k+1, nz-1); + int kb = std::max(k-1, 0); + + // CFL number based on grid spacing interpolated to interface i,j,k-1/2 + Real coef = dtn/dz; //dtn/(0.5*(adz(kb)+adz(k))*dz); + + // Compute cloud ice density in this cell and the ones above/below. + // Since cloud ice is falling, the above cell is u(icrm,upwind), + // this cell is c (center) and the one below is d (downwind). + Real qiu = rho1d_t(kc)*qci_array(i,j,kc); + Real qic = rho1d_t(k )*qci_array(i,j,k ); + Real qid = rho1d_t(kb)*qci_array(i,j,kb); + + // Ice sedimentation velocity depends on ice content. The fiting is + // based on the data by Heymsfield (JAS,2003). -Marat + Real vt_ice = min( 0.4 , 8.66 * pow( (max(0.,qic)+1.e-10) , 0.24) ); // Heymsfield (JAS, 2003, p.2607) + + // Use MC flux limiter in computation of flux correction. + // (MC = monotonized centered difference). + // if (qic.eq.qid) then + Real tmp_phi; + if ( std::abs(qic-qid) < 1.0e-25 ) { // when qic, and qid is very small, qic_qid can still be zero + // even if qic is not equal to qid. so add a fix here +++mhwang + tmp_phi = 0.; + } else { + Real tmp_theta = (qiu-qic)/(qic-qid+1.0e-20); + tmp_phi = max(0., min(0.5*(1.+tmp_theta), min(2., 2.*tmp_theta))); + } + + // Compute limited flux. + // Since falling cloud ice is a 1D advection problem, this + // flux-limited advection scheme is monotonic. + fz_array(i,j,k) = -vt_ice*(qic - 0.5*(1.-coef*vt_ice)*tmp_phi*(qic-qid)); + } + }); + + // for (int j=0; j= std::max(0,kmin-1) && k <= kmax ) { + Real coef = dtn/dz; + // The cloud ice increment is the difference of the fluxes. + Real dqi = coef*(fz_array(i,j,k)-fz_array(i,j,k+1)); + // Add this increment to both non-precipitating and total water. + amrex::Gpu::Atomic::Add(&qt_array(i,j,k), dqi); + // Include this effect in the total moisture budget. + amrex::Gpu::Atomic::Add(&qifall_t(k), dqi); + + // The latent heat flux induced by the falling cloud ice enters + // the liquid-ice static energy budget in the same way as the + // precipitation. Note: use latent heat of sublimation. + Real lat_heat = (fac_cond+fac_fus)*dqi; + + // Add divergence of latent heat flux contribution to liquid-ice static potential temperature. + amrex::Gpu::Atomic::Add(&theta_array(i,j,k), -lat_heat); + // Add divergence to liquid-ice static energy budget. + amrex::Gpu::Atomic::Add(&tlatqi_t(k), -lat_heat); + } + }); + +#if 0 + // for (int j=0; j(ny,nx,ncrms) , YAKL_LAMBDA (int j, int i, int icrm) { + amrex::ParallelFor(box2d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + Real coef = dtn/dz; + Real dqi = -coef*fz(i,j,0); + precsfc (i,j) = precsfc (i,j)+dqi; + precssfc(i,j) = precssfc(i,j)+dqi; + }); +#endif + } +} + diff --git a/Source/Microphysics/IceFall.cpp b/Source/Microphysics/SAM/IceFall.cpp~ similarity index 100% rename from Source/Microphysics/IceFall.cpp rename to Source/Microphysics/SAM/IceFall.cpp~ diff --git a/Source/Microphysics/SAM/Init.cpp b/Source/Microphysics/SAM/Init.cpp new file mode 100644 index 000000000..385eff7f8 --- /dev/null +++ b/Source/Microphysics/SAM/Init.cpp @@ -0,0 +1,245 @@ + +#include +#include "Microphysics.H" +#include "IndexDefines.H" +#include "PlaneAverage.H" +#include "EOS.H" +#include "TileNoZ.H" + +using namespace amrex; + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @param[in] geom Geometry associated with these MultiFabs and grids + * @param[in] dt_advance Timestep for the advance + */ +void SAM::Init (const MultiFab& cons_in, MultiFab& qmoist, + const BoxArray& grids, + const Geometry& geom, + const Real& dt_advance) + { + m_geom = geom; + m_gtoe = grids; + + auto dz = m_geom.CellSize(2); + auto lowz = m_geom.ProbLo(2); + + dt = dt_advance; + + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar::NumVars; ++ivar) { + mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); + mic_fab_vars[ivar]->setVal(0.); + } + + // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled + // The ghost cells of these arrays aren't filled in the boundary condition calls for the state + qmoist.setVal(0.); + + for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { + + const auto& box3d = mfi.tilebox(); + + const auto& lo = amrex::lbound(box3d); + const auto& hi = amrex::ubound(box3d); + + nlev = box3d.length(2); + zlo = lo.z; + zhi = hi.z; + + // parameters + accrrc.resize({zlo}, {zhi}); + accrsi.resize({zlo}, {zhi}); + accrsc.resize({zlo}, {zhi}); + coefice.resize({zlo}, {zhi}); + evaps1.resize({zlo}, {zhi}); + evaps2.resize({zlo}, {zhi}); + accrgi.resize({zlo}, {zhi}); + accrgc.resize({zlo}, {zhi}); + evapg1.resize({zlo}, {zhi}); + evapg2.resize({zlo}, {zhi}); + evapr1.resize({zlo}, {zhi}); + evapr2.resize({zlo}, {zhi}); + + // data (input) + rho1d.resize({zlo}, {zhi}); + pres1d.resize({zlo}, {zhi}); + tabs1d.resize({zlo}, {zhi}); + gamaz.resize({zlo}, {zhi}); + zmid.resize({zlo}, {zhi}); + + // output + qifall.resize({zlo}, {zhi}); + tlatqi.resize({zlo}, {zhi}); + } + + auto accrrc_t = accrrc.table(); + auto accrsi_t = accrsi.table(); + auto accrsc_t = accrsc.table(); + auto coefice_t = coefice.table(); + auto evaps1_t = evaps1.table(); + auto evaps2_t = evaps2.table(); + auto accrgi_t = accrgi.table(); + auto accrgc_t = accrgc.table(); + auto evapg1_t = evapg1.table(); + auto evapg2_t = evapg2.table(); + auto evapr1_t = evapr1.table(); + auto evapr2_t = evapr2.table(); + + auto rho1d_t = rho1d.table(); + auto pres1d_t = pres1d.table(); + auto tabs1d_t = tabs1d.table(); + + auto gamaz_t = gamaz.table(); + auto zmid_t = zmid.table(); + + Real gam3 = erf_gammafff(3.0 ); + Real gamr1 = erf_gammafff(3.0+b_rain ); + Real gamr2 = erf_gammafff((5.0+b_rain)/2.0); + // Real gamr3 = erf_gammafff(4.0+b_rain ); + Real gams1 = erf_gammafff(3.0+b_snow ); + Real gams2 = erf_gammafff((5.0+b_snow)/2.0); + // Real gams3 = erf_gammafff(4.0+b_snow ); + Real gamg1 = erf_gammafff(3.0+b_grau ); + Real gamg2 = erf_gammafff((5.0+b_grau)/2.0); + // Real gamg3 = erf_gammafff(4.0+b_grau ); + + amrex::MultiFab qv(qmoist, amrex::make_alias, 0, 1); + amrex::MultiFab qc(qmoist, amrex::make_alias, 1, 1); + amrex::MultiFab qi(qmoist, amrex::make_alias, 2, 1); + + // Get the temperature, density, theta, qt and qp from input + for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { + auto states_array = cons_in.array(mfi); + auto qc_array = qc.array(mfi); + auto qi_array = qi.array(mfi); + + auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi); + auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi); + auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi); + auto temp_array = mic_fab_vars[MicVar::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + // Get pressure, theta, temperature, density, and qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + rho_array(i,j,k) = states_array(i,j,k,Rho_comp); + theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); + qt_array(i,j,k) = states_array(i,j,k,RhoQt_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQp_comp)/states_array(i,j,k,Rho_comp); + qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); + temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp)); + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; + }); + } + + // calculate the plane average variables + PlaneAverage cons_ave(&cons_in, m_geom, m_axis); + cons_ave.compute_averages(ZDir(), cons_ave.field()); + + // get host variable rho, and rhotheta + int ncell = cons_ave.ncell_line(); + + Gpu::HostVector rho_h(ncell), rhotheta_h(ncell); + cons_ave.line_average(Rho_comp, rho_h); + cons_ave.line_average(RhoTheta_comp, rhotheta_h); + + // copy data to device + Gpu::DeviceVector rho_d(ncell), rhotheta_d(ncell); + Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, rhotheta_h.begin(), rhotheta_h.end(), rhotheta_d.begin()); + Gpu::streamSynchronize(); + + Real* rho_dptr = rho_d.data(); + Real* rhotheta_dptr = rhotheta_d.data(); + + Real gOcp = m_gOcp; + + amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { + Real pressure = getPgivenRTh(rhotheta_dptr[k]); + rho1d_t(k) = rho_dptr[k]; + pres1d_t(k) = pressure/100.; + tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],rhotheta_dptr[k]); + zmid_t(k) = lowz + (k+0.5)*dz; + gamaz_t(k) = gOcp*zmid_t(k); + }); + + // This fills qv + Diagnose(); + +#if 0 + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) { + fluxbmk(l,j,i) = 0.0; + fluxtmk(l,j,i) = 0.0; + }); + + amrex::ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int k, int l, int i0) { + mkwle (l,k) = 0.0; + mkwsb (l,k) = 0.0; + mkadv (l,k) = 0.0; + mkdiff(l,k) = 0.0; + }); +#endif + + if(round(gam3) != 2) { + std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; + std::exit(-1); + } + + amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { + + Real pratio = sqrt(1.29 / rho1d_t(k)); + Real rrr1=393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5); + Real rrr2=std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k)); + Real estw = 100.0*erf_esatw(tabs1d_t(k)); + Real esti = 100.0*erf_esati(tabs1d_t(k)); + + // accretion by snow: + Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0)); + Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); + accrsi_t(k) = coef1 * coef2 * esicoef; + accrsc_t(k) = coef1 * esccoef; + coefice_t(k) = coef2; + + // evaporation of snow: + coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); + evaps1_t(k) = 0.65*4.0*nzeros/sqrt(PI*rhos*nzeros)/(coef1+coef2)/sqrt(rho1d_t(k)); + evaps2_t(k) = 0.49*4.0*nzeros*gams2*sqrt(a_snow/(muelq*rrr1))/pow((PI*rhos*nzeros) , ((5.0+b_snow)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_snow)/8.0))*sqrt(pratio); + + // accretion by graupel: + coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0)); + coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); + accrgi_t(k) = coef1 * coef2 * egicoef; + accrgc_t(k) = coef1 * egccoef; + + // evaporation of graupel: + coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); + evapg1_t(k) = 0.65*4.0*nzerog/sqrt(PI*rhog*nzerog)/(coef1+coef2)/sqrt(rho1d_t(k)); + evapg2_t(k) = 0.49*4.0*nzerog*gamg2*sqrt(a_grau/(muelq*rrr1))/pow((PI * rhog * nzerog) , ((5.0+b_grau)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_grau)/8.0))*sqrt(pratio); + + // accretion by rain: + accrrc_t(k)= 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3+b_rain)/4.))* erccoef; + + // evaporation of rain: + coef1 =(lcond/(tabs1d_t(k)*rv)-1.0)*lcond/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq * rrr2 * estw); + evapr1_t(k) = 0.78 * 2.0 * PI * nzeror / + sqrt(PI * rhor * nzeror) / (coef1+coef2) / sqrt(rho1d_t(k)); + evapr2_t(k) = 0.31 * 2.0 * PI * nzeror * gamr2 * 0.89 * sqrt(a_rain/(muelq*rrr1))/ + pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); + }); +} diff --git a/Source/Microphysics/Init.cpp b/Source/Microphysics/SAM/Init.cpp~ similarity index 100% rename from Source/Microphysics/Init.cpp rename to Source/Microphysics/SAM/Init.cpp~ diff --git a/Source/Microphysics/SAM/Make.package b/Source/Microphysics/SAM/Make.package new file mode 100644 index 000000000..f55996515 --- /dev/null +++ b/Source/Microphysics/SAM/Make.package @@ -0,0 +1,9 @@ +CEXE_sources += Cloud.cpp +CEXE_sources += Init.cpp +CEXE_sources += Diagnose.cpp +CEXE_sources += Update.cpp +CEXE_sources += IceFall.cpp +CEXE_sources += Precip.cpp +CEXE_sources += PrecipFall.cpp +CEXE_headers += SAM.H + diff --git a/Source/Microphysics/SAM/Precip.cpp b/Source/Microphysics/SAM/Precip.cpp new file mode 100644 index 000000000..aa58a6ee4 --- /dev/null +++ b/Source/Microphysics/SAM/Precip.cpp @@ -0,0 +1,148 @@ +/* + * this file is modified from precip_proc from samxx + */ +#include "Microphysics.H" + +using namespace amrex; + +/** + * Compute Precipitation-related Microphysics quantities. + */ +void SAM::Precip () { + + Real powr1 = (3.0 + b_rain) / 4.0; + Real powr2 = (5.0 + b_rain) / 8.0; + Real pows1 = (3.0 + b_snow) / 4.0; + Real pows2 = (5.0 + b_snow) / 8.0; + Real powg1 = (3.0 + b_grau) / 4.0; + Real powg2 = (5.0 + b_grau) / 8.0; + + auto accrrc_t = accrrc.table(); + auto accrsc_t = accrsc.table(); + auto accrsi_t = accrsi.table(); + auto accrgc_t = accrgc.table(); + auto accrgi_t = accrgi.table(); + auto coefice_t = coefice.table(); + auto evapr1_t = evapr1.table(); + auto evapr2_t = evapr2.table(); + auto evaps1_t = evaps1.table(); + auto evaps2_t = evaps2.table(); + auto evapg1_t = evapg1.table(); + auto evapg2_t = evapg2.table(); + auto pres1d_t = pres1d.table(); + + auto qt = mic_fab_vars[MicVar::qt]; + auto qp = mic_fab_vars[MicVar::qp]; + auto qn = mic_fab_vars[MicVar::qn]; + auto tabs = mic_fab_vars[MicVar::tabs]; + + Real dtn = dt; + + // get the temperature, dentisy, theta, qt and qp from input + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi); + auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi); + auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + //------- Autoconversion/accretion + Real omn, omp, omg, qcc, qii, autor, autos, accrr, qrr, accrcs, accris, + qss, accrcg, accrig, tmp, qgg, dq, qsatt, qsat; + + if (qn_array(i,j,k)+qp_array(i,j,k) > 0.0) { + omn = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); + omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); + + if (qn_array(i,j,k) > 0.0) { + qcc = qn_array(i,j,k) * omn; + qii = qn_array(i,j,k) * (1.0-omn); + + if (qcc > qcw0) { + autor = alphaelq; + } else { + autor = 0.0; + } + + if (qii > qci0) { + autos = betaelq*coefice_t(k); + } else { + autos = 0.0; + } + + accrr = 0.0; + if (omp > 0.001) { + qrr = qp_array(i,j,k) * omp; + accrr = accrrc_t(k) * std::pow(qrr, powr1); + } + + accrcs = 0.0; + accris = 0.0; + + if (omp < 0.999 && omg < 0.999) { + qss = qp_array(i,j,k) * (1.0-omp)*(1.0-omg); + tmp = pow(qss, pows1); + accrcs = accrsc_t(k) * tmp; + accris = accrsi_t(k) * tmp; + } + accrcg = 0.0; + accrig = 0.0; + if (omp < 0.999 && omg > 0.001) { + qgg = qp_array(i,j,k) * (1.0-omp)*omg; + tmp = pow(qgg, powg1); + accrcg = accrgc_t(k) * tmp; + accrig = accrgi_t(k) * tmp; + } + qcc = (qcc+dtn*autor*qcw0)/(1.0+dtn*(accrr+accrcs+accrcg+autor)); + qii = (qii+dtn*autos*qci0)/(1.0+dtn*(accris+accrig+autos)); + dq = dtn *(accrr*qcc + autor*(qcc-qcw0)+(accris+accrig)*qii + (accrcs+accrcg)*qcc + autos*(qii-qci0)); + dq = std::min(dq,qn_array(i,j,k)); + qt_array(i,j,k) = qt_array(i,j,k) - dq; + qp_array(i,j,k) = qp_array(i,j,k) + dq; + qn_array(i,j,k) = qn_array(i,j,k) - dq; + + } else if(qp_array(i,j,k) > qp_threshold && qn_array(i,j,k) == 0.0) { + + qsatt = 0.0; + if(omn > 0.001) { + erf_qsatw(tabs_array(i,j,k),pres1d_t(k),qsat); + qsatt = qsatt + omn*qsat; + } + if(omn < 0.999) { + erf_qsati(tabs_array(i,j,k),pres1d_t(k),qsat); + qsatt = qsatt + (1.-omn)*qsat; + } + dq = 0.0; + if(omp > 0.001) { + qrr = qp_array(i,j,k) * omp; + dq = dq + evapr1_t(k)*sqrt(qrr) + evapr2_t(k)*pow(qrr,powr2); + } + if(omp < 0.999 && omg < 0.999) { + qss = qp_array(i,j,k) * (1.0-omp)*(1.0-omg); + dq = dq + evaps1_t(k)*sqrt(qss) + evaps2_t(k)*pow(qss,pows2); + } + if(omp < 0.999 && omg > 0.001) { + qgg = qp_array(i,j,k) * (1.0-omp)*omg; + dq = dq + evapg1_t(k)*sqrt(qgg) + evapg2_t(k)*pow(qgg,powg2); + } + dq = dq * dtn * (qt_array(i,j,k) / qsatt-1.0); + dq = std::max(-0.5*qp_array(i,j,k),dq); + qt_array(i,j,k) = qt_array(i,j,k) - dq; + qp_array(i,j,k) = qp_array(i,j,k) + dq; + + } else { + qt_array(i,j,k) = qt_array(i,j,k) + qp_array(i,j,k); + qp_array(i,j,k) = 0.0; + } + } + dq = qp_array(i,j,k); + qp_array(i,j,k) = max(0.0,qp_array(i,j,k)); + qt_array(i,j,k) = qt_array(i,j,k) + (dq-qp_array(i,j,k)); + }); + } +} + + diff --git a/Source/Microphysics/Precip.cpp b/Source/Microphysics/SAM/Precip.cpp~ similarity index 100% rename from Source/Microphysics/Precip.cpp rename to Source/Microphysics/SAM/Precip.cpp~ diff --git a/Source/Microphysics/SAM/PrecipFall.cpp b/Source/Microphysics/SAM/PrecipFall.cpp new file mode 100644 index 000000000..86cf8daf3 --- /dev/null +++ b/Source/Microphysics/SAM/PrecipFall.cpp @@ -0,0 +1,307 @@ +#include "ERF_Constants.H" +#include "Microphysics.H" +#include "TileNoZ.H" + +using namespace amrex; + +/** + * Compute positive definite monotonic advection with a non-oscillatory option + * and gravitational sedimentation. + * + * Code modified from SAMXX, the C++ version of the SAM code. + * + * @param[in] hydro_type Type selection for the precipitation advection hydrodynamics scheme (0-3) + */ +void SAM::PrecipFall(int hydro_type) { + + Real constexpr eps = 1.e-10; + bool constexpr nonos = true; + +/* + Real gam3 = erf_gammafff(3.0 ); + Real gamr1 = erf_gammafff(3.0+b_rain ); + Real gamr2 = erf_gammafff((5.0+b_rain)/2.0); + Real gams1 = erf_gammafff(3.0+b_snow ); + Real gams2 = erf_gammafff((5.0+b_snow)/2.0); + Real gamg1 = erf_gammafff(3.0+b_grau ); + Real gamg2 = erf_gammafff((5.0+b_grau)/2.0); +*/ + + Real gamr3 = erf_gammafff(4.0+b_rain ); + Real gams3 = erf_gammafff(4.0+b_snow ); + Real gamg3 = erf_gammafff(4.0+b_grau ); + + Real vrain = a_rain*gamr3/6.0/pow((PI*rhor*nzeror),crain); + Real vsnow = a_snow*gams3/6.0/pow((PI*rhos*nzeros),csnow); + Real vgrau = a_grau*gamg3/6.0/pow((PI*rhog*nzerog),cgrau); + + Real dt_advance = dt; + int nz = nlev; + + auto qp = mic_fab_vars[MicVar::qp]; + auto omega = mic_fab_vars[MicVar::omega]; + auto tabs = mic_fab_vars[MicVar::tabs]; + auto theta = mic_fab_vars[MicVar::theta]; + + auto ba = tabs->boxArray(); + auto dm = tabs->DistributionMap(); + auto ngrow = tabs->nGrowVect(); + + MultiFab mx; + MultiFab mn; + MultiFab lfac; + MultiFab www; + MultiFab fz; + MultiFab wp; + MultiFab tmp_qp; + + mx.define(ba,dm, 1, ngrow); + mn.define(ba,dm, 1, ngrow); + lfac.define(ba, dm, 1, ngrow); + www.define(ba, dm, 1, ngrow); + fz.define(ba, dm, 1, ngrow); + wp.define(ba, dm, 1, ngrow); + tmp_qp.define(ba, dm, 1, ngrow); + + TableData irho; + TableData iwmax; + TableData rhofac; + + irho.resize({zlo},{zhi}); + iwmax.resize({zlo},{zhi}); + rhofac.resize({zlo},{zhi}); + + auto irho_t = irho.table(); + auto iwmax_t = iwmax.table(); + auto rhofac_t = rhofac.table(); + auto rho1d_t = rho1d.table(); + + auto dz = m_geom.CellSize(2); + + ParallelFor(nz, [=] AMREX_GPU_DEVICE (int k) noexcept { + rhofac_t(k) = std::sqrt(1.29/rho1d_t(k)); + irho_t(k) = 1.0/rho1d_t(k); + Real wmax = dz/dt_advance; // Velocity equivalent to a cfl of 1.0. + iwmax_t(k) = 1.0/wmax; + }); + + // Add sedimentation of precipitation field to the vert. vel. + MultiFab prec_cfl_fab; + prec_cfl_fab.define(tabs->boxArray(),tabs->DistributionMap(), 1, tabs->nGrowVect()); + + for ( MFIter mfi(lfac, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto lfac_array = lfac.array(mfi); + auto omega_array = omega->array(mfi); + auto qp_array = qp->array(mfi); + auto tabs_array = tabs->array(mfi); + auto wp_array = wp.array(mfi); + auto www_array = www.array(mfi); + auto fz_array = fz.array(mfi); + auto prec_cfl_array = prec_cfl_fab.array(mfi); + + const auto& box3d = mfi.tilebox(); + + const Real fac_cond = m_fac_cond; + const Real fac_sub = m_fac_sub; + const Real fac_fus = m_fac_fus; + + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (hydro_type == 0) { + lfac_array(i,j,k) = fac_cond; + } + else if (hydro_type == 1) { + lfac_array(i,j,k) = fac_sub; + } + else if (hydro_type == 2) { + lfac_array(i,j,k) = fac_cond + (1.0-omega_array(i,j,k))*fac_fus; + } + else if (hydro_type == 3) { + lfac_array(i,j,k) = 0.0; + } + Real tmp = term_vel_qp(i,j,k,qp_array(i,j,k), + vrain, vsnow, vgrau, rho1d_t(k), + tabs_array(i,j,k)); + wp_array(i,j,k)=rhofac_t(k)*tmp; + tmp = wp_array(i,j,k)*iwmax_t(k); + prec_cfl_array(i,j,k) = tmp; + wp_array(i,j,k) = -wp_array(i,j,k)*rho1d_t(k)*dt_advance/dz; + if (k == 0) { + fz_array(i,j,nz-1) = 0.0; + www_array(i,j,nz-1) = 0.0; + lfac_array(i,j,nz-1) = 0.0; + } + }); + } + + auto const& cfl_arrays = prec_cfl_fab.const_arrays(); + Real prec_cfl = ParReduce(TypeList{}, TypeList{}, + prec_cfl_fab, IntVect(0), + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept + ->GpuTuple + { + return { cfl_arrays[box_no](i,j,k) }; + }); + + // If maximum CFL due to precipitation velocity is greater than 0.9, + // take more than one advection step to maintain stability. + int nprec; + if (prec_cfl > 0.9) { + nprec = static_cast(std::ceil(prec_cfl/0.9)); + for (MFIter mfi(wp, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto wp_array = wp.array(mfi); + const auto& box3d = mfi.tilebox(); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) { + // wp already includes factor of dt, so reduce it by a + // factor equal to the number of precipitation steps. + wp_array(i,j,k) = wp_array(i,j,k)/Real(nprec); + }); + } + } else { + nprec = 1; + } + +//std::cout << "precipfall: nprec= " << nprec << std::endl; + +#ifdef ERF_FIXED_SUBCYCLE + nprec = 4; +#endif + + for(int iprec = 1; iprec<=nprec; iprec++) { + for ( MFIter mfi(tmp_qp, TileNoZ()); mfi.isValid(); ++mfi) { + auto qp_array = qp->array(mfi); + auto tabs_array = tabs->array(mfi); + auto theta_array = theta->array(mfi); + auto tmp_qp_array = tmp_qp.array(mfi); + auto mx_array = mx.array(mfi); + auto mn_array = mn.array(mfi); + auto fz_array = fz.array(mfi); + auto wp_array = wp.array(mfi); + auto lfac_array = lfac.array(mfi); + auto www_array = www.array(mfi); + + const auto& box3d = mfi.tilebox(); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + tmp_qp_array(i,j,k) = qp_array(i,j,k); // Temporary array for qp in this column + }); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (nonos) { + int kc=min(nz-1,k+1); + int kb=max(0,k-1); + mx_array(i,j,k) = max(tmp_qp_array(i,j,kb), max(tmp_qp_array(i,j,kc), tmp_qp_array(i,j,k))); + mn_array(i,j,k) = min(tmp_qp_array(i,j,kb), min(tmp_qp_array(i,j,kc), tmp_qp_array(i,j,k))); + } + // Define upwind precipitation flux + fz_array(i,j,k) = tmp_qp_array(i,j,k)*wp_array(i,j,k); + }); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int kc = min(k+1, nz-1); + tmp_qp_array(i,j,k) = tmp_qp_array(i,j,k)-(fz_array(i,j,kc)-fz_array(i,j,k))*irho_t(k); //Update temporary qp + }); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + // Also, compute anti-diffusive correction to previous + // (upwind) approximation to the flux + int kb=max(0,k-1); + // The precipitation velocity is a cell-centered quantity, + // since it is computed from the cell-centered + // precipitation mass fraction. Therefore, a reformulated + // anti-diffusive flux is used here which accounts for + // this and results in reduced numerical diffusion. + www_array(i,j,k) = 0.5*(1.0+wp_array(i,j,k)*irho_t(k))*(tmp_qp_array(i,j,kb)*wp_array(i,j,kb) - + tmp_qp_array(i,j,k)*wp_array(i,j,k)); // works for wp(k)<0 + }); + + if (nonos) { + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int kc=min(nz-1,k+1); + int kb=max(0,k-1); + mx_array(i,j,k) = max(tmp_qp_array(i,j,kb),max(tmp_qp_array(i,j,kc), max(tmp_qp_array(i,j,k), mx_array(i,j,k)))); + mn_array(i,j,k) = min(tmp_qp_array(i,j,kb),min(tmp_qp_array(i,j,kc), min(tmp_qp_array(i,j,k), mn_array(i,j,k)))); + kc = min(nz-1,k+1); + mx_array(i,j,k) = rho1d_t(k)*(mx_array(i,j,k)-tmp_qp_array(i,j,k))/(pn(www_array(i,j,kc)) + + pp(www_array(i,j,k))+eps); + mn_array(i,j,k) = rho1d_t(k)*(tmp_qp_array(i,j,k)-mn_array(i,j,k))/(pp(www_array(i,j,kc)) + + pn(www_array(i,j,k))+eps); + }); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int kb=max(0,k-1); + // Add limited flux correction to fz(k). + fz_array(i,j,k) = fz_array(i,j,k) + pp(www_array(i,j,k))*std::min(1.0,std::min(mx_array(i,j,k), mn_array(i,j,kb))) - + pn(www_array(i,j,k))*std::min(1.0,std::min(mx_array(i,j,kb),mn_array(i,j,k))); // Anti-diffusive flux + }); + } + + // Update precipitation mass fraction and liquid-ice static + // energy using precipitation fluxes computed in this column. + ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int kc=min(k+1, nz-1); + // Update precipitation mass fraction. + // Note that fz is the total flux, including both the + // upwind flux and the anti-diffusive correction. +// Real flagstat = 1.0; + qp_array(i,j,k) = qp_array(i,j,k) - (fz_array(i,j,kc) - fz_array(i,j,k))*irho_t(k); +// Real tmp = -(fz_array(i,j,kc)-fz_array(i,j,k))*irho_t(k)*flagstat; // For qp budget +// +// NOTE qpfall, tlat, and precflux,...are output diagnostic variables, not sure whether we need to calculate these variables here? +// Please correct me!!! by xyuan@anl.gov +// +// amrex::Gpu::Atomic::Add(&qpfall_t(k), tmp); + Real lat_heat = -(lfac_array(i,j,kc)*fz_array(i,j,kc)-lfac_array(i,j,k)*fz_array(i,j,k))*irho_t(k); + amrex::Gpu::Atomic::Add(&theta_array(i,j,k), -lat_heat); +// amrex::Gpu::Atomic::Add(&tlat_t(k), -lat_heat); +// tmp = fz_array(i,j,k)*flagstat; +// amrex::Gpu::Atomic::Add(&precflux_t(k), -tmp); +// yakl::atomicAdd(precflux(k,icrm),-tmp); + if (k == 0) { +// precsfc_t(i,i) = precsfc_t(i,j) - fz_array(i,j,0)*flagstat; +// precssfc_t(i,j) = precssfc_t(i,j) - fz_array(i,j,0)*(1.0-omega_array(i,j,0))*flagstat; +// prec_xy_t(i,j) = prec_xy_t(i,j) - fz_array(i,j,0)*flagstat; + } + }); + + if (iprec < nprec) { + // Re-compute precipitation velocity using new value of qp. + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real tmp = term_vel_qp(i,j,k,qp_array(i,j,k), + vrain, vsnow, vgrau, rho1d_t(k), + tabs_array(i,j,k)); + wp_array(i,j,k) = rhofac_t(k)*tmp; + // Decrease precipitation velocity by factor of nprec + wp_array(i,j,k) = -wp_array(i,j,k)*rho1d_t(k)*dt_advance/dz/nprec; + // Note: Don't bother checking CFL condition at each + // substep since it's unlikely that the CFL will + // increase very much between substeps when using + // monotonic advection schemes. + if (k == 0) { + fz_array(i,j,nz-1) = 0.0; + www_array(i,j,nz-1) = 0.0; + lfac_array(i,j,nz-1) = 0.0; + } + }); + } + } + } // iprec loop +} + +/** + * Wrapper for PrecipFall which computes the temporary variable Omega, needed by the precipitation advection scheme. + */ +void SAM::MicroPrecipFall() { + + for ( MFIter mfi(*(mic_fab_vars[MicVar::omega]), TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto omega_array = mic_fab_vars[MicVar::omega]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + omega_array(i,j,k) = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + }); + } + PrecipFall(2); +} diff --git a/Source/Microphysics/PrecipFall.cpp b/Source/Microphysics/SAM/PrecipFall.cpp~ similarity index 100% rename from Source/Microphysics/PrecipFall.cpp rename to Source/Microphysics/SAM/PrecipFall.cpp~ diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H new file mode 100644 index 000000000..c1d785cfc --- /dev/null +++ b/Source/Microphysics/SAM/SAM.H @@ -0,0 +1,177 @@ +/* + * Implementation 1-moment microphysics model + * NOTE: this model is based on the SAM code, and the Klemp's paper + * 1): Joseph, Klemp, the simulation of three-dimensional convective storm dynamics, + * Journal of the atmospheric sciences, vol35, p1070 + * 2): Marat Khairoutdinov and David Randall, cloud resolving modeling of the ARM summer 1997 IOP: + * model formulation, results, unvertainties, and sensitivities, Journal of the atmospheric sciences, vol60, p607 + */ +#ifndef SAM_H +#define SAM_H + +#include +#include +#include + +#include +#include +#include +#include + +#include "ERF_Constants.H" +#include "Microphysics_Utils.H" +#include "IndexDefines.H" +#include "DataStruct.H" + +namespace MicVar { + enum { + // independent variables + qt = 0, + qp, + theta, // liquid/ice water potential temperature + tabs, // temperature + rho, // density + pres, // pressure + // derived variables + qr, // rain water + qv, // water vapor + qn, // cloud condensate (liquid+ice), initial to zero + qci, // cloud ice + qcl, // cloud water + qpl, // precip rain + qpi, // precip ice + // temporary variable + omega, + NumVars + }; +} + +// +// use MultiFab for 3D data, but table for 1D data +// +class SAM { + + using FabPtr = std::shared_ptr; + + public: + // constructor + SAM () {} + + // destructor + virtual ~SAM () = default; + + // Set up for first time + virtual void + define (SolverChoice& sc) + { + docloud = sc.do_cloud; + doprecip = sc.do_precip; + m_fac_cond = lcond / sc.c_p; + m_fac_fus = lfus / sc.c_p; + m_fac_sub = lsub / sc.c_p; + m_gOcp = CONST_GRAV / sc.c_p; + m_axis = sc.ave_plane; + } + + // cloud physics + void Cloud (); + + // ice physics + void IceFall (); + + // precip + void Precip (); + + // precip fall + void PrecipFall (int hydro_type); + + // micro interface for precip fall + void MicroPrecipFall (); + + // diagnose + void Diagnose (); + + // process microphysics + void Proc (); + + // init + virtual void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance); + + // update ERF variables + virtual void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist); + + // wrapper to do all the updating + virtual void + Advance ( ) + { + this->Cloud(); + this->Diagnose(); + this->IceFall(); + this->Precip(); + this->MicroPrecipFall(); + } + + private: + // geometry + amrex::Geometry m_geom; + // valid boxes on which to evolve the solution + amrex::BoxArray m_gtoe; + + // timestep + amrex::Real dt; + + // number of vertical levels + int nlev, zlo, zhi; + + // plane average axis + int m_axis; + + // model options + bool docloud, doprecip; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // microphysics parameters/coefficients + amrex::TableData accrrc; + amrex::TableData accrsi; + amrex::TableData accrsc; + amrex::TableData coefice; + amrex::TableData evaps1; + amrex::TableData evaps2; + amrex::TableData accrgi; + amrex::TableData accrgc; + amrex::TableData evapg1; + amrex::TableData evapg2; + amrex::TableData evapr1; + amrex::TableData evapr2; + + // vertical plane average data + amrex::TableData rho1d; + amrex::TableData pres1d; + amrex::TableData tabs1d; + amrex::TableData qt1d; + amrex::TableData qv1d; + amrex::TableData qn1d; + + // independent variables + amrex::Array mic_fab_vars; + + amrex::TableData gamaz; + amrex::TableData zmid; // mid value of vertical coordinate in physical domain + + // data (output) + amrex::TableData qifall; + amrex::TableData tlatqi; +}; +#endif diff --git a/Source/Microphysics/SAM/SAM.H~ b/Source/Microphysics/SAM/SAM.H~ new file mode 100644 index 000000000..09e538320 --- /dev/null +++ b/Source/Microphysics/SAM/SAM.H~ @@ -0,0 +1,163 @@ +/* + * Implementation 1-moment microphysics model + * NOTE: this model is based on the SAM code, and the Klemp's paper + * 1): Joseph, Klemp, the simulation of three-dimensional convective storm dynamics, + * Journal of the atmospheric sciences, vol35, p1070 + * 2): Marat Khairoutdinov and David Randall, cloud resolving modeling of the ARM summer 1997 IOP: + * model formulation, results, unvertainties, and sensitivities, Journal of the atmospheric sciences, vol60, p607 + */ +#ifndef MICROPHYSICS_H +#define MICROPHYSICS_H + +#include +#include +#include + +#include +#include +#include +#include + +#include "ERF_Constants.H" +#include "Microphysics_Utils.H" +#include "IndexDefines.H" +#include "DataStruct.H" + +namespace MicVar { + enum { + // independent variables + qt = 0, + qp, + theta, // liquid/ice water potential temperature + tabs, // temperature + rho, // density + pres, // pressure + // derived variables + qr, // rain water + qv, // water vapor + qn, // cloud condensate (liquid+ice), initial to zero + qci, // cloud ice + qcl, // cloud water + qpl, // precip rain + qpi, // precip ice + // temporary variable + omega, + NumVars + }; +} + +// +// use MultiFab for 3D data, but table for 1D data +// +class Microphysics { + + using FabPtr = std::shared_ptr; + + public: + // constructor + Microphysics () {} + + // Set up for first time + void define (SolverChoice& sc) + { + docloud = sc.do_cloud; + doprecip = sc.do_precip; + m_fac_cond = lcond / sc.c_p; + m_fac_fus = lfus / sc.c_p; + m_fac_sub = lsub / sc.c_p; + m_gOcp = CONST_GRAV / sc.c_p; + m_axis = sc.ave_plane; + } + + // destructor + ~Microphysics () = default; + + // cloud physics + void Cloud (); + + // ice physics + void IceFall (); + + // precip + void Precip (); + + // precip fall + void PrecipFall (int hydro_type); + + // micro interface for precip fall + void MicroPrecipFall (); + + // init + void Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance); + + // update ERF variables + void Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist); + + // diagnose + void Diagnose (); + + // process microphysics + void Proc (); + + private: + // geometry + amrex::Geometry m_geom; + // valid boxes on which to evolve the solution + amrex::BoxArray m_gtoe; + + // timestep + amrex::Real dt; + + // number of vertical levels + int nlev, zlo, zhi; + + // plane average axis + int m_axis; + + // model options + bool docloud, doprecip; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // microphysics parameters/coefficients + amrex::TableData accrrc; + amrex::TableData accrsi; + amrex::TableData accrsc; + amrex::TableData coefice; + amrex::TableData evaps1; + amrex::TableData evaps2; + amrex::TableData accrgi; + amrex::TableData accrgc; + amrex::TableData evapg1; + amrex::TableData evapg2; + amrex::TableData evapr1; + amrex::TableData evapr2; + + // vertical plane average data + amrex::TableData rho1d; + amrex::TableData pres1d; + amrex::TableData tabs1d; + amrex::TableData qt1d; + amrex::TableData qv1d; + amrex::TableData qn1d; + + // independent variables + amrex::Array mic_fab_vars; + + amrex::TableData gamaz; + amrex::TableData zmid; // mid value of vertical coordinate in physical domain + + // data (output) + amrex::TableData qifall; + amrex::TableData tlatqi; +}; +#endif diff --git a/Source/Microphysics/SAM/Update.cpp b/Source/Microphysics/SAM/Update.cpp new file mode 100644 index 000000000..3c99d2802 --- /dev/null +++ b/Source/Microphysics/SAM/Update.cpp @@ -0,0 +1,60 @@ + +#include "Microphysics.H" +#include "IndexDefines.H" +#include "TileNoZ.H" + +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc, qi, qr, qs, qg + */ +void SAM::Update (amrex::MultiFab& cons, + amrex::MultiFab& qmoist) +{ + // copy multifab data to qc, qv, and qi + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qv], 0, 0, 1, mic_fab_vars[MicVar::qv]->nGrowVect()); // vapor + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qcl], 0, 1, 1, mic_fab_vars[MicVar::qcl]->nGrowVect()); // cloud water + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qci], 0, 2, 1, mic_fab_vars[MicVar::qci]->nGrowVect()); // cloud ice + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qpl], 0, 3, 1, mic_fab_vars[MicVar::qpl]->nGrowVect()); // rain + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qpi], 0, 4, 1, mic_fab_vars[MicVar::qpi]->nGrowVect()); // snow + + // Don't need to copy this since it is filled below + // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qpi], 0, 5, 1, mic_fab_vars[MicVar::qci]->nGrowVect()); // graupel + + amrex::MultiFab qgraup_mf(qmoist, amrex::make_alias, 5, 1); + + // Get the temperature, density, theta, qt and qp from input + for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar::theta]->array(mfi); + auto qt_arr = mic_fab_vars[MicVar::qt]->array(mfi); + auto qp_arr = mic_fab_vars[MicVar::qp]->array(mfi); + auto qpl_arr = mic_fab_vars[MicVar::qpl]->array(mfi); + auto qpi_arr = mic_fab_vars[MicVar::qpi]->array(mfi); + + auto qgraup_arr= qgraup_mf.array(mfi); + + const auto& box3d = mfi.tilebox(); + + // get potential total density, temperature, qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQt_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); + states_arr(i,j,k,RhoQp_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + + // Graupel == precip total - rain - snow (but must be >= 0) + qgraup_arr(i,j,k) = std::max(0.0, qp_arr(i,j,k)-qpl_arr(i,j,k)-qpi_arr(i,j,k)); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); + qmoist.FillBoundary(m_geom.periodicity()); +} + + diff --git a/Source/Microphysics/Update.cpp b/Source/Microphysics/SAM/Update.cpp~ similarity index 100% rename from Source/Microphysics/Update.cpp rename to Source/Microphysics/SAM/Update.cpp~ diff --git a/Source/TimeIntegration/ERF_advance_microphysics.cpp b/Source/TimeIntegration/ERF_advance_microphysics.cpp index 090a01935..9c0e330c2 100644 --- a/Source/TimeIntegration/ERF_advance_microphysics.cpp +++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp @@ -11,13 +11,7 @@ void ERF::advance_microphysics (int lev, grids[lev], Geom(lev), dt_advance); - - micro.Cloud(); - micro.Diagnose(); - micro.IceFall(); - micro.Precip(); - micro.MicroPrecipFall(); - + micro.Advance(); micro.Update(cons, qmoist[lev]); } #endif