diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
index de7ac75..9592764 100644
--- a/.github/ISSUE_TEMPLATE/bug_report.md
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -1,27 +1,30 @@
 ---
-name: Bug report
-about: Create a report to help us improve
+name: Problem/bug report
+about: Create a report for incorrect/unexpected/erroring behavior
 title: ''
 labels: ''
 assignees: ''
 
 ---
 
-**Describe the bug**
-A clear and concise description of what the bug is.
+**Describe the problem**
+
+A clear and concise description of what the problem is, and what you were expecting to happen instead.
 
 **Minimal Working Example**
+
 Please provide a piece of code that leads to the bug you encounter.
 
-If the code is **runnable**, it will help us identify the problem faster.
+This code should be **runnable**, **self-contained**, and minimal. Remove all irrelevant dependencies
+and actively try to reduce the code to its smallest possible version.
+This will help us identify and fix the problem faster.
 
 **Package versions**
 
 Please provide the versions you use. To do this, run the code:
 ```julia
-using Pkg
-Pkg.status([
-    "Package1", "Package2"]; # etc.
+import Pkg
+Pkg.status(["Package1", "Package2"]; # etc.
     mode = PKGMODE_MANIFEST
 )
 ```
diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md
index 91bb1de..c39ff6d 100644
--- a/.github/ISSUE_TEMPLATE/feature_request.md
+++ b/.github/ISSUE_TEMPLATE/feature_request.md
@@ -1,6 +1,6 @@
 ---
 name: Feature request
-about: Suggest an idea for this project
+about: Suggest a new feature for this project
 title: ''
 labels: ''
 assignees: ''
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
index 6c2390a..7a9c79f 100644
--- a/.github/workflows/CompatHelper.yml
+++ b/.github/workflows/CompatHelper.yml
@@ -1,24 +1,44 @@
 name: CompatHelper
-
 on:
   schedule:
-    - cron: '00 * * * *'
-
+    - cron: 0 0 * * *
+  workflow_dispatch:
+permissions:
+  contents: write
+  pull-requests: write
 jobs:
   CompatHelper:
-    runs-on: ${{ matrix.os }}
-    strategy:
-      matrix:
-        julia-version: [1]
-        julia-arch: [x86]
-        os: [ubuntu-latest]
+    runs-on: ubuntu-latest
     steps:
-      - uses: julia-actions/setup-julia@latest
+      - name: Check if Julia is already available in the PATH
+        id: julia_in_path
+        run: which julia
+        continue-on-error: true
+      - name: Install Julia, but only if it is not already available in the PATH
+        uses: julia-actions/setup-julia@v1
         with:
-          version: ${{ matrix.julia-version }}
-      - name: Pkg.add("CompatHelper")
-        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
-      - name: CompatHelper.main()
+          version: '1'
+          arch: ${{ runner.arch }}
+        if: steps.julia_in_path.outcome != 'success'
+      - name: "Add the General registry via Git"
+        run: |
+          import Pkg
+          ENV["JULIA_PKG_SERVER"] = ""
+          Pkg.Registry.add("General")
+        shell: julia --color=yes {0}
+      - name: "Install CompatHelper"
+        run: |
+          import Pkg
+          name = "CompatHelper"
+          uuid = "aa819f21-2bde-4658-8897-bab36330d9b7"
+          version = "3"
+          Pkg.add(; name, uuid, version)
+        shell: julia --color=yes {0}
+      - name: "Run CompatHelper"
+        run: |
+          import CompatHelper
+          CompatHelper.main()
+        shell: julia --color=yes {0}
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-        run: julia -e 'using CompatHelper; CompatHelper.main()'
+          COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
index 623860f..f49313b 100644
--- a/.github/workflows/TagBot.yml
+++ b/.github/workflows/TagBot.yml
@@ -12,4 +12,4 @@ jobs:
       - uses: JuliaRegistries/TagBot@v1
         with:
           token: ${{ secrets.GITHUB_TOKEN }}
-          ssh: ${{ secrets.DOCUMENTER_KEY }}
\ No newline at end of file
+          ssh: ${{ secrets.DOCUMENTER_KEY }}
diff --git a/.github/workflows/doccleanup.yml b/.github/workflows/doccleanup.yml
index 22d9043..bc29462 100644
--- a/.github/workflows/doccleanup.yml
+++ b/.github/workflows/doccleanup.yml
@@ -23,4 +23,4 @@ jobs:
               git push --force origin gh-pages-new:gh-pages
             fi
         env:
-            PRNUM: ${{ github.event.number }}
\ No newline at end of file
+            PRNUM: ${{ github.event.number }}
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index 33e080c..724cecd 100644
--- a/.github/workflows/documentation.yml
+++ b/.github/workflows/documentation.yml
@@ -22,5 +22,5 @@ jobs:
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
-        run: julia --project=docs/ docs/make.jl
\ No newline at end of file
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key for using TagBot
+        run: julia --project=docs/ docs/make.jl
diff --git a/.github/workflows/ci.yml b/.github/workflows/test_coverage.yml
similarity index 98%
rename from .github/workflows/ci.yml
rename to .github/workflows/test_coverage.yml
index 6278b77..5e9a78b 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/test_coverage.yml
@@ -1,4 +1,4 @@
-name: CI
+name: CI - test & coverage
 on:
   pull_request:
     branches:
diff --git a/Project.toml b/Project.toml
index 6c306a0..43c2b71 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ConceptualClimateModels"
 uuid = "3e640dcb-3446-4ed7-aadb-0366b77312ec"
 authors = ["Datseris <datseris.george@gmail.com>"]
-version = "0.1.3"
+version = "0.2.0"
 
 [deps]
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
diff --git a/docs/make.jl b/docs/make.jl
index af3128b..a73caf7 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,6 +1,7 @@
 cd(@__DIR__)
 
 using ConceptualClimateModels
