Skip to content

Commit ab14894

Browse files
haakon-ehaakon-e
committed
remove dt limiting from p3 melt, nuc
1 parent fbd7881 commit ab14894

File tree

4 files changed

+88
-99
lines changed

4 files changed

+88
-99
lines changed

‎docs/src/plots/P3ImmersionFreezing.jl‎

Lines changed: 12 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,6 @@ p = FT(800 * 1e2)
3030
dd = CMP.DesertDust(FT)
3131
il = CMP.Illite(FT)
3232

33-
# model time step (for limiting)
34-
dt = FT(1)
35-
3633
# plot data
3734
RH_range = range(0.8, stop = 1.2, length = 1000)
3835
T1 = FT(273.15 - 15)
@@ -41,22 +38,23 @@ T2 = FT(273.15 - 35)
4138
#! format: off
4239

4340
# limiters to not nucleate more mass and number than we have in liquid phase
44-
max_dLdt_T1 = [qₗ* p2ρ(T1, RH) / dt for RH in RH_range]
45-
max_dLdt_T2 = [qₗ* p2ρ(T2, RH) / dt for RH in RH_range]
46-
max_dNdt = [Nₗ / dt for RH in RH_range]
41+
# (implicitly dt = 1)
42+
max_dLdt_T1 = @. qₗ * p2ρ(T1, RH_range)
43+
max_dLdt_T2 = @. qₗ * p2ρ(T2, RH_range)
44+
max_dNdt = fill(Nₗ, size(RH_range))
4745

48-
dLdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dLdt for RH in RH_range]
49-
dNdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dNdt for RH in RH_range]
46+
dLdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dLdt for RH in RH_range]
47+
dNdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dNdt for RH in RH_range]
5048

51-
dLdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dLdt for RH in RH_range]
52-
dNdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dNdt for RH in RH_range]
49+
dLdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dLdt for RH in RH_range]
50+
dNdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dNdt for RH in RH_range]
5351

5452

55-
dLdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dLdt for RH in RH_range]
56-
dNdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dNdt for RH in RH_range]
53+
dLdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dLdt for RH in RH_range]
54+
dNdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dNdt for RH in RH_range]
5755

58-
dLdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dLdt for RH in RH_range]
59-
dNdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dNdt for RH in RH_range]
56+
dLdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dLdt for RH in RH_range]
57+
dNdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dNdt for RH in RH_range]
6058

6159
# plotting
6260
fig = PL.Figure(size = (1500, 500), fontsize=22, linewidth=3)

‎docs/src/plots/P3Melting.jl‎

Lines changed: 58 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -23,73 +23,73 @@ Nᵢ = FT(2e5)
2323
ΔT_range = range(1e-4, stop = 0.025, length = 1000)
2424

2525
# limiters to not melt more mass and number than we have
26-
max_dLdt = [Lᵢ / dt for ΔT in ΔT_range]
27-
max_dNdt = [Nᵢ / dt for ΔT in ΔT_range]
26+
max_dLdt = Lᵢ / dt
27+
max_dNdt = Nᵢ / dt
2828

2929
#! format: off
3030

31-
ρₐ1 = FT(1.2)
32-
Fᵣ1 = FT(0.8)
33-
ρᵣ1 = FT(800)
34-
state = P3.get_state(params; F_rim = Fᵣ1, ρ_rim = ρᵣ1, L_ice = Lᵢ, N_ice = Nᵢ)
31+
ρₐ = FT(1.2)
32+
Fᵣ = FT(0.8)
33+
ρᵣ = FT(800)
34+
label1 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
35+
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
3536
logλ = P3.get_distribution_logλ(state)
36-
dLdt1 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
37-
dNdt1 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
37+
melt1 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
38+
dLdt1 = getfield.(melt1, :dLdt)
39+
dNdt1 = getfield.(melt1, :dNdt)
3840

39-
Fᵣ2 = FT(0.2)
40-
state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ1, L_ice = Lᵢ, N_ice = Nᵢ)
41+
Fᵣ = FT(0.2)
42+
label2 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
43+
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
4144
logλ = P3.get_distribution_logλ(state)
42-
dLdt2 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
43-
dNdt2 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
45+
melt2 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
46+
dLdt2 = getfield.(melt2, :dLdt)
47+
dNdt2 = getfield.(melt2, :dNdt)
4448

45-
ρᵣ2 = FT(200)
46-
state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ2, L_ice = Lᵢ, N_ice = Nᵢ)
49+
ρᵣ = FT(200)
50+
label3 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
51+
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
4752
logλ = P3.get_distribution_logλ(state)
48-
dLdt3 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
49-
dNdt3 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
53+
melt3 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
54+
dLdt3 = getfield.(melt3, :dLdt)
55+
dNdt3 = getfield.(melt3, :dNdt)
5056

51-
ρₐ2 = FT(0.5)
52-
state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ2, L_ice = Lᵢ, N_ice = Nᵢ)
57+
ρₐ = FT(0.5)
58+
label4 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
59+
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
5360
logλ = P3.get_distribution_logλ(state)
54-
dLdt4 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt, state, logλ).dLdt for ΔT in ΔT_range]
55-
dNdt4 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt, state, logλ).dNdt for ΔT in ΔT_range]
61+
melt4 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
62+
dLdt4 = getfield.(melt4, :dLdt)
63+
dNdt4 = getfield.(melt4, :dNdt)
5664

