From 7b6a6d4214e08ef799a88fbed8346005fc674ce4 Mon Sep 17 00:00:00 2001 From: Milan Date: Tue, 30 May 2023 20:41:58 -0400 Subject: [PATCH 1/8] SpeedyCondensation implemented --- src/SpeedyWeather.jl | 6 +- src/abstract_types.jl | 1 + src/dynamics/atmospheres.jl | 5 +- src/dynamics/constants.jl | 2 +- src/dynamics/horizontal_diffusion.jl | 2 +- src/dynamics/implicit.jl | 4 +- src/dynamics/initial_conditions.jl | 96 +++++++++-------- src/dynamics/models.jl | 10 +- src/dynamics/orography.jl | 4 +- src/dynamics/planets.jl | 2 +- src/dynamics/spectral_grid.jl | 4 +- src/dynamics/time_integration.jl | 4 +- src/dynamics/vertical_coordinates.jl | 8 +- src/output/output.jl | 14 ++- src/physics/boundary_layer.jl | 2 +- src/physics/column_variables.jl | 8 +- src/physics/convection.jl | 22 ---- src/physics/large_scale_condensation.jl | 138 ++++++++++++++++++------ src/physics/shortwave_radiation.jl | 2 +- src/physics/temperature_relaxation.jl | 4 +- src/physics/tendencies.jl | 2 +- src/physics/thermodynamics.jl | 4 +- src/physics/vertical_diffusion.jl | 4 +- 23 files changed, 212 insertions(+), 136 deletions(-) diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index d68e33222..c4ba991b7 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -1,7 +1,6 @@ module SpeedyWeather # STRUCTURE -import Parameters: @with_kw using DocStringExtensions # NUMERICS @@ -95,6 +94,9 @@ export NoBoundaryLayer, export NoVerticalDiffusion, VerticalLaplacian +# Large scale condensation +export SpeedyCondensation + # EXPORT STRUCTS export DynamicsConstants, SpectralTransform, @@ -170,7 +172,7 @@ include("dynamics/models.jl") # # PHYSICS # include("physics/convection.jl") -# include("physics/large_scale_condensation.jl") +include("physics/large_scale_condensation.jl") # include("physics/longwave_radiation.jl") # include("physics/shortwave_radiation.jl") diff --git a/src/abstract_types.jl b/src/abstract_types.jl index 55a4208b2..5f7dea0c6 100644 --- a/src/abstract_types.jl +++ b/src/abstract_types.jl @@ -32,6 +32,7 @@ abstract type BoundaryLayerDrag{NF} <: AbstractParameterization{NF} end abstract type TemperatureRelaxation{NF} <: AbstractParameterization{NF} end abstract type VerticalDiffusion{NF} <: AbstractParameterization{NF} end abstract type AbstractThermodynamics{NF} <: AbstractParameterization{NF} end +abstract type AbstractCondensation{NF} <: AbstractParameterization{NF} end # INPUT/OUTPUT abstract type AbstractFeedback end diff --git a/src/dynamics/atmospheres.jl b/src/dynamics/atmospheres.jl index 5c543ef3c..00ee3d692 100644 --- a/src/dynamics/atmospheres.jl +++ b/src/dynamics/atmospheres.jl @@ -4,7 +4,7 @@ Create a struct `EarthAtmosphere<:AbstractPlanet`, with the following physical/c characteristics. Note that `radius` is not part of it as this should be chosen in `SpectralGrid`. Keyword arguments are $(TYPEDFIELDS)""" -@with_kw struct EarthAtmosphere <: AbstractAtmosphere +Base.@kwdef struct EarthAtmosphere <: AbstractAtmosphere # ATMOSPHERE "molar mass of dry air [g/mol]" mol_mass_dry_air::Float64 = 28.9649 @@ -50,6 +50,9 @@ $(TYPEDFIELDS)""" "start of the stratosphere in sigma coordinates" σ_tropopause::Float64 = 0.2 + "top of the planetary boundary layer in sigma coordinates" + σ_boundary_layer::Float64 = 0.95 + "scale height for pressure [km]" scale_height::Float64 = 7.5 diff --git a/src/dynamics/constants.jl b/src/dynamics/constants.jl index 80da11c68..2a17edab7 100644 --- a/src/dynamics/constants.jl +++ b/src/dynamics/constants.jl @@ -1,7 +1,7 @@ """ Struct holding constants needed at runtime for the dynamical core in number format NF. $(TYPEDFIELDS)""" -@with_kw struct DynamicsConstants{NF<:AbstractFloat} <: AbstractDynamicsConstants{NF} +Base.@kwdef struct DynamicsConstants{NF<:AbstractFloat} <: AbstractDynamicsConstants{NF} # PHYSICAL CONSTANTS "Radius of Planet [m]" radius::NF diff --git a/src/dynamics/horizontal_diffusion.jl b/src/dynamics/horizontal_diffusion.jl index 1ea907bc8..192e72f82 100644 --- a/src/dynamics/horizontal_diffusion.jl +++ b/src/dynamics/horizontal_diffusion.jl @@ -7,7 +7,7 @@ layers. Furthermore the power can be decreased above the `tapering_σ` to `power_stratosphere` (default 2). For Barotropic, ShallowWater, the default non-adaptive constant-time scale hyper diffusion is used. Options are $(TYPEDFIELDS)""" -@with_kw struct HyperDiffusion{NF} <: HorizontalDiffusion{NF} +Base.@kwdef struct HyperDiffusion{NF} <: HorizontalDiffusion{NF} # DIMENSIONS "spectral resolution" trunc::Int diff --git a/src/dynamics/implicit.jl b/src/dynamics/implicit.jl index 1dc3fe7c0..21ce0b651 100644 --- a/src/dynamics/implicit.jl +++ b/src/dynamics/implicit.jl @@ -8,7 +8,7 @@ initialize!(I::NoImplicit,dt::Real,::DiagnosticVariables,::ModelSetup) = nothing Struct that holds various precomputed arrays for the semi-implicit correction to prevent gravity waves from amplifying in the shallow water model. $(TYPEDFIELDS)""" -@with_kw struct ImplicitShallowWater{NF<:AbstractFloat} <: AbstractImplicit{NF} +Base.@kwdef struct ImplicitShallowWater{NF<:AbstractFloat} <: AbstractImplicit{NF} # DIMENSIONS trunc::Int @@ -130,7 +130,7 @@ end Struct that holds various precomputed arrays for the semi-implicit correction to prevent gravity waves from amplifying in the primitive equation model. $(TYPEDFIELDS)""" -@with_kw struct ImplicitPrimitiveEq{NF<:AbstractFloat} <: AbstractImplicit{NF} +Base.@kwdef struct ImplicitPrimitiveEq{NF<:AbstractFloat} <: AbstractImplicit{NF} # DIMENSIONS "spectral resolution" diff --git a/src/dynamics/initial_conditions.jl b/src/dynamics/initial_conditions.jl index c0c08f701..c2ed983bd 100644 --- a/src/dynamics/initial_conditions.jl +++ b/src/dynamics/initial_conditions.jl @@ -20,7 +20,7 @@ function allocate(::Type{PrognosticVariables},spectral_grid::SpectralGrid{Model} return zeros(PrognosticVariables{NF},Model,trunc,nlev) end -@with_kw struct StartFromRest <: InitialConditions +Base.@kwdef struct StartFromRest <: InitialConditions pressure_on_orography::Bool = false end @@ -40,7 +40,7 @@ end """Start with random vorticity as initial conditions $(TYPEDFIELDS)""" -@with_kw struct StartWithVorticity <: InitialConditions +Base.@kwdef struct StartWithVorticity <: InitialConditions "Power law the vorticity should be spectrally distributed by" power_law::Float64 = -3 @@ -76,7 +76,7 @@ $(TYPEDSIGNATURES) Create a struct that contains all parameters for the Galewsky et al, 2004 zonal jet intitial conditions for the shallow water model. Default values as in Galewsky. $(TYPEDFIELDS)""" -@with_kw struct ZonalJet <: InitialConditions +Base.@kwdef struct ZonalJet <: InitialConditions "jet latitude [˚N]" latitude::Float64 = 45 @@ -196,7 +196,7 @@ $(TYPEDSIGNATURES) Create a struct that contains all parameters for the Jablonowski and Williamson, 2006 intitial conditions for the primitive equation model. Default values as in Jablonowski. $(TYPEDFIELDS)""" -@with_kw struct ZonalWind <: InitialConditions +Base.@kwdef struct ZonalWind <: InitialConditions "conversion from σ to Jablonowski's ηᵥ-coordinates" η₀::Float64 = 0.252 @@ -337,7 +337,12 @@ function initial_conditions!( progn::PrognosticVariables{NF}, # PRESSURE (constant everywhere) lnp₀ = log(pres_ref*100) # logarithm of reference surface pressure, *100 for [hPa] to [Pa] progn.surface.timesteps[1].pres[1] = norm_sphere*lnp₀ - pressure_on_orography && pressure_on_orography!(progn,model) + + lnpₛ = ones(Grid{NF},nlat_half) + lnpₛ .= pressure_on_orography ? pressure_on_orography!(progn,model) : lnp₀ + + # HUMIDITY + initialize_humidity!(progn,lnpₛ,model) end """ @@ -345,7 +350,7 @@ Restart from a previous SpeedyWeather.jl simulation via the restart file restart Applies interpolation in the horizontal but not in the vertical. restart.jld2 is identified by $(TYPEDFIELDS)""" -@with_kw struct StartFromFile <: InitialConditions +Base.@kwdef struct StartFromFile <: InitialConditions "path for restart file" path::String = pwd() @@ -447,40 +452,45 @@ function pressure_on_orography!(progn::PrognosticVariables, return lnp_grid # return grid for use in initialize_humidity! end -# """Initialize specific humidity in spectral space.""" -# function initialize_humidity!( humid::AbstractArray{Complex{NF},3},# spectral specific humidity -# pres_surf_grid::AbstractMatrix, # log of surf pressure (grid space) -# P::Parameters, # Parameters struct -# G::GeoSpectral) where NF # Geospectral struct - -# lmax,mmax,nlev = size(humid) # of size lmax+1, mmax+1, nlev -# lmax, mmax = lmax-1, mmax-1 # hence correct with -1 for 0-based l,m -# (; nlon, nlat, n_stratosphere_levels ) = P -# (; σ_levels_full ) = G.geometry - -# # reference saturation water vapour pressure [Pa] -# # relative humidity reference [1] -# (; water_pres_ref, relhumid_ref ) = P -# humid_ref = relhumid_ref*0.622*water_pres_ref # reference specific humidity [Pa] - -# # scale height [km], scale height for spec humidity [km] -# (; scale_height, scale_height_humid ) = P -# scale_height_ratio = scale_height/scale_height_humid # ratio of scale heights [1] - -# # Specific humidity at the surface (grid space) -# humid_surf_grid = humid_ref*exp.(scale_height_ratio*pres_surf_grid) -# humid_surf = spectral(humid_surf_grid,one_more_l=true) -# # spectral_truncation!(humid_surf,P.trunc) # set the lmax+1 row to zero - -# # stratospheric humidity zero -# fill!(view(humid,:,:,1:n_stratosphere_levels),0) - -# # Specific humidity at tropospheric levels -# for k in n_stratosphere_levels+1:nlev -# for m in 1:mmax+1 -# for l in m:lmax+1 -# humid[l,m,k] = humid_surf[l,m]*σ_levels_full[k]^scale_height_ratio -# end -# end -# end -# end \ No newline at end of file +function initialize_humidity!( progn::PrognosticVariables, + pres_surf_grid::AbstractGrid, + model::PrimitiveDry) + return nothing +end + +function initialize_humidity!( progn::PrognosticVariables, + pres_surf_grid::AbstractGrid, + model::PrimitiveWet) + + # reference saturation water vapour pressure [Pa] + # relative humidity reference [1] + (;water_pres_ref, relhumid_ref, R_dry, R_vapour) = model.atmosphere + gas_ratio = R_dry/R_vapour + humid_ref = relhumid_ref*gas_ratio*water_pres_ref # reference specific humidity [Pa] + + # ratio of scale heights [1], scale height [km], scale height for spec humidity [km] + (;scale_height, scale_height_humid, σ_tropopause) = model.atmosphere + scale_height_ratio = scale_height/scale_height_humid + + (;nlev, σ_levels_full) = model.geometry + n_stratosphere_levels = findfirst(σ->σ>=σ_tropopause,σ_levels_full) + + # Specific humidity at the surface (grid space) + humid_surf_grid = zero(pres_surf_grid) + @. humid_surf_grid = humid_ref*exp(scale_height_ratio*pres_surf_grid) + + humid_surf = spectral(humid_surf_grid,model.spectral_transform) + spectral_truncation!(humid_surf) + + q0 = 3e-3 + q1 = -1e-3 + + # Specific humidity at tropospheric levels (stratospheric humidity remains zero) + for k in n_stratosphere_levels+1:nlev + # for lm in eachharmonic(humid_surf,progn.layers[1].timesteps[1].humid) + # progn.layers[k].timesteps[1].humid[lm] = humid_surf[lm]*σ_levels_full[k]^scale_height_ratio + # end + progn.layers[k].timesteps[1].humid[1] = model.spectral_transform.norm_sphere*q0*σ_levels_full[k]^scale_height_ratio + progn.layers[k].timesteps[1].humid[3] = model.spectral_transform.norm_sphere*q1*σ_levels_full[k]^scale_height_ratio + end +end \ No newline at end of file diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index 23494c7ce..3c16dbae4 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -19,7 +19,7 @@ $(SIGNATURES) The BarotropicModel struct holds all other structs that contain precalculated constants, whether scalars or arrays that do not change throughout model integration. $(TYPEDFIELDS)""" -@with_kw struct BarotropicModel{NF<:AbstractFloat, D<:AbstractDevice} <: Barotropic +Base.@kwdef struct BarotropicModel{NF<:AbstractFloat, D<:AbstractDevice} <: Barotropic "dictates resolution for many other components" spectral_grid::SpectralGrid = SpectralGrid() @@ -70,7 +70,7 @@ $(SIGNATURES) The ShallowWaterModel struct holds all other structs that contain precalculated constants, whether scalars or arrays that do not change throughout model integration. $(TYPEDFIELDS)""" -@with_kw struct ShallowWaterModel{NF<:AbstractFloat, D<:AbstractDevice} <: ShallowWater +Base.@kwdef struct ShallowWaterModel{NF<:AbstractFloat, D<:AbstractDevice} <: ShallowWater "dictates resolution for many other components" spectral_grid::SpectralGrid = SpectralGrid() @@ -125,7 +125,7 @@ $(SIGNATURES) The PrimitiveDryModel struct holds all other structs that contain precalculated constants, whether scalars or arrays that do not change throughout model integration. $(TYPEDFIELDS)""" -@with_kw struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveDry +Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveDry "dictates resolution for many other components" spectral_grid::SpectralGrid = SpectralGrid() @@ -192,7 +192,7 @@ $(SIGNATURES) The PrimitiveDryModel struct holds all other structs that contain precalculated constants, whether scalars or arrays that do not change throughout model integration. $(TYPEDFIELDS)""" -@with_kw struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveWet +Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveWet "dictates resolution for many other components" spectral_grid::SpectralGrid = SpectralGrid() @@ -209,6 +209,7 @@ $(TYPEDFIELDS)""" boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid) temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid) static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid) + large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid) # vertical_diffusion::VerticalDiffusion{NF} = VerticalLaplacian(spectral_grid) # NUMERICS @@ -248,6 +249,7 @@ function initialize!(model::PrimitiveWet) initialize!(model.boundary_layer_drag,model) initialize!(model.temperature_relaxation,model) initialize!(model.static_energy_diffusion,model) + initialize!(model.large_scale_condensation,model) # initialize!(model.vertical_diffusion,model) prognostic_variables = initial_conditions(model) diff --git a/src/dynamics/orography.jl b/src/dynamics/orography.jl index 93e9ae9ce..da2033c76 100644 --- a/src/dynamics/orography.jl +++ b/src/dynamics/orography.jl @@ -28,7 +28,7 @@ initialize!(::NoOrography,::AbstractPlanet,::SpectralTransform,::Geometry) = not """Zonal ridge orography after Jablonowski and Williamson, 2006. $(TYPEDFIELDS)""" -@with_kw struct ZonalRidge{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} <: AbstractOrography{NF,Grid} +Base.@kwdef struct ZonalRidge{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} <: AbstractOrography{NF,Grid} "conversion from σ to Jablonowski's ηᵥ-coordinates" η₀::Float64 = 0.252 @@ -100,7 +100,7 @@ end """Earth's orography read from file, with smoothing. $(TYPEDFIELDS)""" -@with_kw struct EarthOrography{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} <: AbstractOrography{NF,Grid} +Base.@kwdef struct EarthOrography{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} <: AbstractOrography{NF,Grid} # OPTIONS "path to the folder containing the orography file, pkg path default" diff --git a/src/dynamics/planets.jl b/src/dynamics/planets.jl index 124d99a3f..353a78800 100644 --- a/src/dynamics/planets.jl +++ b/src/dynamics/planets.jl @@ -5,7 +5,7 @@ characteristics. Note that `radius` is not part of it as this should be chosen in `SpectralGrid`. Keyword arguments are $(TYPEDFIELDS) """ -@with_kw struct Earth <: AbstractPlanet +Base.@kwdef struct Earth <: AbstractPlanet "angular frequency of Earth's rotation [rad/s]" rotation::Float64 = 7.29e-5 diff --git a/src/dynamics/spectral_grid.jl b/src/dynamics/spectral_grid.jl index eb48ee4c2..8c7491839 100644 --- a/src/dynamics/spectral_grid.jl +++ b/src/dynamics/spectral_grid.jl @@ -9,7 +9,7 @@ $(TYPEDFIELDS) `nlat_half` and `npoints` should not be chosen but are derived from `trunc`, `Grid` and `dealiasing`.""" -@with_kw struct SpectralGrid{Model<:ModelSetup} +Base.@kwdef struct SpectralGrid{Model<:ModelSetup} "number format used throughout the model" NF::Type{<:AbstractFloat} = DEFAULT_NF @@ -76,7 +76,7 @@ Construct Geometry struct containing parameters and arrays describing an iso-lat and the vertical levels. Pass on `SpectralGrid` to calculate the following fields $(TYPEDFIELDS) """ -@with_kw struct Geometry{NF<:AbstractFloat} <: AbstractGeometry{NF} # NF: Number format +Base.@kwdef struct Geometry{NF<:AbstractFloat} <: AbstractGeometry{NF} # NF: Number format "SpectralGrid that defines spectral and grid resolution" spectral_grid::SpectralGrid diff --git a/src/dynamics/time_integration.jl b/src/dynamics/time_integration.jl index 48a29bb7a..b3a933f87 100644 --- a/src/dynamics/time_integration.jl +++ b/src/dynamics/time_integration.jl @@ -2,7 +2,7 @@ Clock struct keeps track of the model time, how many days to integrate for and how many time steps this takes $(TYPEDFIELDS).""" -@with_kw mutable struct Clock +Base.@kwdef mutable struct Clock "current model time" time::DateTime = DateTime(2000,1,1) @@ -27,7 +27,7 @@ end Leapfrog time stepping defined by the following fields $(TYPEDFIELDS) """ -@with_kw struct Leapfrog{NF} <: TimeStepper{NF} +Base.@kwdef struct Leapfrog{NF} <: TimeStepper{NF} # DIMENSIONS "spectral resolution (max degree of spherical harmonics)" diff --git a/src/dynamics/vertical_coordinates.jl b/src/dynamics/vertical_coordinates.jl index ba6d73f4c..f38e2ae06 100644 --- a/src/dynamics/vertical_coordinates.jl +++ b/src/dynamics/vertical_coordinates.jl @@ -1,8 +1,8 @@ -@with_kw struct NoVerticalCoordinates <: VerticalCoordinates +Base.@kwdef struct NoVerticalCoordinates <: VerticalCoordinates nlev::Int = 1 end -@with_kw struct SigmaCoordinates <: VerticalCoordinates +Base.@kwdef struct SigmaCoordinates <: VerticalCoordinates nlev::Int = 8 σ_half::Vector{Float64} = default_sigma_coordinates(nlev) @@ -19,7 +19,7 @@ Default coefficients A,K,C,Q,B,M,ν are fitted to the old L31 configuration at E Following the notation of [https://en.wikipedia.org/wiki/Generalised_logistic_function](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 +Base.@kwdef struct GenLogisticCoefs A::Float64 = -0.283 # obtained from a fit in /input_data/vertical_coordinate/vertical_resolution.ipynb K::Float64 = 0.871 C::Float64 = 0.414 @@ -63,7 +63,7 @@ function sigma_okay(nlev::Integer,σ_half::AbstractVector) end #TODO -@with_kw struct SigmaPressureCoordinates <: VerticalCoordinates +Base.@kwdef struct SigmaPressureCoordinates <: VerticalCoordinates nlev::Int = 8 A::Vector{Float64} = default_hybrid_coordinates(:A,nlev) B::Vector{Float64} = default_hybrid_coordinates(:B,nlev) diff --git a/src/output/output.jl b/src/output/output.jl index faf60a4a7..e06ed28dc 100644 --- a/src/output/output.jl +++ b/src/output/output.jl @@ -2,7 +2,7 @@ Number of mantissa bits to keep for each prognostic variable when compressed for netCDF and .jld2 data output. $(TYPEDFIELDS)""" -@with_kw struct Keepbits +Base.@kwdef struct Keepbits u::Int = 7 v::Int = 7 vor::Int = 5 @@ -37,7 +37,7 @@ Base.@kwdef mutable struct OutputWriter{NF<:Union{Float32,Float64},Model<:ModelS # FILE OPTIONS output::Bool = false # output to netCDF? path::String = pwd() # path to output folder - id::Union{String,Int} = "" # run identification number/string + id::String = "0001" # run identification string "????" run_path::String = "" # will be determined in initalize! filename::String = "output.nc" # name of the output netcdf file write_restart::Bool = true # also write restart file if output==true? @@ -201,7 +201,7 @@ function initialize!( # GET RUN ID, CREATE FOLDER # get new id only if not already specified - output.id = output.id == "" ? get_run_id(output.path) : output.id + output.id = get_run_id(output.path,output.id) output.run_path = create_output_folder(output.path,output.id) feedback.id = output.id # synchronize with feedback struct @@ -250,8 +250,12 @@ Checks existing `run_????` folders in `path` to determine a 4-digit `id` number by counting up. E.g. if folder run_0001 exists it will return the string "0002". Does not create a folder for the returned run id. """ -function get_run_id(path::String) - # pull list of existing run_???? folders via readdir +function get_run_id(path::String,id::String) + # if run_???? folder doesn't exist yet don't change the id + run_id = string("run_",run_id_to_string(id)) + !(run_id in readdir(path)) && return id + + # otherwise pull list of existing run_???? folders via readdir pattern = r"run_\d\d\d\d" # run_???? in regex runlist = filter(x->startswith(x,pattern),readdir(path)) runlist = filter(x->endswith( x,pattern),runlist) diff --git a/src/physics/boundary_layer.jl b/src/physics/boundary_layer.jl index 900775dff..cfe719b7f 100644 --- a/src/physics/boundary_layer.jl +++ b/src/physics/boundary_layer.jl @@ -16,7 +16,7 @@ end """Linear boundary layer drag Following Held and Suarez, 1996 BAMS $(TYPEDFIELDS)""" -@with_kw struct LinearDrag{NF<:AbstractFloat} <: BoundaryLayerDrag{NF} +Base.@kwdef struct LinearDrag{NF<:AbstractFloat} <: BoundaryLayerDrag{NF} # PARAMETERS σb::Float64 = 0.7 # sigma coordinate below which linear drag is applied time_scale::Float64 = 24 # [hours] time scale for linear drag coefficient at σ=1 (=1/kf in HS96) diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl index 978841296..ede3c9445 100644 --- a/src/physics/column_variables.jl +++ b/src/physics/column_variables.jl @@ -64,7 +64,7 @@ function write_column_tendencies!( D::DiagnosticVariables, layer.tendencies.humid_tend_grid[ij] = C.humid_tend[k] end - # D.surface.precip_large_scale[ij] = C.precip_large_scale + D.surface.precip_large_scale[ij] = C.precip_large_scale # D.surface.precip_convection[ij] = C.precip_convection return nothing @@ -75,7 +75,7 @@ end Set the accumulators (tendencies but also vertical sums and similar) back to zero for `column` to be reused at other grid points.""" -function reset_column!(column::ColumnVariables) +function reset_column!(column::ColumnVariables{NF}) where NF # set tendencies to 0 for += accumulation fill!(column.u_tend,0) @@ -101,8 +101,8 @@ function reset_column!(column::ColumnVariables) # fill!(column.net_flux_humid, 0) # fill!(column.net_flux_dry_static_energy, 0) - # # Large-scale condensation - # column.precip_large_scale = zero(NF) + # Large-scale condensation + column.precip_large_scale = zero(NF) return nothing end diff --git a/src/physics/convection.jl b/src/physics/convection.jl index 13af21ef7..5935cd12f 100644 --- a/src/physics/convection.jl +++ b/src/physics/convection.jl @@ -204,12 +204,6 @@ end # entrainment_profile /= sum(entrainment_profile) # Normalise # entrainment_profile *= max_entrainment # fraction of max entrainment - # # PARAMETRIZATIONS - # # Large-scale condensation (occurs when relative humidity exceeds a given threshold) - # RH_thresh_pbl_lsc::NF # Relative humidity threshold for LSC in PBL - # RH_thresh_range_lsc::NF # Vertical range of relative humidity threshold - # RH_thresh_max_lsc ::NF # Maximum relative humidity threshold - # humid_relax_time_lsc::NF # Relaxation time for humidity (hours) # # Convection # pres_thresh_cnv::NF # Minimum (normalised) surface pressure for the occurrence of convection @@ -223,22 +217,6 @@ end # "For computing saturation vapour pressure" # magnus_coefs::Coefficients = MagnusCoefs{NF}() - # # Large-Scale Condensation (from table B10) - # "Index of atmospheric level at which large-scale condensation begins" - # k_lsc::Int = 2 - - # "Relative humidity threshold for boundary layer" - # RH_thresh_pbl_lsc::Float64 = 0.95 - - # "Vertical range of relative humidity threshold" - # RH_thresh_range_lsc::Float64 = 0.1 - - # "Maximum relative humidity threshold" - # RH_thresh_max_lsc::Float64 = 0.9 - - # "Relaxation time for humidity (hours)" - # humid_relax_time_lsc::Float64 = 4.0 - # # Convection # "Minimum (normalised) surface pressure for the occurrence of convection" # pres_thresh_cnv::Float64 = 0.8 diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index a7d2a4f10..8d25aa15d 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -1,47 +1,123 @@ -function large_scale_condensation!( column::ColumnVariables{NF}, - model::PrimitiveEquation, - ) where {NF<:AbstractFloat} +""" +Large scale condensation as in Fortran SPEEDY with default values from therein. +$(TYPEDFIELDS)""" +@kwdef struct SpeedyCondensation{NF<:AbstractFloat} <: AbstractCondensation{NF} - (; gravity, RH_thresh_max_lsc, RH_thresh_range_lsc, RH_thresh_pbl_lsc, humid_relax_time_lsc ) = model.constants - (; σ_levels_full, σ_levels_thick, n_stratosphere_levels ) = model.geometry - (; cₚ, alhc ) = model.parameters + "number of vertical levels" + nlev::Int - (; humid, pres ) = column # prognostic variables: specific humidity, surface pressure - (; temp_tend, humid_tend ) = column # tendencies to write into - (; sat_humid ) = column # intermediate variable - (; nlev ) = column + # "Index of atmospheric level at which large-scale condensation begins" + # k_lsc::Int = 2 - # 1. Tendencies of humidity and temperature due to large-scale condensation - for k in eachlayer(column)[n_stratosphere_levels:end] # top to bottom, skip stratospheric levels + "Relative humidity threshold for boundary layer" + threshold_boundary_layer::Float64 = 0.95 + + "Vertical range of relative humidity threshold" + threshold_range::Float64 = 0.1 + + "Maximum relative humidity threshold [1]" + threshold_max::Float64 = 0.9 + + "Relaxation time for humidity [hrs]" + time_scale::Float64 = 4.0 + + # precomputed arrays + n_stratosphere_levels::Base.RefValue{Int} = Ref(0) + humid_tend_max::Vector{NF} = zeros(NF,nlev) + relative_threshold::Vector{NF} = zeros(NF,nlev) +end + +SpeedyCondensation(SG::SpectralGrid;kwargs...) = SpeedyCondensation{SG.NF}(nlev=SG.nlev;kwargs...) - # Relative humidity threshold for condensation (Formula 24) +""" +$(TYPEDSIGNATURES) +Initialize the SpeedyCondensation scheme.""" +function initialize!(scheme::SpeedyCondensation,model::PrimitiveEquation) + + (;nlev,threshold_boundary_layer,threshold_range,threshold_max) = scheme + (;humid_tend_max, relative_threshold, time_scale) = scheme + (;σ_levels_full) = model.geometry + (;σ_tropopause,σ_boundary_layer) = model.atmosphere + + scheme.n_stratosphere_levels[] = findfirst(σ->σ>=σ_tropopause,σ_levels_full) + + for k in 1:nlev + # the relative humidity threshold above which condensation occurs per layer σₖ² = σ_levels_full[k]^2 - RH_threshold = RH_thresh_max_lsc + RH_thresh_range_lsc * (σₖ² - 1) - if k == nlev - RH_threshold = max(RH_threshold, RH_thresh_pbl_lsc) + relative_threshold[k] = threshold_max + threshold_range * (σₖ² - 1) + if k >= σ_boundary_layer + relative_threshold[k] = max(relative_threshold[k], threshold_boundary_layer) 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σₖ² / 3600humid_relax_time_lsc - humid_threshold = RH_threshold * sat_humid[k] # Specific humidity threshold for condensation + humid_tend_max[k] = 10σₖ² / 3600time_scale # [hrs] → [s] + end +end + +""" +$(TYPEDSIGNATURES) +No condensation in a PrimitiveDry model.""" +function large_scale_condensation!( + column::ColumnVariables, + model::PrimitiveDry, +) + return nothing +end + +"""Function barrier only.""" +function large_scale_condensation!( + column::ColumnVariables, + model::PrimitiveWet, +) + large_scale_condensation!(column,model.large_scale_condensation, + model.geometry,model.constants,model.atmosphere) +end + + +function large_scale_condensation!( + column::ColumnVariables{NF}, + scheme::SpeedyCondensation, + geometry::Geometry, + constants::DynamicsConstants, + atmosphere::AbstractAtmosphere, + ) where NF + + (;relative_threshold,humid_tend_max) = scheme + time_scale = convert(NF,3600*scheme.time_scale) + n_stratosphere_levels = scheme.n_stratosphere_levels[] + + (;humid, pres) = column # prognostic variables: specific humidity, surface pressure + (;temp_tend, humid_tend) = column # tendencies to write into + (;sat_humid) = column # intermediate variable, calculate in thermodynamics! + (;nlev) = column + pₛ = pres[end] # surface pressure + + (;gravity) = constants + (;σ_levels_thick) = geometry + latent_heat = convert(NF, atmosphere.latent_heat_condensation/atmosphere.cₚ) + + # 1. Tendencies of humidity and temperature due to large-scale condensation + for k in n_stratosphere_levels+1:nlev # top to bottom, skip stratospheric levels + + # Specific humidity threshold for condensation + humid_threshold = relative_threshold[k] * sat_humid[k] if humid[k] > humid_threshold # accumulate in tendencies (nothing is added if humidity not above threshold) - humid_tend[k] += -(humid[k] - humid_threshold) / humid_relax_time_lsc # Formula 22 - temp_tend[k] += -alhc / cₚ * min(humid_tend[k], humid_tend_max*pres) # Formula 23 + humid_tend_k = -(humid[k] - humid_threshold) / time_scale # Formula 22 + temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23 # If there is large-scale condensation at a level higher (i.e. smaller k) than - # the cloud-top previously diagnosed due to convection, then increase the - # cloud-top. - column.cloud_top = min(column.cloud_top, k) # Page 7 (last sentence) - end - end + # the cloud-top previously diagnosed due to convection, then increase the cloud-top + column.cloud_top = min(column.cloud_top, k) # Page 7 (last sentence) + + # 2. Precipitation due to large-scale condensation + Δpₖ_g = pₛ*σ_levels_thick[k]/gravity # Formula 4 + column.precip_large_scale -= Δpₖ_g * humid_tend_k # Formula 25 - # 2. Precipitation due to large-scale condensation - pₛ = pres[end] # surface pressure - for k in eachlayer(column)[n_stratosphere_levels:end] # top to bottom, skip stratosphere - Δpₖ = pₛ*σ_levels_thick[k] # Formula 4 # Correct index? - column.precip_large_scale += -1 / gravity * Δpₖ * humid_tend[k] # Formula 25 + # only write into humid_tend now to allow humid_tend != 0 before this scheme is called + humid_tend[k] += humid_tend_k + end end -end +end \ No newline at end of file diff --git a/src/physics/shortwave_radiation.jl b/src/physics/shortwave_radiation.jl index 8ccc3e5e8..ed1c19f6e 100644 --- a/src/physics/shortwave_radiation.jl +++ b/src/physics/shortwave_radiation.jl @@ -1,7 +1,7 @@ """ Parameters for radiation parameterizations. """ -@with_kw struct RadiationCoefs{NF<:Real} <: Coefficients +Base.@kwdef struct RadiationCoefs{NF<:Real} <: Coefficients epslw::NF = 0.05 # Fraction of blackbody spectrum absorbed/emitted by PBL only emisfc::NF = 0.98 # Longwave surface emissivity diff --git a/src/physics/temperature_relaxation.jl b/src/physics/temperature_relaxation.jl index d2e8cda20..4e6a5188b 100644 --- a/src/physics/temperature_relaxation.jl +++ b/src/physics/temperature_relaxation.jl @@ -16,7 +16,7 @@ end """ Struct that defines the temperature relaxation from Held and Suarez, 1996 BAMS $(TYPEDFIELDS)""" -@with_kw struct HeldSuarez{NF<:AbstractFloat} <: TemperatureRelaxation{NF} +Base.@kwdef struct HeldSuarez{NF<:AbstractFloat} <: TemperatureRelaxation{NF} # DIMENSIONS "number of latitude rings" nlat::Int @@ -123,7 +123,7 @@ end """$(TYPEDSIGNATURES) HeldSuarez-like temperature relaxation, but towards the Jablonowski temperature profile with increasing temperatures in the stratosphere.""" -@with_kw struct JablonowskiRelaxation{NF<:AbstractFloat} <: TemperatureRelaxation{NF} +Base.@kwdef struct JablonowskiRelaxation{NF<:AbstractFloat} <: TemperatureRelaxation{NF} # DIMENSIONS nlat::Int nlev::Int diff --git a/src/physics/tendencies.jl b/src/physics/tendencies.jl index 4d55fd81d..a8b620c20 100644 --- a/src/physics/tendencies.jl +++ b/src/physics/tendencies.jl @@ -42,7 +42,7 @@ function parameterization_tendencies!( # Calculate parametrizations (order of execution is important!) # convection!(column,model) - # large_scale_condensation!(column,model) + large_scale_condensation!(column,model) # clouds!(column, model) # shortwave_radiation!(column,model) # longwave_radiation!(column,model) diff --git a/src/physics/thermodynamics.jl b/src/physics/thermodynamics.jl index 4ed63e2fc..ae8fae42c 100644 --- a/src/physics/thermodynamics.jl +++ b/src/physics/thermodynamics.jl @@ -6,7 +6,7 @@ Parameters for computing saturation vapour pressure using the August-Roche-Magnu where T is in Kelvin and i = 1,2 for saturation with respect to water and ice, respectively. $(TYPEDFIELDS)""" -@with_kw struct MagnusCoefs{NF<:AbstractFloat} +Base.@kwdef struct MagnusCoefs{NF<:AbstractFloat} "Saturation vapour pressure at 0°C" e₀::NF = 6.108 @@ -18,7 +18,7 @@ $(TYPEDFIELDS)""" C₂::NF = 21.875 end -@with_kw struct Thermodynamics{NF} <: AbstractThermodynamics{NF} +Base.@kwdef struct Thermodynamics{NF} <: AbstractThermodynamics{NF} magnus_coefs::MagnusCoefs{NF} = MagnusCoefs{NF}() mol_ratio::NF latent_heat_condensation::NF diff --git a/src/physics/vertical_diffusion.jl b/src/physics/vertical_diffusion.jl index 8d3470cd1..6dffeb212 100644 --- a/src/physics/vertical_diffusion.jl +++ b/src/physics/vertical_diffusion.jl @@ -15,7 +15,7 @@ end Diffusion of dry static energy: A relaxation towards a reference gradient of static energy wrt to geopotential, see Fortran SPEEDY documentation. $(TYPEDFIELDS)""" -@with_kw struct StaticEnergyDiffusion{NF<:AbstractFloat} <: VerticalDiffusion{NF} +Base.@kwdef struct StaticEnergyDiffusion{NF<:AbstractFloat} <: VerticalDiffusion{NF} "time scale [hrs] for strength" time_scale::Float64 = 6 @@ -60,7 +60,7 @@ function static_energy_diffusion!( column::ColumnVariables{NF}, end end -# @with_kw struct VerticalLaplacian{NF<:Real} <: VerticalDiffusion{NF} +# Base.@kwdef struct VerticalLaplacian{NF<:Real} <: VerticalDiffusion{NF} # time_scale::NF = 10.0 # [hours] time scale to control the strength of vertical diffusion # height_scale::NF = 100.0 # [m] scales for Δσ so that time_scale is sensible From e6355fd1ce9f892e30980079fc09265bcc4efad8 Mon Sep 17 00:00:00 2001 From: Milan Date: Wed, 31 May 2023 12:49:07 -0400 Subject: [PATCH 2/8] humidity initial conditions --- src/dynamics/initial_conditions.jl | 21 ++++++++++----------- src/dynamics/models.jl | 2 -- src/output/output.jl | 5 +++-- src/physics/large_scale_condensation.jl | 2 +- 4 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src/dynamics/initial_conditions.jl b/src/dynamics/initial_conditions.jl index c2ed983bd..d89ffd9c2 100644 --- a/src/dynamics/initial_conditions.jl +++ b/src/dynamics/initial_conditions.jl @@ -464,7 +464,7 @@ function initialize_humidity!( progn::PrognosticVariables, # reference saturation water vapour pressure [Pa] # relative humidity reference [1] - (;water_pres_ref, relhumid_ref, R_dry, R_vapour) = model.atmosphere + (;water_pres_ref, relhumid_ref, R_dry, R_vapour, pres_ref) = model.atmosphere gas_ratio = R_dry/R_vapour humid_ref = relhumid_ref*gas_ratio*water_pres_ref # reference specific humidity [Pa] @@ -477,20 +477,19 @@ function initialize_humidity!( progn::PrognosticVariables, # Specific humidity at the surface (grid space) humid_surf_grid = zero(pres_surf_grid) - @. humid_surf_grid = humid_ref*exp(scale_height_ratio*pres_surf_grid) - + # @. humid_surf_grid = humid_ref*(exp(pres_surf_grid)/(pres_ref*100))^scale_height_ratio + q_ref = 20e-3 # kg/kg at the surface + @. humid_surf_grid .= q_ref + scale_coslat²!(humid_surf_grid,model.geometry) + humid_surf = spectral(humid_surf_grid,model.spectral_transform) spectral_truncation!(humid_surf) - q0 = 3e-3 - q1 = -1e-3 - # Specific humidity at tropospheric levels (stratospheric humidity remains zero) + a = model.spectral_transform.norm_sphere for k in n_stratosphere_levels+1:nlev - # for lm in eachharmonic(humid_surf,progn.layers[1].timesteps[1].humid) - # progn.layers[k].timesteps[1].humid[lm] = humid_surf[lm]*σ_levels_full[k]^scale_height_ratio - # end - progn.layers[k].timesteps[1].humid[1] = model.spectral_transform.norm_sphere*q0*σ_levels_full[k]^scale_height_ratio - progn.layers[k].timesteps[1].humid[3] = model.spectral_transform.norm_sphere*q1*σ_levels_full[k]^scale_height_ratio + for lm in eachharmonic(humid_surf,progn.layers[1].timesteps[1].humid) + progn.layers[k].timesteps[1].humid[lm] = humid_surf[lm]*σ_levels_full[k]^scale_height_ratio + end end end \ No newline at end of file diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index 3c16dbae4..796609d24 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -162,7 +162,6 @@ end has(::Type{<:PrimitiveDry}, var_name::Symbol) = var_name in (:vor, :div, :temp, :pres) default_concrete_model(::Type{PrimitiveDry}) = PrimitiveDryModel -default_concrete_model(::Type{PrimitiveEquation}) = PrimitiveDryModel """ $(TYPEDSIGNATURES) @@ -231,7 +230,6 @@ end has(::Type{<:PrimitiveWet}, var_name::Symbol) = var_name in (:vor, :div, :temp, :pres, :humid) default_concrete_model(::Type{PrimitiveWet}) = PrimitiveWetModel -default_concrete_model(::Type{PrimitiveEquation}) = PrimitiveDryModel """ $(TYPEDSIGNATURES) diff --git a/src/output/output.jl b/src/output/output.jl index e06ed28dc..67717f062 100644 --- a/src/output/output.jl +++ b/src/output/output.jl @@ -136,7 +136,8 @@ function initialize!( (;startdate) = output time_string = "seconds since $(Dates.format(startdate, "yyyy-mm-dd HH:MM:0.0"))" dim_time = NcDim("time",0,unlimited=true) - var_time = NcVar("time",dim_time,t=Int64,atts=Dict("units"=>time_string,"long_name"=>"time")) + var_time = NcVar("time",dim_time,t=Int64,atts=Dict("units"=>time_string,"long_name"=>"time", + "standard_name"=>"time","calendar"=>"proleptic_gregorian")) # DEFINE NETCDF DIMENSIONS SPACE (;input_Grid, output_Grid, nlat_half, nlev) = output @@ -192,7 +193,7 @@ function initialize!( atts=Dict("long_name"=>"temperature","units"=>"degC","missing_value"=>missing_value, "_FillValue"=>missing_value)), humid = NcVar("humid",[dim_lon,dim_lat,dim_lev,dim_time],t=output_NF,compress=compression_level, - atts=Dict("long_name"=>"specific humidity","units"=>"1","missing_value"=>missing_value, + atts=Dict("long_name"=>"specific humidity","units"=>"kg/kg","missing_value"=>missing_value, "_FillValue"=>missing_value)), orog = NcVar("orog",[dim_lon,dim_lat],t=output_NF,compress=compression_level, atts=Dict("long_name"=>"orography","units"=>"m","missing_value"=>missing_value, diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index 8d25aa15d..48246b5f1 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -1,7 +1,7 @@ """ Large scale condensation as in Fortran SPEEDY with default values from therein. $(TYPEDFIELDS)""" -@kwdef struct SpeedyCondensation{NF<:AbstractFloat} <: AbstractCondensation{NF} +Base.@kwdef struct SpeedyCondensation{NF<:AbstractFloat} <: AbstractCondensation{NF} "number of vertical levels" nlev::Int From f703f096889f6c012ec97572cb9dae18f8572187 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 1 Jun 2023 15:22:26 -0400 Subject: [PATCH 3/8] Magnus formula [hPa] -> [Pa] --- src/physics/thermodynamics.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/thermodynamics.jl b/src/physics/thermodynamics.jl index ae8fae42c..aa7a7b8a5 100644 --- a/src/physics/thermodynamics.jl +++ b/src/physics/thermodynamics.jl @@ -7,8 +7,8 @@ where T is in Kelvin and i = 1,2 for saturation with respect to water and ice, respectively. $(TYPEDFIELDS)""" Base.@kwdef struct MagnusCoefs{NF<:AbstractFloat} - "Saturation vapour pressure at 0°C" - e₀::NF = 6.108 + "Saturation vapour pressure at 0°C [Pa]" + e₀::NF = 610.8 "0°C in Kelvin" T₀::NF = 273.16 From 2cd86499462d24c552fbbe7406b93745908888ca Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 1 Jun 2023 18:08:17 -0400 Subject: [PATCH 4/8] ITS RAINING --- src/SpeedyWeather.jl | 2 +- src/dynamics/atmospheres.jl | 3 ++ src/dynamics/constants.jl | 7 ++- src/dynamics/horizontal_diffusion.jl | 7 +-- src/dynamics/initial_conditions.jl | 6 +-- src/dynamics/models.jl | 2 +- src/output/output.jl | 68 ++++++++++++++++++------- src/physics/column_variables.jl | 7 +-- src/physics/large_scale_condensation.jl | 9 ++-- 9 files changed, 76 insertions(+), 35 deletions(-) diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index c4ba991b7..f75de466c 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -79,7 +79,7 @@ export StartFromFile, StartFromRest, ZonalJet, ZonalWind, - StartWithVorticity + StartWithRandomVorticity # EXPORT TEMPERATURE RELAXATION SCHEMES export NoTemperatureRelaxation, diff --git a/src/dynamics/atmospheres.jl b/src/dynamics/atmospheres.jl index 00ee3d692..7b62aa3d4 100644 --- a/src/dynamics/atmospheres.jl +++ b/src/dynamics/atmospheres.jl @@ -24,6 +24,9 @@ Base.@kwdef struct EarthAtmosphere <: AbstractAtmosphere "specific gas constant for water vapour [J/kg/K]" R_vapour::Float64 = 1000*R_gas/mol_mass_vapour + "water density [kg/m³]" + water_density::Float64 = 1000 + "latent heat of condensation [J/g] for consistency with specific humidity [g/Kg], also called alhc" latent_heat_condensation::Float64 = 2501 diff --git a/src/dynamics/constants.jl b/src/dynamics/constants.jl index 2a17edab7..61b004be1 100644 --- a/src/dynamics/constants.jl +++ b/src/dynamics/constants.jl @@ -31,6 +31,9 @@ Base.@kwdef struct DynamicsConstants{NF<:AbstractFloat} <: AbstractDynamicsConst "= R_dry/cₚ, gas const for air over heat capacity" κ::NF + "water density [kg/m³]" + water_density::NF + "coriolis frequency [1/s] (scaled by radius as is vorticity) = 2Ω*sin(lat)*radius" f_coriolis::Vector{NF} @@ -63,7 +66,7 @@ function DynamicsConstants( spectral_grid::SpectralGrid, geometry::Geometry) # PHYSICAL CONSTANTS - (;R_dry, R_vapour, lapse_rate, cₚ) = atmosphere + (;R_dry, R_vapour, lapse_rate, cₚ, water_density) = atmosphere (;ΔT_stratosphere, σ_tropopause, temp_ref) = atmosphere (;NF, radius) = spectral_grid (;rotation, gravity) = planet @@ -102,7 +105,7 @@ function DynamicsConstants( spectral_grid::SpectralGrid, # This implies conversion to NF return DynamicsConstants{NF}(; radius,rotation,gravity,layer_thickness, - R_dry,R_vapour,μ_virt_temp,cₚ,κ, + R_dry,R_vapour,μ_virt_temp,cₚ,κ,water_density, f_coriolis, σ_lnp_A,σ_lnp_B, Δp_geopot_half, Δp_geopot_full, diff --git a/src/dynamics/horizontal_diffusion.jl b/src/dynamics/horizontal_diffusion.jl index 192e72f82..3d3a51397 100644 --- a/src/dynamics/horizontal_diffusion.jl +++ b/src/dynamics/horizontal_diffusion.jl @@ -260,13 +260,14 @@ function horizontal_diffusion!( progn::PrognosticLayerTimesteps, ∇²ⁿ_implicit = HD.∇²ⁿ_implicit[k] # Primitive equation models diffuse vor and divergence more selective/adaptive - (;vor,div,temp) = progn.timesteps[lf] - (;vor_tend,div_tend,temp_tend) = diagn.tendencies + (;vor,div,temp,humid) = progn.timesteps[lf] + (;vor_tend,div_tend,temp_tend,humid_tend) = diagn.tendencies horizontal_diffusion!(vor_tend,vor,∇²ⁿ,∇²ⁿ_implicit) horizontal_diffusion!(div_tend,div,∇²ⁿ,∇²ⁿ_implicit) - # but use the weaker normal diffusion for temperature + # but use the weaker normal diffusion for temperature, humidity ∇²ⁿ = HD.∇²ⁿ_2D ∇²ⁿ_implicit = HD.∇²ⁿ_2D_implicit horizontal_diffusion!(temp_tend,temp,∇²ⁿ,∇²ⁿ_implicit) + model isa PrimitiveWet && horizontal_diffusion!(humid_tend,humid,∇²ⁿ,∇²ⁿ_implicit) end \ No newline at end of file diff --git a/src/dynamics/initial_conditions.jl b/src/dynamics/initial_conditions.jl index b835e3209..18263d460 100644 --- a/src/dynamics/initial_conditions.jl +++ b/src/dynamics/initial_conditions.jl @@ -1,5 +1,5 @@ # default initial conditions by model -initial_conditions_default(::Type{<:Barotropic}) = StartWithVorticity() +initial_conditions_default(::Type{<:Barotropic}) = StartWithRandomVorticity() initial_conditions_default(::Type{<:ShallowWater}) = ZonalJet() initial_conditions_default(::Type{<:PrimitiveEquation}) = ZonalWind() @@ -46,7 +46,7 @@ end """Start with random vorticity as initial conditions $(TYPEDFIELDS)""" -Base.@kwdef struct StartWithVorticity <: InitialConditions +Base.@kwdef struct StartWithRandomVorticity <: InitialConditions "Power law the vorticity should be spectrally distributed by" power_law::Float64 = -3 @@ -59,7 +59,7 @@ end $(TYPEDSIGNATURES) Start with random vorticity as initial conditions""" function initial_conditions!( progn::PrognosticVariables{NF}, - initial_conditions::StartWithVorticity, + initial_conditions::StartWithRandomVorticity, model::ModelSetup) where NF lmax = progn.trunc+1 diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index e6af03a39..cc92bc8cc 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -28,7 +28,7 @@ Base.@kwdef struct BarotropicModel{NF<:AbstractFloat, D<:AbstractDevice} <: Baro planet::AbstractPlanet = Earth() atmosphere::AbstractAtmosphere = EarthAtmosphere() forcing::AbstractForcing{NF} = NoForcing(spectral_grid) - initial_conditions::InitialConditions = StartWithVorticity() + initial_conditions::InitialConditions = StartWithRandomVorticity() # NUMERICS time_stepping::TimeStepper{NF} = Leapfrog(spectral_grid) diff --git a/src/output/output.jl b/src/output/output.jl index f19454ef1..f8e486f1e 100644 --- a/src/output/output.jl +++ b/src/output/output.jl @@ -10,6 +10,8 @@ Base.@kwdef struct Keepbits temp::Int = 10 pres::Int = 12 humid::Int = 7 + precip_cond::Int = 7 + precip_conv::Int = 7 end function Base.show(io::IO,K::Keepbits) @@ -57,6 +59,9 @@ Base.@kwdef mutable struct OutputWriter{NF<:Union{Float32,Float64},Model<:ModelS "[OPTION] output frequency, time step [hrs]" output_dt::Float64 = 6 + "actual output time step [sec]" + output_dt_sec::Int = 0 + "[OPTION] which variables to output, u, v, vor, div, pres, temp, humid" output_vars::Vector{Symbol} = default_output_vars(Model) @@ -111,6 +116,8 @@ Base.@kwdef mutable struct OutputWriter{NF<:Union{Float32,Float64},Model<:ModelS const temp::Matrix{NF} = fill(missing_value,nlon,nlat) const pres::Matrix{NF} = fill(missing_value,nlon,nlat) const humid::Matrix{NF} = fill(missing_value,nlon,nlat) + const precip_cond::Matrix{NF} = fill(missing_value,nlon,nlat) + const precip_conv::Matrix{NF} = fill(missing_value,nlon,nlat) end # generator function pulling grid resolution and time stepping from ::SpectralGrid and ::TimeStepper @@ -127,7 +134,7 @@ end default_output_vars(::Type{<:Barotropic}) = [:vor,:u] default_output_vars(::Type{<:ShallowWater}) = [:vor,:u] default_output_vars(::Type{<:PrimitiveDry}) = [:vor,:u,:temp,:pres] -default_output_vars(::Type{<:PrimitiveWet}) = [:vor,:u,:temp,:humid,:pres] +default_output_vars(::Type{<:PrimitiveWet}) = [:vor,:u,:temp,:humid,:pres,:precip_cond,:precip_conv] # print all fields with type <: Number function Base.show(io::IO,O::AbstractOutputWriter) @@ -219,7 +226,13 @@ function initialize!( "_FillValue"=>missing_value)), orog = NcVar("orog",[dim_lon,dim_lat],t=output_NF,compress=compression_level, atts=Dict("long_name"=>"orography","units"=>"m","missing_value"=>missing_value, - "_FillValue"=>missing_value)) + "_FillValue"=>missing_value)), + precip_cond = NcVar("precip_cond",[dim_lon,dim_lat,dim_time],t=output_NF,compress=compression_level, + atts=Dict("long_name"=>"Large-scale precipitation","units"=>"mm/dt","missing_value"=>missing_value, + "_FillValue"=>missing_value)), + precip_conv = NcVar("precip_conv",[dim_lon,dim_lat,dim_time],t=output_NF,compress=compression_level, + atts=Dict("long_name"=>"Convective precipitation","units"=>"mm/dt","missing_value"=>missing_value, + "_FillValue"=>missing_value)), ) # GET RUN ID, CREATE FOLDER @@ -232,7 +245,6 @@ function initialize!( feedback.progress_meter.desc = "Weather is speedy: run $(output.id) " feedback.output = true # if output=true set feedback.output=true too! - # CREATE NETCDF FILE, vector of NcVars for output (; run_path, filename, output_vars) = output vars_out = [all_ncvars[key] for key in keys(all_ncvars) if key in vcat(:time,output_vars)] @@ -244,7 +256,8 @@ function initialize!( # OUTPUT FREQUENCY output.output_every_n_steps = max(1,floor(Int,output.output_dt/time_stepping.Δt_hrs)) - + output.output_dt_sec = output.output_every_n_steps*time_stepping.Δt_sec + # RESET COUNTERS output.timestep_counter = 0 # time step counter output.output_counter = 0 # output step counter @@ -357,7 +370,7 @@ function write_netcdf_variables!( output::OutputWriter, (;output_vars) = output # Vector{Symbol} of variables to output i = output.output_counter - (;u, v, vor, div, pres, temp, humid) = output + (;u, v, vor, div, pres, temp, humid, precip_cond, precip_conv) = output (;output_Grid, interpolator) = output (;quadrant_rotation, matrix_quadrant) = output @@ -367,7 +380,9 @@ function write_netcdf_variables!( output::OutputWriter, if output.as_matrix # resort gridded variables interpolation-free into a matrix - # create (matrix,grid) tuples for simultaneous grid -> matrix conversion + # create (matrix,grid) tuples for simultaneous grid -> matrix conversion + # TODO this currently does the Matrix! conversion to all variables, not just output_vars + # as arrays are always initialised MGs = ((M,G) for (M,G) in zip((u,v,vor,div,temp,humid), (u_grid,v_grid,vor_grid,div_grid,temp_grid,humid_grid)) if length(M) > 0) @@ -408,24 +423,41 @@ function write_netcdf_variables!( output::OutputWriter, end # surface pressure, i.e. interface displacement η - if :pres in output_vars - (; pres_grid ) = diagn.surface + (; pres_grid, precip_large_scale, precip_convection ) = diagn.surface - if output.as_matrix - RingGrids.Matrix!(pres,diagn.surface.pres_grid; quadrant_rotation, matrix_quadrant) - else - RingGrids.interpolate!(output_Grid(pres),pres_grid,interpolator) - end + if output.as_matrix + if :pres in output_vars || :precip_cond in output_vars || :precip_conv in output_vars + MGs = ((M,G) for (M,G) in zip((pres,precip_cond,precip_conv), + (pres_grid, precip_large_scale, precip_convection))) - if Model <: PrimitiveEquation - @. pres = exp(pres)/100 # convert from log(pₛ) to surface pressure pₛ [hPa] + RingGrids.Matrix!(MGs...; quadrant_rotation, matrix_quadrant) end + else + :pres in output_vars && RingGrids.interpolate!(output_Grid(pres),pres_grid,interpolator) + :precip_cond in output_vars && RingGrids.interpolate!(output_Grid(precip_cond),precip_large_scale,interpolator) + :precip_conv in output_vars && RingGrids.interpolate!(output_Grid(precip_conv),precip_convection,interpolator) + end - round!(pres,output.keepbits.pres) - - NetCDF.putvar(output.netcdf_file,"pres",pres,start=[1,1,i],count=[-1,-1,1]) + # after output set precip accumulators back to zero + precip_large_scale .= 0 + precip_convection .= 0 + + # convert from [m/s] to [mm/Δt], where Δt is the output time step + precip_cond *= 1000*output.output_dt_sec + precip_conv *= 1000*output.output_dt_sec + + if Model <: PrimitiveEquation + @. pres = exp(pres)/100 # convert from log(pₛ) to surface pressure pₛ [hPa] end + :pres in output_vars && round!(pres,output.keepbits.pres) + :precip_cond in output_vars && round!(precip_cond,output.keepbits.precip_cond) + :precip_conv in output_vars && round!(precip_conv,output.keepbits.precip_conv) + + :pres in output_vars && NetCDF.putvar(output.netcdf_file,"pres",pres,start=[1,1,i],count=[-1,-1,1]) + :precip_cond in output_vars && NetCDF.putvar(output.netcdf_file,"precip_cond",precip_cond,start=[1,1,i],count=[-1,-1,1]) + :precip_conv in output_vars && NetCDF.putvar(output.netcdf_file,"precip_conv",precip_conv,start=[1,1,i],count=[-1,-1,1]) + return nothing end diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl index ede3c9445..0f66464c2 100644 --- a/src/physics/column_variables.jl +++ b/src/physics/column_variables.jl @@ -64,8 +64,9 @@ function write_column_tendencies!( D::DiagnosticVariables, layer.tendencies.humid_tend_grid[ij] = C.humid_tend[k] end - D.surface.precip_large_scale[ij] = C.precip_large_scale - # D.surface.precip_convection[ij] = C.precip_convection + # accumulate (set back to zero when netcdf output) + D.surface.precip_large_scale[ij] += C.precip_large_scale + D.surface.precip_convection[ij] += C.precip_convection return nothing end @@ -97,7 +98,7 @@ function reset_column!(column::ColumnVariables{NF}) where NF # column.cloud_top = column.nlev+1 # column.conditional_instability = false # column.activate_convection = false - # column.precip_convection = zero(NF) + column.precip_convection = zero(NF) # fill!(column.net_flux_humid, 0) # fill!(column.net_flux_dry_static_energy, 0) diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index 48246b5f1..d092ccf27 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -93,7 +93,7 @@ function large_scale_condensation!( (;nlev) = column pₛ = pres[end] # surface pressure - (;gravity) = constants + (;gravity, water_density) = constants (;σ_levels_thick) = geometry latent_heat = convert(NF, atmosphere.latent_heat_condensation/atmosphere.cₚ) @@ -112,9 +112,10 @@ function large_scale_condensation!( # the cloud-top previously diagnosed due to convection, then increase the cloud-top column.cloud_top = min(column.cloud_top, k) # Page 7 (last sentence) - # 2. Precipitation due to large-scale condensation - Δpₖ_g = pₛ*σ_levels_thick[k]/gravity # Formula 4 - column.precip_large_scale -= Δpₖ_g * humid_tend_k # Formula 25 + # 2. Precipitation due to large-scale condensation [kg/m²/s] /ρ for [m/s] + # += for vertical integral + Δpₖ_g = pₛ*σ_levels_thick[k]/gravity/water_density # Formula 4 + column.precip_large_scale += -Δpₖ_g * humid_tend_k # Formula 25 # only write into humid_tend now to allow humid_tend != 0 before this scheme is called humid_tend[k] += humid_tend_k From b13d25ad5fae9968756b392d71cd961cc21dc517 Mon Sep 17 00:00:00 2001 From: Milan Date: Sat, 3 Jun 2023 09:41:06 -0400 Subject: [PATCH 5/8] precip units in [m] to [mm] for output --- src/SpeedyWeather.jl | 2 +- src/output/output.jl | 6 ++--- src/physics/column_variables.jl | 34 ++++++++++++------------- src/physics/define_column.jl | 30 +++++++++++----------- src/physics/large_scale_condensation.jl | 14 ++++++---- 5 files changed, 44 insertions(+), 42 deletions(-) diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index 3ca257fcb..194890462 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -168,6 +168,7 @@ include("physics/thermodynamics.jl") include("physics/boundary_layer.jl") include("physics/temperature_relaxation.jl") include("physics/vertical_diffusion.jl") +include("physics/large_scale_condensation.jl") include("physics/pretty_printing.jl") # MODELS @@ -175,7 +176,6 @@ include("dynamics/models.jl") # # PHYSICS # include("physics/convection.jl") -include("physics/large_scale_condensation.jl") # include("physics/longwave_radiation.jl") # include("physics/shortwave_radiation.jl") diff --git a/src/output/output.jl b/src/output/output.jl index f8e486f1e..7a62bad00 100644 --- a/src/output/output.jl +++ b/src/output/output.jl @@ -442,9 +442,9 @@ function write_netcdf_variables!( output::OutputWriter, precip_large_scale .= 0 precip_convection .= 0 - # convert from [m/s] to [mm/Δt], where Δt is the output time step - precip_cond *= 1000*output.output_dt_sec - precip_conv *= 1000*output.output_dt_sec + # convert from [m] to [mm] within output time step (e.g. 6hours) + precip_cond *= 1000 + precip_conv *= 1000 if Model <: PrimitiveEquation @. pres = exp(pres)/100 # convert from log(pₛ) to surface pressure pₛ [hPa] diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl index 0f66464c2..33edc17f5 100644 --- a/src/physics/column_variables.jl +++ b/src/physics/column_variables.jl @@ -1,6 +1,5 @@ """ - get_column!(C,D,ij,G) - +$(TYPEDSIGNATURES) Update `C::ColumnVariables` by copying the prognostic variables from `D::DiagnosticVariables` at gridpoint index `ij`. Provide `G::Geometry` for coordinate information.""" function get_column!( C::ColumnVariables, @@ -31,6 +30,7 @@ function get_column!( C::ColumnVariables, C.humid[k] = layer.grid_variables.humid_grid[ij] # as well as geopotential (not actually prognostic though) + # TODO geopot on the grid is currently not computed in dynamics C.geopot[k] = layer.grid_variables.geopot_grid[ij] end end @@ -47,8 +47,7 @@ function get_column!( C::ColumnVariables, end """ - write_column_tendencies!(D,C,ij,G) - +$(TYPEDSIGNATURES) Write the parametrization tendencies from `C::ColumnVariables` into the horizontal fields of tendencies stored in `D::DiagnosticVariables` at gridpoint index `ij`.""" function write_column_tendencies!( D::DiagnosticVariables, @@ -72,27 +71,26 @@ function write_column_tendencies!( D::DiagnosticVariables, end """ - reset_column!(column::ColumnVariables) - +$(TYPEDSIGNATURES) Set the accumulators (tendencies but also vertical sums and similar) back to zero for `column` to be reused at other grid points.""" function reset_column!(column::ColumnVariables{NF}) where NF # set tendencies to 0 for += accumulation - fill!(column.u_tend,0) - fill!(column.v_tend,0) - fill!(column.temp_tend,0) - fill!(column.humid_tend,0) + column.u_tend .= 0 + column.v_tend .= 0 + column.temp_tend .= 0 + column.humid_tend .= 0 # set fluxes to 0 for += accumulation - fill!(column.flux_u_upward,0) - fill!(column.flux_u_downward,0) - fill!(column.flux_v_upward,0) - fill!(column.flux_v_downward,0) - fill!(column.flux_humid_upward,0) - fill!(column.flux_humid_downward,0) - fill!(column.flux_temp_upward,0) - fill!(column.flux_temp_downward,0) + column.flux_u_upward .= 0 + column.flux_u_downward .= 0 + column.flux_v_upward .= 0 + column.flux_v_downward .= 0 + column.flux_humid_upward .= 0 + column.flux_humid_downward .= 0 + column.flux_temp_upward .= 0 + column.flux_temp_downward .= 0 # # Convection # column.cloud_top = column.nlev+1 diff --git a/src/physics/define_column.jl b/src/physics/define_column.jl index 548f101c2..36ef159bf 100644 --- a/src/physics/define_column.jl +++ b/src/physics/define_column.jl @@ -17,14 +17,14 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV latd::NF = 0 # latitude, needed for shortwave radiation # PROGNOSTIC VARIABLES - const u::Vector{NF} = zeros(NF,nlev) # zonal velocity - const v::Vector{NF} = zeros(NF,nlev) # meridional velocity - const temp::Vector{NF} = zeros(NF,nlev) # temperature - const humid::Vector{NF} = zeros(NF,nlev) # specific humidity + const u::Vector{NF} = zeros(NF,nlev) # zonal velocity [m/s] + const v::Vector{NF} = zeros(NF,nlev) # meridional velocity [m/s] + const temp::Vector{NF} = zeros(NF,nlev) # temperature [K] + const humid::Vector{NF} = zeros(NF,nlev) # specific humidity [kg/kg] # (log) pressure per layer, surface is prognostic, last element here, but precompute other layers too - const ln_pres::Vector{NF} = zeros(NF,nlev+1) # logarithm of pressure [log(hPa)] - const pres::Vector{NF} = zeros(NF,nlev+1) # pressure [hPa] + const ln_pres::Vector{NF} = zeros(NF,nlev+1) # logarithm of pressure [log(Pa)] + const pres::Vector{NF} = zeros(NF,nlev+1) # pressure [Pa] # TENDENCIES to accumulate the parametrizations into const u_tend::Vector{NF} = zeros(NF,nlev) # zonal velocity [m] @@ -49,8 +49,8 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV const flux_humid_downward::Vector{NF} = zeros(NF,nlev+1) # THERMODYNAMICS - const sat_humid::Vector{NF} = zeros(NF,nlev) # Saturation specific humidity - const sat_vap_pres::Vector{NF} = zeros(NF,nlev) # Saturation vapour pressure + const sat_humid::Vector{NF} = zeros(NF,nlev) # Saturation specific humidity [kg/kg] + const sat_vap_pres::Vector{NF} = zeros(NF,nlev) # Saturation vapour pressure [Pa] const dry_static_energy::Vector{NF} = zeros(NF,nlev) # Dry static energy const moist_static_energy::Vector{NF} = zeros(NF,nlev) # Moist static energy @@ -62,18 +62,18 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV sat_moist_static_energy_half::Vector{NF} = zeros(NF,nlev) # Saturation moist static energy interpolated to half-levels # Convection - conditional_instability::Bool = false # Whether a conditional instability exists in this column (condition 1) - activate_convection::Bool = false # Whether convection should be activated in this column (condition 2) - cloud_top::Int = nlev+1 # Top-of-convection layer - excess_humidity::NF = 0 # Excess humidity due to convection - cloud_base_mass_flux::NF = 0 # Mass flux at the top of the PBL - precip_convection::NF = 0 # Precipitation due to convection + conditional_instability::Bool = false # in this column? (condition 1) + activate_convection::Bool = false # in this column? (condition 2) + cloud_top::Int = nlev+1 # Top-of-convection layer + excess_humidity::NF = 0 # due to convection + cloud_base_mass_flux::NF = 0 # Mass flux at the top of the PBL + precip_convection::NF = 0 # Precipitation due to convection [m] net_flux_humid::Vector{NF} = zeros(NF,nlev) # Net fluxes of moisture in this column net_flux_dry_static_energy::Vector{NF} = zeros(NF,nlev) # Net fluxes of dry static energy in this column entrainment_profile::Vector{NF} = zeros(NF,nlev) # Entrainment coefficients # Large-scale condensation - precip_large_scale::NF = 0 # Precipitation due to large-scale condensation + precip_large_scale::NF = 0 # precipitation [m] # Longwave radiation ## New vars in radlw_down! diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index d092ccf27..8533c4978 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -71,7 +71,7 @@ function large_scale_condensation!( model::PrimitiveWet, ) large_scale_condensation!(column,model.large_scale_condensation, - model.geometry,model.constants,model.atmosphere) + model.geometry,model.constants,model.atmosphere,model.time_stepping) end @@ -81,6 +81,7 @@ function large_scale_condensation!( geometry::Geometry, constants::DynamicsConstants, atmosphere::AbstractAtmosphere, + time_stepping::TimeStepper, ) where NF (;relative_threshold,humid_tend_max) = scheme @@ -94,6 +95,9 @@ function large_scale_condensation!( pₛ = pres[end] # surface pressure (;gravity, water_density) = constants + (;Δt_sec) = time_stepping + Δt_gρ = Δt_sec/gravity/water_density + (;σ_levels_thick) = geometry latent_heat = convert(NF, atmosphere.latent_heat_condensation/atmosphere.cₚ) @@ -106,16 +110,16 @@ function large_scale_condensation!( if humid[k] > humid_threshold # accumulate in tendencies (nothing is added if humidity not above threshold) humid_tend_k = -(humid[k] - humid_threshold) / time_scale # Formula 22 - temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23 + temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23 # If there is large-scale condensation at a level higher (i.e. smaller k) than # the cloud-top previously diagnosed due to convection, then increase the cloud-top - column.cloud_top = min(column.cloud_top, k) # Page 7 (last sentence) + column.cloud_top = min(column.cloud_top, k) # Page 7 (last sentence) # 2. Precipitation due to large-scale condensation [kg/m²/s] /ρ for [m/s] # += for vertical integral - Δpₖ_g = pₛ*σ_levels_thick[k]/gravity/water_density # Formula 4 - column.precip_large_scale += -Δpₖ_g * humid_tend_k # Formula 25 + Δpₖ_g = pₛ*σ_levels_thick[k]*Δt_gρ # Formula 4 *Δt for [m] of rain during Δt + column.precip_large_scale += -Δpₖ_g * humid_tend_k # Formula 25, unit [m] # only write into humid_tend now to allow humid_tend != 0 before this scheme is called humid_tend[k] += humid_tend_k From 9617a2c7a9035ecea17081b8d80cf3bd287c9dfc Mon Sep 17 00:00:00 2001 From: Milan Date: Sat, 3 Jun 2023 10:00:48 -0400 Subject: [PATCH 6/8] docstrings --- src/physics/large_scale_condensation.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index 8533c4978..4f4fee1cd 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -6,9 +6,6 @@ Base.@kwdef struct SpeedyCondensation{NF<:AbstractFloat} <: AbstractCondensation "number of vertical levels" nlev::Int - # "Index of atmospheric level at which large-scale condensation begins" - # k_lsc::Int = 2 - "Relative humidity threshold for boundary layer" threshold_boundary_layer::Float64 = 0.95 @@ -74,7 +71,12 @@ function large_scale_condensation!( model.geometry,model.constants,model.atmosphere,model.time_stepping) end - +""" +$(TYPEDSIGNATURES) +Large-scale condensation for a `column` by relaxation back to a reference +relative humidity if larger than that. Calculates the tendencies for +specific humidity and temperature and integrates the large-scale +precipitation vertically for output.""" function large_scale_condensation!( column::ColumnVariables{NF}, scheme::SpeedyCondensation, @@ -96,13 +98,13 @@ function large_scale_condensation!( (;gravity, water_density) = constants (;Δt_sec) = time_stepping - Δt_gρ = Δt_sec/gravity/water_density + pₛΔt_gρ = pₛ*Δt_sec/gravity/water_density # precompute constant (;σ_levels_thick) = geometry latent_heat = convert(NF, atmosphere.latent_heat_condensation/atmosphere.cₚ) # 1. Tendencies of humidity and temperature due to large-scale condensation - for k in n_stratosphere_levels+1:nlev # top to bottom, skip stratospheric levels + @inbounds for k in n_stratosphere_levels+1:nlev # top to bottom, skip stratospheric levels # Specific humidity threshold for condensation humid_threshold = relative_threshold[k] * sat_humid[k] @@ -114,12 +116,12 @@ function large_scale_condensation!( # If there is large-scale condensation at a level higher (i.e. smaller k) than # the cloud-top previously diagnosed due to convection, then increase the cloud-top - column.cloud_top = min(column.cloud_top, k) # Page 7 (last sentence) + column.cloud_top = min(column.cloud_top, k) # Page 7 (last sentence) # 2. Precipitation due to large-scale condensation [kg/m²/s] /ρ for [m/s] # += for vertical integral - Δpₖ_g = pₛ*σ_levels_thick[k]*Δt_gρ # Formula 4 *Δt for [m] of rain during Δt - column.precip_large_scale += -Δpₖ_g * humid_tend_k # Formula 25, unit [m] + ΔpₖΔt_gρ = σ_levels_thick[k]*pₛΔt_gρ # Formula 4 *Δt for [m] of rain during Δt + column.precip_large_scale += -ΔpₖΔt_gρ * humid_tend_k # Formula 25, unit [m] # only write into humid_tend now to allow humid_tend != 0 before this scheme is called humid_tend[k] += humid_tend_k From 41461aedee05768085f8708f6ea385da206ffa6a Mon Sep 17 00:00:00 2001 From: Milan Date: Sun, 4 Jun 2023 18:02:47 -0400 Subject: [PATCH 7/8] plot(::AbstractGrid) --- src/RingGrids/RingGrids.jl | 9 +++-- src/RingGrids/grids_general.jl | 1 + src/RingGrids/interpolation.jl | 50 +++++++++++++++++----------- src/RingGrids/show.jl | 26 +++++++++++++++ src/SpeedyWeather.jl | 3 +- src/dynamics/prognostic_variables.jl | 23 ++----------- 6 files changed, 69 insertions(+), 43 deletions(-) diff --git a/src/RingGrids/RingGrids.jl b/src/RingGrids/RingGrids.jl index f59e27d2c..f2c67bde4 100644 --- a/src/RingGrids/RingGrids.jl +++ b/src/RingGrids/RingGrids.jl @@ -1,7 +1,10 @@ module RingGrids -using DocStringExtensions +# DOCUMENTATION AND VISUALISATION +using DocStringExtensions +import UnicodePlots +# NUMERICS import Statistics: mean import FastGaussQuadrature @@ -55,14 +58,16 @@ export interpolate, update_locator, update_locator! +export plot + include("utility_functions.jl") include("grids_general.jl") -include("show.jl") include("full_grids.jl") include("octahedral.jl") include("healpix.jl") include("octahealpix.jl") include("quadrature_weights.jl") include("interpolation.jl") +include("show.jl") end \ No newline at end of file diff --git a/src/RingGrids/grids_general.jl b/src/RingGrids/grids_general.jl index 33dd02574..53dec4916 100644 --- a/src/RingGrids/grids_general.jl +++ b/src/RingGrids/grids_general.jl @@ -133,6 +133,7 @@ function get_latdlonds(Grid::Type{<:AbstractGrid},nlat_half::Integer) end get_latd(grid::Grid) where {Grid<:AbstractGrid} = get_latd(Grid,grid.nlat_half) +get_lond(grid::Grid) where {Grid<:AbstractGrid} = get_lond(Grid,grid.nlat_half) function get_latd(Grid::Type{<:AbstractGrid},nlat_half::Integer) colat = get_colat(Grid,nlat_half) diff --git a/src/RingGrids/interpolation.jl b/src/RingGrids/interpolation.jl index a618f284a..a2d10a105 100644 --- a/src/RingGrids/interpolation.jl +++ b/src/RingGrids/interpolation.jl @@ -407,15 +407,11 @@ function find_lon_indices( λ::NF, # longitude to find incides for (0˚... end """ - N,S = average_on_poles( A::AbstractGrid, - rings::Vector{<:UnitRange}) - -N,S are the interpolated values of A onto the north/south pole, by averaging all values on the -northern/southern-most rings respectively.""" -average_on_poles(A::AbstractGrid{NF},rings::Vector{<:UnitRange}) where NF = average_on_poles(NF,A,rings) - -function average_on_poles( ::Type{NF}, - A::AbstractGrid, +$(TYPEDSIGNATURES) +Computes the average at the North and South pole from a given grid `A` and it's precomputed +ring indices `rings`. The North pole average is an equally weighted average of all grid points +on the northern-most ring. Similar for the South pole.""" +function average_on_poles( A::AbstractGrid{NF}, rings::Vector{<:UnitRange{<:Integer}} ) where {NF<:AbstractFloat} @@ -424,6 +420,19 @@ function average_on_poles( ::Type{NF}, return convert(NF,A_northpole), convert(NF,A_southpole) end +""" +$(TYPEDSIGNATURES) +Method for `A::Abstract{T<:Integer}` which rounds the averaged values +to return the same number format `NF`.""" +function average_on_poles( A::AbstractGrid{NF}, + rings::Vector{<:UnitRange{<:Integer}} + ) where {NF<:Integer} + + A_northpole = mean(view(A,rings[1])) # average of all grid points around the north pole + A_southpole = mean(view(A,rings[end])) # same for south pole + return round(NF,A_northpole), round(NF,A_southpole) +end + """ $(TYPEDSIGNATURES) The bilinear average of a,b,c,d which are values at grid points @@ -446,19 +455,20 @@ longitude/x-coordinate. See schematic: 0...............1 # fraction of distance Δcd between c,d ``` ^ fraction of distance Δy between a-b and c-d.""" -function anvil_average( a::NF, # top left value - b::NF, # top right value - c::NF, # bottom left value - d::NF, # bottom right value - Δab::Real, # fraction of distance between a,b ∈ [0,1) - Δcd::Real, # fraction of distance between c,d ∈ [0,1) - Δy::Real, # fraction of distance between ab,cd ∈ [0,1) - ) where NF +function anvil_average( + a, # top left value + b, # top right value + c, # bottom left value + d, # bottom right value + Δab, # fraction of distance between a,b ∈ [0,1) + Δcd, # fraction of distance between c,d ∈ [0,1) + Δy, # fraction of distance between ab,cd ∈ [0,1) + ) # the type of the weights is ::Real, but the following is written such that # always NF (the type of the data values a,b,c,d) is returned - ab_average = a + (b-a)*convert(NF,Δab) # a for Δab=0, b for Δab=1 - cd_average = c + (d-c)*convert(NF,Δcd) # c for Δab=0, b for Δab=1 - abcd_average = ab_average + (cd_average-ab_average)*convert(NF,Δy) + ab_average = a + (b-a)*Δab # a for Δab=0, b for Δab=1 + cd_average = c + (d-c)*Δcd # c for Δab=0, b for Δab=1 + abcd_average = ab_average + (cd_average-ab_average)*Δy return abcd_average end \ No newline at end of file diff --git a/src/RingGrids/show.jl b/src/RingGrids/show.jl index 85900130f..6d3545d2d 100644 --- a/src/RingGrids/show.jl +++ b/src/RingGrids/show.jl @@ -4,4 +4,30 @@ function Base.array_summary(io::IO, grid::AbstractGrid, inds::Tuple{Vararg{Base.OneTo}}) print(io, Base.dims2string(length.(inds)), ", $(get_nlat(grid))-ring ") Base.showarg(io, grid, true) +end + +function plot(A::AbstractGrid;title::String="$(get_nlat(A))-ring $(typeof(A))") + A_full = interpolate(full_grid(typeof(A)),A.nlat_half,A) + plot(A_full;title) +end + +function plot(A::AbstractFullGrid;title::String="$(get_nlat(A))-ring $(typeof(A))") + + A_matrix = Matrix(A) + nlon,nlat = size(A_matrix) + A_view = view(A_matrix,:,nlat:-1:1) + + plot_kwargs = pairs(( xlabel="˚E", + xfact=360/(nlon-1), + ylabel="˚N", + yfact=180/(nlat-1), + yoffset=-90, + title=title, + colormap=:viridis, + compact=true, + colorbar=true, + width=60, + height=30)) + + UnicodePlots.heatmap(A_view';plot_kwargs...) end \ No newline at end of file diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index 194890462..63ba8ffa4 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -66,7 +66,8 @@ export LowerTriangularMatrix, OctahedralGaussianGrid, OctahedralClenshawGrid, HEALPixGrid, - OctaHEALPixGrid + OctaHEALPixGrid, + plot export Leapfrog diff --git a/src/dynamics/prognostic_variables.jl b/src/dynamics/prognostic_variables.jl index b23e3efc3..f8becf44e 100644 --- a/src/dynamics/prognostic_variables.jl +++ b/src/dynamics/prognostic_variables.jl @@ -354,24 +354,7 @@ get_humidity(progn::PrognosticVariables; kwargs...) = get_var(progn, :humid; kwa get_pressure(progn::PrognosticVariables; lf::Integer=1) = progn.surface.timesteps[lf].pres function Base.show(io::IO, P::PrognosticVariables) - - ζ = P.layers[end].timesteps[1].vor # create a view on vorticity - ζ_grid = Matrix(gridded(ζ)) # to grid space - ζ_grid = ζ_grid[:,end:-1:1] # flip latitudes - - nlon,nlat = size(ζ_grid) - - plot_kwargs = pairs(( xlabel="˚E", - xfact=360/(nlon-1), - ylabel="˚N", - yfact=180/(nlat-1), - yoffset=-90, - title="Surface relative vorticity", - colormap=:viridis, - compact=true, - colorbar=true, - width=60, - height=30)) - - print(io,UnicodePlots.heatmap(ζ_grid';plot_kwargs...)) + ζ = P.layers[end].timesteps[1].vor # create a view on surface relative vorticity + ζ_grid = gridded(ζ) # to grid space + print(io,plot(ζ_grid,title="Surface relative vorticity")) end \ No newline at end of file From 97e35640e0319ce4720437c98a3b0353e2294eea Mon Sep 17 00:00:00 2001 From: Milan Date: Sun, 4 Jun 2023 21:53:10 -0400 Subject: [PATCH 8/8] @example for ringgrids doc --- docs/make.jl | 4 +- docs/src/ringgrids.md | 328 ++++++++++----------------------- src/RingGrids/interpolation.jl | 5 +- 3 files changed, 99 insertions(+), 238 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 99e492963..a31b57dd8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,8 +1,8 @@ using Documenter, SpeedyWeather makedocs( - format = Documenter.HTML( - prettyurls = get(ENV, "CI", nothing) == "true"), + format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing)=="true", + ansicolor=true), sitename="SpeedyWeather.jl", authors="M Klöwer and SpeedyWeather contributors", modules=[SpeedyWeather], diff --git a/docs/src/ringgrids.md b/docs/src/ringgrids.md index 1eba388d9..2b57f73b4 100644 --- a/docs/src/ringgrids.md +++ b/docs/src/ringgrids.md @@ -19,89 +19,62 @@ The following explanation of how to use these can be mostly applied to any of th there are certain functions that are not defined, e.g. the full grids can be trivially converted to a `Matrix` but not the `OctahedralGaussianGrid`. +!!! note "What is a ring?" + We use the term *ring*, short for *iso-latitude ring*, to refer to a sequence of grid points + that all share the same latitude. A latitude-longitude grid is a ring grid, as it organises + its grid-points into rings. However, other grids, like the + [cubed-sphere](https://en.wikipedia.org/wiki/Quadrilateralized_spherical_cube) + are not based on iso-latitude rings. SpeedyWeather.jl only works with ring grids + because its a requirement for the [Spherical Harmonic Transform](@ref) to be efficient. + See [Grids](@ref). + ## Creating data on a RingGrid -Every grid in RingGrids has a `data` field, which is a vector containing the data on the grid. +Every `grid` in RingGrids has a `grid.data` field, which is a vector containing the data on the grid. The grid points are unravelled west to east then north to south, meaning that it starts at 90˚N and 0˚E then walks eastward for 360˚ before jumping on the next latitude ring further south, this way circling around the sphere till reaching the south pole. This may also be called _ring order_. -Data on a 96x48 `Matrix` which follows this ring order can be put on a `FullGaussianGrid` like so -```julia -julia> map -96×48 Matrix{Float32}: - -0.219801 -0.519598 1.2066 … 0.81304 -1.16023 1.0353 - -0.55615 -1.05712 -0.227948 -2.06369 1.10353 1.60918 - 0.446913 -0.856431 1.58896 -1.0894 -0.894315 0.632353 - 0.445915 -0.107201 -0.577785 0.574784 -0.825049 1.29677 - 1.194 -0.353374 1.30581 0.465554 0.358457 -0.726567 - 1.28693 1.43997 0.691283 … -0.330544 -0.267588 0.181308 - ⋮ ⋱ ⋮ - -0.432703 0.17233 0.89222 0.888913 1.32787 -0.248779 - -0.404498 0.127172 -0.64237 0.127979 -1.55253 -2.00749 - -0.857746 -0.433251 -0.468293 1.09825 -0.291169 1.07452 - 0.375367 -0.218278 0.492855 -0.287976 0.878996 -1.19745 - -0.0619525 -0.129129 -1.35502 … 0.0824819 0.481736 0.845638 +Data in a `Matrix` which follows this ring order can be put on a `FullGaussianGrid` like so +```@example ringgrids +using SpeedyWeather.RingGrids +map = randn(Float32,8,4) +``` -julia> grid = FullGaussianGrid(map) -4608-element, 48-ring FullGaussianGrid{Float32}: - -0.21980117 - -0.5561496 - 0.44691312 - 0.4459149 - 1.1940043 - 1.2869298 - ⋮ - -0.24877885 - -2.007495 - 1.0745221 - -1.197454 - 0.84563845 +```@example ringgrids +grid = FullGaussianGrid(map) ``` -A full Gaussian grid has always ``2N`` x ``N`` grid points, but a `FullClenshawGrid` has ``2N`` x ``N-1``, if those dimensions don't match, the creation will throw an error. To reobtain the data from a grid, you can access its `data` field -which returns a normal `Vector` -```julia -julia> grid.data -4608-element Vector{Float32}: - -0.21980117 - -0.5561496 - 0.44691312 - 0.4459149 - 1.1940043 - 1.2869298 - ⋮ - -0.24877885 - -2.007495 - 1.0745221 - -1.197454 - 0.84563845 +A full Gaussian grid has always ``2N`` x ``N`` grid points, but a `FullClenshawGrid` has ``2N`` x ``N-1``, +if those dimensions don't match, the creation will throw an error. To reobtain the data from a grid, +you can access its `data` field which returns a normal `Vector` +```@example ringgrids +grid.data ``` Which can be reshaped to reobtain `map` from above. Alternatively you can `Matrix(grid)` to do this in one step -```julia -julia> map == Matrix(FullGaussianGrid(map)) -true +```@example ringgrids +map == Matrix(FullGaussianGrid(map)) ``` You can also use `zeros`,`ones`,`rand`,`randn` to create a grid, whereby `nlat_half`, i.e. the number of latitude rings on one hemisphere, Equator included, is used as a resolution parameter and here as a second argument. -```julia -julia> nlat_half = 4 -julia> grid = randn(OctahedralGaussianGrid{Float16},nlat_half) -208-element, 8-ring OctahedralGaussianGrid{Float16}: - -1.868 - 0.493 - 0.3142 - 1.871 - 1.349 - 0.623 - ⋮ - 1.064 - 0.4346 - -0.641 - 0.1445 - 0.3909 +```@example ringgrids +nlat_half = 4 +grid = randn(OctahedralGaussianGrid{Float16},nlat_half) ``` and any element type `T` can be used for `OctahedralGaussianGrid{T}` and similar for other grid types. +## Visualising RingGrid data + +As only the full grids can be reshaped into a matrix, the underyling data structure of any `AbstractGrid` +is a vector. As shown in the examples above, one can therefore inspect the data as if it was a vector. +But as that data has, through its `<:AbstractGrid` type, all the geometric information available to plot +it on a map, RingGrids also exports `plot` function, +based on [UnicodePlots](https://github.com/JuliaPlots/UnicodePlots.jl)' `heatmap`. +```@example ringgrids +nlat_half = 24 +grid = randn(OctahedralGaussianGrid,nlat_half) +plot(grid) +``` + ## Indexing RingGrids All RingGrids have a single index `ij` which follows the ring order. While this is obviously not super @@ -109,20 +82,12 @@ exciting here are some examples how to make better use of the information that t We obtain the latitudes of the rings of a grid by calling `get_latd` (`get_lond` is only defined for full grids, or use `get_latdlonds` for latitudes, longitudes per grid point not per ring) -```julia -julia> latd = get_latd(grid) -8-element Vector{Float64}: - 73.79921362856324 - 52.81294318999426 - 31.704091745007943 - 10.569882312576093 - -10.569882312576093 - -31.704091745007943 - -52.81294318999426 - -73.79921362856324 +```@example ringgrids +grid = randn(OctahedralClenshawGrid,5) +latd = get_latd(grid) ``` Now we could calculate Coriolis and add it on the grid as follows -```julia +```@example ringgrids rotation = 7.29e-5 # angular frequency of Earth's rotation [rad/s] coriolis = 2rotation*sind.(latd) # vector of coriolis parameters per latitude ring @@ -139,7 +104,7 @@ end loop over all in-ring indices `i` (changing longitudes) to do something on the grid. Something similar can be done to scale/unscale with the cosine of latitude for example. We can always loop over all grid-points like so -```julia +```@example ringgrids for ij in eachgridpoint(grid) grid[ij] end @@ -154,53 +119,16 @@ above. This information generally can also be used to interpolate our data from request an interpolated value on some coordinates. Using our data on `grid` which is an `OctahedralGaussianGrid` from above we can use the `interpolate` function to get it onto a `FullGaussianGrid` (or any other grid for purpose) -```julia -julia> grid -208-element, 8-ring OctahedralGaussianGrid{Float16}: - -1.868 - 0.493 - 0.3142 - 1.871 - 1.349 - 0.623 - ⋮ - 1.064 - 0.4346 - -0.641 - 0.1445 - 0.3909 - -julia> interpolate(FullGaussianGrid,grid) -128-element, 8-ring FullGaussianGrid{Float64}: - -1.8681640625 - 0.4482421875 - 1.0927734375 - 1.4794921875 - 0.623046875 - -0.6435546875 - ⋮ - -0.57763671875 - 1.064453125 - 0.16552734375 - -0.248291015625 - 0.329345703125 +```@example ringgrids +grid = randn(OctahedralGaussianGrid{Float32},4) ``` -By default this will linearly interpolate (it's an anvil interpolator, see below) onto a grid with the same +```@example ringgrids +interpolate(FullGaussianGrid,grid) +``` +By default this will linearly interpolate (it's an [Anvil interpolator](@ref), see below) onto a grid with the same `nlat_half`, but we can also coarse-grain or fine-grain by specifying `nlat_half` directly as 2nd argument -```julia -julia> interpolate(FullGaussianGrid,6,grid) -288-element, 12-ring FullGaussianGrid{Float64}: - -1.248046875 - 0.08984375 - 0.2763671875 - 0.76513671875 - 1.1767578125 - ⋮ - 0.26416015625 - -0.295166015625 - -0.271728515625 - 0.0511474609375 - 0.0814208984375 +```@example ringgrids +interpolate(FullGaussianGrid,6,grid) ``` So we got from an `8-ring OctahedralGaussianGrid{Float16}` to a `12-ring FullGaussianGrid{Float64}`, so it did a conversion from `Float16` to `Float64` on the fly too, because the default precision is `Float64` unless @@ -208,19 +136,13 @@ specified. `interpolate(FullGaussianGrid{Float16},6,grid)` would have interpolat `Float16`. One can also interpolate onto a give cordinate ˚N, ˚E like so -```julia -julia> interpolate(30.0,10.0,grid) -1-element Vector{Float16}: - 0.9297 +```@example ringgrids +interpolate(30.0,10.0,grid) ``` we interpolated the data from `grid` onto 30˚N, 10˚E. To do this simultaneously for many coordinates they can be packed into a vector too -```julia -julia> interpolate([30.0,40.0,50.0],[10.0,10.0,10.0],grid) -3-element Vector{Float16}: - 0.9297 - 0.08887 - -0.929 +```@example ringgrids +interpolate([30.0,40.0,50.0],[10.0,10.0,10.0],grid) ``` which returns the data on `grid` at 30˚N, 40˚N, 50˚N, and 10˚E respectively. Note how the interpolation here retains the element type of `grid`. @@ -247,138 +169,74 @@ to grid B many times, you can just reuse the same interpolator. If you want to c of the output grid but its total number of points remain constants then you can update the locator inside the interpolator and only else you will need to create a new interpolator. Let's look at this in practice. Say we have two grids an want to interpolate between them -```julia -julia> grid_in = rand(HEALPixGrid,4) -julia> grid_out = zeros(FullClenshawGrid,6) -julia> interp = RingGrids.interpolator(grid_out,grid_in) +```@example ringgrids +grid_in = rand(HEALPixGrid,4) +grid_out = zeros(FullClenshawGrid,6) +interp = RingGrids.interpolator(grid_out,grid_in) ``` Now we have created an interpolator `interp` which knows about the geometry where to interpolate *from* and the coordinates there to interpolate *to*. It is also initialized, meaning it has precomputed the indices to of `grid_in` that are supposed to be used. It just does not know about the data of `grid_in` (and neither of `grid_out` which will be overwritten anyway). We can now do -```julia -julia> interpolate!(grid_out,grid_in,interp); -julia> grid_out -264-element, 11-ring FullClenshawGrid{Float64}: - 0.47698810225785837 - 0.49923033302273034 - 0.5214725637876022 - 0.5437147945524742 - ⋮ - 0.6277318221906577 - 0.5934538182075797 - 0.6009226488782581 - 0.6083914795489366 +```@example ringgrids +interpolate!(grid_out,grid_in,interp) +grid_out ``` -which is identical to `interpolate(grid_out,grid_in)` but you can reuse `interp` with more data. +which is identical to `interpolate(grid_out,grid_in)` but you can reuse `interp` for other data. The interpolation can also handle various element types (the interpolator `interp` does not have to be updated for this either) -```julia -julia> grid_out = zeros(FullClenshawGrid{Float16},6); -julia> interpolate!(grid_out,grid_in,interp) -julia> grid_out -264-element, 11-ring FullClenshawGrid{Float16}: - 0.477 - 0.4993 - 0.5215 - 0.544 - 0.5493 - 0.555 - ⋮ - 0.662 - 0.628 - 0.5933 - 0.601 - 0.6084 +```@example ringgrids +grid_out = zeros(FullClenshawGrid{Float16},6); +interpolate!(grid_out,grid_in,interp) +grid_out ``` -and we have converted data from a `HEALPixGrid{Float64}` (which is always default if not specified) +and we have converted data from a `HEALPixGrid{Float64}` (`Float64` is always default if not specified) to a `FullClenshawGrid{Float16}` including the type conversion Float64-Float16 on the fly. Technically there are three data types and their combinations possible: The input data will come with a type, the output array has an element type and the interpolator has precomputed weights with a given type. Say we want to go from Float16 data on an `OctahedralGaussianGrid` to Float16 on a `FullClenshawGrid` but using Float32 precision for the interpolation itself, we would do this by -```julia -julia> grid_in = randn(OctahedralGaussianGrid{Float16},24) -julia> grid_out = zeros(FullClenshawGrid{Float16},24) -julia> interp = RingGrids.interpolator(Float32,grid_out,grid_in) -julia> interpolate!(grid_out,grid_in,interp) -julia> grid_out -4512-element, 47-ring FullClenshawGrid{Float16}: - -0.954 - -0.724 - -0.4941 - -0.264 - -0.03433 - 0.1796 - ⋮ - -0.5703 - -0.3481 - -0.07666 - 0.1958 - 0.467 +```@example ringgrids +grid_in = randn(OctahedralGaussianGrid{Float16},24) +grid_out = zeros(FullClenshawGrid{Float16},24) +interp = RingGrids.interpolator(Float32,grid_out,grid_in) +interpolate!(grid_out,grid_in,interp) +grid_out ``` As a last example we want to illustrate a situation where we would always want to interplate onto 10 coordinates, but their locations may change. In order to avoid recreating an interpolator object we would do (`AnvilInterpolator` is described in [Anvil interpolator](@ref)) -```julia -julia> npoints = 10 # number of coordinates to interpolate onto -julia> interp = AnvilInterpolator(Float32,HEALPixGrid,24,npoints) +```@example ringgrids +npoints = 10 # number of coordinates to interpolate onto +interp = AnvilInterpolator(Float32,HEALPixGrid,24,npoints) ``` with the first argument being the number format used during interpolation, then the input grid type, its resolution in terms of `nlat_half` and then the number of points to interpolate onto. However, `interp` is not yet initialized as it does not know about the destination coordinates yet. Let's define them, but note that we already decided there's only 10 of them above. -```julia -julia> latds = collect(0.0:5.0:45.0) -10-element Vector{Float64}: - 0.0 - 5.0 - ⋮ - 40.0 - 45.0 - -julia> londs = collect(-10.0:2.0:8.0) -10-element Vector{Float64}: - -10.0 - -8.0 - -6.0 - ⋮ - 6.0 - 8.0 +```@example ringgrids +latds = collect(0.0:5.0:45.0) +londs = collect(-10.0:2.0:8.0) +nothing # hide ``` now we can update the locator inside our interpolator as follows -```julia -julia> RingGrids.update_locator!(interp,latds,londs) +```@example ringgrids +RingGrids.update_locator!(interp,latds,londs) ``` With data matching the input from above, a `nlat_half=24` HEALPixGrid, and allocate 10-element output vector -```julia -julia> output_vec = zeros(10) -julia> grid_input = rand(HEALPixGrid,24) +```@example ringgrids +output_vec = zeros(10) +grid_input = rand(HEALPixGrid,24) +nothing # hide ``` we can use the interpolator as follows -```julia -julia> interpolate!(output_vec,grid_input,interp) -10-element Vector{Float64}: - 0.3182548251299291 - 0.7499448926757676 - 0.25733825675836064 - ⋮ - 0.2949249541923441 - 0.6690698461409016 - 0.6159433564856793 +```@example ringgrids +interpolate!(output_vec,grid_input,interp) ``` which is the approximately the same as doing it directly without creating an interpolator first and updating its locator -```julia -julia> interpolate(latds,londs,grid_input) -10-element Vector{Float64}: - 0.31825482404891603 - 0.7499448923165136 - 0.25733824520344434 - ⋮ - 0.294924962125593 - 0.6690698486360254 - 0.6159433558779497 +```@example ringgrids +interpolate(latds,londs,grid_input) ``` but allows for a reuse of the interpolator. Note that the two output arrays are not exactly identical because we manually set our interpolator `interp` to use `Float32` for the interplation whereas the default is `Float64`. diff --git a/src/RingGrids/interpolation.jl b/src/RingGrids/interpolation.jl index a2d10a105..ea108f379 100644 --- a/src/RingGrids/interpolation.jl +++ b/src/RingGrids/interpolation.jl @@ -174,7 +174,10 @@ function interpolator( Aout::AbstractGrid, end ## FUNCTIONS -interpolate(latd::Real,lond::Real,A::AbstractGrid) = interpolate([latd],[lond],A) +function interpolate(latd::Real,lond::Real,A::AbstractGrid) + Ai = interpolate([latd],[lond],A) + return Ai[1] +end function interpolate( latds::Vector{NF}, # latitudes to interpolate onto (90˚N...-90˚N) londs::Vector{NF}, # longitudes to interpolate into (0˚...360˚E)