Merge pull request #3743 from CliMA/aj/ml_number_of_droplets

calibratable number_conc for r_eff

Author: trontrytel

.buildkite/Manifest-v1.11.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.4"
 manifest_format = "2.0"
-project_hash = "04945e2a7116cfacb2d57340f477dbe23e208227"
+project_hash = "dbb001c3c134c9962d9a283b94c060c2d2eb8b6a"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"
@@ -395,9 +395,9 @@ version = "0.2.12"
 
 [[deps.ClimaParams]]
 deps = ["TOML"]
-git-tree-sha1 = "bf8fdee01cf5ea7e5e9fa00ac32367f65448484e"
+git-tree-sha1 = "8bf3c3ca9073ae455682eb4f011f94a86a0f88a4"
 uuid = "5c42b081-d73a-476f-9059-fd94b934656c"
-version = "0.10.23"
+version = "0.10.24"
 
 [[deps.ClimaReproducibilityTests]]
 deps = ["OrderedCollections", "PrettyTables"]

Project.toml

@@ -45,7 +45,7 @@ AtmosphericProfilesLibrary = "0.1.7"
 ClimaComms = "0.6.6"
 ClimaCore = "0.14.24"
 ClimaDiagnostics = "0.2.12"
-ClimaParams = "0.10.23"
+ClimaParams = "0.10.24"
 ClimaTimeSteppers = "0.8.2"
 ClimaUtilities = "0.1.22"
 CloudMicrophysics = "0.22.9"

src/callbacks/callbacks.jl

@@ -206,18 +206,53 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
                 max(cloud_fraction, eps(FT))
             @. ᶜfrac = cloud_fraction
             # RRTMGP needs effective radius in microns
+
+            seasalt_aero_conc = p.scratch.ᶜtemp_scalar
+            dust_aero_conc = p.scratch.ᶜtemp_scalar_2
+            SO4_aero_conc = p.scratch.ᶜtemp_scalar_3
+            @. seasalt_aero_conc = 0
+            @. dust_aero_conc = 0
+            @. SO4_aero_conc = 0
+            # Get aerosol mass concentrations if available
+            seasalt_names = [:SSLT01, :SSLT02, :SSLT03, :SSLT04, :SSLT05]
+            dust_names = [:DST01, :DST02, :DST03, :DST04, :DST05]
+            SO4_names = [:SO4]
+            if :prescribed_aerosols_field in propertynames(p.tracers)
+                aerosol_field = p.tracers.prescribed_aerosols_field
+                for aerosol_name in propertynames(aerosol_field)
+                    if aerosol_name in seasalt_names
+                        data = getproperty(aerosol_field, aerosol_name)
+                        @. seasalt_aero_conc += data
+                    elseif aerosol_name in dust_names
+                        data = getproperty(aerosol_field, aerosol_name)
+                        @. dust_aero_conc += data
+                    elseif aerosol_name in SO4_names
+                        data = getproperty(aerosol_field, aerosol_name)
+                        @. SO4_aero_conc += data
+                    end
+                end
+            end
+
             @. ᶜreliq = ifelse(
                 cloud_liquid_water_content > FT(0),
                 CM.CloudDiagnostics.effective_radius_Liu_Hallet_97(
                     cmc.liquid,
                     Y.c.ρ,
                     cloud_liquid_water_content / max(eps(FT), cloud_fraction),
-                    cmc.N_cloud_liquid_droplets,
+                    ml_N_cloud_liquid_droplets(
+                        (cmc,),
+                        dust_aero_conc,
+                        seasalt_aero_conc,
+                        SO4_aero_conc,
+                        cloud_liquid_water_content /
+                        max(eps(FT), cloud_fraction),
+                    ),
                     FT(0),
                     FT(0),
                 ) * m_to_um_factor,
                 FT(0),
             )
+
             @. ᶜreice = ifelse(
                 cloud_ice_water_content > FT(0),
                 CM.CloudDiagnostics.effective_radius_const(cmc.ice) *

src/parameterized_tendencies/microphysics/microphysics_wrappers.jl

@@ -34,6 +34,31 @@ function limit(q, dt, n::Int)
     return q / float(dt) / n
 end
 
+"""
+    ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq)
+
+ - cmc - a struct with cloud and aerosol parameters
+ - c_dust, c_seasalt, c_SO4 - dust, seasalt and ammonium sulfate mass concentrations [kg/kg]
+ - q_liq - liquid water specific humidity
+
+Returns the liquid cloud droplet number concentration diagnosed based on the
+aerosol loading and cloud liquid water.
+"""
+function ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq)
+    # We can also add w, T, RH, w' ...
+    # Also consider lookind only at around cloud base height
+    (; α_dust, α_seasalt, α_SO4, α_q_liq) = cmc.aml
+    (; c₀_dust, c₀_seasalt, c₀_SO4) = cmc.aml
+    N₀ = cmc.N_cloud_liquid_droplets
+
+    FT = eltype(N₀)
+    return N₀ +
+           α_dust * log(max(c_dust, eps(FT)) / c₀_dust) +
+           α_seasalt * log(max(c_seasalt, eps(FT)) / c₀_seasalt) +
+           α_SO4 * log(max(c_SO4, eps(FT)) / c₀_SO4) +
+           α_q_liq * log(max(q_liq, eps(FT)))
+end
+
 """
     cloud_sources(cm_params, thp, ts, dt)
 

src/parameters/create_parameters.jl

@@ -121,6 +121,7 @@ cloud_parameters(toml_dict::CP.AbstractTOMLDict) = (;
         "prescribed_cloud_droplet_number_concentration",
         "ClimaAtmos",
     ).prescribed_cloud_droplet_number_concentration,
+    aml = aerosol_ml_parameters(toml_dict),
 )
 
 microphys_1m_parameters(::Type{FT}) where {FT <: AbstractFloat} =
@@ -147,6 +148,19 @@ function vert_diff_parameters(toml_dict)
     return CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos")
 end
 
+function aerosol_ml_parameters(toml_dict)
+    name_map = (;
+        :prescribed_cloud_droplet_number_concentration => :N₀,
+        :dust_calibration_coefficient => :α_dust,
+        :seasalt_calibration_coefficient => :α_seasalt,
+        :ammonium_sulfate_calibration_coefficient => :α_SO4,
+        :liquid_water_specific_humidity_calibration_coefficent => :α_q_liq,
+        :reference_dust_aerosol_mass_concentration => :c₀_dust,
+        :reference_seasalt_aerosol_mass_concentration => :c₀_seasalt,
+        :reference_ammonium_sulfate_mass_concentration => :c₀_SO4,
+    )
+    return CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos")
+end
 
 to_svec(x::AbstractArray) = SA.SVector{length(x)}(x)
 to_svec(x) = x