Merge pull request #611 from CliMA/aj/delete_pp_where_possible

Wip on deleting PP from microphysics interface

Author: trontrytel

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.25.0"
+version = "0.26.0"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
@@ -32,5 +32,5 @@ MLJ = "0.20"
 QuadGK = "2.9"
 RootSolvers = "0.3, 0.4"
 SpecialFunctions = "1, 2"
-Thermodynamics = "0.12.13"
+Thermodynamics = "0.12.14"
 julia = "1.6"

docs/src/plots/Microphysics1M_plots.jl

@@ -64,10 +64,9 @@ MK.set_theme!(MK.theme_minimal())
 
 # autoconversion rate figure
 T = 273.15
-q_range = TDI.PP.(q_tot, 0.0, q_ice_range)
 fig = MK.Figure()
 ax = MK.Axis(fig[1, 1]; xlabel = "q_liq or q_ice [g/kg]", ylabel = "autoconversion rate [1/s]", limits)
-q_ice_to_q_sno_rate = T -> CM1.conv_q_ice_to_q_sno.(ice, aps, tps, q_range, q_rai, q_sno, ρ_air, T)
+q_ice_to_q_sno_rate = T -> CM1.conv_q_ice_to_q_sno.(ice, aps, tps, q_tot, FT(0), q_ice_range, q_rai, q_sno, ρ_air, T)
 MK.lines!(q_liq_range * 1e3, CM1.conv_q_liq_to_q_rai.(rain.acnv1M, q_liq_range), label = "Rain")
 MK.lines!(q_ice_range * 1e3, q_ice_to_q_sno_rate(T - 5), label = "Snow T = −5°C")
 MK.lines!(q_ice_range * 1e3, q_ice_to_q_sno_rate(T - 10), label = "Snow T = −10°C")
@@ -144,15 +143,27 @@ q_vap = 0.15 * q_sat
 q_ice = 0.0
 q_liq = q_tot - q_vap - q_ice
 q_sno = 0.0
-q = TDI.PP(q_tot, q_liq, q_ice)
 R = TDI.Rₘ(tps, q_tot, q_liq, q_ice)
 ρ = p / R / T
 
 fig = MK.Figure()
 ax = MK.Axis(fig[1, 1]; xlabel = "q_rain [g/kg]", ylabel = "rain evaporation rate [1/s]")
 MK.xlims!(ax, 0, q_max * 1e3)
 MK.lines!(
-    q_rain_range * 1e3, CM1.evaporation_sublimation.(rain, Blk1MVel.rain, aps, tps, Ref(q), q_rain_range, q_sno, ρ, T),
+    q_rain_range * 1e3,
+    CM1.evaporation_sublimation.(
+        rain,
+        Blk1MVel.rain,
+        aps,
+        tps,
+        q_tot,
+        q_liq .- q_rain_range,
+        q_ice,
+        q_rain_range,
+        q_sno,
+        ρ,
+        T,
+    ),
     label = "ClimateMachine",
 )
 MK.lines!(q_rain_range * 1e3, rain_evap_empirical.(tps, q_rain_range, q_tot, q_liq, T, p, ρ), label = "empirical")
@@ -173,12 +184,25 @@ q_snow_range = range(1e-8, stop = 5e-3, length = 100)
 q_tot = 15e-3
 q_vap = 0.15 * q_sat
 q_liq = 0.0
+q_lcl = 0.0
 q_ice = q_tot - q_vap - q_liq
 q_rai = 0.0
-q = TDI.PP(q_tot, q_liq, q_ice)
 R = TDI.Rₘ(tps, q_tot, q_liq, q_ice)
 ρ = p / R / T
-rate = CM1.evaporation_sublimation.(snow, Blk1MVel.snow, aps, tps, Ref(q), q_rai, q_snow_range, ρ, T)
+rate =
+    CM1.evaporation_sublimation.(
+        snow,
+        Blk1MVel.snow,
+        aps,
+        tps,
+        q_tot,
+        q_lcl,
+        q_ice .- q_snow_range,
+        q_rai,
+        q_snow_range,
+        ρ,
+        T,
+    )
 MK.lines!(q_snow_range * 1e3, rate, label = "T < 0°C")
 
 # example values 2
