microphysics ahhhhhhh

Author: szy21

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -754,18 +754,20 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
             # 0-moment microphysics: sink of q_tot from precipitation removal
             if microphysics_model isa Microphysics0Moment
                 @. S_q_totʲ_prev_level = q_tot_0M_precipitation_sources(
-                    thermo_params,
                     microphys_0m_params,
                     dt,
                     q_totʲ_prev_level,
-                    tsʲ_prev_level,
+                    TD.total_specific_humidity(thermo_params, tsʲ_prev_level),
+                    TD.liquid_specific_humidity(thermo_params, tsʲ_prev_level),
+                    TD.ice_specific_humidity(thermo_params, tsʲ_prev_level),
                 )
                 # 1-moment microphysics: cloud water (liquid and ice) and
                 # precipitation (rain and snow) sources. q_tot is constant, because
                 # all the species are considered a part of the working fluid.
             elseif moisture_model isa NonEquilMoistModel &&
                    microphysics_model isa Microphysics1Moment
                 # Rain formation from the updrafts
+                Tʲ_prev_level = @. lazy(TD.air_temperature(thermo_params, tsʲ_prev_level))
                 compute_precipitation_sources!(
                     Sᵖ_prev_level,
                     Sᵖ_snow_prev_level,
@@ -779,7 +781,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     q_iceʲ_prev_level,
                     q_raiʲ_prev_level,
                     q_snoʲ_prev_level,
-                    tsʲ_prev_level,
+                    Tʲ_prev_level,
                     dt,
                     microphys_1m_params,
                     thermo_params,
@@ -795,7 +797,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     q_iceʲ_prev_level,
                     q_raiʲ_prev_level,
                     q_snoʲ_prev_level,
-                    tsʲ_prev_level,
+                    Tʲ_prev_level,
                     dt,
                     microphys_1m_params,
                     thermo_params,
@@ -810,7 +812,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     q_raiʲ_prev_level,
                     q_snoʲ_prev_level,
                     ρʲ_prev_level,
-                    TD.air_temperature(thermo_params, tsʲ_prev_level),
+                    Tʲ_prev_level,
                     dt,
                 )
                 @. S_q_iceʲ_prev_level += cloud_sources(
@@ -822,7 +824,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     q_raiʲ_prev_level,
                     q_snoʲ_prev_level,
                     ρʲ_prev_level,
-                    TD.air_temperature(thermo_params, tsʲ_prev_level),
+                    Tʲ_prev_level,
                     dt,
                 )
 
@@ -986,14 +988,19 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     mseʲ_prev_level,
                 )
             if microphysics_model isa Microphysics0Moment
+                ᶜTʲ = @. lazy(TD.air_temperature(thermo_params, tsʲ_prev_level))
+                ᶜq_liq_raiʲ = @. lazy(TD.liquid_specific_humidity(thermo_params, tsʲ_prev_level))
+                ᶜq_ice_snoʲ = @. lazy(TD.ice_specific_humidity(thermo_params, tsʲ_prev_level))
                 @. ρaʲu³ʲ_datamse += microphysics_sources(
                     local_geometry_halflevel.J,
                     local_geometry_prev_level.J,
                     ρaʲ_prev_level,
                     S_q_totʲ_prev_level *
                     e_tot_0M_precipitation_sources_helper(
                         thermo_params,
-                        tsʲ_prev_level,
+                        ᶜTʲ,
+                        ᶜq_liq_raiʲ,
+                        ᶜq_ice_snoʲ,
                         Φ_prev_level,
                     ),
                 )
@@ -1403,19 +1410,19 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
     t,
     microphysics_model::Microphysics0Moment,
 )
-    thermo_params = CAP.thermodynamics_params(p.params)
     microphys_0m_params = CAP.microphysics_0m_params(p.params)
     (; dt) = p
-    (; ᶜts, ᶜSqₜᵖ⁰) = p.precomputed
+    (; ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, ᶜSqₜᵖ⁰) = p.precomputed
 
     # Environment precipitation sources (to be applied to grid mean)
     ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
     @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources(
-        thermo_params,
         microphys_0m_params,
         dt,
         ᶜq_tot,
-        ᶜts,
+        ᶜq_tot_safe,
+        ᶜq_liq_rai,
+        ᶜq_ice_sno,
     )
     return nothing
 end
