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 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/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); }); }