Skip to content

Commit

Permalink
Split LookUp table and Pade cloud optics tables.
Browse files Browse the repository at this point in the history
  • Loading branch information
@sriharshakandala
sriharshakandala committed Dec 4, 2023
1 parent 40486d9 commit b9a063a
Show file tree
Hide file tree
Showing 6 changed files with 184 additions and 151 deletions.
162 changes: 89 additions & 73 deletions src/optics/CloudOptics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,7 @@ end
This function computed the TwoSteam cloud optics properties using either the
lookup table method or pade method.
"""
function compute_cld_props(lkp_cld, as, cld_mask, glay, gcol, ibnd, igpt)
use_lut = lkp_cld.use_lut
function compute_cld_props(lkp_cld::LookUpCld, as, cld_mask, glay, gcol, ibnd, igpt)
re_liq = as.cld_r_eff_liq[glay, gcol]
re_ice = as.cld_r_eff_ice[glay, gcol]
ice_rgh = as.ice_rgh
Expand All @@ -75,83 +74,99 @@ function compute_cld_props(lkp_cld, as, cld_mask, glay, gcol, ibnd, igpt)
FT = eltype(re_liq)
τl, τl_ssa, τl_ssag = FT(0), FT(0), FT(0)
τi, τi_ssa, τi_ssag = FT(0), FT(0), FT(0)
if use_lut # use LUT interpolation
(;
lut_extliq,
lut_ssaliq,
lut_asyliq,
lut_extice,
lut_ssaice,
lut_asyice,
radliq_lwr,
radliq_upr,
radice_lwr,
radice_upr,
nsize_liq,
nsize_ice,
) = lkp_cld
Δr_liq = (radliq_upr - radliq_lwr) / FT(nsize_liq - 1)
Δr_ice = (radice_upr - radice_lwr) / FT(nsize_ice - 1)
# cloud liquid particles
if cld_path_liq > eps(FT)
loc = Int(max(min(unsafe_trunc(Int, (re_liq - radliq_lwr) / Δr_liq) + 1, nsize_liq - 1), 1))
fac = (re_liq - radliq_lwr - (loc - 1) * Δr_liq) / Δr_liq
fc1 = FT(1) - fac
τl = (fc1 * lut_extliq[loc, ibnd] + fac * lut_extliq[loc + 1, ibnd]) * cld_path_liq
τl_ssa = (fc1 * lut_ssaliq[loc, ibnd] + fac * lut_ssaliq[loc + 1, ibnd]) * τl
τl_ssag = (fc1 * lut_asyliq[loc, ibnd] + fac * lut_asyliq[loc + 1, ibnd]) * τl_ssa
end
# cloud ice particles
if cld_path_ice > eps(FT)
loc = Int(max(min(unsafe_trunc(Int, (re_ice - radice_lwr) / Δr_ice) + 1, nsize_ice - 1), 1))
fac = (re_ice - radice_lwr - (loc - 1) * Δr_ice) / Δr_ice
fc1 = FT(1) - fac
τi = (fc1 * lut_extice[loc, ibnd, ice_rgh] + fac * lut_extice[loc + 1, ibnd, ice_rgh]) * cld_path_ice
τi_ssa = (fc1 * lut_ssaice[loc, ibnd, ice_rgh] + fac * lut_ssaice[loc + 1, ibnd, ice_rgh]) * τi
τi_ssag = (fc1 * lut_asyice[loc, ibnd, ice_rgh] + fac * lut_asyice[loc + 1, ibnd, ice_rgh]) * τi_ssa
end
else # use pade interpolation
(;
pade_extliq,
pade_ssaliq,
pade_asyliq,
pade_extice,
pade_ssaice,
pade_asyice,
pade_sizreg_extliq,
pade_sizreg_ssaliq,
pade_sizreg_asyliq,
pade_sizreg_extice,
pade_sizreg_ssaice,
pade_sizreg_asyice,
) = lkp_cld
m_ext, m_ssa_g = 3, 3
n_ext, n_ssa_g = 3, 2
# Finds index into size regime table
# This works only if there are precisely three size regimes (four bounds) and it's
# previously guaranteed that size_bounds(1) <= size <= size_bounds(4)
if cld_path_liq > eps(FT)
irad = Int(min(floor((re_liq - pade_sizreg_extliq[2]) / pade_sizreg_extliq[3]) + 2, 3))
τl = pade_eval(ibnd, re_liq, irad, m_ext, n_ext, pade_extliq) * cld_path_liq
(;
lut_extliq,
lut_ssaliq,
lut_asyliq,
lut_extice,
lut_ssaice,
lut_asyice,
radliq_lwr,
radliq_upr,
radice_lwr,
radice_upr,
nsize_liq,
nsize_ice,
) = lkp_cld
Δr_liq = (radliq_upr - radliq_lwr) / FT(nsize_liq - 1)
Δr_ice = (radice_upr - radice_lwr) / FT(nsize_ice - 1)
# cloud liquid particles
if cld_path_liq > eps(FT)
loc = Int(max(min(unsafe_trunc(Int, (re_liq - radliq_lwr) / Δr_liq) + 1, nsize_liq - 1), 1))
fac = (re_liq - radliq_lwr - (loc - 1) * Δr_liq) / Δr_liq
fc1 = FT(1) - fac
τl = (fc1 * lut_extliq[loc, ibnd] + fac * lut_extliq[loc + 1, ibnd]) * cld_path_liq
τl_ssa = (fc1 * lut_ssaliq[loc, ibnd] + fac * lut_ssaliq[loc + 1, ibnd]) * τl
τl_ssag = (fc1 * lut_asyliq[loc, ibnd] + fac * lut_asyliq[loc + 1, ibnd]) * τl_ssa
end
# cloud ice particles
if cld_path_ice > eps(FT)
loc = Int(max(min(unsafe_trunc(Int, (re_ice - radice_lwr) / Δr_ice) + 1, nsize_ice - 1), 1))
fac = (re_ice - radice_lwr - (loc - 1) * Δr_ice) / Δr_ice
fc1 = FT(1) - fac
τi = (fc1 * lut_extice[loc, ibnd, ice_rgh] + fac * lut_extice[loc + 1, ibnd, ice_rgh]) * cld_path_ice
τi_ssa = (fc1 * lut_ssaice[loc, ibnd, ice_rgh] + fac * lut_ssaice[loc + 1, ibnd, ice_rgh]) * τi
τi_ssag = (fc1 * lut_asyice[loc, ibnd, ice_rgh] + fac * lut_asyice[loc + 1, ibnd, ice_rgh]) * τi_ssa
end

irad = Int(min(floor((re_liq - pade_sizreg_ssaliq[2]) / pade_sizreg_ssaliq[3]) + 2, 3))
τl_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_ssaliq))) * τl
τ = τl + τi
τ_ssa = τl_ssa + τi_ssa
τ_ssag = (τl_ssag + τi_ssag) / max(eps(FT), τ_ssa)
τ_ssa /= max(eps(FT), τ)

irad = Int(min(floor((re_liq - pade_sizreg_asyliq[2]) / pade_sizreg_asyliq[3]) + 2, 3))
τl_ssag = pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_asyliq) * τl_ssa
end
return (τ, τ_ssa, τ_ssag)
end

function compute_cld_props(lkp_cld::PadeCld, as, cld_mask, glay, gcol, ibnd, igpt)
re_liq = as.cld_r_eff_liq[glay, gcol]
re_ice = as.cld_r_eff_ice[glay, gcol]
ice_rgh = as.ice_rgh
cld_path_liq = as.cld_path_liq[glay, gcol]
cld_path_ice = as.cld_path_ice[glay, gcol]
FT = eltype(re_liq)
τl, τl_ssa, τl_ssag = FT(0), FT(0), FT(0)
τi, τi_ssa, τi_ssag = FT(0), FT(0), FT(0)

if cld_path_ice > eps(FT)
irad = Int(min(floor((re_ice - pade_sizreg_extice[2]) / pade_sizreg_extice[3]) + 2, 3))
(;
pade_extliq,
pade_ssaliq,
pade_asyliq,
pade_extice,
pade_ssaice,
pade_asyice,
pade_sizreg_extliq,
pade_sizreg_ssaliq,
pade_sizreg_asyliq,
pade_sizreg_extice,
pade_sizreg_ssaice,
pade_sizreg_asyice,
) = lkp_cld
m_ext, m_ssa_g = 3, 3
n_ext, n_ssa_g = 3, 2
# Finds index into size regime table
# This works only if there are precisely three size regimes (four bounds) and it's
# previously guaranteed that size_bounds(1) <= size <= size_bounds(4)
if cld_path_liq > eps(FT)
irad = Int(min(floor((re_liq - pade_sizreg_extliq[2]) / pade_sizreg_extliq[3]) + 2, 3))
τl = pade_eval(ibnd, re_liq, irad, m_ext, n_ext, pade_extliq) * cld_path_liq

τi = pade_eval(ibnd, re_ice, irad, m_ext, n_ext, pade_extice, ice_rgh) * cld_path_ice
irad = Int(min(floor((re_liq - pade_sizreg_ssaliq[2]) / pade_sizreg_ssaliq[3]) + 2, 3))
τl_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_ssaliq))) * τl

irad = Int(min(floor((re_ice - pade_sizreg_ssaice[2]) / pade_sizreg_ssaice[3]) + 2, 3))
τi_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_ssaice, ice_rgh))) * τi
irad = Int(min(floor((re_liq - pade_sizreg_asyliq[2]) / pade_sizreg_asyliq[3]) + 2, 3))
τl_ssag = pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_asyliq) * τl_ssa
end

irad = Int(min(floor((re_ice - pade_sizreg_asyice[2]) / pade_sizreg_asyice[3]) + 2, 3))
τi_ssag = pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_asyice, ice_rgh) * τi_ssa
end
if cld_path_ice > eps(FT)
irad = Int(min(floor((re_ice - pade_sizreg_extice[2]) / pade_sizreg_extice[3]) + 2, 3))

τi = pade_eval(ibnd, re_ice, irad, m_ext, n_ext, pade_extice, ice_rgh) * cld_path_ice

irad = Int(min(floor((re_ice - pade_sizreg_ssaice[2]) / pade_sizreg_ssaice[3]) + 2, 3))
τi_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_ssaice, ice_rgh))) * τi

irad = Int(min(floor((re_ice - pade_sizreg_asyice[2]) / pade_sizreg_asyice[3]) + 2, 3))
τi_ssag = pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_asyice, ice_rgh) * τi_ssa
end

τ = τl + τi
Expand All @@ -162,6 +177,7 @@ function compute_cld_props(lkp_cld, as, cld_mask, glay, gcol, ibnd, igpt)
return (τ, τ_ssa, τ_ssag)
end


"""
pade_eval(
ibnd,
Expand Down
131 changes: 74 additions & 57 deletions src/optics/LookUpTables.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ module LookUpTables
using DocStringExtensions
using Adapt

export AbstractLookUp, LookUpLW, LookUpSW, LookUpCld
export AbstractLookUp, LookUpLW, LookUpSW, LookUpCld, PadeCld

"""
AbstractLookUp{FT}
Expand Down Expand Up @@ -798,7 +798,7 @@ These are used to determine the optical properties of ice and water cloud togeth
# Fields
$(DocStringExtensions.FIELDS)
"""
struct LookUpCld{B, FT, FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp{FT}
struct LookUpCld{FT, FTA1D, FTA2D, FTA3D} <: AbstractLookUp{FT}
"number of bands"
nband::Int
"number of ice roughness types"
Expand All @@ -809,12 +809,6 @@ struct LookUpCld{B, FT, FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp{FT}
nsize_ice::Int
"number of size regimes"
nsizereg::Int
"number of extinction coefficients for pade approximation"
ncoeff_ext::Int
"number of ssa/g coefficients for pade approximation"
ncoeff_ssa_g::Int
"number of size regime boundaries for pade interpolation"
nbound::Int
"pair = 2"
pair::Int
"liquid particle size lower bound for LUT interpolation"
Expand All @@ -841,43 +835,16 @@ struct LookUpCld{B, FT, FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp{FT}
lut_ssaice::FTA3D
"LUT ice asymmetry parameter (`nsize_ice, nband, nrghice`)"
lut_asyice::FTA3D
"pade coefficients for liquid extinction (`nband, nsizereg, ncoeff_ext`)"
pade_extliq::FTA3D
"pade coefficients for liquid single scattening albedo (`nband, nsizereg, ncoeff_ssa_g`)"
pade_ssaliq::FTA3D
"pade coefficients for liquid asymmetry paramter (`nband, nsizereg, ncoeff_ssa_g`)"
pade_asyliq::FTA3D
"pade coefficients for ice extinction (`nband, nsizereg, ncoeff_ext, nrghice`)"
pade_extice::FTA4D
"pade coefficients for ice single scattening albedo (`nband, nsizereg, ncoeff_ssa_g, nrghice`)"
pade_ssaice::FTA4D
"pade coefficients for ice asymmetry paramter (`nband, nsizereg, ncoeff_ssa_g, nrghice`)"
pade_asyice::FTA4D
"pade size regime boundaries for liquid extinction coefficient for pade interpolation (`nbound`) μm"
pade_sizreg_extliq::FTA1D
"pade size regime boundaries for liquid single scattering albedo for pade interpolation (`nbound`) μm"
pade_sizreg_ssaliq::FTA1D
"pade size regime boundaries for liquid asymmetry parameter for pade interpolation (`nbound`) μm"
pade_sizreg_asyliq::FTA1D
"pade size regime boundaries for ice extinction coefficient for pade interpolation (`nbound`) μm"
pade_sizreg_extice::FTA1D
"pade size regime boundaries for ice single scattering albedo for pade interpolation (`nbound`) μm"
pade_sizreg_ssaice::FTA1D
"pade size regime boundaries for ice asymmetry parameter for pade interpolation (`nbound`) μm"
pade_sizreg_asyice::FTA1D
"beginning and ending wavenumber for each band (`2, nband`) cm⁻¹"
bnd_lims_wn::FTA2D
"use LUT (default) or pade interpolation"
use_lut::B
end
Adapt.@adapt_structure LookUpCld

function LookUpCld(ds, ::Type{FT}, ::Type{DA}, use_lut::Bool = true) where {FT <: AbstractFloat, DA}
function LookUpCld(ds, ::Type{FT}, ::Type{DA}) where {FT <: AbstractFloat, DA}

FTA1D = DA{FT, 1}
FTA2D = DA{FT, 2}
FTA3D = DA{FT, 3}
FTA4D = DA{FT, 4}

nband = Int(ds.dim["nband"])
nrghice = Int(ds.dim["nrghice"])
Expand Down Expand Up @@ -907,6 +874,76 @@ function LookUpCld(ds, ::Type{FT}, ::Type{DA}, use_lut::Bool = true) where {FT <
lut_ssaice = FTA3D(Array(ds["lut_ssaice"]))
lut_asyice = FTA3D(Array(ds["lut_asyice"]))

bnd_lims_wn = FTA2D(Array(ds["bnd_limits_wavenumber"]))

return (LookUpCld{FT, FTA1D, FTA2D, FTA3D}(
nband,
nrghice,
nsize_liq,
nsize_ice,
nsizereg,
pair,
radliq_lwr,
radliq_upr,
radliq_fac,
radice_lwr,
radice_upr,
radice_fac,
lut_extliq,
lut_ssaliq,
lut_asyliq,
lut_extice,
lut_ssaice,
lut_asyice,
bnd_lims_wn,
))
end

struct PadeCld{FT, FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp{FT}
"number of extinction coefficients for pade approximation"
ncoeff_ext::Int
"number of ssa/g coefficients for pade approximation"
ncoeff_ssa_g::Int
"number of size regime boundaries for pade interpolation"
nbound::Int
"pade coefficients for liquid extinction (`nband, nsizereg, ncoeff_ext`)"
pade_extliq::FTA3D
"pade coefficients for liquid single scattening albedo (`nband, nsizereg, ncoeff_ssa_g`)"
pade_ssaliq::FTA3D
"pade coefficients for liquid asymmetry paramter (`nband, nsizereg, ncoeff_ssa_g`)"
pade_asyliq::FTA3D
"pade coefficients for ice extinction (`nband, nsizereg, ncoeff_ext, nrghice`)"
pade_extice::FTA4D
"pade coefficients for ice single scattening albedo (`nband, nsizereg, ncoeff_ssa_g, nrghice`)"
pade_ssaice::FTA4D
"pade coefficients for ice asymmetry paramter (`nband, nsizereg, ncoeff_ssa_g, nrghice`)"
pade_asyice::FTA4D
"pade size regime boundaries for liquid extinction coefficient for pade interpolation (`nbound`) μm"
pade_sizreg_extliq::FTA1D
"pade size regime boundaries for liquid single scattering albedo for pade interpolation (`nbound`) μm"
pade_sizreg_ssaliq::FTA1D
"pade size regime boundaries for liquid asymmetry parameter for pade interpolation (`nbound`) μm"
pade_sizreg_asyliq::FTA1D
"pade size regime boundaries for ice extinction coefficient for pade interpolation (`nbound`) μm"
pade_sizreg_extice::FTA1D
"pade size regime boundaries for ice single scattering albedo for pade interpolation (`nbound`) μm"
pade_sizreg_ssaice::FTA1D
"pade size regime boundaries for ice asymmetry parameter for pade interpolation (`nbound`) μm"
pade_sizreg_asyice::FTA1D
"beginning and ending wavenumber for each band (`2, nband`) cm⁻¹"
bnd_lims_wn::FTA2D
end

function PadeCld(ds, ::Type{FT}, ::Type{DA}) where {FT <: AbstractFloat, DA}
FTA1D = DA{FT, 1}
FTA2D = DA{FT, 2}
FTA3D = DA{FT, 3}
FTA4D = DA{FT, 4}

ncoeff_ext = Int(ds.dim["ncoeff_ext"])
ncoeff_ssa_g = Int(ds.dim["ncoeff_ssa_g"])
nbound = Int(ds.dim["nbound"])

pade_extliq = FTA3D(Array(ds["pade_extliq"]))
pade_ssaliq = FTA3D(Array(ds["pade_ssaliq"]))
pade_asyliq = FTA3D(Array(ds["pade_asyliq"]))
Expand All @@ -925,28 +962,10 @@ function LookUpCld(ds, ::Type{FT}, ::Type{DA}, use_lut::Bool = true) where {FT <

bnd_lims_wn = FTA2D(Array(ds["bnd_limits_wavenumber"]))

return (LookUpCld{Bool, FT, FTA1D, FTA2D, FTA3D, FTA4D}(
nband,
nrghice,
nsize_liq,
nsize_ice,
nsizereg,
return (PadeCld{FT, FTA1D, FTA2D, FTA3D, FTA4D}(
ncoeff_ext,
ncoeff_ssa_g,
nbound,
pair,
radliq_lwr,
radliq_upr,
radliq_fac,
radice_lwr,
radice_upr,
radice_fac,
lut_extliq,
lut_ssaliq,
lut_asyliq,
lut_extice,
lut_ssaice,
lut_asyice,
pade_extliq,
pade_ssaliq,
pade_asyliq,
Expand All @@ -960,9 +979,7 @@ function LookUpCld(ds, ::Type{FT}, ::Type{DA}, use_lut::Bool = true) where {FT <
pade_sizreg_ssaice,
pade_sizreg_asyice,
bnd_lims_wn,
use_lut,
))

end

end
4 changes: 2 additions & 2 deletions src/optics/Optics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ end
gcol::Int,
igpt::Int,
lkp::LookUpLW{FT},
lkp_cld::Union{LookUpCld,Nothing} = nothing,
lkp_cld::Union{LookUpCld,PadeCld,Nothing} = nothing,
) where {FT<:AbstractFloat}
Computes optical properties for the longwave problem.
Expand Down Expand Up @@ -172,7 +172,7 @@ end
gcol::Int,
igpt::Int,
lkp::LookUpSW{FT},
lkp_cld::Union{LookUpCld,Nothing} = nothing,
lkp_cld::Union{LookUpCld,PadeCld,Nothing} = nothing,
) where {FT<:AbstractFloat}
Computes optical properties for the shortwave problem.
Expand Down
Loading

0 comments on commit b9a063a

Please sign in to comment.