Merge #1244

1244: Add 2-moment, prescribed Nd autoconversions and accretions r=trontrytel a=trontrytel Co-authored-by: lynnyang <[email protected]> # PULL REQUEST ## Purpose and Content This PR makes use of 2-moment, prescribed cloud droplet number concentration microphysics schemes in the new release of CloudMicrophysics.jl. On the TC.jl side we are creating a rain model that accepts the prescribed droplet number concentration as argument. The default scheme is unchanged. The new keyword options are: `--rain_formation_scheme` (valid options: `clima_1m_default`, `KK2000`, `B1994`, `TC1980`, `LD2004`) `--prescribed_Nd` (number of droplets per m3) See the 2-moment CloudMicrophysics.jl [docs](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics2M/) for details. ## Benefits and Risks We will now have more microphysics options. (Uncertain if its a benefit or risk). Co-authored-by: Anna Jaruga <[email protected]>
CliMA · Jul 28, 2022 · d9acc64 · d9acc64 · trontrytel · Jul 28, 2022
2 parents 2ce0936 + 209f180
commit d9acc64
Show file tree

Hide file tree

Showing 12 changed files with 189 additions and 23 deletions.
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -161,6 +161,11 @@ steps:
         command: "julia --color=yes --project=integration_tests integration_tests/driver.jl --case Rico --thermo_covariance_model prognostic --skip_tests true --suffix _prognostic_covarinaces"
         artifact_paths: "Output.Rico.01_prognostic_covarinaces/stats/comparison/*"
 
+      - label: ":partly_sunny: Rico with KK2000 microphysics"
+        command: "julia --color=yes --project=integration_tests integration_tests/driver.jl --case Rico --skip_tests true --precipitation_model clima_1m --rain_formation_scheme KK2000 --prescribed_Nd 2e8 --suffix _two_moment_micro_KK2000"
+        artifact_paths: "Output.Rico.01_prognostic_covarinaces/stats/comparison/*"
+
+
   - group: "Sphere cases"
     steps:
 

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "TurbulenceConvection"
 uuid = "8e072fc4-01f8-44fb-b9dc-f9336c367e6b"
 authors = ["Climate Modeling Alliance"]
-version = "0.34.2"
+version = "0.34.3"
 
 [deps]
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"

diff --git a/driver/generate_namelist.jl b/driver/generate_namelist.jl
@@ -399,6 +399,8 @@ function Rico(namelist_defaults)
     namelist["forcing"]["coriolis"] = 4.5e-5
 
     namelist["microphysics"]["precipitation_model"] = "clima_1m"
+    namelist["microphysics"]["rain_formation_scheme"] = "clima_1m_default"
+    namelist["microphysics"]["prescribed_Nd"] = 1e8
     namelist["microphysics"]["precip_fraction_model"] = "prescribed" # "prescribed" or "cloud_cover"
     namelist["microphysics"]["prescribed_precip_frac_value"] = 1.0
     namelist["microphysics"]["precip_fraction_limiter"] = 0.3
@@ -440,6 +442,8 @@ function TRMM_LBA(namelist_defaults)
     namelist["time_stepping"]["dt_min"] = 1.0
 
     namelist["microphysics"]["precipitation_model"] = "clima_1m" # "cutoff"
+    namelist["microphysics"]["rain_formation_scheme"] = "clima_1m_default"
+    namelist["microphysics"]["prescribed_Nd"] = 1e8
     namelist["microphysics"]["precip_fraction_model"] = "prescribed" # "prescribed" or "cloud_cover"
     namelist["microphysics"]["prescribed_precip_frac_value"] = 1.0
     namelist["microphysics"]["precip_fraction_limiter"] = 0.3
@@ -537,6 +541,8 @@ function DYCOMS_RF02(namelist_defaults)
     namelist["time_stepping"]["dt_min"] = 1.0
 
     namelist["microphysics"]["precipitation_model"] = "clima_1m" #"cutoff"
+    namelist["microphysics"]["rain_formation_scheme"] = "clima_1m_default"
+    namelist["microphysics"]["prescribed_Nd"] = 1e8
     namelist["microphysics"]["precip_fraction_model"] = "prescribed" # "prescribed" or "cloud_cover"
     namelist["microphysics"]["prescribed_precip_frac_value"] = 1.0
     namelist["microphysics"]["precip_fraction_limiter"] = 0.3

diff --git a/driver/main.jl b/driver/main.jl
@@ -12,6 +12,7 @@ import CloudMicrophysics as CM
 import CloudMicrophysics.MicrophysicsNonEq as CMNe
 import CloudMicrophysics.Microphysics0M as CM0
 import CloudMicrophysics.Microphysics1M as CM1
+import CloudMicrophysics.Microphysics2M as CM2
 
 import ClimaCore as CC
 import SciMLBase