+using ConceptualClimateModels.GlobalMeanEBM
 
 import Downloads
 Downloads.download(
@@ -12,7 +13,9 @@ include("build_docs_with_style.jl")
 pages =  [
     "Introduction" => "index.md",
     "Tutorial" => "tutorial.md",
-    "Predefined processes" => "processes.md",
+    "Submodules" => [
+        "Global mean EBM" => "submodules/globalmeanebm.md",
+    ],
     "Examples" => "examples.md",
     "References" => "references.md",
 ]
@@ -24,7 +27,7 @@ bib = CitationBibliography(
     style=:authoryear
 )
 
-build_docs_with_style(pages, ConceptualClimateModels, ProcessBasedModelling;
+build_docs_with_style(pages, ConceptualClimateModels, GlobalMeanEBM, ProcessBasedModelling;
     authors = "George Datseris <datseris.george@gmail.com>",
     bib, warnonly = [:doctest, :missing_docs, :cross_references],
     repo = Remotes.GitHub("JuliaDynamics", "ConceptualClimateModels.jl")
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 5ba1014..c11f31e 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -18,7 +18,7 @@ We will combine the processes:
 
 ```@example MAIN
 using ConceptualClimateModels
-using ConceptualClimateModels.CCMV
+using ConceptualClimateModels.GlobalMeanEBM
 
 budyko_processes = [
     BasicRadiationBalance(),
@@ -30,7 +30,7 @@ budyko_processes = [
     # absorbed solar radiation has a default process
 ]
 
-budyko = processes_to_coupledodes(budyko_processes)
+budyko = processes_to_coupledodes(budyko_processes, GlobalMeanEBM)
 println(dynamical_system_summary(budyko))
 ```
 
@@ -44,10 +44,10 @@ For setting up the continuation we leverage the integration with DynamicalSystem
 ```@example MAIN
 using DynamicalSystems
 
-grid = physically_plausible_grid(budyko)
+grid = plausible_grid(budyko)
 mapper = AttractorsViaRecurrences(budyko, grid)
 rfam = RecurrencesFindAndMatch(mapper)
-sampler = physically_plausible_ic_sampler(budyko)
+sampler = plausible_ic_sampler(budyko)
 sampler() # randomly sample initial conditions
 ```
 
@@ -61,7 +61,7 @@ index = :ε_0
 Now we perform the continuation versus the effective emissivity, to approximate increasing or decreasing the strength of the greenhouse effect:
 ```@example MAIN
 εrange = 0.3:0.01:0.8
-fractions_curves, attractors_info = continuation(
+fractions_curves, attractors_info = global_continuation(
     rfam, εrange, index, sampler;
     samples_per_parameter = 1000,
     show_progress = false,
@@ -85,7 +85,7 @@ dependent variable instead, like so:
 
 ```@example MAIN
 using ConceptualClimateModels
-using ConceptualClimateModels.CCMV
+using ConceptualClimateModels.GlobalMeanEBM
 
 @variables η1(t)
 @parameters η1_0 = 2.0 # starting value for η1 parameter
@@ -96,7 +96,7 @@ processes = [
     η1 ~ η1_0 + r_η*t, # this symbolic variable has its own equation!
 ]
 
-stommel = processes_to_coupledodes(processes)
+stommel = processes_to_coupledodes(processes, GlobalMeanEBM)
 println(dynamical_system_summary(stommel))
 ```
 
@@ -106,15 +106,15 @@ estimate its bifurcation diagram using the same steps as the above example.
 ```@example MAIN
 using DynamicalSystems
 
-grid = physically_plausible_grid(stommel)
+grid = plausible_grid(stommel)
 mapper = AttractorsViaRecurrences(stommel, grid;
     consecutive_recurrences = 1000, attractor_locate_steps = 10,
 )
 rfam = RecurrencesFindAndMatch(mapper)
-sampler = physically_plausible_ic_sampler(stommel)
+sampler = plausible_ic_sampler(stommel)
 
 ηrange = 2.0:0.01:4.0
-fractions_curves, attractors_info = continuation(
+fractions_curves, attractors_info = global_continuation(
     rfam, ηrange, η1_0, sampler;
     samples_per_parameter = 1000,
     show_progress = false,
@@ -159,14 +159,3 @@ fig
 As you can see from the figure, depending on the rate the system either "tracks"
 the fixed point of high ΔΤ or it collapses down to the small ΔT branch.
 This happens because the system crosses the unstable manifold of the lower branch [Datseris2022; Chap. 12](@cite).
-To visualize the unstable manifold we could use BifurcationKit.jl, however,
-it is very inconvenient to do so, because BifurcationKit.jl does not provide most of the conveniences
-that DynamicalSystems.jl does. For example, it does not integrate well enough with
-DifferentialEquations.jl (to allow passing `ODEProblem` which is created by `DynamicalSystem`).
-It also does not allow indexing parameters by their symbolic bindings.
-Lastly, it does not work with models generated with ModelingToolkit.jl so we would
-have to re-write all the equations that the chosen `processes` made for us.
-
-## Glacial oscillations
-
-Coming soon!
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index 92e0769..46f42cd 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -5,29 +5,14 @@ ConceptualClimateModels
 ```
 
 To get started with ConceptualClimateModels.jl see the [tutorial](@ref tutorial).
-Predefined processes that can be part of a model are in the
-[predefined processes](@ref predefined_processes) page.
+Predefined processes that can be part of a model are in individual module pages
+(see "predefined processes" on the side bar).
 See the [examples](@id examples) for a couple of applications.
 
 ## [Asking questions](@id ask_questions)
 
-There are three options for asking questions:
-
-1. As a new post in the official [Julia discourse](https://discourse.julialang.org/) and ask a question under the category Specific Domains > Modelling & Simulations, also using `dynamical-systems` as a tag.
-2. As a message in our channel `#dynamics-bridged` in the [Julia Slack](https://julialang.org/slack/) workplace.
-3. By opening an issue directly on the GitHub page of DynamicalSystems.jl while providing a Minimal Working Example. This is the most useful approach when you encounter unexpected behaviour.
+See <https://github.com/JuliaDynamics/.github/blob/main/SUPPORT.md>.
 
 ## Contributing
 
-There are many ways to contribute to ConceptualClimateModels.jl:
-
-1. Just *use it*! Share it with your colleagues if it was useful for you, and report
-   unexpected behaviour if you find any.
-2. Suggest processes that you think should be included in our library. This should be
-   done by opening a new GitHub issue that describes the process and gives references to papers using the method.
-3. Contribute code by adding new documentation examples.
-4. Contribute code by implementing new processes! That is by far the most impactful
-   way to contribute to the library.
-
-Contributed code must be documented per the standards of
-[DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contributors_guide/#Documentation-string-style).
\ No newline at end of file
+See <https://github.com/JuliaDynamics/.github/blob/main/CONTRIBUTING.md>.
diff --git a/docs/src/processes.md b/docs/src/processes.md
deleted file mode 100644
index bf95b35..0000000
--- a/docs/src/processes.md
+++ /dev/null
@@ -1,99 +0,0 @@
-# [Predefined processes](@id 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 ConceptualClimateModels.jl.
-
-## Temperature
-
-```@docs
-BasicRadiationBalance
-```
-
-## Temperature difference
-
-```@docs
-ΔTLinearRelaxation
-ΔTStommelModel
-```
-
-## Longwave radiation
-
-```@docs
-LinearOLR
-LinearClearSkyOLR
-EmissivityStefanBoltzmanOLR
-EmissivityFeedbackTanh
-EmissivitySellers1969
-SoedergrenClearSkyEmissivity
-```
-
-## Shortwave radiation
-
-```@docs
-DirectAlbedoAddition
-CoAlbedoProduct
-SeparatedClearAllSkyAlbedo
-```
-
-## Ice/snow
-
-```@docs
-IceAlbedoFeedback
-```
-
-## Water vapor
-
-```@docs
-saturation_vapor_pressure
-```
-
-## Insolation
-
-```@docs
-AstronomicalForcingDeSaedeleer
-```
-
-## Forcings
-
-```@docs
-CO2Forcing
-```
-
-## Clouds
-
-```@docs
-CloudAlbedoExponential
-CloudAlbedoLinear
-BudykoOLR
-```
-
-## [Generic processes](@id generic_processes)
-
-Processes that do not depend on any particular physical concept and instead provide
-a simple way to create new processes for a given climate variable:
-
-```@docs
-ParameterProcess
-TimeDerivative
-ExpRelaxation
-AdditionProcess
-TanhProcess
-```
-
-## [Default processes](@id default_processes)
-
-The list of default processes that are used by default in [`processes_to_coupledodes`](@ref) if one does not explicitly provide a list of default processes are:
-
-```@setup MAIN
-using ConceptualClimateModels
-struct ShowFile
-    file::String
-end
-function Base.show(io::IO, ::MIME"text/plain", f::ShowFile)
-    write(io, read(f.file))
-end
-```
-```@example MAIN
-ShowFile(joinpath(dirname(pathof(ConceptualClimateModels)), "default.jl")) # hide
-```
\ No newline at end of file
diff --git a/docs/src/submodules/globalmeanebm.md b/docs/src/submodules/globalmeanebm.md
new file mode 100644
index 0000000..87ef923
--- /dev/null
+++ b/docs/src/submodules/globalmeanebm.md
@@ -0,0 +1,83 @@
+# GlobalMeanEBM
+
+```@docs
+GlobalMeanEBM
+```
+
+## Default variables
+
+```@example MAIN
+using ConceptualClimateModels
+GlobalMeanEBM.globalmeanebm_variables
+```
+
+## Temperature
+
+```@docs
+BasicRadiationBalance
+```
+
+## Temperature difference
+
+```@docs
+ΔTLinearRelaxation
+ΔTStommelModel
+```
+
+## Longwave radiation
+
+```@docs
+LinearOLR
+LinearClearSkyOLR
+EmissivityStefanBoltzmanOLR
+EmissivityFeedbackTanh
+EmissivitySellers1969
+SoedergrenClearSkyEmissivity
+```
+
+## Shortwave radiation
+
+```@docs
+DirectAlbedoAddition
+CoAlbedoProduct
+SeparatedClearAllSkyAlbedo
+```
+
+## Ice/snow
+
+```@docs
+IceAlbedoFeedback
+```
+
+## Water vapor
+
+```@docs
+saturation_vapor_pressure
+```
+
+## Insolation
+
+```@docs
+AstronomicalForcingDeSaedeleer
+```
+
+## Forcings
+
+```@docs
+CO2Forcing
+```
+
+## Clouds
+
+```@docs
+CloudAlbedoExponential
+CloudAlbedoLinear
+BudykoOLR
+```
+
+## Default processes
+
+```@example MAIN
+using ConceptualClimateModels.GlobalMeanEBM
+default_processes_eqs(GlobalMeanEBM)
+```
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
index e0b1c32..a22cf44 100644
--- a/docs/src/tutorial.md
+++ b/docs/src/tutorial.md
@@ -1,5 +1,7 @@
 # [Tutorial](@id tutorial)
 
+## Terminology
+
 ConceptualClimateModels.jl follows a process-based modelling approach
 to make differential equation systems from _processes_.
 A _process_ is simply a particular _equation_ defining the dynamics of a climate
@@ -7,12 +9,12 @@ variable, while also explicitly defining _which_ variable the equation defines.
 A vector of processes is composed by the user, and given to the main function [`processes_to_coupledodes`](@ref) which bundles them into a system of equations
 that creates a dynamical system. The dynamical system can then be further analyzed
 in terms of stability properties, multistability, tipping, periodic (or not)
-behavior, and many more aspects, via the DynamicalSystems.jl library (see the examples).
+behavior, chaos, and many more aspects, via the DynamicalSystems.jl library (see the examples).
 
 Note the distinction: a _process_ is _not_ the climate variable
 (such as "clouds" or  "insolation"); rather it is the _exact equation_ that
 defines the behavior of the climate variable, and could itself utilize
-many other already existing climate variables.
+many other already existing climate variables, or introduce new ones.
 In terminology of climate science a _process_ is a generalization
 of the term "parameterization".
 Many different processes may _describe_ the behavior of a particular variable and
@@ -21,7 +23,7 @@ process.
 
 
 !!! note "Familiarity with DynamicalSystems.jl and ModelingToolkit.jl"
-    ConceptualClimateModels.jl builds on [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) for building the equations representing the climate model, and it builds on [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/) to further analyze the models. Familiarity with either package is good to have, and will allow you to faster and better understand the concepts discussed here. Nevertheless familiarity is actually optional as the steps required to use ConceptualClimateModels.jl are all explained in this tutorial.
+    ConceptualClimateModels.jl uses [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) for building the equations representing the climate model via **symbolic expressions**. It then uses [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/) to further analyze the models. Familiarity with either package is good to have (but not mandatory), and will allow you to faster and better understand the concepts discussed here.
 
 
 ## Introductory example
@@ -48,16 +50,16 @@ To create this model with ConceptualClimateModels.jl while providing the least i
 
 ```@example MAIN
 using ConceptualClimateModels
-using ConceptualClimateModels.CCMV
+using ConceptualClimateModels.GlobalMeanEBM # submodule for global mean models
 
 processes = [
     BasicRadiationBalance(),
     LinearOLR(),
     ParameterProcess(α),
-    CO2Forcing(), # note that for default CO2 value this is zero forcing
+    CO2Forcing(), # for default CO2 valu this is zero forcing
 ]
 
-ds = processes_to_coupledodes(processes)
+ds = processes_to_coupledodes(processes, GlobalMeanEBM)
 println(dynamical_system_summary(ds))
 ```
 
@@ -86,17 +88,15 @@ parameters(mtk)
 ```
 
 The symbolic variables and parameters can also be used to query or alter the
-dynamical system. For example, we can obtain, or alter, any parameter by providing the
-symbolic parameter index.
-There are multiple ways to obtain the symbolic index provided we know its name.
+dynamical system. In the equations above we say that the symbolic variable `CO2` was equated to the parameter `CO2_0`.
+There are multiple ways to obtain a symbolic index provided we know its name.
 First, we can re-create the symbolic parameter:
 
 ```@example MAIN
-index = first(@parameters CO2_0)
+index = only(@parameters CO2_0)
 ```
 
-Second, we can use the retrieved MTK model and access its `CO2_0` parameter:
-
+Second, we can use the retrieved MTK model and access its `CO2_0` field, which will return the symbolic variable:
 
 ```@example MAIN
 index = mtk.CO2_0
@@ -123,24 +123,31 @@ set_parameter!(ds, index, 800)
 
 Similarly, we can obtain or alter values corresponding to the dynamic variables,
 or observed functions of the state of the system, using their symbolic indices.
-For example we can obtain the value corresponding to symbolic variable ``T`` by:
+For example we can obtain the value corresponding to symbolic variable ``T`` at the current state of the dynamical system with:
 
 ```@example MAIN
 observe_state(ds, T) # binding `T` already exists in scope
 ```
 
-or obtain the ``OLR`` (outgoing longwave radiation)
+We didn't have to create `T` because when we did `using ConceptualClimateModels.GlobalMeanEBM`, `T` was brought into scope.
+Similarly, we can obtain the ``OLR`` (outgoing longwave radiation), using the `Symbol` representation of the symbolic variable:
 
 ```@example MAIN
 observe_state(ds, :OLR)
 ```
 
+The advantage of using the symbolic variable `OLR` itself instead of its `Symbol` representation, is that symbolic variables can do arbitrary symbolic computations. For example let's say that when we created the equations of the system we didn't have defined a variable representing the expression ``T^2/OLR``. That's not a problem, we can query the expression nevertheless:
+
+```@example MAIN
+observe_state(ds, T^2/OLR) # symbolic expression
+```
+
 Let's unpack the steps that led to this level of automation.
 
 ## Processes
 
 A process is conceptually an _equation the defines a climate variable or observable_.
-All processes that composed the system are then composed into a set of differential equations via [`processes_to_coupledodes`](@ref) (or [`processes_to_mtkmodel`](@ref)) that represent the climate model.
+All processes that define the system are then composed into a set of differential equations via [`processes_to_coupledodes`](@ref) (or [`processes_to_mtkmodel`](@ref)) that represent the climate model.
 
 For example, the process
 ```@example MAIN
@@ -166,7 +173,7 @@ This allows very easily exchanging the way processes are represented by equation
 OLR_process = EmissivityStefanBoltzmanOLR()
 lhs(OLR_process) ~ rhs(OLR_process)
 ```
-then we would have used a Stefan-Boltzmann grey-atmosphere representation for the outgoing longwave radiation.
+then we would have used a Stefan-Boltzmann grey-atmosphere representation for the outgoing longwave radiation. Notice how this introduced an additional variable ``\varepsilon``.
 
 ## Default processes
 
@@ -179,14 +186,26 @@ processes = [
     CO2Forcing(),
 ]
 ```
-there was no process that defined the absorbed solar radiation ``ASR``!
+there was no process that defined for example the absorbed solar radiation ``ASR``!
 
-Well, ConceptualClimateModels.jl has a list of [predefined processes](@ref predefined_processes) that are automatically included in every call to [`processes_to_coupledodes`](@ref).
-The default processes for the ASR is $ASR = S(1-\alpha)$ with $S$ the solar constant.
-The function [`processes_to_coupledodes`](@ref) goes through all processes the user provided and identifies variables that themselves do not have a process.
-It then checks the list of default processes and attempt to assign one to these variables.
+Well, ConceptualClimateModels.jl allows the concept of **default processes**.
+The package exports some **submodules**, and each submodule targets a particular application of conceptual climate models. In this Tutorial we are using the [`GlobalMeanEBM`](@ref) submodule, which provides functionality to model global mean energy balance models.
+
+Each submodule defines and exports its own **list of predefined symbolic variables**.
+When we wrote
+```julia
+using ConceptualClimateModels.GlobalMeanEBM
+```
+we brought into scope all the variables that this (sub)module defines and exports, such as `T, α, OLR`.
 
-If there are no default processes, it makes the variables themselves parameters with the same name but with a subscript 0.
+Each (sub)module and provides a list of **predefined processes for its predefined symbolic variables**.
+These predefined default processes are loaded automatically when we provide the (sub)module as a second argument to `processes_to_coupledodes`, which we did above.
+In this (sub)module, the default process for the $ASR$ is $ASR = S(1-\alpha)$ with $S$ the solar constant.
+
+The function [`processes_to_coupledodes`](@ref) goes through all processes the user provided and identifies variables that themselves do not have a process.
+It then checks the list of default processes and attempts to assign one to these variables.
+If there are no default processes for some variables, it makes the variables themselves parameters with the same name but with a subscript 0 if the variables
+have a default value.
 
 For example, let's assume that we completely remove default processes and we don't specify
 a process for the albedo ``α``:
@@ -195,14 +214,16 @@ a process for the albedo ``α``:
 processes = [
     BasicRadiationBalance(),
     LinearOLR(),
-    CO2Forcing(), # note that for default CO2 values this is zero forcing
+    CO2Forcing(), # for default CO2 value this is zero forcing
     ASR ~ S*(1-α), # add the processes for ASR, but not for S or α
 ]
+```
 
-# note the empty list as 2nd argument, which is otherwise
-# the default processes. Notice also that we make an MTK model
-# (which is the step _before_ converting to a dynamical system)
-mtk = processes_to_mtkmodel(processes, [])
+which we'll give to `processes_to_mtkmodel`. Notice how there is no second argument given, which would normally be the (sub)module to obtain default processes from.
+
+_side note: we use `process_to_mtkmodel` here, as we don't care yet for making a `DynamicalSystem`; just the symbolic model representation_
+```@example MAIN
+mtk = processes_to_mtkmodel(processes)
 # we access the equations directly from the model
 equations(mtk)
 ```
@@ -222,19 +243,17 @@ default_value(mtk.α_0)
 When this automation occurs a warning is thrown:
 
 ```
-┌ Warning: Variable α was introduced in process of variable ASR(t).
-│ However, a process for α was not provided,
+┌ Warning: Variable α(t) was introduced in process of variable ASR(t).
+│ However, a process for α(t) was not provided,
 │ and there is no default process for it either.
 │ Since it has a default value, we make it a parameter by adding a process:
-│ `ParameterProcess(α)`.
-└ @ ProcessBasedModelling ...\ProcessBasedModelling\src\make.jl:65
+└ `ParameterProcess(α)`.
 ```
 
-[`ParameterProcess`](@ref) is the most trivial process: it simply means that the corresponding variable does not have any physical process and rather it is a system
-parameter.
+[`ParameterProcess`](@ref) is the most trivial process: it simply means that the corresponding variable does not have any actual process describing it and rather it is a system parameter.
 
 This automation does not occur if there is no default value.
-For example, variables that can never be dynamic state variables, such as $ASR$, do not have a default value. If we have not assigned a process for $ASR$, the system construction would error instead:
+For example, the $ASR$ variable does not have a default value. If we have not assigned a process for $ASR$, the system construction would error instead:
 
 ```julia
 processes = [
@@ -243,25 +262,25 @@ processes = [
     CO2Forcing(),
 ]
 
-mtk = processes_to_mtkmodel(processes, [])
+mtk = processes_to_mtkmodel(processes)
 ```
 
 ```
-ERROR: ArgumentError: Variable ASR was introduced in process of
-variable T(t). However, a process for ASR was not provided,
-there is no default process for ASR, and ASR doesn't have a default value.
-Please provide a process for variable ASR.
+ERROR: ArgumentError: Variable ASR(t) was introduced in process of
+variable T(t). However, a process for ASR(t) was not provided,
+there is no default process for ASR(t), and ASR(t) doesn't have a default value.
+Please provide a process for variable ASR(t).
 ```
 
-These warnings and errors are always perfectly informative.
+These warnings and errors are always "perfectly" informative.
 They tell us exactly which variable does not have a process, and exactly which other process introduced the process-less variable first.
-This makes the modelling experience stress-free, especially when large and complex
+This drastically improves the modelling experience, especially when large and complex
 models are being created.
 
 ## Adding your own processes
 
-ConceptualClimateModels.jl provides an increasing list of
-[predefined processes](@ref predefined_processes) that you can use out of
+Each of the submodules of ConceptualClimateModels.jl provides an increasing list of
+predefined processes that you can use out of
 the box to compose climate models. The predefined processes all come from existing
 literature and cite their source via BiBTeX.
 
@@ -279,31 +298,20 @@ x_process = x ~ 0.5*T^2 # x is just a function of temperature
 A more re-usable approach however is to create a function that generates a process
 or create a new process type as we describe in [making new processes](@ref new_processes).
 
-## [Premade symbolic variable instances](@id global_vars)
+## A note on symbolic variable instances
 
-You might be wondering, when we wrote the equation `ASR ~ S*(1-α)` for the $ASR$ process,
-or when we wrote `x ~ 0.5 * T^2`, where did the variable bindings `ASR, S, α` come from?
-For convenience, ConceptualClimateModels.jl defines some symbolic variables
-for typical climate quantities and assigns default processes to them.
-we brought all of these into scope when we did
+Recall that when we wrote
 ```julia
-using ConceptualClimateModels.CCMV
+using ConceptualClimateModels.GlobalMeanEBM
 ```
-where `CCMV` (standing for ConceptualClimateModels Variables) is a submodule
-that defines and exports the variables. We list all of these
-[below](@ref list_vars). These default variables are used throughout the library as the
-default variables in [predefined processes](@ref predefined_processes).
-When going through documentation strings of [predefined processes](@ref predefined_processes),
-such as [`BasicRadiationBalance`](@ref),
-you will notice that the function call signatures are like:
+at the start of this tutorial, we brought into scope variables that this (sub)module defines and exports, such as `T, α, OLR`. They are listed on the (sub)module's documentation page, and are used in that module's default processes.
 
+In some predefined processes documentation you will notice call signatures like
 ```julia
 BasicRadiationBalance(; T, f, kwargs...)
 ```
 There are keywords that do not have an assignment like `T, f` above.
-These always represent climate
-variables, never parameters, and for the variables they use the [existing predefined
-climate variables](@ref list_vars).
+This means that use the (sub)module's predefined variables.
 
 Crucially, these default variables are _symbolic variables_. They are defined as
 ```julia
@@ -313,6 +321,7 @@ Crucially, these default variables are _symbolic variables_. They are defined as
 end
 ```
 which means that expressions that involve them result in symbolic expressions,
+for example
 ```@example MAIN
 A2 = 0.5
 B2 = 0.5
@@ -322,11 +331,11 @@ OLR2 = A2 + B2*T
 In contrast, if we did instead
 ```@example MAIN
 T2 = 0.5 # _not_ symbolic!
-OLR2 = A2 + B2*T2
+OLR3 = A2 + B2*T2
 ```
-This `OLR2` is _not_ a symbolic expression and _cannot_ be used to represent a process.
+This `OLR3` is _not_ a symbolic expression and _cannot_ be used to represent a process.
 
-You can use your own variables for any of the [predefined processes](@ref predefined_processes)
+You can use your own variables in any predefined processes.
 You can define them by doing
 ```@example MAIN
 @variables begin
@@ -339,7 +348,8 @@ process = BasicRadiationBalance(T = T1_tropics)
 ```
 
 Defining variables with the extra `bounds, description` annotations is
-useful for integrating with the rest of the functionality of the library.
+useful for integrating with the rest of the functionality of the library,
+and therefore it is strongly recommended.
 
 !!! warn "Custom variables need to be assigned everywhere!"
     Above we assigned `T1_tropics` as the temperature variable.
@@ -347,45 +357,17 @@ useful for integrating with the rest of the functionality of the library.
     `OLR` variable by also providing the processes
     `LinearOLR(T = T1_tropics)` (for example).
 
-### [List of premade variables](@id list_vars)
 
-The premade variables are not exported by default.
-To bring them into global scope you have to do:
+## Default values, limits, etc.
 
-```@example MAIN
-using ConceptualClimateModels
-using ConceptualClimateModels.CCMV
-```
-and then use them,
-```@example MAIN
-T, q, OLR
-```
-
-All the premade variables are:
-```@example MAIN
-PREDEFINED_CCM_VARIABLES
-```
-
-### Default values, limits, etc.
-
-All premade variables have a default value, a description, and plausible physical bounds.
+All predefined variables that could be dynamic variables (i.e., could have a time derivative applied to them) have a default value, a description, and plausible physical bounds.
 
 To obtain the default value, use `default_value(x)`. For the description,
-use `getdescription(x)`. For the bounds, use:
-
-```@docs
-physically_plausible_limits(::Any)
-```
-
-e.g.,
-
-```@example MAIN
-physically_plausible_limits(T)
-```
+use `getdescription(x)`. For the bounds, see [`plausible_limits`](@ref).
 
 ## Automatically named parameters
 
-The majority of [predefined processes](@ref predefined_processes) create symbolic
+The majority of predefined processes create symbolic
 parameters that are automatically named based on the variables governing the processes.
 This default behaviour can be altered in two ways.
 
@@ -427,17 +409,17 @@ parameters(mtk)
 
 ## Integration with DynamicalSystems.jl
 
-ConceptualClimateModels.jl integrates with [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/) by providing initial condition
+ConceptualClimateModels.jl integrates with [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/). It provides a summary function  [`dynamical_system_summary`](@ref), and also provides initial condition
 sampling to use when e.g., finding attractors and their basin fractions with
-`DynamicalSystems.basins_fractions`, and with the function [`dynamical_system_summary`](@ref).
+`DynamicalSystems.basins_fractions`, via the function [`plausible_limits`](@ref).
 Moreover, since all dynamical systems generated by ConceptualClimateModels.jl have
 symbolic bindings, one can use the symbolic variables in e.g., [interactive GUI exploration](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/visualizations/#Interactive-or-animated-trajectory-evolution)
 or to access or set the parameters of the system.
 
 ```@docs
-physically_plausible_limits(::DynamicalSystem)
-physically_plausible_ic_sampler
-physically_plausible_grid
+plausible_limits
+plausible_ic_sampler
+plausible_grid
 dynamical_system_summary
 ```
 
@@ -452,8 +434,7 @@ processes_to_mtkmodel
 
 To make a new processes you can:
 
-1. Create a _function_ that given some keyword arguments (including which symbolic
-   variables to use) uses one of the existing [generic processes](@ref generic_processes)
+1. Create a _function_ that given some keyword arguments uses one of the existing [generic processes](@ref generic_processes)
    to make and return a process instance. Or, it can return an equation directly,
    provided it satisfies the format of [`processes_to_mtkmodel`](@ref).
    For an example of this, see the source code of [`SeparatedClearAllSkyAlbedo`](@ref)
@@ -466,4 +447,18 @@ To make a new processes you can:
 
 ```@docs
 Process
+register_default_process!
+```
+
+## [Generic processes](@id generic_processes)
+
+Processes that do not depend on any particular physical concept and instead provide
+a simple way to create new processes for a given climate variable:
+
+```@docs
+ParameterProcess
+TimeDerivative
+ExpRelaxation
+AdditionProcess
+SigmoidProcess
 ```
diff --git a/src/ConceptualClimateModels.jl b/src/ConceptualClimateModels.jl
index 7a93181..6b8985d 100644
--- a/src/ConceptualClimateModels.jl
+++ b/src/ConceptualClimateModels.jl
@@ -14,26 +14,10 @@ using ProcessBasedModelling: rhs, lhs_variable, lhs
 
 import NaNMath # so we can use them for sqrt, log, etc. in MTK
 
-include("constants.jl")
-
-include("variables.jl")
-using .CCMV
-
-include("statespace.jl")
-
-# include all processes first
-for (root, dirs, files) in walkdir(joinpath(@__DIR__, "processes"))
-    for file in files
-        include(joinpath(root, file))
-    end
-end
-
-include("default.jl")
-export DEFAULT_CCM_PROCESSES
-
 include("processes_advanced.jl")
+include("statespace.jl")
 include("dynamicalsystems.jl")
 
-export DEFAULT_PROCESSES
+include("GlobalMeanEBM/GlobalMeanEBM.jl")
 
 end # module
\ No newline at end of file
diff --git a/src/GlobalMeanEBM/GlobalMeanEBM.jl b/src/GlobalMeanEBM/GlobalMeanEBM.jl
new file mode 100644
index 0000000..0f468d9
--- /dev/null
+++ b/src/GlobalMeanEBM/GlobalMeanEBM.jl
@@ -0,0 +1,47 @@
+"""
+    GlobalMeanEBM
+
+Submodule of `ConceptualClimateModels` which provides processes
+useful in creating global mean energy balance models.
+It is inspired by Budyko-Sellers-Ghil type models (without the latitudinal dependence).
+The mean does not have to be global, it can be hemispheric, or for other smaller regions.
+"""
+module GlobalMeanEBM
+using ConceptualClimateModels
+import NaNMath
+
+include("variables.jl")
+
+# include all processes
+for (root, dirs, files) in walkdir(joinpath(@__DIR__))
+    for file in files
+        file in ("variables.jl", "GlobalMeanEBM.jl") && continue
+        include(joinpath(root, file))
+    end
+end
+
+function __init__()
+    # note; all variables that do not have a process here
+    # become parameters via `ParameterProcess` by default.
+    register_default_process!.(
+        [
+            BasicRadiationBalance(),
+            # shortwave
+            ASR ~ S*(1 - α)*solar_constant,
+            IceAlbedoFeedback(),
+            DirectAlbedoAddition(), # Albedo uses fact that cloud albedo is defined as additive
+            # longwave
+            BudykoOLR(),
+            CO2Forcing(), # for default CO2 values this is zero forcing
+            AbsoluteHumidityIGLCRH(),
+            # misc
+            ΔTLinearRelaxation(),
+        ],
+        Ref(GlobalMeanEBM)
+    )
+end
+
+
+end # Module
+
+export GlobalMeanEBM
\ No newline at end of file
diff --git a/src/processes/albedo.jl b/src/GlobalMeanEBM/albedo.jl
similarity index 100%
rename from src/processes/albedo.jl
rename to src/GlobalMeanEBM/albedo.jl
diff --git a/src/processes/clouds/albedo.jl b/src/GlobalMeanEBM/clouds/albedo.jl
similarity index 100%
rename from src/processes/clouds/albedo.jl
rename to src/GlobalMeanEBM/clouds/albedo.jl
diff --git a/src/processes/clouds/longwave.jl b/src/GlobalMeanEBM/clouds/longwave.jl
similarity index 100%
rename from src/processes/clouds/longwave.jl
rename to src/GlobalMeanEBM/clouds/longwave.jl
diff --git a/src/constants.jl b/src/GlobalMeanEBM/constants.jl
similarity index 100%
rename from src/constants.jl
rename to src/GlobalMeanEBM/constants.jl
diff --git a/src/processes/forcing.jl b/src/GlobalMeanEBM/forcing.jl
similarity index 100%
rename from src/processes/forcing.jl
rename to src/GlobalMeanEBM/forcing.jl
diff --git a/src/processes/ice.jl b/src/GlobalMeanEBM/ice.jl
similarity index 92%
rename from src/processes/ice.jl
rename to src/GlobalMeanEBM/ice.jl
index 40e35c9..e1061db 100644
--- a/src/processes/ice.jl
+++ b/src/GlobalMeanEBM/ice.jl
@@ -37,14 +37,14 @@ end
 ProcessBasedModelling.lhs_variable(a::IceAlbedoFeedback) = a.α_ice
 ProcessBasedModelling.timescale(a::IceAlbedoFeedback) = a.τ
 
+# TODO: Re-write this now that sigmoid allows more freedom.
 function ProcessBasedModelling.rhs(a::IceAlbedoFeedback)
     y = a.α_ice
     suffixes = ["max", "min", "Tscale", "Tfreeze"]
     values = (a.max, a.min, a.Tscale, a.Tfreeze)
     mtk_pars = map((val, s) -> new_derived_named_parameter(y, val, s), values, suffixes)
     Tref = mtk_pars[4] - mtk_pars[3]/2
-    # don't use `TanhProcess`, use the function instead
-    iproc = ExpRelaxation(y, tanh_expression(a.T, mtk_pars[1:3]..., Tref), a.τ)
+    iproc = ExpRelaxation(y, ConceptualClimateModels.sigmoid_expression(a.T, mtk_pars[1:3]..., Tref), a.τ)
     return ProcessBasedModelling.rhs(iproc)
 end
 # we don't need to extend `lhs`, `ProcessBasedModelling` takes care of that.
\ No newline at end of file
diff --git a/src/processes/insolation.jl b/src/GlobalMeanEBM/insolation.jl
similarity index 100%
rename from src/processes/insolation.jl
rename to src/GlobalMeanEBM/insolation.jl
diff --git a/src/processes/longwave/emissivity.jl b/src/GlobalMeanEBM/longwave/emissivity.jl
similarity index 98%
rename from src/processes/longwave/emissivity.jl
rename to src/GlobalMeanEBM/longwave/emissivity.jl
index b391eb7..c856226 100644
--- a/src/processes/longwave/emissivity.jl
+++ b/src/GlobalMeanEBM/longwave/emissivity.jl
@@ -10,7 +10,7 @@ similar to [`EmissivitySellers1969`](@ref).
 In essence this is a [`TanhProcess`](@ref) with the given keywords as parameters.
 """
 function EmissivityFeedbackTanh(; ε=ε, T=T, left = 0.5, right = 0.4, rate = 0.5, Tref = 288.0)
-    return TanhProcess(ε, T, left, right, rate, Tref)
+    return SigmoidProcess(ε, T, left, right, rate, Tref)
 end
 
 """
diff --git a/src/processes/longwave/olr.jl b/src/GlobalMeanEBM/longwave/olr.jl
similarity index 100%
rename from src/processes/longwave/olr.jl
rename to src/GlobalMeanEBM/longwave/olr.jl
diff --git a/src/processes/temperature/gmt.jl b/src/GlobalMeanEBM/temperature/gmt.jl
similarity index 100%
rename from src/processes/temperature/gmt.jl
rename to src/GlobalMeanEBM/temperature/gmt.jl
diff --git a/src/processes/temperature/tempdiff.jl b/src/GlobalMeanEBM/temperature/tempdiff.jl
similarity index 100%
rename from src/processes/temperature/tempdiff.jl
rename to src/GlobalMeanEBM/temperature/tempdiff.jl
diff --git a/src/processes/thermodynamics.jl b/src/GlobalMeanEBM/thermodynamics.jl
similarity index 100%
rename from src/processes/thermodynamics.jl
rename to src/GlobalMeanEBM/thermodynamics.jl
diff --git a/src/variables.jl b/src/GlobalMeanEBM/variables.jl
similarity index 75%
rename from src/variables.jl
rename to src/GlobalMeanEBM/variables.jl
index f07c571..732630b 100644
--- a/src/variables.jl
+++ b/src/GlobalMeanEBM/variables.jl
@@ -1,15 +1,10 @@
-function physically_plausible_limits end
-
-module CCMV
-
-using ConceptualClimateModels
 # These are quantities that change with time. Hence, they are either state variables
 # or observables of the state variables. They have given names and are utilized
 # in the rest of the files to make the equations that compose the dynamical system.
 # In the final generated energy balance model it is not necessary that all of these
 # will exist in the equations.
 
-const PREDEFINED_CCM_VARIABLES = @variables begin
+globalmeanebm_variables = @variables begin
     (T(t) = 290.0),      [bounds = (200.0, 350.0), description = "temperature, in Kelvin"]
     (S(t) = 1.0),        [bounds = (0.8, 1.2), description = "insolation, normalized to units of the solar constant"]
     (f(t) = 0.0),        [bounds = (-0.1, 0.1), description = "external forcing, in W/m²"]
@@ -28,20 +23,4 @@ const PREDEFINED_CCM_VARIABLES = @variables begin
     (ASR(t)), [description = "absorved shortwave radiation, in W/m²"]
 end
 
-export PREDEFINED_CCM_VARIABLES
 export T, S, f, α, α_ice, α_cloud, ΔT, ΔS, ε, ℓ, C, CO2, OLR, ASR, q
-
-function ConceptualClimateModels.physically_plausible_limits(var::String)::Tuple{Float64, Float64}
-    if var[1] == 'T'
-        return (200, 350)
-    elseif var[1] == 'α' || var[1] == 'ε' || var[1] == 'C'
-        return (0, 1)
-    else
-        error("""
-        Unpsecified plausible physical limits for $(var): it has no defined bounds or
-        a default variable. You need to redefine the variable to have either of the two.
-        """)
-    end
-end
-
-end
\ No newline at end of file
diff --git a/src/processes/water_vapor.jl b/src/GlobalMeanEBM/water_vapor.jl
similarity index 100%
rename from src/processes/water_vapor.jl
rename to src/GlobalMeanEBM/water_vapor.jl
diff --git a/src/default.jl b/src/default.jl
deleted file mode 100644
index 8edc28a..0000000
--- a/src/default.jl
+++ /dev/null
@@ -1,18 +0,0 @@
-# note; all variables that do not have a process here
-# become parameters via `ParameterProcess` by default.
-DEFAULT_CCM_PROCESSES = [
-    BasicRadiationBalance(),
-    # shortwave
-    ASR ~ S*(1 - α)*solar_constant,
-    IceAlbedoFeedback(),
-    DirectAlbedoAddition(), # Albedo uses fact that cloud albedo is defined as additive
-    S ~ 1, # don't make insolation a parameter by default
-    # longwave
-    BudykoOLR(),
-    CO2Forcing(), # for default CO2 values this is zero forcing
-    ParameterProcess(CO2),
-    AbsoluteHumidityIGLCRH(),
-    # misc
-    C ~ default_value(C),
-    ΔTLinearRelaxation(),
-]
diff --git a/src/dynamicalsystems.jl b/src/dynamicalsystems.jl
index d67a946..af52fb6 100644
--- a/src/dynamicalsystems.jl
+++ b/src/dynamicalsystems.jl
@@ -5,7 +5,7 @@ export all_equations
 DEFAULT_DIFFEQ = DynamicalSystemsBase.DEFAULT_DIFFEQ
 
 """
-    processes_to_coupledodes(processes, default = DEFAULT_PROCESSES; kw...)
+    processes_to_coupledodes(processes [, default]; kw...)
 
 Convert a given `Vector` of processes to a `DynamicalSystem`, in particular `CoupledODEs`.
 All processes represent symbolic equations managed by ModelingToolkit.jl.
@@ -27,7 +27,7 @@ or see the online [Tutorial](@ref).
   This accelerates parameter access, assuming all parameters are of the same type.
 - `kw...`: all other keywords are propagated to `processes_to_mtkmodel`.
 """
-function processes_to_coupledodes(proc, default = DEFAULT_CCM_PROCESSES;
+function processes_to_coupledodes(proc, default = [];
         diffeq = DEFAULT_DIFFEQ, inplace = nothing, split::Bool = false, kwargs...
     )
     sys = processes_to_mtkmodel(proc, default; kwargs...)
@@ -55,20 +55,20 @@ end
 Return a printable/writable string containing a summary of `ds`,
 which outlines its current status and lists all symbolic
 equations and parameters that constitute the system, if a referrence to a
-ModelingToolkit.jl exists in `ds`.
+ModelingToolkit.jl model exists in `ds`.
 """
 function dynamical_system_summary(ebm::DynamicalSystem)
     if !DynamicalSystemsBase.has_referrenced_model(ebm)
         summary = sprint(show, MIME"text/plain"(), ebm)
         return summary
     end
-    # Else, use symbolic stuff
+    # Else, use symbolic stuff; first print the first 4 lines
     summary = keepfirstlines(sprint(show, MIME"text/plain"(), ebm), 4)*"\n"
     model = referrenced_sciml_model(ebm)
     summary *= "\nwith equations:\n"*skipfirstline(sprint(show, MIME"text/plain"(), all_equations(model)))
-    # TODO: here add dump metadata
-#     See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref),
-# [`ModelingToolkit.dump_unknowns`](@ref).
+    # We can use
+    # [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref),
+    # [`ModelingToolkit.dump_unknowns`](@ref), as well; but this works fine.
     summary *= "\n\nand parameter values:\n"*skipfirstline(sprint(show, MIME"text/plain"(), named_current_parameters(ebm)))
     return summary
 end
diff --git a/src/processes_advanced.jl b/src/processes_advanced.jl
index 170f40a..4b7b798 100644
--- a/src/processes_advanced.jl
+++ b/src/processes_advanced.jl
@@ -1,27 +1,28 @@
-# TODO: Rename to sigmoid
-# TODO: Logistic function is faster in Julia on windows 10
-export TanhProcess
+# TODO: Logistic function is faster than tanh in Julia on windows 10
+export SigmoidProcess
 
 """
-    TanhProcess(variable, driver, left, right, scale, reference) <: Process
-    TanhProcess(variable, driver; left, right, scale, reference) <: Process
+    SigmoidProcess <: Process
+    SigmoidProcess(variable, driver; left, right, scale, reference)
 
-A common process for when a `variable` has a tanh-dependence on a `driver` variable.
+A common process for when a `variable` has a sigmoidal dependence on a `driver` variable.
 The rest of the input arguments should be real numbers or `@parameter` named parameters.
 
-The process creates the expression:
+The process creates a sigmoidal relationship based on the `tanh` function:
 ```
 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`.
+i.e., the variable goes from value `left` to value `right` as `driver` increases
+over a range of `scale` (centered at `reference`).
+Instead of `reference` you may provide `start` or `finish` keywords,
+which make `reference = start + scale/2` or `reference = finish - scale/2` respectively.
 
 If the values given to the parameters of the expression are real numbers, they become
 named parameters prefixed with the name of `variable`, then the name of the `driver`,
-and then `_tanh_left`, `_tanh_right`, `_tanh_rate` and `_tanh_ref` respectively.
+and then `_sigmoid_left`, `_sigmoid_right`, `_sigmoid_rate` and `_sigmoid_ref` respectively.
 Use `LiteralParameter` for parameters you do not wish to rename.
 """
-struct TanhProcess <: Process
+struct SigmoidProcess <: Process
     variable
     driver_variable
     left
@@ -29,21 +30,32 @@ struct TanhProcess <: Process
     scale
     reference
 end
-TanhProcess(variable, driver; left=0, right=1, scale=1, reference=0) =
-TanhProcess(variable, driver, left, right, scale, reference)
 
-function ProcessBasedModelling.rhs(p::TanhProcess)
+function SigmoidProcess(variable, driver; left=0, right=1, scale=1,
+    reference=nothing, start=nothing, finish = nothing)
+    if count(!isnothing, (reference, start, finish)) > 1
+        error("Exactly one of `reference, start, finish` must be given!")
+    end
+    if !isnothing(start)
+        reference = start + scale/2
+    elseif !isnothing(finish)
+        reference = finish - scale/2
+    end
+    return SigmoidProcess(variable, driver, left, right, scale, reference)
+end
+
+function ProcessBasedModelling.rhs(p::SigmoidProcess)
     y = p.variable
     x = p.driver_variable
     # Create the names for everything
     (; left, right, scale, reference) = p
     values = left, right, scale, reference
     xname = string(ModelingToolkit.getname(x))
-    suffixes = (xname*"_") .* ["tanh_left", "tanh_right", "tanh_scale", "tanh_ref"]
+    suffixes = (xname*"_") .* ["sigmoid_left", "sigmoid_right", "sigmoid_scale", "sigmoid_ref"]
     mtk_vars = map((val, s) -> new_derived_named_parameter(y, val, s; prefix = false), values, suffixes)
     # pass them to the regular tanh Julia function
-    return tanh_expression(x, mtk_vars...)
+    return sigmoid_expression(x, mtk_vars...)
 end
-function tanh_expression(T, left, right, scale, reference)
+function sigmoid_expression(T, left, right, scale, reference)
     return left + (right - left)*(1 + tanh(2(T - reference)/(scale)))*0.5
 end
diff --git a/src/statespace.jl b/src/statespace.jl
index 8548261..f19382d 100644
--- a/src/statespace.jl
+++ b/src/statespace.jl
@@ -1,55 +1,56 @@
 # This function is only meaningful for dynamic variables!
 """
-    physically_plausible_limits(x)
+    plausible_limits(x::Num)
 
 Return a tuple (min, max) of plausible limiting values for the variable `x`.
-If the variable does not have defined `bounds` metadata, then the default value ± 20% is used.
-If there is no default value, a heuristic is tried, and an error is thrown if it fails.
+If `x` is defined with the `bounds` metadata, this is returned as the limits.
+Else, if `x` has a default value, the limits are this value ± 20%.
+Else, if there is no default value, an error is thrown.
 """
-function physically_plausible_limits(var)
+function plausible_limits(var)
     if ModelingToolkit.hasbounds(var)
-        return getbounds(var)
+        return ModelingToolkit.getbounds(var)
     elseif !isnothing(default_value(var))
         return (0.8default_value(var), 1.2default_value(var))
     else
-        return CCMV.physically_plausible_limits(string(ModelingToolkit.getname(var)))
+        return throw(ArgumentError("""
+        Variable $(var) does not have any defined limits.
+        Define it with a `bounds` metadata or a default value.
+        """))
     end
 end
 
 """
-    physically_plausible_limits(ds::DynamicalSystem [, idxs])
+    plausible_limits(ds::DynamicalSystem [, idxs])
 
 Return a vector of limits (min, max) for each dynamic state variable in `ds`.
 Optionally provide the `idxs` of the variables to use as a vector of `Symbol`s
 or a vector of symbolic variables present in the referrenced MTK model of `ds`.
+
+See also [`plausible_grid`](@ref), [`plausible_ic_sampler`](@ref).
 """
-function physically_plausible_limits(ds::DynamicalSystem, idxs = nothing)
+function plausible_limits(ds::DynamicalSystem, idxs = nothing)
     model = referrenced_sciml_model(ds)
     vars = ModelingToolkit.unknowns(model)
-    if idxs isa Nothing
-        vars = ModelingToolkit.unknowns(model)
-    elseif idxs isa Vector{Symbol}
+    if idxs isa Vector{Symbol}
         allvars = ModelingToolkit.unknowns(model)
         allnames = ModelingToolkit.getname.(allvars)
         j = [findfirst(isequal(i), allnames) for i in idxs]
         vars = allvars[j]
-    else
-        vars = idxs
     end
-    minmax = physically_plausible_limits.(vars)
-    return minmax
+    return plausible_limits.(vars)
 end
 
 """
-    physically_plausible_ic_sampler(ds::DynamicalSystem)
+    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`.
+[`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)
+function plausible_ic_sampler(ds::DynamicalSystem)
+    minmax = plausible_limits(ds)
     mini = getindex.(minmax, 1)
     maxi = getindex.(minmax, 2)
     sampler, = statespace_sampler(HRectangle(mini, maxi))
@@ -57,15 +58,15 @@ function physically_plausible_ic_sampler(ds::DynamicalSystem)
 end
 
 """
-    physically_plausible_grid(ds::DynamicalSystem, n = 101)
+    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`.
+[`plausible_limits`](@ref)(`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)
+function plausible_grid(ds::DynamicalSystem, n = 101)
+    minmax = plausible_limits(ds)
     if n isa Integer
         ranges = [range(mm[1], mm[2]; length = n) for mm in minmax]
     elseif n isa AbstractVector
@@ -74,6 +75,4 @@ function physically_plausible_grid(ds::DynamicalSystem, n = 101)
     return Tuple(ranges)
 end
 
-export physically_plausible_ic_sampler,
-    physically_plausible_limits,
-    physically_plausible_grid
\ No newline at end of file
+export plausible_ic_sampler, plausible_limits, plausible_grid
\ No newline at end of file
diff --git a/test/bifkit.md b/test/bifkit.md
deleted file mode 100644
index 978d4e9..0000000
--- a/test/bifkit.md
+++ /dev/null
@@ -1,48 +0,0 @@
-### Using BifurcationKit.jl
-
-Because for this example having _unstable_
-fixed point branches is crucial, we will perform the alternative continuation
-provided by BifurcationKit.jl. BifurcationKit.jl does not provide most of the conveniences
-that DynamicalSystems.jl does. For example, it does not integrate well enough with
-DifferentialEquations.jl (to allow passing `ODEProblem` which is created by `DynamicalSystem`).
-It also does not allow indexing parameters by their symbolic bindings.
-
-Due to this, we need to change the dynamical system generation to disable [parameter splitting](https://github.com/SciML/ModelingToolkit.jl/issues/2482#issuecomment-1959330130),
-so that the parameter container becomes a simple vector of numbers
-that can be indexed by the integers:
-
-```@example MAIN
-stommel = processes_to_coupledodes(processes; split = false, inplace = true)
-```
-
-We now need to find out which integer corresponds the the parameters of interest,
-which we do by seeing the list:
-
-```@example MAIN
-mtk = referrenced_sciml_model(stommel)
-pnames = parameters(mtk)
-hcat(eachindex(pnames), pnames)
-```
-
-from where we obtain that index 1 corresponds to `r_η`.
-
-From here, we follow the [introductory tutorial](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/gettingstarted/) of BifurcationKit.jl
-and first make a `BifurcationProblem`:
-
-```@example MAIN
-import BifurcationKit
-
-vf = dynamic_rule(stommel) # vector field, the dynamic rule
-u = current_state(stommel) # initial guess for bifurcation
-p = current_parameters(stommel) # parameter container
-i = (BifurcationKit.@lens _[1]) # BifurcationKit.jl doesn't allow normal indices, only some weird "lenses"
-
-prob = BifurcationKit.BifurcationProblem(vf, u, p, i)
-```
-
-Great, that took a while but now we can finally perform the continuation
-
-```@example MAIN
-opts = BifurcationKit.ContinuationPar(p_min = 0.0, p_max = 5.0)
-br = BifurcationKit.continuation(prob, BifurcationKit.PALC(), opts)
-```
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 1934d7c..6343f44 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,5 @@
 using ConceptualClimateModels
-using ConceptualClimateModels.CCMV
+using ConceptualClimateModels.GlobalMeanEBM
 using Test
 
 @testset "all processes" begin
@@ -32,12 +32,12 @@ p1 = [
     EmissivityLogConcentration(ε = ε4),
 ]
 
-ds = processes_to_coupledodes(p1)
+ds = processes_to_coupledodes(p1, GlobalMeanEBM)
 mtk = referrenced_sciml_model(ds)
 eqs = all_equations(mtk)
 
 # This also tests a bunch of default processes
-for var in (T, ε, ε1, ε2, ε3, ε4, q, :ε1_T_tanh_ref, mtk.ε_0)
+for var in (T, ε, ε1, ε2, ε3, ε4, q, :ε1_T_sigmoid_ref, mtk.ε_0)
     test_symbolic_var(eqs, var)
 end
 
@@ -68,7 +68,7 @@ p2 = [
     ΔTStommelModel(ΔT = ΔT1, ΔS = ΔS1),
 ]
 
-ds = processes_to_coupledodes(p2)
+ds = processes_to_coupledodes(p2, GlobalMeanEBM)
 mtk = referrenced_sciml_model(ds)
 for var in (T, C, OLR1, OLR2, :α_bg, a1, a2, a3, a4, ΔS1, ΔT)
     if has_symbolic_var(mtk, var)
@@ -85,8 +85,6 @@ end
 
 @testset "dynamical systems integration" begin
     using DynamicalSystems
-    using ConceptualClimateModels
-    using ConceptualClimateModels.CCMV
 
     budyko_processes = [
         BasicRadiationBalance(),
@@ -98,12 +96,12 @@ end
         # absorbed solar radiation has a default process
     ]
 
-    budyko = processes_to_coupledodes(budyko_processes)
+    budyko = processes_to_coupledodes(budyko_processes, GlobalMeanEBM)
 
-    grid = physically_plausible_grid(budyko)
+    grid = plausible_grid(budyko)
     mapper = AttractorsViaRecurrences(budyko, grid)
     rfam = RecurrencesFindAndMatch(mapper)
-    sampler = physically_plausible_ic_sampler(budyko)
+    sampler = plausible_ic_sampler(budyko)
 
     set_parameter!(budyko, :ε_0, 0.5)