@@ -190,12 +214,25 @@ q_snow_range = range(1e-8, stop = 5e-3, length = 100)
 q_tot = 15e-3
 q_vap = 0.15 * q_sat
 q_liq = 0.0
+q_lcl = 0.0
 q_ice = q_tot - q_vap - q_liq
 q_rai = 0.0
-q = TDI.PP(q_tot, q_liq, q_ice)
 R = TDI.Rₘ(tps, q_tot, q_liq, q_ice)
 ρ = p / R / T
-rate = CM1.evaporation_sublimation.(snow, Blk1MVel.snow, aps, tps, Ref(q), q_rai, q_snow_range, ρ, T)
+rate =
+    CM1.evaporation_sublimation.(
+        snow,
+        Blk1MVel.snow,
+        aps,
+        tps,
+        q_tot,
+        q_lcl,
+        q_ice .- q_snow_range,
+        q_rai,
+        q_snow_range,
+        ρ,
+        T,
+    )
 MK.lines!(q_snow_range * 1e3, rate, label = "T > 0°C")
 
 MK.axislegend(ax; position = :rt)

src/Microphysics0M.jl

@@ -11,8 +11,8 @@ module Microphysics0M
 
 import CloudMicrophysics.Parameters as CMP
 
-# Only needed for the wrapper for calling remove_precipitation with
-# TD.PhasePartition as an argument. Can be dropped if we drop this pattern.
+# Needed for the wrapper for calling the 0-moment remove_precipitation
+# with TD.PhasePartition in ClimaAtmos
 import CloudMicrophysics.ThermodynamicsInterface as TDI
 
 export remove_precipitation