5765
# plotting
58-
fig = Makie.Figure(size = (1500, 500), fontsize=22, linewidth=3)
59-
60-
ax1 = Makie.Axis(fig[1, 1]; yscale = log10)
61-
ax2 = Makie.Axis(fig[1, 2]; yscale = log10)
62-
63-
ax1.xlabel = "T [C]"
64-
ax1.ylabel = "ice mass melting rate [g/m3/s]"
65-
ax2.xlabel = "T [C]"
66-
ax2.ylabel = "ice number melting rate [1/cm3/s]"
67-
68-
l_max_dLdt = Makie.lines!(ax1, ΔT_range, max_dLdt * 1e3, color = :thistle)
69-
l_max_dNdt = Makie.lines!(ax2, ΔT_range, max_dNdt * 1e-6, color = :thistle)
70-
71-
l_dLdt1 = Makie.lines!(ax1, ΔT_range, dLdt1 * 1e3, color = :skyblue)
72-
l_dNdt1 = Makie.lines!(ax2, ΔT_range, dNdt1 * 1e-6, color = :skyblue)
73-
74-
l_dLdt2 = Makie.lines!(ax1, ΔT_range, dLdt2 * 1e3, color = :blue3)
75-
l_dNdt2 = Makie.lines!(ax2, ΔT_range, dNdt2 * 1e-6, color = :blue3)
76-
77-
l_dLdt3 = Makie.lines!(ax1, ΔT_range, dLdt3 * 1e3, color = :orchid)
78-
l_dNdt3 = Makie.lines!(ax2, ΔT_range, dNdt3 * 1e-6, color = :orchid)
79-
80-
l_dLdt4 = Makie.lines!(ax1, ΔT_range, dLdt4 * 1e3, color = :purple)
81-
l_dNdt4 = Makie.lines!(ax2, ΔT_range, dNdt4 * 1e-6, color = :purple)
82-
83-
Makie.Legend(
84-
fig[1, 3],
85-
[l_max_dNdt, l_dNdt1, l_dNdt2, l_dNdt3, l_dNdt4],
86-
[
87-
"limit",
88-
"ρₐ=1.2 kg/m3, Fᵣ=0.8, ρᵣ=800kg/m3",
89-
"ρₐ=1.2 kg/m3, Fᵣ=0.2, ρᵣ=800kg/m3",
90-
"ρₐ=1.2 kg/m3, Fᵣ=0.2, ρᵣ=200kg/m3",
91-
"ρₐ=0.5 kg/m3, Fᵣ=0.2, ρᵣ=200kg/m3",
92-
],
93-
framevisible = false,
94-
)
95-
Makie.save("P3_ice_melt.svg", fig)
66+
Makie.with_theme(Makie.theme_minimal(), fontsize = 22, linewidth = 3) do
67+
fig = Makie.Figure(size = (800, 600))
68+
69+
ax1 = Makie.Axis(fig[1, 1];
70+
xlabel = "T [°C]", ylabel = "ice mass melting rate [g/m³/s]",
71+
title = "P3 ice melting",
72+
)
73+
ax2 = Makie.Axis(fig[1, 1];
74+
xlabel = "T [°C]", ylabel = "ice number melting rate [1/cm³/s]",
75+
yaxisposition = :right, rightspinevisible = true,
76+
)
77+
78+
l_max_dLdt = Makie.hlines!(ax1, max_dLdt * 1e3; color = :gray, linewidth = 1)
79+
l_max_dNdt = Makie.hlines!(ax2, max_dNdt * 1e-6; color = :gray, linewidth = 1, label = "limit")
80+
81+
l_dLdt1 = Makie.lines!(ax1, ΔT_range, dLdt1 * 1e3; color = :skyblue)
82+
l_dNdt1 = Makie.lines!(ax2, ΔT_range, dNdt1 * 1e-6; color = :skyblue, label = label1)
83+
84+
l_dLdt2 = Makie.lines!(ax1, ΔT_range, dLdt2 * 1e3; color = :blue3)
85+
l_dNdt2 = Makie.lines!(ax2, ΔT_range, dNdt2 * 1e-6; color = :blue3, label = label2)
86+
87+
l_dLdt3 = Makie.lines!(ax1, ΔT_range, dLdt3 * 1e3; color = :orchid)
88+
l_dNdt3 = Makie.lines!(ax2, ΔT_range, dNdt3 * 1e-6; color = :orchid, label = label3)
89+
90+
l_dLdt4 = Makie.lines!(ax1, ΔT_range, dLdt4 * 1e3; color = :purple)
91+
l_dNdt4 = Makie.lines!(ax2, ΔT_range, dNdt4 * 1e-6; color = :purple, label = label4)
92+
93+
Makie.axislegend(ax2; position = :rb, framevisible = false)
94+
Makie.save("P3_ice_melt.svg", fig)
95+
end

