Merge pull request #3631 from CliMA/ck/rrtmgp_aero

add AOD to diagnostics

Author: szy21

.buildkite/Manifest-v1.11.toml

@@ -1983,9 +1983,9 @@ version = "1.11.0"
 
 [[deps.RRTMGP]]
 deps = ["Adapt", "Artifacts", "ClimaComms", "DocStringExtensions", "Random"]
-git-tree-sha1 = "d258e9495bb6a1b7ad2b28683c9382c0053dafc1"
+git-tree-sha1 = "d2318b99ce2d40d563d63c401f2bd264a76e0f41"
 uuid = "a01a1ee8-cea4-48fc-987c-fc7878d79da1"
-version = "0.20.1"
+version = "0.21.0"
 weakdeps = ["CUDA", "ClimaParams"]
 
     [deps.RRTMGP.extensions]

Project.toml

@@ -59,7 +59,7 @@ Logging = "1"
 NCDatasets = "0.14.2"
 NullBroadcasts = "0.1"
 NVTX = "0.3"
-RRTMGP = "0.20.1"
+RRTMGP = "0.21"
 Random = "1"
 SciMLBase = "2.12"
 StaticArrays = "1.7"

config/longrun_configs/amip_target_diagedmf.yml

@@ -34,7 +34,7 @@ t_end: "120days"
 toml: [toml/longrun_aquaplanet_diagedmf.toml]
 netcdf_output_at_levels: true
 diagnostics:
-  - short_name: [mmrso4, mmrdust, mmrss, o3]
+  - short_name: [mmrso4, mmrdust, mmrss, o3, od550aer, odsc550aer]
     reduction_time: average
     period: "1months"
 

config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml

@@ -23,6 +23,6 @@ reproducibility_test: true
 prescribed_aerosols: ["CB1", "CB2", "SO4"]
 toml: [toml/sphere_aquaplanet_1M.toml]
 diagnostics:
-  - short_name: [edt, evu, mmrso4, mmrbcpo, mmrbcpi, o3]
+  - short_name: [edt, evu, mmrso4, mmrbcpo, mmrbcpi, o3, od550aer, odsc550aer]
     reduction_time: average
     period: "1days"

src/diagnostics/radiation_diagnostics.jl

@@ -942,3 +942,88 @@ add_diagnostic_variable!(
     comments = "In-cloud ratio of the third moment over the second moment of the particle size distribution. Set to zero outside of clouds.",
     compute! = compute_reffcli!,
 )
