Soil moisture and vegetation

src/SpeedyWeather.jl

@@ -26,6 +26,9 @@ import BitInformation: round, round!
 import UnicodePlots
 import ProgressMeter
 
+# to avoid a `using Dates` to pass on DateTime arguments
+export DateTime
+
 # EXPORT MONOLITHIC INTERFACE TO SPEEDY
 export  run_speedy,
         run_speedy!,
@@ -183,7 +186,7 @@ include("physics/pretty_printing.jl")
 
 # OCEAN AND LAND
 include("dynamics/ocean.jl")
-include("dynamics/land.jl")
+include("physics/land.jl")
 
 # MODELS
 include("dynamics/models.jl")

src/abstract_types.jl

@@ -21,7 +21,6 @@ abstract type InitialConditions end
 abstract type AbstractOrography{NF,Grid} end
 abstract type AbstractLandSeaMask{NF,Grid} end
 abstract type AbstractAlbedo{NF,Grid} end
-abstract type AbstractVegetation{NF,Grid} end
 
 # ATMOSPHERIC COLUMN FOR PARAMETERIZATIONS
 abstract type AbstractColumnVariables{NF} end

src/dynamics/diagnostic_variables.jl

@@ -138,6 +138,7 @@ Base.@kwdef struct SurfaceVariables{NF<:AbstractFloat,Grid<:AbstractGrid{NF}}
 
     precip_large_scale::Grid = zeros(Grid,nlat_half)    # large scale precipitation (for output)
     precip_convection::Grid = zeros(Grid,nlat_half)     # convective precipitation (for output)
+    soil_moisture_availability::Grid = zeros(Grid,nlat_half)
 end
 
 # generator function based on a SpectralGrid

src/dynamics/models.jl

@@ -55,7 +55,7 @@ $(TYPEDSIGNATURES)
 Calls all `initialize!` functions for components of `model`,
 except for `model.output` and `model.feedback` which are always called
 at in `time_stepping!`."""
-function initialize!(model::Barotropic)
+function initialize!(model::Barotropic;time::DateTime = DEFAULT_DATE)
     (;spectral_grid) = model
 
     spectral_grid.nlev > 1 && @warn "Only nlev=1 supported for BarotropicModel, \
@@ -68,6 +68,7 @@ function initialize!(model::Barotropic)
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid,model)
     initialize!(prognostic_variables,model.initial_conditions,model)
+    prognostic_variables.clock.time = time       # set the time
 
     diagnostic_variables = DiagnosticVariables(spectral_grid,model)
     return Simulation(prognostic_variables,diagnostic_variables,model)
@@ -115,7 +116,7 @@ $(TYPEDSIGNATURES)
 Calls all `initialize!` functions for components of `model`,
 except for `model.output` and `model.feedback` which are always called
 at in `time_stepping!` and `model.implicit` which is done in `first_timesteps!`."""
-function initialize!(model::ShallowWater)
+function initialize!(model::ShallowWater;time::DateTime = DEFAULT_DATE)
     (;spectral_grid) = model
 
     spectral_grid.nlev > 1 && @warn "Only nlev=1 supported for ShallowWaterModel, \
@@ -129,6 +130,7 @@ function initialize!(model::ShallowWater)
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid,model)
     initialize!(prognostic_variables,model.initial_conditions,model)
+    prognostic_variables.clock.time = time       # set the time
 
     diagnostic_variables = DiagnosticVariables(spectral_grid,model)
     return Simulation(prognostic_variables,diagnostic_variables,model)
@@ -153,7 +155,7 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
     # BOUNDARY CONDITIONS
     land_sea_mask::AbstractLandSeaMask{NF} = LandSeaMask(spectral_grid)
     ocean::AbstractOcean{NF} = SeasonalOceanClimatology(spectral_grid)
