From 7c03e4b5ddc4a79e135434cd8f944435d4330feb Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 10 Oct 2023 12:58:18 -0700 Subject: [PATCH 01/23] Changes to EOS.H with moisture --- Source/EOS.H | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/Source/EOS.H b/Source/EOS.H index ad53ddb30..189303e7e 100644 --- a/Source/EOS.H +++ b/Source/EOS.H @@ -13,10 +13,11 @@ * @params[in] rhotheta density times potential temperature theta */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta) +amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0) { - amrex::Real p_loc = p_0 * std::pow(R_d * rhotheta * ip_0, Gamma); - return p_loc / (R_d * rho); + // rho and rhotheta are dry values. We should be using moist values with moisture + amrex::Real p_loc = p_0 * std::pow(R_d * rhotheta * (1.0 + R_v/R_d*qv) * ip_0, Gamma); + return p_loc / (R_d * rho * (1.0 + R_v/R_d*qv) ); } /** @@ -27,9 +28,9 @@ amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta) * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp) +amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp, const amrex::Real qv=0.0) { - amrex::Real p_loc = rho * R_d * T; + amrex::Real p_loc = rho * R_d * T * (1.0 + R_v/R_d*qv); return T * std::pow((p_0/p_loc),rdOcp); } @@ -65,9 +66,9 @@ amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv = 0.) * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp) +amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0) { - return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * theta); + return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * theta * (1.0 + R_v/R_d*qv) ); } /** From 4a20716be52d3676381e54a635c71384f1bf7a37 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 10 Oct 2023 13:29:23 -0700 Subject: [PATCH 02/23] Change all functions in EOS.H to work with moisture --- Source/EOS.H | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/Source/EOS.H b/Source/EOS.H index 189303e7e..886674bba 100644 --- a/Source/EOS.H +++ b/Source/EOS.H @@ -15,8 +15,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0) { - // rho and rhotheta are dry values. We should be using moist values with moisture + // rho and rhotheta are dry values. We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) amrex::Real p_loc = p_0 * std::pow(R_d * rhotheta * (1.0 + R_v/R_d*qv) * ip_0, Gamma); + // p = rho_d * R_d * T_v (not T) + // T_v = T * (1 + R_v/R_d*qv) return p_loc / (R_d * rho * (1.0 + R_v/R_d*qv) ); } @@ -30,7 +33,9 @@ amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp, const amrex::Real qv=0.0) { + // p = rho_d * R_d * T_moist amrex::Real p_loc = rho * R_d * T * (1.0 + R_v/R_d*qv); + // theta_d = T * (p0/p)^(R_d/C_p) return T * std::pow((p_0/p_loc),rdOcp); } @@ -68,6 +73,8 @@ amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv = 0.) AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0) { + // We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * theta * (1.0 + R_v/R_d*qv) ); } @@ -79,9 +86,11 @@ amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getdPdRgivenConstantTheta(const amrex::Real rho, const amrex::Real theta) +amrex::Real getdPdRgivenConstantTheta(const amrex::Real rho, const amrex::Real theta, const amrex::Real qv=0.0) { - return Gamma * p_0 * std::pow( (R_d * theta * ip_0), Gamma) * std::pow(rho, Gamma-1.0) ; + // We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) + return Gamma * p_0 * std::pow( (R_d * theta * (1.0 + R_v/R_d*qv) * ip_0), Gamma) * std::pow(rho, Gamma-1.0) ; } /** @@ -103,10 +112,12 @@ amrex::Real getExnergivenP(const amrex::Real P, const amrex::Real rdOcp) * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp) +amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=0.0 ) { // Exner function pi in terms of (rho theta) - return std::pow(R_d * rhotheta * ip_0, Gamma * rdOcp); + // We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) + return std::pow(R_d * rhotheta * (1.0 + R_v/R_d*qv) * ip_0, Gamma * rdOcp); } /** @@ -115,11 +126,12 @@ amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp * @params[in] p pressure */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getRhoThetagivenP(const amrex::Real p) +amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0) { // diagnostic relation for the full pressure // see https://erf.readthedocs.io/en/latest/theory/NavierStokesEquations.html - return std::pow(p*std::pow(p_0, Gamma-1), iGamma) * iR_d; + // For cases with mositure theta = theta_m / (1 + R_v/R_d*qv) + return std::pow(p*std::pow(p_0, Gamma-1), iGamma) * iR_d / (1.0 + R_v/R_d*qv) ; } #endif From 5e6005609d1c19d87255a107a710ac91ab437a0d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:35:29 -0800 Subject: [PATCH 03/23] Adding qp for buoyancy type 1 --- Source/TimeIntegration/ERF_make_buoyancy.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index 7ad36127c..b3c461a14 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -173,9 +173,9 @@ void make_buoyancy (Vector& S_data, amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k ) - + qi_data(i,j,k )) - r0_arr(i,j,k ); + + qi_data(i,j,k ) + cell_data(i,j,k,RhoQp_comp)/cell_data(i,j,k,Rho_comp)) - r0_arr(i,j,k ); Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1) - + qi_data(i,j,k-1)) - r0_arr(i,j,k-1); + + qi_data(i,j,k-1) + + cell_data(i,j,k-1,RhoQp_comp)/cell_data(i,j,k-1,Rho_comp)) - r0_arr(i,j,k-1); buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); }); } // mfi From 5582efc27d497b6f62dcc42d2627c45c87b15b81 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:48:46 -0800 Subject: [PATCH 04/23] Updating docs for buoyancy term --- Docs/sphinx_doc/theory/Buoyancy.rst | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 9e5f0d819..36cc012f8 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,6 +7,28 @@ .. _Buoyancy: +Density of the mixture +======================== + +The total density in a given cell is given by +.. math:: + \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + +where `m_a` is the mass of dry air, `m_v` is the mass of water vapor, `m_c` is the mass of liquid water, and `m_p` is the mass of precipitate. +From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component + +.. math:: + q_i \equiv = \frac{m_i}{m_a}. + +Using this we can write + +.. math:: + \rho = m_a\frac{(1 + q_v + q_c + q_p)}{V} + = \rho_d(1 + q_v + q_c + q_p), + +where `\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. + + Buoyancy ========= From ee798397b0959701d99f47fbe3f41fc5818ca2b6 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:51:04 -0800 Subject: [PATCH 05/23] Updating docs for buoyancy term --- Docs/sphinx_doc/theory/Buoyancy.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 36cc012f8..2005be0dc 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -11,6 +11,7 @@ Density of the mixture ======================== The total density in a given cell is given by + .. math:: \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, From 38c286126925d1cfbd4697e9a18766735d16e0bb Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:55:29 -0800 Subject: [PATCH 06/23] Update docs for buoyancy --- Docs/sphinx_doc/theory/Buoyancy.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 2005be0dc..1eb960c17 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -15,7 +15,7 @@ The total density in a given cell is given by .. math:: \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, -where `m_a` is the mass of dry air, `m_v` is the mass of water vapor, `m_c` is the mass of liquid water, and `m_p` is the mass of precipitate. +where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component .. math:: @@ -27,7 +27,7 @@ Using this we can write \rho = m_a\frac{(1 + q_v + q_c + q_p)}{V} = \rho_d(1 + q_v + q_c + q_p), -where `\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. +where :math:`\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. Buoyancy From fb48acc89f673a07be777772eeb6484fc97e217d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 16:09:37 -0800 Subject: [PATCH 07/23] Avoiding extra divide --- Source/TimeIntegration/ERF_make_buoyancy.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index b3c461a14..6e3b7f39e 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -173,9 +173,9 @@ void make_buoyancy (Vector& S_data, amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k ) - + qi_data(i,j,k ) + cell_data(i,j,k,RhoQp_comp)/cell_data(i,j,k,Rho_comp)) - r0_arr(i,j,k ); + + qi_data(i,j,k )) + cell_data(i,j,k,RhoQp_comp) - r0_arr(i,j,k ); Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1) - + qi_data(i,j,k-1) + + cell_data(i,j,k-1,RhoQp_comp)/cell_data(i,j,k-1,Rho_comp)) - r0_arr(i,j,k-1); + + qi_data(i,j,k-1)) + cell_data(i,j,k-1,RhoQp_comp) - r0_arr(i,j,k-1); buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); }); } // mfi From 926ad1ef197ebe9ee702dda3ee586f36b1e4b6b7 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 16:34:23 -0800 Subject: [PATCH 08/23] Update buoyancy docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 1eb960c17..204d262e3 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -13,7 +13,7 @@ Density of the mixture The total density in a given cell is given by .. math:: - \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component From fe0b6d67c114498e0881550959ff31fa08633972 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 10:26:09 -0800 Subject: [PATCH 09/23] Some corrections to buoyancy docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 204d262e3..87df660f0 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,19 +7,25 @@ .. _Buoyancy: +Buoyancy +========= + +ERF has three options for how to define the buoyancy force. Even in the absence of moisture these +expressions are not equivalent. + Density of the mixture -======================== +----------------------- -The total density in a given cell is given by +The total density in a cell containing air, water vapor, liquid water and precipitates is given by .. math:: - \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component .. math:: - q_i \equiv = \frac{m_i}{m_a}. + q_i \equiv \frac{m_i}{m_a}. Using this we can write @@ -30,12 +36,6 @@ Using this we can write where :math:`\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. -Buoyancy -========= - -ERF has three options for how to define the buoyancy force. Even in the absence of moisture these -expressions are not equivalent. - Type 1 ------ @@ -47,8 +47,10 @@ One version of the buoyancy force is expressed simply as .. math:: \rho^\prime = \rho_{total} - \rho_0 -where the full density :math:`\rho_{total}` is the sum of dry and moist components and :math:`\rho_0` is the base state density -for dry air only. +where the total density :math:`\rho_{total} = \rho_d(1 + q_v + q_c + q_p)` is the sum of dry and moist components and :math:`\rho_0` is the total density +for the background state. For eg., a usual scenario is that of a background state that contains only air and vapor and no cloud water or precipitates. For such a state, +the total background density :math:`\rho_0 = \rho_{d_0}(1 + q_{v_0})`, where :math:`\rho_{d_0}` and :math:`q_{v_0}` are the background dry density and vapor mixing ratio respectively. +As a check, we observe that :math:`\rho^\prime_0 = 0`, which means that the background state is not buoyant. Type 2 ------ From a89a48f0162e973d50ca4ed2f0a1b330f1404c80 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 10:33:44 -0800 Subject: [PATCH 10/23] Some corrections to Docs in buoyancy --- Docs/sphinx_doc/theory/Buoyancy.rst | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index de71f28b9..87df660f0 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,29 +7,6 @@ .. _Buoyancy: -Density of the mixture -======================== - -The total density in a given cell is given by - -.. math:: - \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, - -where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. -From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component - -.. math:: - q_i \equiv = \frac{m_i}{m_a}. - -Using this we can write - -.. math:: - \rho = m_a\frac{(1 + q_v + q_c + q_p)}{V} - = \rho_d(1 + q_v + q_c + q_p), - -where :math:`\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. - - Buoyancy ========= From ab100b4f3fd1bc30bf048f95b0ba2431f740821b Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 12:20:04 -0800 Subject: [PATCH 11/23] Correcting docs issue --- Docs/sphinx_doc/theory/Buoyancy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 87df660f0..6213c379f 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -19,7 +19,7 @@ Density of the mixture The total density in a cell containing air, water vapor, liquid water and precipitates is given by .. math:: - \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component From d277d3390729f47103a51baa405b22d3c41874ea Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 12:22:43 -0800 Subject: [PATCH 12/23] Correcting docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 87df660f0..6213c379f 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -19,7 +19,7 @@ Density of the mixture The total density in a cell containing air, water vapor, liquid water and precipitates is given by .. math:: - \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component From 051117a4c3df86f352df0d709c43be1c0c48222d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:07:49 -0800 Subject: [PATCH 13/23] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index 42684f0cf..eadc410e8 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -7,6 +7,21 @@ .. _Microphysics: +Kessler Microphysics model +=========================== +Governing equations for the microphysical quantities for Kessler microphysics from `gabervsek2012dry`_ are + +.. math:: + \frac{\partial q_v}{\partial t} = -C_c + E_c + E_r + \frac{\partial q_c}{\partial t} = C_c - E_c - (A_c + K_c) + \frac{\partial q_p}{\partial t} = \frac{1}{\overline{\rho}}\frac{\partial}{\partial z}(\overline{\rho}Vq_p) + (A_c + K_c) - E_r + \frac{\partial q_t}{\partial t} = \frac{\partial q_v}{\partial t} + \frac{\partial q_c}{\partial t} + = E_r - (A_c + K_c) +where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in \cite{klemp1978simulation}, and is given below. Note that in all the equations, $p$ is specified in millibars and $\overline{\rho}$ is specified in g cm$^{-3}$. + +.. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml + + Single Moment Microphysics Model =================================== The conversion rates among the moist hydrometeors are parameterized assuming that From 16e76d76d1b0e9c591489b5b96a921a2b84950ab Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:10:32 -0800 Subject: [PATCH 14/23] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index eadc410e8..d18b21d7f 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -13,10 +13,14 @@ Governing equations for the microphysical quantities for Kessler microphysics fr .. math:: \frac{\partial q_v}{\partial t} = -C_c + E_c + E_r +.. math:: \frac{\partial q_c}{\partial t} = C_c - E_c - (A_c + K_c) +.. math:: \frac{\partial q_p}{\partial t} = \frac{1}{\overline{\rho}}\frac{\partial}{\partial z}(\overline{\rho}Vq_p) + (A_c + K_c) - E_r +.. math:: \frac{\partial q_t}{\partial t} = \frac{\partial q_v}{\partial t} + \frac{\partial q_c}{\partial t} = E_r - (A_c + K_c) + where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in \cite{klemp1978simulation}, and is given below. Note that in all the equations, $p$ is specified in millibars and $\overline{\rho}$ is specified in g cm$^{-3}$. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml From f33cdd8060d3bb8e5ed26dc46c554ff8dd6d3baa Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:14:35 -0800 Subject: [PATCH 15/23] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index d18b21d7f..f747ace71 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -21,9 +21,13 @@ Governing equations for the microphysical quantities for Kessler microphysics fr \frac{\partial q_t}{\partial t} = \frac{\partial q_v}{\partial t} + \frac{\partial q_c}{\partial t} = E_r - (A_c + K_c) -where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in \cite{klemp1978simulation}, and is given below. Note that in all the equations, $p$ is specified in millibars and $\overline{\rho}$ is specified in g cm$^{-3}$. +where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, +:math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of +rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in `klemp1978simulation`_, and is given +below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm:math:`^{-3}`. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml +.. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml Single Moment Microphysics Model From a0bb74dddf0c5de225692a763a2bb01bcfaf0aec Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:16:11 -0800 Subject: [PATCH 16/23] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index f747ace71..a644d86d4 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -24,7 +24,7 @@ Governing equations for the microphysical quantities for Kessler microphysics fr where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in `klemp1978simulation`_, and is given -below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm:math:`^{-3}`. +below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml .. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml From 56957bd18dd54149bf21dce22ba35eebb3710874 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 30 Nov 2023 09:30:58 -0800 Subject: [PATCH 17/23] Updating inputs --- Exec/SquallLine_2D/inputs_moisture | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 2e2757fe5..7ede7c674 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -68,9 +68,9 @@ erf.les_type = "None" erf.molec_diff_type = "ConstantAlpha" #erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities -erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.dynamicViscosity = 100.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s erf.alpha_T = 00.0 # [m^2/s] -erf.alpha_C = 25.0 +erf.alpha_C = 100.0 erf.moisture_model = "Kessler" erf.use_moist_background = true From 14073876bc2c1f404de1303a68ed1cd111d9631d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Sat, 2 Dec 2023 16:38:09 -0800 Subject: [PATCH 18/23] Adding inputs for Gabersek test case: --- Exec/SquallLine_2D/inputs_moisture_Gabersek | 91 +++++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 Exec/SquallLine_2D/inputs_moisture_Gabersek diff --git a/Exec/SquallLine_2D/inputs_moisture_Gabersek b/Exec/SquallLine_2D/inputs_moisture_Gabersek new file mode 100644 index 000000000..ad5b642b4 --- /dev/null +++ b/Exec/SquallLine_2D/inputs_moisture_Gabersek @@ -0,0 +1,91 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 18000 +stop_time = 90000.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 2048 1024 2048 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = -60000. 0. 0. +geometry.prob_hi = 60000. 400. 20000. +amr.n_cell = 384 4 128 # dx=dy=dz=100 m + +# periodic in x to match WRF setup +# - as an alternative, could use symmetry at x=0 and outflow at x=25600 +geometry.is_periodic = 1 1 0 +#xlo.type = "Outflow" +#xhi.type = "Outflow" +zlo.type = "SlipWall" +zhi.type = "Outflow" + +erf.sponge_strength = 2.0 +#erf.use_zhi_sponge_damping = true +erf.zhi_sponge_start = 12000.0 + +erf.sponge_density = 1.2 +erf.sponge_x_velocity = 0.0 +erf.sponge_y_velocity = 0.0 +erf.sponge_z_velocity = 0.0 + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 0.5 # fixed time step [s] -- Straka et al 1993 +erf.fixed_fast_dt = 0.25 # fixed time step [s] -- Straka et al 1993 +#erf.no_substepping = 1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +amr.check_file = chk # root name of checkpoint file +amr.check_int = 1000 # number of timesteps between checkpoints +#amr.restart = chk09000 + +# PLOTFILES +erf.plot_file_1 = plt # root name of plotfile +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens + +# SOLVER CHOICE +erf.use_gravity = true +erf.buoyancy_type = 1 +erf.use_coriolis = false +erf.use_rayleigh_damping = false + +#erf.les_type = "Smagorinsky" +erf.Cs = 0.25 +erf.les_type = "None" + +# +# diffusion coefficient from Straka, K = 75 m^2/s +# +erf.molec_diff_type = "ConstantAlpha" +#erf.molec_diff_type = "Constant" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 100.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 00.0 # [m^2/s] +erf.alpha_C = 50.0 + +erf.moisture_model = "Kessler" +erf.use_moist_background = true + +erf.moistscal_horiz_adv_string = "Centered_2nd" +erf.moistscal_vert_adv_string = "Centered_2nd" + +# PROBLEM PARAMETERS (optional) +prob.z_tr = 12000.0 +prob.height = 1200.0 +prob.theta_0 = 300.0 +prob.theta_tr = 343.0 +prob.T_tr = 213.0 +prob.x_c = 0.0 +prob.z_c = 2000.0 +prob.x_r = 10000.0 +prob.z_r = 1500.0 +prob.theta_c = 3.0 From 59439a708494a34b86ba99a961dd6cb2fc9981ea Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 6 Dec 2023 22:22:56 -0800 Subject: [PATCH 19/23] Correcting VisIt plot issue --- Source/Microphysics/Kessler/Kessler.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 529bb03e9..6c32a0ff2 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -197,8 +197,12 @@ void Kessler::AdvanceKessler() { } + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; + if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; - + if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; + //dq_sed = 0.0; + 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; From 398ce20be4f55489bc5ebe29f3206d01606b4e8a Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Mon, 11 Dec 2023 09:42:10 -0800 Subject: [PATCH 20/23] Reverting changes --- Exec/SquallLine_2D/inputs_moisture | 4 +- Exec/SquallLine_2D/inputs_moisture_Gabersek | 91 --------------------- 2 files changed, 2 insertions(+), 93 deletions(-) delete mode 100644 Exec/SquallLine_2D/inputs_moisture_Gabersek diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 7ede7c674..2e2757fe5 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -68,9 +68,9 @@ erf.les_type = "None" erf.molec_diff_type = "ConstantAlpha" #erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities -erf.dynamicViscosity = 100.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s erf.alpha_T = 00.0 # [m^2/s] -erf.alpha_C = 100.0 +erf.alpha_C = 25.0 erf.moisture_model = "Kessler" erf.use_moist_background = true diff --git a/Exec/SquallLine_2D/inputs_moisture_Gabersek b/Exec/SquallLine_2D/inputs_moisture_Gabersek deleted file mode 100644 index ad5b642b4..000000000 --- a/Exec/SquallLine_2D/inputs_moisture_Gabersek +++ /dev/null @@ -1,91 +0,0 @@ -# ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 18000 -stop_time = 90000.0 - -amrex.fpe_trap_invalid = 1 - -fabarray.mfiter_tile_size = 2048 1024 2048 - -# PROBLEM SIZE & GEOMETRY -geometry.prob_lo = -60000. 0. 0. -geometry.prob_hi = 60000. 400. 20000. -amr.n_cell = 384 4 128 # dx=dy=dz=100 m - -# periodic in x to match WRF setup -# - as an alternative, could use symmetry at x=0 and outflow at x=25600 -geometry.is_periodic = 1 1 0 -#xlo.type = "Outflow" -#xhi.type = "Outflow" -zlo.type = "SlipWall" -zhi.type = "Outflow" - -erf.sponge_strength = 2.0 -#erf.use_zhi_sponge_damping = true -erf.zhi_sponge_start = 12000.0 - -erf.sponge_density = 1.2 -erf.sponge_x_velocity = 0.0 -erf.sponge_y_velocity = 0.0 -erf.sponge_z_velocity = 0.0 - -# TIME STEP CONTROL -erf.use_native_mri = 1 -erf.fixed_dt = 0.5 # fixed time step [s] -- Straka et al 1993 -erf.fixed_fast_dt = 0.25 # fixed time step [s] -- Straka et al 1993 -#erf.no_substepping = 1 - -# DIAGNOSTICS & VERBOSITY -erf.sum_interval = 1 # timesteps between computing mass -erf.v = 1 # verbosity in ERF.cpp -amr.v = 1 # verbosity in Amr.cpp - -# REFINEMENT / REGRIDDING -amr.max_level = 0 # maximum level number allowed - -# CHECKPOINT FILES -amr.check_file = chk # root name of checkpoint file -amr.check_int = 1000 # number of timesteps between checkpoints -#amr.restart = chk09000 - -# PLOTFILES -erf.plot_file_1 = plt # root name of plotfile -erf.plot_int_1 = 100 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens - -# SOLVER CHOICE -erf.use_gravity = true -erf.buoyancy_type = 1 -erf.use_coriolis = false -erf.use_rayleigh_damping = false - -#erf.les_type = "Smagorinsky" -erf.Cs = 0.25 -erf.les_type = "None" - -# -# diffusion coefficient from Straka, K = 75 m^2/s -# -erf.molec_diff_type = "ConstantAlpha" -#erf.molec_diff_type = "Constant" -erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities -erf.dynamicViscosity = 100.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s -erf.alpha_T = 00.0 # [m^2/s] -erf.alpha_C = 50.0 - -erf.moisture_model = "Kessler" -erf.use_moist_background = true - -erf.moistscal_horiz_adv_string = "Centered_2nd" -erf.moistscal_vert_adv_string = "Centered_2nd" - -# PROBLEM PARAMETERS (optional) -prob.z_tr = 12000.0 -prob.height = 1200.0 -prob.theta_0 = 300.0 -prob.theta_tr = 343.0 -prob.T_tr = 213.0 -prob.x_c = 0.0 -prob.z_c = 2000.0 -prob.x_r = 10000.0 -prob.z_r = 1500.0 -prob.theta_c = 3.0 From bb7f3a8e535269acfc00142600a49bf3f36f9ae9 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Mon, 11 Dec 2023 09:45:34 -0800 Subject: [PATCH 21/23] Reverting changes --- Source/Microphysics/Kessler/Kessler.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index e8e13e50e..64fe35639 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -196,6 +196,7 @@ void Kessler::AdvanceKessler() { dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k)); } + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; From 56dac2853de933ebf9110a7aaf19f1722188093b Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Mon, 11 Dec 2023 10:45:12 -0800 Subject: [PATCH 22/23] FastEddy microphysics clean up --- 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..e86b9b5d5 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); }); } From 91257bc7dc1794a69b8236f1e26676696d2b8c7f Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Mon, 11 Dec 2023 10:47:40 -0800 Subject: [PATCH 23/23] FastEddy microphysics clean up --- Source/Microphysics/FastEddy/FastEddy.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp index e86b9b5d5..40f04132d 100644 --- a/Source/Microphysics/FastEddy/FastEddy.cpp +++ b/Source/Microphysics/FastEddy/FastEddy.cpp @@ -50,12 +50,12 @@ void FastEddy::AdvanceFE() { // 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)); + 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)); + 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;