Fix reference in Xiyuan's ocean chemistry example 3

- Add Sarmiento & Gruber (2006) to docs/src/paleo_references.bib
  (this is needed so that the documentation builds)
- Reuse config_ex2.yaml and modify parameters
- Add more explicit refs for choice of modern steady-state values
(- also add an unrelated bug fix for ReactionReservoirAtm for consistency)

docs/src/paleo_references.bib

@@ -14,6 +14,15 @@ @article{Clarkson2018
 	pages = {2918--2923},
 }
 
+@book{SarmientoGruber2006,
+	address = {Princeton},
+	title = {Ocean biogeochemical dynamics},
+	isbn = {978-0-691-01707-5},
+	publisher = {Princeton University Press},
+	author = {Sarmiento, Jorge Louis and Gruber, Nicolas},
+	year = {2006},
+	note = {OCLC: ocm60651167},
+}
 
 @book{Zeebe2001,
 	title = {{CO2} in {Seawater}: {Equilibrium}, {Kinetics}, {Isotopes}},

examples/ocean_chemistry/README.md

@@ -174,18 +174,10 @@ show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL
 
 This example shows how ocean(TAlk-DIC)-air(CO2) exchange in modern state [Example 3 Minimal modern earth ocean-air](@ref).
 
-This configuration is the same as `Example 2 Air-sea exchange`. What is different is that we set `TAlk` and `K_0` to be modern state value [Sarmiento2006](@cite) at S=35, T=25℃.
+This configuration .yaml file is the same as `Example 2 Air-sea exchange`. What is different is that we set `TAlk`, `DIC` and `K_0` to be modern state value (from [SarmientoGruber2006](@cite)) using [`PB.set_parameter_value!`](https://paleotoolkit.github.io/PALEOboxes.jl/stable/Solver%20API/#PALEOboxes.set_parameter_value!) and [`PB.set_variable_attribute!`](https://paleotoolkit.github.io/PALEOboxes.jl/stable/Solver%20API/#PALEOboxes.set_variable_attribute!)
 
-### yaml configuration file
-The model configuration (file `examples/ocean_chemistry/config_ex3.yaml`) contains three Reservoirs `DIC`, `TAlk` and `CO2`. Following `reservoirs` [Example 4 Transfer between Domains](@ref), we use `ReactionFluxTarget` and `ReactionFluxTransfer` to transfer `CO2_airsea_exchange` between `DIC` reservoir and `CO2` reservoir:
-```@eval
-str = read("../../../../examples/ocean_chemistry/config_ex3.yaml", String)
-str = """```julia
-      $str
-      ```"""
-import Markdown
-Markdown.parse(str)
-```
+NB: this is not an accurate model for the modern system ! See [PALEOocean.jl](https://github.com/PALEOtoolkit/PALEOocean.jl) for more complete models,
+including representations of ocean spatial structure and circulation, and more detailed parameterisations of carbonate chemistry ([ReactionCO2SYS](https://paleotoolkit.github.io/PALEOaqchem.jl/dev/PALEOaqchem%20Reactions/#PALEOaqchem.CarbChem.ReactionCO2SYS)) and air-sea exchange ([ReactionAirSeaCO2](https://paleotoolkit.github.io/PALEOocean.jl/dev/PALEOocean_Reactions/#PALEOocean.Oceansurface.AirSeaExchange.ReactionAirSeaCO2))
 
 ### Run script
 The script to run the model (file `examples/ocean_chemistry/run_ex3.jl`) contains:
@@ -197,7 +189,7 @@ str = """```julia
 import Markdown
 Markdown.parse(str)
 ```
-and produces output showing the change, if `TAlk_conc` increase, how the carbonic acid and pH change in ocean and CO2 change in the air:
+and produces output showing an approximate modern steady state:
 ```@example ex3
 include("../../../../examples/ocean_chemistry/run_ex3.jl") # hide
 plot(paleorun.output, ["ocean.TAlk_conc", "ocean.DIC_conc"],                                           (cell=1,); ylabel="TAlk, DIC conc (mol m-3)") # hide
@@ -218,12 +210,5 @@ savefig("ex3_plot5.svg"); nothing  # hide
 ![](ex3_plot4.svg)
 ![](ex3_plot5.svg)
 
-### Displaying model structure and variables
-
-All metadata for model variables can be displayed with `PB.show_variables`:
-```@example ex3
-show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL
-# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code
-```
 
 For more information and cooperation, please communicate with us!

examples/ocean_chemistry/reactions_ReservoirAtm.jl

@@ -51,12 +51,30 @@ function PB.register_methods!(rj::ReactionReservoirAtm)
     @info "register_methods! ReactionReservoirAtm $(PB.fullname(rj)) field_data=$(rj.pars.field_data[])"
     PB.setfrozen!(rj.pars.field_data)
 
+    # callback function to store Variable norm during setup
+    function setup_callback(m, attribute_value, v, vdata)
+        v.localname == "R" || error("setup_callback unexpected Variable $(PB.fullname(v))")
+        if attribute_value == :norm_value
+            m.reaction.norm_value = PB.value_ad(PB.get_total(vdata[]))
+        end
+        return nothing
+    end
+
+
     if rj.pars.const[]
         R = PB.VarPropScalar(      "R", "mol", "scalar constant reservoir", attributes=(:field_data=>rj.pars.field_data[],))
         R_sms = PB.VarTarget(     "R_sms", "mol yr-1", "constant reservoir source-sinks", attributes=(:field_data=>rj.pars.field_data[],))
+        PB.add_method_setup_initialvalue_vars_default!(
+            rj, 
+            [R],
+            filterfn = v->true,  # set filterfn to force setup even if R is constant, not a state Variable
+            force_initial_norm_value=true, # setup :norm_value, :initial_value to get norm_value callback, even though R is not a state Variable
+            setup_callback=setup_callback
+        )  
     else
         R = PB.VarStateExplicitScalar("R", "mol", "scalar reservoir", attributes=(:field_data=>rj.pars.field_data[],))
         R_sms = PB.VarDerivScalar(     "R_sms", "mol yr-1", "scalar reservoir source-sinks", attributes=(:field_data=>rj.pars.field_data[],))
+        PB.add_method_setup_initialvalue_vars_default!(rj, [R], setup_callback=setup_callback)
     end
     PB.setfrozen!(rj.pars.const)
 
@@ -74,16 +92,7 @@ function PB.register_methods!(rj::ReactionReservoirAtm)
             PB.VarPropScalar(     "R_delta", "per mil", "scalar atmosphere reservoir isotope delta"))
     end
 
-    # callback function to store Variable norm during setup
-    function setup_callback(m, attribute_value, v, vdata)
-        v.localname == "R" || error("setup_callback unexpected Variable $(PB.fullname(v))")
-        if attribute_value == :norm_value
-            m.reaction.norm_value = PB.value_ad(PB.get_total(vdata[]))
-        end
-        return nothing
-    end
-    # set filterfn to force setup even if R is constant, not a state Variable
-    PB.add_method_setup_initialvalue_vars_default!(rj, [R], filterfn = v->true, setup_callback=setup_callback)  
+
 
     PB.add_method_do!(
         rj,

examples/ocean_chemistry/run_ex3.jl

@@ -11,7 +11,18 @@ include(joinpath(@__DIR__,"reactions_ReservoirAtm.jl")) # @__DIR__ so still runs
 #####################################################
 # Create model
 
-model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex3.yaml"), "Minimal_modern_AirSea")
+model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex2.yaml"), "Minimal_Alk_pH_AirSea")
+
+# set to ~modern (1990s) values from Sarmiento 2006 for global mean ocean surface
+
+# K0 (NB: Table A.3 gives global mean ocean surface temperature  17.88 C)
+PB.set_parameter_value!(model, "oceansurface", "solve_AirSea_Exchange", "K_0", 38.3) # mol m-3 atm-1, Table 3.2.3 at S=35. T=15℃
+# PB.set_parameter_value!(model, "oceansurface", "solve_AirSea_Exchange", "K_0", 33.11) # mol m-3 atm-1, Table 3.2.3 at S=35. T=20℃
+# PB.set_parameter_value!(model, "oceansurface", "solve_AirSea_Exchange", "K_0", 29.6) # mol m-3 atm-1, from Sarmiento 2006, at S=35. T=25℃
+
+PB.set_parameter_value!(model, "global", "add_Alk", "perturb_totals", [0.0, 0.0]) # don't add any Alk
+PB.set_variable_attribute!(model, "ocean.TAlk", :initial_value, 2.308)  # mol m-3 ~ modern value global mean ocean surface (Table A.3)
+PB.set_variable_attribute!(model, "ocean.DIC", :initial_value, 2.026)  # mol m-3 ~ modern value global mean ocean surface  (Table A.3)
 
 #########################################################
 # Initialize