From eb9129ebd361a6fbe1a10e60d6200ea752c6ee27 Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Tue, 18 Jun 2024 11:54:12 +0100 Subject: [PATCH] Add ICBM global ocean and Black Sea examples from Romeniello & Derry (2010) examples/romglb/ examples/blacksea/ See README.md for details including how to download data files. These use 'ReactionOceanTransportRomanielloICBM' which reads the Matlab .mat files containing transport matrix and geometry for the ICBM models of global ocean and Black sea from: - Romaniello & Derry (2010) 'An intermediate-complexity model for simulating marine biogeochemistry in deep time: Validation against the modern global ocean', GGG, - Romaniello & Derry (2010) 'Validation of an intermediate-complexity model for simulating marine biogeochemistry under anoxic conditions in the modern Black Sea', GGG, --- .gitignore | 1 + docs/src/PALEOocean_Reactions.md | 3 +- docs/src/paleo_references.bib | 23 + .../PALEO_examples_blacksea_O2_only.jl | 58 + .../blacksea/PALEO_examples_blacksea_P_O2.jl | 59 + .../PALEO_examples_blacksea_P_O2_SO4.jl | 58 + .../blacksea/PALEO_examples_blacksea_cfg.yaml | 543 +++++++++ examples/blacksea/README.md | 41 + .../blacksea/config_ocean_blacksea_expts.jl | 23 + examples/blacksea/plot_ocean_blacksea.jl | 92 ++ examples/blacksea/runtests.jl | 112 ++ examples/doc_order.txt | 2 + .../romglb/PALEO_examples_romglb_O2_only.jl | 60 + examples/romglb/PALEO_examples_romglb_P_O2.jl | 60 + .../PALEO_examples_romglb_P_O2_S_Carb_open.jl | 77 ++ .../romglb/PALEO_examples_romglb_P_O2_open.jl | 64 + .../romglb/PALEO_examples_romglb_cfg.yaml | 1045 +++++++++++++++++ examples/romglb/README.md | 56 + examples/romglb/SedimentationRate_dev.jl | 108 ++ examples/romglb/Test xarray netcdf.ipynb | 224 ++++ examples/romglb/config_ocean_romglb_expts.jl | 26 + examples/romglb/plot_ocean_romglb.jl | 185 +++ examples/romglb/runtests.jl | 212 ++++ src/ocean/Ocean.jl | 2 + src/ocean/OceanTransportRomaniello2010.jl | 489 ++++++++ 25 files changed, 3622 insertions(+), 1 deletion(-) create mode 100644 examples/blacksea/PALEO_examples_blacksea_O2_only.jl create mode 100644 examples/blacksea/PALEO_examples_blacksea_P_O2.jl create mode 100644 examples/blacksea/PALEO_examples_blacksea_P_O2_SO4.jl create mode 100644 examples/blacksea/PALEO_examples_blacksea_cfg.yaml create mode 100644 examples/blacksea/README.md create mode 100644 examples/blacksea/config_ocean_blacksea_expts.jl create mode 100644 examples/blacksea/plot_ocean_blacksea.jl create mode 100644 examples/blacksea/runtests.jl create mode 100644 examples/romglb/PALEO_examples_romglb_O2_only.jl create mode 100644 examples/romglb/PALEO_examples_romglb_P_O2.jl create mode 100644 examples/romglb/PALEO_examples_romglb_P_O2_S_Carb_open.jl create mode 100644 examples/romglb/PALEO_examples_romglb_P_O2_open.jl create mode 100644 examples/romglb/PALEO_examples_romglb_cfg.yaml create mode 100644 examples/romglb/README.md create mode 100644 examples/romglb/SedimentationRate_dev.jl create mode 100644 examples/romglb/Test xarray netcdf.ipynb create mode 100644 examples/romglb/config_ocean_romglb_expts.jl create mode 100644 examples/romglb/plot_ocean_romglb.jl create mode 100644 examples/romglb/runtests.jl create mode 100644 src/ocean/OceanTransportRomaniello2010.jl diff --git a/.gitignore b/.gitignore index 29e748e..eb7dec6 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,4 @@ LocalPreferences.toml *.nc *.zip examples/shelf1D/S2P3_transport_20240614 +examples/romglb/romaniello2010_transport diff --git a/docs/src/PALEOocean_Reactions.md b/docs/src/PALEOocean_Reactions.md index 52435a5..6381f4f 100644 --- a/docs/src/PALEOocean_Reactions.md +++ b/docs/src/PALEOocean_Reactions.md @@ -9,8 +9,9 @@ CurrentModule = PALEOocean.Ocean OceanNoTransport.ReactionOceanNoTransport OceanTransport3box.ReactionOceanTransport3box OceanTransport6box.ReactionOceanTransport6box -OceanTransportTMM.ReactionOceanTransportTMM OceanTransportColumn.ReactionOceanTransportColumn +OceanTransportRomaniello2010.ReactionOceanTransportRomanielloICBM +OceanTransportTMM.ReactionOceanTransportTMM ``` ## Vertical Transport diff --git a/docs/src/paleo_references.bib b/docs/src/paleo_references.bib index 428c946..7ee1616 100644 --- a/docs/src/paleo_references.bib +++ b/docs/src/paleo_references.bib @@ -151,6 +151,29 @@ @article{Ridgwell2007 pages = {87--104}, } +@article{Romaniello2010a, +author = {Romaniello, Stephen J. and Derry, Louis A}, +doi = {10.1029/2009GC002712}, +journal = {Geochemistry Geophysics Geosystems}, +month = {aug}, +number = {8}, +pages = {1--18}, +title = {{Validation of an intermediate-complexity model for simulating marine biogeochemistry under anoxic conditions in the modern Black Sea}}, +volume = {11}, +year = {2010} +} + +@article{Romaniello2010, +author = {Romaniello, Stephen J. and Derry, Louis A}, +doi = {10.1029/2009GC002711}, +journal = {Geochemistry Geophysics Geosystems}, +month = {aug}, +number = {8}, +pages = {1--33}, +title = {{An intermediate-complexity model for simulating marine biogeochemistry in deep time: Validation against the modern global ocean}}, +volume = {11}, +year = {2010} +} @article{Sarmiento1984, author = {Sarmiento, Jorge L and Toggweiler, J. R.}, diff --git a/examples/blacksea/PALEO_examples_blacksea_O2_only.jl b/examples/blacksea/PALEO_examples_blacksea_O2_only.jl new file mode 100644 index 0000000..d7274b5 --- /dev/null +++ b/examples/blacksea/PALEO_examples_blacksea_O2_only.jl @@ -0,0 +1,58 @@ +using Logging +using DiffEqBase +using Sundials + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOocean + + +global_logger(ConsoleLogger(stderr,Logging.Info)) + +include("config_ocean_blacksea_expts.jl") +include("plot_ocean_blacksea.jl") + +# abiotic atmosphere-ocean, O2 only +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_blacksea_cfg.yaml"), "blacksea_abiotic_O2", + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "../romglb/romaniello2010_transport") # assume zip file has been downloaded and unpacked in this subfolder + ) +) + +config_ocean_blacksea_expts(model, ["baseline"]); tspan=(0, 1e4) + +initial_state, modeldata = PALEOmodel.initialize!(model) + +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# No carbonate system, so can integrate as an ODE +PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, + solvekwargs=(reltol=1e-5,), +) + + +######################################## +# Plot output +######################################## + +# individual plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# assemble plots onto screens with 6 subplots +gr(size=(1200, 900)) +pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm))) + +plot_totals(paleorun.output; species=["T", "O2"], pager=pager) +plot_airsea(paleorun.output; species=["O2"], pager=pager) +plot_ocean_tracers( + paleorun.output; + tracers=["insol", "T_conc", "O2_conc"], + tcol=[-Inf, 10.0, 100.0, 1000.0, 1e4, 1e5], + pager=pager +) +pager(:newpage) # flush output diff --git a/examples/blacksea/PALEO_examples_blacksea_P_O2.jl b/examples/blacksea/PALEO_examples_blacksea_P_O2.jl new file mode 100644 index 0000000..d072563 --- /dev/null +++ b/examples/blacksea/PALEO_examples_blacksea_P_O2.jl @@ -0,0 +1,59 @@ +using Logging +using DiffEqBase +using Sundials + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOocean + + +global_logger(ConsoleLogger(stderr,Logging.Info)) + +include("config_ocean_blacksea_expts.jl") +include("plot_ocean_blacksea.jl") + +# abiotic atmosphere-ocean, O2 only +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_blacksea_cfg.yaml"), "blacksea_P_O2", + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "../romglb/romaniello2010_transport") # assume zip file has been downloaded and unpacked in this subfolder + ) +) + +config_ocean_blacksea_expts(model, ["baseline"]); tspan=(0, 1e5) + +initial_state, modeldata = PALEOmodel.initialize!(model) + +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# No carbonate system, so can integrate as an ODE +PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, + solvekwargs=(reltol=1e-5,), +) + + +######################################## +# Plot output +######################################## + +# individual plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# assemble plots onto screens with 6 subplots +gr(size=(1200, 900)) +pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm))) + +plot_totals(paleorun.output; species=["T", "O2", "P"], pager=pager) +plot_airsea(paleorun.output; species=["O2"], pager=pager) +pager(:newpage) +plot_ocean_tracers( + paleorun.output; + tracers=["insol", "T_conc", "O2_conc", "P_conc"], + tcol=[-Inf, 10.0, 100.0, 1000.0, 1e4, 1e5], + pager=pager +) +pager(:newpage) # flush output diff --git a/examples/blacksea/PALEO_examples_blacksea_P_O2_SO4.jl b/examples/blacksea/PALEO_examples_blacksea_P_O2_SO4.jl new file mode 100644 index 0000000..968cc08 --- /dev/null +++ b/examples/blacksea/PALEO_examples_blacksea_P_O2_SO4.jl @@ -0,0 +1,58 @@ +using Logging +using DiffEqBase +using Sundials + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOocean + + +global_logger(ConsoleLogger(stderr,Logging.Info)) + +include("config_ocean_blacksea_expts.jl") +include("plot_ocean_blacksea.jl") + +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_blacksea_cfg.yaml"), "blacksea_P_O2_SO4", + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "../romglb/romaniello2010_transport") # assume zip file has been downloaded and unpacked in this subfolder + ) +) + +config_ocean_blacksea_expts(model, ["baseline"]); tspan=(0, 1e5) + +initial_state, modeldata = PALEOmodel.initialize!(model) + +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# No carbonate system, so can integrate as an ODE +PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, + solvekwargs=(reltol=1e-5,), +) + + +######################################## +# Plot output +######################################## + +# individual plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# assemble plots onto screens with 6 subplots +gr(size=(1200, 900)) +pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm))) + +plot_totals(paleorun.output; species=["T", "O2", "P", "S"], pager=pager) +plot_airsea(paleorun.output; species=["O2"], pager=pager) +pager(:newpage) +plot_ocean_tracers( + paleorun.output; + tracers=["insol", "T_conc", "O2_conc", "P_conc", "SO4_conc", "H2S_conc"], + tcol=[-Inf, 10.0, 100.0, 1000.0, 1e4, 1e5], + pager=pager +) +pager(:newpage) # flush output diff --git a/examples/blacksea/PALEO_examples_blacksea_cfg.yaml b/examples/blacksea/PALEO_examples_blacksea_cfg.yaml new file mode 100644 index 0000000..dddd5ba --- /dev/null +++ b/examples/blacksea/PALEO_examples_blacksea_cfg.yaml @@ -0,0 +1,543 @@ +#################################### +# Abiotic atmosphere/ocean, O2 only +# Bosphorus outflow is added back to oceansurface to form a closed system +#################################### +blacksea_abiotic_O2: + parameters: + # CIsotope: ScalarData + matdir: ../romglb/romaniello2010_transport # folder with data files from Romaniello (2010) GGG + domains: + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["O2"] + + fluxBosphorusOutflow: + # scalar domain + reactions: + outflowtarget: + class: ReactionFluxTarget + parameters: + fluxlist: ["T", "O2"] + + global: + # scalar domain + + reactions: + + atm: + + + reactions: + constant_pO2: + class: ReactionConst + + parameters: + constnames: ["pO2atm"] + variable_attributes: + pO2atm:initial_value: [0.21] # mol/mol + + oceansurface: + reactions: + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston: 4.8 # m d-1 + + surface_insol: + class: ReactionConst + parameters: + constnames: [surface_insol] + variable_attributes: + surface_insol:initial_value: [350.0] # W m-2 downwelling radiative flux + + transfer_AtmtoOceansurface: + 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 + + transfer_BosphorusOutflow: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + input_fluxes: fluxBosphorusOutflow.flux_$fluxname$ + output_fluxes: ocean.oceansurface.$fluxname$_sms + + ocean: + reactions: + transportromblacksea: + class: ReactionOceanTransportRomanielloICBM + parameters: + matdir: external%matdir + circname: Black_Sea + bosph_outflow: 6.04e11 # m^3 yr^-1 + bosph_inflow: 3.05e11 + variable_links: + # the easy way to create a closed system: add back Bosphorus outflow to surface ocean + # bosph_outflow_*: ocean.oceansurface.*_sms + # source Bosphorus inflow from surface ocean to create closed system + bosph_inflow_*_conc: ocean.oceansurface.*_conc + bosph_inflow_*_sms: ocean.oceansurface.*_sms + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_T: # passive test tracer to check conservation + class: ReactionReservoirTotal + variable_links: + R*: T* + variable_attributes: + R:initial_value: 1.0 # concentration m-3 + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.04 + + + oceanfloor: + reactions: + +################################# +# Atmosphere/ocean, O2 and P only +# Closed system (outflow -> surface box, and surface box -> inflow) +################################# + +blacksea_P_O2: + parameters: + CIsotope: ScalarData + SIsotope: ScalarData + matdir: ../romglb/romaniello2010_transport # folder with data files from Romaniello (2010) GGG + domains: + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["O2"] + + fluxOceanfloor: + reactions: + particulatetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: particulateflux_ + fluxlist: ["Corg::CIsotope", "N", "P", "Ccarb::CIsotope"] + + solutetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: soluteflux_ + fluxlist: ["DIC::CIsotope", "TAlk", "O2", "P"] + + global: + # scalar domain + + reactions: + total_O2: + class: ReactionSum + parameters: + vars_to_add: [ocean.O2] + variable_links: + sum: total_O2 + atm: + + + reactions: + constant_pO2: + class: ReactionConst + + parameters: + constnames: ["pO2atm"] + variable_attributes: + pO2atm:initial_value: [0.21] # mol/mol + + oceansurface: + reactions: + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston: 4.8 # m d-1 + + surface_insol: + class: ReactionConst + parameters: + constnames: [surface_insol] + variable_attributes: + surface_insol:initial_value: [350.0] # W m-2 downwelling radiative flux + + transfer_AtmtoOceansurface: + 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 + + ocean: + reactions: + transportromglb: + class: ReactionOceanTransportRomanielloICBM + parameters: + matdir: external%matdir + circname: Black_Sea + bosph_outflow: 6.04e11 # m^3 yr^-1 + bosph_inflow: 3.05e11 + variable_links: + # add back Bosphorus outflow to surface ocean to create a closed system + bosph_outflow_*: ocean.oceansurface.*_sms + # source Bosphorus inflow from surface ocean to create closed system + bosph_inflow_*_conc: ocean.oceansurface.*_conc + bosph_inflow_*_sms: ocean.oceansurface.*_sms + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_P: + class: ReactionReservoirTotal + variable_links: + R*: P* + variable_attributes: + R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) + + reservoir_T: # passive test tracer to check conservation + class: ReactionReservoirTotal + variable_links: + R*: T* + variable_attributes: + R:initial_value: 1.0 # concentration m-3 + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.0 # zero opacity as we just want to prescribe insol in surface boxes + + bioprod: + class: ReactionBioProdMMPop + parameters: + depthlimit: -10.0 # surface cells only + + rCorgPO4: 106.0 + rNPO4: 16.0 + rCcarbCorg: 0.25 + rCcarbCorg_fixed: true + + nuDOM: 0.0 # no DOM pool + + k_poptype: Constant + k_uPO4: 2.615e-3 # 1.91e-6 mol / kg-sw / yr * 1027 kg m-3 * 4/3 + + k_nuttype: PO4MM + k_KPO4: 0.21567e-3 # 0.21e-6 mol / kg-sw * 1027 kg m-3 + + k_lightlim: linear + k_Irel: 1.0 + + variable_links: + partprod_*: export_* + + biopump: + class: ReactionExportDirectColumn + + parameters: + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] + transportfloor: true + exportfunction: SumExp + input_frac: [0.9354, 0.0646] # 2-G model of Ridgwell & Hargreaves (2007) + sumexp_scale: [550.5195, 1e6] + + reminocean: + class: ReactionReminO2 + + parameters: + + variable_links: + soluteflux_*: "*_sms" + + + oceanfloor: + reactions: + + reminoceanfloor: + class: ReactionReminO2 + + parameters: + + variable_links: + remin*: particulateflux* + + soluteflux_*: fluxOceanfloor.soluteflux_* + + transferparticulate_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: particulateflux_$fluxname$ + + transfersolute_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ + output_fluxes: ocean.oceanfloor.$fluxname$_sms + + sedcrust: + + +################################# +# Atmosphere/ocean, O2, SO4/H2S and P only +################################# + +blacksea_P_O2_SO4: + parameters: + CIsotope: ScalarData + SIsotope: ScalarData + matdir: ../romglb/romaniello2010_transport # folder with data files from Romaniello (2010) GGG + domains: + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["O2"] + + fluxOceanfloor: + reactions: + particulatetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: particulateflux_ + fluxlist: ["Corg::CIsotope", "N", "P", "Ccarb::CIsotope"] + + solutetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: soluteflux_ + fluxlist: ["O2", "P", "SO4", "H2S"] + + + + global: + # scalar domain + + reactions: + total_O2: + class: ReactionSum + parameters: + vars_to_add: [ocean.O2_total] + variable_links: + sum: total_O2 + + total_S: + class: ReactionSum + parameters: + vars_to_add: [ocean.SO4_total, ocean.H2S_total] + variable_links: + sum: total_S + + atm: + + + reactions: + constant_pO2: + class: ReactionConst + + parameters: + constnames: ["pO2atm"] + variable_attributes: + pO2atm:initial_value: [0.21] # mol/mol + + oceansurface: + reactions: + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston: 4.8 # m d-1 + + surface_insol: + class: ReactionConst + parameters: + constnames: [surface_insol] + variable_attributes: + surface_insol:initial_value: [350.0] # W m-2 downwelling radiative flux + + transfer_AtmtoOceansurface: + 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 + + ocean: + reactions: + transportromglb: + class: ReactionOceanTransportRomanielloICBM + parameters: + matdir: external%matdir + circname: Black_Sea + bosph_outflow: 6.04e11 # m^3 yr^-1 + bosph_inflow: 3.05e11 + variable_links: + # add back Bosphorus outflow to surface ocean to create a closed system + bosph_outflow_*: ocean.oceansurface.*_sms + # source Bosphorus inflow from surface ocean to create closed system + bosph_inflow_*_conc: ocean.oceansurface.*_conc + bosph_inflow_*_sms: ocean.oceansurface.*_sms + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_SO4: + class: ReactionReservoirTotal + parameters: + field_data: external%SIsotope + variable_links: + R*: SO4* + variable_attributes: + R:initial_value: 28756.0e-3 # concentration mol m-3 = 28e-3 mol/kg * 1027 kg m-3 + R:norm_value: 1000.0e-3 # for scaling only + + reservoir_H2S: + class: ReactionReservoirTotal + parameters: + field_data: external%SIsotope + variable_links: + R*: H2S* + variable_attributes: + R:initial_value: 1e-6 # concentration mol m-3 ~ 1e-9 mol/kg * 1027 kg m-3 + R:norm_value: 1.0e-3 # for scaling only + + reservoir_P: + class: ReactionReservoirTotal + variable_links: + R*: P* + variable_attributes: + R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) + + reservoir_T: # passive test tracer to check conservation + class: ReactionReservoirTotal + variable_links: + R*: T* + variable_attributes: + R:initial_value: 1.0 # concentration m-3 + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.0 # zero opacity as we just want to prescribe insol in surface boxes + + bioprod: + class: ReactionBioProdMMPop + parameters: + depthlimit: -10.0 # surface cells only + + rCorgPO4: 106.0 + rNPO4: 16.0 + rCcarbCorg: 0.25 + rCcarbCorg_fixed: true + + nuDOM: 0.0 # no DOM pool + + k_poptype: Constant + k_uPO4: 2.615e-3 # 1.91e-6 mol / kg-sw / yr * 1027 kg m-3 * 4/3 + + k_nuttype: PO4MM + k_KPO4: 0.21567e-3 # 0.21e-6 mol / kg-sw * 1027 kg m-3 + + k_lightlim: linear + k_Irel: 1.0 + + variable_links: + partprod_*: export_* + + biopump: + class: ReactionExportDirectColumn + + parameters: + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] + transportfloor: true + exportfunction: SumExp + input_frac: [0.9354, 0.0646] # 2-G model of Ridgwell & Hargreaves (2007) + sumexp_scale: [550.5195, 1e6] + + reminocean: + class: ReactionReminO2_SO4 + + parameters: + + variable_links: + soluteflux_*: "*_sms" + + redox_H2S_O2: + class: ReactionRedoxH2S_O2 + + parameters: + R_H2S_O2: 3.65e3 # (mol m-3) yr-1 + + + oceanfloor: + reactions: + + reminoceanfloor: + class: ReactionReminO2_SO4 + + parameters: + + variable_links: + remin*: particulateflux* + + O2_conc: ocean.oceanfloor.O2_conc + + soluteflux_*: fluxOceanfloor.soluteflux_* + + transferparticulate_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: particulateflux_$fluxname$ + + transfersolute_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ + output_fluxes: ocean.oceanfloor.$fluxname$_sms + + sedcrust: + + + + diff --git a/examples/blacksea/README.md b/examples/blacksea/README.md new file mode 100644 index 0000000..016e4d1 --- /dev/null +++ b/examples/blacksea/README.md @@ -0,0 +1,41 @@ +# Black Sea intermediate-complexity box model + +Black sea examples and test cases, +using the 1-column model for the Black sea circulation from +[Romaniello2010a](@cite). + +## Installation + +These examples use `ReactionOceanTransportRomanielloICBM` to read Matlab data files with the +Black sea circulation from [Romaniello2010a](@cite). + +The Matlab datafiles are available as a zip file from , +generated from the Matlab model code available as Supplementary Information to [Romaniello2010](@cite). + +The examples assume the zip file has been downloaded and unpacked to subfolder `../romglb/romaniello2010_transport` +(NB: a subfolder of `romglb`, not this folder `blacksea`, as the zip file contains both global and Black sea data files) + +## O2 only air-sea exchange and transport test + + include("PALEO_examples_blacksea_O2_only.jl") + +## P, O2 with parameterized export production + + include("PALEO_examples_blacksea_P_O2.jl") + +Phosphorus and oxygen, with a parameterization of export production based on light and nutrient availability. + +Configured to create a closed circulation by returning Bosphorus outflow back to surface box, +and sourcing Bosphorus inflow from surface box. + +No burial fluxes (ie oceanfloor phosphorus flux is recycled into water column), +so the Black sea is effectively a closed system for phosphorus. + + +## P, O2, S with parameterized export production + + include("PALEO_examples_blacksea_P_O2_SO4.jl") + +Phosphorus, oxygen, sulphur, with a parameterization of +export production based on light and nutrient availability. + diff --git a/examples/blacksea/config_ocean_blacksea_expts.jl b/examples/blacksea/config_ocean_blacksea_expts.jl new file mode 100644 index 0000000..6a0145f --- /dev/null +++ b/examples/blacksea/config_ocean_blacksea_expts.jl @@ -0,0 +1,23 @@ +function config_ocean_blacksea_expts( + model, expts +) + for expt in expts + if expt == "baseline" + # baseline configuration + + elseif length(expt) == 5 && expt[1] == "setpar" + # generic parameter set (setpar, , , , matdir, + ) + ) + + tspan=(0, 1e5) + + initial_state, modeldata = PALEOmodel.initialize!(model) + paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + + PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, solvekwargs=(reltol=1e-9,) + ) + + println("conservation checks:") + conschecks = [ + ("ocean", "T_total", 1e-6 ), + ] + for (domname, varname, rtol) in conschecks + startval, endval = PB.get_data(paleorun.output, domname*"."*varname)[[1, end]] + println(" check $domname.$varname $startval $endval $rtol") + @test isapprox(startval, endval, rtol=rtol) + end + + println("check values at end of run:") + checkvals = [ + ] + for (domname, varname, checkval, rtol) in checkvals + outputval = PB.get_data(paleorun.output, domname*"."*varname)[end] + println(" check $domname.$varname $outputval $checkval $rtol") + @test isapprox(outputval, checkval, rtol=rtol) + end + +end + +!("P_O2" in skipped_testsets) && @testset "P_O2" begin + + model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_blacksea_cfg.yaml"), "blacksea_P_O2", + modelpars=Dict( + "matdir"=>matdir, + ) + ) + + tspan=(0, 1e5) + + initial_state, modeldata = PALEOmodel.initialize!(model) + paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + + PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, solvekwargs=(reltol=1e-9,) + ) + + ncoutfile = tempname(cleanup=true) + PALEOmodel.OutputWriters.save_netcdf(paleorun.output, ncoutfile; check_ext=false) + ncoutput = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), ncoutfile; check_ext=false) + + for output in (paleorun.output, ncoutput) + println("conservation checks:") + conschecks = [ + ("ocean", "T_total", 1e-6 ), + ("ocean", "P_total", 1e-6 ), + ] + for (domname, varname, rtol) in conschecks + startval, endval = PB.get_data(output, domname*"."*varname)[[1, end]] + println(" check $domname.$varname $startval $endval $rtol") + @test isapprox(startval, endval, rtol=rtol) + end + + println("check values at end of run:") + checkvals = [ + ] + for (domname, varname, checkval, rtol) in checkvals + outputval = PB.get_data(output, domname*"."*varname)[end] + println(" check $domname.$varname $outputval $checkval $rtol") + @test isapprox(outputval, checkval, rtol=rtol) + end + end + +end + + +end diff --git a/examples/doc_order.txt b/examples/doc_order.txt index dd86548..8afad92 100644 --- a/examples/doc_order.txt +++ b/examples/doc_order.txt @@ -1,6 +1,8 @@ ocean3box PTBClarkson2014 shelf1D +blacksea +romglb mitgcm skeleton_configuration transport_examples diff --git a/examples/romglb/PALEO_examples_romglb_O2_only.jl b/examples/romglb/PALEO_examples_romglb_O2_only.jl new file mode 100644 index 0000000..ba6cf77 --- /dev/null +++ b/examples/romglb/PALEO_examples_romglb_O2_only.jl @@ -0,0 +1,60 @@ +using Logging +using DiffEqBase +using Sundials + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOocean + + +global_logger(ConsoleLogger(stderr,Logging.Info)) + +include("../atmreservoirreaction.jl") +include("config_ocean_romglb_expts.jl") +include("plot_ocean_romglb.jl") + +# abiotic atmosphere-ocean, O2 only +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_abiotic_O2"; + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "romaniello2010_transport") # assume zip file has been downloaded and unpacked in this subfolder + ) +) + +config_ocean_romglb_expts(model, ["baseline"]); tspan=(0, 1e5) + +initial_state, modeldata = PALEOmodel.initialize!(model) + +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# No carbonate system, so can integrate as an ODE +PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, + solvekwargs=(reltol=1e-5,), +) + + +######################################## +# Plot output +######################################## + +# individual plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# assemble plots onto screens with 6 subplots +gr(size=(1200, 900)) +pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm))) + +plot_totals(paleorun.output; species=["T", "O2"], pager=pager) +plot_airsea(paleorun.output; species=["O2"], pager=pager) +pager(:newpage) +plot_ocean_tracers( + paleorun.output; + tracers=["insol", "T_conc", "O2_conc"], + tcol=[-Inf, 10.0, 100.0, 1000.0, 1e4, 1e5], + pager=pager +) +pager(:newpage) # flush output diff --git a/examples/romglb/PALEO_examples_romglb_P_O2.jl b/examples/romglb/PALEO_examples_romglb_P_O2.jl new file mode 100644 index 0000000..2d12d81 --- /dev/null +++ b/examples/romglb/PALEO_examples_romglb_P_O2.jl @@ -0,0 +1,60 @@ +using Logging +using DiffEqBase +using Sundials + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOocean + + +global_logger(ConsoleLogger(stderr,Logging.Info)) + + +include("../atmreservoirreaction.jl") +include("config_ocean_romglb_expts.jl") +include("plot_ocean_romglb.jl") + +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_P_O2", + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "romaniello2010_transport") # assume zip file has been downloaded and unpacked in this subfolder + ) +) + +config_ocean_romglb_expts(model, ["baseline"]); tspan=(0, 1e5) + +initial_state, modeldata = PALEOmodel.initialize!(model) + +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# No carbonate system, so can integrate as an ODE +PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, + solvekwargs=(reltol=1e-5,), +) + + +######################################## +# Plot output +######################################## + +# individual plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# assemble plots onto screens with 6 subplots +gr(size=(1200, 900)) +pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm))) + +plot_totals(paleorun.output; species=["O2", "P"], pager=pager) +plot_airsea(paleorun.output; species=["O2"], pager=pager) +pager(:newpage) +plot_ocean_tracers( + paleorun.output; + tracers=["insol", "O2_conc", "P_conc"], + tcol=[-Inf, 10.0, 100.0, 1000.0, 1e4, 1e5], + pager=pager +) +pager(:newpage) # flush output diff --git a/examples/romglb/PALEO_examples_romglb_P_O2_S_Carb_open.jl b/examples/romglb/PALEO_examples_romglb_P_O2_S_Carb_open.jl new file mode 100644 index 0000000..d475b84 --- /dev/null +++ b/examples/romglb/PALEO_examples_romglb_P_O2_S_Carb_open.jl @@ -0,0 +1,77 @@ +using Logging +using DiffEqBase +using Sundials + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOocean +import PALEOcopse + +global_logger(ConsoleLogger(stderr,Logging.Info)) + +include("../atmreservoirreaction.jl") +include("SedimentationRate_dev.jl") +include("config_ocean_romglb_expts.jl") +include("plot_ocean_romglb.jl") + +# Atmosphere/ocean, Ccarb, Corg, S and P burial +# Test two options for solving carbonate chemistry +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_P_O2_S_Carb_open"; + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "romaniello2010_transport"), # assume zip file has been downloaded and unpacked in this subfolder + # Option 1: add pHfree, TAlk to state variables and apply constraint on TAlk + "TAlkStateExplicit"=>true, + # Option 2: TAlk is an implicit variable, pHfree is a state variable + # "TAlkStateExplicit"=>false, + ), +) + +# additional configuration to set shelf and deep carbonate burial +config_ocean_romglb_expts( + model, + [ + "set_carbonate_config_79box", # configure shelf and deep carbonate burial for modern Earth + ], +) +tspan=(0,1e5) # constraint + +initial_state, modeldata = PALEOmodel.initialize!(model) + +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# DAE solver required when using carbonate system +PALEOmodel.ODE.integrateDAEForwardDiff( + paleorun, initial_state, modeldata, tspan, + solvekwargs=( + reltol=1e-5, + dtmax=1e5, + ), +) + + +######################################## +# Plot output +######################################## + +# individual plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# assemble plots onto screens with 6 subplots +gr(size=(1200, 900)) +pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm))) + +plot_totals(paleorun.output; species=["P", "S"], pager=pager) +plot_airsea(paleorun.output; pager=pager) +pager(:newpage) +plot_ocean_tracers( + paleorun.output; + tracers=["insol", "O2_conc", "P_conc", "H2S_conc", "DIC_conc", "TAlk_conc", "pHtot", "OmegaAR"], + tcol=[-Inf, 10.0, 100.0, 1000.0, 1e4, 1e5], + pager=pager +) +plot_ocean_burial(paleorun.output; pager=pager) +pager(:newpage) # flush output diff --git a/examples/romglb/PALEO_examples_romglb_P_O2_open.jl b/examples/romglb/PALEO_examples_romglb_P_O2_open.jl new file mode 100644 index 0000000..e17953d --- /dev/null +++ b/examples/romglb/PALEO_examples_romglb_P_O2_open.jl @@ -0,0 +1,64 @@ +using Logging +using DiffEqBase +using Sundials + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOocean + +global_logger(ConsoleLogger(stderr,Logging.Info)) + + +include("../atmreservoirreaction.jl") +include("SedimentationRate_dev.jl") +include("config_ocean_romglb_expts.jl") +include("plot_ocean_romglb.jl") + +# Atmosphere/ocean, O2 and P only, open system with P and Corg burial +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_P_O2_open", + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "romaniello2010_transport") # assume zip file has been downloaded and unpacked in this subfolder + ) +) + +config_ocean_romglb_expts(model, ["baseline"]); tspan=(0, 1e4) + +initial_state, modeldata = PALEOmodel.initialize!(model) + +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# No carbonate system, so can integrate as an ODE +PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, + solvekwargs=(reltol=1e-5,), +) + + + +######################################## +# Plot output +######################################## + +# individual plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# assemble plots onto screens with 6 subplots +gr(size=(1200, 900)) +pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm))) + + +plot_ocean_tracers( + paleorun.output; + tracers=["insol", "O2_conc", "P_conc"], + tcol=[-Inf, 10.0, 100.0, 1000.0, 1e4, 1e5], + pager=pager +) +pager(:newpage) # flush output +plot_totals(paleorun.output; species=["O2", "P"], pager=pager) +plot_airsea(paleorun.output; species=["O2"], pager=pager) +plot_ocean_burial_CorgP(paleorun.output; pager=pager) +pager(:newpage) # flush output \ No newline at end of file diff --git a/examples/romglb/PALEO_examples_romglb_cfg.yaml b/examples/romglb/PALEO_examples_romglb_cfg.yaml new file mode 100644 index 0000000..bdddec1 --- /dev/null +++ b/examples/romglb/PALEO_examples_romglb_cfg.yaml @@ -0,0 +1,1045 @@ +#################################### +# Abiotic atmosphere/ocean, O2 only +#################################### +romglb_abiotic_O2: + parameters: + # CIsotope: ScalarData + matdir: romaniello2010_transport # folder with data files from Romaniello (2010) GGG + domains: + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["O2"] + + global: + # scalar domain + + reactions: + total_O2: + class: ReactionSum + parameters: + vars_to_add: [atm.O2, ocean.O2] + variable_links: + sum: total_O2 + atm: + + + reactions: + reservoir_O2: + class: ReactionReservoirAtm + + variable_links: + R*: O2* + pRatm: pO2atm + pRnorm: pO2PAL + variable_attributes: + R:norm_value: 3.7e19 # present-day atmospheric level + R:initial_value: 3.7e19 + + transfer_AtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + transfer_multiplier: -1.0 + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: $fluxname$_sms + + oceansurface: + reactions: + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston: 4.8 # m d-1 + + surface_insol: + class: ReactionConst + parameters: + constnames: [surface_insol] + variable_attributes: + surface_insol:initial_value: [70.0, 350.0, 350.0] # W m-2 downwelling radiative flux + + transfer_AtmtoOceansurface: + 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 + + ocean: + reactions: + transportromglb: + class: ReactionOceanTransportRomanielloICBM + parameters: + matdir: external%matdir + circname: Global_79_Box + + reservoir_O2: + class: ReactionReservoir + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_T: + class: ReactionReservoirTotal + variable_links: + R*: T* + variable_attributes: + R:initial_value: 1.0 # concentration m-3 + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.04 + + + + oceanfloor: + reactions: + +################################# +# Atmosphere/ocean, O2 and P only +################################# + +romglb_P_O2: + parameters: + CIsotope: ScalarData + SIsotope: ScalarData + matdir: romaniello2010_transport # folder with data files from Romaniello (2010) GGG + domains: + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["O2"] + + fluxOceanfloor: + reactions: + particulatetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: particulateflux_ + fluxlist: ["Corg::CIsotope", "N", "P", "Ccarb::CIsotope"] + + solutetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: soluteflux_ + fluxlist: ["DIC::CIsotope", "TAlk", "O2", "P"] + + global: + # scalar domain + + reactions: + total_O2: + class: ReactionSum + parameters: + vars_to_add: [atm.O2, ocean.O2] + variable_links: + sum: total_O2 + atm: + + + reactions: + reservoir_O2: + class: ReactionReservoirAtm + + variable_links: + R*: O2* + pRatm: pO2atm + pRnorm: pO2PAL + variable_attributes: + R:norm_value: 3.7e19 # present-day atmospheric level + R:initial_value: 3.7e19 + + transfer_AtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + transfer_multiplier: -1.0 + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: $fluxname$_sms + + oceansurface: + reactions: + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston: 4.8 # m d-1 + + surface_insol: + class: ReactionConst + parameters: + constnames: [surface_insol] + variable_attributes: + surface_insol:initial_value: [70.0, 350.0, 350.0] # W m-2 downwelling radiative flux + + transfer_AtmtoOceansurface: + 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 + + ocean: + reactions: + transportromglb: + class: ReactionOceanTransportRomanielloICBM + parameters: + matdir: external%matdir + circname: Global_79_Box + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_P: + class: ReactionReservoirTotal + variable_links: + R*: P* + variable_attributes: + R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.0 # zero opacity as we just want to prescribe insol in surface boxes + + bioprod: + class: ReactionBioProdMMPop + parameters: + depthlimit: -10.0 # surface cells only + + rCorgPO4: 106.0 + rNPO4: 16.0 + rCcarbCorg: 0.25 + rCcarbCorg_fixed: true + + nuDOM: 0.0 # no DOM pool + + k_poptype: Constant + k_uPO4: 2.615e-3 # 1.91e-6 mol / kg-sw / yr * 1027 kg m-3 * 4/3 + + k_nuttype: PO4MM + k_KPO4: 0.21567e-3 # 0.21e-6 mol / kg-sw * 1027 kg m-3 + + k_lightlim: linear + k_Irel: 1.0 + + variable_links: + partprod_*: export_* + + biopump: + class: ReactionExportDirectColumn + + parameters: + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] + transportfloor: true + exportfunction: SumExp + input_frac: [0.9354, 0.0646] # 2-G model of Ridgwell & Hargreaves (2007) + sumexp_scale: [550.5195, 1e6] + + reminocean: + class: ReactionReminO2 + + parameters: + + variable_links: + soluteflux_*: "*_sms" + + + oceanfloor: + reactions: + + reminoceanfloor: + class: ReactionReminO2 + + parameters: + + variable_links: + remin*: particulateflux* + soluteflux_*: fluxOceanfloor.soluteflux_* + + transferparticulate_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: particulateflux_$fluxname$ + + transfersolute_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ + output_fluxes: ocean.oceanfloor.$fluxname$_sms + + sedcrust: + + +#################################################################################################### +# Atmosphere/ocean, O2 and P only, open system with P and Corg burial +# This is configured to run to a steady-state, with restoring fluxes balancing burial +# Atmosphere O2 is restored to a constant value (so restoring flux will balance net O2 production from Corg burial) +# Ocean P is restored to constant values (so restoring fluxe will balance P burial) +################################################################################################# + +romglb_P_O2_open: + parameters: + CIsotope: ScalarData + SIsotope: ScalarData + matdir: romaniello2010_transport # folder with data files from Romaniello (2010) GGG + domains: + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["O2"] + + fluxRtoOcean: + + reactions: + transfer: + class: ReactionFluxTarget + parameters: + fluxlist: ["P"] # for restoring flux + + fluxOceanfloor: + reactions: + particulatetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: particulateflux_ + fluxlist: ["Corg::CIsotope", "N", "P", "Ccarb::CIsotope"] + + solutetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: soluteflux_ + fluxlist: ["DIC::CIsotope", "TAlk", "O2", "P"] + + fluxOceanBurial: + reactions: + transfer: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["Corg::CIsotope", "P", "Pauth", "PFe", "Porg"] + + global: + # scalar domain + + reactions: + total_O2: + class: ReactionSum + parameters: + vars_to_add: [atm.O2, ocean.O2] + variable_links: + sum: total_O2 + atm: + + reactions: + reservoir_O2: + class: ReactionReservoirAtm + + variable_links: + R*: O2* + pRatm: pO2atm + pRnorm: pO2PAL + variable_attributes: + R:norm_value: 3.7e19 # present-day atmospheric level + R:initial_value: 3.7e19 + + restoreO2: + class: ReactionRestore + + parameters: + RequiredLevel: 3.7e19 # PAL + trestore: 100.0 + source_only: false + variable_links: + WatchLevel: O2 + RestoringFlux: O2_sms + RestoringApplied: O2_restoring + + + transfer_AtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + transfer_multiplier: -1.0 + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: $fluxname$_sms + + oceansurface: + reactions: + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston: 4.8 # m d-1 + + surface_insol: + class: ReactionConst + parameters: + constnames: [surface_insol] + variable_attributes: + surface_insol:initial_value: [70.0, 350.0, 350.0] # W m-2 downwelling radiative flux + + transfer_AtmtoOceansurface: + 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 + + transfer_RtoOcean: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + input_fluxes: fluxRtoOcean.flux_$fluxname$ + output_fluxes: ocean.oceansurface.$fluxname$_sms + + ocean: + reactions: + transportromglb: + class: ReactionOceanTransportRomanielloICBM + parameters: + matdir: romaniello2010_transport # folder with data files from Romaniello (2010) GGG + circname: Global_79_Box + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_P: + class: ReactionReservoirTotal + variable_links: + R*: P* + variable_attributes: + R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.0 # zero opacity as we just want to prescribe insol in surface boxes + + bioprod: + class: ReactionBioProdMMPop + parameters: + depthlimit: -10.0 # surface cells only + + rCorgPO4: 106.0 + rNPO4: 16.0 + rCcarbCorg: 0.25 + rCcarbCorg_fixed: true + + nuDOM: 0.0 # no DOM pool + + k_poptype: Constant + k_uPO4: 2.615e-3 # 1.91e-6 mol / kg-sw / yr * 1027 kg m-3 * 4/3 + + k_nuttype: PO4MM + k_KPO4: 0.21567e-3 # 0.21e-6 mol / kg-sw * 1027 kg m-3 + + k_lightlim: linear + k_Irel: 1.0 + + variable_links: + partprod_*: export_* + + biopump: + class: ReactionExportDirectColumn + + parameters: + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] + transportfloor: true + exportfunction: SumExp + input_frac: [0.9354, 0.0646] # 2-G model of Ridgwell & Hargreaves (2007) + sumexp_scale: [550.5195, 1e6] + + reminocean: + class: ReactionReminO2 + + parameters: + + variable_links: + soluteflux_*: "*_sms" + + restoreP: + class: ReactionRestore + + parameters: + RequiredLevel: 3.466e15 # mol 1.570e18 m^3 * 2.208e-3 mol m-3 + trestore: 100.0 + source_only: true + variable_links: + WatchLevel: P_total + RestoringFlux: fluxRtoOcean.flux_P + + + oceanfloor: + reactions: + + sedrate: + class: ReactionSedimentationRate_dev + parameters: + sed_rate_function: Tromp1995 + + sedBEcorgP: + class: ReactionBurialEffCorgP + parameters: + burial_eff_function: Ozaki2011 + + # overall normalization for Corg burial to give 2.5e12 mol C yr-1 + BECorgNorm: 0.145 # 2.5e12/17.22e12 + + # P:Corg burial ratios to give 4e+10 mol P yr-1 + BPorgCorg: [4.0e-3] # prescribed P:Corg + BPFeCorg: [4.0e-3] # prescribed P:Corg + BPauthCorg: [8.0e-3] # prescribed P:Corg + variable_links: + reminflux_Corg: remin_Corg # -> reminoceanfloor + reminflux_N: remin_N # -> reminoceanfloor + reminflux_P: remin_P # -> reminoceanfloor + + reminoceanfloor: + class: ReactionReminO2 + + parameters: + + variable_links: + remin_Corg: remin_Corg + remin_N: remin_N + remin_P: remin_P + remin_Ccarb: particulateflux_Ccarb # no burial Reaction, so remin everything + + soluteflux_*: fluxOceanfloor.soluteflux_* + + transferparticulate_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: particulateflux_$fluxname$ + + transfersolute_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ + output_fluxes: ocean.oceanfloor.$fluxname$_sms + + sedcrust: + + + +################################################################################################ +# Atmosphere/ocean, Ccarb, Corg, and P burial +# This is configured to run to a steady-state, with restoring fluxes balancing burial +# Atmosphere O2 is restored to a constant value (so restoring flux will balance net O2 production from Corg burial) +# Ocean P is restored to a constant value (so restoring fluxes will balance P burial) +# No pyrite or gypsum burial fluxes (so ocean is a closed system to sulphur) +# Ocean TAlk and atmosphere pCO2 are restored to constant values (so restoring fluxes will balance CaCO3 burial) +########################################################################################################### + +romglb_P_O2_S_Carb_open: + parameters: + CIsotope: ScalarData + SIsotope: ScalarData + # Define the solution method for pH/TAlk + # Option 1: add pHfree, TAlk to state variables and apply constraint on TAlk + TAlkStateExplicit: true + # Option 2: TAlk is an implicit variable, pHfree is a state variable + # TAlkStateExplicit: false + + matdir: romaniello2010_transport # folder with data files from Romaniello (2010) GGG + domains: + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["O2", "CO2::CIsotope"] + + fluxRtoOcean: + + reactions: + transfer: + class: ReactionFluxTarget + parameters: + fluxlist: ["P", "TAlk", "SO4::SIsotope"] # for restoring fluxes + + fluxOceanfloor: + reactions: + particulatetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: particulateflux_ + fluxlist: ["Corg::CIsotope", "N", "P", "Ccarb::CIsotope"] + + solutetarget: + class: ReactionFluxTarget + parameters: + flux_totals: true + target_prefix: soluteflux_ + fluxlist: ["DIC::CIsotope", "TAlk", "O2", "P", "SO4::SIsotope", "H2S::SIsotope"] + + + fluxOceanBurial: + reactions: + transfer: + class: ReactionFluxTarget + parameters: + flux_totals: true + fluxlist: ["Corg::CIsotope", "Ccarb::CIsotope", "GYP::SIsotope", "PYR::SIsotope", "P", "Pauth", "PFe", "Porg"] + + global: + + reactions: + force_RHOSFW: # constant 1.0 + class: ReactionForceInterp + variable_links: + F: RHOSFW + + total_S: + class: ReactionSum + parameters: + vars_to_add: [ocean.SO4_total, ocean.H2S_total] + variable_links: + sum: total_S + + total_O2: + class: ReactionSum + parameters: + vars_to_add: [atm.O2, ocean.O2] + variable_links: + sum: total_O2 + + atm: + + + reactions: + reservoir_O2: + class: ReactionReservoirAtm + + variable_links: + R*: O2* + pRatm: pO2atm + pRnorm: pO2PAL + variable_attributes: + R:norm_value: 3.7e19 # present-day atmospheric level + R:initial_value: 3.7e19 + + reservoir_CO2: + class: ReactionReservoirAtm + parameters: + field_data: external%CIsotope + + variable_links: + R*: CO2* + pRatm: pCO2atm + pRnorm: pCO2PAL + variable_attributes: + R:norm_value: 4.956e16 # pre ind 280e-6 + R:initial_value: 4.956e16 # pre ind 280e-6 + + restoreO2: + class: ReactionRestore + + parameters: + RequiredLevel: 3.7e19 # PAL + trestore: 100.0 + source_only: false + variable_links: + WatchLevel: O2 + RestoringFlux: O2_sms + RestoringApplied: O2_restoring + + restoreCO2: + class: ReactionRestore + + parameters: + RequiredLevel: 4.956e16 # pre ind 280e-6 + trestore: 100.0 + source_only: true + variable_links: + WatchLevel: CO2 + RestoringFlux: CO2_sms + RestoringApplied: CO2_restoring + + transfer_AtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + transfer_multiplier: -1.0 + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: $fluxname$_sms + + oceansurface: + reactions: + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston: 4.8 # m d-1 + + airsea_CO2: + class: ReactionAirSeaCO2 + parameters: + piston: 4.8 # m d-1 + + + variable_links: + # Xatm_delta: atm.CO2_delta + # Xocean_delta: ocean.oceansurface.DIC_delta + + surface_insol: + class: ReactionConst + parameters: + constnames: [surface_insol] + variable_attributes: + surface_insol:initial_value: [70.0, 350.0, 350.0] # W m-2 downwelling radiative flux + + transfer_RtoOcean: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + input_fluxes: fluxRtoOcean.flux_$fluxname$ + output_fluxes: ocean.oceansurface.$fluxname$_sms + + 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 + + ocean: + reactions: + transportromglb: + class: ReactionOceanTransportRomanielloICBM + parameters: + matdir: romaniello2010_transport # folder with data files from Romaniello (2010) GGG + circname: Global_79_Box + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_P: + class: ReactionReservoirTotal + variable_links: + R*: P* + variable_attributes: + R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) + + reservoir_SO4: + class: ReactionReservoirTotal + parameters: + field_data: external%SIsotope + variable_links: + R*: SO4* + variable_attributes: + R:initial_value: 28756.0e-3 # concentration mol m-3 = 28e-3 mol/kg * 1027 kg m-3 + R:initial_delta: 0.0 + R:norm_value: 1000.0e-3 # for scaling only + + reservoir_H2S: + class: ReactionReservoirTotal + parameters: + field_data: external%SIsotope + variable_links: + R*: H2S* + variable_attributes: + R:initial_value: 1e-6 # concentration mol m-3 ~ 1e-9 mol/kg * 1027 kg m-3 + R:initial_delta: 0.0 + R:norm_value: 1.0e-3 # for scaling only + + reservoir_Ca: + class: ReactionReservoirConst + + variable_links: + R*: Ca* + + variable_attributes: + R_conc:initial_value: 10.56 # mol m-3 modern Ca (10.28e-3*1027 kg m-3) + + reservoir_DIC: + class: ReactionReservoirTotal + parameters: + field_data: external%CIsotope + variable_links: + R*: DIC* + variable_attributes: + R:initial_value: 2208.1e-3 # 2150e-6 mol kg-1 * 1027 kg m-3 + R:initial_delta: -1.0 + R:norm_value: 1000.0e-3 # for scaling only + + # TAlk / pHfree configuration: + # either: + # - reservoir_TAlk + TAlk_constraint_H_primary_species + # or + # - TAlk_total_H_primary_species + # + reservoir_TAlk: + enabled: external%TAlkStateExplicit + class: ReactionReservoirTotal + + variable_links: + R*: TAlk* + variable_attributes: + R:initial_value: 2403.2e-3 # 2340e-6 mol kg-1 * 1027 kg m-3 + R:norm_value: 1000.0e-3 # for scaling only + + # [H+] primary species (as pHfree) for algebraic constraint on TAlk + TAlk_constraint_H_primary_species: + enabled: external%TAlkStateExplicit + class: ReactionConstraintReservoir + variable_links: + Primary_pconc: pHfree + Primary_conc: H_conc + R*: TAlk* + parameters: + primary_total_stoich: 0.0 # ReactionCO2SYS adds H to TAlk_calc + primary_variable: p_concentration # provide pHfree as state variable to solver + constraint_variable: amount # provide TAlk_constraint (mol) as algebraic constraint to solver + variable_attributes: + Primary_pconc%initial_value: 8.0 + Primary_pconc%norm_value: 1.0 + R_constraint%norm_value: 1.0 + + # TA as total and H (as pH) as primary species + TAlk_total_H_primary_species: + disabled: external%TAlkStateExplicit + class: ReactionImplicitReservoir + variable_links: + Primary_pconc: pHfree + Primary_conc: H_conc + R*: TAlk* + parameters: + primary_total_stoich: 0.0 # ReactionCO2SYS adds H to TAlk_calc + primary_variable: p_concentration # provide solver with -log10(H_conc) + total_variable: amount # provide TAlk (mol) to solver + total: true # add R_total = sum(R) + variable_attributes: + Primary_pconc%initial_value: 8.0 + Primary_pconc%norm_value: 1.0 + R_conc%norm_value: 1.0 + + carbchem: + class: ReactionCO2SYS + parameters: + components: ["Ci", "B", "S", "F", "Omega", "H2S"] + defaultconcs: ["TF", "TB"] #, "Ca"] + solve_pH: speciationTAlk + outputs: ["pCO2", "xCO2dryinp", "CO2", "CO3", "OmegaCA", "OmegaAR"] + variable_links: + TCi_conc: DIC_conc + TS_conc: SO4_conc + CO2: CO2_conc + CO3: CO3_conc + TH2S_conc: H2S_conc + pCO2: pCO2 + OmegaAR: OmegaAR + OmegaCA: OmegaCA + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.0 # zero opacity as we just want to prescribe insol in surface boxes + + bioprod: + class: ReactionBioProdMMPop + parameters: + depthlimit: -10.0 # surface cells only + + rCorgPO4: 106.0 + rNPO4: 16.0 + rCcarbCorg: 0.25 + rCcarbCorg_fixed: true + + nuDOM: 0.0 # no DOM pool + + k_poptype: Constant + k_uPO4: 2.615e-3 # 1.91e-6 mol / kg-sw / yr * 1027 kg m-3 * 4/3 + + k_nuttype: PO4MM + k_KPO4: 0.21567e-3 # 0.21e-6 mol / kg-sw * 1027 kg m-3 + + k_lightlim: linear + k_Irel: 1.0 + + variable_links: + Prod_Corg_total: Prod_Corg_total + Prod_Ccarb_total: Prod_Ccarb_total + partprod_*: export_* + + biopumpCorg: + class: ReactionExportDirectColumn + + parameters: + fluxlist: ["P", "N", "Corg::CIsotope"] + transportfloor: true + exportfunction: SumExp + input_frac: [0.9354, 0.0646] # 2-G model of Ridgwell & Hargreaves (2007) + sumexp_scale: [550.5195, 1e6] + + biopumpCcarb: + class: ReactionExportDirectColumn + + parameters: + fluxlist: ["Ccarb::CIsotope"] + transportfloor: true + exportfunction: SumExp + input_frac: [1.0] # single inorganic carbon fraction + sumexp_scale: [2000.0] + + remin: + class: ReactionReminO2_SO4 + + variable_links: + soluteflux_*: "*_sms" + + redox_H2S_O2: + class: ReactionRedoxH2S_O2 + + parameters: + R_H2S_O2: 3.65e3 # (mol m-3) yr-1 + + restoreP: + class: ReactionRestore + + parameters: + RequiredLevel: 3.466e15 # mol 1.570e18 m^3 * 2.208e-3 mol m-3 + trestore: 100.0 + source_only: true + variable_links: + WatchLevel: P_total + RestoringFlux: fluxRtoOcean.flux_P + + restoreTAlk: + class: ReactionRestore + + parameters: + RequiredLevel: 3.773e18 # mol 1.570e18 m^3 * 2403.2e-3 mol m-3 + trestore: 100.0 + source_only: false + variable_links: + WatchLevel: TAlk_total + RestoringFlux: fluxRtoOcean.flux_TAlk + + + oceanfloor: + reactions: + sedrate: + class: ReactionSedimentationRate_dev + parameters: + sed_rate_function: Tromp1995 + + sedBEcorgP: + class: ReactionBurialEffCorgP + parameters: + burial_eff_function: Ozaki2011 + + # overall normalization for Corg burial to give 2.5e12 mol C yr-1 + BECorgNorm: 0.145 # 2.5e12/17.22e12 + + # P:Corg burial ratios to give 4e+10 mol P yr-1 + BPorgCorg: [4.0e-3] # prescribed P:Corg + BPFeCorg: [4.0e-3] # prescribed P:Corg + BPauthCorg: [8.0e-3] # prescribed P:Corg + variable_links: + reminflux_Corg: remin_Corg # -> reminoceanfloor + reminflux_N: remin_N # -> reminoceanfloor + reminflux_P: remin_P # -> reminoceanfloor + + shelfcarb: + class: ReactionShelfCarb + parameters: + shelfareanorm: [.NaN] # requires configuration to set only low lat surface box + carbsedshallow: 0.987e12 # for 6.66e12 mol C yr-1 + + deepcarb: + class: ReactionBurialEffCarb + parameters: + burial_eff_function: OmegaCA + hascarbseddeep: [false] # requires configuration to set deep boxes + m0: 2.0 + m1: 1.8 # guesses to get 12.77e12 mol C yr-1 + variable_links: + particulateflux_Ccarb: particulateflux_Ccarb + + + reminoceanfloor: + class: ReactionReminO2_SO4 + + parameters: + + variable_links: + # remin_Ccarb: particulateflux_Ccarb # handled by deepcarb + + O2_conc: ocean.oceanfloor.O2_conc + # SO4_delta: ocean.oceanfloor.SO4_delta + + soluteflux_*: fluxOceanfloor.soluteflux_* + + sfw: + class: ReactionSeafloorWeathering + + parameters: + k_sfw : 3e12 # mol yr-1 rate normalisation + f_sfw_force : None # tectonic forcing + f_sfw_opt : sfw_noT # + sfw_distribution_method : Depth + sfw_depthrange: [-3000.0, -1000.0] # guess depth range to distribute sfw flux + + transferparticulate_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: particulateflux_$fluxname$ + + transfersolute_fluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ + output_fluxes: ocean.oceanfloor.$fluxname$_sms + + sedcrust: + + diff --git a/examples/romglb/README.md b/examples/romglb/README.md new file mode 100644 index 0000000..b68ebe2 --- /dev/null +++ b/examples/romglb/README.md @@ -0,0 +1,56 @@ +# Intermediate-complexity global ocean + +Intermediate-complexity global ocean examples and test cases, +using the 3-column box model for the global ocean circulation from +[Romaniello2010](@cite). + +## Installation + +These examples use `ReactionOceanTransportRomanielloICBM` to read Matlab data files with the +3-column box model circulation from [Romaniello2010](@cite). + +The Matlab datafiles are available as a zip file from , +generated from the Matlab model code available as Supplementary Information to [Romaniello2010](@cite). + +The examples assume the zip file has been downloaded and unpacked to subfolder `romaniello2010_transport` + +## O2 only air-sea exchange and transport test + + include("PALEO_examples_romglb_O2_only.jl") + +## P, O2 with parameterized export production + + include("PALEO_examples_romglb_P_O2.jl") + +Phosphorus and oxygen, with a parameterization of +export production based on light and nutrient availability. + +No burial fluxes (ie oceanfloor phosphorus flux is recycled into water column), +so the ocean is effectively a closed system for phosphorus. + +## P, O2 with organic carbon and phosphorus burial + + include("PALEO_examples_romglb_P_O2_open.jl") + +Phosphorus and oxygen, with a parameterization of +export production based on light and nutrient availability. + +Burial efficiency parameterization for burial fluxes of organic carbon and phosphorus, +with ocean phosphorus restored to modern level. + +## P, O2, S, DIC with organic carbon, phosphorus, and carbonate burial + + include("PALEO_examples_romglb_P_O2_S_Carb_open.jl") + +Phosphorus, oxygen, sulphur and DIC, with a parameterization of +export production based on light and nutrient availability. + +Burial efficiency parameterization for burial fluxes of organic carbon and phosphorus, +with ocean phosphorus restored to modern level. + +Burial efficiency parameterisation of carbonate burial, with shelf carbonate burial +based on saturation state and shelf area, and deep carbonate burial based on oceanfloor +flux and saturation state. Ocean TAlk and atmosphere pCO2 are restored to constant values +(so restoring fluxes will balance CaCO3 burial). + +No pyrite or gypsum burial (so ocean is a closed system for sulphur). \ No newline at end of file diff --git a/examples/romglb/SedimentationRate_dev.jl b/examples/romglb/SedimentationRate_dev.jl new file mode 100644 index 0000000..23ad854 --- /dev/null +++ b/examples/romglb/SedimentationRate_dev.jl @@ -0,0 +1,108 @@ +module SedimentationRate_dev + + +import PALEOboxes as PB +using PALEOboxes.DocStrings + +using SpecialFunctions +using Interpolations + +import Infiltrator + + + +""" + sedRate_Tromp1995(depth) -> sedimentation_rate + +`sedimentation_rate` (m/yr) at `depth` (m, -ve), [Tromp1995](@cite) +eg used by [Ozaki2011](@cite) + +# Examples: +Check value at depth 1000m + ```jldoctest +julia> round(Main.SedimentationRate_dev.sedRate_Tromp1995(-1000.0), sigdigits=5) +0.00034152 +julia> round(Main.SedimentationRate_dev.Burial.sedRate_Tromp1995(-100.0), sigdigits=5) +0.002369 +``` +""" +function sedRate_Tromp1995(depth) + # eqn (3) sed rate: (NB: log = log10) + + sedimentation_rate = NaN*one(depth) + + if (depth < -5400) + sedimentation_rate = zero(depth) + elseif depth < 0 + sedimentation_rate = 1e-2*10^(erfcinv(-depth/2700)-2.1) + end + + return sedimentation_rate +end + +""" + ReactionSedimentationRate_dev + +Sedimentation rate parameterized from water depth + +# Parameters +$(PARS) + +# Methods and Variables +$(METHODS_SETUP) +""" +Base.@kwdef mutable struct ReactionSedimentationRate_dev{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + PB.ParString("sed_rate_function", "Tromp1995", allowed_values=["Tromp1995"], + description="sedimentation rate parameterisation"), + ) + +end + +function PB.register_methods!(rj::ReactionSedimentationRate_dev) + + sr_Tromp1995(pars, vars, i) = sedRate_Tromp1995(vars.zlower[i]) + + if rj.pars.sed_rate_function[] == "Tromp1995" + sr_func = sr_Tromp1995 + else + error("configuration error: unknown sed_rate_function $(rj.pars.sed_rate_function[])") + end + PB.setfrozen!(rj.pars.sed_rate_function) + + vars = [ + PB.VarDepStateIndep("ocean.oceanfloor.zlower", "m", "depth of lower surface of box (m)"), + PB.VarPropStateIndep("sedimentation_rate", "m yr-1", "sedimentation rate"), + ] + + PB.add_method_setup!( + rj, + setup_sedimentation_rate, + (PB.VarList_namedtuple(vars), ), + p = sr_func, + ) + return nothing +end + +function setup_sedimentation_rate( + m::PB.ReactionMethod, + pars, + (vars, ), + cellrange::PB.AbstractCellRange, + attribute_name +) + attribute_name == :setup || return nothing + + sr_func = m.p + + @inbounds for i in cellrange.indices + vars.sedimentation_rate[i] = sr_func(pars, vars, i) + end + + return nothing +end + + +end diff --git a/examples/romglb/Test xarray netcdf.ipynb b/examples/romglb/Test xarray netcdf.ipynb new file mode 100644 index 0000000..3edd7d3 --- /dev/null +++ b/examples/romglb/Test xarray netcdf.ipynb @@ -0,0 +1,224 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "68667024", + "metadata": {}, + "source": [ + "# Example reading netcdf output using python xarray\n", + "\n", + "Generate file 'test.nc' with:\n", + "\n", + " julia> include(\"PALEO_examples_romglb_P_O2_S_Carb_open.jl\")\n", + " julia> PALEOmodel.OutputWriters.save_netcdf(run.output, \"test.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d640c5a0", + "metadata": {}, + "outputs": [], + "source": [ + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9a1bd15", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5ddcd9d", + "metadata": {}, + "outputs": [], + "source": [ + "os.getcwd()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "401e5aef", + "metadata": {}, + "outputs": [], + "source": [ + "# read top-level dataset\n", + "ds = xr.open_dataset(\"test.nc\")\n", + "ds.attrs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b80a015", + "metadata": {}, + "outputs": [], + "source": [ + "# read ocean group (= PALEO Domain ocean) from netcdf file\n", + "ds_ocean = xr.open_dataset(\"test.nc\", group=\"ocean\")\n", + "\n", + "# attach z coordinates (this is not automatic!)\n", + "ds_ocean = ds_ocean.set_coords([\"zmid\", \"zlower\", \"zupper\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ef1c0ef", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "ds_ocean" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11dcd1af", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# timeseries for a scalar variable\n", + "print(ds_ocean[\"O2_total\"])\n", + "ds_ocean[\"O2_total\"].plot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "853eb848", + "metadata": {}, + "outputs": [], + "source": [ + "# split dataset into columns\n", + "cells_in_column = ds_ocean[\"cells_in_column\"]\n", + "\n", + "print(\"Icolumn comment:\", ds_ocean[\"Icolumns\"].attrs[\"comment\"]) # zero-based indices of cells from top to bottom ordered by columns\n", + "Icolumns = ds_ocean[\"Icolumns\"].values\n", + "\n", + "ds_ocean_columns = []\n", + "colstartidx = 0\n", + "\n", + "for cidx, clength in enumerate(cells_in_column.values):\n", + " ccells = Icolumns[np.arange(colstartidx, colstartidx + clength)]\n", + " ds_ocean_columns.append(ds_ocean.isel(cells=ccells, columns=cidx))\n", + " colstartidx += clength\n", + " \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68d49cf3", + "metadata": {}, + "outputs": [], + "source": [ + "ds_ocean_columns[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47459078", + "metadata": {}, + "outputs": [], + "source": [ + "ds_ocean_columns[0].columnnames.values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aae8816e", + "metadata": {}, + "outputs": [], + "source": [ + "print(ds_ocean_columns[0][\"O2_conc\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "799531d6", + "metadata": {}, + "outputs": [], + "source": [ + "ds_ocean_columns[0][\"O2_conc\"].dims" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bdc14198", + "metadata": {}, + "outputs": [], + "source": [ + "ds_ocean_columns[0][\"O2_conc\"].coords" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd95a3fa", + "metadata": {}, + "outputs": [], + "source": [ + "np.shape(ds_ocean_columns[0][\"O2_conc\"].values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e87c215", + "metadata": {}, + "outputs": [], + "source": [ + "ds_ocean_columns[0][\"O2_conc\"].sel(tmodel=1e5).plot(y=\"zmid\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec691d08", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/romglb/config_ocean_romglb_expts.jl b/examples/romglb/config_ocean_romglb_expts.jl new file mode 100644 index 0000000..154ab04 --- /dev/null +++ b/examples/romglb/config_ocean_romglb_expts.jl @@ -0,0 +1,26 @@ + +function config_ocean_romglb_expts( + model, expts; +) + for expt in expts + if expt == "baseline" + # baseline configuration + + elseif expt == "set_carbonate_config_79box" + shelfareanorm = zeros(79) + shelfareanorm[14] = 1.0 # low lat (:gyre) surface box only + r_shelfcarb = PB.get_reaction(model, "oceanfloor", "shelfcarb") + PB.setvalue!(r_shelfcarb.pars.shelfareanorm, shelfareanorm) + + hascarbseddeep = ones(Bool, 79) + hascarbseddeep[[1, 14, 47]] .= false # disable deposition in surface boxes ('shelves') + r_deepcarb = PB.get_reaction(model, "oceanfloor", "deepcarb") + PB.setvalue!(r_deepcarb.pars.hascarbseddeep, hascarbseddeep) + + else + error("unknown expt ", expt) + end + end + + return model +end diff --git a/examples/romglb/plot_ocean_romglb.jl b/examples/romglb/plot_ocean_romglb.jl new file mode 100644 index 0000000..9303850 --- /dev/null +++ b/examples/romglb/plot_ocean_romglb.jl @@ -0,0 +1,185 @@ +function plot_totals( + output; + species=["T", "P", "O2", "S"], + pager=PALEOmodel.DefaultPlotPager(), + plotargs=NamedTuple(), +) + # conservation checks + + for sp in species + if PB.has_variable(output, "global.total_$sp") + pager(plot(title="Total $sp (global)", output, ["global.total_$sp"]; ylabel="$sp (mol)",), plotargs...) + else # ocean only + pager(plot(title="Total $sp (ocean)", output, ["ocean.$(sp)_total"]; ylabel="$sp (mol)",), plotargs...) + end + end + + return nothing +end + +function plot_airsea( + output; + columns=[:hlat, :gyre, :upw], + species=["O2", "CO2"], + pager=PALEOmodel.DefaultPlotPager(), + plotargs=NamedTuple(), +) + for sp in species + pager((plot(title="Air-sea $sp", output, "fluxAtmtoOceansurface.flux_$sp", (cell=columns,); plotargs...); + plot!(output, "fluxAtmtoOceansurface.flux_total_$sp"; ylabel="$sp flux (mol yr-1)", plotargs...))) + end + + return nothing +end + +function plot_ocean_tracers( + output; + columns=[:hlat, :gyre, :upw], + colskip=0, # number of blank panels to add (use to format pages) + tracers=[ + "insol", "DIC_conc", "TAlk_conc", "temp", "pHtot", "O2_conc", "SO4_conc", "H2S_conc", "CH4_conc","P_conc", "Ca_conc", + "H2S_delta", "SO4_delta", "CH4_delta", "OmegaAR"], + tcol=[-Inf, 10.0, 100.0, 1000.0, Inf], # model time for column plots + pager=PALEOmodel.DefaultPlotPager(), + plotargs=NamedTuple(), +) + for tr in tracers + for col in columns + pager(plot(title="$tr $col", output, "ocean.$tr", ( tmodel=tcol, column=col); + swap_xy=true, labelattribute=:filter_records, plotargs...)) + end + for i in 1:colskip + pager(:skip) + end + end + + return nothing +end + +function plot_ocean_tracers_ts( + output; + columns=[:gyre, :upw], + colskip=0, # number of blank panels to add (use to format pages) + tracers=["O2_conc", "H2S_conc"], + pager=PALEOmodel.DefaultPlotPager(), + plotargs=NamedTuple(), +) + bodgedplotargs = Dict(pairs(plotargs)) + if get(bodgedplotargs, :xflip, false) + @warn "plot_ocean_tracers_ts: replacing xflip with mult_x_coord=-1.0" + bodgedplotargs[:mult_x_coord]=-1.0 + delete!(bodgedplotargs, :xflip) + # fix up xlim if present + xlim = get(bodgedplotargs, :xlim, nothing) + if !isnothing(xlim) + bodgedplotargs[:xlim] = (-xlim[2], -xlim[1]) + end + end + + for tr in tracers + for col in columns + pager( + heatmap(title="Ocean $tr :$col", output, "ocean.$tr", (column=col,); bodgedplotargs...), + ) + end + for i in 1:colskip + pager(:skip) + end + end + + return nothing +end + +function plot_ocean_floor( + output; + columns=[:hlat, :gyre, :upw], + colskip=0, + tracers=["flys"], + tcol=1e12, + pager=PALEOmodel.DefaultPlotPager(), + plotargs=NamedTuple(), +) + for tr in tracers + for col in columns + p = plot(title="$tr $col"; plotargs...) + for t in tcol + t_a = PALEOmodel.get_array(output, "oceanfloor.$tr", tmodel=t, column=col,) + # add zdepth dimension + zdepth = PB.NamedDimension("zdepth", PALEOmodel.get_array(output, "oceanfloor.zfloor", tmodel=t, column=col,).values) + t_a = PB.FieldArray(tr, t_a.values, (zdepth, ), t_a.attributes) + plot!(title="$tr $col", t_a; swap_xy=true, labelattribute=:filter_records, plotargs...) + end + pager(p) + end + for i in 1:colskip + pager(:skip) + end + end + + return nothing +end + +function plot_ocean_burial_CorgP( + output; + pager=PALEOmodel.DefaultPlotPager(), + plotargs=NamedTuple(), +) + + pager( + plot(title="P input output", output, ["fluxRtoOcean.flux_P", "fluxOceanBurial.flux_total_P"]; + ylabel="P flux (mol P yr-1)", plotargs...), + plot(title="Production burial", output, ["ocean.Prod_Corg_total", "fluxOceanBurial.flux_total_Corg"], + ylabel="C flux (mol C yr-1)", plotargs...), + ( + plot(title="O2 Corg burial", output, ["fluxOceanBurial.flux_total_Corg"], + ylabel="O2, C flux (mol yr-1)", ylim=(-5e12, 5e12), plotargs...); + plot!(-1*PALEOmodel.get_array(output, "atm.O2_restoring")); + plot!(-1*PALEOmodel.get_array(output, "fluxAtmtoOceansurface.flux_total_O2")); + ) + ) + + return nothing +end + +function plot_ocean_burial( + output; + ocean_only=true, # extra restoring flux contributions + pager=PALEOmodel.DefaultPlotPager(), + plotargs=NamedTuple(), +) + + pager( + plot(title="P input output", output, ["fluxRtoOcean.flux_P", "fluxOceanBurial.flux_total_P"]; + ylabel="P flux (mol P yr-1)", plotargs...), + plot(title="Production burial", output, ["ocean.Prod_Corg_total", "ocean.Prod_Ccarb_total", "fluxOceanBurial.flux_total_Corg", "oceanfloor.deep_Ccarb_total"]; + ylabel="C flux (mol C yr-1)", plotargs...), + plot(title="S input output", output, ["fluxRtoOcean.flux_SO4", "fluxOceanBurial.flux_total_GYP", "fluxOceanBurial.flux_total_PYR"]; + ylabel="S flux (mol S yr-1)", plotargs...), + plot(title="TAlk input", output, ["fluxRtoOcean.flux_TAlk"]; + ylabel="flux (mol yr-1)", plotargs...), + ) + + if ocean_only + pager( + plot(title="Corg Pyr burial O2 flux", output, ["fluxOceanBurial.flux_total_Corg", "fluxOceanBurial.flux_total_PYR", + "atm.O2_restoring", "fluxAtmtoOceansurface.flux_total_O2"]; + ylabel="flux (mol yr-1)", plotargs...), + plot(title="C input output", output, ["fluxOceanBurial.flux_total_Corg", "fluxOceanBurial.flux_total_Ccarb", + "oceanfloor.shelf_Ccarb_total", "oceanfloor.deep_Ccarb_total", "oceanfloor.sfw_total", "atm.CO2_restoring"]; + ylabel="flux (mol yr-1)", plotargs...), + ) + else + pager( + plot(title="Corg Pyr burial O2 flux", output, ["fluxOceanBurial.flux_total_Corg", "fluxOceanBurial.flux_total_PYR", "fluxAtmtoOceansurface.flux_total_O2"]; + ylabel="flux (mol yr-1)", ylim=(-1e13, 1e13), plotargs...), + plot(title="C input output", output, ["fluxOceanBurial.flux_total_Corg", "fluxOceanBurial.flux_total_Ccarb", + "oceanfloor.shelf_Ccarb_total", "oceanfloor.deep_Ccarb_total", "oceanfloor.sfw_total"]; + ylabel="flux (mol yr-1)", plotargs...), + ) + end + + return nothing +end + + + diff --git a/examples/romglb/runtests.jl b/examples/romglb/runtests.jl new file mode 100644 index 0000000..a929171 --- /dev/null +++ b/examples/romglb/runtests.jl @@ -0,0 +1,212 @@ +using Test +using Logging +using DiffEqBase +using Sundials +import DataFrames + +import PALEOboxes as PB + +import PALEOocean +import PALEOmodel + + +@testset "romglb examples" begin + +skipped_testsets = [ + # "O2_only", + # "P_O2", + # "P_O2_S_Carb_open", + # "P_O2_S_Carb_open_implicit", # requires PALEOaqchem > v0.3.9 +] + +matdir = joinpath(@__DIR__, "romaniello2010_transport") # assume zip file has been downloaded and unpacked in this subfolder + +include("../atmreservoirreaction.jl") +include("SedimentationRate_dev.jl") + +include("config_ocean_romglb_expts.jl") + +!("O2_only" in skipped_testsets) && @testset "O2_only" begin + + model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_abiotic_O2"; + modelpars=Dict( + "matdir"=>matdir, + ) + ) + tspan=(0,1e5) + + initial_state, modeldata = PALEOmodel.initialize!(model) + paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + + PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, solvekwargs=(reltol=1e-9,) + ) + + println("conservation checks:") + conschecks = [ + ("global", "total_O2", 1e-6 ), + ] + for (domname, varname, rtol) in conschecks + startval, endval = PB.get_data(paleorun.output, domname*"."*varname)[[1, end]] + println(" check $domname.$varname $startval $endval $rtol") + @test isapprox(startval, endval, rtol=rtol) + end + + println("check values at end of run:") + checkvals = [ + ("atm", "pO2PAL", 0.995817, 1e-4), + ("atm", "O2_sms", 1000.0, 2.0), + ] + for (domname, varname, checkval, rtol) in checkvals + outputval = PB.get_data(paleorun.output, domname*"."*varname)[end] + println(" check $domname.$varname $outputval $checkval $rtol") + @test isapprox(outputval, checkval, rtol=rtol) + end + +end + +!("P_O2" in skipped_testsets) && @testset "P_O2" begin + + model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_P_O2"; + modelpars=Dict( + "matdir"=>matdir, + ) + ) + tspan=(0,1e5) + + initial_state, modeldata = PALEOmodel.initialize!(model) + paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + + PALEOmodel.ODE.integrateForwardDiff( + paleorun, initial_state, modeldata, tspan, solvekwargs=(reltol=1e-9,) + ) + + println("conservation checks:") + conschecks = [ + ("global", "total_O2", 1e-6 ), + ("ocean", "P_total", 1e-6 ), + ] + for (domname, varname, rtol) in conschecks + startval, endval = PB.get_data(paleorun.output, domname*"."*varname)[[1, end]] + println(" check $domname.$varname $startval $endval $rtol") + @test isapprox(startval, endval, rtol=rtol) + end + + println("check values at end of run:") + checkvals = [ + ("atm", "pO2PAL", 1.00155, 1e-4), + ("atm", "O2_sms", 1000.0, 2.0), + ] + for (domname, varname, checkval, rtol) in checkvals + outputval = PB.get_data(paleorun.output, domname*"."*varname)[end] + println(" check $domname.$varname $outputval $checkval $rtol") + @test isapprox(outputval, checkval, rtol=rtol) + end + +end + +!("P_O2_S_Carb_open" in skipped_testsets) && @testset "P_O2_S_Carb_open" begin + + model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_P_O2_S_Carb_open"; + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "romaniello2010_transport"), # assume zip file has been downloaded and unpacked in this subfolder + # Option 1: add pHfree, TAlk to state variables and apply constraint on TAlk + "TAlkStateExplicit"=>true, + # Option 2: TAlk is an implicit variable, pHfree is a state variable + # "TAlkStateExplicit"=>false, + ), + ) + + # additional configuration to set shelf and deep carbonate burial + config_ocean_romglb_expts( + model, + [ + "set_carbonate_config_79box", # configure shelf and deep carbonate burial for modern Earth + ], + ) + + tspan=(0,1e5) + + initial_state, modeldata = PALEOmodel.initialize!(model) + run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + + PALEOmodel.ODE.integrateDAEForwardDiff( + run, initial_state, modeldata, tspan, + solvekwargs=(reltol=1e-6,) + # solvekwargs=(reltol=1e-5,) + ) + + ncoutfile = tempname(cleanup=true) + PALEOmodel.OutputWriters.save_netcdf(paleorun.output, ncoutfile; check_ext=false) + ncoutput = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), ncoutfile; check_ext=false) + + for output in (paleorun.output, ncoutput) + + println("check values at end of run:") + checkvals = [ + ("fluxRtoOcean", "flux_P", 4.17847e10, 1e-4), # restoring flux + ("atm", "O2_restoring", -2.6115e12, 1e-4), + ("fluxOceanBurial", "flux_total_P", 4.17847e10, 1e-4), + ("fluxOceanBurial", "flux_total_Corg", 2.61155e12, 1e-4), + ("fluxOceanBurial", "flux_total_Ccarb", 2.3240e13, 1e-4), + ] + for (domname, varname, checkval, rtol) in checkvals + outputval = PB.get_data(output, domname*"."*varname)[end] + println(" check $domname.$varname $outputval $checkval $rtol") + @test isapprox(outputval, checkval, rtol=rtol) + end + end + +end + +!("P_O2_S_Carb_open_implicit" in skipped_testsets) && @testset "P_O2_S_Carb_open_implicit" begin + + model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_P_O2_S_Carb_open"; + modelpars=Dict( + "matdir"=>joinpath(@__DIR__, "romaniello2010_transport"), # assume zip file has been downloaded and unpacked in this subfolder + # Option 1: add pHfree, TAlk to state variables and apply constraint on TAlk + # "TAlkStateExplicit"=>true, + # Option 2: TAlk is an implicit variable, pHfree is a state variable + "TAlkStateExplicit"=>false, + ), + ) + + # additional configuration to set shelf and deep carbonate burial + config_ocean_romglb_expts( + model, + [ + "set_carbonate_config_79box", # configure shelf and deep carbonate burial for modern Earth + ], + ) + + tspan=(0,1e5) + + initial_state, modeldata = PALEOmodel.initialize!(model) + run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + + PALEOmodel.ODE.integrateDAEForwardDiff( + run, initial_state, modeldata, tspan, solvekwargs=(reltol=1e-5,) + ) + + println("check values at end of run:") + checkvals = [ + ("fluxRtoOcean", "flux_P", 4.17847e10, 1e-4), # restoring flux + ("atm", "O2_restoring", -2.6115e12, 1e-4), + ("fluxOceanBurial", "flux_total_P", 4.17847e10, 1e-4), + ("fluxOceanBurial", "flux_total_Corg", 2.61155e12, 1e-4), + ("fluxOceanBurial", "flux_total_Ccarb", 2.3240e13, 1e-4), + ] + for (domname, varname, checkval, rtol) in checkvals + outputval = PB.get_data(run.output, domname*"."*varname)[end] + println(" check $domname.$varname $outputval $checkval $rtol") + @test isapprox(outputval, checkval, rtol=rtol) + end + +end + + +end diff --git a/src/ocean/Ocean.jl b/src/ocean/Ocean.jl index 678a725..580d277 100644 --- a/src/ocean/Ocean.jl +++ b/src/ocean/Ocean.jl @@ -14,6 +14,8 @@ include("OceanTransportTMM.jl") include("OceanTransportColumn.jl") +include("OceanTransportRomaniello2010.jl") + include("BioProd.jl") include("VerticalTransport.jl") diff --git a/src/ocean/OceanTransportRomaniello2010.jl b/src/ocean/OceanTransportRomaniello2010.jl new file mode 100644 index 0000000..9833dc4 --- /dev/null +++ b/src/ocean/OceanTransportRomaniello2010.jl @@ -0,0 +1,489 @@ +module OceanTransportRomaniello2010 + +import SparseArrays +import MAT # Matlab file access + +import PALEOboxes as PB +using PALEOboxes.DocStrings +import ..Ocean + +import Infiltrator # Julia debugger + +""" + ReactionOceanTransportRomanielloICBM + +Modern global and Black Sea ocean transport from [Romaniello2010a](@cite) [Romaniello2010](@cite), + +# Global configurations (`circname = Global_79_Box`, `circname = Global_13_Box`): + +The `ocean` Domain consists of three columns, :hlat, :gyre, :upw + +Ocean box indices are: + + :hlat :gyre :upw + --------------------------- + | 1 | 14 | 47 | + | | | | + | | | | + | | | | + | 13 | 46 | 79 | + ------------------- + + :hlat :gyre :upw + --------------------------- + | 1 | 4 | 9 | + | | | | + | | | | + | | | | + | 3 | 8 | 13 | + ------------------- + +The `oceansurface` Domain has three boxes, the `oceanfloor` Domain has 1 box per ocean box. + +# Black Sea configuration (`circname = Black_Sea`) + +The `ocean` Domain consists of a single column. The `oceansurface` Domain has one box, the `oceanfloor` Domain has 1 box per ocean box. + +Bosphorus outflow flux (Parameter `bosph_outflow`) is represented by Variables with default linking to a `fluxBosphorusOutflow` Domain: + + bosph_outflow_ --> fluxBosphorusOutflow.flux_ + +The .yaml configuration file can be used to configure a closed system by adding outflow flux back into the ocean surface: + + variable_links: + bosph_outflow_*: ocean.oceansurface.*_sms + +Bosphorus inflow (plume) flux is represented by Variables with default linking to concentrations and source - sink fluxes in a `BosphorusInflow` Domain: + + bosph_inflow__conc --> BosphorusInflow._conc + bosph_inflow__sms --> BosphorusInflow._sms # will be -ve for flux into Black Sea + +the .yaml configuration file can be used to configure a closed system by sourcing inflow flux from ocean surface: + + variable_links: + bosph_inflow_*_conc: ocean.oceansurface.*_conc + bosph_inflow_*_sms: ocean.oceansurface.*_sms + + +# Implementation +Reads Matlab .mat files created from the published SI with Matlab commands: + + >> circ_global_79_box = Global_79_Box_ICBM_Params + >> circ_global_13_box = Global_13_Box_ICBM_Params + >> circ_black_sea = Black_Sea_ICBM_Params + >> save('romaniello_global79','-struct', 'circ_global_79_box', '-v6') + >> save('romaniello_global13','-struct', 'circ_global_13_box', '-v6') + >> save('romaniello_blacksea','-struct', 'circ_black_sea', '-v6') + +# Parameters +$(PARS) + +# Methods and Variables +$(METHODS_SETUP) +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionOceanTransportRomanielloICBM{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + PB.ParString("matdir", "romaniello2010_transport", + description="folder with Romaniello (2010) transport and geometry data files"), + PB.ParString("circname", "Global_79_Box", allowed_values=["Global_79_Box", "Global_13_Box", "Black_Sea"], + description="transport matrix from Romaniello(2010) G3"), + + PB.ParBool("set_temp", true, + description="true to provide and set ocean.temp variable to (very approximate) values"), + + PB.ParBool("temp_trackglobal", false, + description="track global temperature (apply offset of global temp -15C"), + + PB.ParDouble("bosph_outflow", 6.04e11, units="m^3 yr^-1", + description="Bosphorus outflow (Black_Sea only)"), + + PB.ParDouble("bosph_inflow", 3.05e11, units="m^3 yr^-1", + description="Bosphorus inflow (plume, Black_Sea only)"), + ) + + + temp_oceanC::Vector{Float64} = Float64[] # ocean temperature (degrees C) + + grid_ocean = nothing + grid_oceansurface = nothing + grid_oceanfloor = nothing + + + "Transport matrix: Units yr^{-1} + so dc/dt = trspt_dtm * c [yr^-1] + where c is column vector of tracer concentrations" + trspt_dtm::SparseArrays.SparseMatrixCSC{Float64,Int64} = SparseArrays.spzeros(0, 0) + + "Transpose of trspt_dtm" + trspt_dtm_tr::SparseArrays.SparseMatrixCSC{Float64,Int64} = SparseArrays.spzeros(0, 0) + + "Black_Sea only: fraction of Bosphorus inflow (plume) flux to each box" + plume_frac::SparseArrays.SparseVector{Float64, Int64} = SparseArrays.spzeros(0) + + rom_data = nothing # raw data from .mat file +end + + + +function read_datafiles(rj::ReactionOceanTransportRomanielloICBM) + + matdir = rj.pars.matdir[] # directory containing .mat files + + # rom keys: + # rom[] indices of boxes for this column, ordered surface to floor + # rom["depths"] mid depths (m) cells + # rom[_bnd] indices of horizontal surfaces for this column (length of column + 1) + # rom["bnd_depths"][rom[_bnd]] depths (m) of horizontal surfaces (length of column + 1) + # rom["Hyps"][rom[_bnd]] area (m^2) of horizontal surfaces (length of column + 1) + # rom["T"] transport matrix + + if rj.pars.circname[] == "Global_13_Box" + matfilename = joinpath(matdir, "romaniello_global13.mat") + @info "read_datafiles: $(fullname(rj)) reading transport matrix from file $(matfilename)" + rj.rom_data = MAT.matread(matfilename) + rj.grid_ocean = PB.Grids.UnstructuredColumnGrid( + # domain=rj.domain, + ncells=rj.rom_data["nboxes"], + columnnames=[:hlat, :gyre, :upw], + Icolumns=[Int.(vec(rj.rom_data["hlat"])), Int.(vec(rj.rom_data["gyre"])), Int.(vec(rj.rom_data["upw"]))]) + elseif rj.pars.circname[] == "Global_79_Box" + matfilename = joinpath(matdir, "romaniello_global79.mat") + @info "read_datafiles: $(PB.fullname(rj)) reading transport matrix from file $(matfilename)" + rj.rom_data = MAT.matread(matfilename) + rj.grid_ocean = PB.Grids.UnstructuredColumnGrid( + # domain=rj.domain, + ncells=rj.rom_data["nboxes"], + columnnames=[:hlat, :gyre, :upw], + Icolumns=[Int.(vec(rj.rom_data["hlat"])), Int.(vec(rj.rom_data["gyre"])), Int.(vec(rj.rom_data["upw"]))]) + elseif rj.pars.circname[] == "Black_Sea" + matfilename = joinpath(matdir, "romaniello_blacksea.mat") + @info "read_datafiles: $(PB.fullname(rj)) reading transport matrix from file $(matfilename)" + rj.rom_data = MAT.matread(matfilename) + # single column, no name + rj.grid_ocean = PB.Grids.UnstructuredColumnGrid( + # domain=rj.domain, + ncells=rj.rom_data["nboxes"], + columnnames=[:-], + Icolumns=[collect(1:rj.rom_data["nboxes"])]) + else + error("unknown circname '$(rj.pars.circname[])") + end + + return nothing +end + + +function PB.set_model_geometry(rj::ReactionOceanTransportRomanielloICBM, model::PB.Model) + + read_datafiles(rj) + + ocean_cells = rj.grid_ocean.ncells + + # set subdomain mappings + isurf=[col[1] for col in rj.grid_ocean.Icolumns] # first index in each column is surface box + ifloor=collect(1:ocean_cells) # every box has an oceanfloor box under it + + PB.Grids.set_subdomain!(rj.grid_ocean, "oceansurface", PB.Grids.BoundarySubdomain(isurf), true) + @info " set ocean.oceansurface Subdomain size=$(length(isurf))" + PB.Grids.set_subdomain!(rj.grid_ocean, "oceanfloor", PB.Grids.BoundarySubdomain(ifloor), true) + @info " set ocean.oceanfloor Subdomain size=$(length(ifloor))" + + rj.grid_oceansurface = PB.Grids.UnstructuredVectorGrid( + # domain=PB.get_domain(model, "oceansurface"), + ncells=length(isurf), + cellnames=Dict(rj.grid_ocean.columnnames[i]=>i for i in 1:length(rj.grid_ocean.columnnames))) + PB.Grids.set_subdomain!(rj.grid_oceansurface, "ocean", PB.Grids.InteriorSubdomain(ocean_cells, isurf), true) + + rj.grid_oceanfloor = PB.Grids.UnstructuredColumnGrid( + # domain=PB.get_domain(model, "oceanfloor"), + ncells=length(ifloor), + columnnames=rj.grid_ocean.columnnames, + Icolumns=rj.grid_ocean.Icolumns) + PB.Grids.set_subdomain!(rj.grid_oceanfloor, "ocean", PB.Grids.InteriorSubdomain(ocean_cells, ifloor), true) + + Ocean.set_model_domains(model, rj.grid_ocean, rj.grid_oceansurface, rj.grid_oceanfloor) + + return nothing +end + +function PB.register_methods!(rj::ReactionOceanTransportRomanielloICBM) + + physvars = [ + PB.VarProp("sal", "psu", "Ocean salinity"), + PB.VarProp("rho", "kg m^-3", "physical ocean density"), + # no length check as create and set oceansurface Variable from atm Domain + PB.VarProp("oceansurface.open_area_fraction","","fraction of area open to atmosphere", attributes=(:check_length=>false,)), + ] + + PB.add_method_setup!( + rj, + do_setup_grid, + ( + PB.VarList_namedtuple(Ocean.grid_vars_all), + PB.VarList_namedtuple(physvars), + ), + ) + + if rj.pars.set_temp[] + tempvars = [ + PB.VarDepScalar("(global.TEMP)", "K", "global mean temperature"), + PB.VarProp("temp", "Kelvin", "Ocean temperature"), + ] + + PB.add_method_do!( + rj, + do_temperature, + (PB.VarList_namedtuple(tempvars), ), + ) + end + + return nothing +end + + +function do_setup_grid( + m::PB.ReactionMethod, + (grid_vars, physvars), + cellrange::PB.AbstractCellRange, + attribute_name +) + attribute_name == :setup || return + + rj = m.reaction + # rom keys: + # rom[] indices of boxes for this column, ordered surface to floor + # rom["depths"] mid depths (m) cells + # rom[_bnd] indices of horizontal surfaces for this column (length of column + 1) + # rom["bnd_depths"][rom[_bnd]] depths (m) of horizontal surfaces (length of column + 1) + # rom["Hyps"][rom[_bnd]] area (m^2) of horizontal surfaces (length of column + 1) + # rom["T"] transport matrix + + for icl in 1:length(rj.grid_ocean.Icolumns) + iboxes = rj.grid_ocean.Icolumns[icl] # indices for this column + if rj.pars.circname[] == "Black_Sea" + ibnd = collect(1:length(iboxes)+1) + else + ibnd = Int.(rj.rom_data[String(rj.grid_ocean.columnnames[icl])*"_bnd"]) + end + + grid_vars.volume[iboxes] .= rj.rom_data["V"][iboxes]/1000.0 # convert to m^3 + + grid_vars.Abox[iboxes] .= rj.rom_data["Hyps"][ibnd[1:end-1]] + grid_vars.Asurf[icl] = rj.rom_data["Hyps"][ibnd[1]] + # assume column closed at bottom box (Romaniello hypsometry has small term here) + grid_vars.Afloor[iboxes[1:end-1]] .= grid_vars.Abox[iboxes[1:end-1]] - grid_vars.Abox[iboxes[2:end]] + grid_vars.Afloor[iboxes[end]] = grid_vars.Abox[iboxes[end]] + + # NB: there is an inconsistency in Matlab (rj.rom_data) vs PALEO conventions: volume != Abox * (zupper - zlower) + grid_vars.zupper[iboxes] .= -rj.rom_data["bnd_depths"][ibnd[1:end-1]] + grid_vars.zmid[iboxes] .= -rj.rom_data["depths"][iboxes] + # grid_vars.zlower[iboxes] .= grid_vars.zupper[iboxes]-grid_vars.volume[iboxes]./grid_vars.Abox[iboxes] + grid_vars.zlower[iboxes] .= -rj.rom_data["bnd_depths"][ibnd[2:end]] + + grid_vars.zfloor[iboxes] .= grid_vars.zlower[iboxes] + # m^2 (k_nbox,k_nbox) Area(i, j) of horizontal contact between lower surface of box i and box j + # rj.grid_vars.Aoverlap + end + grid_vars.volume_total[] = sum(grid_vars.volume) + grid_vars.Afloor_total[] = sum(grid_vars.Afloor) + + # attach coordinates to grid for output visualisation etc + empty!(rj.domain.grid.z_coords) + push!(rj.domain.grid.z_coords, PB.FixedCoord("zmid", grid_vars.zmid, PB.get_variable(m, "zmid").attributes)) + push!(rj.domain.grid.z_coords, PB.FixedCoord("zlower", grid_vars.zlower, PB.get_variable(m, "zlower").attributes)) + push!(rj.domain.grid.z_coords, PB.FixedCoord("zupper", grid_vars.zupper, PB.get_variable(m, "zupper").attributes)) + + # constant density + grid_vars.rho_ref .= 1027 + + grid_vars.pressure .= -grid_vars.zmid # pressure(dbar) ~ depth (m) + + # read transport matrix + rj.trspt_dtm = copy(rj.rom_data["T"]) + + + if rj.pars.circname[] == "Black_Sea" + + # Add Bosphorus outflow back into box 1 so trspt_dtm represents a closed system + # (the Matlab model included this flux in the T matrix - we will explicitly add it in separately if needed) + # yr-1 m^3 yr-1 m^3 + rj.trspt_dtm[1, 1] += rj.rom_data["Bosph_Outflow"]/grid_vars.volume[1] + + # Calculate the distribution of Bosphorus inflow plume to each box + # (the Matlab model supplies this as Tplume, a transport matrix) + # We want plume_frac, the fraction of inflow flux to each box, where sum(plume_frac) == 1.0 + # m^3 yr-1 = m^3 * yr-1 + plume_flux = sum(grid_vars.volume.*rj.rom_data["Tplume"][:, 1]) + # unitless = m^3 * yr-1 / m3 yr-1 + rj.plume_frac = grid_vars.volume.*rj.rom_data["Tplume"][:, 1] ./ plume_flux + + # check conservation + # for i in 1:rj.grid_ocean.ncells + # println("cell ", i, " sum of flux with conc in cell i = 1 mol m-3: ", sum(rj.trspt_dtm[:, i].*grid_vars.volume)) + # end + end + + rj.trspt_dtm_tr = SparseArrays.sparse(transpose(rj.trspt_dtm)) + + # set salinity + physvars.sal .= 35.0 + # constant density + physvars.rho .= 1027 + + physvars.open_area_fraction .= 1.0 + + # optionally set a very approximate default ocean temperature + if rj.pars.set_temp[] + rj.temp_oceanC = Vector{Float64}(undef, rj.grid_ocean.ncells) + rj.temp_oceanC .= 2.0 # all interior boxes at 2C + if rj.pars.circname[] == "Black_Sea" + rj.temp_oceanC[rj.domain.grid.subdomains["oceansurface"].indices] .= 7.5 # guessed 7.5C surface temp + else + rj.temp_oceanC[rj.domain.grid.subdomains["oceansurface"].indices] .= [2.5, 25.0, 25.0] # guessed + end + end + + return nothing +end + + + +function PB.register_dynamic_methods!(rj::ReactionOceanTransportRomanielloICBM) + + (transport_conc_vars, transport_sms_vars, transport_input_vars) = + Ocean.find_transport_vars(rj.domain, add_transport_input_vars=true) + + PB.add_method_do!( + rj, + do_transport, + ( + PB.VarList_namedtuple(PB.VarDep.(Ocean.grid_vars_ocean)), + PB.VarList_components(transport_conc_vars), + PB.VarList_components(transport_sms_vars), + PB.VarList_components(transport_input_vars), + ), + preparefn=Ocean.prepare_transport + ) + + if rj.pars.circname[] == "Black_Sea" + tracernames = [v.localname[1:end-5] for v in transport_conc_vars] + + bosph_outflow_vars = [ + PB.VarContribScalar("bosph_outflow_$n"=>"fluxBosphorusOutflow.flux_$n", "mol yr-1", "Bosphorus outflow flux") + for n in tracernames + ] + PB.add_method_do!( + rj, + do_bosph_outflow, + ( + PB.VarList_namedtuple(PB.VarDep.(Ocean.grid_vars_ocean)), + PB.VarList_components(transport_conc_vars), + PB.VarList_components(transport_sms_vars), + PB.VarList_components(bosph_outflow_vars), + ), + ) + + bosph_inflow_conc_vars = [ + PB.VarDepScalar("bosph_inflow_$(n)_conc"=>"BosphorusInflow.$(n)_conc", "mol m-3", "Bosphorus inflow concentration") + for n in tracernames + ] + bosph_inflow_sms_vars = [ + PB.VarContribScalar("bosph_inflow_$(n)_sms"=>"BosphorusInflow.$(n)_sms", "mol yr-1", "Bosphorus inflow source - sink (-ve for flux to Black Sea)") + for n in tracernames + ] + PB.add_method_do!( + rj, + do_bosph_inflow, + ( + PB.VarList_components(transport_sms_vars), + PB.VarList_components(bosph_inflow_conc_vars), + PB.VarList_components(bosph_inflow_sms_vars), + ), + ) + end + + return nothing +end + + +function do_transport( + m::PB.ReactionMethod, + (grid_vars, transport_conc_components, transport_sms_components, transport_input_components, buffer), + cellrange::PB.AbstractCellRange, + deltat +) + rj = m.reaction + + Ocean.do_transport_tr( + grid_vars, transport_conc_components, transport_sms_components, transport_input_components, buffer, + rj.trspt_dtm_tr, + cellrange) + return nothing +end + +function do_bosph_outflow( + m::PB.ReactionMethod, + pars, + (grid_vars, transport_conc_components, transport_sms_components, bosph_outflow_components), + cellrange::PB.AbstractCellRange, + deltat +) + + errmsg="do_bosph_outflow: components length mismatch transport_conc_components, transport_sms_components, bosph_outflow_components (check :field_data (ScalarData, IsotopeLinear etc) match)" + for (conc, sms, outflow) in PB.IteratorUtils.zipstrict( + transport_conc_components, transport_sms_components, bosph_outflow_components; errmsg=errmsg + ) + for i in cellrange.indices + if i == 1 # surface box + # mol yr-1 = m^3 yr-1 * mol m-3 + flux = pars.bosph_outflow[]*conc[1] + sms[1] -= flux + outflow[] += flux + end + end + end + + return nothing +end + +function do_bosph_inflow( + m::PB.ReactionMethod, + pars, + (transport_sms_components, bosph_inflow_conc_components, bosph_inflow_sms_components), + cellrange::PB.AbstractCellRange, + deltat +) + rj = m.reaction + + errmsg = "do_bosph_inflow: components length mismatch transport_sms_components, bosph_inflow_conc_components, bosph_inflow_sms_components (check :field_data (ScalarData, IsotopeLinear etc) match)" + for (sms, inflow_conc, inflow_sms) in PB.IteratorUtils.zipstrict( + transport_sms_components, bosph_inflow_conc_components, bosph_inflow_sms_components; errmsg=errmsg + ) + for i in cellrange.indices + # mol yr-1 = m^3 yr-1 * mol m-3 + flux = pars.bosph_inflow[]*inflow_conc[]*rj.plume_frac[i] + sms[i] += flux + inflow_sms[] -= flux + end + end + + return nothing +end + +function do_temperature(m::PB.ReactionMethod, pars, (tempvars, ), cellrange::PB.AbstractCellRange, deltat) + rj = m.reaction + + # Set temperature + if pars.temp_trackglobal[] + tempvars.temp .= rj.temp_oceanC .- 15.0 .+ tempvars.TEMP[] # temperature (K) + else + tempvars.temp .= rj.temp_oceanC .+ PB.Constants.k_CtoK # temperature (K) + end + + return nothing +end + +end # module