diff --git a/box/Alpert_Knopf_2016_forward.jl b/box/Alpert_Knopf_2016_forward.jl
index f71bf849dc..aa407d2e66 100644
--- a/box/Alpert_Knopf_2016_forward.jl
+++ b/box/Alpert_Knopf_2016_forward.jl
@@ -46,7 +46,7 @@ Aj_sorted = sort(Aj_unsorted, rev = true)
 # initial condition for the problem
 IC = [T_initial, r_l, FT(N₀), FT(N_ice)]
 # additional model parameters
-p_def = (; tps, r_l, 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))
diff --git a/box/box.jl b/box/box.jl
index 37d59fc3c1..6e7f6035f1 100644
--- a/box/box.jl
+++ b/box/box.jl
@@ -11,12 +11,13 @@ ODE problem definitions
 function box_model(dY, Y, p, t)
 
     # Get simulation parameters
-    (; tps, r_l, 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
 
@@ -45,13 +46,13 @@ ODE problem definitions
 function box_model_with_probability(dY, Y, p, t)
 
     # Get simulation parameters
-    (; tps, const_dt, r_l, 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
-    r = Y[2]                      # total surface area converted to radius
+    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
-    drj_dt = FT(0)
+    dr_l_dt = FT(0)
     dN_ice_dt = FT(0)
     dN_liq_dt = FT(0)
     if N_liq > 0
@@ -79,18 +80,17 @@ function box_model_with_probability(dY, Y, p, t)
                 Aj_sorted[j] = FT(0)
             end
         end
-        A = N_liq * 4 * FT(π) * r^2
         Aj_sum = sum(Aj_sorted)
         rj_sum = sqrt(Aj_sum / 4 / FT(π))
-        drj_dt = (rj_sum - r * N_liq) / const_dt / N_liq
-        #drj_dt = - sqrt(abs(Aj_sum - A) / 4 / FT(π)) / const_dt / N_liq
+        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] = drj_dt         # radius
+    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