diff --git a/input_data/albedo.nc b/input_data/albedo.nc
new file mode 100644
index 000000000..82d3abe8d
Binary files /dev/null and b/input_data/albedo.nc differ
diff --git a/input_data/sea_ice.nc b/input_data/sea_ice.nc
index 72bd0c49f..5ebf7aee3 100644
Binary files a/input_data/sea_ice.nc and b/input_data/sea_ice.nc differ
diff --git a/input_data/snow.nc b/input_data/snow.nc
index f37de794b..7fcccb742 100644
Binary files a/input_data/snow.nc and b/input_data/snow.nc differ
diff --git a/input_data/soil.nc b/input_data/soil.nc
deleted file mode 100644
index a40f7b3cb..000000000
Binary files a/input_data/soil.nc and /dev/null differ
diff --git a/input_data/soil_moisture.nc b/input_data/soil_moisture.nc
new file mode 100644
index 000000000..0a1619c2d
Binary files /dev/null and b/input_data/soil_moisture.nc differ
diff --git a/input_data/surface.nc b/input_data/surface.nc
deleted file mode 100644
index 007a09102..000000000
Binary files a/input_data/surface.nc and /dev/null differ
diff --git a/input_data/vegetation.nc b/input_data/vegetation.nc
new file mode 100644
index 000000000..bd78d364c
Binary files /dev/null and b/input_data/vegetation.nc differ
diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl
index 777662b09..0cf01a018 100644
--- a/src/SpeedyWeather.jl
+++ b/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")
diff --git a/src/abstract_types.jl b/src/abstract_types.jl
index 714c01ab8..248354a8f 100644
--- a/src/abstract_types.jl
+++ b/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
diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl
index e24bb1593..d48657756 100644
--- a/src/dynamics/diagnostic_variables.jl
+++ b/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
diff --git a/src/dynamics/land.jl b/src/dynamics/land.jl
deleted file mode 100644
index dc8f609ad..000000000
--- a/src/dynamics/land.jl
+++ /dev/null
@@ -1,110 +0,0 @@
-abstract type AbstractLand{NF,Grid} end
-
-function Base.show(io::IO,L::AbstractLand)
-    println(io,"$(typeof(L)) <: AbstractLand")
-    keys = propertynames(L)
-    print_fields(io,L,keys)
-end
-
-Base.@kwdef struct SeasonalLandClimatology{NF,Grid<:AbstractGrid{NF}} <: AbstractLand{NF,Grid}
-
-    "number of latitudes on one hemisphere, Equator included"
-    nlat_half::Int
-
-    # OPTIONS
-    "Time step used to update land surface temperatures"
-    Δt::Dates.Day = Dates.Day(3)
-
-    "path to the folder containing the land-sea mask file, pkg path default"
-    path::String = "SpeedyWeather.jl/input_data"
-
-    "filename of sea surface temperatures"
-    file::String = "land_surface_temperature.nc"
-
-    "Grid the sea surface temperature file comes on"
-    file_Grid::Type{<:AbstractGrid} = FullGaussianGrid
-
-    "The missing value in the data respresenting land"
-    missing_value::NF = NF(NaN)
-
-    # to be filled from file
-    "Monthly land surface temperatures [K], interpolated onto Grid"
-    monthly_temperature::Vector{Grid} = [zeros(Grid,nlat_half) for _ in 1:12]
-end
-
-# generator function
-function SeasonalLandClimatology(SG::SpectralGrid;kwargs...)
-    (;NF,Grid,nlat_half) = SG
-    return SeasonalLandClimatology{NF,Grid{NF}}(;nlat_half,kwargs...)
-end
-
-function initialize!(land::SeasonalLandClimatology{NF,Grid}) where {NF,Grid}
-    # LOAD NETCDF FILE
-    if land.path == "SpeedyWeather.jl/input_data"
-        path = joinpath(@__DIR__,"../../input_data",land.file)
-    else
-        path = joinpath(land.path,land.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["lst"].attrib["_FillValue"])
-    lst = land.file_Grid(zeros(NF_file,npoints))
-    interp = RingGrids.interpolator(NF,land.monthly_temperature[1],lst)
-
-    # interpolate and store in land
-    for month in 1:12
-        lst_this_month = ncfile["lst"][:,:,month]
-        ij = 0
-        for j in 1:ny
-            for i in 1:nx
-                ij += 1
-                x = lst_this_month[i,j]
-                lst[ij] = ismissing(x) ? land.missing_value : x
-            end
-        end
-        interpolate!(land.monthly_temperature[month],lst,interp)
-    end
-    return nothing
-end
-
-function initialize!(   land::PrognosticVariablesLand,
-                        time::DateTime,
-                        model::PrimitiveEquation)
-    land_timestep!(land,time,model,initialize=true)
-end
-
-# function barrier
-function land_timestep!(land::PrognosticVariablesLand,
-                        time::DateTime,
-                        model::PrimitiveEquation;
-                        initialize::Bool = false)
-    land_timestep!(land,time,model.land;initialize)
-end
-
-function land_timestep!(land::PrognosticVariablesLand{NF},
-                        time::DateTime,
-                        land_model::SeasonalLandClimatology;
-                        initialize::Bool = false) where NF
-
-    # escape immediately if Δt of land model hasn't passed yet
-    # unless the land hasn't been initialized yet
-    initialize || (time - land.time) < land_model.Δt && return nothing
-
-    # otherwise update land prognostic variables:
-    land.time = time
-    this_month = Dates.month(time)
-    next_month = (this_month % 12) + 1      # mod for dec 12 -> jan 1
-
-    # linear interpolation weight between the two months
-    # TODO check whether this shifts the climatology by 1/2 a month
-    weight = convert(NF,Dates.days(time-Dates.firstdayofmonth(time))/Dates.daysinmonth(time))
-
-    (;monthly_temperature) = land_model
-    @. land.land_surface_temperature = (1-weight) * monthly_temperature[this_month] +
-                                          weight  * monthly_temperature[next_month]
-
-    return nothing
-end
diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl
index d9c864759..9f51332d5 100644
--- a/src/dynamics/models.jl
+++ b/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)
diff --git a/src/dynamics/ocean.jl b/src/dynamics/ocean.jl
index f759e5383..5fa5d418e 100644
--- a/src/dynamics/ocean.jl
+++ b/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
diff --git a/src/dynamics/prognostic_variables.jl b/src/dynamics/prognostic_variables.jl
index 694f1c46d..d7d6c66ce 100644
--- a/src/dynamics/prognostic_variables.jl
+++ b/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
diff --git a/src/dynamics/time_integration.jl b/src/dynamics/time_integration.jl
index eb39e517a..c668d17fc 100644
--- a/src/dynamics/time_integration.jl
+++ b/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
diff --git a/src/physics/boundary_layer.jl b/src/physics/boundary_layer.jl
index cfe719b7f..c6fb27b0e 100644
--- a/src/physics/boundary_layer.jl
+++ b/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
 
diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl
index d4a485e66..9a1989ce7 100644
--- a/src/physics/column_variables.jl
+++ b/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."""
diff --git a/src/physics/define_column.jl b/src/physics/define_column.jl
index 77968dcef..efa5a10bc 100644
--- a/src/physics/define_column.jl
+++ b/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
diff --git a/src/physics/land.jl b/src/physics/land.jl
new file mode 100644
index 000000000..e5053c719
--- /dev/null
+++ b/src/physics/land.jl
@@ -0,0 +1,290 @@
+## LAND TEMPERATURE
+abstract type AbstractLand{NF,Grid} end
+
+function Base.show(io::IO,L::AbstractLand)
+    println(io,"$(typeof(L)) <: AbstractLand")
+    keys = propertynames(L)
+    print_fields(io,L,keys)
+end
+
+Base.@kwdef struct SeasonalLandTemperature{NF,Grid<:AbstractGrid{NF}} <: AbstractLand{NF,Grid}
+
+    "number of latitudes on one hemisphere, Equator included"
+    nlat_half::Int
+
+    # OPTIONS
+    "Time step used to update land surface temperatures"
+    Δt::Dates.Day = Dates.Day(3)
+
+    "path to the folder containing the land temperature file, pkg path default"
+    path::String = "SpeedyWeather.jl/input_data"
+
+    "filename of land surface temperatures"
+    file::String = "land_surface_temperature.nc"
+
+    "variable name in netcdf file"
+    varname::String = "lst"
+
+    "Grid the land surface temperature file comes on"
+    file_Grid::Type{<:AbstractGrid} = FullGaussianGrid
+
+    "The missing value in the data respresenting ocean"
+    missing_value::NF = NF(NaN)
+
+    # to be filled from file
+    "Monthly land surface temperatures [K], interpolated onto Grid"
+    monthly_temperature::Vector{Grid} = [zeros(Grid,nlat_half) for _ in 1:12]
+end
+
+# generator function
+function SeasonalLandTemperature(SG::SpectralGrid;kwargs...)
+    (;NF,Grid,nlat_half) = SG
+    return SeasonalLandTemperature{NF,Grid{NF}}(;nlat_half,kwargs...)
+end
+
+function initialize!(land::SeasonalLandTemperature{NF,Grid}) where {NF,Grid}
+    load_monthly_climatology!(land.monthly_temperature,land)
+end
+
+function initialize!(   land::PrognosticVariablesLand,
+                        time::DateTime,
+                        model::PrimitiveEquation)
+    land_timestep!(land,time,model.land,initialize=true)
+    model isa PrimitiveWet && soil_timestep!(land,time,model.soil)
+end
+
+# function barrier
+function land_timestep!(land::PrognosticVariablesLand,
+                        time::DateTime,
+                        model::PrimitiveEquation;
+                        initialize::Bool = false)
+    # the time step is dictated by the land "model" 
+    executed = land_timestep!(land,time,model.land;initialize)
+    executed && model isa PrimitiveWet && soil_timestep!(land,time,model.soil)
+end
+
+function land_timestep!(land::PrognosticVariablesLand{NF},
+                        time::DateTime,
+                        land_model::SeasonalLandTemperature;
+                        initialize::Bool = false) where NF
+
+    # escape immediately if Δt of land model hasn't passed yet
+    # unless the land hasn't been initialized yet
+    initialize || (time - land.time) < land_model.Δt && return false    # = executed
+
+    # otherwise update land prognostic variables:
+    land.time = time
+    this_month = Dates.month(time)
+    next_month = (this_month % 12) + 1      # mod for dec 12 -> jan 1
+
+    # linear interpolation weight between the two months
+    # TODO check whether this shifts the climatology by 1/2 a month
+    weight = convert(NF,Dates.days(time-Dates.firstdayofmonth(time))/Dates.daysinmonth(time))
+
+    (;monthly_temperature) = land_model
+    @. land.land_surface_temperature = (1-weight) * monthly_temperature[this_month] +
+                                          weight  * monthly_temperature[next_month]
+
+    return true # = executed
+end
+
+## SOIL MOISTURE
+abstract type AbstractSoil{NF,Grid} end
+
+function Base.show(io::IO,S::AbstractSoil)
+    println(io,"$(typeof(S)) <: AbstractSoil")
+    keys = propertynames(S)
+    print_fields(io,S,keys)
+end
+
+Base.@kwdef struct SeasonalSoilMoisture{NF,Grid<:AbstractGrid{NF}} <: AbstractSoil{NF,Grid}
+
+    "number of latitudes on one hemisphere, Equator included"
+    nlat_half::Int
+
+    # OPTIONS
+    "Depth of top soil layer [m]"
+    D_top::Float64 = 0.07
+
+    "Depth of root layer [m]"
+    D_root::Float64 = 0.21
+
+    "Soil wetness at field capacity [volume fraction]"
+    W_cap::Float64 = 0.3
+
+    "Soil wetness at wilting point [volume fraction]"
+    W_wilt::Float64 = 0.17
+
+    # READ CLIMATOLOGY FROM FILE
+    "path to the folder containing the soil moisture file, pkg path default"
+    path::String = "SpeedyWeather.jl/input_data"
+
+    "filename of soil moisture"
+    file::String = "soil_moisture.nc"
+
+    "variable name in netcdf file"
+    varname_layer1::String = "swl1"
+    varname_layer2::String = "swl2"
+
+    "Grid the soil moisture file comes on"
+    file_Grid::Type{<:AbstractGrid} = FullGaussianGrid
+
+    "The missing value in the data respresenting ocean"
+    missing_value::NF = NF(NaN)
+
+    # to be filled from file
+    "Monthly soil moisture volume fraction [1], top layer, interpolated onto Grid"
+    monthly_soil_moisture_layer1::Vector{Grid} = [zeros(Grid,nlat_half) for _ in 1:12]
+
+    "Monthly soil moisture volume fraction [1], 2nd layer, interpolated onto Grid"
+    monthly_soil_moisture_layer2::Vector{Grid} = [zeros(Grid,nlat_half) for _ in 1:12]
+end
+
+# generator function
+function SeasonalSoilMoisture(SG::SpectralGrid;kwargs...)
+    (;NF,Grid,nlat_half) = SG
+    return SeasonalSoilMoisture{NF,Grid{NF}}(;nlat_half,kwargs...)
+end
+
+function initialize!(soil::SeasonalSoilMoisture{NF,Grid}) where {NF,Grid}
+    load_monthly_climatology!(soil.monthly_soil_moisture_layer1,soil,varname=soil.varname_layer1)
+    load_monthly_climatology!(soil.monthly_soil_moisture_layer2,soil,varname=soil.varname_layer2)
+end
+
+function soil_timestep!(land::PrognosticVariablesLand{NF},
+                        time::DateTime,
+                        soil_model::SeasonalSoilMoisture) where NF
+
+    this_month = Dates.month(time)
+    next_month = (this_month % 12) + 1      # mod for dec 12 -> jan 1
+
+    # linear interpolation weight between the two months
+    # TODO check whether this shifts the climatology by 1/2 a month
+    weight = convert(NF,Dates.days(time-Dates.firstdayofmonth(time))/Dates.daysinmonth(time))
+
+    (;monthly_soil_moisture_layer1, monthly_soil_moisture_layer2) = soil_model
+    @. land.soil_moisture_layer1 = (1-weight) * monthly_soil_moisture_layer1[this_month] +
+                                        weight  * monthly_soil_moisture_layer1[next_month]
+    @. land.soil_moisture_layer2 = (1-weight) * monthly_soil_moisture_layer2[this_month] +
+                                        weight  * monthly_soil_moisture_layer2[next_month]
+
+    return nothing
+end
+
+## SOIL MOISTURE
+abstract type AbstractVegetation{NF,Grid} end
+
+function Base.show(io::IO,A::AbstractVegetation)
+    println(io,"$(typeof(A)) <: AbstractVegetation")
+    keys = propertynames(A)
+    print_fields(io,A,keys)
+end
+
+Base.@kwdef struct VegetationClimatology{NF,Grid<:AbstractGrid{NF}} <: AbstractVegetation{NF,Grid}
+
+    "number of latitudes on one hemisphere, Equator included"
+    nlat_half::Int
+
+    # OPTIONS
+    "Combine high and low vegetation factor, a in high + a*low [1]"
+    low_veg_factor::Float64 = 0.8
+
+    "path to the folder containing the soil moisture file, pkg path default"
+    path::String = "SpeedyWeather.jl/input_data"
+
+    "filename of soil moisture"
+    file::String = "vegetation.nc"
+
+    "variable name in netcdf file"
+    varname_vegh::String = "vegh"
+    varname_vegl::String = "vegl"
+
+    "Grid the soil moisture file comes on"
+    file_Grid::Type{<:AbstractGrid} = FullGaussianGrid
+
+    "The missing value in the data respresenting ocean"
+    missing_value::NF = NF(NaN)
+
+    # to be filled from file
+    "High vegetation cover [1], interpolated onto Grid"
+    high_cover::Grid = zeros(Grid,nlat_half)
+
+    "Low vegetation cover [1], interpolated onto Grid"
+    low_cover::Grid = zeros(Grid,nlat_half)
+end
+
+# generator function
+function VegetationClimatology(SG::SpectralGrid;kwargs...)
+    (;NF,Grid,nlat_half) = SG
+    return VegetationClimatology{NF,Grid{NF}}(;nlat_half,kwargs...)
+end
+
+function initialize!(vegetation::VegetationClimatology)
+
+    # LOAD NETCDF FILE
+    if vegetation.path == "SpeedyWeather.jl/input_data"
+        path = joinpath(@__DIR__,"../../input_data",vegetation.file)
+    else
+        path = joinpath(vegetation.path,vegetation.file)
+    end
+    ncfile = NCDataset(path)
+
+    # high and low vegetation cover
+    vegh = vegetation.file_Grid(ncfile[vegetation.varname_vegh][:,:])
+    vegl = vegetation.file_Grid(ncfile[vegetation.varname_vegl][:,:])
+
+    # interpolate onto grid
+    high_vegetation_cover = vegetation.high_cover
+    low_vegetation_cover = vegetation.low_cover
+    interpolator = RingGrids.interpolator(Float32,high_vegetation_cover,vegh)
+    interpolate!(high_vegetation_cover,vegh,interpolator)
+    interpolate!(low_vegetation_cover,vegl,interpolator)
+end
+
+# function barrier
+function soil_moisture_availability!(
+    diagn::SurfaceVariables,
+    land::PrognosticVariablesLand,
+    model::PrimitiveWet)
+    soil_moisture_availability!(diagn,land,model.soil,model.vegetation)
+end
+
+function soil_moisture_availability!(
+    ::SurfaceVariables,
+    ::PrognosticVariablesLand,
+    ::PrimitiveDry)
+    return nothing
+end
+
+function soil_moisture_availability!(
+    diagn::SurfaceVariables{NF},
+    land::PrognosticVariablesLand,
+    soil::AbstractSoil,
+    vegetation::AbstractVegetation) where NF
+
+    (;soil_moisture_availability) = diagn
+    (;soil_moisture_layer1, soil_moisture_layer2) = land
+    (;high_cover, low_cover) = vegetation
+    
+    D_top = convert(NF,soil.D_top)
+    D_root = convert(NF,soil.D_root)
+    W_cap = convert(NF,soil.W_cap)
+    W_wilt = convert(NF,soil.W_wilt)
+    low_veg_factor = convert(NF,vegetation.low_veg_factor)
+
+    # precalculate
+    r = 1/(D_top*W_cap + D_root*(W_cap - W_wilt))
+
+    @inbounds for ij in eachgridpoint(soil_moisture_availability,
+                            soil_moisture_layer1,
+                            soil_moisture_layer2,
+                            high_cover,low_cover)
+        
+        # Fortran SPEEDY source/land_model.f90 line 111 origin unclear
+        veg = max(0,high_cover[ij] + low_veg_factor*low_cover[ij])
+
+        # Fortran SPEEDY documentation eq. 51
+        soil_moisture_availability[ij] = r*(D_top*soil_moisture_layer1[ij] +
+            veg*D_root*max(soil_moisture_layer2[ij]- W_wilt,0))
+    end
+end
diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl
index 3f13a74f7..af67f9033 100644
--- a/src/physics/surface_fluxes.jl
+++ b/src/physics/surface_fluxes.jl
@@ -7,7 +7,6 @@ function surface_fluxes!(column::ColumnVariables,model::PrimitiveEquation)
     surface_wind_stress!(column,model.surface_wind)
 
     # now call other heat and humidity fluxes
-    # soil_moisture!(column,model.vegetation)
     sensible_heat_flux!(column,model.surface_heat_flux,model.constants)
     evaporation!(column,model)
 end
@@ -116,8 +115,8 @@ function sensible_heat_flux!(   column::ColumnVariables{NF},
     flux_sea  = ρ*heat_exchange_sea*V₀*cₚ*(T_skin_sea  - T)
 
     # mix fluxes for fractional land-sea mask
-    land_available = ~isnan(T_skin_land)
-    sea_available = ~isnan(T_skin_sea)
+    land_available = isfinite(T_skin_land)
+    sea_available = isfinite(T_skin_sea)
 
     if land_available && sea_available
         column.flux_temp_upward[end] += land_fraction*flux_land + (1-land_fraction)*flux_sea
@@ -137,59 +136,11 @@ function sensible_heat_flux!(   column::ColumnVariables{NF},
 end
 
 Base.@kwdef struct SurfaceEvaporation{NF<:AbstractFloat} <: AbstractEvaporation{NF}
-
-    nlat_half::Int
-
     moisture_exchange_land::Float64 = 1.2e-3    # for neutral stability
     moisture_exchange_sea::Float64 = 0.9e-3
-
-    "Depth of top soil layer [m]"
-    D_top::Float64 = 0.07
-
-    "Depth of root layer [m]"
-    D_root::Float64 = 0.21
-
-    "Soil wetness at field capacity [volume fraction]"
-    W_cap::Float64 = 0.3
-
-    "Soil wetness at wilting point [volume fraction]"
-    W_wilt::Float64 = 0.17
-
-    # FILE OPTIONS
-    "path to the folder containing the soil moisture and vegetation file, pkg path default"
-    path::String = "SpeedyWeather.jl/input_data"
-
-    "filename of that netcdf file"
-    file::String = "soil_moisture.nc"
-
-    # "Grid the soil and vegetation file comes on"
-    # file_Grid::Type{<:AbstractGrid} = FullGaussianGrid
-
-    # # arrays to be initialised
-    # "Soil moisture availability index"
-    # soil_moisture_availability::Grid = zeros(Grid{NF},nlat_half)
 end
 
-SurfaceEvaporation(SG::SpectralGrid;kwargs...) = SurfaceEvaporation{SG.NF}(nlat_half=SG.nlat_half;kwargs...)
-
-function initialize!(evaporation::SurfaceEvaporation)
-
-    # LOAD NETCDF FILE
-    if evaporation.path == "SpeedyWeather.jl/input_data"
-        path = joinpath(@__DIR__,"../../input_data",evaporation.file)
-    else
-        path = joinpath(evaporation.path,evaporation.file)
-    end
-    ncfile = NCDataset(path)
-
-    # high resolution land-sea mask
-    lsm_highres = evaporation.file_Grid(ncfile["lsm"][:,:])
-
-    # average onto grid cells of the model
-    RingGrids.grid_cell_average!(land_sea_mask.land_sea_mask,lsm_highres)
-end
-
-SurfaceEvaporation(SG::SpectralGrid;kwargs...) = SurfaceEvaporation{SG.NF}(nlat_half=SG.nlat_half;kwargs...)
+SurfaceEvaporation(SG::SpectralGrid;kwargs...) = SurfaceEvaporation{SG.NF}(;kwargs...)
 
 # don't do anything for dry core
 function evaporation!(  column::ColumnVariables,
@@ -210,9 +161,7 @@ function evaporation!(  column::ColumnVariables{NF},
     (;skin_temperature_sea, skin_temperature_land, pres, humid) = column
     moisture_exchange_land = convert(NF,evaporation.moisture_exchange_land)
     moisture_exchange_sea  = convert(NF,evaporation.moisture_exchange_sea)
-
-    # α = soil_moisture.availability[column.ij]
-    α = convert(NF,0.2)
+    α = column.soil_moisture_availability
 
     # SATURATION HUMIDITY OVER LAND AND OCEAN
     surface_pressure = pres[end]
@@ -226,8 +175,8 @@ function evaporation!(  column::ColumnVariables{NF},
     flux_land = ρ*moisture_exchange_land*V₀*α*max(sat_humid_land  - humid[end],0)
 
     # mix fluxes for fractional land-sea mask
-    land_available = ~isnan(skin_temperature_land)
-    sea_available = ~isnan(skin_temperature_sea)
+    land_available = isfinite(skin_temperature_land) && isfinite(α)
+    sea_available = isfinite(skin_temperature_sea)
 
     if land_available && sea_available
         column.flux_humid_upward[end] += land_fraction*flux_land + (1-land_fraction)*flux_sea
diff --git a/src/run_speedy.jl b/src/run_speedy.jl
index b9b359d10..ff5dafcce 100644
--- a/src/run_speedy.jl
+++ b/src/run_speedy.jl
@@ -1,19 +1,14 @@
 """
 $(TYPEDSIGNATURES)
-Run a SpeedyWeather.jl `simulation`. The `simulation.model` is assumed to be initialized,
-otherwise use `initialize=true` as keyword argument."""
+Run a SpeedyWeather.jl `simulation`. The `simulation.model` is assumed to be initialized."""
 function run!(  simulation::Simulation;
                 n_days::Real = 10,
-                startdate::Union{Nothing,DateTime} = nothing,
                 output::Bool = false)
     
     (;prognostic_variables, diagnostic_variables, model) = simulation
     (;clock) = prognostic_variables
 
-    # set the clock
-    if typeof(startdate) == DateTime
-        clock.time = startdate
-    end
+    # set the clock's enddate
     clock.n_days = n_days
     initialize!(clock,model.time_stepping)