Merge pull request #338 from SpeedyWeather/mk/condensation

SpeedyCondensation implemented

src/SpeedyWeather.jl

@@ -1,7 +1,6 @@
 module SpeedyWeather
 
 # STRUCTURE
-import Parameters: @with_kw
 using DocStringExtensions
 
 # NUMERICS
@@ -98,6 +97,9 @@ export  NoBoundaryLayer,
 export  NoVerticalDiffusion,
         VerticalLaplacian
 
+# Large scale condensation
+export  SpeedyCondensation
+
 # EXPORT STRUCTS
 export  DynamicsConstants,
         SpectralTransform,
@@ -166,14 +168,14 @@ 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
 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")
 

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

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
@@ -24,6 +24,9 @@ $(TYPEDFIELDS)"""
     "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
 
@@ -50,6 +53,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
 

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
@@ -31,6 +31,9 @@ $(TYPEDFIELDS)"""
     "= 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,

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
@@ -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

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"

src/dynamics/initial_conditions.jl

@@ -26,7 +26,7 @@ function allocate(
     return zeros(PrognosticVariables{NF},Model,trunc,nlev)
 end
 
-@with_kw struct StartFromRest <: InitialConditions 
+Base.@kwdef struct StartFromRest <: InitialConditions 
     pressure_on_orography::Bool = false
 end
 
@@ -46,10 +46,9 @@ end
 
 """Start with random vorticity as initial conditions
 $(TYPEDFIELDS)"""
-@with_kw struct StartWithRandomVorticity <: InitialConditions
-
-    "Power law the vorticity should be spectrally distributed by"
-    power_law::Float64 = -3
+Base.@kwdef struct StartWithRandomVorticity <: InitialConditions
+    "Power of the spectral distribution k^power"
+    power::Float64 = -3
 
     "(approximate) amplitude in [1/s], used as standard deviation of spherical harmonic coefficients"
     amplitude::Float64 = 1e-5
@@ -63,7 +62,7 @@ function initial_conditions!(   progn::PrognosticVariables{NF},
                                 model::ModelSetup) where NF
 
     lmax = progn.trunc+1
-    power = initial_conditions.power_law + 1    # +1 as power is summed of orders m
+    power = initial_conditions.power + 1    # +1 as power is summed of orders m
     ξ = randn(Complex{NF},lmax,lmax)*convert(NF,initial_conditions.amplitude)
 
     for progn_layer in progn.layers
@@ -82,7 +81,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
 
@@ -202,7 +201,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
 
@@ -343,15 +342,20 @@ 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
 
 """
 Restart from a previous SpeedyWeather.jl simulation via the restart file restart.jld2
 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()
 
@@ -453,40 +457,44 @@ 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
+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, pres_ref) = 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(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)
+
+    # 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
+    end
+end

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()
 
@@ -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)
@@ -192,7 +191,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 +208,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
@@ -230,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)
@@ -248,6 +247,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)

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"