‎src/BulkMicrophysicsTendencies.jl‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -719,10 +719,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
719719
T_freeze = TDI.TD.Parameters.T_freeze(tps)
720720
if T > T_freeze && L_ice > zero(L_ice)
721721
state = CMP3.P3State(p3, L_ice, N_ice, F_rim, ρ_rim)
722-
# TODO: Using a function that takes dt as an argument is not compatible with the current API.
723-
# We should use a function that doesn't take dt as an argument.
724-
dt_dummy = 1000 * one(T) # P3 uses dt for limiting, we'll limit later
725-
melt = CMP3.ice_melt(vel, aps, tps, T, ρ, dt_dummy, state, logλ)
722+
melt = CMP3.ice_melt(vel, aps, tps, T, ρ, state, logλ)
726723

727724
# Melting converts ice to rain
728725
dq_ice_dt -= melt.dLdt / ρ

‎src/P3_processes.jl‎

Lines changed: 17 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,28 +1,25 @@
11
"""
22
het_ice_nucleation(pdf_c, p3, tps, q_lcl, N, T, ρₐ, p, aerosol)
33
4-
- aerosol - aerosol parameters (supported types: desert dust, illite, kaolinite)
5-
- tps - thermodynamics parameters
6-
- q_lcl - cloud liquid water specific content
7-
- N_lcl - cloud droplet number concentration
8-
- RH - relative humidity
9-
- T - temperature
10-
- ρₐ - air density
11-
- dt - model time step
4+
# Arguments
5+
- `aerosol`: aerosol parameters (supported types: desert dust, illite, kaolinite)
6+
- `tps`: thermodynamics parameters
7+
- `q_lcl`: cloud liquid water specific content
8+
- `N_lcl`: cloud droplet number concentration
9+
- `RH`: relative humidity
10+
- `T`: temperature
11+
- `ρₐ`: air density
12+
- `dt`: model time step
1213
1314
Returns a named tuple with ice number concentration and ice content
1415
hetergoeneous freezing rates from cloud droplets.
1516
"""
1617
function het_ice_nucleation(
1718
aerosol::Union{CMP.DesertDust, CMP.Illite, CMP.Kaolinite},
1819
tps::TDI.PS,
19-
q_lcl::FT,
20-
N_lcl::FT,
21-
RH::FT,
22-
T::FT,
23-
ρₐ::FT,
24-
dt::FT,
25-
) where {FT}
20+
q_lcl, N_lcl, RH, T, ρₐ, #dt,
21+
)
22+
FT = eltype(tps)
2623
#TODO - Also consider rain freezing
2724

2825
# Immersion freezing nucleation rate coefficient
@@ -40,22 +37,21 @@ function het_ice_nucleation(
4037
dNdt = max(0, dNdt)
4138
dLdt = max(0, dLdt)
4239
# ... and dont exceed the available number and mass of water droplets
43-
dNdt = min(dNdt, N_lcl / dt)
44-
dLdt = min(dLdt, q_lcl * ρₐ / dt)
40+
# dNdt = min(dNdt, N_lcl / dt)
41+
# dLdt = min(dLdt, q_lcl * ρₐ / dt)
4542

4643
return (; dNdt, dLdt)
4744
end
4845

4946
"""
50-
ice_melt(velocity_params::CMP.Chen2022VelType, aps, tps, Tₐ, ρₐ, dt, state, logλ; ∫kwargs...)
47+
ice_melt(velocity_params, aps, tps, Tₐ, ρₐ, state, logλ; ∫kwargs...)
5148
5249
# Arguments
5350
- `velocity_params`: [`CMP.Chen2022VelType`](@ref)
5451
- `aps`: [`CMP.AirProperties`](@ref)
5552
- `tps`: thermodynamics parameters
5653
- `Tₐ`: temperature (K)
5754
- `ρₐ`: air density
58-
- `dt`: model time step (for limiting the tendnecy)
5955
- `state`: a [`P3State`](@ref) object
6056
- `logλ`: the log of the slope parameter [log(1/m)]
6157
@@ -65,7 +61,8 @@ end
6561
Returns the melting rate of ice (QIMLT in Morrison and Mildbrandt (2015)).
6662
"""
6763
function ice_melt(
68-
velocity_params::CMP.Chen2022VelType, aps::CMP.AirProperties, tps::TDI.PS, Tₐ, ρₐ, dt, state::P3State, logλ;
64+
velocity_params, aps::CMP.AirProperties, tps::TDI.PS,
65+
Tₐ, ρₐ, state::P3State, logλ;
6966
∫kwargs...,
7067
)
7168
# Note: process not dependent on `F_liq`
@@ -93,9 +90,6 @@ function ice_melt(
9390
# compute change of N_ice proportional to change in L
9491
dNdt = N_ice / L_ice * dLdt
9592

96-
# ... and don't exceed the available number and mass of water droplets
97-
dNdt = min(dNdt, N_ice / dt) # TODO: Apply limiters in CA.jl
98-
dLdt = min(dLdt, L_ice / dt)
9993
return (; dNdt, dLdt)
10094
end
10195

0 commit comments

Comments
 (0)