diff --git a/docs/src/paleo_references.bib b/docs/src/paleo_references.bib index 835041b..081397c 100644 --- a/docs/src/paleo_references.bib +++ b/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}}, diff --git a/examples/ocean_chemistry/README.md b/examples/ocean_chemistry/README.md index 254beeb..381b2aa 100644 --- a/examples/ocean_chemistry/README.md +++ b/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! \ No newline at end of file diff --git a/examples/ocean_chemistry/config_ex3.yaml b/examples/ocean_chemistry/config_ex3.yaml deleted file mode 100644 index 813f198..0000000 --- a/examples/ocean_chemistry/config_ex3.yaml +++ /dev/null @@ -1,207 +0,0 @@ -Minimal_modern_AirSea: - - domains: - - - global: - - reactions: - - sum_C: - - class: ReactionSum - - parameters: - - vars_to_add: ["ocean.DIC_total", "atm.CO2"] - - variable_links: - - sum: C_total - - add_Alk: - # from PALEOboxes reaction catalog: constant input flux of Alk - class: ReactionFluxPerturb - parameters: - # linear interpolation for input flux - set to constant - perturb_times: [-1e30, 1e30] - # ocean volume = 3.6e14 * 3688.0 = 1.33e18 m^3 - # DIC in ocean = 2 * 1.33e18 = 2.66e18 mol - # set Alk flux so Alk = DIC after 100 yr-1 = 2.66e18/100 = 2.66e16 mol yr-1 - perturb_totals: [0.0, 0.0] - variable_links: - F: ocean.TAlk_sms # add flux to TAlk reservoir NB: works for single-box ocean only !! - - - - # -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - - ocean: - - reactions: - - oceantrasport: - - class: ReactionOceanNoTransport - - parameters: - - area: [3.6e14] - depth: [3688.0] - - - - reservoir_TA: - - class: ReactionReservoirTotal - - variable_links: - R*: TAlk* - variable_attributes: - R:initial_value: 2.4e0 - R:norm_value: 2.4e0 - - - - reservoir_DIC: - class: ReactionReservoirTotal - - variable_links: - R*: DIC* - variable_attributes: - R:initial_value: 2.0e0 - R:norm_value: 2.0e0 - - constant_B: - class: ReactionReservoirConst - variable_links: - R*: B - variable_attributes: - R_conc:initial_value: 0.427 # mol m-3 contemporary value - - solve_Alk_pH: - - class: Reaction_Alk_pH - - parameters: - - K_1: 1.4e-6 - K_2: 1.2e-9 - K_w: 6.0e-14 - K_B: 2.5e-9 - - variable_links: - density: rho_ref # OceanNoTransport provides density as rho_ref - - -# -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - - oceansurface: - - # air-sea exchange reactions - - reactions: - # calculate air-sea exchange flux - solve_AirSea_Exchange: - - class: Reaction_Min_AirSeaExchange - - parameters: - - K_0: 29.6 # mol m-3 atm-1, from Sarmiento 2006, at S=35. T=25℃, -log10(K0)=1.54 - vpiston: 1138.8 # m yr-1 - - variable_links: - - pXatm: atm.pCO2atm - X_aq_conc: ocean.oceansurface.CO2_aq_conc - - area: Asurf # ocean surface area - - X_airsea_exchange: fluxAtmtoOceansurface.flux_CO2 - - - # apply air-sea flux to ocean surface - transfer_fluxAtmtoOceansurface: - - class: ReactionFluxTransfer - - parameters: - - transfer_matrix: Identity - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: ocean.oceansurface.$fluxname$_sms - - variable_links: - - output_CO2: ocean.oceansurface.DIC_sms - -# -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - - oceanfloor: - -# -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - - atm: - - - reactions: - - reservoir_CO2_atm: - - class: ReactionReservoirAtm - - parameters: - - moles1atm: 1.77e20 - - - variable_links: - - R*: CO2* - pRatm: pCO2atm - pRnorm: pCO2PAL - - variable_attributes: - - R:norm_value: 4.956e16 # 280e-6 ppm - R:initial_value: 4.956e16 - - transfer_AtmtoOceansurface: - - class: ReactionFluxTransfer - - parameters: - - transfer_matrix: Distribute - transfer_multiplier: -1.0 - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: $fluxname$_sms - - variable_links: - - output_CO2: atm.CO2_sms - - - - - -# -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - - fluxAtmtoOceansurface: - - reactions: - - fluxtarget: - - class: ReactionFluxTarget - - parameters: - - flux_totals: true - - fluxlist: ["CO2"] - - - - diff --git a/examples/ocean_chemistry/reactions_ReservoirAtm.jl b/examples/ocean_chemistry/reactions_ReservoirAtm.jl index c92c4c0..14d3bc6 100644 --- a/examples/ocean_chemistry/reactions_ReservoirAtm.jl +++ b/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, diff --git a/examples/ocean_chemistry/run_ex3.jl b/examples/ocean_chemistry/run_ex3.jl index f3eee29..9cc3e80 100644 --- a/examples/ocean_chemistry/run_ex3.jl +++ b/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