Add aerosol optical depth to AerosolState

NEWS.md

@@ -4,6 +4,11 @@ RRTMGP.jl Release Notes
 main
 ------
 
+v0.21.0
+------
+- Aerosol optical depth was added to the `AerosolState`
+  PR [#567](https://github.com/CliMA/RRTMGP.jl/pulls/567)
+
 v0.20.1
 ------
 - Clip effective radius by the look-up table range

Project.toml

@@ -1,7 +1,7 @@
 name = "RRTMGP"
 uuid = "a01a1ee8-cea4-48fc-987c-fc7878d79da1"
 authors = ["Climate Modeling Alliance"]
-version = "0.20.1"
+version = "0.21.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/optics/AtmosphericStates.jl

@@ -111,11 +111,15 @@ end
 Adapt.@adapt_structure CloudState
 
 """
-    AerosolState{B, D}
+    AerosolState{A, B, D}
 
 Aerosol state, used to compute optical properties.
 """
-struct AerosolState{B, D}
+struct AerosolState{A, B, D}
+    "shortwave aerosol optical depth"
+    aod_sw_ext::A
+    "shortwave aerosol optical depth (scattering component)"
+    aod_sw_sca::A
     "aerosol mask, = true if any aerosol is present"
     aero_mask::B
     "aerosol size (microns)"

src/optics/LookUpTables.jl

@@ -859,6 +859,8 @@ struct LookUpAerosolMerra{D, D1, D2, D3, D4, W} <: AbstractLookUp
     organic_carbon::D2
     "beginning and ending wavenumber for each band (`2, nband`) cm⁻¹"
     bnd_lims_wn::W
+    "Band number index corresponding to 550 nm"
+    iband_550nm::Int
 end
 Adapt.@adapt_structure LookUpAerosolMerra
 
@@ -873,7 +875,13 @@ function LookUpAerosolMerra(ds, ::Type{FT}, ::Type{DA}) where {FT <: AbstractFlo
     black_carbon = DA{FT}(Array(ds["aero_bcar_tbl"]))
     organic_carbon_rh = DA{FT}(Array(ds["aero_ocar_rh_tbl"]))
     organic_carbon = DA{FT}(Array(ds["aero_ocar_tbl"]))
-    bnd_lims_wn = DA{FT}(Array(ds["bnd_limits_wavenumber"]))
+    bnd_lims_wn = Array(ds["bnd_limits_wavenumber"])
+
+    iband_550nm = findfirst(1:size(bnd_lims_wn, 2)) do i
+        1 / (bnd_lims_wn[2, i] * 100) ≤ 550 * 1e-9 ≤ 1 / (bnd_lims_wn[1, i] * 100) # bnd_lims_wn is in cm⁻¹
+    end
+    iband_550nm = isnothing(iband_550nm) ? 0 : iband_550nm
+
     # map aerosol name to aerosol idx
     idx_aerosol = Dict(
         "dust1" => 1,
@@ -903,7 +911,8 @@ function LookUpAerosolMerra(ds, ::Type{FT}, ::Type{DA}) where {FT <: AbstractFlo
         black_carbon,
         organic_carbon_rh,
         organic_carbon,
-        bnd_lims_wn,
+        DA{FT}(bnd_lims_wn),
+        iband_550nm,
     ),
     idx_aerosol,
     idx_aerosize

src/optics/Optics.jl

@@ -238,12 +238,26 @@ Computes optical properties for the longwave problem.
             )
         end
         if !isnothing(lkp_aero)
+            aod_sw_ext = nothing
+            aod_sw_sca = nothing
+            iband_550nm = nothing
             aero_mask = view(as.aerosol_state.aero_mask, :, gcol)
             aero_size = view(as.aerosol_state.aero_size, :, :, gcol)
             aero_mass = view(as.aerosol_state.aero_mass, :, :, gcol)
             rel_hum = AtmosphericStates.getview_rel_hum(as, gcol)
 
-            add_aerosol_optics_1scalar!(τ, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd)
+            add_aerosol_optics_1scalar!(
+                τ,
+                aod_sw_ext,
+                aod_sw_sca,
+                aero_mask,
+                aero_size,
+                aero_mass,
+                rel_hum,
+                lkp_aero,
+                ibnd,
+                iband_550nm,
+            )
         end
     end
     return nothing
@@ -325,12 +339,28 @@ end
         )
     end
     if !isnothing(lkp_aero)
+        aod_sw_ext = nothing
+        aod_sw_sca = nothing
+        iband_550nm = nothing
         aero_mask = view(as.aerosol_state.aero_mask, :, gcol)
         aero_size = view(as.aerosol_state.aero_size, :, :, gcol)
         aero_mass = view(as.aerosol_state.aero_mass, :, :, gcol)
         rel_hum = AtmosphericStates.getview_rel_hum(as, gcol)
 
-        add_aerosol_optics_2stream!(τ, ssa, g, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd)
+        add_aerosol_optics_2stream!(
+            τ,
+            ssa,
+            g,
+            aod_sw_ext,
+            aod_sw_sca,
+            aero_mask,
+            aero_size,
+            aero_mass,
+            rel_hum,
+            lkp_aero,
+            ibnd,
+            iband_550nm,
+        )
     end
     return nothing
 end
@@ -420,6 +450,9 @@ end
         )
     end
     if !isnothing(lkp_aero)
