Relative Humidity -> Relative Super Saturation

CliMA · Dec 19, 2019 · c87d8f9 · c87d8f9
1 parent 5d962ae
commit c87d8f9
Show file tree

Hide file tree

Showing 4 changed files with 169 additions and 115 deletions.
diff --git a/src/Common/MoistThermodynamics/MoistThermodynamics.jl b/src/Common/MoistThermodynamics/MoistThermodynamics.jl
@@ -826,7 +826,6 @@ function saturation_adjustment_q_tot_θ_liq_ice_given_pressure(θ_liq_ice::FT, p
     T_2 = bound_upper_temperature(T_1, T_2)
     sol = find_zero(
       T -> liquid_ice_pottemp_sat(T, air_density(T, p, PhasePartition(q_tot)), q_tot) - θ_liq_ice,
-      # T -> liquid_ice_pottemp_sat_given_pressure(T, p, q_tot, tol, maxiter) - θ_liq_ice,
       T_1, T_2, SecantMethod(), CompactSolution(), tol, maxiter)
     if !sol.converged
       error("saturation_adjustment_q_tot_θ_liq_ice_given_pressure did not converge")
@@ -1027,27 +1026,6 @@ The saturated liquid ice potential temperature where
 liquid_ice_pottemp_sat(T::FT, ρ::FT, q_tot::FT) where {FT<:Real} =
     liquid_ice_pottemp(T, ρ, PhasePartition_equil(T, ρ, q_tot))
 
-"""
-    liquid_ice_pottemp_sat_given_pressure(T, p, q_tot)
-
-The saturated liquid ice potential temperature where
-
- - `T` temperature
- - `p` pressure
- - `q_tot` total specific humidity
- - `tol` tolerance for non-linear equation solve
- - `maxiter` maximum iterations for non-linear equation solve
-"""
-function liquid_ice_pottemp_sat_given_pressure(T::FT, p::FT, q_tot::FT, tol::FT, maxiter::Int) where {FT<:Real}
-  sol = find_zero(ρ_ -> ρ_ - air_density(T, p, PhasePartition_equil(T, ρ_, q_tot)),
-    FT(0.1), FT(1.2), SecantMethod(), CompactSolution(), tol, maxiter)
-  ρ = sol.root
-  if !sol.converged
-    error("liquid_ice_pottemp_sat_given_pressure did not converge")
-  end
-  return liquid_ice_pottemp(T, ρ, PhasePartition_equil(T, ρ, q_tot))
-end
-
 """
     liquid_ice_pottemp_sat(ts::ThermodynamicState)
 

diff --git a/src/Common/MoistThermodynamics/states.jl b/src/Common/MoistThermodynamics/states.jl
@@ -78,7 +78,7 @@ end
 function PhaseEquil(e_int::FT,
                     ρ::FT,
                     q_tot::FT,
-                    tol::FT=FT(1e-2),
+                    tol::FT=FT(1e-1),
                     maxiter::Int=3,
                     sat_adjust::F=saturation_adjustment) where {FT<:Real,F}
     return PhaseEquil{FT}(e_int, ρ, q_tot, sat_adjust(e_int, ρ, q_tot, tol, maxiter))
@@ -118,8 +118,8 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
 function LiquidIcePotTempSHumEquil(θ_liq_ice::FT,
                                    ρ::FT,
                                    q_tot::FT,
-                                   tol::FT=FT(1e-2),
-                                   maxiter::Int=5
+                                   tol::FT=FT(1e-1),
+                                   maxiter::Int=30
                                    ) where {FT<:Real}
     T = saturation_adjustment_q_tot_θ_liq_ice(θ_liq_ice, ρ, q_tot, tol, maxiter)
     q_pt = PhasePartition_equil(T, ρ, q_tot)
@@ -141,8 +141,8 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
 function LiquidIcePotTempSHumEquil_given_pressure(θ_liq_ice::FT,
                                                   p::FT,
                                                   q_tot::FT,
-                                                  tol::FT=FT(1e-2),
-                                                  maxiter::Int=10) where {FT<:Real}
+                                                  tol::FT=FT(1e-1),
+                                                  maxiter::Int=30) where {FT<:Real}
     T = saturation_adjustment_q_tot_θ_liq_ice_given_pressure(θ_liq_ice, p, q_tot, tol, maxiter)
     ρ = air_density(T, p, PhasePartition(q_tot))
     q = PhasePartition_equil(T, ρ, q_tot)
@@ -205,8 +205,8 @@ and, optionally
 function LiquidIcePotTempSHumNonEquil(θ_liq_ice::FT,
                                       ρ::FT,
                                       q_pt::PhasePartition{FT},
-                                      tol::FT=FT(1e-2),
-                                      maxiter::Int=10
+                                      tol::FT=FT(1e-1),
+                                      maxiter::Int=5
                                       ) where {FT<:Real}
     T = air_temperature_from_liquid_ice_pottemp_non_linear(θ_liq_ice, ρ, tol, maxiter, q_pt)
     e_int = internal_energy(T, q_pt)
@@ -265,9 +265,13 @@ Note that the output vectors are of size ``n*n_RH``, and they
 should span the input arguments to all of the constructors.
 """
 function tested_convergence_range(FT, n)
-  n_RH = 30
-  z_range  = range(FT(0), stop=FT(25*10^3), length=n)
-  RH       = range(FT(0), stop=FT(1.2), length=n_RH)
+  n_RS1 = 10
+  n_RS2 = 20
+  n_RS = n_RS1+n_RS2
+  z_range  = range(FT(0), stop=FT(2.5e4), length=n)
+  relative_sat1 = range(FT(0), stop=FT(1), length=n_RS1)
+  relative_sat2 = range(FT(1), stop=FT(1.02), length=n_RS2)
+  relative_sat = [relative_sat1...,relative_sat2...]
   T_min = FT(150)
   T_surface = FT(350)
 
@@ -276,14 +280,14 @@ function tested_convergence_range(FT, n)
           getindex.(args, 2),
           getindex.(args, 3)
 
-  p  = collect(Iterators.flatten([p  for RH_ in 1:n_RH]))
-  ρ  = collect(Iterators.flatten([ρ  for RH_ in 1:n_RH]))
-  T  = collect(Iterators.flatten([T  for RH_ in 1:n_RH]))
-  RH = collect(Iterators.flatten([RH for RH_ in 1:n]))
+  p  = collect(Iterators.flatten([p  for RS in 1:n_RS]))
+  ρ  = collect(Iterators.flatten([ρ  for RS in 1:n_RS]))
+  T  = collect(Iterators.flatten([T  for RS in 1:n_RS]))
+  relative_sat = collect(Iterators.flatten([relative_sat for RS in 1:n]))
 
   # Additional variables
   q_sat = q_vap_saturation.(T, ρ)
-  q_tot = min.(RH .*q_sat, FT(1))
+  q_tot = min.(relative_sat .*q_sat, FT(1))
   q_pt = PhasePartition_equil.(T, ρ, q_tot)
   e_int = internal_energy.(T, q_pt)
   θ_liq_ice = liquid_ice_pottemp.(T, ρ, q_pt)

diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl
@@ -90,7 +90,7 @@ function compute_thermo!(FT, state, k, ijk, ev, e, z, zvals, thermoQ)
 
     e_int = e_tot - 1//2 * (u^2 + v^2 + w^2) - grav * z
 
-    ts = PhaseEquil(convert(FT, e_int), state.ρ, q_tot)
+    ts = PhaseEquil(convert(FT, e_int), state.ρ, q_tot, FT(1e-2), 3)
     Phpart = PhasePartition(ts)
 
     th = thermo_vars(thermoQ[ijk,e])

diff --git a/test/Common/MoistThermodynamics/runtests.jl b/test/Common/MoistThermodynamics/runtests.jl
@@ -5,6 +5,8 @@ MT = MoistThermodynamics
 using CLIMA.PlanetParameters
 using LinearAlgebra
 
+float_types = [Float32, Float64]
+
 @testset "moist thermodynamics - correctness" begin
   FT = Float64
   # ideal gas law
@@ -124,85 +126,155 @@ using LinearAlgebra
 end
 
 
+@testset "moist thermodynamics - default behavior accuracy" begin
+  # Input arguments should be accurate within machine precision
+  # Temperature is approximated via saturation adjustment, and should be within a physical tolerance
+
+  for FT in float_types
+    rtol = FT(1e-2)
+    e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_convergence_range(FT, 50)
+
+    # PhaseEquil
+    ts_exact = PhaseEquil.(e_int, ρ, q_tot, FT(1e-4), 100)
+    ts = PhaseEquil.(e_int, ρ, q_tot)
+    # Should be machine accurate (because ts contains `e_int`,`ρ`,`q_tot`):
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(PhasePartition.(ts_exact),:tot))
+    @test all(internal_energy.(ts) .≈ internal_energy.(ts_exact))
+    @test all(air_density.(ts) .≈ air_density.(ts_exact))
+    # Approximate (temperature must be computed via saturation adjustment):
+    @test all(isapprox.(air_pressure.(ts), air_pressure.(ts_exact), rtol=rtol))
+    @test all(isapprox.(air_temperature.(ts), air_temperature.(ts_exact), rtol=rtol))
+
+    # PhaseEquil
+    ts_exact = PhaseEquil.(e_int, ρ, q_tot, FT(1e-4), 100, MT.saturation_adjustment_SecantMethod)
+    ts = PhaseEquil.(e_int, ρ, q_tot, FT(1e-1), 30, MT.saturation_adjustment_SecantMethod) # Needs to be in sync with default
+    # Should be machine accurate (because ts contains `e_int`,`ρ`,`q_tot`):
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(PhasePartition.(ts_exact),:tot))
+    @test all(internal_energy.(ts) .≈ internal_energy.(ts_exact))
+    @test all(air_density.(ts) .≈ air_density.(ts_exact))
+    # Approximate (temperature must be computed via saturation adjustment):
+    @test all(isapprox.(air_pressure.(ts), air_pressure.(ts_exact), rtol=rtol))
+    @test all(isapprox.(air_temperature.(ts), air_temperature.(ts_exact), rtol=rtol))
+
+    # LiquidIcePotTempSHumEquil
+    ts_exact = LiquidIcePotTempSHumEquil.(θ_liq_ice, ρ, q_tot, FT(1e-3), 40)
+    ts = LiquidIcePotTempSHumEquil.(θ_liq_ice, ρ, q_tot)
+    # Should be machine accurate:
+    @test all(air_density.(ts) .≈ air_density.(ts_exact))
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(PhasePartition.(ts_exact),:tot))
+    # Approximate (temperature must be computed via saturation adjustment):
+    @test all(isapprox.(internal_energy.(ts), internal_energy.(ts_exact), rtol=rtol))
+    @test all(isapprox.(liquid_ice_pottemp.(ts), liquid_ice_pottemp.(ts_exact), rtol=rtol))
+    @test all(isapprox.(air_temperature.(ts), air_temperature.(ts_exact), rtol=rtol))
+
+    # LiquidIcePotTempSHumEquil_given_pressure
+    ts_exact = LiquidIcePotTempSHumEquil_given_pressure.(θ_liq_ice, p, q_tot, FT(1e-3), 40)
+    ts = LiquidIcePotTempSHumEquil_given_pressure.(θ_liq_ice, p, q_tot)
+    # Should be machine accurate:
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(PhasePartition.(ts_exact),:tot))
+    # Approximate (temperature must be computed via saturation adjustment):
+    @test all(isapprox.(air_density.(ts), air_density.(ts_exact), rtol=rtol))
+    @test all(isapprox.(internal_energy.(ts), internal_energy.(ts_exact), rtol=rtol))
+    @test all(isapprox.(liquid_ice_pottemp.(ts), liquid_ice_pottemp.(ts_exact), rtol=rtol))
+    @test all(isapprox.(air_temperature.(ts), air_temperature.(ts_exact), rtol=rtol))
+
+    # LiquidIcePotTempSHumNonEquil
+    ts_exact = LiquidIcePotTempSHumNonEquil.(θ_liq_ice, ρ, q_pt, FT(1e-3), 40)
+    ts = LiquidIcePotTempSHumNonEquil.(θ_liq_ice, ρ, q_pt)
+    # Should be machine accurate:
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(PhasePartition.(ts_exact),:tot))
+    @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(PhasePartition.(ts_exact),:liq))
+    @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(PhasePartition.(ts_exact),:ice))
+    @test all(air_density.(ts) .≈ air_density.(ts_exact))
+    # Approximate (temperature must be computed via non-linear solve):
+    @test all(isapprox.(internal_energy.(ts), internal_energy.(ts_exact), rtol=rtol))
+    @test all(isapprox.(liquid_ice_pottemp.(ts), liquid_ice_pottemp.(ts_exact), rtol=rtol))
+    @test all(isapprox.(air_temperature.(ts), air_temperature.(ts_exact), rtol=rtol))
+
+  end
+
+end
+
 @testset "moist thermodynamics - constructor consistency" begin
 
   # Make sure `ThermodynamicState` arguments are returned unchanged
 
