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