Merge pull request #496 from SpeedyWeather/mk/ocean

docs/make.jl

@@ -37,6 +37,7 @@ makedocs(
                 "Parameterizations"=>"parameterizations.md",
                 "Orography"=>"orography.md",
                 "Land-Sea Mask"=>"land_sea_mask.md",
+                "Ocean"=>"ocean.md",
                 "Callbacks"=>"callbacks.md",
             ],
             "RingGrids"=>"ringgrids.md",

docs/src/callbacks.md

@@ -225,7 +225,7 @@ global mean surface temperatures in Kelvin on every time step while the model ra
 model.callbacks[:temperature].temp
 ```
 
-## Intrusive callbacks
+## [Intrusive callbacks](@id intrusive_callbacks)
 
 In the sections above, callbacks were introduced as a tool to define custom
 diagnostics or simulation output. This is the simpler and recommended way of using 
@@ -243,6 +243,10 @@ Another example would be to switch on/off certain model components over time.
 If these components are implemented as *mutable* struct then one could define
 a callback that weakens their respective strength parameter over time.
 
+As an example of a callback that changes the model components see
+
+- Millenium flood: [Time-dependent land-sea mask](@ref)
+
 Changing the diagnostic variables, however, will not have any effect. All of
 them are treated as work arrays, meaning that their state is completely
 overwritten on every time step.  Changing the prognostic variables in spectral space
@@ -315,11 +319,11 @@ function SpeedyWeather.callback!(
     diagn::DiagnosticVariables,
     model::ModelSetup,
 )
-    # scheduled callbacks start with this line to execute only when scheduled!
+    # scheduled callbacks start with this line to execute only when scheduled!
     # else escape immediately
     isscheduled(callback.schedule, progn.clock) || return nothing
 
-    # Just print the North Pole surface temperature to screen
+    # Just print the North Pole surface temperature to screen
     (;time) = progn.clock
     temp_at_north_pole = diagn.layers[end].grid_variables.temp_grid[1]
     @info "North pole has a temperature of $temp_at_north_pole on $time."

docs/src/land_sea_mask.md

@@ -0,0 +1,183 @@
+# Land-sea mask
+
+The following describes how a custom land-sea mask can be defined.
+SpeedyWeather uses a _fractional_ land-sea mask, i.e. for every grid-point
+
+- 1 indicates land
+- 0 indicates ocean
+- a value in between indicates a grid-cell partially covered by ocean and land
+
+The land-sea mask determines solely how to weight the surface fluxes
+coming from land or from the ocean. For the sensible heat fluxes this uses
+land and sea surface temperatures and weights the respective fluxes
+proportional to the fractional mask. Similar for evaporation.
+You can therefore define an ocean on top of a mountain, or a land without
+heat fluxes when the land-surface temperature is not defined, i.e. `NaN`.
+Let ``F_L, F_S`` be the fluxes coming from land and sea, respectively.
+Then the land-sea mask ``a \in [0,1]`` weights them as follows for the
+total flux ``F``
+
+```math
+F = aF_L + (1-a)F_S
+```
+
+but ``F=F_L`` if the sea flux is NaN (because the ocean temperature is not defined)
+and ``F=F_S`` if the land flux is NaN (because the land temperature or soil moisture
+is not defined, for sensible heat fluxes or evaporation), and ``F=0`` if both fluxes
+are NaN.
+
+Setting the land-sea mask to ocean therefore will disable any fluxes that
+may come from land, and vice versa. However, with an ocean-everywhere land-sea mask
+you must also define sea surface temperatures everywhere, otherwise the fluxes
+in those regions will be zero.
+
+## Manual land-sea mask
+
+You can create the default land-sea mask as follows
+
+```@example landseamask
+using SpeedyWeather
+spectral_grid = SpectralGrid(trunc=31, nlev=8)
+land_sea_mask = LandSeaMask(spectral_grid)
+```
+
+which will automatically interpolate the land-sea mask onto grid and resolution
+as defined in `spectral_grid` at initialization. The actual mask is in
+`land_sea_mask.land_sea_mask` and you can visualise it with
+
+```@example landseamask
+model = PrimitiveWetModel(;spectral_grid, land_sea_mask)
+simulation = initialize!(model)     # triggers also initialization of model.land_sea_mask
+plot(land_sea_mask.mask)
+```
+
+Now before you run a simulation you could manually change the land-sea mask by
+
+```@example landseamask
+# unpack, this is a flat copy, changing it will also change the mask inside model
+(; mask) = land_sea_mask
+
+# ocean everywhere, or
+mask .= 0    
+
+# random land-sea mask, or
+for i in eachindex(mask)
+    mask[i] = rand()     
+end
+
+# ocean only between 10˚S and 10˚N
+for (j, ring) in enumerate(RingGrids.eachring(mask))
+    for ij in ring
+        mask[ij] = abs(model.geometry.latd[j]) > 10 ? 1 : 0
+    end
+end
+```
+
+And now you can run the simulation as usual with `run!(simulation)`. Most useful
+for the generation of custom land-sea masks in this manual way is probably the
+`model.geometry` component which has all sorts of coordinates like `latd`
+(latitudes in degrees on rings) or `latds, londs` (latitude, longitude in degrees
+for every grid point).
+
+## Earth's land-sea mask
+
+The Earth's [`LandSeaMask`](@ref) has itself the option to load another
+land-sea mask from file, but you also have to specify the grid that mask
+from files comes on. It will then attempt to read it via `NCDatasets`
+and interpolate onto the model grid.
+
+## AquaPlanetMask
+
+Predefined is also the [`AquaPlanetMask`](@ref) which can be created as
+```@example landseamask
+land_sea_mask = AquaPlanetMask(spectral_grid)
+```
+and is equivalent to using Earth's [`LandSeaMask`](@ref) but setting
+the entire mask to zero afterwards `land_sea_mask.mask .= 0`.
+
+## Custom land-sea mask
+
+Every (custom) land-sea mask has to be a subtype of `AbstractLandSeaMask`.
+A custom land-sea mask has to be defined as a new type (`struct` or `mutable struct`)
+
+```julia
+CustomMask{NF<:AbstractFloat, Grid<:AbstractGrid{NF}} <: AbstractLandSeaMask{NF, Grid}
+```
+
+and needs to have at least a field called `mask::Grid` that uses a `Grid` as defined
+by the spectral grid object, so of correct size and with the number format `NF`.
+All `AbstractLandSeaMask` have a convenient generator function to be used like
+`mask = CustomMask(spectral_grid, option=argument)`, but you may add your own or customize by
+defining `CustomMask(args...)` which should return an instance of type `CustomMask{NF, Grid}`
+with parameters matching the spectral grid. Then the initialize function has to be extended for
+that new mask
+
+```julia
+initialize!(mask::CustomMask, model::PrimitiveEquation)
+```
+
+which generally is used to tweak the mask.mask grid as you like, using
+any other options you have included in `CustomMask` as fields or anything else (preferably read-only,
+because this is only to initialize the land-sea mask, nothing else) from `model`. You can
+for example read something from file, set some values manually, or use coordinates from `model.geometry`.
+
+## Time-dependent land-sea mask
+
+It is possible to define an [intrusive callback](@ref intrusive_callbacks) to change the
+land-sea mask during integration. The grid in `model.land_sea_mask.mask`
+is mutable, meaning you can change the values of grid points in-place but not replace
+the entire mask or change its size. If that mask is changed, this will be reflected
+in all relevant model components. For example, we can define a callback that
+floods the entire planet at the beginning of the 21st century as
+
+```@example landseamask
+Base.@kwdef struct MilleniumFlood <: SpeedyWeather.AbstractCallback
+    schedule::Schedule = Schedule(DateTime(2000,1,1))
+end
+
+# initialize the schedule
+function SpeedyWeather.initialize!(
+    callback::MilleniumFlood,
+    progn::PrognosticVariables,
+    diagn::DiagnosticVariables,
+    model::ModelSetup,
+)
+    initialize!(callback.schedule, progn.clock)
+end
+
+function SpeedyWeather.callback!(
+    callback::MilleniumFlood,
+    progn::PrognosticVariables,
+    diagn::DiagnosticVariables,
+    model::ModelSetup,
+)
+    # escape immediately if not scheduled yet
+    isscheduled(callback.schedule, progn.clock) || return nothing
+
+    # otherwise set the entire land-sea mask to ocean
+    model.land_sea_mask.mask .= 0
+    @info "Everything flooded on $(progn.clock.time)"
+end
+
+# nothing needs to be done when finishing
+SpeedyWeather.finish!(::MilleniumFlood, args...) = nothing
+```
+
+Note that the flooding will take place only at the start of the 21st century,
+last indefinitely, but not if the model integration period does not cover that
+exact event, see [Schedules](@ref). Initializing a model a few days earlier
+would then have this `MilleniumFlood` take place
+
+```@example landseamask
+land_sea_mask = LandSeaMask(spectral_grid)      # start with Earth's land-sea mask
+model = PrimitiveWetModel(;spectral_grid, land_sea_mask)
+add!(model, MilleniumFlood())   # or MilleniumFlood(::DateTime) for any non-default date
+model.feedback.verbose = false # hide
+
+simulation = initialize!(model, time=DateTime(1999,12,29))
+run!(simulation, period=Day(5))
+plot(model.land_sea_mask.mask)
+```
+
+And the land-sea mask has succesfully been set to ocean everywhere at the start
+of the 21st century.

docs/src/ocean.md

@@ -0,0 +1,2 @@
+# Ocean
+

docs/src/parameterizations.md

@@ -1,3 +1 @@
-# Parameterizations
-
-More to follow ...
+# Parameterizations

docs/src/radiation.md

@@ -0,0 +1 @@
+# Radiation

docs/src/setups.md

@@ -36,8 +36,8 @@ There are other options to create a planet but they are irrelevant for the
 barotropic vorticity equations. We also want to specify the initial conditions,
 randomly distributed vorticity is already defined
 ```@example barotropic_setup
