Enhancements on Clausius-Clapeyron section in Docs (#54)

* format docstring

* enhance docstring + code to avoid confusion of ℒᵛ(T=0) and ℒ₀

* cleaner

* cleaner

* homogenize docs, docstring and code

* clarify bc

* avoid confusing subscript 0 with reference temperature

* add comma

* use pᵛ⁺ notation

* plot saturation vapor pressure

* Update CI workflow to include new Go versions

Author: navidcy

.github/workflows/CI.yml

@@ -40,8 +40,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - "min"
-          - "1"
+          - "1.10"
+          - "1.11"
+          - "1.12"
         os:
           - ubuntu-latest
           - macOS-latest

docs/src/thermodynamics.md

@@ -269,14 +269,14 @@ heats, the latent heat of a phase transition is linear in temperature.
 For example, for phase change from vapor to liquid,
 
 ```math
-ℒˡ(T) = ℒˡ₀ + \big ( \underbrace{cᵖᵛ - cᵖˡ}_{≡Δcˡ} \big ) T ,
+ℒˡ(T) = ℒˡ_{0K} + \big ( \underbrace{cᵖᵛ - cᵖˡ}_{≡Δcˡ} \big ) T ,
 ```
 
-where ``ℒˡ₀`` is the latent heat at ``T = 0``, with ``T`` in Kelvin.
-Integrate that to get
+where ``ℒˡ_{0K}`` is the latent heat at absolute zero, ``T = 0 \; \mathrm{K}``.
+By integrating from the triple-point temperature ``Tᵗʳ`` for which ``p(Tᵗʳ) = pᵗʳ``, we get
 
 ```math
-pᵛ⁺(T) = pᵗʳ \left ( \frac{T}{Tᵗʳ} \right )^{Δcˡ / Rᵛ} \exp \left [ \frac{ℒˡ₀}{Rᵛ} \left (\frac{1}{Tᵗʳ} - \frac{1}{T} \right ) \right ] .
+pᵛ⁺(T) = pᵗʳ \left ( \frac{T}{Tᵗʳ} \right )^{Δcˡ / Rᵛ} \exp \left [ \frac{ℒˡ_{0K}}{Rᵛ} \left (\frac{1}{Tᵗʳ} - \frac{1}{T} \right ) \right ] .
 ```
 
 Consider parameters for liquid water,
@@ -292,6 +292,29 @@ or water ice,
 water_ice = CondensedPhase(latent_heat=2834000, heat_capacity=2108)
 ```
 
+The saturation vapor pressure is
+
+```@example
+using Breeze
+using Breeze.Thermodynamics: saturation_vapor_pressure
+
+thermo = AtmosphereThermodynamics()
+
+T = collect(200:0.1:320)
+pᵛˡ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, thermo.liquid) for Tⁱ in T]
+pᵛⁱ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, thermo.solid) for Tⁱ in T]
+pᵛⁱ⁺[T .> thermo.triple_point_temperature] .= NaN
+
+using CairoMakie
+
+fig = Figure()
+ax = Axis(fig[1, 1], xlabel="Temperature (ᵒK)", ylabel="Saturation vapor pressure pᵛ⁺ (Pa)", yscale = log10, xticks=200:20:320)
+lines!(ax, T, pᵛˡ⁺, label="vapor pressure over liquid")
+lines!(ax, T, pᵛⁱ⁺, linestyle=:dash, label="vapor pressure over ice")
+axislegend(ax, position=:rb)
+fig
+```
+
 The saturation specific humidity is
 
 ```math

src/Thermodynamics/atmosphere_thermodynamics.jl

@@ -79,7 +79,8 @@ Likewise, during deposition, water molecules in the gas phase cluster into ice c
 Arguments
 =========
 - `FT`: Float type to use (defaults to Oceananigans.defaults.FloatType)
-- `latent_heat`: Difference between the internal energy of the gaseous phase at the energy_reference_temperature.
+- `latent_heat`: Difference between the internal energy of the gaseous phase at
+  the `energy_reference_temperature`.
 - `heat_capacity`: Heat capacity of the phase of matter.
 """
 function CondensedPhase(FT = Oceananigans.defaults.FloatType; latent_heat, heat_capacity)

src/Thermodynamics/vapor_saturation.jl

@@ -1,39 +1,37 @@
 """
     saturation_vapor_pressure(T, thermo, phase::CondensedPhase)
 
-Compute the saturation vapor pressure over a liquid surface by integrating
+Compute the saturation vapor pressure ``pᵛ⁺`` over a liquid surface by integrating
 the Clausius-Clapeyron relation,
 
 ```math
-dp/dT = ℒᵛ / (Rᵛ T^2)
+dpᵛ⁺ / dT = pᵛ⁺ ℒᵛ(T) / (Rᵛ T^2)
 ```
 
-which integrates to the expression
+where ``ℒᵛ(T) = ℒᵛ(T=0) + Δcᵖ T``, with ``Δcᵖ ≡ (cᵖᵛ - cᵖˡ)``.
 
-```math
-p(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{aᵛ} \\exp \\left [ bᵛ (1/Tᵗʳ - 1/T) \\right ]
-```
-
-where
+The saturation vapor pressure ``pᵛ⁺`` is obtained after integrating the above from
+the triple point, i.e., ``p(Tᵗʳ) = pᵗʳ`` to get
 
 ```math
-aᵛ ≡ (cᵖᵛ - cᵖˡ) / Rᵛ \\\\
-bᵛ ≡ ℒ₀ / Rᵛ - aᵛ T₀
+pᵛ⁺(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{Δcᵖ / Rᵛ} \\exp \\left [ (1/Tᵗʳ - 1/T) ℒᵛ(T=0) / Rᵛ \\right ]
 ```
 """
 @inline function saturation_vapor_pressure(T, thermo, phase::CondensedPhase)
-    ℒ₀ = phase.latent_heat
+    ℒ₀ = phase.latent_heat # at thermo.energy_reference_temperature
     cᵖˡ = phase.heat_capacity
     T₀ = thermo.energy_reference_temperature
     Tᵗʳ = thermo.triple_point_temperature
     pᵗʳ = thermo.triple_point_pressure
     cᵖᵛ = thermo.vapor.heat_capacity
     Rᵛ = vapor_gas_constant(thermo)
 
-    aᵛ = (cᵖᵛ - cᵖˡ) / Rᵛ
-    bᵛ = ℒ₀ / Rᵛ - aᵛ * T₀
+    Δcᵖ = cᵖᵛ - cᵖˡ
+
+    # latent heat at T = 0 ᵒK assuming temperature-independent specific heats
+    ℒ₀ₖ = ℒ₀ - Δcᵖ * T₀
 
-    return pᵗʳ * (T / Tᵗʳ)^aᵛ * exp(bᵛ * (1/Tᵗʳ - 1/T))
+    return pᵗʳ * (T / Tᵗʳ)^(Δcᵖ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * ℒ₀ₖ / Rᵛ)
 end
 
 # Over a liquid surface