Skip to content

Commit

Permalink
Split LookUp table and Pade cloud optics tables. (#408)
Browse files Browse the repository at this point in the history
  • Loading branch information
@sriharshakandala
sriharshakandala authored Dec 5, 2023
1 parent 40486d9 commit e53af28
Show file tree
Hide file tree
Showing 8 changed files with 222 additions and 204 deletions.
1 change: 1 addition & 0 deletions docs/src/optics/LookUpTables.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,5 +8,6 @@ CurrentModule = RRTMGP.LookUpTables
LookUpLW
LookUpSW
LookUpCld
PadeCld
AbstractLookUp
```
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,
) = lkp_cld
nsize_liq = LookUpTables.get_nsize_liq(lkp_cld)
nsize_ice = LookUpTables.get_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
37 changes: 10 additions & 27 deletions src/optics/GasOptics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ function compute_col_gas_kernel!(
end
"""
compute_interp_fractions(
lkp::AbstractLookUp{FT},
lkp::AbstractLookUp,
vmr,
p_lay,
t_lay,
Expand All @@ -50,16 +50,7 @@ end
compute interpolation fractions for binary species parameter, pressure and temperature.
"""
@inline function compute_interp_fractions(
lkp::AbstractLookUp{FT},
vmr,
p_lay,
t_lay,
tropo,
ibnd,
glay,
gcol,
) where {FT <: AbstractFloat}
@inline function compute_interp_fractions(lkp::AbstractLookUp, vmr, p_lay, t_lay, tropo, ibnd, glay, gcol)
jftemp = compute_interp_frac_temp(lkp, t_lay, glay, gcol)
jfpress = compute_interp_frac_press(lkp, p_lay, tropo, glay, gcol)
jfη, col_mix = compute_interp_frac_η(lkp, vmr, tropo, jftemp[1], ibnd, glay, gcol)
Expand All @@ -68,15 +59,15 @@ end

"""
compute_interp_frac_temp(
lkp::AbstractLookUp{FT},
lkp::AbstractLookUp,
t_lay,
glay,
gcol,
) where {FT<:AbstractFloat}
compute interpolation fraction for temperature.
"""
@inline function compute_interp_frac_temp(lkp::AbstractLookUp{FT}, t_lay, glay, gcol) where {FT <: AbstractFloat}
@inline function compute_interp_frac_temp(lkp::AbstractLookUp, t_lay, glay, gcol)
(; Δ_t_ref, n_t_ref, t_ref) = lkp

@inbounds jtemp = loc_lower(t_lay, Δ_t_ref, n_t_ref, t_ref)
Expand Down Expand Up @@ -110,32 +101,24 @@ end

"""
compute_interp_frac_η(
lkp::AbstractLookUp{FT},
lkp::AbstractLookUp,
vmr,
tropo,
jtemp,
ibnd,
glay,
gcol,
) where {FT<:AbstractFloat}
)
Compute interpolation fraction for binary species parameter.
"""
@inline function compute_interp_frac_η(
lkp::AbstractLookUp{FT},
vmr,
tropo,
jtemp,
ibnd,
glay,
gcol,
) where {FT <: AbstractFloat}
@inline function compute_interp_frac_η(lkp::AbstractLookUp, vmr, tropo, jtemp, ibnd, glay, gcol)
(; n_η, key_species, vmr_ref) = lkp
ig = view(key_species, :, tropo, ibnd)

vmr1 = get_vmr(vmr, ig[1], glay, gcol)
vmr2 = get_vmr(vmr, ig[2], glay, gcol)

FT = eltype(vmr1)
itemp = 1
@inbounds η_half = vmr_ref[tropo, ig[1] + 1, jtemp + itemp - 1] / vmr_ref[tropo, ig[2] + 1, jtemp + itemp - 1]
col_mix1 = vmr1 + η_half * vmr2
Expand All @@ -161,7 +144,7 @@ end

"""
compute_τ_ssa_lw_src!(
lkp::AbstractLookUp{FT},
lkp::AbstractLookUp,
vmr,
col_dry,
igpt,
Expand All @@ -176,7 +159,7 @@ Compute optical thickness, single scattering albedo, asymmetry parameter
and longwave sources whenever applicable.
"""
@inline function compute_τ_ssa_lw_src!(
lkp::AbstractLookUp{FT},
lkp::AbstractLookUp,
vmr,
col_dry,
igpt,
Expand Down
Loading

0 comments on commit e53af28

Please sign in to comment.