Update docs and source

README.md

@@ -1,18 +1,36 @@
 # EnergyBalanceModels.jl
 
-Any quantity that may potentially be a dynamic (state) variable, or an explicit function of time, is defined in `src/variables.jl` and must have a default value when defined. They should also be exported and in the rest of the source code they are accessed from the module-level scope. _TODO: Example here_.
-Variables that may be dynamic (state) variables must also obtain limits in the `physicall_plausible_limits` function which is in the same file.
-Variables that can never be dynamic but are guaranteed to be observed, such as the outgoing longwave radiation, should be preferably defined in the same file but without a default value or limits.
+[![docsdev](https://img.shields.io/badge/docs-dev-lightblue.svg)](https://juliadynamics.github.io/ProcessBasedModelling,jl/dev/)
+[![docsstable](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliadynamics.github.io/EnergyBalanceModels.jl/stable/)
+[![CI](https://github.com/JuliaDynamics/EnergyBalanceModels.jl/workflows/CI/badge.svg)](https://github.com/JuliaDynamics/EnergyBalanceModels.jl/actions?query=workflow%3ACI)
+[![codecov](https://codecov.io/gh/JuliaDynamics/EnergyBalanceModels.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaDynamics/EnergyBalanceModels.jl)
+[![Package Downloads](https://shields.io/endpoint?url=https://pkgs.genieframework.com/api/v1/badge/ProcessBasedModelling)](https://pkgs.genieframework.com?packages=ProcessBasedModelling)
 
-# TODO:
-explain that all processes by default use the global variables,
-which is why docstrings are written this way.
+EnergyBalanceModels.jl is a Julia package for creating and analyzing conceptual
+models of climate, such as energy balance models or climate tipping models.
+Such conceptual models are simplified representation of basic climate components,
+and the processes that connect them, such as flows of energy or mass.
+Within this context such models are typically coupled ordinary differential
+equations (with partial or stochastic DEs also being possible).
 
-EnergyBalanceModels
+EnergyBalanceModels.jl accelerates both modelling and analysis aspects of working
+with such models by:
 
-"""
-module defining various equations and models representing energy balance
-processes. It utilizes ModelingToolkit.jl to create a flexible framework
-for adding or remove dynamical state variables and adding or removing
-equations that represent various processes in the earth system.
-"""
+- Building upon [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/)
+  for creating equations from _symbolic expressions_.
+- Utilizing [ProcessBasedModelling.jl](https://github.com/JuliaDynamics/ProcessBasedModelling.jl?tab=readme-ov-file)
+  to provide a field-specific framework that allows easily testing different physical
+  hypotheses regarding how climate variables couple to each
+  other, or how climate processes are represented by equations.
+- Offering many predefined processes from current literature and ongoing research.
+- Being easy to extend with more climate variables or physical processes.
+- Allowing the straightforward coupling of different conceptual models with each other.
+- Automating the naming of custom parameters relating to existing climate processes.
+- Integrating with [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/)
+  to automate the start-up phase of typical nonlinear dynamics based workflows.
+
+To install it, run `import Pkg; Pkg.add("EnergyBalanceModels")`.
+
+All further information is provided in the documentation, which you can either find
+[online](https://juliadynamics.github.io/EnergyBalanceModels.jl/stable/) or build
+locally by running the `docs/make.jl` file.

docs/make.jl

@@ -11,6 +11,8 @@ include("build_docs_with_style.jl")
 
 pages =  [
     "Introduction" => "index.md",
+    "Predefined processes" => "processes.md",
+    "References" => "references.md",
 ]
 
 using DocumenterCitations

docs/src/index.md

@@ -1,69 +1,70 @@
 # EnergyBalanceModels.jl
 
-# ## Variables and default values
 
-...
 
-# ## Temperature
+## [Premade symbolic variable instances](@id global_vars)
 
-```@docs
-BasicRadiationBalance
-ΔTLinearRelaxation
-```
+For convenience, EnergyBalanceModels.jl defines and exports some symbolic variables
+that we [list below](@ref list_vars). These are used throughout the library as the
+default variables in [implemented processes](@id processes).
+When going through documentation strings of processes, such as [`BasicRadiationBalance`](@ref),
+you will notice that the function call signature is like:
 
-# ## Longwave
+```julia
+BasicRadiationBalance(; T=T, kwargs...)
+```
 
-```@docs
-LinearOLR
-LinearClearSkyOLR
-EmissivityStefanBoltzmanOLR
-EmissivityFeedbackTanh
-EmissivitySellers1969
-SoedergrenClearSkyEmissivity
+This `T=T` means that the keyword argument `T`, which represents the
+"temperature variable", takes the value of `EnergyBalanceModels.T`,
+which itself is a premade symbolic variable instance that is exported by
+EnergyBalanceModels.jl. You can pass in your own variables instead, by doing
+```julia
+@variables begin
+    (T1_tropics(t) = 290.0), [bounds = (200.0, 350.0), description = "temperature in tropical box 1, in Kelvin"]
+end
+BasicRadiationBalance(; T=T1_tropics, kwargs...)
 ```
+_(you would also need to give `T1_tropics` to all other processes that utilize temperature)_
 
-# ## Shortwave
+Defining variables with the extra `bounds, description` annotations is
+useful for integrating with the rest of the functionality of the library.
 
-```@docs
-DirectAlbedoAddition
-CoAlbedoProduct
-SeparatedClearAllSkyAlbedo
+### [List of premade variables](@id list_vars)
+
+```@example MAIN
+using EnergyBalanceModels
+PREDEFINED_EBM_VARIABLES
 ```
 
-# ## Ice/snow
+## Default values, limits, etc.
+
+All [exported symbolic variable instances](@ref) are defiled with a default value and have plausible physical limits that can be obtained with the following function:
 
 ```@docs
-IceAlbedoFeedback
+physically_plausible_limits(::Any)
 ```
 
-# ## Insolation
+e.g.,
 
-```@docs
-AstronomicalForcingDeSaedeleer
+```@example MAIN
+physically_plausible_limits(T)
 ```
 
-# ## Forcings
 
-```@docs
-CO2Forcing
-```
+## Integration with DynamicalSystems.jl
+
+EnergyBalanceModels.jl integrates with [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/) in many ways.
 
-# ## Clouds
 
 ```@docs
-CloudAlbedoExponential
-CloudAlbedoLinear
-EmissivityCumulativeSubtraction
-BudykoOLR
+physically_plausible_limits(::DynamicalSystem)
+physically_plausible_ic_sampler
+physically_plausible_grid
 ```
 
-# ## Utils
+# ## Utilities
 
 ```@docs
 TanhProcess
 ```
 
-# ## References
-
-```@bibliography
-```

docs/src/processes.md

@@ -0,0 +1,57 @@
+# Predefined processes
+
+Predefined processes are provided in this page.
+Those extracted by the literature cite their according resource.
+Note that by default all processes utilize the globally-exported [predefined variables](@ref global_vars) of EnergyBalanceModels.jl.
+
+## Temperature
+
+```@docs
+BasicRadiationBalance
+ΔTLinearRelaxation
+```
+
+## Longwave radiation
+
+```@docs
+LinearOLR
+LinearClearSkyOLR
+EmissivityStefanBoltzmanOLR
+EmissivityFeedbackTanh
+EmissivitySellers1969
+SoedergrenClearSkyEmissivity
+```
+
+## Shortwave radiation
+
+```@docs
+DirectAlbedoAddition
+CoAlbedoProduct
+SeparatedClearAllSkyAlbedo
+```
+
+## Ice/snow
+
+```@docs
+IceAlbedoFeedback
+```
+
+## Insolation
+
+```@docs
+AstronomicalForcingDeSaedeleer
+```
+
+## Forcings
+
+```@docs
+CO2Forcing
+```
+
+## Clouds
+
+```@docs
+CloudAlbedoExponential
+CloudAlbedoLinear
+BudykoOLR
+```

docs/src/references.md

@@ -0,0 +1,4 @@
+# References
+
+```@bibliography
+```

src/processes/clouds/albedo.jl

@@ -0,0 +1,37 @@
+export CloudAlbedoExponential, CloudAlbedoLinear
+
+"""
+    CloudAlbedoExponential(
+        α_cloud = α_cloud, C = C, a = 2.499219232848238, b = 17.596369331717433
+    )
+
+Create the equation `α_cloud ~ sinh(a*C)/b` relating cloud albedo
+to cloud fraction `C`. This equation is exponential and not linear, as in observations.
+[Engstrom2015](@cite) (and also [Bender2017](@cite)) discuss this exponential relation
+in detail, and provide as explanation that cloud effective albedo increases with latitude
+(due to solar zenith changes) while cloud fraction also increases with latitude.
+
+Note that here however we modify the equation `α_cloud ~ exp(a*C - b)` of [Engstrom2015](@cite)
+to utilize the hyperbolic sine, so that `α_cloud = 0` when `C = 0` as is physically
+necessary. Then, `a, b` are extracted by fitting CERES data, using as `α_cloud` the
+energetically consistent cloud albedo as defined by [Datseris2021](@cite),
+further yearly averaged and within latitudes (-60, 60) as in [Bender2017](@cite).
+This albedo can be directly added to the clear sky albedo to produce the planetary albedo.
+"""
+function CloudAlbedoExponential(;
+        α_cloud = α_cloud, C = C, a = 2.499219232848238, b = 17.596369331717433
+    )
+    return α_cloud ~ sinh(a*C)/b
+end
+
+"""
+    CloudAlbedoLinear(; α_cloud = α_cloud, C = C, a = 0.2252861764703435)
+
+Same as in [`CloudAlbedoExponential`](@ref), but now the linear form `α_cloud ~ a*C`
+is returned, with `a` fitted from CERES data in the same way.
+"""
+function CloudAlbedoLinear(;
+        α_cloud = α_cloud, C = C, a = 0.2252861764703435
+    )
+    return α_cloud ~ a*C
+end

src/processes/clouds/longwave.jl

@@ -0,0 +1,24 @@
+export BudykoOLR
+
+"""
+    BudykoOLR(; T=T, C=C,
+        A = -461.8068, B = 2.58978, Ac = -377.22741, Bc = 1.536171
+    )
+
+Create the equation `OLR ~ A + B*T - C*(Ac + Bc*T)` for the dependence of OLR on
+both temperature and cloud fraction (in 0-1). This is the same
+as Eq. (1) of [Budyko1969](@cite). However, here `T` is expected in Kelvin,
+and the coefficients have been extracted by fitting into CERES data in the same
+way as in [`LinearOLR`](@ref).
+"""
+function BudykoOLR(; T=T, C=C,
+        A = -461.8068, B = 2.58978, Ac = -377.22741, Bc = 1.536171
+    )
+    @parameters begin
+        Budyko_A = A
+        Budyko_B = B
+        Budyko_Ac = Ac
+        Budyko_Bc = Bc
+    end
+    return OLR ~ Budyko_A + Budyko_B*T - C*(Budyko_Ac + Budyko_Bc*T)
+end

src/processes_advanced.jl

@@ -12,7 +12,7 @@ The rest of the input arguments should be real numbers or `@parameter` named par
 
 The process creates the expression:
 ```
-variable ~ left + (right - mvap)*(1 + tanh(2(driver - reference)/scale))*0.5
+variable ~ left + (right - left)*(1 + tanh(2(driver - reference)/scale))*0.5
 ```
 i.e., a tanh formula that goes from value `left` to value `right` as a function of `driver`
 over a range of `scale` being centered at `reference`.

src/statespace.jl

@@ -6,14 +6,20 @@ The states need to be named as the limits are deduced from the function
 `physically_plausible_limits(varname::String)`.
 """
 function physically_plausible_limits(ds::DynamicalSystem)
-    sys = referrenced_sciml_model(ds)
-    vars = ModelingToolkit.states(sys)
-    varstrings = @. string(ModelingToolkit.getname(vars))
-    minmax = physically_plausible_limits.(varstrings)
+    model = referrenced_sciml_model(ds)
+    vars = ModelingToolkit.states(model)
+    minmax = physically_plausible_limits.(vars)
     return minmax
 end
 
+"""
+    physically_plausible_ic_sampler(ds::DynamicalSystem)
 
+Return a `sampler` that can be called as a 0-argument function `sampler()`, which
+yields random initial conditions within the hyperrectangle defined by the
+[`physically_plausible_limits`](@ref) of `ds`.
+The `sampler` is useful to give to e.g., `DynamicalSystems.basins_fractions`.
+"""
 function physically_plausible_ic_sampler(ds::DynamicalSystem)
     minmax = physically_plausible_limits(ds)
     mini = getindex.(minmax, 1)
@@ -22,6 +28,14 @@ function physically_plausible_ic_sampler(ds::DynamicalSystem)
     return sampler
 end
 
+"""
+    physically_plausible_grid(ds::DynamicalSystem, n = 101)
+
+Return a `grid` that is a tuple of `range` objects that each spans the
+[`physically_plausible_limits`](@ref) of `ds`.
+`n` can be an integer or a vector of integers (for each dimension).
+The resulting `grid` is useful to give to `DynamicalSystems.AttractorsViaRecurrences`.
+"""
 function physically_plausible_grid(ds::DynamicalSystem, n = 101)
     minmax = physically_plausible_limits(ds)
     if n isa Integer