Update docs

dev_notes/todo.txt

@@ -3,10 +3,16 @@
 	*** High pri ***
 
 
-		There is a pretty bad order-of-module load error in nivafjord_lake. The model can be wrong or not load if airsea_fpv (or airsea) is loaded in the outer and not the inner model!!!!!
-			The energy loss in the simplytox version seems unrelated. What causes it to happen here and not in the other version though?? Seems unrelated to the river too.
-
-
+		There is a pretty bad order-of-module-load error in nivafjord_lake. The model compilation crashes if airsea_fpv if is loaded in the outer and not the inner model!!!!!
+			It doesn't detect index set dependencies correctly! Sort error? Should still have been detected earlier.
+			It doesn't know that surf.ced needs to index over basin_idx even though it looks up stab, which again depends on top_water.temp.
+			But it depends on when the module was loaded, which it should not at all!
+
+
+		Need some form of way to write parameter checks at the start.
+			- Landscape units sum to 1
+			- NIVAFjord areas are decreasing with depth (causes errors all the time when two layers have the same area!)
+			- etc.
 
 		Add forums on github page? (github Discussions)
 
@@ -20,7 +26,12 @@
 			Maybe factor out common phyto module for EasyChem and FjordChem.
 			Also for O2 (found a difficulty with that earlier, but maybe we could resolve it?)
 
-			NIVAFjord!
+			NIVAFjord
+
+				Go to the simpler version from MyLake of parametrizing vertical mixing
+					Mathematically the same, just removes unnecessary parameters
+					Allows us to compare with MyLake setups.
+					See guide chapter 07
 
 				Finish biochemistry.
 					Do we need separate NH4? Do we need Si?

docs/existingmodels/autogen/auxiliary.md

