Skip to content

Commit

Permalink
Use ifelse in liquid_fraction, plus a comment
Browse files Browse the repository at this point in the history
  • Loading branch information
glwagner committed Jan 18, 2024
1 parent d73a834 commit 17c8f2d
Showing 1 changed file with 21 additions and 20 deletions.
41 changes: 21 additions & 20 deletions src/relations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -968,7 +968,7 @@ end

# we may be hitting a slow path:
# https://stackoverflow.com/questions/14687665/very-slow-stdpow-for-bases-very-close-to-1
pow_hack(x, y) = exp(y * log(x))
@inline pow_hack(x, y) = exp(y * log(x))

function saturation_vapor_pressure(
param_set::APS,
Expand Down Expand Up @@ -1290,24 +1290,26 @@ is a function that is 1 above `T_freeze` and goes to zero below `T_icenuc`.
function liquid_fraction(
param_set::APS,
T::FT,
::Type{phase_type},
::Type{State},
q::PhasePartition{FT} = q_pt_0(FT),
) where {FT <: Real, phase_type <: ThermodynamicState}
_T_freeze::FT = TP.T_freeze(param_set)
_T_icenuc::FT = TP.T_icenuc(param_set)

if T > _T_freeze
return FT(1)
elseif T > _T_icenuc
_pow_icenuc::FT = TP.pow_icenuc(param_set)
if _pow_icenuc == 1
return (T - _T_icenuc) / (_T_freeze - _T_icenuc)
else
return ((T - _T_icenuc) / (_T_freeze - _T_icenuc))^_pow_icenuc
end
else
return FT(0)
end
) where State <: ThermodynamicState

Tᶠ = TP.T_freeze(param_set) # freezing temperature
Tⁿ = TP.T_icenuc(param_set) # temperature of homogeneous ice nucleation
α = TP.pow_icenuc(param_set) # power law partial ice nucleation parameter

# Model for cloud water condensate probabilty when temperature is below
# freezing (all liquid condensate), but above the temperature of homogeneous
# ice nucleation (all ice condensate). In it's simplest form, this is simply
# a number that varies between 0 and 1. For example, see figure 6 in
# Hu et al., https://doi.org/10.1029/2009JD012384, JGR (2010).
λᵖ = ((T - Tⁿ) / (Tᶠ - Tⁿ))^α

above_freezing = T > Tᶠ
supercooled_liquid = (T Tᶠ) & (T > Tⁿ)

return ifelse(above_freezing, one(FT),
ifelse(supercooled_liquid, λᵖ, zero(FT)))
end

function liquid_fraction(
Expand Down Expand Up @@ -2776,8 +2778,7 @@ relative_humidity(
PhasePartition(param_set, ts),
)

relative_humidity(param_set::APS, ts::AbstractPhaseDry{FT}) where {FT <: Real} =
FT(0)
@inline relative_humidity(::APS, ::AbstractPhaseDry{FT}) where FT = zero(FT)

"""
total_specific_enthalpy(e_tot, R_m, T)
Expand Down

0 comments on commit 17c8f2d

Please sign in to comment.