Skip to content

Commit

Permalink
Add type restrictions for debugging
Browse files Browse the repository at this point in the history
  • Loading branch information
charleskawczynski committed Nov 22, 2024
1 parent 7e00bf3 commit 4b49acb
Showing 1 changed file with 40 additions and 34 deletions.
74 changes: 40 additions & 34 deletions src/parameterized_tendencies/radiation/held_suarez.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
##### Held-Suarez
#####

import Thermodynamics as TD
import Thermodynamics.Parameters as TDP
import ClimaCore.Spaces as Spaces
import ClimaCore.Fields as Fields

Expand Down Expand Up @@ -32,8 +34,6 @@ function held_suarez_ΔT_y_T_equator(
return ΔT_y, T_equator
end

height_factor(σ, σ_b) = max(0, (σ - σ_b) / (1 - σ_b))

struct HeldSuarezForcingParams{FT}
day::FT
σ_b::FT
Expand All @@ -49,40 +49,39 @@ end
Base.Broadcast.broadcastable(x::HeldSuarezForcingParams) = tuple(x)

function compute_ΔρT(
thermo_params,
ts_surf,
ρ,
p,
lat,
z_surface,
thermo_params::TDP.ThermodynamicsParameters,
ts_surf::TD.ThermodynamicState,
ρ::FT,
p::FT,
lat::FT,
z_surface::FT,
s::HeldSuarezForcingParams,
)
) where {FT}
σ = compute_σ(thermo_params, z_surface, p, ts_surf, s)
k_a = 1 / (40 * s.day)
k_s = 1 / (4 * s.day)

φ = deg2rad(lat)
return
(k_a + (k_s - k_a) * height_factor(σ, s.σ_b) * abs2(abs2(cos(φ)))) *
ρ *
( # ᶜT - ᶜT_equil
p /* s.R_d) - max(
s.T_min,
(
s.T_equator - ΔT_y * abs2(sin(φ)) -
s.Δθ_z * log(p / s.p_ref_theta) * abs2(cos(φ))
) * fast_pow(p / s.p_ref_theta, s.κ_d),
)
)
return (k_a + (k_s - k_a) * height_factor(σ, s.σ_b) * abs2(abs2(cos(φ)))) *
ρ *
( # ᶜT - ᶜT_equil
p /* s.R_d) - max(
s.T_min,
(
s.T_equator - ΔT_y * abs2(sin(φ)) -
s.Δθ_z * log(p / s.p_ref_theta) * abs2(cos(φ))
) * fast_pow(p / s.p_ref_theta, s.κ_d),
)
)
end

function compute_σ(
thermo_params,
z_surface,
p,
ts_surf,
thermo_params::TDP.ThermodynamicsParameters,
z_surface::FT,
p::FT,
ts_surf::TD.ThermodynamicState,
s::HeldSuarezForcingParams,
)
) where {FT}
p / (
s.MSLP * exp(
-s.grav * z_surface / s.R_d /
Expand All @@ -91,6 +90,16 @@ function compute_σ(
)
end

height_factor::FT, σ_b::FT) where {FT} = max(0, (σ - σ_b) / (1 - σ_b))
height_factor(
thermo_params::TDP.ThermodynamicsParameters,
z_surface::FT,
p::FT,
ts_surf::TD.ThermodynamicState,
s::HeldSuarezForcingParams,
) where {FT} =
height_factor(compute_σ(thermo_params, z_surface, p, ts_surf, s), s.σ_b)

function forcing_tendency!(Yₜ, Y, p, t, ::HeldSuarezForcing)
(; params) = p
(; ᶜp, sfc_conditions) = p.precomputed
Expand Down Expand Up @@ -131,14 +140,11 @@ function forcing_tendency!(Yₜ, Y, p, t, ::HeldSuarezForcing)
@. Yₜ.c.uₕ -=
(
k_f * height_factor(
compute_σ(
thermo_params,
z_surface,
ᶜp,
sfc_conditions.ts,
hs_params,
),
σ_b,
thermo_params,
z_surface,
ᶜp,
sfc_conditions.ts,
hs_params,
)
) * Y.c.uₕ
@. Yₜ.c.ρe_tot -=
Expand Down

0 comments on commit 4b49acb

Please sign in to comment.