-using Random # hide
-Random.seed!(1234) # hide
+using Random # hide
+Random.seed!(1234) # hide
 initial_conditions = StartWithRandomVorticity()
 ```
 By default, the power of vorticity is spectrally distributed with ``k^{-3}``, ``k`` being the
@@ -124,7 +124,7 @@ Vorticity `vor` is stored as a lon x lat x vert x time array, we may want to loo
 which is the end of the previous simulation (time=6days) which we didn't store output for.
 ```@example galewsky_setup
 t = 1
-vor = Matrix{Float32}(ds["vor"][:, :, 1, t]) # convert from Matrix{Union{Missing, Float32}} to Matrix{Float32}
+vor = Matrix{Float32}(ds["vor"][:, :, 1, t]) # convert from Matrix{Union{Missing, Float32}} to Matrix{Float32}
 lat = ds["lat"][:]
 lon = ds["lon"][:]
 
@@ -211,7 +211,7 @@ simulation = initialize!(model)
 model.feedback.verbose = false # hide
 run!(simulation, period=Day(20))   # discard first 20 days   
 run!(simulation, period=Day(20), output=true)
-nothing # hide
+nothing # hide
 ```
 
 We want to simulate polar jet streams in the shallow water model. We add a `JetStreamForcing`
@@ -250,8 +250,8 @@ nothing # hide
 
 Setup script to copy and paste:
 ```@example gravity_wave_setup
