add advection, sedimentation and diffusion for 2M; next add a test

Author: sajjadazimi

src/cache/precipitation_precomputed_quantities.jl

@@ -104,7 +104,7 @@ function set_precipitation_velocities!(
     moisture_model::NonEquilMoistModel,
     precip_model::Microphysics2Moment,
 )
-    (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwNₗ, ᶜwNᵢ, ᶜwNᵣ, ᶜwNₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
+    (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwnₗ, ᶜwnᵢ, ᶜwnᵣ, ᶜwnₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
     (; q_liq, q_ice, q_rai, q_sno) = p.precomputed.ᶜspecific
     (; ᶜΦ) = p.core
 
@@ -115,21 +115,21 @@ function set_precipitation_velocities!(
 
     # compute the precipitation terminal velocity [m/s]
     # TODO sedimentation of snow is based on the 1M scheme
-    @. ᶜwNᵣ = getindex(CM2.rain_terminal_velocity(
+    @. ᶜwnᵣ = getindex(CM2.rain_terminal_velocity(
         cm2p.sb,
         cm2p.tv,
         q_rai,
         Y.c.ρ,
-        Y.c.N_rai,
+        Y.c.ρn_rai,
     ), 1)
     @. ᶜwᵣ = getindex(CM2.rain_terminal_velocity(
         cm2p.sb,
         cm2p.tv,
         q_rai,
         Y.c.ρ,
-        Y.c.N_rai,
+        Y.c.ρn_rai,
     ), 2)
-    @. ᶜwNₛ = CM1.terminal_velocity(
+    @. ᶜwnₛ = CM1.terminal_velocity(
         cm1p.ps,
         cm1p.tv.snow,
         Y.c.ρ,
@@ -144,7 +144,7 @@ function set_precipitation_velocities!(
     # compute sedimentation velocity for cloud condensate [m/s]
     # TODO sedimentation velocities of cloud condensates are based 
     # on the 1M scheme.
-    @. ᶜwNₗ = CMNe.terminal_velocity(
+    @. ᶜwnₗ = CMNe.terminal_velocity(
         cm1c.liquid,
         cm1c.Ch2022.rain,
         Y.c.ρ,
@@ -156,7 +156,7 @@ function set_precipitation_velocities!(
         Y.c.ρ,
         q_liq,
     )
-    @. ᶜwNᵢ = CMNe.terminal_velocity(
+    @. ᶜwnᵢ = CMNe.terminal_velocity(
         cm1c.ice,
         cm1c.Ch2022.small_ice,
         Y.c.ρ,
@@ -355,35 +355,34 @@ function set_precipitation_cache!(Y, p, ::Microphysics2Moment, _)
     (; dt) = p
     (; ᶜts) = p.precomputed
     (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed
-    (; ᶜSNₗᵖ, ᶜSNᵢᵖ, ᶜSNᵣᵖ, ᶜSNₛᵖ) = p.precomputed
+    (; ᶜSnₗᵖ, ᶜSnᵢᵖ, ᶜSnᵣᵖ, ᶜSnₛᵖ) = p.precomputed
 
-    (; q_liq, q_rai, q_sno) = p.precomputed.ᶜspecific
+    (; q_liq, q_rai, n_liq, n_rai) = p.precomputed.ᶜspecific
 
     ᶜSᵖ = p.scratch.ᶜtemp_scalar
     ᶜS₂ᵖ = p.scratch.ᶜtemp_scalar_2
 
     # get thermodynamics and microphysics params
     (; params) = p
-    cm1p = CAP.microphysics_1m_params(params)
-    cm2p = CAP.microphysics_2m_params(params)
+    cmp = CAP.microphysics_2m_params(params)
     thp = CAP.thermodynamics_params(params)
 
     # compute warm precipitation sources on the grid mean (based on SB2006 2M scheme)
     compute_warm_precipitation_sources_2M!(
         ᶜSᵖ,
         ᶜS₂ᵖ,
-        ᶜSNₗᵖ,
-        ᶜSNᵣᵖ,
+        ᶜSnₗᵖ,
+        ᶜSnᵣᵖ,
         ᶜSqₗᵖ,
         ᶜSqᵣᵖ,
         Y.c.ρ,
-        Y.c.N_liq,
-        Y.c.N_rai,
+        n_liq,
+        n_rai,
         q_liq,
         q_rai,
         ᶜts,
         dt,
-        cm2p,
+        cmp,
         thp,
     )
 

src/cache/precomputed_quantities.jl

@@ -126,35 +126,35 @@ function precomputed_quantities(Y, atmos)
     )
     sedimentation_quantities =
         atmos.moisture_model isa NonEquilMoistModel ?
-        (; ᶜwₗ = similar(Y.c, FT), ᶜwᵢ = similar(Y.c, FT), ᶜwNₗ = similar(Y.c, FT), ᶜwNᵢ = similar(Y.c, FT)) : (;)
+        (; ᶜwₗ = similar(Y.c, FT), ᶜwᵢ = similar(Y.c, FT), ᶜwnₗ = similar(Y.c, FT), ᶜwnᵢ = similar(Y.c, FT)) : (;)
     if atmos.precip_model isa Microphysics0Moment
         precipitation_quantities =
             (; ᶜS_ρq_tot = similar(Y.c, FT), ᶜS_ρe_tot = similar(Y.c, FT))
-        elseif atmos.precip_model isa Microphysics1Moment
-            precipitation_quantities = (;
-                ᶜwᵣ = similar(Y.c, FT),
-                ᶜwₛ = similar(Y.c, FT),
-                ᶜSqₗᵖ = similar(Y.c, FT),
-                ᶜSqᵢᵖ = similar(Y.c, FT),
-                ᶜSqᵣᵖ = similar(Y.c, FT),
-                ᶜSqₛᵖ = similar(Y.c, FT),
-            )
-        elseif atmos.precip_model isa Microphysics2Moment
-            precipitation_quantities = (;
-                ᶜwᵣ = similar(Y.c, FT),
-                ᶜwₛ = similar(Y.c, FT),
-                ᶜSqₗᵖ = similar(Y.c, FT),
-                ᶜSqᵢᵖ = similar(Y.c, FT),
-                ᶜSqᵣᵖ = similar(Y.c, FT),
-                ᶜSqₛᵖ = similar(Y.c, FT),
-                ᶜwNᵣ = similar(Y.c, FT),
-                ᶜwNₛ = similar(Y.c, FT),
-                ᶜSNₗᵖ = similar(Y.c, FT),
-                ᶜSNᵢᵖ = similar(Y.c, FT),
-                ᶜSNᵣᵖ = similar(Y.c, FT),
-                ᶜSNₛᵖ = similar(Y.c, FT),
-            )
-        else
+    elseif atmos.precip_model isa Microphysics1Moment
+        precipitation_quantities = (;
+            ᶜwᵣ = similar(Y.c, FT),
+            ᶜwₛ = similar(Y.c, FT),
+            ᶜSqₗᵖ = similar(Y.c, FT),
+            ᶜSqᵢᵖ = similar(Y.c, FT),
+            ᶜSqᵣᵖ = similar(Y.c, FT),
+            ᶜSqₛᵖ = similar(Y.c, FT),
+        )
+    elseif atmos.precip_model isa Microphysics2Moment
+        precipitation_quantities = (;
+            ᶜwᵣ = similar(Y.c, FT),
+            ᶜwₛ = similar(Y.c, FT),
+            ᶜSqₗᵖ = similar(Y.c, FT),
+            ᶜSqᵢᵖ = similar(Y.c, FT),
+            ᶜSqᵣᵖ = similar(Y.c, FT),
+            ᶜSqₛᵖ = similar(Y.c, FT),
+            ᶜwnᵣ = similar(Y.c, FT),
+            ᶜwnₛ = similar(Y.c, FT),
+            ᶜSnₗᵖ = similar(Y.c, FT),
+            ᶜSnᵢᵖ = similar(Y.c, FT),
+            ᶜSnᵣᵖ = similar(Y.c, FT),
+            ᶜSnₛᵖ = similar(Y.c, FT),
+        )
+    else
         precipitation_quantities = (;)
     end
     precipitation_sgs_quantities =

src/initial_conditions/atmos_state.jl

@@ -130,10 +130,10 @@ precip_variables(ls, ::Microphysics1Moment) = (;
 precip_variables(ls, ::Microphysics2Moment) = (;
     ρq_rai = ls.ρ * ls.precip_state.q_rai,
     ρq_sno = ls.ρ * ls.precip_state.q_sno,
-    N_rai = ls.precip_state.N_rai,
-    N_sno = ls.precip_state.N_sno,
-    N_liq = ls.precip_state.N_liq,
-    N_ice = ls.precip_state.N_ice,
+    ρn_rai = ls.ρ * ls.precip_state.n_rai,
+    ρn_sno = ls.ρ * ls.precip_state.n_sno,
+    ρn_liq = ls.ρ * ls.precip_state.n_liq,
+    ρn_ice = ls.ρ * ls.precip_state.n_ice,
 )
 
 # We can use paper-based cases for LES type configurations (no TKE)

src/parameterized_tendencies/microphysics/cloud_condensate.jl

@@ -47,14 +47,14 @@ function cloud_condensate_tendency!(
 )
     (; ᶜts) = p.precomputed
     (; params, dt) = p
-    (; q_rai, q_sno) = p.precomputed.ᶜspecific
+    (; q_rai, q_sno, n_liq, n_ice) = p.precomputed.ᶜspecific
     FT = eltype(params)
     thp = CAP.thermodynamics_params(params)
     cmc = CAP.microphysics_cloud_params(params)
 
     @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources(cmc.liquid, thp, ᶜts, q_rai, dt)
     @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources(cmc.ice, thp, ᶜts, q_sno, dt)
 
-    @. Yₜ.c.N_liq += aerosol_activation_sources(cmc.liquid, thp, ᶜts, Y.c.ρ, q_rai, Y.c.N_liq, cmc.N_cloud_liquid_droplets, dt)
-    @. Yₜ.c.N_ice += aerosol_activation_sources(cmc.ice, thp, ᶜts, Y.c.ρ, q_sno, Y.c.N_ice, cmc.N_cloud_liquid_droplets, dt)
+    @. Yₜ.c.ρn_liq += Y.c.ρ * aerosol_activation_sources(cmc.liquid, thp, ᶜts, q_rai, n_liq, cmc.N_cloud_liquid_droplets, dt)
+    @. Yₜ.c.ρn_ice += Y.c.ρ * aerosol_activation_sources(cmc.ice, thp, ᶜts, q_sno, n_ice, cmc.N_cloud_liquid_droplets, dt)
 end

src/parameterized_tendencies/microphysics/microphysics_wrappers.jl

@@ -62,7 +62,7 @@ function ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq)
 end
 
 """
-    cloud_mass_sources(cm_params, thp, ts, dt)
+    cloud_sources(cm_params, thp, ts, dt)
 
  - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate
  - thp - Thermodynamics parameters struct
@@ -336,10 +336,9 @@ end
  - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate
  - thp - Thermodynamics parameters struct
  - ts - thermodynamics state
- - ρ - air density
  - qₚ - precipitation (rain or snow) specific humidity
- - N_dp - droplets (liquid or ice) number concentration
- _ N_dp_prescribed - droplets (liquid or ice) prescribed number concentration
+ - n_dp - number concentration droplets (liquid or ice) per mass
+ _ n_dp_prescribed - prescribed number concentration of droplets (liquid or ice) per mass
  - dt - model time step
 
 Returns the activation rate. #TODO This function temporarily computes activation rate 
@@ -349,30 +348,29 @@ function aerosol_activation_sources(
     cm_params::Union{CMP.CloudLiquid{FT}, CMP.CloudIce{FT}},
     thp,
     ts,
-    ρ,
     qₚ,
-    N_dp,
-    N_dp_prescribed,
+    n_dp,
+    n_dp_prescribed,
     dt,
 ) where {FT}
     r_dp = FT(2e-6) # 2 μm
-    m_dp = 4 / 3 * π * r_dp^3
-    S = ρ * cloud_sources(cm_params, thp, ts, qₚ, dt) / m_dp
+    m_dp = 4 / 3 * π * r_dp^3 * 1000 # ρ_water = 1000 kg/m^3
+    Sn = cloud_sources(cm_params, thp, ts, qₚ, dt) / m_dp
 
     return ifelse(
-        S > FT(0),
-        triangle_inequality_limiter(S, limit((N_dp_prescribed - N_dp), dt, 2)),
-        -triangle_inequality_limiter(abs(S), limit(N_dp, dt, 2)),
+        Sn > FT(0),
+        triangle_inequality_limiter(Sn, limit((n_dp_prescribed - n_dp), dt, 2)),
+        -triangle_inequality_limiter(abs(Sn), limit(n_dp, dt, 2)),
         )
 end
 
 """
-    compute_warm_precipitation_sources_2M!(Sᵖ, S₂ᵖ, SNₗᵖ, SNᵣᵖ, Sqₗᵖ, Sqᵣᵖ, ρ, Nₗ, Nᵣ, qₗ, qᵣ, dt, sb, thp)
+    compute_warm_precipitation_sources_2M!(Sᵖ, S₂ᵖ, Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ, ρ, nₗ, nᵣ, qₗ, qᵣ, dt, sb, thp)
 
  - Sᵖ, S₂ᵖ - temporary containters to help compute precipitation source terms
- - SNₗᵖ, SNᵣᵖ, Sqₗᵖ, Sqᵣᵖ - cached storage for precipitation source terms
+ - Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ - cached storage for precipitation source terms
  - ρ - air density
- - Nₗ, Nᵣ - cloud liquid and rain number concentration
+ - nₗ, nᵣ - cloud liquid and rain number concentration per mass
  - qₗ, qᵣ - cloud liquid and rain specific humidity
  - ts - thermodynamic state (see td package for details)
  - dt - model time step
@@ -384,13 +382,13 @@ Seifert and Beheng (2006) scheme.
 function compute_warm_precipitation_sources_2M!(
     Sᵖ,
     S₂ᵖ,
-    SNₗᵖ,
-    SNᵣᵖ,
+    Snₗᵖ,
+    Snᵣᵖ,
     Sqₗᵖ,
     Sqᵣᵖ,
     ρ,
-    Nₗ,
-    Nᵣ,
+    nₗ,
+    nᵣ,
     qₗ,
     qᵣ,
     ts,
@@ -400,71 +398,71 @@ function compute_warm_precipitation_sources_2M!(
 )
 
     FT = eltype(thp)
-    @. SNₗᵖ = ρ * FT(0)
-    @. SNᵣᵖ = ρ * FT(0)
+    @. Snₗᵖ = ρ * FT(0)
+    @. Snᵣᵖ = ρ * FT(0)
     @. Sqₗᵖ = ρ * FT(0)
     @. Sqᵣᵖ = ρ * FT(0)
 
     # auto-conversion (mass)
     @. Sᵖ = triangle_inequality_limiter(
-        CM2.autoconversion(mp.sb.acnv, mp.sb.pdf_c, qₗ, qᵣ, ρ, Nₗ).dq_rai_dt,
+        CM2.autoconversion(mp.sb.acnv, mp.sb.pdf_c, qₗ, qᵣ, ρ, ρ * nₗ).dq_rai_dt,
         limit(qₗ, dt, 5),
     )
     @. Sqₗᵖ -= Sᵖ
     @. Sqᵣᵖ += Sᵖ
 
     # auto-conversion (number) and liquid self-collection
     @. Sᵖ = triangle_inequality_limiter(
-        CM2.autoconversion(mp.sb.acnv, mp.sb.pdf_c, qₗ, qᵣ, ρ, Nₗ).dN_liq_dt,
-        limit(Nₗ, dt, 10),
+        CM2.autoconversion(mp.sb.acnv, mp.sb.pdf_c, qₗ, qᵣ, ρ, ρ * nₗ).dN_liq_dt / ρ,
+        limit(nₗ, dt, 10),
     )
     @. S₂ᵖ = -triangle_inequality_limiter(
-        -CM2.liquid_self_collection(mp.sb.acnv, mp.sb.pdf_c, qₗ, ρ, Sᵖ),
-        limit(Nₗ, dt, 5)
+        -CM2.liquid_self_collection(mp.sb.acnv, mp.sb.pdf_c, qₗ, ρ, Sᵖ) / ρ,
+        limit(nₗ / ρ, dt, 5)
     )
-    @. SNₗᵖ += Sᵖ
-    @. SNₗᵖ += S₂ᵖ
-    @. SNᵣᵖ -= 0.5*Sᵖ
+    @. Snₗᵖ += Sᵖ
+    @. Snₗᵖ += S₂ᵖ
+    @. Snᵣᵖ -= 0.5*Sᵖ
 
     # rain self-collection and breakup
     @. Sᵖ = -triangle_inequality_limiter(
-        -CM2.rain_self_collection(mp.sb.pdf_r, mp.sb.self, qᵣ, ρ, Nᵣ),
-        limit(Nᵣ, dt, 5),
+        -CM2.rain_self_collection(mp.sb.pdf_r, mp.sb.self, qᵣ, ρ, ρ * nᵣ) / ρ,
+        limit(nᵣ, dt, 5),
     )
     @. S₂ᵖ = triangle_inequality_limiter(
-        CM2.rain_breakup(mp.sb.pdf_r, mp.sb.brek, qᵣ, ρ, Nᵣ, Sᵖ),
-        limit(Nᵣ, dt, 5),
+        CM2.rain_breakup(mp.sb.pdf_r, mp.sb.brek, qᵣ, ρ, ρ * nᵣ, Sᵖ) / ρ,
+        limit(nᵣ, dt, 5),
     )
-    @. SNᵣᵖ += Sᵖ
-    @. SNᵣᵖ += S₂ᵖ
+    @. Snᵣᵖ += Sᵖ
+    @. Snᵣᵖ += S₂ᵖ
 
     # accretion (mass)
     @. Sᵖ = triangle_inequality_limiter(
-        CM2.accretion(mp.sb, qₗ, qᵣ, ρ, Nₗ).dq_rai_dt,
+        CM2.accretion(mp.sb, qₗ, qᵣ, ρ, ρ * nₗ).dq_rai_dt,
         limit(qₗ, dt, 5),
     )
     @. Sqₗᵖ -= Sᵖ
     @. Sqᵣᵖ += Sᵖ
 
     # accretion (number)
     @. Sᵖ = -triangle_inequality_limiter(
-        -CM2.accretion(mp.sb, qₗ, qᵣ, ρ, Nₗ).dN_liq_dt,
-        limit(Nₗ, dt, 5),
+        -CM2.accretion(mp.sb, qₗ, qᵣ, ρ, ρ * nₗ).dN_liq_dt / ρ,
+        limit(nₗ, dt, 5),
     )
-    @. SNₗᵖ += Sᵖ
+    @. Snₗᵖ += Sᵖ
 
     # evaporation (mass)
     @. Sᵖ = -triangle_inequality_limiter(
-        -CM2.rain_evaporation(mp.sb, mp.aps, thp, PP(thp, ts), qᵣ, ρ, Nᵣ, Tₐ(thp, ts)).evap_rate_1,
+        -CM2.rain_evaporation(mp.sb, mp.aps, thp, PP(thp, ts), qᵣ, ρ, ρ * nᵣ, Tₐ(thp, ts)).evap_rate_1,
         limit(qᵣ, dt, 5),
     )
     @. Sqᵣᵖ += Sᵖ
 
     # evaporation (number)
     @. Sᵖ = -triangle_inequality_limiter(
-        -CM2.rain_evaporation(mp.sb, mp.aps, thp, PP(thp, ts), qᵣ, ρ, Nᵣ, Tₐ(thp, ts)).evap_rate_0,
-        limit(Nᵣ, dt, 5),
+        -CM2.rain_evaporation(mp.sb, mp.aps, thp, PP(thp, ts), qᵣ, ρ, ρ * nᵣ, Tₐ(thp, ts)).evap_rate_0 / ρ,
+        limit(nᵣ, dt, 5),
     )
-    @. SNᵣᵖ += Sᵖ
+    @. Snᵣᵖ += Sᵖ
 
 end