Update to the new CloudMicrophysics interface

CliMA · Oct 30, 2023 · 4160a99 · 4160a99
1 parent 06e610b
commit 4160a99
Show file tree

Hide file tree

Showing 7 changed files with 35 additions and 142 deletions.
diff --git a/Project.toml b/Project.toml
@@ -53,7 +53,7 @@ CUDA = "4"
 ClimaComms = "0.5.5"
 ClimaCore = "0.10.55"
 ClimaTimeSteppers = "0.7.14"
-CloudMicrophysics = "0.14.0"
+CloudMicrophysics = "0.15.0"
 Colors = "0.12"
 Dierckx = "0.5"
 DiffEqBase = "6"

diff --git a/examples/Project.toml b/examples/Project.toml
@@ -68,6 +68,7 @@ ClimaCorePlots = "0.2"
 ClimaCoreSpectra = "0.1"
 ClimaCoreTempestRemap = "0.3"
 ClimaCoreVTK = "0.7"
+CloudMicrophysics = "0.15"
 ForwardDiff = "0.10"
 Glob = "1"
 JET = "0.7, 0.8"

diff --git a/perf/Project.toml b/perf/Project.toml
@@ -76,6 +76,7 @@ ClimaCorePlots = "0.2"
 ClimaCoreSpectra = "0.1"
 ClimaCoreTempestRemap = "0.3"
 ClimaCoreVTK = "0.7"
+CloudMicrophysics = "0.15"
 DiffEqNoiseProcess = "5"
 ForwardDiff = "0.10"
 Glob = "1"

diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl
@@ -144,127 +144,3 @@ function precipitation_tendency!(
     end
     return nothing
 end
-
-
-#####
-##### 1-Moment coupled to sgs
-#####
-# TODO: move 1-moment microphysics cache / tendency here
-function precipitation_cache(Y, precip_model::Microphysics1Moment)
-    FT = Spaces.undertype(axes(Y.c))
-
-    return (;
-        ᶜS_ρe_tot = similar(Y.c, FT),
-        ᶜS_ρq_tot = similar(Y.c, FT),
-        ᶜS_ρq_rai = similar(Y.c, FT),
-        ᶜS_ρq_sno = similar(Y.c, FT),
-        ᶜS_q_rai_evap = similar(Y.c, FT),
-        ᶜS_q_sno_melt = similar(Y.c, FT),
-        ᶜS_q_sno_sub_dep = similar(Y.c, FT),
-        ᶜq_rai = similar(Y.c, FT),
-        ᶜq_sno = similar(Y.c, FT),
-    )
-end
-
-"""
-Computes the rain and snow advection (down) tendency
-"""
-function precipitation_advection_tendency!(
-    Yₜ,
-    Y,
-    p,
-    colidx,
-    ::Microphysics1Moment,
-)
-    FT = Spaces.undertype(axes(Y.c))
-    (; params) = p
-
-    ρ_c = Y.c.ρ[colidx]
-
-    # helper to calculate the rain velocity
-    # TODO: assuming w_gm = 0
-    # TODO: verify translation
-
-    ρq_rai = Y.c.ρq_rai[colidx]
-    ρq_sno = Y.c.ρq_sno[colidx]
-
-    RB = Operators.RightBiasedC2F(; top = Operators.SetValue(FT(0)))
-    ᶜdivᵥ = Operators.DivergenceF2C(; bottom = Operators.Extrapolate())
-    wvec = Geometry.WVector
-    microphys_params = CAP.microphysics_params(params)
-    rain_type = CM.CommonTypes.RainType()
-    snow_type = CM.CommonTypes.SnowType()
-    velo_type = CM.CommonTypes.Blk1MVelType()
-
-    # TODO - some positivity limiters are needed
-
-    # TODO: need to add horizontal advection + vertical velocity of air
-
-    # TODO: use correct advection operators
-    if !isnothing(Yₜ)
-        @. Yₜ.c.ρq_rai[colidx] += ᶜdivᵥ(
-            wvec(
-                RB(
-                    ρq_rai * CM1.terminal_velocity(
-                        microphys_params,
-                        rain_type,
-                        velo_type,
-                        ρ_c,
-                        ρq_rai / ρ_c,
-                    ),
-                ),
-            ),
-        )
-
-        @. Yₜ.c.ρq_sno[colidx] += ᶜdivᵥ(
-            wvec(
-                RB(
-                    ρq_sno * CM1.terminal_velocity(
-                        microphys_params,
-                        snow_type,
-                        velo_type,
-                        ρ_c,
-                        ρq_sno / ρ_c,
-                    ),
-                ),
-            ),
-        )
-    end
-    return nothing
-end
-
-function precipitation_tendency!(
-    Yₜ,
-    Y,
-    p,
-    t,
-    colidx,
-    precip_model::Microphysics1Moment,
-)
-    (; ᶜS_ρe_tot, ᶜS_ρq_tot, ᶜS_ρq_rai, ᶜS_ρq_sno) = p.precipitation
-    compute_precipitation_cache!(
-        Y,
-        p,
-        colidx,
-        precip_model,
-        p.atmos.turbconv_model,
-    )
-
-    if !isnothing(Yₜ)
-        @. Yₜ.c.ρ[colidx] += ᶜS_ρq_tot[colidx]
-        @. Yₜ.c.ρq_tot[colidx] += ᶜS_ρq_tot[colidx]
-        @. Yₜ.c.ρq_rai[colidx] += ᶜS_ρq_rai[colidx]
-        @. Yₜ.c.ρq_sno[colidx] += ᶜS_ρq_sno[colidx]
-
-        if :ρe_tot in propertynames(Y.c)
-            @. Yₜ.c.ρe_tot[colidx] += ᶜS_ρe_tot[colidx]
-        else
-            error(
-                "1-moment microphysics can only be coupled to ρe_tot energy variable",
-            )
-        end
-    end
-
-    precipitation_advection_tendency!(Yₜ, Y, p, colidx, precip_model)
-    return nothing
-end
diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl
@@ -43,11 +43,13 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP
     updraft_number::Int
 end
 
-Base.@kwdef struct ClimaAtmosParameters{FT, TP, RP, IP, MPP, SFP, TCP} <: ACAP
+Base.@kwdef struct ClimaAtmosParameters{FT, TP, RP, IP, MPP, WP, SFP, TCP} <:
+                   ACAP
     thermodynamics_params::TP
     rrtmgp_params::RP
     insolation_params::IP
     microphysics_params::MPP
+    water_params::WP
     surface_fluxes_params::SFP
     turbconv_params::TCP
     Omega::FT
@@ -94,7 +96,7 @@ for var in [:molmass_ratio, :R_d, :R_v, :cp_d, :cv_v, :cv_l, :cv_d]
 end
 
 # Forwarding CloudMicrophysics parameters
-ρ_cloud_liq(ps::ACAP) = CM.Parameters.ρ_cloud_liq(microphysics_params(ps))
+ρ_cloud_liq(ps::ACAP) = ps.water_params.ρw
 
 # Forwarding SurfaceFluxes parameters
 von_karman_const(ps::ACAP) =

diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
@@ -13,14 +13,6 @@ function create_parameter_set(config::AtmosConfig)
         aliases = string.(fieldnames(param_struct))
         aliases = setdiff(aliases, string.(propertynames(subparam_structs)))
         pairs = CP.get_parameter_values!(toml_dict, aliases)
-        # Workaround for setting τ_precip = dt
-        if parsed_args["override_τ_precip"] &&
-           param_struct == CM.Parameters.CloudMicrophysicsParameters
-            pairs = (;
-                pairs...,
-                τ_precip = FT(CA.time_to_seconds(parsed_args["dt"])),
-            )
-        end
         return param_struct{FT, typeof.(values(subparam_structs))...}(;
             pairs...,
             subparam_structs...,
@@ -30,15 +22,35 @@ function create_parameter_set(config::AtmosConfig)
     (; toml_dict, parsed_args) = config
     FT = CP.float_type(toml_dict)
 
+    # EDMF parameters
+    turbconv_params =
+        create_parameter_struct(CAP.TurbulenceConvectionParameters)
+    # Thermodynamics.jl parameters
     thermo_params =
         create_parameter_struct(TD.Parameters.ThermodynamicsParameters)
+    # Radiation parameters
     rrtmgp_params = create_parameter_struct(RP.RRTMGPParameters)
+    # Insolation.jl parameters
     insolation_params = create_parameter_struct(IP.InsolationParameters)
-    microphys_params = create_parameter_struct(
-        CM.Parameters.CloudMicrophysicsParameters;
-        subparam_structs = (; thermo_params),
-    )
+    # Water properties parameters (from CloudMicrophysics.jl)
+    water_params = CM.Parameters.WaterProperties(FT, toml_dict)
+
+    # Microphysics scheme parameters (from CloudMicrophysics.jl)
+    # TODO - repeating the logic from solver/model_getters.jl...
+    if parsed_args["override_τ_precip"]
+        toml_dict["τ_precip"]["value"] =
+            FT(CA.time_to_seconds(parsed_args["dt"]))
+    end
+    precip_model = parsed_args["precip_model"]
+    microphys_params = if precip_model == nothing || precip_model == "nothing"
+        nothing
+    elseif precip_model == "0M"
+        CM.Parameters.Parameters0M(FT, toml_dict)
+    else
+        error("Invalid precip_model $(precip_model)")
+    end
 
+    # SurfaceFluxes.jl parameters
     aliases = [
         "Pr_0_Businger",
         "a_m_Businger",
@@ -56,20 +68,20 @@ function create_parameter_set(config::AtmosConfig)
         γ = pairs.γ_Businger,
     )
     ufp = UF.BusingerParams{FT}(; pairs...)
-
     surf_flux_params = create_parameter_struct(
         SF.Parameters.SurfaceFluxesParameters;
         subparam_structs = (; ufp, thermo_params),
     )
-    turbconv_params =
-        create_parameter_struct(CAP.TurbulenceConvectionParameters;)
+
+    # Create the big ClimaAtmos parameters struct
     return create_parameter_struct(
         CAP.ClimaAtmosParameters;
         subparam_structs = (;
             thermodynamics_params = thermo_params,
             rrtmgp_params,
             insolation_params,
             microphysics_params = microphys_params,
+            water_params,
             surface_fluxes_params = surf_flux_params,
             turbconv_params,
         ),

diff --git a/test/Project.toml b/test/Project.toml
@@ -73,6 +73,7 @@ ClimaCorePlots = "0.2"
 ClimaCoreSpectra = "0.1"
 ClimaCoreTempestRemap = "0.3"
 ClimaCoreVTK = "0.7"
+CloudMicrophysics = "0.15"
 Dierckx = "0.5"
 DiffEqNoiseProcess = "5"
 ForwardDiff = "0.10"