-using Random # hide
-Random.seed!(1234) # hide
+using Random # hide
+Random.seed!(1234) # hide
 using SpeedyWeather
 spectral_grid = SpectralGrid(trunc=127, nlev=1)
 time_stepping = SpeedyWeather.Leapfrog(spectral_grid, Δt_at_T31=Minute(30))
@@ -361,6 +361,70 @@ nothing # hide
 ```
 ![Jablonowski pyplot](jablonowski.png)
 
+## Aquaplanet
+
+```@example aquaplanet
+using SpeedyWeather
+
+# components
+spectral_grid = SpectralGrid(trunc=31, nlev=5)
+ocean = AquaPlanet(spectral_grid, temp_equator=302, temp_poles=273)
+land_sea_mask = AquaPlanetMask(spectral_grid)
+orography = NoOrography(spectral_grid)
+
+# create model, initialize, run
+model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography)
+simulation = initialize!(model)
+model.feedback.verbose = false # hide
+run!(simulation, period=Day(50), output=true)
+nothing # hide
+```
+
+Here we have defined an aquaplanet simulation by
+- creating an `ocean::AquaPlanet`. This will use constant sea surface temperatures that only vary with latitude.
+- creating a `land_sea_mask::AquaPlanetMask` this will use a land-sea mask with `false`=ocean everywhere.
+- creating an `orography::NoOrography` which will have no orography and zero surface geopotential.
+
+All passed on to the model constructor for a `PrimitiveWetModel`, we have now a model with humidity
+and physics parameterization as they are defined by default (typing `model` will give you an overview
+of its components). We could have change the `model.land` and `model.vegetation` components too,
+but given the land-sea masks masks those contributions to the surface fluxes anyway, this is not
+necessary. Note that neither sea surface temperature, land-sea mask
+or orography have to agree. It is possible to have an ocean on top of a mountain.
+For an ocean grid-cell that is (partially) masked by the land-sea mask, its value will
+be (fractionally) ignored in the calculation of surface fluxes (potentially leading
+to a zero flux depending on land surface temperatures).
+
+Now with the following we visualize the surface humidity after the 50 days of
+simulation. We use 50 days as without mountains it takes longer for the initial conditions to
+become unstable. The surface humidity shows small-scale patches in the tropics, which is a result
+of the convection scheme, causing updrafts and downdrafts in both humidity and temperature.
+
+```@example aquaplanet
+using PythonPlot, NCDatasets
+ioff() # hide
+
+id = model.output.id
+ds = NCDataset("run_$id/output.nc")
+timestep = ds.dim["time"]   # last time step
+surface = ds.dim["lev"]     # surface layer
+humid = Matrix{Float32}(ds["humid"][:, :, surface, timestep])
+lat = ds["lat"][:]
+lon = ds["lon"][:]
+
+fig, ax = subplots(1, 1, figsize=(10, 6))
+q = ax.pcolormesh(lon, lat, humid')
+ax.set_xlabel("longitude")
+ax.set_ylabel("latitude")
+ax.set_title("Surface humidity [kg/kg]")
+colorbar(q)
+tight_layout() # hide
+savefig("aquaplanet.png", dpi=70) # hide
+nothing # hide
+```
+![Aquaplanet pyplot](aquaplanet.png)
+
+
 
 ## References
 

docs/src/surface_fluxes.md

@@ -1,3 +1,3 @@
-# Surface fluxes
+# Surface fluxes
 
 More to follow ... 

src/physics/column_variables.jl

@@ -32,16 +32,16 @@ function get_column!(   C::ColumnVariables,
     C.lond = geometry.londs[ij]
     C.ij = ij                       # grid-point index
     C.jring = jring                 # ring index j of column, used to index latitude vectors
-    C.land_fraction = land_sea_mask.land_sea_mask[ij]
+    C.land_fraction = land_sea_mask.mask[ij]
     C.orography = orography.orography[ij]
     C.surface_geopotential = C.orography * planet.gravity
 
     # pressure [Pa] or [log(Pa)]
-    lnpₛ = D.surface.pres_grid[ij]          # logarithm of surf pressure used in dynamics
+    lnpₛ = D.surface.pres_grid[ij]          # logarithm of surf pressure used in dynamics
     pₛ = exp(lnpₛ)                          # convert back here
     C.ln_pres .= ln_σ_levels_full .+ lnpₛ   # log pressure on every level ln(p) = ln(σ) + ln(pₛ)
-    C.pres[1:end-1] .= σ_levels_full.*pₛ    # pressure on every level p = σ*pₛ
-    C.pres[end] = pₛ                        # last element is surface pressure pₛ
+    C.pres[1:end-1] .= σ_levels_full.*pₛ    # pressure on every level p = σ*pₛ
+    C.pres[end] = pₛ                        # last element is surface pressure pₛ
 
     @inbounds for (k, layer) = enumerate(D.layers)   # read out prognostic variables on grid
         C.u[k] = layer.grid_variables.u_grid[ij]