Parametrization of Large-Scale Condensation

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,
@@ -48,13 +48,17 @@ module SpeedyWeather
     include("horizontal_diffusion.jl")      # defines HorizontalDiffusion
 
     include("models.jl")                    # defines ModelSetups
-    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

@@ -9,7 +9,7 @@ With keywords such that default values can be changed at creation.
 
     # 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
 
@@ -47,14 +47,21 @@ With keywords such that default values can be changed at creation.
     # 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)
+    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
@@ -85,8 +92,8 @@ With keywords such that default values can be changed at creation.
     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         ::Array{NF,3}   # Humidity tendencies due to physics
+    temp_tend          ::Array{NF,3}   # Temperature tendencies due to physics
+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         = zeros(NF,nlon,nlat,nlev)  # Humidity tendencies due to physics
+    temp_tend          = zeros(NF,nlon,nlat,nlev)  # Temperature tendencies due to physics
+
+    return ParametrizationVariables(sat_vap_pressure,
+                                    sat_spec_humidity,
+                                    cloud_top,
+                                    precip_large_scale,
+                                    humid_tend,
+                                    temp_tend,
+                                    )
+end
+
 """Struct holding intermediate quantities that are used and shared when calculating tendencies"""
 struct IntermediateVariables{NF<:AbstractFloat}
-    
+
     ### VORTICITY INVERSION
     stream_function ::Array{Complex{NF},3}
     coslat_u        ::Array{Complex{NF},3}
@@ -99,16 +134,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
 
 """
@@ -123,7 +158,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)
@@ -138,28 +173,28 @@ 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(   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,
@@ -169,18 +204,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,55 @@
+"""
+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 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
+
+    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] = 0.622 * sat_vap_pressure[i, j, k] /
+            (pres[i, j] * σ_levels_full[k] - (1 - 0.622) * sat_vap_pressure[i, j, k])
+    end
+
+    return nothing
+end

src/large_scale_condensation.jl

@@ -0,0 +1,43 @@
+function get_large_scale_condensation_tendencies!(
+    Diag::DiagnosticVariables{NF},
+    M   ::ModelSetup,
+) where {NF<:AbstractFloat}
+    @unpack gravity, RH_thresh_max, RH_thresh_range, RH_thresh_boundary, humid_relax_time = M.constants
+    @unpack cp, alhc = M.parameters
+    @unpack nlon, nlat, nlev, σ_levels_full, σ_levels_thick = M.geospectral.geometry
+    @unpack temp_grid, humid_grid, pres_surf_grid = Diag.grid_variables
+    @unpack temp_tend, humid_tend, sat_spec_humidity, sat_vap_pressure, cloud_top, precip_large_scale = Diag.parametrization_variables
+
+    pres = exp.(pres_surf_grid)  # Normalised surface pressure - TODO(alistair): check pressure units throughout
+
+    get_saturation_specific_humidity!(sat_spec_humidity, sat_vap_pressure, temp_grid, pres, M)
+
+    # 1. Tendencies of humidity and temperature due to large-scale condensation
+    for k = 2:nlev
+        σₖ = σ_levels_full[k]
+        RH_threshold = RH_thresh_max + RH_thresh_range * (σₖ^2 - 1)  # Relative humidity threshold for condensation (Formula 24)
+        if k == nlev
+            RH_threshold = max(RH_threshold, RH_thresh_boundary)
+        end
+
+        # Impose a maximum heating rate to avoid grid-point storm instability
+        humid_tend_max = 10σₖ^2 / 3600humid_relax_time  # This formula does not appear in the documentation
+
+        for j = 1:nlat, i = 1:nlon
+            humid_threshold = RH_threshold * sat_spec_humidity[i, j, k]  # Specific humidity threshold for condensation
+            if humid_grid[i, j, k] > humid_threshold
+                humid_tend[i, j, k] += -(humid_grid[i, j, k] - humid_threshold) / humid_relax_time  # Formula 22
+                temp_tend[i, j, k] += -alhc / cp * min(humid_tend[i, j, k], humid_tend_max * pres[i, j])  # Formula 23
+                cloud_top[i, j] = min(k, cloud_top[i, j])  # Page 7 (last sentence)
+            end
+        end
+    end
+
+    # 2. Precipitation due to large-scale condensation
+    for k = 2:nlev
+        Δpₖ = pres * σ_levels_thick[k]  # Formula 4
+        for j = 1:nlat, i = 1:nlon
+            precip_large_scale[i, j] += -1 / gravity * Δpₖ[i, j] * humid_tend[i, j, k]  # Formula 25
+        end
+    end
+end