@@ -41,7 +42,7 @@ function DiffEqNoiseProcess.wiener_randn!(rng::Random.AbstractRNG, rand_vec::CC.
     parent(rand_vec.face) .= Random.randn.()
 end
 
-struct Simulation1d{IONT, P, A, C, F, R, SM, EDMF, PM, D, TIMESTEPPING, STATS, PS, SRS, LESDK}
+struct Simulation1d{IONT, P, A, C, F, R, SM, EDMF, PM, RFM, D, TIMESTEPPING, STATS, PS, SRS, LESDK}
     io_nt::IONT
     prog::P
     aux::A
@@ -51,6 +52,7 @@ struct Simulation1d{IONT, P, A, C, F, R, SM, EDMF, PM, D, TIMESTEPPING, STATS, P
     surf_params::SM
     edmf::EDMF
     precip_model::PM
+    rain_formation_model::RFM
     diagnostics::D
     TS::TIMESTEPPING
     Stats::STATS
@@ -120,7 +122,29 @@ function Simulation1d(namelist)
         error("Invalid precip_name $(precip_name)")
     end
 
-    edmf = TC.EDMFModel(FTD, namelist, precip_model)
+    if precip_name == "clima_1m"
+        rain_formation_name =
+            TC.parse_namelist(namelist, "microphysics", "rain_formation_scheme"; default = "clima_1m_default")
+        prescribed_Nd = TC.parse_namelist(namelist, "microphysics", "prescribed_Nd"; default = 1e8)
+
+        rain_formation_model = if rain_formation_name == "clima_1m_default"
+            TC.Clima1M_default()
+        elseif rain_formation_name == "KK2000"
+            TC.KK2000(prescribed_Nd)
+        elseif rain_formation_name == "B1994"
+            TC.B1994(prescribed_Nd)
+        elseif rain_formation_name == "TC1980"
+            TC.TC1980(prescribed_Nd)
+        elseif rain_formation_name == "LD2004"
+            TC.LD2004(prescribed_Nd)
+        else
+            error("Invalid rain_formation_name $(rain_formation_name)")
+        end
+    else
+        rain_formation_model = TC.NoRainFormation()
+    end
+
+    edmf = TC.EDMFModel(FTD, namelist, precip_model, rain_formation_model)
     if isbits(edmf)
         @info "edmf = \n$(summary(edmf))"
     else
@@ -204,6 +228,7 @@ function Simulation1d(namelist)
         surf_params,
         edmf,
         precip_model,
+        rain_formation_model,
         diagnostics,
         TS,
         Stats,
@@ -317,6 +342,7 @@ function solve_args(sim::Simulation1d)
         calibrate_io,
         edmf = sim.edmf,
         precip_model = sim.precip_model,
+        rain_formation_model = sim.rain_formation_model,
         param_set = sim.param_set,
         aux = aux,
         io_nt = sim.io_nt,

diff --git a/integration_tests/cli_options.jl b/integration_tests/cli_options.jl
@@ -53,23 +53,35 @@ function parse_commandline()
         default = ""
         "--n_up"           # Number of updrafts
         arg_type = Int
-        "--moisture_model" # Moisture model (equilibrium or nonequilibrium)
+        "--moisture_model"
+        help = "Cloud condensate formation model [`equilibrium` (default), `nonequilibrium`]"
         arg_type = String
-        "--precipitation_model" # Precipitation model (None, cutoff or clima_1m)
+        "--precipitation_model"
+        help = "Precipitation model [`None` (default), `cutoff` or `clima_1m`]"
         arg_type = String
-        "--precip_fraction_model" # Precipitation model (prescribed or cloud_cover)
+        "--rain_formation_scheme"
+        help = "Rain autoconversion and accretion scheme [`clima_1m_default` (default) , `KK2000`, `B1994`, `TC1980`, `LD2004`]"
         arg_type = String
-        "--prescribed_precip_frac_value" # Value of the precipitation fraction, if prescribed
+        "--prescribed_Nd"
+        help = "Prescribed cloud droplet number concentration. Valid when rain_formation_scheme is `KK2000`, `B1994`, `TC1980` or `LD2004`]"
         arg_type = Float64
-        "--precip_fraction_limiter" # Minimum precipitation fraction, if diagnostic
+        "--precip_fraction_model"
+        help = "Assumed (constant with height) precipitation fraction [`prescribed` (default), `cloud_cover`]"
+        arg_type = String
+        "--prescribed_precip_frac_value"
+        help = "Value of the precipitation fraction, if prescribed."
+        arg_type = Float64
+        "--precip_fraction_limiter"
+        help = "Minimum precipitation fraction, if diagnostic."
         arg_type = Float64
-        "--thermo_covariance_model" # covariance model (prognostic or diagnostic)
+        "--thermo_covariance_model"
+        help = "The type of equation for the sgs covariance [`prognostic`, `diagnostic` (default)]"
         arg_type = String
         "--float_type"
         arg_type = String
         default = "Float64"
         "--config"
-        help = "Spatial configuration [`sphere` (default), `column`]"
+        help = "Spatial configuration [`sphere`, `column` (default)]"
         arg_type = String
         default = "column"
         "--set_src_seed"

diff --git a/integration_tests/overwrite_namelist.jl b/integration_tests/overwrite_namelist.jl
@@ -19,6 +19,8 @@ function overwrite_namelist!(namelist, parsed_args)
     "n_up"                         => (nl, pa, key) -> (nl["turbulence"]["EDMF_PrognosticTKE"]["updraft_number"] = pa[key]),
     "moisture_model"               => (nl, pa, key) -> (nl["thermodynamics"]["moisture_model"] = pa[key]),
     "precipitation_model"          => (nl, pa, key) -> (nl["microphysics"]["precipitation_model"] = pa[key]),
+    "rain_formation_scheme"        => (nl, pa, key) -> (nl["microphysics"]["rain_formation_scheme"] = pa[key]),
+    "prescribed_Nd"                => (nl, pa, key) -> (nl["microphysics"]["prescribed_Nd"] = pa[key]),
     "precip_fraction_model"        => (nl, pa, key) -> (nl["microphysics"]["precip_fraction_model"] = pa[key]),
     "prescribed_precip_frac_value" => (nl, pa, key) -> (nl["microphysics"]["prescribed_precip_frac_value"] = pa[key]),
     "precip_fraction_limiter"      => (nl, pa, key) -> (nl["microphysics"]["precip_fraction_limiter"] = pa[key]),

diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl
@@ -4,6 +4,7 @@ function microphysics(
     state::State,
     edmf::EDMFModel,
     precip_model::AbstractPrecipitationModel,
+    rain_formation_model::AbstractRainFormationModel,
     Δt::Real,
     param_set::APS,
 )
@@ -34,6 +35,7 @@ function microphysics(
         mph = precipitation_formation(
             param_set,
             precip_model,
+            rain_formation_model,
             prog_pr.q_rai[k],
             prog_pr.q_sno[k],
             aux_en.area[k],
@@ -72,7 +74,7 @@ function microphysics(
     return nothing
 end
 
-function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt::Real)
+function quad_loop(en_thermo::SGSQuadrature, precip_model, rain_formation_model, vars, param_set, Δt::Real)
 
     env_len = 8
     src_len = 8
@@ -168,8 +170,18 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt:
             q_ice_en = TD.ice_specific_humidity(thermo_params, ts)
             T = TD.air_temperature(thermo_params, ts)
             # autoconversion and accretion
-            mph =
-                precipitation_formation(param_set, precip_model, q_rai, q_sno, subdomain_area, ρ_c, Δt, ts, precip_frac)
+            mph = precipitation_formation(
+                param_set,
+                precip_model,
+                rain_formation_model,
+                q_rai,
+                q_sno,
+                subdomain_area,
+                ρ_c,
+                Δt,
+                ts,
+                precip_frac,
+            )
 
             # environmental variables
             inner_env[i_ql] += q_liq_en * weights[m_h] * sqpi_inv
@@ -231,6 +243,7 @@ function microphysics(
     state::State,
     edmf::EDMFModel,
     precip_model,
+    rain_formation_model,
     Δt::Real,
     param_set::APS,
 )
@@ -284,7 +297,7 @@ function microphysics(
                 precip_frac = precip_fraction,
                 zc = FT(grid.zc[k].z),
             )
-            outer_env, outer_src = quad_loop(en_thermo, precip_model, vars, param_set, Δt)
+            outer_env, outer_src = quad_loop(en_thermo, precip_model, rain_formation_model, vars, param_set, Δt)
 
             # update environmental variables
 
@@ -339,6 +352,7 @@ function microphysics(
             mph = precipitation_formation(
                 param_set,
                 precip_model,
+                rain_formation_model,
                 prog_pr.q_rai[k],
                 prog_pr.q_sno[k],
                 aux_en.area[k],

diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl
@@ -49,6 +49,7 @@ function compute_precipitation_formation_tendencies(
     state::State,
     edmf::EDMFModel,
     precip_model::AbstractPrecipitationModel,
+    rain_formation_model::AbstractRainFormationModel,
     Δt::Real,
     param_set::APS,
 )
@@ -87,6 +88,7 @@ function compute_precipitation_formation_tendencies(
             mph = precipitation_formation(
                 param_set,
                 precip_model,
+                rain_formation_model,
                 prog_pr.q_rai[k],
                 prog_pr.q_sno[k],
                 aux_up[i].area[k],

diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl
@@ -18,6 +18,7 @@ import CloudMicrophysics as CM
 import CloudMicrophysics.MicrophysicsNonEq as CMNe
 import CloudMicrophysics.Microphysics0M as CM0
 import CloudMicrophysics.Microphysics1M as CM1
+import CloudMicrophysics.Microphysics2M as CM2
 import UnPack
 import Random
 import StochasticDiffEq as SDE

diff --git a/src/microphysics_coupling.jl b/src/microphysics_coupling.jl
@@ -47,6 +47,7 @@ Computes the tendencies to qt and θ_liq_ice due to precipitation formation
 function precipitation_formation(
     param_set::APS,
     precip_model::AbstractPrecipitationModel,
+    rain_formation_model::AbstractRainFormationModel,
     qr::FT,
     qs::FT,
     area::FT,
@@ -111,19 +112,82 @@ function precipitation_formation(
             # The saturation adjustment scheme prevents using the
             # 1-moment snow autoconversion rate that assumes
             # that the supersaturation is present in the domain.
-            S_qt_rain = -min(q.liq / Δt, α_acnv * CM1.conv_q_liq_to_q_rai(microphys_params, q.liq))
+            if rain_formation_model isa Clima1M_default
+                S_qt_rain = -min(q.liq / Δt, α_acnv * CM1.conv_q_liq_to_q_rai(microphys_params, q.liq))
+            elseif rain_formation_model isa KK2000
+                S_qt_rain =
+                    -min(
+                        q.liq / Δt,
+                        α_acnv * CM2.conv_q_liq_to_q_rai_KK2000(
+                            microphys_params,
+                            q.liq,
+                            ρ,
+                            N_d = rain_formation_model.prescribed_Nd,
+                        ),
+                    )
+            elseif rain_formation_model isa B1994
+                S_qt_rain =
+                    -min(
+                        q.liq / Δt,
+                        α_acnv * CM2.conv_q_liq_to_q_rai_B1994(
+                            microphys_params,
+                            q.liq,
+                            ρ,
+                            N_d = rain_formation_model.prescribed_Nd,
+                        ),
+                    )
+            elseif rain_formation_model isa TC1980
+                S_qt_rain =
+                    -min(
+                        q.liq / Δt,
+                        α_acnv * CM2.conv_q_liq_to_q_rai_TC1980(
+                            microphys_params,
+                            q.liq,
+                            ρ,
+                            N_d = rain_formation_model.prescribed_Nd,
+                        ),
+                    )
+            elseif rain_formation_model isa LD2004
+                S_qt_rain =
+                    -min(
+                        q.liq / Δt,
+                        α_acnv * CM2.conv_q_liq_to_q_rai_LD2004(
+                            microphys_params,
+                            q.liq,
+                            ρ,
+                            N_d = rain_formation_model.prescribed_Nd,
+                        ),
+                    )
+            else
+                error("Unrecognized rain formation model")
+            end
             S_qt_snow = -min(q.ice / Δt, α_acnv * CM1.conv_q_ice_to_q_sno_no_supersat(microphys_params, q.ice))
             qr_tendency -= S_qt_rain
             qs_tendency -= S_qt_snow
             qt_tendency += S_qt_rain + S_qt_snow
             ql_tendency += S_qt_rain
             qi_tendency += S_qt_snow
+
             θ_liq_ice_tendency -= 1 / Π_m / c_pm * (L_v0 * S_qt_rain + L_s0 * S_qt_snow)
 
             # accretion cloud water + rain
-            S_qr =
-                min(q.liq / Δt, α_accr * CM1.accretion(microphys_params, liq_type, rain_type, q.liq, qr, ρ)) *
-                precip_fraction
+            if rain_formation_model isa Clima1M_default
+                S_qr =
+                    min(q.liq / Δt, α_accr * CM1.accretion(microphys_params, liq_type, rain_type, q.liq, qr, ρ)) *
+                    precip_fraction
+            elseif rain_formation_model isa KK2000
+                S_qr = min(q.liq / Δt, α_accr * CM2.accretion_KK2000(microphys_params, q.liq, qr, ρ)) * precip_fraction
+            elseif rain_formation_model isa B1994
+                S_qr = min(q.liq / Δt, α_accr * CM2.accretion_B1994(microphys_params, q.liq, qr, ρ)) * precip_fraction
+            elseif rain_formation_model isa TC1980
+                S_qr = min(q.liq / Δt, α_accr * CM2.accretion_TC1980(microphys_params, q.liq, qr)) * precip_fraction
+            elseif rain_formation_model isa LD2004
+                S_qr =
+                    min(q.liq / Δt, α_accr * CM1.accretion(microphys_params, liq_type, rain_type, q.liq, qr, ρ)) *
+                    precip_fraction
+            else
+                error("Unrecognized rain formation model")
+            end
             qr_tendency += S_qr
             qt_tendency -= S_qr
             ql_tendency -= S_qr