-  FT = Float64
-  e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_convergence_range(FT, 50)
-
-  # PhaseDry
-  ts = PhaseDry.(e_int, ρ)
-  @test all(internal_energy.(ts) .≈ e_int)
-  @test all(air_density.(ts) .≈ ρ)
-
-  # PhaseEquil
-  ts = PhaseEquil.(e_int, ρ, q_tot, FT(1e-2), 15, Ref(MT.saturation_adjustment_SecantMethod))
-  @test all(internal_energy.(ts) .≈ e_int)
-  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot)
-  @test all(air_density.(ts) .≈ ρ)
-
-  ts = PhaseEquil.(e_int, ρ, q_tot, FT(1e-2), 15)
-  @test all(internal_energy.(ts) .≈ e_int)
-  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot)
-  @test all(air_density.(ts) .≈ ρ)
-
-  # PhaseNonEquil
-  ts = PhaseNonEquil.(e_int, ρ, q_pt)
-  @test all(internal_energy.(ts) .≈ e_int)
-  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt, :tot))
-  @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt, :liq))
-  @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt, :ice))
-  @test all(air_density.(ts) .≈ ρ)
-
-  # air_temperature_from_liquid_ice_pottemp_given_pressure-liquid_ice_pottemp inverse
-  θ_liq_ice_ = liquid_ice_pottemp_given_pressure.(T, p, q_pt)
-  @test all(air_temperature_from_liquid_ice_pottemp_given_pressure.(θ_liq_ice_, p, q_pt) .≈ T)
-
-  # liquid_ice_pottemp-air_temperature_from_liquid_ice_pottemp_given_pressure inverse
-  T = air_temperature_from_liquid_ice_pottemp_given_pressure.(θ_liq_ice, p, q_pt)
-  @test all(liquid_ice_pottemp_given_pressure.(T, p, q_pt) .≈ θ_liq_ice)
-
-  # Accurate but expensive `LiquidIcePotTempSHumNonEquil` constructor (Non-linear temperature from θ_liq_ice)
-  T_non_linear = air_temperature_from_liquid_ice_pottemp_non_linear.(θ_liq_ice, ρ, 1e-3, 10, q_pt)
-  e_int_ = internal_energy.(T_non_linear, q_pt)
-  ts = PhaseNonEquil.(e_int_, ρ, q_pt)
-  @test all(T_non_linear .≈ air_temperature.(ts))
-  @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts))
-
-  # LiquidIcePotTempSHumEquil
-  ts = LiquidIcePotTempSHumEquil.(θ_liq_ice, ρ, q_tot, 1e-4, 10)
-  @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice)
-  @test all(air_density.(ts) .≈ ρ)
-  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot)
-
-  # The LiquidIcePotTempSHumEquil_given_pressure constructor
-  # passes the consistency test within sufficient physical precision,
-  # however, it fails to satisfy the consistency test within machine
-  # precision for the input pressure.
-
-  # LiquidIcePotTempSHumEquil_given_pressure
-  ts = LiquidIcePotTempSHumEquil_given_pressure.(θ_liq_ice, p, q_tot, FT(1e-4), 20)
-  @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice)
-  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt, :tot))
-  # @test all(air_pressure.(ts) .≈ p) # TODO: would be nice to get working
-  @test max(abs.(air_pressure.(ts) .- p)...) < 1100 # Best we can do right now
-
-  # LiquidIcePotTempSHumNonEquil_given_pressure
-  ts = LiquidIcePotTempSHumNonEquil_given_pressure.(θ_liq_ice, p, q_pt)
-  @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice)
-  @test all(air_pressure.(ts) .≈ p)
-  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt,:tot))
-  @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt,:liq))
-  @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt,:ice))
-
-  # LiquidIcePotTempSHumNonEquil
-  ts = LiquidIcePotTempSHumNonEquil.(θ_liq_ice, ρ, q_pt, FT(1e-3))
-  @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts))
-  @test all(air_density.(ts) .≈ ρ)
-  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt,:tot))
-  @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt,:liq))
-  @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt,:ice))
+  for FT in float_types
+    rtol = FT(1e-2)
+    e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_convergence_range(FT, 50)
+
+    # PhaseDry
+    ts = PhaseDry.(e_int, ρ)
+    @test all(internal_energy.(ts) .≈ e_int)
+    @test all(air_density.(ts) .≈ ρ)
+
+    # PhaseEquil
+    ts = PhaseEquil.(e_int, ρ, q_tot, FT(1e-1), 30, Ref(MT.saturation_adjustment_SecantMethod))
+    @test all(internal_energy.(ts) .≈ e_int)
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot)
+    @test all(air_density.(ts) .≈ ρ)
+
+    ts = PhaseEquil.(e_int, ρ, q_tot)
+    @test all(internal_energy.(ts) .≈ e_int)
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot)
+    @test all(air_density.(ts) .≈ ρ)
+
+    # PhaseNonEquil
+    ts = PhaseNonEquil.(e_int, ρ, q_pt)
+    @test all(internal_energy.(ts) .≈ e_int)
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt, :tot))
+    @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt, :liq))
+    @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt, :ice))
+    @test all(air_density.(ts) .≈ ρ)
+
+    # air_temperature_from_liquid_ice_pottemp_given_pressure-liquid_ice_pottemp inverse
+    θ_liq_ice_ = liquid_ice_pottemp_given_pressure.(T, p, q_pt)
+    @test all(air_temperature_from_liquid_ice_pottemp_given_pressure.(θ_liq_ice_, p, q_pt) .≈ T)
+
+    # liquid_ice_pottemp-air_temperature_from_liquid_ice_pottemp_given_pressure inverse
+    T = air_temperature_from_liquid_ice_pottemp_given_pressure.(θ_liq_ice, p, q_pt)
+    @test all(liquid_ice_pottemp_given_pressure.(T, p, q_pt) .≈ θ_liq_ice)
+
+    # Accurate but expensive `LiquidIcePotTempSHumNonEquil` constructor (Non-linear temperature from θ_liq_ice)
+    T_non_linear = air_temperature_from_liquid_ice_pottemp_non_linear.(θ_liq_ice, ρ, FT(1e-3), 10, q_pt)
+    e_int_ = internal_energy.(T_non_linear, q_pt)
+    ts = PhaseNonEquil.(e_int_, ρ, q_pt)
+    @test all(T_non_linear .≈ air_temperature.(ts))
+    @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts))
+
+    # LiquidIcePotTempSHumEquil
+    ts = LiquidIcePotTempSHumEquil.(θ_liq_ice, ρ, q_tot, FT(1e-3), 40)
+    @test all(isapprox.(liquid_ice_pottemp.(ts), θ_liq_ice, atol=1e-1))
+    @test all(isapprox.(air_density.(ts), ρ, rtol=rtol))
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot)
+
+    # The LiquidIcePotTempSHumEquil_given_pressure constructor
+    # passes the consistency test within sufficient physical precision,
+    # however, it fails to satisfy the consistency test within machine
+    # precision for the input pressure.
+
+    # LiquidIcePotTempSHumEquil_given_pressure
+    ts = LiquidIcePotTempSHumEquil_given_pressure.(θ_liq_ice, p, q_tot, FT(1e-3), 40)
+    @test all(isapprox.(liquid_ice_pottemp.(ts), θ_liq_ice, atol=1e-1))
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt, :tot))
+    @test all(isapprox.(air_pressure.(ts), p, atol=FT(MSLP)*2e-2))
+
+    # LiquidIcePotTempSHumNonEquil_given_pressure
+    ts = LiquidIcePotTempSHumNonEquil_given_pressure.(θ_liq_ice, p, q_pt)
+    @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice)
+    @test all(air_pressure.(ts) .≈ p)
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt,:tot))
+    @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt,:liq))
+    @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt,:ice))
+
+    # LiquidIcePotTempSHumNonEquil
+    ts = LiquidIcePotTempSHumNonEquil.(θ_liq_ice, ρ, q_pt, FT(1e-3))
+    @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts))
+    @test all(air_density.(ts) .≈ ρ)
+    @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt,:tot))
+    @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt,:liq))
+    @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt,:ice))
+  end
 
 end
 