@@ -1430,7 +1437,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
     cloud_params = CAP.microphysics_cloud_params(p.params)
     (; dt) = p
 
-    (; ᶜts, ᶜT, ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed
+    (; ᶜT, ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed
     ᶜSᵖ = p.scratch.ᶜtemp_scalar
     ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2
 
@@ -1448,7 +1455,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
         specific.(Y.c.ρq_ice, Y.c.ρ),
         specific.(Y.c.ρq_rai, Y.c.ρ),
         specific.(Y.c.ρq_sno, Y.c.ρ),
-        ᶜts,
+        ᶜT,
         dt,
         microphys_1m_params,
         thermo_params,
@@ -1464,7 +1471,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
         specific.(Y.c.ρq_ice, Y.c.ρ),
         specific.(Y.c.ρq_rai, Y.c.ρ),
         specific.(Y.c.ρq_sno, Y.c.ρ),
-        ᶜts,
+        ᶜT,
         dt,
         microphys_1m_params,
         thermo_params,

src/cache/precipitation_precomputed_quantities.jl

@@ -56,7 +56,7 @@ function set_precipitation_velocities!(
     microphysics_model::Microphysics1Moment,
     _,
 )
-    (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
+    (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜT, ᶜu) = p.precomputed
     (; ᶜΦ) = p.core
     cmc = CAP.microphysics_cloud_params(p.params)
     cmp = CAP.microphysics_1m_params(p.params)
@@ -99,10 +99,10 @@ function set_precipitation_velocities!(
         ) / Y.c.ρ
     @. ᶜwₕhₜ =
         Geometry.WVector(
-            ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) +
-            ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) +
-            ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) +
-            ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))),
+            ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) +
+            ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) +
+            ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) +
+            ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))),
         ) / Y.c.ρ
     return nothing
 end
@@ -284,7 +284,7 @@ function set_precipitation_velocities!(
     microphysics_model::Microphysics2Moment,
     _,
 )
-    (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₙₗ, ᶜwₙᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
+    (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₙₗ, ᶜwₙᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜT, ᶜu) = p.precomputed
     (; ᶜΦ) = p.core
 
     cmc = CAP.microphysics_cloud_params(p.params)
@@ -359,10 +359,10 @@ function set_precipitation_velocities!(
         ) / Y.c.ρ
     @. ᶜwₕhₜ =
         Geometry.WVector(
-            ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) +
-            ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) +
-            ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) +
-            ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))),
+            ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) +
+            ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) +
+            ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) +
+            ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))),
         ) / Y.c.ρ
     return nothing
 end
@@ -555,7 +555,7 @@ function set_precipitation_velocities!(
     Y, p, ::NonEquilMoistModel, ::Microphysics2MomentP3,
 )
     ## liquid quantities (2M warm rain)
-    (; ᶜwₗ, ᶜwᵣ, ᶜwnₗ, ᶜwnᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed
+    (; ᶜwₗ, ᶜwᵣ, ᶜwnₗ, ᶜwnᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜT, ᶜu) = p.precomputed
     (; ᶜΦ) = p.core
 
     (; ρ, ρq_liq, ρn_liq, ρq_rai, ρn_rai) = Y.c
@@ -606,9 +606,9 @@ function set_precipitation_velocities!(
     @. ᶜwₜqₜ = Geometry.WVector(ᶜwₗ * ρq_liq + ᶜwᵢ * ρq_ice + ᶜwᵣ * ρq_rai) / ρ
     @. ᶜwₕhₜ =
         Geometry.WVector(
-            ᶜwₗ * ρq_liq * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) +
-            ᶜwᵢ * ρq_ice * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) +
-            ᶜwᵣ * ρq_rai * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))),
+            ᶜwₗ * ρq_liq * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) +
+            ᶜwᵢ * ρq_ice * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) +
+            ᶜwᵣ * ρq_rai * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))),
         ) / ρ
     return nothing
 end
@@ -624,22 +624,23 @@ sources from the sub-domains.
 set_precipitation_cache!(Y, p, _, _) = nothing
 function set_precipitation_cache!(Y, p, ::Microphysics0Moment, _)
     (; params, dt) = p
