From f83644c7e430e84b07190a6d927da75eafc3be32 Mon Sep 17 00:00:00 2001
From: Mahesh Natarajan <mn556@cornell.edu>
Date: Mon, 11 Dec 2023 11:36:23 -0800
Subject: [PATCH] Fast eddy microphysics clean up (#1336)
* Changes to EOS.H with moisture
* Change all functions in EOS.H to work with moisture
* Adding qp for buoyancy type 1
* Updating docs for buoyancy term
* Updating docs for buoyancy term
* Update docs for buoyancy
* Avoiding extra divide
* Update buoyancy docs
* Some corrections to buoyancy docs
* Some corrections to Docs in buoyancy
* Correcting docs issue
* Correcting docs
* adding Kessler microphysics docs
* adding Kessler microphysics docs
* adding Kessler microphysics docs
* adding Kessler microphysics docs
* Updating inputs
* Adding inputs for Gabersek test case:
* Correcting VisIt plot issue
* Reverting changes
* Reverting changes
* FastEddy microphysics clean up
* FastEddy microphysics clean up
---------
Co-authored-by: Mahesh Natarajan <mahesh@dogora.dhcp.lbl.gov>
Co-authored-by: Ann Almgren <asalmgren@lbl.gov>
---
Source/Microphysics/FastEddy/Diagnose_FE.cpp | 33 ---
Source/Microphysics/FastEddy/FastEddy.H | 15 +-
Source/Microphysics/FastEddy/FastEddy.cpp | 201 +++----------------
Source/Microphysics/FastEddy/Init_FE.cpp | 183 +----------------
Source/Microphysics/FastEddy/Update_FE.cpp | 26 +--
5 files changed, 40 insertions(+), 418 deletions(-)
diff --git a/Source/Microphysics/FastEddy/Diagnose_FE.cpp b/Source/Microphysics/FastEddy/Diagnose_FE.cpp
index b3c1307c1..369a5c5d5 100644
--- a/Source/Microphysics/FastEddy/Diagnose_FE.cpp
+++ b/Source/Microphysics/FastEddy/Diagnose_FE.cpp
@@ -5,37 +5,4 @@
* from the existing Microphysics variables.
*/
void FastEddy::Diagnose () {
-
- auto qt = mic_fab_vars[MicVar_FE::qt];
- auto qp = mic_fab_vars[MicVar_FE::qp];
- auto qv = mic_fab_vars[MicVar_FE::qv];
- auto qn = mic_fab_vars[MicVar_FE::qn];
- auto qcl = mic_fab_vars[MicVar_FE::qcl];
- auto qci = mic_fab_vars[MicVar_FE::qci];
- auto qpl = mic_fab_vars[MicVar_FE::qpl];
- auto qpi = mic_fab_vars[MicVar_FE::qpi];
- auto tabs = mic_fab_vars[MicVar_FE::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);
-
- const auto& box3d = mfi.tilebox();
-
- amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
- qv_array(i,j,k) = std::max(0.0,qt_array(i,j,k) - qn_array(i,j,k));
- amrex::Real omn = 1.0;
- 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 = 1.0;;
- qpl_array(i,j,k) = qp_array(i,j,k)*omp;
- qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp);
- });
- }
}
diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H
index df54a1736..a60eb888d 100644
--- a/Source/Microphysics/FastEddy/FastEddy.H
+++ b/Source/Microphysics/FastEddy/FastEddy.H
@@ -20,22 +20,11 @@
namespace MicVar_FE {
enum {
// independent variables
- qt = 0,
- qp,
+ qv = 0,
+ qc,
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
};
}
diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp
index e3575ed99..40f04132d 100644
--- a/Source/Microphysics/FastEddy/FastEddy.cpp
+++ b/Source/Microphysics/FastEddy/FastEddy.cpp
@@ -13,201 +13,58 @@ using namespace amrex;
*/
void FastEddy::AdvanceFE() {
- auto qt = mic_fab_vars[MicVar_FE::qt];
- auto qp = mic_fab_vars[MicVar_FE::qp];
- auto qn = mic_fab_vars[MicVar_FE::qn];
- auto tabs = mic_fab_vars[MicVar_FE::tabs];
+ auto tabs = mic_fab_vars[MicVar_FE::tabs];
- auto qcl = mic_fab_vars[MicVar_FE::qcl];
- auto theta = mic_fab_vars[MicVar_FE::theta];
- auto qv = mic_fab_vars[MicVar_FE::qv];
- auto rho = mic_fab_vars[MicVar_FE::rho];
-
- auto dz = m_geom.CellSize(2);
- auto domain = m_geom.Domain();
- int k_lo = domain.smallEnd(2);
- int k_hi = domain.bigEnd(2);
-
- MultiFab fz;
- auto ba = tabs->boxArray();
- auto dm = tabs->DistributionMap();
- fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells
-
- for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
- auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
- auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
- auto fz_array = fz.array(mfi);
- const Box& tbz = mfi.tilebox();
-
- ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-
- Real rho_avg, qp_avg;
-
- if (k==k_lo) {
- rho_avg = rho_array(i,j,k);
- qp_avg = qp_array(i,j,k);
- } else if (k==k_hi+1) {
- rho_avg = rho_array(i,j,k-1);
- qp_avg = qp_array(i,j,k-1);
- } else {
- rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
- qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
- }
-
- qp_avg = std::max(0.0, qp_avg);
-
- Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s
-
- fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
-
- /*if(k==0){
- fz_array(i,j,k) = 0;
- }*/
- });
- }
-
- Real dtn = dt;
-
- /*for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-
- auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi);
- auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi);
- auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
- auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi);
-
- const auto& box3d = mfi.tilebox();
-
- ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-
- qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
- qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
- qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));
-
- if(qt_array(i,j,k) == 0.0){
- qv_array(i,j,k) = 0.0;
- qn_array(i,j,k) = 0.0;
- }
- });
- }
-
-
- for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
- auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
- auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
-
- const auto& box3d = mfi.tilebox();
- auto fz_array = fz.array(mfi);
-
- ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
- Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn;
- qp_array(i,j,k) = qp_array(i,j,k) + dq_sed;
- });
- }*/
-
-
-
-
- // get the temperature, dentisy, theta, qt and qp from input
+ // get the temperature, dentisy, theta, qt and qc from input
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
- auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi);
- auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi);
- auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi);
- auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
-
- auto theta_array = theta->array(mfi);
- auto qv_array = qv->array(mfi);
- auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
+ auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi);
+ auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi);
+ auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi);
+ auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi);
+ auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
const auto& box3d = mfi.tilebox();
- auto fz_array = fz.array(mfi);
-
// Expose for GPU
Real d_fac_cond = m_fac_cond;
ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
- qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
- qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
- qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));
-
- if(qt_array(i,j,k) == 0.0){
- qv_array(i,j,k) = 0.0;
- qn_array(i,j,k) = 0.0;
- }
+ qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
//------- Autoconversion/accretion
- Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat;
+ Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat;
Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0;
erf_qsatw(tabs_array(i,j,k), pressure, qsat);
- // If there is precipitating water (i.e. rain), and the cell is not saturated
- // then the rain water can evaporate leading to extraction of latent heat, hence
- // reducing temperature and creating negative buoyancy
-
- dq_clwater_to_rain = 0.0;
- dq_rain_to_vapor = 0.0;
- dq_vapor_to_clwater = 0.0;
- dq_clwater_to_vapor = 0.0;
-
+ // If there is precipitating water (i.e. rain), and the cell is not saturated
+ // then the rain water can evaporate leading to extraction of latent heat, hence
+ // reducing temperature and creating negative buoyancy
- //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
- Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
+ dq_vapor_to_clwater = 0.0;
+ dq_clwater_to_vapor = 0.0;
- // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
- if(qv_array(i,j,k) > qsat){
- dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
- }
+ //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
+ Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
-
- // If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and
- // reducing temperature
- if(qv_array(i,j,k) < qsat and qn_array(i,j,k) > 0.0){
- dq_clwater_to_vapor = std::min(qn_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
- }
-
- if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) {
- Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046);
- dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/
- (5.4e5 + 2.55e6/(pressure*qsat))*dtn;
- // The negative sign is to make this variable (vapor formed from evaporation)
- // a poistive quantity (as qv/qs < 1)
- dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor});
-
- // Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature
- }
-
- // If there is cloud water present then do accretion and autoconversion to rain
-
- if (qn_array(i,j,k) > 0.0) {
- qcc = qn_array(i,j,k);
-
- autor = 0.0;
- if (qcc > qcw0) {
- autor = 0.001;
- }
-
- accrr = 0.0;
- accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875);
- dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001));
-
- // If the amount of change is more than the amount of qc present, then dq = qc
- dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k));
+ // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
+ if(qv_array(i,j,k) > qsat){
+ dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
+ }
+ // If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and
+ // reducing temperature
+ if(qv_array(i,j,k) < qsat and qc_array(i,j,k) > 0.0){
+ dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
}
+ qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor;
+ qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor;
- Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn;
-
- qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain;
- qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor;
- qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;
-
- theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);
+ theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor);
- qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
- qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
- qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));
+ qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
+ qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
});
}
diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp
index 033d56136..7708a1243 100644
--- a/Source/Microphysics/FastEddy/Init_FE.cpp
+++ b/Source/Microphysics/FastEddy/Init_FE.cpp
@@ -27,9 +27,6 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
m_geom = geom;
m_gtoe = grids;
- auto dz = m_geom.CellSize(2);
- auto lowz = m_geom.ProbLo(2);
-
dt = dt_advance;
// initialize microphysics variables
@@ -41,7 +38,7 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
// 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.);
+ qmoist.setVal(0.);
for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
@@ -53,83 +50,17 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
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 qv_array_from_moist = qv.array(mfi);
- auto qc_array = qc.array(mfi);
- auto qi_array = qi.array(mfi);
- auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi);
- auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
- auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi);
auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi);
+ auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi);
auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi);
auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi);
- auto pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi);
const auto& box3d = mfi.tilebox();
@@ -137,113 +68,9 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
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,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
- qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
- qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
- qv_array(i,j,k) = qv_array_from_moist(i,j,k);
+ qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
+ qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k));
- 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<Real> 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<Real> 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/FastEddy/Update_FE.cpp b/Source/Microphysics/FastEddy/Update_FE.cpp
index a157279f7..07d9a86c3 100644
--- a/Source/Microphysics/FastEddy/Update_FE.cpp
+++ b/Source/Microphysics/FastEddy/Update_FE.cpp
@@ -13,17 +13,6 @@
void FastEddy::Update (amrex::MultiFab& cons,
amrex::MultiFab& qmoist)
{
- // copy multifab data to qc, qv, and qi
- amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qv], 0, 0, 1, mic_fab_vars[MicVar_FE::qv]->nGrowVect()); // vapor
- amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qcl], 0, 1, 1, mic_fab_vars[MicVar_FE::qcl]->nGrowVect()); // cloud water
- amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qci], 0, 2, 1, mic_fab_vars[MicVar_FE::qci]->nGrowVect()); // cloud ice
- amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpl], 0, 3, 1, mic_fab_vars[MicVar_FE::qpl]->nGrowVect()); // rain
- amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpi], 0, 4, 1, mic_fab_vars[MicVar_FE::qpi]->nGrowVect()); // snow
-
- // Don't need to copy this since it is filled below
- // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpi], 0, 5, 1, mic_fab_vars[MicVar_FE::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) {
@@ -31,22 +20,15 @@ void FastEddy::Update (amrex::MultiFab& cons,
auto rho_arr = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi);
- auto qt_arr = mic_fab_vars[MicVar_FE::qt]->array(mfi);
- auto qp_arr = mic_fab_vars[MicVar_FE::qp]->array(mfi);
-
- auto qgraup_arr= qgraup_mf.array(mfi);
-
+ auto qv_arr = mic_fab_vars[MicVar_FE::qv]->array(mfi);
+ auto qc_arr = mic_fab_vars[MicVar_FE::qc]->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,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k);
- states_arr(i,j,k,RhoQ2_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) = 0.0;//
+ states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k);
+ states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k);
});
}