Parametrization of Large-Scale Condensation (#82)

* Add humidity functions

* Remove prematurely elaborate docstrings

* Switch iteration order

* Add inbounds block

* Fix typo

* Move params into struct

* Add large-scale condensation parametrization

* Make coefficient structs parametric

* Remove @inbounds

* Remove unreachable branch

* Revert renaming of file

* Fix include

* Fix params unpacking

* Add ParametrizationVariables struct

* Add tests

* Tidy up

* Remove Distributions dependency

* Don't export humidity coefficients

* Rename variables

* Add comments

* Rename variables and add docstrings

* Add comment

* Parametrise k for large-scale condensation

* Allocate lsc tendencies explicitly

* Update docstring

* Remove test-specific dependencies

* Define and convert ratio outside loop

* Fix NF in tests

* Update docstring

src/SpeedyWeather.jl

@@ -2,30 +2,30 @@ module SpeedyWeather
 
     # STRUCTURE
     import Parameters: @with_kw, @unpack
-    
+
     # NUMERICS
     import FastGaussQuadrature
     import AssociatedLegendrePolynomials
     import FFTW
     import Primes
     import LinearAlgebra
 
-    # INPUT OUTPUT 
+    # INPUT OUTPUT
     import Dates: Dates, DateTime
     import Printf: @sprintf
     import NetCDF: NetCDF, NcFile, NcDim, NcVar
     import BitInformation: round, round!
     import UnicodePlots
     import ProgressMeter
-    
+
     # EXPORT MAIN INTERFACE TO SPEEDY
     export run_speedy, initialize_speedy
 
     # EXPORT STRUCTS
     export Parameters, GenLogisticCoefs,
         GeoSpectral, Boundaries, Constants, Geometry, SpectralTransform,
         PrognosticVariables, DiagnosticVariables
-    
+
     # EXPORT SPECTRAL FUNCTIONS
     export  spectral, gridded,
         spectral_truncation, spectral_interpolation,
@@ -45,17 +45,21 @@ module SpeedyWeather
     include("boundaries.jl")                # defines Boundaries
     include("horizontal_diffusion.jl")      # defines HorizontalDiffusion
     include("models.jl")                    # defines ModelSetups
-    
+
     include("prognostic_variables.jl")      # defines PrognosticVariables
     include("diagnostic_variables.jl")      # defines DiagnosticVariables
 
-    include("run_speedy.jl")                
+    include("run_speedy.jl")
     include("tendencies_parametrizations.jl")
     include("tendencies_dynamics.jl")
     include("tendencies.jl")
     include("feedback.jl")                  # defines Feedback
-    include("output.jl")                    
+    include("output.jl")
+
+    # PHYSICS
+    include("humidity.jl")
+    include("large_scale_condensation.jl")
 
     include("time_integration.jl")
     include("pretty_printing.jl")
-end
+end

src/constants.jl

@@ -25,6 +25,13 @@ Struct holding the parameters needed at runtime in number format NF.
 
     # DIFFUSION AND DRAG
     drag_strat::NF               # drag [1/s] for zonal wind in the stratosphere
+
+    # PARAMETRIZATIONS
+    # Large-scale condensation (occurs when relative humidity exceeds a given threshold)
+    RH_thresh_boundary::NF  # Relative humidity threshold for boundary layer
+    RH_thresh_range::NF     # Vertical range of relative humidity threshold
+    RH_thresh_max ::NF      # Maximum relative humidity threshold
+    humid_relax_time::NF    # Relaxation time for humidity (hours)
 end
 
 """
@@ -34,11 +41,14 @@ function Constants(P::Parameters)
 
     # PHYSICAL CONSTANTS
     @unpack radius_earth, rotation_earth, gravity, akap, R_gas = P
-    
+
     # TIME INTEGRATION CONSTANTS
     @unpack robert_filter, williams_filter = P
     @unpack trunc, Δt_at_T85, n_days, output_dt = P
 
+    # PARAMETRIZATION CONSTANTS
+    @unpack RH_thresh_boundary, RH_thresh_range, RH_thresh_max, humid_relax_time = P  # Large-scale condensation
+
     Δt_min_at_trunc = Δt_at_T85*(85/trunc)      # scale time step Δt to specified resolution
     Δt      = round(Δt_min_at_trunc*60)         # convert time step Δt from minutes to whole seconds
     Δt_sec  = convert(Int,Δt)                   # encode time step Δt [s] as integer
@@ -54,11 +64,12 @@ function Constants(P::Parameters)
     # SCALING
     Δt_unscaled = Δt        # [s] not scaled
     Δt /= radius_earth      # scale with Earth's radius
-    
+
     # This implies conversion to NF
     return Constants{P.NF}( radius_earth,rotation_earth,gravity,akap,R_gas,
                             Δt,Δt_unscaled,Δt_sec,Δt_hrs,
                             robert_filter,williams_filter,n_timesteps,
                             output_every_n_steps, n_outputsteps,
-                            drag_strat)
-end  
+                            drag_strat, RH_thresh_boundary, RH_thresh_range,
+                            RH_thresh_max, humid_relax_time)
+end

src/default_parameters.jl

@@ -12,7 +12,7 @@ The default values of the keywords define the default model setup.
 
     # RESOLUTION
     trunc::Int=31                       # spectral truncation
-    nlon::Int=roundup_fft(3*trunc+1)    # number of longitudes 
+    nlon::Int=roundup_fft(3*trunc+1)    # number of longitudes
     nlat::Int=nlon÷2                    # number of latitudes
     nlev::Int=8                         # number of vertical levels
 
@@ -50,14 +50,22 @@ The default values of the keywords define the default model setup.
     # DIFFUSION AND DRAG
     diffusion_power::Real=4                 # Power n of Laplacian in horizontal diffusion ∇²ⁿ
     diffusion_time::Real=2.4                # Diffusion time scale [hrs] for temperature and vorticity
-    diffusion_time_div::Real=diffusion_time # Diffusion time scale [hrs] for divergence           
+    diffusion_time_div::Real=diffusion_time # Diffusion time scale [hrs] for divergence
     diffusion_time_strat::Real=12           # Diffusion time scale [hrs] for extra ∇² in the stratosphere
-    damping_time_strat::Real=24*30          # Damping time [hrs] for drag on zonal-mean wind in the stratosphere     
+    damping_time_strat::Real=24*30          # Damping time [hrs] for drag on zonal-mean wind in the stratosphere
 
     # PARAMETRIZATIONS
-    seasonal_cycle::Bool=true   # Seasonal cycle?
-    n_shortwave::Int=3          # Compute shortwave radiation every n steps
-    sppt_on::Bool=false         # Turn on SPPT?
+    seasonal_cycle::Bool=true                  # Seasonal cycle?
+    n_shortwave::Int=3                         # Compute shortwave radiation every n steps
+    sppt_on::Bool=false                        # Turn on SPPT?
+    magnus_coefs::MagnusCoefs = MagnusCoefs()  # For computing saturation vapour pressure
+
+    # Large-Scale Condensation (from table B10)
+    k_lsc::Int = 2                   # Index of atmospheric level at which large-scale condensation begins
+    RH_thresh_boundary::Real = 0.95  # Relative humidity threshold for boundary layer
+    RH_thresh_range::Real = 0.1      # Vertical range of relative humidity threshold
+    RH_thresh_max::Real = 0.9        # Maximum relative humidity threshold
+    humid_relax_time::Real = 4.0     # Relaxation time for humidity (hours)
 
     # TIME STEPPING
     Δt_at_T85::Real=20          # time step in minutes for T85, scale linearly for specified trunc
@@ -88,8 +96,8 @@ The default values of the keywords define the default model setup.
     out_path::String=pwd()      # path to output folder
     output_vars::Vector{String}=["u","v","T","humid","logp0"]
     compression_level::Int=3    # 1=low but fast, 9=high but slow
-    keepbits::Int=10            # mantissa bits to keep for every variable 
+    keepbits::Int=10            # mantissa bits to keep for every variable
 
     # TODO assert not allowed parameter values
     @assert α in [0,0.5,1] "Only semi-implicit α = 0, 0.5 or 1 allowed."
-end
+end

src/diagnostic_variables.jl

@@ -66,9 +66,44 @@ function GridVariables(G::GeoSpectral{NF}) where NF
                         u_grid,v_grid,temp_grid_anomaly)
 end
 
+"""
+Struct holding quantities calculated from the physical parameterisations. All quantities
+are in grid-point space.
+"""
+struct ParametrizationVariables{NF<:AbstractFloat}
+    sat_vap_pressure   ::Array{NF,3}   # Saturation vapour pressure
+    sat_spec_humidity  ::Array{NF,3}   # Saturation specific humidity
+    cloud_top          ::Array{Int,2}  # Cloud-top
+    precip_large_scale ::Array{NF,2}   # Large-scale precipitation
+    humid_tend_lsc     ::Array{NF,3}   # Humidity tendencies due to large-scale condensation
+    temp_tend_lsc      ::Array{NF,3}   # Temperature tendencies due to large-scale condensation
+end
+
+"""
+Generator function for the ParametrizationVariables struct. Initialises with zeros.
+"""
+function ParametrizationVariables(G::GeoSpectral{NF}) where NF
+    @unpack nlon, nlat, nlev = G.geometry
+
+    sat_vap_pressure   = zeros(NF,nlon,nlat,nlev)  # Saturation vapour pressure
+    sat_spec_humidity  = zeros(NF,nlon,nlat,nlev)  # Saturation specific humidity
+    cloud_top          = zeros(Int,nlon,nlat)      # Cloud-top
+    precip_large_scale = zeros(NF,nlon,nlat)       # Large-scale precipitation
+    humid_tend_lsc     = zeros(NF,nlon,nlat,nlev)  # Humidity tendencies due to large-scale condensation
+    temp_tend_lsc      = zeros(NF,nlon,nlat,nlev)  # Temperature tendencies due to large-scale condensation
+
+    return ParametrizationVariables(sat_vap_pressure,
+                                    sat_spec_humidity,
+                                    cloud_top,
+                                    precip_large_scale,
+                                    humid_tend_lsc,
+                                    temp_tend_lsc,
+                                    )
+end
+
 """Struct holding intermediate quantities that are used and shared when calculating tendencies"""
 struct IntermediateVariables{NF<:AbstractFloat}
-    
+
     ### VORTICITY INVERSION
     velocity_potential  ::Array{Complex{NF},3}
     stream_function     ::Array{Complex{NF},3}
@@ -100,16 +135,16 @@ struct IntermediateVariables{NF<:AbstractFloat}
     puv        ::Array{NF,3} #(ug -umean)*px + (vg -vmean)*py
 
     ###------Defined in zonal_wind_tendency!()
-    sigma_u ::Array{NF,3}  #some quantity used for later calculations 
+    sigma_u ::Array{NF,3}  #some quantity used for later calculations
 
     ###------Defined in vor_div_tendency_and_corrections!()
     L2_velocity_complex ::Array{Complex{NF},2} # -laplacian(0.5*(u**2+v**2))
 
-    ###-----Defined in tendencies.jl/get_spectral_tendencies!() 
+    ###-----Defined in tendencies.jl/get_spectral_tendencies!()
     vertical_mean_divergence ::Array{Complex{NF},2}
     sigdtc ::Array{Complex{NF},3} # what is this quantity, physically?
     dumk ::Array{Complex{NF},3} #ditto
-    spectral_geopotential ::Array{Complex{NF},3} #This should probably go elsewhere 
+    spectral_geopotential ::Array{Complex{NF},3} #This should probably go elsewhere
 end
 
 """
@@ -125,7 +160,7 @@ function IntermediateVariables(G::GeoSpectral{NF}) where NF
     stream_function = zeros(Complex{NF},lmax+1,mmax+1,nlev)
     coslat_u = zeros(Complex{NF},lmax+2,mmax+1,nlev)
     coslat_v = zeros(Complex{NF},lmax+2,mmax+1,nlev)
-    
+
     # VORTICITY ADVECTION
     uω_grid  = zeros(NF,nlon,nlat,nlev)
     vω_grid  = zeros(NF,nlon,nlat,nlev)
@@ -140,29 +175,29 @@ function IntermediateVariables(G::GeoSpectral{NF}) where NF
 
     # one more l for recursion in meridional gradients
     # X,Y gradient of the surface pressure in spectral space
-    pres_surf_gradient_spectral_x = zeros(Complex{NF},lmax+2,mmax+1)  
+    pres_surf_gradient_spectral_x = zeros(Complex{NF},lmax+2,mmax+1)
     pres_surf_gradient_spectral_y = zeros(Complex{NF},lmax+2,mmax+1)
 
     # X,Y gradient of the surface pressure in grid space
-    pres_surf_gradient_grid_x = zeros(NF,nlon,nlat)  
+    pres_surf_gradient_grid_x = zeros(NF,nlon,nlat)
     pres_surf_gradient_grid_y = zeros(NF,nlon,nlat)
 
     sigma_tend  = zeros(NF,nlon,nlat,nlev+1)
     sigma_m     = zeros(NF,nlon,nlat,nlev+1)
     puv         = zeros(NF,nlon,nlat,nlev)
     sigma_u     = zeros(NF,nlon,nlat,nlev+1)
 
-    L2_velocity_complex         = zeros(Complex{NF},lmax+2,mmax+1)  
+    L2_velocity_complex         = zeros(Complex{NF},lmax+2,mmax+1)
 
-    vertical_mean_divergence    = zeros(Complex{NF},lmax+2,mmax+1)  
-    sigdtc                      = zeros(Complex{NF},lmax+2,mmax+1,nlev+1)  
-    dumk                        = zeros(Complex{NF},lmax+2,mmax+1,nlev+1)  
+    vertical_mean_divergence    = zeros(Complex{NF},lmax+2,mmax+1)
+    sigdtc                      = zeros(Complex{NF},lmax+2,mmax+1,nlev+1)
+    dumk                        = zeros(Complex{NF},lmax+2,mmax+1,nlev+1)
     spectral_geopotential       = zeros(Complex{NF},lmax+2,mmax+1,nlev)
 
     return IntermediateVariables(   velocity_potential, stream_function,
                                     coslat_u, coslat_v,
                                     uω_grid,vω_grid,
-                                    uω,vω,∂uω_∂lon,∂vω_∂lat,         
+                                    uω,vω,∂uω_∂lon,∂vω_∂lat,
                                     u_mean,v_mean,div_mean,
                                     pres_surf_gradient_spectral_x,pres_surf_gradient_spectral_y,
                                     pres_surf_gradient_grid_x,pres_surf_gradient_grid_y,
@@ -172,18 +207,21 @@ end
 
 """Struct holding the diagnostic variables."""
 struct DiagnosticVariables{NF<:AbstractFloat}
-    tendencies             ::Tendencies{NF}
-    grid_variables         ::GridVariables{NF}
-    intermediate_variables ::IntermediateVariables{NF}
+    tendencies                ::Tendencies{NF}
+    grid_variables            ::GridVariables{NF}
+    intermediate_variables    ::IntermediateVariables{NF}
+    parametrization_variables ::ParametrizationVariables{NF}
 end
 
 """Generator function for Diagnostic Variables """
 function DiagnosticVariables(G::GeoSpectral)
-    tendencies             = Tendencies(G)
-    grid_variables         = GridVariables(G)
-    intermediate_variables = IntermediateVariables(G)
+    tendencies                = Tendencies(G)
+    grid_variables            = GridVariables(G)
+    intermediate_variables    = IntermediateVariables(G)
+    parametrization_variables = ParametrizationVariables(G)
     return DiagnosticVariables( tendencies,
                                 grid_variables,
-                                intermediate_variables)
+                                intermediate_variables,
+                                parametrization_variables,
+                                )
 end
-

src/humidity.jl

@@ -0,0 +1,59 @@
+"""
+Compute the saturation vapour pressure as a function of temperature using the
+    August-Roche-Magnus formula,
+
+    eᵢ(T) = e₀ * exp(Cᵢ * (T - T₀) / (T - Tᵢ)),
+
+    where T is in Kelvin and i = 1,2 for saturation with respect to water and ice,
+    respectively.
+"""
+function get_saturation_vapour_pressure!(
+    sat_vap_pressure ::Array{NF,3},
+    temp_grid        ::Array{NF,3},
+    model            ::ModelSetup,
+    ) where {NF<:AbstractFloat}
+    @unpack nlon, nlat, nlev = model.geospectral.geometry
+    @unpack e₀, T₀, C₁, C₂, T₁, T₂ = model.parameters.magnus_coefs
+
+    for k = 1:nlev, j = 1:nlat, i = 1:nlon
+        if temp_grid[i, j, k] > T₀
+            # Saturation vapour pressure over water
+            sat_vap_pressure[i, j, k] = e₀ * exp(C₁ * (temp_grid[i, j, k] - T₀) / (temp_grid[i, j, k] - T₁))
+        else
+            # Saturation vapour pressure over ice
+            sat_vap_pressure[i, j, k] = e₀ * exp(C₂ * (temp_grid[i, j, k] - T₀) / (temp_grid[i, j, k] - T₂))
+        end
+    end
+
+    return nothing
+end
+
+"""
+Compute the saturation specific humidity according to the formula,
+
+    0.622 * e / (p - (1 - 0.622) * e),
+
+    where e is the saturation vapour pressure, p is the pressure, and 0.622 is the ratio of
+    the molecular weight of water to dry air.
+"""
+function get_saturation_specific_humidity!(
+    sat_spec_humidity ::Array{NF,3},
+    sat_vap_pressure  ::Array{NF,3},
+    temp_grid         ::Array{NF,3},
+    pres              ::Array{NF,2},
+    model             ::ModelSetup,
+) where {NF<:AbstractFloat}
+    @unpack nlon, nlat, nlev, σ_levels_full = model.geospectral.geometry
+
+    mol_ratio = convert(NF, 0.622)
+    one_minus_mol_ratio = convert(NF, 1 - mol_ratio)
+
+    get_saturation_vapour_pressure!(sat_vap_pressure, temp_grid, model)
+
+    for k = 1:nlev, j = 1:nlat, i = 1:nlon
+        sat_spec_humidity[i, j, k] = mol_ratio * sat_vap_pressure[i, j, k] /
+            (pres[i, j] * σ_levels_full[k] - one_minus_mol_ratio * sat_vap_pressure[i, j, k])
+    end
+
+    return nothing
+end