-    land::AbstractLand{NF} = SeasonalLandClimatology(spectral_grid)
+    land::AbstractLand{NF} = SeasonalLandTemperature(spectral_grid)
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
@@ -189,7 +191,7 @@ $(TYPEDSIGNATURES)
 Calls all `initialize!` functions for components of `model`,
 except for `model.output` and `model.feedback` which are always called
 at in `time_stepping!` and `model.implicit` which is done in `first_timesteps!`."""
-function initialize!(model::PrimitiveDry)
+function initialize!(model::PrimitiveDry;time::DateTime = DEFAULT_DATE)
     (;spectral_grid) = model
 
     # numerics (implicit is initialized later)
@@ -209,10 +211,12 @@ function initialize!(model::PrimitiveDry)
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid,model)
     initialize!(prognostic_variables,model.initial_conditions,model)
-
-    (;time) = prognostic_variables.clock
-    initialize!(prognostic_variables.ocean,time,model)
-    initialize!(prognostic_variables.land,time,model)
+    (;clock) = prognostic_variables
+    clock.time = time       # set the time
+
+    # initialize ocean and land and synchronize clocks
+    initialize!(prognostic_variables.ocean,clock.time,model)
+    initialize!(prognostic_variables.land,clock.time,model)
 
     diagnostic_variables = DiagnosticVariables(spectral_grid,model)
     return Simulation(prognostic_variables,diagnostic_variables,model)
@@ -237,7 +241,9 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
     # BOUNDARY CONDITIONS
     land_sea_mask::AbstractLandSeaMask{NF} = LandSeaMask(spectral_grid)
     ocean::AbstractOcean{NF} = SeasonalOceanClimatology(spectral_grid)
-    land::AbstractLand{NF} = SeasonalLandClimatology(spectral_grid)
+    land::AbstractLand{NF} = SeasonalLandTemperature(spectral_grid)
+    soil::AbstractSoil{NF} = SeasonalSoilMoisture(spectral_grid)
+    vegetation::AbstractVegetation{NF} = VegetationClimatology(spectral_grid)
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
@@ -276,7 +282,7 @@ $(TYPEDSIGNATURES)
 Calls all `initialize!` functions for components of `model`,
 except for `model.output` and `model.feedback` which are always called
 at in `time_stepping!` and `model.implicit` which is done in `first_timesteps!`."""
-function initialize!(model::PrimitiveWet)
+function initialize!(model::PrimitiveWet;time::DateTime = DEFAULT_DATE)
     (;spectral_grid) = model
 
     # numerics (implicit is initialized later)
@@ -287,6 +293,8 @@ function initialize!(model::PrimitiveWet)
     initialize!(model.land_sea_mask)
     initialize!(model.ocean)
     initialize!(model.land)
+    initialize!(model.soil)
+    initialize!(model.vegetation)
 
     # parameterizations
     initialize!(model.boundary_layer_drag,model)
@@ -297,10 +305,12 @@ function initialize!(model::PrimitiveWet)
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid,model)
     initialize!(prognostic_variables,model.initial_conditions,model)
+    (;clock) = prognostic_variables
+    clock.time = time       # set the time
 
-    (;time) = prognostic_variables.clock
-    initialize!(prognostic_variables.ocean,time,model)
-    initialize!(prognostic_variables.land,time,model)
+    # initialize ocean and land and synchronize clocks
+    initialize!(prognostic_variables.ocean,clock.time,model)
+    initialize!(prognostic_variables.land,clock.time,model)
 
     diagnostic_variables = DiagnosticVariables(spectral_grid,model)
     return Simulation(prognostic_variables,diagnostic_variables,model)

src/dynamics/ocean.jl

@@ -21,6 +21,9 @@ Base.@kwdef struct SeasonalOceanClimatology{NF,Grid<:AbstractGrid{NF}} <: Abstra
     "filename of sea surface temperatures"
     file::String = "sea_surface_temperature.nc"
 
+    "Variable name in netcdf file"
+    varname::String = "sst"
+
     "Grid the sea surface temperature file comes on"
     file_Grid::Type{<:AbstractGrid} = FullGaussianGrid
 
@@ -39,33 +42,42 @@ function SeasonalOceanClimatology(SG::SpectralGrid;kwargs...)
 end
 
 function initialize!(ocean::SeasonalOceanClimatology{NF,Grid}) where {NF,Grid}
+    load_monthly_climatology!(ocean.monthly_temperature,ocean)
+end
+
+function load_monthly_climatology!( 
+    monthly::Vector{Grid},
+    scheme;
+    varname::String = scheme.varname
+) where {Grid<:AbstractGrid}
+
     # LOAD NETCDF FILE
-    if ocean.path == "SpeedyWeather.jl/input_data"
-        path = joinpath(@__DIR__,"../../input_data",ocean.file)
+    if scheme.path == "SpeedyWeather.jl/input_data"
+        path = joinpath(@__DIR__,"../../input_data",scheme.file)
     else
-        path = joinpath(ocean.path,ocean.file)
+        path = joinpath(scheme.path,scheme.file)
     end
     ncfile = NCDataset(path)
 
     # create interpolator from grid in file to grid used in model
     nx, ny = ncfile.dim["lon"], ncfile.dim["lat"]
     npoints = nx*ny
-    NF_file = typeof(ncfile["sst"].attrib["_FillValue"])
-    sst = ocean.file_Grid(zeros(NF_file,npoints))
-    interp = RingGrids.interpolator(NF,ocean.monthly_temperature[1],sst)
+    NF_file = typeof(ncfile[varname].attrib["_FillValue"])
+    grid = scheme.file_Grid(zeros(NF_file,npoints))
+    interp = RingGrids.interpolator(Float32,monthly[1],grid)
 
