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 [#](https://github.com/CliMA/RRTMGP.jl/pulls/)
+
 v0.20.0
 ------
 - Add five size bins of dust and sea salt aerosols.

src/optics/AtmosphericStates.jl

@@ -111,11 +111,13 @@ 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}
+    "aerosol optical depth"
+    aod::A
     "aerosol mask, = true if any aerosol is present"
     aero_mask::B
     "aerosol size (microns)"

src/optics/Optics.jl

@@ -238,12 +238,13 @@ Computes optical properties for the longwave problem.
             )
         end
         if !isnothing(lkp_aero)
+            aod = view(as.aerosol_state.aod, :, 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)
             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, τ, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd)
         end
     end
     return nothing
@@ -325,12 +326,13 @@ end
         )
     end
     if !isnothing(lkp_aero)
+        aod = view(as.aerosol_state.aod, :, 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)
         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!(aod, τ, ssa, g, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd)
     end
     return nothing
 end
@@ -420,12 +422,14 @@ end
         )
     end
     if !isnothing(lkp_aero)
+        aod = view(as.aerosol_state.aod, :, 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)
         rel_hum = AtmosphericStates.getview_rel_hum(as, gcol)
 
         add_aerosol_optics_2stream!(
+            aod,
             τ,
             ssa,
             g,

src/optics/aerosol_optics.jl

@@ -1,5 +1,6 @@
 """
     add_aerosol_optics_1scalar!(
+        aod,
         τ,
         aero_mask,
         aero_size,
@@ -12,7 +13,7 @@
 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, τ, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd)
     FT = eltype(τ)
     @inbounds begin
         nlay = length(τ)
@@ -21,6 +22,7 @@ to the exising `OneScalar` optics properties.
                 τ_aero, τ_ssa_aero, τ_ssag_aero =
                     compute_lookup_aerosol(lkp_aero, ibnd, aero_mass, aero_size, rel_hum[glay], glay)
                 τ[glay] += (τ_aero - τ_ssa_aero)
+                aod[glay] += τ_aero
             end
         end
     end
@@ -29,6 +31,7 @@ end
 
 """
     add_aerosol_optics_2stream!(
+        aod,
         τ,
         ssa,
         g,
@@ -42,9 +45,12 @@ end
     )
 
 This function computes the `TwoStream` aerosol optics properties and adds them
-to the exising `TwoStream` optics properties.
+to the exising `TwoStream` optics properties:
+
+ - `aod` aerosol optical depth
 """
 @inline function add_aerosol_optics_2stream!(
+    aod,
     τ,
     ssa,
     g,
@@ -69,6 +75,7 @@ to the exising `TwoStream` optics properties.
                     τ_aero, ssa_aero, g_aero = delta_scale(τ_aero, ssa_aero, g_aero)
                 end
                 τ[glay], ssa[glay], g[glay] = increment_2stream(τ[glay], ssa[glay], g[glay], τ_aero, ssa_aero, g_aero)
+                aod[glay] += τ_aero
             end
         end
     end

test/read_all_sky_with_aerosols.jl

@@ -103,6 +103,7 @@ 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(zeros(FT, nlay, ncol)),
         DA{Bool}(undef, nlay, ncol),
         DA(repeat(aero_size, 1, 1, nrepeat)[:, :, 1:ncol]),
         DA(repeat(aero_mass, 1, 1, nrepeat)[:, :, 1:ncol]),