Fix bug in saturation adjustment solver for MoistAirBuoyancy (#120)

This PR fixes a bug in the saturation adjustment solver that only gets hit under certain conditions: when the secant iteration accesses unsaturated states during the course of finding a saturated state. Before this PR, we assume that during adjustment the state was saturated. However, this assumption does not hold if "intermediate" (eg test) states that are accessed during the iteration are not saturated. The result of the bug is a spurious change in the total moisture fraction. It is fixed in this PR.

Author: glwagner

src/MoistAirBuoyancies.jl

@@ -203,22 +203,15 @@ Solution of ``r(T) = 0`` is found via the [secant method](https://en.wikipedia.o
     q₁ = MoistureMassFractions(qᵛ⁺₁, qˡ₁, zero(qˡ₁))
     𝒰₁ = with_moisture(𝒰₀, q₁)
 
-    # We generate a second guess simply by adding 1 K to T₁...
-
-    # NOTE: We could also generate a second guess to start a secant iteration
+    # We generate a second guess to start a secant iteration
     # by applying the potential temperature assuming a liquid fraction
     # associated with T₁. This should represent an _overestimate_,
     # since ``qᵛ⁺₁(T₁)`` underestimates the saturation specific humidity,
     # and therefore qˡ₁ is overestimated. This is similar to an approach
     # used in Pressel et al 2015. However, it doesn't work for large liquid fractions.
-    T₂ = T₁ + 1
-
-    #=
     ℒˡᵣ = thermo.liquid.reference_latent_heat
     cᵖᵐ = mixture_heat_capacity(q₁, thermo)
     T₂ = T₁ + ℒˡᵣ * qˡ₁ / cᵖᵐ
-    =#
-
     𝒰₂ = adjust_state(𝒰₁, T₂, thermo)
 
     # Initialize saturation adjustment
@@ -267,7 +260,8 @@ end
     qᵛ⁺ = adjustment_saturation_specific_humidity(T, 𝒰₀, thermo)
     qᵗ = total_moisture_mass_fraction(𝒰₀)
     qˡ = max(0, qᵗ - qᵛ⁺)
-    q₁ = MoistureMassFractions(qᵛ⁺, qˡ, zero(qˡ))
+    qᵛ = qᵗ - qˡ
+    q₁ = MoistureMassFractions(qᵛ, qˡ, zero(qˡ))
     return with_moisture(𝒰₀, q₁)
 end