@@ -47,6 +47,9 @@ remove_precipitation(
 ###
 ### Wrappers for calling with TD.PhasePartition
 ###
+### For now leaving the PhasePartition wrapper because I'm not sure how to get
+### rid of equilibrium thermo state in the Atmos model.
+###
 remove_precipitation(params::CMP.Parameters0M, q::TDI.TD.PhasePartition) =
     remove_precipitation(params, q.liq, q.ice)
 remove_precipitation(params::CMP.Parameters0M, q::TDI.TD.PhasePartition, q_vap_sat) =

src/Microphysics1M.jl

@@ -646,24 +646,4 @@ function snow_melt(
     return snow_melt_rate
 end
 
-###
-### Wrappers for calling with TD.PhasePartition
-###
-
-conv_q_ice_to_q_sno(
-    ice_params::CMP.CloudIce, aps::CMP.AirProperties, tps::TDI.PS,
-    q::TDI.TD.PhasePartition, q_rai, q_sno, ρ, T,
-) = conv_q_ice_to_q_sno(ice_params, aps, tps, TDI.q_(q, q_rai, q_sno)..., q_rai, q_sno, ρ, T)
-
-evaporation_sublimation(
-    rain_params::CMP.Rain, vel::CMP.Blk1MVelTypeRain, aps::CMP.AirProperties, tps::TDI.PS,
-    q::TDI.TD.PhasePartition, q_rai, q_sno, ρ, T,
-) = evaporation_sublimation(rain_params, vel, aps, tps, TDI.q_(q, q_rai, q_sno)..., q_rai, q_sno, ρ, T)
-
-evaporation_sublimation(
-    snow_params::CMP.Snow, vel::CMP.Blk1MVelTypeSnow, aps::CMP.AirProperties, tps::TDI.PS,
-    q::TDI.TD.PhasePartition, q_rai, q_sno, ρ, T,
-) = evaporation_sublimation(snow_params, vel, aps, tps, TDI.q_(q, q_rai, q_sno)..., q_rai, q_sno, ρ, T)
-
-
 end #module Microphysics1M.jl

src/Microphysics2M.jl

@@ -932,14 +932,4 @@ function accretion((; accr)::CMP.TC1980{FT}, q_liq, q_rai) where {FT}
     return A * q_liq * q_rai
 end
 
-###
-### Wrappers for calling with TD.PhasePartition
-###
-
-rain_evaporation(
-    sb_params::CMP.SB2006, air_params::CMP.AirProperties, tps::TDI.PS,
-    q::TDI.TD.PhasePartition,
-    q_rai, q_sno, ρ, N_rai, T,
-) = rain_evaporation(sb_params, air_params, tps, q.tot, q.liq, q.ice, q_rai, q_sno, ρ, N_rai, T)
-
 end # module

src/MicrophysicsNonEq.jl

@@ -170,18 +170,4 @@ function terminal_velocity(
     return fall_w
 end
 
-###
-### Wrappers fro calling with TD.PhasePartition
-###
-conv_q_vap_to_q_liq_ice_MM2015(
-    prs::CMP.CloudLiquid, tps::TDI.PS, q::TDI.TD.PhasePartition, q_rai, q_sno, ρ, T,
-    #) = conv_q_vap_to_q_liq_ice_MM2015(prs, tps, TDI.q_(q, q_rai, q_sno)..., ρ, T)
-) = conv_q_vap_to_q_liq_ice_MM2015(prs, tps, q.tot, q.liq - q_rai, q.ice - q_sno, q_rai, q_sno..., ρ, T)
-
-conv_q_vap_to_q_liq_ice_MM2015(
-    prs::CMP.CloudIce, tps::TDI.PS, q::TDI.TD.PhasePartition, q_rai, q_sno, ρ, T,
-) = conv_q_vap_to_q_liq_ice_MM2015(prs, tps, q.tot, q.liq - q_rai, q.ice - q_sno, q_rai, q_sno..., ρ, T)
-#) = conv_q_vap_to_q_liq_ice_MM2015(prs, tps, TDI.q_(q, q_rai, q_sno)..., ρ, T)
-
-
 end #module MicrophysicsNonEq.jl

src/ThermodynamicsInterface.jl

@@ -2,48 +2,29 @@ module ThermodynamicsInterface
 
 import Thermodynamics as TD
 const PS = TD.Parameters.ThermodynamicsParameters
-const PP = TD.PhasePartition
-
-
-# Constants
+const PP = TD.PhasePartition  # Still used in Atmos EquilMoist model rain formation
 
+###
+### Constants
+###
 grav(tps::PS) = TD.Parameters.grav(tps) # Needed in parcel model
 
-
-# Constants for moist air
-
 Rᵥ(tps::PS) = TD.Parameters.R_v(tps)
 Rd(tps::PS) = TD.Parameters.R_d(tps)
 Rd_over_Rv(tps::PS) = 1 / TD.Parameters.molmass_ratio(tps)
-# Gas constant for mois air
-Rₘ(tps::PS, qₜ, qₗ, qᵢ) = TD.gas_constant_air(tps, PP(qₜ, qₗ, qᵢ)) #TODO - PP
+
+# Gas constant for moist air
+Rₘ(tps::PS, qₜ, qₗ, qᵢ) = TD.gas_constant_air(tps, qₜ, qₗ, qᵢ)
 
 Lᵥ(tps::PS, T) = TD.latent_heat_vapor(tps, T)
 Lₛ(tps::PS, T) = TD.latent_heat_sublim(tps, T)
 Lf(tps::PS, T) = TD.latent_heat_fusion(tps, T)
 
 cpₘ(tps::PS, qₜ, qₗ, qᵢ) = TD.cp_m(tps, qₜ, qₗ, qᵢ)
 
-# Supersaturations
-
-saturation_vapor_pressure_over_liquid(tps::PS, T) =
-    TD.saturation_vapor_pressure(tps, T, TD.Liquid())
-saturation_vapor_pressure_over_ice(tps::PS, T) =
-    TD.saturation_vapor_pressure(tps, T, TD.Ice())
-
-# (only used in tests)
-saturation_vapor_specific_content_over_liquid(tps::PS, T, ρ) =
-    TD.q_vap_saturation_generic(tps, T, ρ, TD.Liquid())
-saturation_vapor_specific_content_over_ice(tps::PS, T, ρ) =
-    TD.q_vap_saturation_generic(tps, T, ρ, TD.Ice())
-
-supersaturation_over_liquid(tps::PS, qₜ, qₗ, qᵢ, ρ, T) =
-    TD.supersaturation(tps, PP(qₜ, qₗ, qᵢ), ρ, T, TD.Liquid()) # TODO - PP
-supersaturation_over_ice(tps::PS, qₜ, qₗ, qᵢ, ρ, T) =
-    TD.supersaturation(tps, PP(qₜ, qₗ, qᵢ), ρ, T, TD.Ice()) # TODO - PP
-
-
-# Utility functions
+###
+### Utility functions
+###
 
 # Get vapor specific content from total water, total liquid water and total ice
 # water specific contents.
@@ -53,17 +34,41 @@ q_vap(q_tot, q_liq, q_ice) = q_tot - q_liq - q_ice
 q_vap(q_tot, q_liq, q_ice, q_rai, q_sno) = q_tot - q_liq - q_ice - q_rai - q_sno
 
 # Get specific content from partial pressure
-p2q(tps::PS, T, ρ, p) = TD.q_vap_saturation_from_density(tps, T, ρ, p)
+p2q(tps::PS, T, ρ, pᵥ) = TD.q_vap_saturation_from_density(tps, T, ρ, pᵥ)
 
+# Get partial pressure from specific content
+q2p(tps::PS, T, ρ, qᵥ) = qᵥ * ρ * Rᵥ(tps) * T
 
 # Get air density from temperature, pressure and water content
 # (only used in tests)
 air_density(tps::PS, T, p, q_tot, q_liq, q_ice) =
     TD.air_density(tps, T, p, PP(q_tot, q_liq, q_ice))
 
-# Get a tuple containing total water, cloud liquid water and cloud ice specific
-# contents from TD.PhasePartition and rain and snow specific contents.
-# Assumes that q::PhasePartition = (q_tot, q_cloud_liq + q_rai, q_cloud_ice + q_sno)
-q_(q::PP, q_rai, q_sno) = (q.tot, q.liq - q_rai, q.ice - q_sno)
+###
+### Supersaturations
+###
+saturation_vapor_pressure_over_liquid(tps::PS, T) =
+    TD.saturation_vapor_pressure(tps, T, TD.Liquid())
+saturation_vapor_pressure_over_ice(tps::PS, T) =
+    TD.saturation_vapor_pressure(tps, T, TD.Ice())
+
+# (only used in tests)
+saturation_vapor_specific_content_over_liquid(tps::PS, T, ρ) =
+    TD.q_vap_saturation_generic(tps, T, ρ, TD.Liquid())
+saturation_vapor_specific_content_over_ice(tps::PS, T, ρ) =
+    TD.q_vap_saturation_generic(tps, T, ρ, TD.Ice())
+
+function supersaturation_over_liquid(tps::PS, qₜ, qₗ, qᵢ, ρ, T)
+    pᵥ_sat = saturation_vapor_pressure_over_liquid(tps, T)
+    qᵥ = q_vap(qₜ, qₗ, qᵢ)
+    pᵥ = q2p(tps, T, ρ, qᵥ)
+    return pᵥ / pᵥ_sat - 1
+end
+function supersaturation_over_ice(tps::PS, qₜ, qₗ, qᵢ, ρ, T)
+    pᵥ_sat = saturation_vapor_pressure_over_ice(tps, T)
+    qᵥ = q_vap(qₜ, qₗ, qᵢ)
+    pᵥ = q2p(tps, T, ρ, qᵥ)
+    return pᵥ / pᵥ_sat - 1
+end
 
 end

test/aerosol_activation_emulators.jl

@@ -300,7 +300,9 @@ function test_emulator(
     w = FT(0.5)    # vertical velocity m/s
     p_vs = TDI.saturation_vapor_pressure_over_liquid(tps, T)
     q_vs = 1 / (1 - 1 / TDI.Rd_over_Rv(tps) * (p_vs - p) / p_vs)
-    q = TDI.TD.PhasePartition(q_vs)
+    q_tot = q_vs
+    q_liq = FT(0)
+    q_ice = FT(0)
 
     # Aerosol size distribution
     salt = CMP.Seasalt(FT)
@@ -323,14 +325,14 @@ function test_emulator(
         machine_name = machine_name,
     )
 
-    TT.@test AA.N_activated_per_mode(mach, ap, ad, aip, tps, T, p, w, q)[1] ≈
-             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[1] rtol =
+    TT.@test AA.N_activated_per_mode(mach, ap, ad, aip, tps, T, p, w, q_tot, q_liq, q_ice)[1] ≈
+             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q_tot, q_liq, q_ice)[1] rtol =
         rtols[1]
-    TT.@test AA.N_activated_per_mode(mach, ap, ad, aip, tps, T, p, w, q)[2] ≈
-             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[2] rtol =
+    TT.@test AA.N_activated_per_mode(mach, ap, ad, aip, tps, T, p, w, q_tot, q_liq, q_ice)[2] ≈
+             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q_tot, q_liq, q_ice)[2] rtol =
         rtols[2]
-    TT.@test AA.total_N_activated(mach, ap, ad, aip, tps, T, p, w, q) ≈
-             AA.total_N_activated(ap, ad, aip, tps, T, p, w, q) rtol = rtols[2]
+    TT.@test AA.total_N_activated(mach, ap, ad, aip, tps, T, p, w, q_tot, q_liq, q_ice) ≈
+             AA.total_N_activated(ap, ad, aip, tps, T, p, w, q_tot, q_liq, q_ice) rtol = rtols[2]
 
 end
 

test/gpu_tests.jl

@@ -117,15 +117,6 @@ end
 
     @inbounds begin
         output[1, i] = CMN.conv_q_vap_to_q_liq_ice_MM2015(
-            liquid,
-            tps,
-            TDI.TD.PhasePartition(qᵥ_sl[i]),
-            FT(0),
-            FT(0),
-            ρ[i],
-            T[i],
-        )
-        output[2, i] = CMN.conv_q_vap_to_q_liq_ice_MM2015(
             liquid,
             tps,
             qᵥ_sl[i],
@@ -136,7 +127,7 @@ end
             ρ[i],
             T[i],
         )
-        output[3, i] = CMN.conv_q_vap_to_q_liq_ice(ice, qᵢ_s[i], qᵢ[i])
+        output[2, i] = CMN.conv_q_vap_to_q_liq_ice(ice, qᵢ_s[i], qᵢ[i])
     end
 end
 
@@ -343,7 +334,6 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        q = TDI.TD.PhasePartition(FT(qt[i]), FT(ql[i]), FT(0))
 
         output[1, i] =
             CM2.autoconversion_and_liquid_self_collection(
@@ -857,7 +847,7 @@ function test_gpu(FT)
     end
 
     TT.@testset "non-equilibrium microphysics kernels" begin
-        dims = (3, 1)
+        dims = (2, 1)
         (; output, ndrange) = setup_output(dims, FT)
 
         ρ = ArrayType([FT(0.8)])
@@ -871,8 +861,7 @@ function test_gpu(FT)
 
         # test that nonequilibrium cloud formation is callable and returns a reasonable value
         TT.@test Array(output)[1] ≈ FT(3.763783850665844e-5)
-        TT.@test Array(output)[2] ≈ FT(3.763783850665844e-5)
-        TT.@test Array(output)[3] ≈ FT(-1e-4)
+        TT.@test Array(output)[2] ≈ FT(-1e-4)
     end
 
     TT.@testset "Chen 2022 terminal velocity kernels" begin