diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index 49a2ddede..427f61dd6 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -2,7 +2,7 @@ module SpeedyWeather # STRUCTURE import Parameters: @with_kw, @unpack - + # NUMERICS import FastGaussQuadrature import AssociatedLegendrePolynomials @@ -10,14 +10,14 @@ module SpeedyWeather 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 @@ -25,7 +25,7 @@ module SpeedyWeather 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 \ No newline at end of file +end diff --git a/src/constants.jl b/src/constants.jl index b2b5781bb..6ea402f91 100644 --- a/src/constants.jl +++ b/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 \ No newline at end of file + drag_strat, RH_thresh_boundary, RH_thresh_range, + RH_thresh_max, humid_relax_time) +end diff --git a/src/default_parameters.jl b/src/default_parameters.jl index fe5c4d74c..1ab6f78dc 100644 --- a/src/default_parameters.jl +++ b/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 \ No newline at end of file +end diff --git a/src/diagnostic_variables.jl b/src/diagnostic_variables.jl index 0d63ef8cf..0105e2bf3 100644 --- a/src/diagnostic_variables.jl +++ b/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,11 +175,11 @@ 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) @@ -152,17 +187,17 @@ function IntermediateVariables(G::GeoSpectral{NF}) where NF 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 - diff --git a/src/humidity.jl b/src/humidity.jl new file mode 100644 index 000000000..fc24f3c81 --- /dev/null +++ b/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 diff --git a/src/large_scale_condensation.jl b/src/large_scale_condensation.jl new file mode 100644 index 000000000..f1211e2d1 --- /dev/null +++ b/src/large_scale_condensation.jl @@ -0,0 +1,48 @@ +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, k_lsc = 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_lsc, humid_tend_lsc, 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 = k_lsc:nlev # Used to be 2:nlev in original Speedy + σₖ = σ_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 + # This formula does not appear in the original Speedy documentation + humid_tend_max = 10σₖ^2 / 3600humid_relax_time + + 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_lsc[i, j, k] = -(humid_grid[i, j, k] - humid_threshold) / humid_relax_time # Formula 22 + temp_tend_lsc[i, j, k] = -alhc / cp * min(humid_tend_lsc[i, j, k], humid_tend_max * pres[i, j]) # Formula 23 + cloud_top[i, j] = min(cloud_top[i, j], k) # Page 7 (last sentence) + else + humid_tend_lsc[i, j, k] = zero(NF) + temp_tend_lsc[i, j, k] = zero(NF) + end + end + end + + # 2. Precipitation due to large-scale condensation + fill!(precip_large_scale, zero(NF)) + for k = k_lsc: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_lsc[i, j, k] # Formula 25 + end + end +end diff --git a/src/parameter_structs.jl b/src/parameter_structs.jl index bc169ea59..3a86a6b03 100644 --- a/src/parameter_structs.jl +++ b/src/parameter_structs.jl @@ -1,3 +1,5 @@ +abstract type Coefficients end + """Coefficients of the generalised logistic function to describe the vertical coordinate. Default coefficients A,K,C,Q,B,M,ν are fitted to the old L31 configuration at ECMWF. See geometry.jl and function vertical_coordinate for more informaiton. @@ -5,12 +7,29 @@ See geometry.jl and function vertical_coordinate for more informaiton. Following the notation of https://en.wikipedia.org/wiki/Generalised_logistic_function (Dec 15 2021). Change default parameters for more/fewer levels in the stratosphere vs troposphere vs boundary layer.""" -@with_kw struct GenLogisticCoefs - A::Real=-0.283 # obtained from a fit in /input_date/vertical_coordinate/vertical_resolution.ipynb - K::Real= 0.871 - C::Real= 0.414 - Q::Real= 6.695 - B::Real=10.336 - M::Real= 0.602 - ν::Real= 5.812 -end \ No newline at end of file +@with_kw struct GenLogisticCoefs{NF<:Real} <: Coefficients + A::NF = -0.283 # obtained from a fit in /input_date/vertical_coordinate/vertical_resolution.ipynb + K::NF = 0.871 + C::NF = 0.414 + Q::NF = 6.695 + B::NF = 10.336 + M::NF = 0.602 + ν::NF = 5.812 +end + +""" +Parameters for computing saturation vapour pressure 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. +""" +@with_kw struct MagnusCoefs{NF<:Real} <: Coefficients + e₀::NF = 6.108 # Saturation vapour pressure at 0°C + T₀::NF = 273.16 # 0°C in Kelvin + T₁::NF = 35.86 + T₂::NF = 7.66 + C₁::NF = 17.269 + C₂::NF = 21.875 +end diff --git a/src/run_speedy.jl b/src/run_speedy.jl index c05055109..e6d03f9d6 100644 --- a/src/run_speedy.jl +++ b/src/run_speedy.jl @@ -11,13 +11,13 @@ function run_speedy(::Type{NF}=Float32; # number format, use Float32 progn_vars,diagn_vars,model_setup = initialize_speedy(NF;kwargs...) # START MODEL INTEGRATION - time_stepping!(progn_vars,diagn_vars,model_setup) + time_stepping!(progn_vars,diagn_vars,model_setup) return progn_vars # return prognostic variables when finished end """ progn_vars, diagn_vars, model_setup = initialize_speedy(NF,kwargs...) - + Initialize the model by returning - `progn_vars`, the initial conditions of the prognostic variables - `diagn_vars`, the preallocated the diagnotic variables (initialised to zero) @@ -43,6 +43,6 @@ function initialize_speedy(::Type{NF}=Float32; # number format, use Float32 prognostic_vars = initial_conditions(M) # initialize prognostic variables diagnostic_vars = DiagnosticVariables(G) # preallocate all diagnostic variables with zeros - + return prognostic_vars, diagnostic_vars, M -end \ No newline at end of file +end diff --git a/src/tendencies_parametrizations.jl b/src/tendencies_parametrizations.jl index b0ac4d466..5b07bfd02 100644 --- a/src/tendencies_parametrizations.jl +++ b/src/tendencies_parametrizations.jl @@ -1,23 +1,16 @@ -struct parametrization_tendencies{T<:AbstractFloat} - -#Just a placeholder - -end - """ -Compute physical parametrization tendencies +Compute physical parametrization tendencies. """ -function parametrization_tendencies!(Prog::PrognosticVariables{NF}, # Prognostic variables - Diag::DiagnosticVariables{NF}, # Diagnostic variables - M - ) where {NF<:AbstractFloat} +function parametrization_tendencies!( + Prog::PrognosticVariables{NF}, + Diag::DiagnosticVariables{NF}, + M, +) where {NF<:AbstractFloat} + get_large_scale_condensation_tendencies!(Diag, M) - - #Calculate + #Calculate #utend: u-wind tendency (gp) #vtend: v-wind tendency (gp) #ttend: temp. tendency (gp) #htend: spec. hum. tendency (gp) - - end diff --git a/test/humidity.jl b/test/humidity.jl new file mode 100644 index 000000000..90a166392 --- /dev/null +++ b/test/humidity.jl @@ -0,0 +1,29 @@ +@testset "humidity.jl" begin + @testset "get_saturation_vapour_pressure!" begin + NF = Float32 + _, diag, model = SpeedyWeather.initialize_speedy(NF) + (;nlon, nlat, nlev) = model.geospectral.geometry + (;sat_vap_pressure) = diag.parametrization_variables + + temp_grid = 200 .+ 150 * rand(NF, nlon, nlat, nlev) # Typical values between 200-350K + + SpeedyWeather.get_saturation_vapour_pressure!(sat_vap_pressure, temp_grid, model) + + @test all(sat_vap_pressure .> 0.0) && all(sat_vap_pressure .< 500.0) + end + + @testset "get_saturation_specific_humidity!" begin + NF = Float32 + _, diag, model = SpeedyWeather.initialize_speedy(NF) + (;nlon, nlat, nlev) = model.geospectral.geometry + (;sat_vap_pressure, sat_spec_humidity) = diag.parametrization_variables + + temp_grid = 200 .+ 150 * rand(NF, nlon, nlat, nlev) # Typical values between 200-350 K + pres_grid = 300 .* 1700 * rand(NF, nlon, nlat) # Typical values between 300-2000 hPa + + SpeedyWeather.get_saturation_specific_humidity!(sat_spec_humidity, sat_vap_pressure, temp_grid, pres_grid, model) + + @test all(isfinite.(sat_spec_humidity)) + @test !any(iszero.(sat_spec_humidity)) + end +end diff --git a/test/large_scale_condensation.jl b/test/large_scale_condensation.jl new file mode 100644 index 000000000..72cce27dc --- /dev/null +++ b/test/large_scale_condensation.jl @@ -0,0 +1,8 @@ +@testset "large_scale_condensation.jl" begin + @testset "get_large_scale_condensation_tendencies!" begin + _, diag, model = SpeedyWeather.initialize_speedy() + + # For now, just check that it runs without errors + SpeedyWeather.get_large_scale_condensation_tendencies!(diag, model) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index c0350f564..2f3b02b7a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,3 +8,7 @@ include("diffusion.jl") include("time_stepping.jl") include("initialize_from_rest.jl") include("run_speedy.jl") + +# PHYSICS +include("humidity.jl") +include("large_scale_condensation.jl")