diff --git a/NEWS.md b/NEWS.md index 1ad3a87b3..c60671606 100644 --- a/NEWS.md +++ b/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 diff --git a/Project.toml b/Project.toml index c6e9ed255..6d860e392 100644 --- a/Project.toml +++ b/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" diff --git a/src/optics/AtmosphericStates.jl b/src/optics/AtmosphericStates.jl index 3bbf096ee..cd90353e2 100644 --- a/src/optics/AtmosphericStates.jl +++ b/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)" diff --git a/src/optics/LookUpTables.jl b/src/optics/LookUpTables.jl index 6fd4c2f4b..46801c490 100644 --- a/src/optics/LookUpTables.jl +++ b/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 diff --git a/src/optics/Optics.jl b/src/optics/Optics.jl index 599998f94..49f176412 100644 --- a/src/optics/Optics.jl +++ b/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 diff --git a/src/optics/aerosol_optics.jl b/src/optics/aerosol_optics.jl index 0dc48b9ed..1302780e0 100644 --- a/src/optics/aerosol_optics.jl +++ b/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 diff --git a/test/all_sky_with_aerosols_utils.jl b/test/all_sky_with_aerosols_utils.jl index 697eee94d..54457f593 100644 --- a/test/all_sky_with_aerosols_utils.jl +++ b/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 diff --git a/test/read_all_sky_with_aerosols.jl b/test/read_all_sky_with_aerosols.jl index 7d5707614..e48a6883a 100644 --- a/test/read_all_sky_with_aerosols.jl +++ b/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]),