Fixes

src/optics/AtmosphericStates.jl

@@ -116,10 +116,10 @@ Adapt.@adapt_structure CloudState
 Aerosol state, used to compute optical properties.
 """
 struct AerosolState{IA, A, B, D}
-    "shortwave aerosol optical depth band numbers"
-    aod_sw_band_number::IA
     "shortwave aerosol optical depth"
-    aod_sw::A
+    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/Optics.jl

@@ -327,14 +327,15 @@ end
         )
     end
     if !isnothing(lkp_aero)
-        aod_sw = nothing
+        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, aod_sw, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd, iband_550nm)
+        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
@@ -425,7 +426,8 @@ end
     end
     if !isnothing(lkp_aero)
         (; iband_550nm) = lkp_aero
-        aod_sw = view(as.aerosol_state.aod_sw, :, gcol)
+        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)
@@ -435,7 +437,8 @@ end
             τ,
             ssa,
             g,
-            aod_sw,
+            aod_sw_ext,
+            aod_sw_sca,
             aero_mask,
             aero_size,
             aero_mass,

src/optics/aerosol_optics.jl

@@ -1,7 +1,8 @@
 """
     add_aerosol_optics_1scalar!(
         τ,
-        aod,
+        aod_ext,
+        aod_sca,
         aero_mask,
         aero_size,
         aero_mass,
@@ -14,19 +15,21 @@
 This function computes the `OneScalar` aerosol optics properties and adds them
 to the exising `OneScalar` optics properties.
 """
-@inline function add_aerosol_optics_1scalar!(τ, aod, aero_mask, aero_size, aero_mass, rel_hum, lkp_aero, ibnd, iband_550nm)
+@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) && (ibnd == iband_550nm)
+        collect_aod = !isnothing(aod_ext) && (ibnd == iband_550nm)
         nlay = length(τ)
-        collect_aod && (aod[] = 0)
+        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[] += τ_aero
+                    aod_ext[] += τ_aero
+                    aod_sca[] += τ_ssa_aero
                 end
             end
         end
@@ -39,7 +42,8 @@ end
         τ,
         ssa,
         g,
-        aod,
+        aod_ext,
+        aod_sca,
         aero_mask,
         aero_size,
         aero_mass,
@@ -53,13 +57,15 @@ end
 This function computes the `TwoStream` aerosol optics properties and adds them
 to the exising `TwoStream` optics properties:
 
- - `aod` aerosol optical depth
+ - `aod_ext` total aerosol optical depth
+ - `aod_sca` scattering component of aerosol optical depth
 """
 @inline function add_aerosol_optics_2stream!(
     τ,
     ssa,
     g,
-    aod,
+    aod_ext,
+    aod_sca,
     aero_mask,
     aero_size,
     aero_mass,
@@ -70,10 +76,11 @@ to the exising `TwoStream` optics properties:
     delta_scaling = false,
 )
     @inbounds begin
-        collect_aod = !isnothing(aod) && (ibnd == iband_550nm)
+        collect_aod = !isnothing(aod_ext) && (ibnd == iband_550nm)
         nlay = length(τ)
         FT = eltype(τ)
-        collect_aod && (aod[] = 0)
+        collect_aod && (aod_ext[] = 0)
+        collect_aod && (aod_sca[] = 0)
         for glay in 1:nlay
             if aero_mask[glay]
                 τ_aero, τ_ssa_aero, τ_ssag_aero =
@@ -85,7 +92,8 @@ to the exising `TwoStream` optics properties:
                 end
                 τ[glay], ssa[glay], g[glay] = increment_2stream(τ[glay], ssa[glay], g[glay], τ_aero, ssa_aero, g_aero)
                 if collect_aod
-                    aod[] += τ_aero
+                    aod_ext[] += τ_aero
+                    aod_sca[] += τ_ssa_aero
                 end
             end
         end

test/all_sky_with_aerosols_utils.jl

@@ -176,7 +176,8 @@ function all_sky_with_aerosols(
     println("L∞ error in flux_up           = $max_err_flux_up_sw")
     println("L∞ error in flux_dn           = $max_err_flux_dn_sw")
     println("L∞ error in flux_net          = $max_err_flux_net_sw")
-    @info as.aerosol_state.aod
+    @info as.aerosol_state.aod_sw_ext
+    @info as.aerosol_state.aod_sw_sca
     println("L∞ relative error in flux_net = $(max_rel_err_flux_net_sw * 100) %\n")
 
     # The reference results for the longwave solver are generated using a non-scattering solver,

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