@@ -13,7 +13,7 @@ Since the modules can be dynamically loaded with different arguments, this docum
 
 See the note on [notation](autogen.html#notation).
 
-The file was generated at 2024-08-22 11:08:06.
+The file was generated at 2024-09-09 12:02:42.
 
 ---
 
@@ -72,9 +72,9 @@ File: [modules/hbv_snow.txt](https://github.com/NIVANorge/Mobius2/tree/main/mode
 
 This is an adaption of the snow module from HBV-Nordic (Sælthun 1995)
 
-[NVE home page](https://www.nve.no/vann-og-vassdrag/vannets-kretslop/analysemetoder-og-modeller/hbv-modellen/)
+[NVE home page](https://www.nve.no/vann-og-vassdrag/vannets-kretsloep/analysemetoder-og-modeller/hbv-modellen/)
 
-[Model description](https://www.uio.no/studier/emner/matnat/geofag/nedlagte-emner/GEO4430/v06/undervisningsmateriale/HBVMOD.PDF)
+[Model description](https://publikasjoner.nve.no/publication/1996/publication1996_07.pdf)
 
 Authors: Magnus D. Norling
 
@@ -428,7 +428,9 @@ Authors: Magnus D. Norling
 | GMT zone | **time_zone** | hr | Only used if sub-daily precision is used. |
 | Elevation | **elev** | m | Only used if Daily average radiation is not provided as a series |
 | **Radiation** | | |  |
-| Use sub-daily precision | **subdaily** |  | Compute hourly averages of solar radiation |
+| Cloud absorption scaling factor | **crs** |  | Used if 'Daily average global radiation' is not provided as a series. Scaling factor for shortwave absorption by clouds |
+| Cloud absorption power factor | **crp** |  | Used if 'Daily average global radiation' is not provided as a series. Power factor for shortwave absorption by clouds |
+| Use sub-daily precision | **subdaily** |  | Compute hourly averages of solar radiation. Only works correctly if the model sampling step is [hr] or lower. |
 
 ### State variables
 
@@ -477,7 +479,7 @@ Unit: W m⁻²
 Value:
 
 $$
-\left(1-0.75\cdot \mathrm{cloud}^{3.4}\right)\cdot \mathrm{g\_rad\_cloud\_free}
+\left(1-\mathrm{crs}\cdot \mathrm{cloud}^{\mathrm{crp}}\right)\cdot \mathrm{g\_rad\_cloud\_free}
 $$
 
 #### **Global radiation**

docs/existingmodels/autogen/easylake.md

@@ -13,7 +13,7 @@ Since the modules can be dynamically loaded with different arguments, this docum
 
 See the note on [notation](autogen.html#notation).
 
-The file was generated at 2024-08-22 11:08:05.
+The file was generated at 2024-09-09 12:02:42.
 
 ---
 

docs/existingmodels/autogen/nivafjord.md

@@ -13,7 +13,7 @@ Since the modules can be dynamically loaded with different arguments, this docum
 
 See the note on [notation](autogen.html#notation).
 
-The file was generated at 2024-08-22 11:08:05.
+The file was generated at 2024-09-09 12:02:42.
 
 ---
 
@@ -96,6 +96,7 @@ Authors: Magnus D. Norling
 | Sea level | **h** | property |
 | Area | **area** | property |
 | Shortwave radiation | **sw** | property |
+|  | **ice_ind** | loc |
 | Fjord vertical | **vert** | connection |
 | Shortwave vertical | **sw_vert** | connection |
 
@@ -112,6 +113,7 @@ Authors: Magnus D. Norling
 | Diminishing rate of additional mixing energy | **hsfact** | m |  |
 | **Layer specific mixing** | | | Distributes like: `layer` |
 | Mixing factor reference | **K0** | m² s⁻¹ | Mixing factor when the B-V frequency is equal to the reference |
+| Mixing factor reference under ice | **K0ice** | m² s⁻¹ |  |
 | **Initial layer physical** | | | Distributes like: `layer` |
 | Initial layer temperature | **init_t** | °C |  |
 | Initial layer salinity | **init_s** |  |  |
@@ -351,7 +353,7 @@ Unit: J
 Value:
 
 $$
-\mathrm{A}\lbrack\mathrm{vert}.\mathrm{top}\rbrack\cdot \sqrt{\frac{\mathrm{stress}^{3}}{\mathrm{layer}.\mathrm{water}.\mathrm{rho}\lbrack\mathrm{vert}.\mathrm{top}\rbrack}}\cdot \mathrm{time}.\mathrm{step\_length\_in\_seconds}
+\;\text{not}\;\mathrm{ice\_ind}\cdot \mathrm{A}\lbrack\mathrm{vert}.\mathrm{top}\rbrack\cdot \sqrt{\frac{\mathrm{stress}^{3}}{\mathrm{layer}.\mathrm{water}.\mathrm{rho}\lbrack\mathrm{vert}.\mathrm{top}\rbrack}}\cdot \mathrm{time}.\mathrm{step\_length\_in\_seconds}
 $$
 
 #### **Sum V above**
@@ -423,7 +425,7 @@ Unit: m s⁻¹
 Value:
 
 $$
-\mathrm{dz\_} = 0.5\cdot \left(\mathrm{dz}+\mathrm{dz}\lbrack\mathrm{vert}.\mathrm{below}\rbrack\right) \\ \href{stdlib.html#basic}{\mathrm{safe\_divide}}\left(\mathrm{K0},\, \mathrm{dz\_}\cdot \left(\frac{\mathrm{Nfreq}}{\mathrm{N0}}\right)^{\mathrm{alpha}}\right)
+\mathrm{K} = \begin{cases}\mathrm{K0ice} & \text{if}\;\mathrm{ice\_ind} \\ \mathrm{K0} & \text{otherwise}\end{cases} \\ \mathrm{dz\_} = 0.5\cdot \left(\mathrm{dz}+\mathrm{dz}\lbrack\mathrm{vert}.\mathrm{below}\rbrack\right) \\ \href{stdlib.html#basic}{\mathrm{safe\_divide}}\left(\mathrm{K},\, \mathrm{dz\_}\cdot \left(\frac{\mathrm{Nfreq}}{\mathrm{N0}}\right)^{\mathrm{alpha}}\right)
 $$
 
 #### **Additional mixing**

docs/existingmodels/autogen/simplycnp.md

@@ -13,7 +13,7 @@ Since the modules can be dynamically loaded with different arguments, this docum
 
 See the note on [notation](autogen.html#notation).
 
-The file was generated at 2024-08-22 11:08:05.
+The file was generated at 2024-09-09 12:02:42.
 
 ---
 

docs/existingmodels/autogen/stdlib.md

@@ -14,7 +14,7 @@ The standard library provides common functions and constants for many models.
 
 See the note on [notation](autogen.html#notation).
 
-The file was generated at 2024-08-22 11:08:06.
+The file was generated at 2024-09-09 12:02:42.
 
 ---
 
@@ -218,7 +218,7 @@ File: [stdlib/physiochemistry.txt](https://github.com/NIVANorge/Mobius2/tree/mai
 
 ### Description
 
-These are simplified functions for computing properties of freshwater at surface pressure. See the SeaWater module for functions that work better in more general cases.
+These are simplified functions for computing properties of freshwater at surface pressure. See the Seawater library for functions that work better in more general cases.
 
 References to be inserted.
 
@@ -524,7 +524,7 @@ File: [stdlib/seawater.txt](https://github.com/NIVANorge/Mobius2/tree/main/stdli
 
 ### Description
 
-This is a library for highly accurate, but more computationally expensive, properties of sea water (taking salinity into account). Se the library "Water utils" for simplified versions of these functions.
+This is a library for highly accurate, but more computationally expensive, properties of sea water (taking salinity into account). Se the library [Water utils](https://nivanorge.github.io/Mobius2/existingmodels/autogen/stdlib.html#water-utils) for simplified freshwater versions of these functions.
 
 The formulas for density are taken from the [Matlab seawater package](http://www.marine.csiro.au/~morgan/seawater).
 
@@ -595,7 +595,7 @@ The implementation is influenced by the one in [SELMA](https://github.com/fabm-m
 **o2_saturation(T : °C, S)** = 
 
 $$
-\mathrm{T\_k} = \left(\left(\mathrm{T}\rightarrow \mathrm{K}\,\right)\Rightarrow 1\right) \\ \left(e^{-135.902+\frac{157570}{\mathrm{T\_k}}-\frac{6.64231\cdot 10^{7}}{\mathrm{T\_k}^{2}}+\frac{1.2438\cdot 10^{10}}{\mathrm{T\_k}^{3}}-\frac{8.62195\cdot 10^{11}}{\mathrm{T\_k}^{4}}-\mathrm{S}\cdot \left(0.017674-\frac{10.754}{\mathrm{T\_k}}+\frac{2140.7}{\mathrm{T\_k}^{2}}\right)}\Rightarrow \mathrm{mmol}\,\mathrm{m}^{-3}\,\right)
+\mathrm{T\_k} = \left(\left(\mathrm{T}\rightarrow \mathrm{K}\,\right)\Rightarrow 1\right) \\ \mathrm{logsat} = -135.902+\frac{157570}{\mathrm{T\_k}}-\frac{6.64231\cdot 10^{7}}{\mathrm{T\_k}^{2}}+\frac{1.2438\cdot 10^{10}}{\mathrm{T\_k}^{3}}-\frac{8.62195\cdot 10^{11}}{\mathrm{T\_k}^{4}}-\mathrm{S}\cdot \left(0.017674-\frac{10.754}{\mathrm{T\_k}}+\frac{2140.7}{\mathrm{T\_k}^{2}}\right) \\ \left(e^{\mathrm{logsat}}\Rightarrow \mathrm{mmol}\,\mathrm{m}^{-3}\,\right)
 $$
 
 **o2_piston_velocity(wind : m s⁻¹, temp : °C)** = 

docs/existingmodels/simply.md

@@ -21,7 +21,7 @@ Examples of potential model uses include:
 3. Exploring the potential response of the system to future environmental change (e.g. climate, land use and management), including potential storm and low-flow dynamics.
 Providing evidence to support decision-making, e.g. to help set water quality and load reduction goals, and advise on means of achieving those goals.
 
-The Simply models originate from \[JacksonBlake17\], where the first version of SimplyP (including what was later named SimplyQ and SimplySed) was presented as a [python program](https://github.com/LeahJB/SimplyP).
+The Simply models originate from \[JacksonBlake17\], where the first version of SimplyP (including what was later named SimplyQ and SimplySed) was implemented as a [python program](https://github.com/LeahJB/SimplyP).
 
 An early version of SimplyC was developed in \[Norling21\].
 

docs/mobius2docs/guide_chapters/07_layered_lake.md

@@ -80,13 +80,13 @@ The mixing coefficient $$K$$ is (as in MyLake) given by
 
 $$
 K_r = \begin{cases}
-K_0 A^{0.56} & \text{If there lake is ice free}\\
+K_0 A^{0.56} & \text{if the lake is ice free}\\
 K_{0,ice}    & \text{otherwise}
-\end{cases},\\
+\end{cases}\\
 K = K_r N^{2\alpha}
 $$
 
-Here $$A$$ is the lake surface area, $$K_0$$, and $$K_{0,ice}$$ are configurable parameters (but with good default values), and $$\alpha=-0.43$$ is a constant.
+Here $$A$$ is the lake surface area, $$K_0$$, and $$K_{0,ice}$$ are configurable parameters (but the default values in the dataset are based on literature), and $$\alpha=-0.43$$ is a constant.
 
 The mixing rate of a layer $$i$$ with the one below it is then
 
@@ -97,20 +97,33 @@ $$
 This is formulated in Mobius2 as
 
 ```python
+var(layer.water.K, [m 2, day-1]) {
+	# Set a to have a dimensionless unit so that we can raise 
+	# it to a power in the empirical formula below.
+	a := (lake.area->[k m 2])=>[],
+	K_ref := {
+		K0*a^0.56   if !lake.ice.ind,
+		K0_ice      otherwise
+	},
+	K_ref*(N2freq=>[])^(2*alpha)
+}
+
 flux(layer.water, vert, [m 3, day-1], "Layer mixing down") { 
 	mdz := 0.5*(dz + dz[vert.below]),
 	K*A[vert.below]/mdz->>
 } @mixing
 ```
 
-The `@mixing` note tells Mobius2 that this flux happens in both directions so that the net exchange of water between the two layers is zero, but that it should still cause mixing of dissolved substances (this includes heat energy). In practice, the ODE system that is generated from this is mathematically equivalent to the finite element discretization of the diffusion equation described in \[SalorantaAndersen07\].
+The `@mixing` note tells Mobius2 that this flux happens in both directions so that the net exchange of water between the two layers is zero, but that it should still cause mixing of dissolved substances (this includes heat energy, see below). In practice, the ODE system that is generated from this is mathematically equivalent to the finite element discretization of the diffusion equation described in \[SalorantaAndersen07\].
 
 ## Surface fluxes, ice and heat
 
-The exchange of heat energy between the lake surface and the atmosphere is entirely taken care of by the AirSea module. It also takes care of ice growth and melt, and it computes evaporation. We won't go into detail about that here. Just note that we model `layer.water.heat` as a dissolved variable, then compute temperature from the heat density (concentration).
+The exchange of heat energy between the lake surface and the atmosphere is entirely taken care of by the AirSea module. It also takes care of ice growth and melt, and it computes evaporation. We won't go into detail about that here. Just note that we model the water heat energy `layer.water.heat` as a "dissolved variable", then compute temperature from the heat density (concentration).
 
 This module also takes care of the water and heat balance connected to direct precipitation inputs to the lake surface and evaporation.
 
+The particular implementation in that module is based on the [GOTM model](https://github.com/gotm-model) and MyLake. We will not explain it further here.
+
 ## Shortwave radiation
 
 Unlike other heat fluxes, shortwave radiation can be passed down through multiple layers. The AirSea module is loaded in a way that specifies the surface incoming shortwave radiation to be passed to the top layer of the lake along the `sw_vert` connection. This is because the location target for the shortwave is specified as `loc(layer.water.heat[sw_vert.top])` (see the load arguments for AirSea in the model).
@@ -145,7 +158,7 @@ var(layer.water.attn, []) {
 }
 ```
 
-The ligh attenuation uses a constant light extinction coefficient `att_c` (parameter). When we add biochemistry, we can make it depend on particle and DOC density in the lake.
+The ligh attenuation uses a constant light extinction coefficient `att_c` (parameter). When we add biochemistry, we will make it depend on particle and DOC density in the lake.
 
 ## Water balance and discharge
 

guide/07/lake_module.txt

@@ -88,13 +88,12 @@ module("Simple layered lake", version(0, 0, 1),
 	} @no_store
 
 	var(layer.water.K, [m 2, day-1]) {
-		# TODO: Is the parametrization by surface area any good?
 		a := (lake.area->[k m 2])=>[],
 		K_ref := {
 			K0*a^0.56   if !lake.ice.ind,
 			K0_ice      otherwise
 		},
-		K_ref*(N2freq=>[])^alpha
+		K_ref*(N2freq=>[])^(2*alpha)
 	}
 
 	flux(layer.water, vert, [m 3, day-1], "Layer mixing down") { 

models/modules/nivafjord/basin.txt

@@ -75,7 +75,7 @@ Authors: Magnus D. Norling
 
 	par_group("Mixing parameters", basin) {	
 		N0     : par_real("Brunt-Väisälä frequency reference", [s-1], 0.8e-2)
-		Nmin   : par_real("Minimum B-V frequency", [s-1], 1e-4, 1e-10, 0.1)
+		Nmin   : par_real("Minimum B-V frequency", [s-1], 0.7e-2, 1e-10, 0.1)
 		alpha  : par_real("Mixing non-linear coefficient", [], 1.4)
 
 		# TODO: Find correct defaults:

models/modules/nivafjord/lake_basin.txt

@@ -32,7 +32,7 @@ This is a module that can be used if a NIVAFjord basin is to be treated as a lak
 
 	flux(layer.water[vert.top], out, [m 3, s-1], "Lake outlet discharge") {
 		excess := max(0, dz[vert.top] - dz0[vert.top]),
-		r1 * excess
+		rl * excess
 	}
 
 }

models/nivafjord_isolated_basin_fpv_model.txt

@@ -10,11 +10,11 @@ It is also currently without horizontal exchange between basins, so it only work
 	# TODO: FPV cover does not dynamically affect wind mixing (not sure to how large a degree that would happen).
 	extend("nivafjord_lake_model.txt")
 
-	fpv       : compartment("FPV", basin_idx)
+	#fpv       : compartment("FPV", basin_idx)
 
-	load("modules/airsea_fpv.txt",
-		module("AirSea FPV", air, basin, fpv, ice, heat, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, indicator,
-			evap, cos_z, loc(basin.freeze_t), loc(basin.area), loc(layer.water[vert.top]), loc(layer.water.heat[sw_vert.top])))
+	#load("modules/airsea_fpv.txt",
+	#	module("AirSea FPV", air, basin, fpv, ice, heat, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, indicator,
+	#		evap, cos_z, loc(basin.freeze_t), loc(basin.area), loc(layer.water[vert.top]), loc(layer.water.heat[sw_vert.top])))
 
 	module("Basin inflow", version(0, 0, 0)) {
 

models/nivafjord_lake_model.txt

@@ -56,10 +56,10 @@ It is almost possible to unify a lot of this one with nivafjord_model.txt . The
 	load("modules/atmospheric.txt",
 		module("Atmospheric", air, temp, wind, g_rad, pressure, a_hum, rho, lwd, cos_z, a_vap, s_vap))
 
-	#fpv : compartment("Floating photovoltaics", basin_idx)
-	#load("modules/airsea_fpv.txt",
-	#	module("AirSea FPV", air, basin, fpv, ice, heat, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, indicator,
-	#		evap, cos_z, loc(basin.freeze_t), loc(basin.area), loc(layer.water[vert.top]), loc(layer.water.heat[sw_vert.top])))
+	fpv : compartment("Floating photovoltaics", basin_idx)
+	load("modules/airsea_fpv.txt",
+		module("AirSea FPV", air, basin, fpv, ice, heat, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, indicator,
+			evap, cos_z, loc(basin.freeze_t), loc(basin.area), loc(layer.water[vert.top]), loc(layer.water.heat[sw_vert.top])))
 
 	load("modules/nivafjord/basin.txt",
 		dims : preamble("NIVAFjord dimensions", basin_idx, layer_idx, layer),

models/nivafjord_lake_simplytox_model.txt

@@ -6,9 +6,6 @@ TODO: This one is incomplete.
 
 Maaybe we want a simpler biochemistry model for this one, not the entire NIVAFjord biochemistry...
 """
-
-
-	extend("nivafjord_lake_model.txt")
 
 	extend("simplytox_model.txt") @exclude(
 		downstream : connection,
@@ -22,6 +19,8 @@ Maaybe we want a simpler biochemistry model for this one, not the entire NIVAFjo
 		wind : property,
 		sol : solver
 	)
+
+	extend("nivafjord_lake_model.txt")
 
 	downstream : connection("Downstream") @directed_graph { river+ } @no_cycles
 

stdlib/seawater.txt

@@ -46,7 +46,7 @@ The implementation used here is influenced by the implementation in [GOTM](https
 
 library("Seawater") {
 """
-This is a library for highly accurate, but more computationally expensive, properties of sea water (taking salinity into account). Se the library "Water utils" for simplified versions of these functions.
+This is a library for highly accurate, but more computationally expensive, properties of sea water (taking salinity into account). Se the library [Water utils](https://nivanorge.github.io/Mobius2/existingmodels/autogen/stdlib.html#water-utils) for simplified freshwater versions of these functions.
 
 The formulas for density are taken from the [Matlab seawater package](http://www.marine.csiro.au/~morgan/seawater).
 
@@ -190,15 +190,14 @@ The implementation is influenced by the one in [SELMA](https://github.com/fabm-m
 
 	o2_saturation : function (T : [deg_c], S : []) {
 		# Formula from Weiss 1970 (also in Selma).
-		T_k := (T -> [K]) => [], # 
-		exp(
-			- 135.90205
+		T_k := (T -> [K]) => [],
+		logsat := - 135.90205
 			+ 1.575701e5  / T_k
 			- 6.642308e7  / T_k^2
 			+ 1.243800e10 / T_k^3
 			- 8.621949e11 / T_k^4
-			- S*(0.017674 - 10.754/T_k + 2140.7/T_k^2)
-		) => [m mol, m-3]
+			- S*(0.017674 - 10.754/T_k + 2140.7/T_k^2),
+		exp(logsat) => [m mol, m-3]
 	}
 
 	o2_piston_velocity : function(wind : [m, s-1], temp : [deg_c]) {