+        (; iband_550nm) = lkp_aero
+        aod_sw_ext = view(as.aerosol_state.aod_sw_ext, gcol)
+        aod_sw_sca = view(as.aerosol_state.aod_sw_sca, gcol)
         aero_mask = view(as.aerosol_state.aero_mask, :, gcol)
         aero_size = view(as.aerosol_state.aero_size, :, :, gcol)
         aero_mass = view(as.aerosol_state.aero_mass, :, :, gcol)
@@ -429,12 +462,15 @@ end
             τ,
             ssa,
             g,
+            aod_sw_ext,
+            aod_sw_sca,
             aero_mask,
             aero_size,
             aero_mass,
             rel_hum,
             lkp_aero,
             ibnd,
+            iband_550nm,
             delta_scaling = true,
         )
     end

src/optics/aerosol_optics.jl

@@ -1,26 +1,47 @@
 """
     add_aerosol_optics_1scalar!(
         τ,
+        aod_ext,
+        aod_sca,
         aero_mask,
         aero_size,
         aero_mass,
         rel_hum,
         lkp_aero,
         ibnd,
+        iband_550nm
     )
 
 This function computes the `OneScalar` aerosol optics properties and adds them
 to the exising `OneScalar` optics properties.
 """
-@inline function add_aerosol_optics_1scalar!(τ, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd)
+@inline function add_aerosol_optics_1scalar!(
+    τ,
+    aod_ext,
+    aod_sca,
+    aero_mask,
+    aero_size,
+    aero_mass,
+    rel_hum,
+    lkp_aero,
+    ibnd,
+    iband_550nm,
+)
     FT = eltype(τ)
     @inbounds begin
+        collect_aod = !isnothing(aod_ext) && (ibnd == iband_550nm)
         nlay = length(τ)
+        collect_aod && (aod_ext[] = 0)
+        collect_aod && (aod_sca[] = 0)
         for glay in 1:nlay
             if aero_mask[glay]
                 τ_aero, τ_ssa_aero, τ_ssag_aero =
                     compute_lookup_aerosol(lkp_aero, ibnd, aero_mass, aero_size, rel_hum[glay], glay)
                 τ[glay] += (τ_aero - τ_ssa_aero)
+                if collect_aod
+                    aod_ext[] += τ_aero
+                    aod_sca[] += τ_ssa_aero
+                end
             end
         end
     end
@@ -32,39 +53,55 @@ end
         τ,
         ssa,
         g,
+        aod_ext,
+        aod_sca,
         aero_mask,
         aero_size,
         aero_mass,
         rel_hum,
         lkp_aero,
-        ibnd;
+        ibnd,
+        iband_550nm;
         delta_scaling = false,
     )
 
 This function computes the `TwoStream` aerosol optics properties and adds them
-to the exising `TwoStream` optics properties.
+to the exising `TwoStream` optics properties:
+
+ - `aod_ext` total aerosol optical depth
+ - `aod_sca` scattering component of aerosol optical depth
 """
 @inline function add_aerosol_optics_2stream!(
     τ,
     ssa,
     g,
+    aod_ext,
+    aod_sca,
     aero_mask,
     aero_size,
     aero_mass,
     rel_hum,
     lkp_aero,
-    ibnd;
+    ibnd,
+    iband_550nm;
     delta_scaling = false,
 )
     @inbounds begin
+        collect_aod = !isnothing(aod_ext) && (ibnd == iband_550nm)
         nlay = length(τ)
         FT = eltype(τ)
+        collect_aod && (aod_ext[] = 0)
+        collect_aod && (aod_sca[] = 0)
         for glay in 1:nlay
             if aero_mask[glay]
                 τ_aero, τ_ssa_aero, τ_ssag_aero =
                     compute_lookup_aerosol(lkp_aero, ibnd, aero_mass, aero_size, rel_hum[glay], glay)
                 g_aero = τ_ssag_aero / max(eps(FT), τ_ssa_aero)
                 ssa_aero = τ_ssa_aero / max(eps(FT), τ_aero)
+                if collect_aod
+                    aod_ext[] += τ_aero
+                    aod_sca[] += τ_ssa_aero
+                end
                 if delta_scaling # delta scaling is applied for shortwave problem
                     τ_aero, ssa_aero, g_aero = delta_scale(τ_aero, ssa_aero, g_aero)
                 end

test/all_sky_with_aerosols_utils.jl

@@ -188,5 +188,9 @@ function all_sky_with_aerosols(
     @test max_err_flux_dn_sw ≤ toler_sw[FT]
     @test max_err_flux_net_sw ≤ toler_sw[FT]
 
+    @test minimum(as.aerosol_state.aod_sw_ext) >= 0
+    @test minimum(as.aerosol_state.aod_sw_sca) >= 0
+    @test minimum(as.aerosol_state.aod_sw_ext .- as.aerosol_state.aod_sw_sca) >= 0
+
     return nothing
 end

test/read_all_sky_with_aerosols.jl

@@ -103,6 +103,8 @@ function setup_allsky_with_aerosols_as(
     # repeat the input data to set problem size to ncols
     nrepeat = Int(cld(ncol, ncol_ref))
     aerosol_state = AerosolState(
+        DA{FT}(undef, ncol),
+        DA{FT}(undef, ncol),
         DA{Bool}(undef, nlay, ncol),
         DA(repeat(aero_size, 1, 1, nrepeat)[:, :, 1:ncol]),
         DA(repeat(aero_mass, 1, 1, nrepeat)[:, :, 1:ncol]),