@@ -222,8 +294,8 @@ end
   ts_eq              = PhaseEquil.(e_int, ρ, q_tot, FT(1e-1), 15)
   ts_T               = TemperatureSHumEquil.(air_temperature.(ts_dry), air_pressure.(ts_dry), q_tot)
   ts_neq             = PhaseNonEquil.(e_int, ρ, q_pt)
-  ts_θ_liq_ice_eq    = LiquidIcePotTempSHumEquil.(θ_liq_ice, ρ, q_tot, FT(1e-2), 15)
-  ts_θ_liq_ice_eq_p  = LiquidIcePotTempSHumEquil_given_pressure.(θ_liq_ice, p, q_tot, FT(1e-2), 20)
+  ts_θ_liq_ice_eq    = LiquidIcePotTempSHumEquil.(θ_liq_ice, ρ, q_tot, FT(1e-3), 40)
+  ts_θ_liq_ice_eq_p  = LiquidIcePotTempSHumEquil_given_pressure.(θ_liq_ice, p, q_tot, FT(1e-3), 40)
   ts_θ_liq_ice_neq   = LiquidIcePotTempSHumNonEquil.(θ_liq_ice, ρ, q_pt)
   ts_θ_liq_ice_neq_p = LiquidIcePotTempSHumNonEquil_given_pressure.(θ_liq_ice, p, q_pt)