box tracks droplet radius instead of area

box/Alpert_Knopf_2016_forward.jl

@@ -21,6 +21,7 @@ include(joinpath(pkgdir(CM), "box", "box.jl"))
 
 # initial conditions following KA16 figure 4 (Cr1 case)
 A_droplet = FT(1e-5) * 1e-4 # INP surface area, m^2
+r_l = sqrt(A_droplet / 4 / FT(π))
 σg = 10
 N₀ = 1000
 N_ice = 0
@@ -42,10 +43,10 @@ end
 Aj_sorted = zeros(N₀)
 Aj_sorted = sort(Aj_unsorted, rev = true)
 
-# initial condition for the ODE problem
-IC = [T_initial, A_sum, FT(N₀), FT(N_ice)]
+# initial condition for the problem
+IC = [T_initial, r_l, FT(N₀), FT(N_ice)]
 # additional model parameters
-p_def = (; tps, A_droplet, aerosol, cooling_rate, N₀)
+p_def = (; tps, aerosol, cooling_rate, N₀)
 
 # run the box model without and with distribution sampling
 sol_cst = run_box(IC, t_0, t_end, (p_def..., const_dt = dt, flag = false))
@@ -82,7 +83,7 @@ MK.lines!(ax2, sol_var[1, :],          ff_var,               color = :green3, li
 idx = (sol_cst[3, :] .> 0)
 MK.lines!(ax3, sol_cst[1, idx], sol_cst[3, idx] .* A_droplet .* 1e4, color = :orange, linewidth = 3, label = "const A")
 idx = (sol_var[2, :] .> 0)
-MK.lines!(ax3, sol_var[1, idx], sol_var[2, idx] .* 1e4,              color = :green3, linewidth = 3, label = "variable A")
+MK.lines!(ax3, sol_var[1, idx], (4 .* FT(π) .* (sol_var[2, idx] .* sol_var[3, idx]).^2) .* 1e4, color = :green3, linewidth = 3, label = "variable A")
 MK.ylims!(ax3, 0.75e-6, 2e-1)
 
 MK.axislegend(ax1, framevisible = false, labelsize = 14, orientation = :horizontal, nbanks = 2, position = :rt)

box/box.jl

@@ -6,17 +6,18 @@ import CloudMicrophysics.Common as CMO
 import CloudMicrophysics.HetIceNucleation as CMI_het
 
 """
-    ODE problem definitions
+ODE problem definitions
 """
 function box_model(dY, Y, p, t)
 
     # Get simulation parameters
-    (; tps, A_droplet, aerosol, cooling_rate) = p
+    (; tps, aerosol, cooling_rate) = p
     # Numerical precision used in the simulation
     FT = eltype(Y)
 
     # Our state vector
     T = Y[1]          # temperature
+    r_l = Y[2]        # droplet radius
     N_liq = Y[3]      # number concentration of existing water droplets
     N_ice = Y[4]      # number concentration of activated ice crystals
 
@@ -28,7 +29,7 @@ function box_model(dY, Y, p, t)
     dN_ice_dt = FT(0)
     dN_liq_dt = FT(0)
     if N_liq > 0
-        dN_ice_dt = J_immer * N_liq * A_droplet
+        dN_ice_dt = J_immer * N_liq * 4 * FT(π) * r_l^2
         dN_liq_dt = -dN_ice_dt
     end
 
@@ -40,18 +41,18 @@ function box_model(dY, Y, p, t)
 end
 
 """
-    ODE problem definitions
+ODE problem definitions
 """
 function box_model_with_probability(dY, Y, p, t)
 
     # Get simulation parameters
-    (; tps, const_dt, A_droplet, aerosol, cooling_rate, Aj_sorted, N₀) = p
+    (; tps, const_dt, aerosol, cooling_rate, Aj_sorted, N₀) = p
     # Numerical precision used in the simulation
     FT = eltype(Y)
 
     # Our state vector
     T = Y[1]                      # temperature
-    A = Y[2]                      # available surface area for freezing
+    r_l = Y[2]                    # total surface area converted to radius
     N_liq = max(FT(0), Y[3])      # number concentration of existing water droplets
     N_ice = max(FT(0), Y[4])      # number concentration of activated ice crystals
 
@@ -61,7 +62,7 @@ function box_model_with_probability(dY, Y, p, t)
 
     # Update the tendecies
     dT_dt = -cooling_rate
-    dAj_dt = FT(0)
+    dr_l_dt = FT(0)
     dN_ice_dt = FT(0)
     dN_liq_dt = FT(0)
     if N_liq > 0
@@ -80,15 +81,17 @@ function box_model_with_probability(dY, Y, p, t)
             end
         end
         Aj_sum = sum(Aj_sorted)
-        dAj_dt = (Aj_sum - A) / const_dt
+        rj_sum = sqrt(Aj_sum / 4 / FT(π))
+        dr_l_dt = (rj_sum - r_l * N_liq) / const_dt / N_liq
+        #dr_l_dt = - sqrt(abs(Aj_sum - A) / 4 / FT(π)) / const_dt / N_liq
         dN_ice_dt = max(FT(0), n_frz / const_dt)
         dN_liq_dt = -dN_ice_dt
     end
 
     # Set tendencies
     dY[1] = dT_dt          # temperature
-    dY[2] = dAj_dt         # total Aj
-    dY[3] = dN_liq_dt      # mumber concentration of droplets
+    dY[2] = dr_l_dt        # droplet radius
+    dY[3] = dN_liq_dt      # number concentration of droplets
     dY[4] = dN_ice_dt      # number concentration of ice activated particles
 
 end