add aerosol radiation

CliMA · Jun 26, 2024 · a4b022e · a4b022e
1 parent 6905a73
commit a4b022e
Show file tree

Hide file tree

Showing 6 changed files with 54 additions and 5 deletions.
diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
@@ -108,6 +108,9 @@ idealized_insolation:
 add_isothermal_boundary_layer:
   help: "Add an isothermal boundary layer above the domain top for radiation [`false`, `true` (default)]"
   value: true
+use_aerosols:
+  help: "Use aerosol radiative effects[`false` (default)]"
+  value: false
 dt_cloud_fraction:
   help: "Time between calling cloud fraction update"
   value: "3hours"

diff --git a/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml b/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml
@@ -10,12 +10,13 @@ approximate_linear_solve_iters: 2
 moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
+use_aerosols: true
 idealized_insolation: false
 rayleigh_sponge: true
 non_orographic_gravity_wave: true
 orographic_gravity_wave: "gfdl_restart"
 surface_setup: "DefaultMoninObukhov"
-prescribed_aerosols: ["CB1", "CB2"]
+prescribed_aerosols: ["CB1", "CB2", "SO4"]
 toml: [toml/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.toml]
 diagnostics:
   - short_name: [edt, evu]

diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
@@ -61,6 +61,7 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
     (; params) = p
     (; idealized_insolation, idealized_h2o, idealized_clouds) = p.radiation
     (; ᶠradiation_flux, radiation_model) = p.radiation
+    (; radiation_mode) = p.atmos
 
     # If we have prescribed aerosols, we need to update them
     for (key, tv) in pairs(p.tracers.prescribed_aerosol_timevaryinginputs)
@@ -85,7 +86,7 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
     @. ᶜT =
         min(max(TD.air_temperature(thermo_params, ᶜts), FT(T_min)), FT(T_max))
 
-    if !(radiation_model.radiation_mode isa RRTMGPI.GrayRadiation)
+    if !(radiation_mode isa RRTMGPI.GrayRadiation)
         ᶜrh = Fields.array2field(
             radiation_model.center_relative_humidity,
             axes(Y.c),
@@ -132,8 +133,8 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
     end
 
     if !idealized_clouds && !(
-        radiation_model.radiation_mode isa RRTMGPI.GrayRadiation ||
-        radiation_model.radiation_mode isa RRTMGPI.ClearSkyRadiation
+        radiation_mode isa RRTMGPI.GrayRadiation ||
+        radiation_mode isa RRTMGPI.ClearSkyRadiation
     )
         ᶜΔz = Fields.local_geometry_field(Y.c).∂x∂ξ.components.data.:9
         ᶜlwp = Fields.array2field(
@@ -157,6 +158,15 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
         @. ᶜfrac = cloud_diagnostics_tuple.cf
     end
 
+    if radiation_mode.use_aerosols && !(radiation_mode isa RRTMGPI.GrayRadiation)
+        ᶜΔz = Fields.local_geometry_field(Y.c).∂x∂ξ.components.data.:9
+        ᶜaero_conc = Fields.array2field(
+            radiation_model.center_aerosol_column_mass_density,
+            axes(Y.c),
+        )
+        @. ᶜaero_conc = p.tracers.prescribed_aerosol_fields.:CB1 * ᶜΔz
+    end
+
     set_surface_albedo!(Y, p, t, p.atmos.surface_albedo)
 
     RRTMGPI.update_fluxes!(radiation_model)

diff --git a/src/parameterized_tendencies/radiation/RRTMGPInterface.jl b/src/parameterized_tendencies/radiation/RRTMGPInterface.jl
@@ -23,18 +23,21 @@ struct ClearSkyRadiation <: AbstractRRTMGPMode
     idealized_insolation::Bool
     idealized_clouds::Bool
     add_isothermal_boundary_layer::Bool
+    use_aerosols::Bool
 end
 struct AllSkyRadiation <: AbstractRRTMGPMode
     idealized_h2o::Bool
     idealized_insolation::Bool
     idealized_clouds::Bool
     add_isothermal_boundary_layer::Bool
+    use_aerosols::Bool
 end
 struct AllSkyRadiationWithClearSkyDiagnostics <: AbstractRRTMGPMode
     idealized_h2o::Bool
     idealized_insolation::Bool
     idealized_clouds::Bool
     add_isothermal_boundary_layer::Bool
+    use_aerosols::Bool
 end
 
 """
@@ -787,6 +790,21 @@ function RRTMGPModel(
             )
         end
 
+        if radiation_mode.use_aerosols
+            aero_type = DA{FT}(undef, nlay, ncol)
+            name = "center_aerosol_type"
+            set_and_save!(aero_type, name, t..., dict)
+            aero_size = DA{FT}(undef, nlay, ncol)
+            name = "center_aerosol_radius"
+            set_and_save!(aero_size, name, t..., dict)
+            aero_mass = DA{FT}(undef, nlay, ncol)
+            name = "center_aerosol_column_mass_density"
+            set_and_save!(aero_mass, name, t..., dict)
+            aerosol_state = RRTMGP.AtmosphericStates.AerosolState(aero_type, aero_size, aero_mass)
+        else
+            aerosol_state = nothing
+        end
+
         as = RRTMGP.AtmosphericStates.AtmosphericState(
             lon,
             lat,
@@ -797,7 +815,7 @@ function RRTMGPModel(
             t_sfc,
             vmr,
             cloud_state,
-            nothing,
+            aerosol_state,
         )
     end
 

diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl
@@ -40,6 +40,9 @@ function radiation_model_cache(
     context = ClimaComms.context(axes(Y.c))
     device = context.device
     (; idealized_h2o, idealized_insolation, idealized_clouds) = radiation_mode
+    if !(radiation_mode isa RRTMGPI.GrayRadiation)
+        (; use_aerosols) = radiation_mode
+    end
     FT = Spaces.undertype(axes(Y.c))
     DA = ClimaComms.array_type(device){FT}
     rrtmgp_params = CAP.rrtmgp_params(params)
@@ -180,6 +183,15 @@ function radiation_model_cache(
                     )
                 end
             end
+
+            if use_aerosols
+                kwargs = (;
+                    kwargs...,
+                    center_aerosol_type = 3,
+                    center_aerosol_radius = 0.2,
+                    center_aerosol_column_mass_density = NaN, # initialized in callback
+                )
+            end
         end
 
         if RRTMGPI.requires_z(interpolation) ||

diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
@@ -205,6 +205,8 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
     @assert idealized_clouds in (true, false)
     add_isothermal_boundary_layer = parsed_args["add_isothermal_boundary_layer"]
     @assert add_isothermal_boundary_layer in (true, false)
+    use_aerosols = parsed_args["use_aerosols"]
+    @assert use_aerosols in (true, false)
     radiation_name = parsed_args["rad"]
     @assert radiation_name in (
         nothing,
@@ -222,6 +224,7 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
             idealized_insolation,
             idealized_clouds,
             add_isothermal_boundary_layer,
+            use_aerosols,
         )
     elseif radiation_name == "gray"
         RRTMGPI.GrayRadiation(
@@ -236,13 +239,15 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
             idealized_insolation,
             idealized_clouds,
             add_isothermal_boundary_layer,
+            use_aerosols,
         )
     elseif radiation_name == "allskywithclear"
         RRTMGPI.AllSkyRadiationWithClearSkyDiagnostics(
             idealized_h2o,
             idealized_insolation,
             idealized_clouds,
             add_isothermal_boundary_layer,
+            use_aerosols,
         )
     elseif radiation_name == "DYCOMS"
         RadiationDYCOMS{FT}()