+
+###
+# Aerosol extinction optical depth (2d)
+###
+compute_od550aer!(out, state, cache, time) =
+    compute_od550aer!(out, state, cache, time, cache.atmos.radiation_mode)
+compute_od550aer!(_, _, _, _, radiation_mode::T) where {T} =
+    error_diagnostic_variable("od550aer", radiation_mode)
+
+function compute_od550aer!(
+    out,
+    state,
+    cache,
+    time,
+    radiation_mode::Union{
+        RRTMGPI.AllSkyRadiationWithClearSkyDiagnostics,
+        RRTMGPI.AllSkyRadiation,
+        RRTMGPI.ClearSkyRadiation,
+    },
+)
+    @assert cache.atmos.radiation_mode.aerosol_radiation "aerosol_radiation must be true to enable aerosol optical depth diagnostics"
+    FT = eltype(state)
+    if isnothing(out)
+        return Fields.array2field(
+            cache.radiation.rrtmgp_model.aod_sw_extinction,
+            axes(Fields.level(state.f, half)),
+        )
+    else
+        out .= Fields.array2field(
+            cache.radiation.rrtmgp_model.aod_sw_extinction,
+            axes(Fields.level(state.f, half)),
+        )
+    end
+end
+
+add_diagnostic_variable!(
+    short_name = "od550aer",
+    long_name = "Ambient Aerosol Optical Thickness at 550nm",
+    standard_name = "atmosphere_optical_thickness_due_to_ambient_aerosol_particles",
+    units = "",
+    comments = "Aerosol optical depth from the ambient aerosols at wavelength 550 nm",
+    compute! = compute_od550aer!,
+)
+
+###
+# Aerosol scattering optical depth (2d)
+###
+compute_odsc550aer!(out, state, cache, time) =
+    compute_odsc550aer!(out, state, cache, time, cache.atmos.radiation_mode)
+compute_odsc550aer!(_, _, _, _, radiation_mode::T) where {T} =
+    error_diagnostic_variable("odsc550aer", radiation_mode)
+
+function compute_odsc550aer!(
+    out,
+    state,
+    cache,
+    time,
+    radiation_mode::Union{
+        RRTMGPI.AllSkyRadiationWithClearSkyDiagnostics,
+        RRTMGPI.AllSkyRadiation,
+        RRTMGPI.ClearSkyRadiation,
+    },
+)
+    @assert cache.atmos.radiation_mode.aerosol_radiation "aerosol_radiation must be true to enable aerosol optical depth diagnostics"
+    FT = eltype(state)
+    if isnothing(out)
+        return Fields.array2field(
+            cache.radiation.rrtmgp_model.aod_sw_scattering,
+            axes(Fields.level(state.f, half)),
+        )
+    else
+        out .= Fields.array2field(
+            cache.radiation.rrtmgp_model.aod_sw_scattering,
+            axes(Fields.level(state.f, half)),
+        )
+    end
+end
+
+add_diagnostic_variable!(
+    short_name = "odsc550aer",
+    long_name = "Ambient Scattering Aerosol Optical Thickness at 550nm",
+    units = "",
+    comments = "Aerosol scattering optical depth from the ambient aerosols at wavelength 550 nm",
+    compute! = compute_odsc550aer!,
+)

src/parameterized_tendencies/radiation/RRTMGPInterface.jl

@@ -824,7 +824,11 @@ function RRTMGPModel(
         end
 
         if radiation_mode.aerosol_radiation
+            aod_sw_ext = DA{FT}(undef, ncol)
+            aod_sw_sca = DA{FT}(undef, ncol)
             aero_mask = DA{Bool}(undef, nlay, ncol)
+            set_and_save!(aod_sw_ext, "aod_sw_extinction", t..., dict)
+            set_and_save!(aod_sw_sca, "aod_sw_scattering", t..., dict)
 
             n_aerosol_sizes = maximum(values(idx_aerosize))
             n_aerosols = length(idx_aerosol)
@@ -882,6 +886,8 @@ function RRTMGPModel(
                 )
             end
             aerosol_state = RRTMGP.AtmosphericStates.AerosolState(
+                aod_sw_ext,
+                aod_sw_sca,
                 aero_mask,
                 aero_size,
                 aero_mass,

src/parameterized_tendencies/radiation/radiation.jl

@@ -262,6 +262,8 @@ function radiation_model_cache(
                 else
                     kwargs = (;
                         kwargs...,
+                        aod_sw_extinction = NaN,
+                        aod_sw_scattering = NaN,
                         # assuming fixed aerosol radius
                         center_dust1_radius = 0.55,
                         center_dust2_radius = 1.4,

test/restart.jl

@@ -246,8 +246,8 @@ function test_restart(test_dict; job_id, comms_ctx, more_ignore = Symbol[])
 
     simulation_restarted = CA.get_simulation(config_should_be_same)
 
-    if pkgversion(CA.RRTMGP) < v"0.21"
-        # Versions of RRTMGP older than 0.21 have a bug and do not set the
+    if pkgversion(CA.RRTMGP) < v"0.22"
+        # Versions of RRTMGP older than 0.22 have a bug and do not set the
         # flux_dn_dir, so that face_clear_sw_direct_flux_dn and
         # face_sw_direct_flux_dn is uninitialized and not deterministic
         rrtmgp_clear_fix =