diff --git a/Source/Microphysics/Null/Make.package~ b/Source/Microphysics/Null/Make.package~ deleted file mode 100644 index 36572d9e5..000000000 --- a/Source/Microphysics/Null/Make.package~ +++ /dev/null @@ -1,2 +0,0 @@ -CEXE_headers += Microphysics.H - diff --git a/Source/Microphysics/Null/NullMoist.H~ b/Source/Microphysics/Null/NullMoist.H~ deleted file mode 100644 index 4d1a3109c..000000000 --- a/Source/Microphysics/Null/NullMoist.H~ +++ /dev/null @@ -1,25 +0,0 @@ -#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~ deleted file mode 100644 index 0e2a37858..000000000 --- a/Source/Microphysics/SAM/Cloud.cpp~ +++ /dev/null @@ -1,134 +0,0 @@ - -#include "Microphysics.H" -#include "IndexDefines.H" -#include "TileNoZ.H" - -using namespace amrex; - -/** - * Compute Cloud-related Microphysics quantities. - */ -void Microphysics::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/SAM/Diagnose.cpp~ b/Source/Microphysics/SAM/Diagnose.cpp~ deleted file mode 100644 index 53f095530..000000000 --- a/Source/Microphysics/SAM/Diagnose.cpp~ +++ /dev/null @@ -1,59 +0,0 @@ -#include "Microphysics.H" - -/** - * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid - * from the existing Microphysics variables. - */ -void Microphysics::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 Microphysics::Proc () { - - MicroPrecipFall(); - - if (docloud) { - Cloud(); - if (doprecip) { - Precip(); - } - } - - Diagnose(); -} diff --git a/Source/Microphysics/SAM/IceFall.cpp~ b/Source/Microphysics/SAM/IceFall.cpp~ deleted file mode 100644 index 09e98076b..000000000 --- a/Source/Microphysics/SAM/IceFall.cpp~ +++ /dev/null @@ -1,171 +0,0 @@ - -#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 Microphysics::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/SAM/Init.cpp~ b/Source/Microphysics/SAM/Init.cpp~ deleted file mode 100644 index a349c59cf..000000000 --- a/Source/Microphysics/SAM/Init.cpp~ +++ /dev/null @@ -1,245 +0,0 @@ - -#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 Microphysics::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/SAM/Precip.cpp~ b/Source/Microphysics/SAM/Precip.cpp~ deleted file mode 100644 index 61f8daaff..000000000 --- a/Source/Microphysics/SAM/Precip.cpp~ +++ /dev/null @@ -1,148 +0,0 @@ -/* - * this file is modified from precip_proc from samxx - */ -#include "Microphysics.H" - -using namespace amrex; - -/** - * Compute Precipitation-related Microphysics quantities. - */ -void Microphysics::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/SAM/PrecipFall.cpp~ b/Source/Microphysics/SAM/PrecipFall.cpp~ deleted file mode 100644 index 4c896b5a9..000000000 --- a/Source/Microphysics/SAM/PrecipFall.cpp~ +++ /dev/null @@ -1,307 +0,0 @@ -#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 Microphysics::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 Microphysics::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/SAM/SAM.H~ b/Source/Microphysics/SAM/SAM.H~ deleted file mode 100644 index 09e538320..000000000 --- a/Source/Microphysics/SAM/SAM.H~ +++ /dev/null @@ -1,163 +0,0 @@ -/* - * 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~ deleted file mode 100644 index a13044380..000000000 --- a/Source/Microphysics/SAM/Update.cpp~ +++ /dev/null @@ -1,60 +0,0 @@ - -#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 Microphysics::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()); -} - -