Skip to content

Commit 1f5f9a8

Browse files
haakon-ehaakon-e
committed
todo p3
1 parent 8ff445c commit 1f5f9a8

7 files changed

Lines changed: 281 additions & 201 deletions

File tree

‎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.png", fig)
95+
end

‎src/BulkMicrophysicsTendencies.jl‎

Lines changed: 66 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ import ..Microphysics1M as CM1
3232
import ..Microphysics2M as CM2
3333
import ..MicrophysicsNonEq as CMNonEq
3434
import ..P3Scheme as CMP3
35+
import ..HetIceNucleation as CM_HetIce
3536
import ...ThermodynamicsInterface as TDI
3637

3738
export MicrophysicsScheme,
@@ -665,6 +666,9 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
665666
vel = mp.ice.terminal_velocity
666667
pdf_c = mp.ice.cloud_pdf
667668
pdf_r = mp.ice.rain_pdf
669+
heterogeneous = mp.ice.heterogeneous
670+
deposition_condfreeze = mp.ice.deposition_condfreeze
671+
668672

669673
# Only compute ice processes if there is ice mass/number present
670674
if (q_ice > zero(q_ice) || n_ice > zero(n_ice))
@@ -684,23 +688,8 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
684688
# Only compute if there is ice present
685689
if L_ice > zero(L_ice) && N_ice > zero(N_ice)
686690
coll = CMP3.bulk_liquid_ice_collision_sources(
687-
p3,
688-
logλ,
689-
L_ice,
690-
N_ice,
691-
F_rim,
692-
ρ_rim,
693-
pdf_c,
694-
pdf_r,
695-
L_lcl,
696-
N_lcl,
697-
L_rai,
698-
N_rai,
699-
aps,
700-
tps,
701-
vel,
702-
ρ,
703-
T,
691+
p3, logλ, L_ice, N_ice, F_rim, ρ_rim, pdf_c, pdf_r,
692+
L_lcl, N_lcl, L_rai, N_rai, aps, tps, vel, ρ, T,
704693
)
705694

706695
# Add collision tendencies
@@ -719,10 +708,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
719708
T_freeze = TDI.TD.Parameters.T_freeze(tps)
720709
if T > T_freeze && L_ice > zero(L_ice)
721710
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λ)
711+
melt = CMP3.ice_melt(vel, aps, tps, T, ρ, state, logλ)
726712

727713
# Melting converts ice to rain
728714
dq_ice_dt -= melt.dLdt / ρ
@@ -733,9 +719,65 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
733719
end
734720
end
735721

736-
# TODO: When ice number concentration is tracked, add dn_ice_dt to return tuple:
737-
# return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt, dq_ice_dt, dn_ice_dt, dq_rim_dt, db_rim_dt)
738-
return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt, dq_ice_dt, dq_rim_dt, db_rim_dt)
722+
# --- -------------------- ---
723+
# --- Ice Nucleation Modes ---
724+
# --- -------------------- ---
725+
# TODO: These need to be rates, ∂N/∂t, not number changes ΔN
726+
727+
# --- Ice Nucleation (Deposition) ---
728+
# Assume nucleated particles have diameter 1μm --> nucleated mass (per particle)
729+
D_nuc_ice = 1e-6 # 1μm
730+
m_nuc = p3.ρ_i * CO.volume_sphere_D(D_nuc_ice) # [kg]
731+
732+
n_dep = CM_HetIce.P3_deposition_N_i(deposition_condfreeze, T) / ρ # [particles / kg air]
733+
m_dep = n_dep * m_nuc # [kg ice / kg air]
734+
dn_ice_dt += n_dep
735+
dq_ice_dt += m_dep
736+
dq_rim_dt += m_dep
737+
db_rim_dt += m_dep / 900
738+
739+
# --- Heterogeneous Ice Nucleation (Immersion freezing) ---
740+
# Assume mass loss is mean condensate mass
741+
m_lcl = q_lcl / n_lcl # mean liquid mass
742+
743+
inpc = CM_HetIce.INP_concentration_mean(heterogeneous, T) / ρ # [particles / kg air]
744+
n_het = max(FT(0), inpc - n_ice) # [particles / kg air]
745+
q_het = n_het * m_lcl # [kg ice / kg air]
746+
747+
dn_ice_dt += n_het
748+
dq_ice_dt += q_het
749+
dq_rim_dt += q_het
750+
db_rim_dt += q_het / 900 # TODO: make this a parameter?
751+
dn_lcl_dt -= n_het
752+
dq_lcl_dt -= q_het
753+
754+
# --- Cloud Droplet Condensation Freezing ---
755+
# ref: `homogeneous_freezing` in `parcel/ParcelTendencies.jl`
756+
# get mean diameter of cloud droplets, then convert to volume
757+
758+
# --- Ice Sublimation / Deposition ---
759+
# Deposition/sublimation of cloud ice
760+
∂ₜq_ice_dep = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(icl, tps, q_tot, q_lcl, q_ice, q_rai, zero(q_ice), ρ, T)
761+
# No ice deposition above freezing (lack of INPs)
762+
∂ₜq_ice_dep = ifelse(T > tps.T_freeze, min(∂ₜq_ice_dep, zero(T)), ∂ₜq_ice_dep)
763+
∂ₜn_ice_dep = ifelse(∂ₜq_ice_dep < 0, (N_ice / q_ice) * ∂ₜq_ice_dep, zero(∂ₜq_ice_dep))
764+
dq_ice_dt += ∂ₜq_ice_dep
765+
dn_ice_dt += ∂ₜn_ice_dep
766+
767+
# --- Ice Self-collection (Aggregation) ---
768+
# TODO: Implement P3 ice self-collection (aggregation)
769+
# This process is currently missing in P3_processes.jl
770+
# S_ice_agg = CMP3.ice_self_collection(p3, q_ice, n_ice, ...)
771+
# dn_ice_dt -= S_ice_agg
772+
773+
# --- Rain Heterogeneous Freezing ---
774+
# TODO: Implement heterogeneous freezing of rain
775+
# This process is currently missing in P3_processes.jl
776+
# S_rai_frz = ...
777+
# dq_rai_dt -= S_rai_frz
778+
# dq_ice_dt += S_rai_frz
779+
780+
return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt, dq_ice_dt, dn_ice_dt, dq_rim_dt, db_rim_dt)
739781
end
740782

741783
end # module BulkMicrophysicsTendencies

0 commit comments

Comments
 (0)