Merge pull request #282 from CliMA/al/melt_to_parc

lat heat of fus in parcel dT/dt and dSl/dt

docs/src/IceNucleationParcel0D.md

@@ -55,13 +55,15 @@ where:
 From the moist adiabatic assumption
 ```math
 \begin{equation}
-\frac{dT}{dt} = \frac{R_a T}{c_p p} \frac{dp}{dt} + \frac{L_v}{c_p} \frac{dq_l}{dt} + \frac{L_s}{c_p} \frac{d q_i}{dt}
+\frac{dT}{dt} = \frac{R_a T}{c_p p} \frac{dp}{dt} + \frac{L_v}{c_p} \frac{dq_{l,vap}}{dt} + \frac{L_s}{c_p} \frac{d q_{i,subl}}{dt} + \frac{L_f}{c_p} \frac{d q_{i,fus}}{dt}
 \end{equation}
 ```
 where:
-- ``q_l`` is the cloud liquid water specific humidity,
-- ``q_i`` is the cloud ice specific humidity,
-- ``L_s`` is the latent heat of sublimation.
+- ``q_{l,vap}`` is the cloud liquid water specific humidity from vaporization/condensation,
+- ``q_{i,subl}`` is the cloud ice specific humidity from sublimation/deposition,
+- ``q_{i,fus}`` is the cloud ice specific humidity from melting/freezing,
+- ``L_s`` is the latent heat of sublimation,
+- ``L_f`` is the latent heat of fusion.
 
 From hydrostatic balance and assuming constant vertical velocity:
 ```math
@@ -73,11 +75,11 @@ where:
 - ``g`` is the gravitational acceleration
 - ``w`` is the constant vertical velocity.
 
-Accounting for conservation of water, i.e. ``\frac{dq_v}{dt} + \frac{dq_l}{dt} + \frac{dq_i}{dt} = 0``,
+Accounting for conservation of water, i.e. ``\frac{dq_v}{dt} = - \frac{dq_{l,vap}}{dt} - \frac{dq_{i,subl}}{dt}``,
 and rearranging the terms
 ```math
 \begin{equation}
-\frac{dS_l}{dt} = a_1 w S_l - \left(a_2 + a_3 \right) S_l \frac{dq_l}{dt} - \left(a_2 + a_4\right) S_l \frac{dq_i}{dt}
+\frac{dS_l}{dt} = a_1 w S_l - \left(a_2 + a_3 \right)S_l  \frac{dq_{l,vap}}{dt} - \left(a_2 + a_4 \right) S_l \frac{dq_{i,subl}}{dt} - a_5 S_l \frac{dq_{i,fus}}{dt}
 \end{equation}
 ```
 where:
@@ -101,6 +103,11 @@ a_3 = \frac{L_v^2}{R_v T^2 c_p}
 a_4 = \frac{L_v L_s}{R_v T^2 c_p}
 \end{equation}
 ```
+```math
+\begin{equation}
+a_5 = \frac{L_v L_f}{R_v T^2 c_p}
+\end{equation}
+```
 
 Saturation ratio over ice can then be related to ``S_l`` by the relation
 ```math

parcel/parcel.jl

@@ -53,6 +53,7 @@ function parcel_model(dY, Y, p, t)
     R_a = TD.gas_constant_air(tps, q)
     cp_a = TD.cp_m(tps, q)
     L_subl = TD.latent_heat_sublim(tps, T)
+    L_fus = TD.latent_heat_fusion(tps, T)
     L_vap = TD.latent_heat_vapor(tps, T)
     ρ_air = TD.air_density(tps, ts)
     e = q_vap * p_a * R_v / R_a
@@ -62,6 +63,7 @@ function parcel_model(dY, Y, p, t)
     a2 = 1 / q_vap
     a3 = L_vap^2 / R_v / T^2 / cp_a
     a4 = L_vap * L_subl / R_v / T^2 / cp_a
+    a5 = L_vap * L_fus / R_v / cp_a / (T^2)
 
     # TODO - we should zero out all tendencies and augemnt them
     # TODO - add immersion, homogeneous, ...
@@ -135,12 +137,16 @@ function parcel_model(dY, Y, p, t)
         end
     end
 
+    dq_liq_dt_vap_to_liq = dql_dt_cond        # from liq-vap transitions
+    dq_ice_dt_vap_to_ice = dqi_dt_new_depo + dqi_dt_depo        # from ice-vap transitions
+    dq_ice_dt_liq_to_ice = dqi_dt_new_immers   # from ice-liq transitions
+
     # Update the tendecies
-    dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_immers
-    dq_liq_dt = dql_dt_cond - dqi_dt_new_immers
-    dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt - (a2 + a4) * S_liq * dq_ice_dt
+    dq_ice_dt = dq_ice_dt_vap_to_ice + dq_ice_dt_liq_to_ice
+    dq_liq_dt = dq_liq_dt_vap_to_liq - dq_ice_dt_liq_to_ice
+    dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt_vap_to_liq - (a2 + a4) * S_liq * dq_ice_dt_vap_to_ice - a5 * S_liq * dq_ice_dt_liq_to_ice
     dp_a_dt = -p_a * grav / R_a / T * w
-    dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt + L_subl / cp_a * dq_ice_dt
+    dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_vap_to_liq + L_fus / cp_a * dq_ice_dt_liq_to_ice + L_subl / cp_a * dq_ice_dt_vap_to_ice
     dq_vap_dt = -dq_ice_dt - dq_liq_dt
 
     # Set tendencies