-    (; ᶜts) = p.precomputed
+    (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed
     (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precomputed
     (; ᶜΦ) = p.core
     cm_params = CAP.microphysics_0m_params(params)
     thermo_params = CAP.thermodynamics_params(params)
     @. ᶜS_ρq_tot =
         Y.c.ρ * q_tot_0M_precipitation_sources(
-            thermo_params,
             cm_params,
             dt,
             Y.c.ρq_tot / Y.c.ρ,
-            ᶜts,
+            ᶜq_tot_safe,
+            ᶜq_liq_rai,
+            ᶜq_ice_sno,
         )
     @. ᶜS_ρe_tot =
         ᶜS_ρq_tot *
-        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, ᶜΦ)
+        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶜΦ)
     return nothing
 end
 function set_precipitation_cache!(
@@ -653,7 +654,7 @@ function set_precipitation_cache!(
     (; ᶜΦ) = p.core
     (; ᶜSqₜᵖ⁰, ᶜSqₜᵖʲs, ᶜρaʲs) = p.precomputed
     (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precomputed
-    (; ᶜts, ᶜtsʲs) = p.precomputed
+    (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, ᶜtsʲs) = p.precomputed
     thermo_params = CAP.thermodynamics_params(p.params)
 
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
@@ -663,15 +664,20 @@ function set_precipitation_cache!(
     @. ᶜS_ρe_tot =
         ᶜSqₜᵖ⁰ *
         ρ *
-        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, ᶜΦ)
+        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶜΦ)
     for j in 1:n
+        ᶜTʲ = @. lazy(TD.air_temperature(thermo_params, ᶜtsʲs.:($$j)))
+        ᶜq_liq_raiʲ = @. lazy(TD.liquid_specific_humidity(thermo_params, ᶜtsʲs.:($$j)))
+        ᶜq_ice_snoʲ = @. lazy(TD.ice_specific_humidity(thermo_params, ᶜtsʲs.:($$j)))
         @. ᶜS_ρq_tot += ᶜSqₜᵖʲs.:($$j) * ᶜρaʲs.:($$j)
         @. ᶜS_ρe_tot +=
             ᶜSqₜᵖʲs.:($$j) *
             ᶜρaʲs.:($$j) *
             e_tot_0M_precipitation_sources_helper(
                 thermo_params,
-                ᶜtsʲs.:($$j),
+                ᶜTʲ,
+                ᶜq_liq_raiʲ,
+                ᶜq_ice_snoʲ,
                 ᶜΦ,
             )
     end
@@ -693,26 +699,34 @@ function set_precipitation_cache!(
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
 
     @. ᶜS_ρq_tot = ᶜSqₜᵖ⁰ * ᶜρa⁰
+    ᶜT⁰ = @. lazy(TD.air_temperature(thermo_params, ᶜts⁰))
+    ᶜq_liq_rai⁰ = @. lazy(TD.liquid_specific_humidity(thermo_params, ᶜts⁰))
+    ᶜq_ice_sno⁰ = @. lazy(TD.ice_specific_humidity(thermo_params, ᶜts⁰))
     @. ᶜS_ρe_tot =
         ᶜSqₜᵖ⁰ *
         ᶜρa⁰ *
-        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts⁰, ᶜΦ)
+        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜT⁰, ᶜq_liq_rai⁰, ᶜq_ice_sno⁰, ᶜΦ)
     for j in 1:n
+        ᶜTʲ = @. lazy(TD.air_temperature(thermo_params, ᶜtsʲs.:($$j)))
+        ᶜq_liq_raiʲ = @. lazy(TD.liquid_specific_humidity(thermo_params, ᶜtsʲs.:($$j)))
+        ᶜq_ice_snoʲ = @. lazy(TD.ice_specific_humidity(thermo_params, ᶜtsʲs.:($$j)))
         @. ᶜS_ρq_tot += ᶜSqₜᵖʲs.:($$j) * Y.c.sgsʲs.:($$j).ρa
         @. ᶜS_ρe_tot +=
             ᶜSqₜᵖʲs.:($$j) *
             Y.c.sgsʲs.:($$j).ρa *
             e_tot_0M_precipitation_sources_helper(
                 thermo_params,
-                ᶜtsʲs.:($$j),
+                ᶜTʲ,
+                ᶜq_liq_raiʲ,
+                ᶜq_ice_snoʲ,
                 ᶜΦ,
             )
     end
     return nothing
 end
 function set_precipitation_cache!(Y, p, ::Microphysics1Moment, _)
     (; dt) = p
-    (; ᶜts, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed
+    (; ᶜT, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed
     (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed
 
     ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
@@ -744,7 +758,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics1Moment, _)
         ᶜq_ice,
         ᶜq_rai,
         ᶜq_sno,
-        ᶜts,
+        ᶜT,
         dt,
         cmp,
         thp,
@@ -761,7 +775,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics1Moment, _)
         ᶜq_ice,
         ᶜq_rai,
         ᶜq_sno,
-        ᶜts,
+        ᶜT,
         dt,
         cmp,
         thp,
@@ -790,8 +804,7 @@ function set_precipitation_cache!(
 end
 function set_precipitation_cache!(Y, p, ::Microphysics2Moment, _)
     (; dt) = p
-    (; ᶜts) = p.precomputed
-    (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed
+    (; ᶜT, ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed
     (; ᶜSnₗᵖ, ᶜSnᵣᵖ) = p.precomputed
 
     ᶜSᵖ = p.scratch.ᶜtemp_scalar
@@ -818,7 +831,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics2Moment, _)
         lazy.(specific.(Y.c.ρq_ice, Y.c.ρ)),
         lazy.(specific.(Y.c.ρq_rai, Y.c.ρ)),
         lazy.(specific.(Y.c.ρq_sno, Y.c.ρ)),
-        ᶜts,
+        ᶜT,
         dt,
         cmp,
         thp,
@@ -894,7 +907,7 @@ function set_precipitation_surface_fluxes!(
     p,
     microphysics_model::Microphysics0Moment,
 )
-    (; ᶜts, ᶜT) = p.precomputed  # assume ᶜts has been updated
+    (; ᶜT) = p.precomputed
     (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precomputed
     (; surface_rain_flux, surface_snow_flux) = p.precomputed
     (; col_integrated_precip_energy_tendency) = p.conservation_check

src/cache/precomputed_quantities.jl

@@ -16,7 +16,10 @@ The following grid-scale quantities are treated implicitly and are precomputed:
     - `ᶜu`: covariant velocity on cell centers
     - `ᶠu`: contravariant velocity on cell faces
     - `ᶜK`: kinetic energy on cell centers
-    - `ᶜts`: thermodynamic state on cell centers
+    - `ᶜT`: air temperature on cell centers
+    - `ᶜq_tot_safe`: total water specific humidity on cell centers
+    - `ᶜq_liq_rai`: liquid water specific humidity on cell centers
+    - `ᶜq_ice_sno`: ice specific humidity on cell centers
     - `ᶜp`: air pressure on cell centers
 If the `turbconv_model` is `PrognosticEDMFX`, there also two SGS versions of
 every quantity except for `ᶜp` (which is shared across all subdomains):
@@ -675,8 +678,8 @@ NVTX.@annotate function set_implicit_precomputed_quantities!(Y, p, t)
         # quantities of the form ᶜρaχ⁰ / ᶜρ⁰ and ᶜρaχʲ / ᶜρʲ to ᶜK, rather than
         # quantities of the form ᶜρaχ⁰ / ᶜρ and ᶜρaχʲ / ᶜρ. However, we cannot
         # compute ᶜρ⁰ and ᶜρʲ without first computing ᶜts⁰ and ᶜtsʲ, both of
-        # which depend on the value of ᶜp, which in turn depends on ᶜts. Since
-        # ᶜts depends on ᶜK, this
+        # which depend on the value of ᶜp, which in turn depends on ᶜT. Since
+        # ᶜT depends on ᶜK, this
         # means that the amount by which ᶜK needs to be incremented is a
         # function of ᶜK itself. So, unless we run a nonlinear solver here, this
         # circular dependency will prevent us from computing the exact value of