-    # interpolate and store in ocean
+    # interpolate and store in monthly
     for month in 1:12
-        sst_this_month = ncfile["sst"][:,:,month]
+        this_month = ncfile[varname][:,:,month]
         ij = 0
         for j in 1:ny
             for i in 1:nx
                 ij += 1
-                x = sst_this_month[i,j]
-                sst[ij] = ismissing(x) ? ocean.missing_value : x
+                x = this_month[i,j]
+                grid[ij] = ismissing(x) ? scheme.missing_value : x
             end
         end
-        interpolate!(ocean.monthly_temperature[month],sst,interp)
+        interpolate!(monthly[month],grid,interp)
     end
     return nothing
 end

src/dynamics/prognostic_variables.jl

@@ -136,11 +136,11 @@ Base.@kwdef mutable struct PrognosticVariablesLand{NF<:AbstractFloat,Grid<:Abstr
     "Snow depth [m]"
     const snow_depth::Grid = zeros(Grid,nlat_half)
 
-    "Soil wetness layer 1, volume fraction [1]"
-    const soil_wetness_layer1::Grid = zeros(Grid,nlat_half)
+    "Soil moisture layer 1, volume fraction [1]"
+    const soil_moisture_layer1::Grid = zeros(Grid,nlat_half)
 
-    "Soil wetness layer 2, volume fraction [1]"
-    const soil_wetness_layer2::Grid = zeros(Grid,nlat_half)
+    "Soil moisture layer 2, volume fraction [1]"
+    const soil_moisture_layer2::Grid = zeros(Grid,nlat_half)
 end
 
 # generator function based on a SpectralGrid

src/dynamics/time_integration.jl

@@ -221,12 +221,17 @@ function timestep!( progn::PrognosticVariables{NF}, # all prognostic variables
     model.feedback.nars_detected && return nothing  # exit immediately if NaRs already present
     (;time) = progn.clock                           # current time
 
-    ocean_timestep!(progn.ocean,time,model)
-    land_timestep!(progn.land,time,model)
-
-    # switch on/off all physics
-    model.physics && parameterization_tendencies!(diagn,progn,time,model)
-    model.physics || zero_tendencies!(diagn)        # set tendencies to zero otherwise
+    if model.physics                                # switch on/off all physics parameterizations
+        # time step ocean (temperature and TODO sea ice) and land (temperature and soil moisture)
+        ocean_timestep!(progn.ocean,time,model)
+        land_timestep!(progn.land,time,model)
+        soil_moisture_availability!(diagn.surface,progn.land,model)
+
+        # calculate all parameterizations
+        parameterization_tendencies!(diagn,progn,time,model)
+    else                                            # set tendencies to zero otherwise for accumulators
+        zero_tendencies!(diagn)
+    end       
 
     dynamics_tendencies!(diagn,progn,model,lf2)         # dynamical core
     implicit_correction!(diagn,model.implicit,progn)    # semi-implicit time stepping corrections

src/physics/boundary_layer.jl

@@ -1,6 +1,8 @@
 """Concrete type that disables the boundary layer drag scheme."""
 struct NoBoundaryLayerDrag{NF} <: BoundaryLayerDrag{NF} end
 
+NoBoundaryLayerDrag(SG::SpectralGrid) = NoBoundaryLayerDrag{SG.NF}()
+
 """NoBoundaryLayer scheme does not need any initialisation."""
 function initialize!(   scheme::NoBoundaryLayerDrag,
                         model::PrimitiveEquation)
@@ -9,8 +11,7 @@ end
 
 """NoBoundaryLayer scheme just passes."""
 function boundary_layer_drag!(  column::ColumnVariables,
-                                scheme::NoBoundaryLayerDrag,
-                                model::PrimitiveEquation)
+                                scheme::NoBoundaryLayerDrag)
     return nothing
 end
 

src/physics/column_variables.jl

@@ -40,6 +40,7 @@ function get_column!(   C::ColumnVariables,
     # TODO skin = surface approximation for now
     C.skin_temperature_sea = P.ocean.sea_surface_temperature[ij]
     C.skin_temperature_land = P.land.land_surface_temperature[ij]
+    C.soil_moisture_availability = D.surface.soil_moisture_availability[ij]
 end
 
 """Recalculate ring index if not provided."""

src/physics/define_column.jl

@@ -53,6 +53,7 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     surface_wind_speed::NF = 0
     skin_temperature_sea::NF = 0
     skin_temperature_land::NF = 0
+    soil_moisture_availability::NF = 0
 
     # THERMODYNAMICS
     surface_air_density::NF = 0