Merge pull request #22 from PALEOtoolkit/add_shelf1D

Add examples/shelf1D

.gitignore

@@ -4,3 +4,6 @@ LocalPreferences.toml
 /docs/src/collated_examples/
 .vscode
 *.jld2
+*.nc
+*.zip
+examples/shelf1D/S2P3_transport_20240614

Project.toml

@@ -9,6 +9,7 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 PALEOaqchem = "673cec3b-17d1-411f-9fcd-71c01c593120"
 PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"

docs/src/PALEOocean_Reactions.md

@@ -10,6 +10,7 @@ OceanNoTransport.ReactionOceanNoTransport
 OceanTransport3box.ReactionOceanTransport3box
 OceanTransport6box.ReactionOceanTransport6box
 OceanTransportTMM.ReactionOceanTransportTMM
+OceanTransportColumn.ReactionOceanTransportColumn
 ```
 
 ## Vertical Transport
@@ -55,4 +56,13 @@ Burial.ReactionBurialEffCarb
 ### Organic carbon and phosphorus burial
 ```@docs
 Burial.ReactionBurialEffCorgP
-```
+```
+
+## Global
+```@meta
+CurrentModule = PALEOocean.Global
+```
+```@docs
+Insolation.ReactionForceInsolationModernEarth
+Insolation.insolMITgcmDIC
+```

docs/src/paleo_references.bib

@@ -1,3 +1,16 @@
+@article{Brock1981,
+author = {Brock, Thomas D},
+doi = {10.1016/0304-3800(81)90011-9},
+journal = {Ecological Modelling},
+month = {nov},
+number = {1-2},
+pages = {1--19},
+title = {{Calculating solar radiation for ecological studies}},
+volume = {14},
+year = {1981}
+}
+
+
 @article{Caldeira1993,
 author = {Caldeira, Ken and Rampino, Michael R},
 doi = {10.1029/93PA01163},

examples/Project.toml

@@ -12,6 +12,7 @@ PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4"
 PALEOcopse = "4a6ed817-0e58-48c6-8452-9e9afc8cb508"
 PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c"
 PALEOocean = "41de04b1-2efd-44ae-92ae-39d71a4fd99b"
+PALEOsediment = "e0a37952-6f01-4236-91ff-62fdc855f67b"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

examples/mitgcm/MITgcm_2deg8_COPDOM.yaml

@@ -203,7 +203,7 @@ PO4MMbase:
         oceansurface:
             reactions:
                 force_par:
-                    class: ReactionForceInsolation
+                    class: ReactionForceInsolationModernEarth
 
                     variable_links:
                         insolation: surface_insol
@@ -618,7 +618,7 @@ PO4MMcarbSCH4:
         oceansurface:
             reactions:
                 force_par:
-                    class: ReactionForceInsolation
+                    class: ReactionForceInsolationModernEarth
 
                     variable_links:
                         insolation: surface_insol            

examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.jl

@@ -16,8 +16,8 @@ model = PB.create_model_from_config(
     joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMcarbSCH4"
 )
 
-# toutputs = [0, 1.0, 10.0] # , 100.0, 1000.0, 1999.5, 2000.0, 2999.5, 3000.0]
-toutputs = [0, 1.0, 10.0, 100.0, 1000.0, 1999.5, 2000.0,]
+toutputs = [0, 1.0, 10.0] # , 100.0, 1000.0, 1999.5, 2000.0, 2999.5, 3000.0]
+# toutputs = [0, 1.0, 10.0, 100.0, 1000.0, 1999.5, 2000.0,]
 
 # start with low oxygen to test marine sulphur system
 PB.set_variable_attribute!(model, "atm", "O2", :initial_value, 0.1*3.71e19)
@@ -28,8 +28,8 @@ tstep = transportMITgcm.pars.Aimp_deltat[]/PB.Constants.k_secpyr
 @info "using tstep=$tstep yr"
 
 # output_filename = "MITgcm_PO4MMcarbSCH42deg8FP64_3000yr_20210210"
-output_filename = "MITgcm_PO4MMcarbSCH42deg8FP64_2000yr_20230422"
-# output_filename = ""
+# output_filename = "MITgcm_PO4MMcarbSCH42deg8FP64_2000yr_20230422"
+output_filename = ""
 
 pickup_output = nothing
 initial_state, modeldata = PALEOmodel.initialize!(model, threadsafe=use_threads, pickup_output=pickup_output)  

examples/mitgcm/MITgcm_2deg8_abiotic.yaml

@@ -104,7 +104,7 @@ abiotic_O2:
         oceansurface:
             reactions:
                 force_par:
-                    class: ReactionForceInsolation
+                    class: ReactionForceInsolationModernEarth
 
                     variable_links:
                         insolation: surface_insol

examples/mitgcm/MITgcm_ECCO_COPDOM.yaml

@@ -171,7 +171,7 @@ PO4MMbase:
         oceansurface:
             reactions:
                 force_par:
-                    class: ReactionForceInsolation
+                    class: ReactionForceInsolationModernEarth
 
                     variable_links:
                         insolation: surface_insol
@@ -550,7 +550,7 @@ PO4MMcarbSCH4:
         oceansurface:
             reactions:
                 force_par:
-                    class: ReactionForceInsolation
+                    class: ReactionForceInsolationModernEarth
 
                     variable_links:
                         insolation: surface_insol            

examples/mitgcm/config_mitgcm_expts.jl

@@ -5,7 +5,6 @@ import PALEOocean
 import PALEOmodel
 using Printf
 
-include("Insolation.jl")
 include("../atmreservoirreaction.jl") # temporary solution to make ReactionReservoirAtm available
 
 

examples/shelf1D/PALEO_examples_shelf1D_O2_only.jl

@@ -0,0 +1,82 @@
+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_shelf1D_expts.jl")
+include("plot_shelf.jl")
+include("../atmreservoirreaction.jl") # ReactionReservoirAtm
+
+transport_dir = "S2P3_transport_20240614" # folder containing S2P3 physical variables output collated to netcdf files
+
+# abiotic atmosphere-ocean, O2 only
+model = PB.create_model_from_config(
+    joinpath(@__DIR__, "PALEO_examples_shelf1D_cfg.yaml"),
+    "shelf1D_abiotic_O2";
+    modelpars=Dict(
+        "phys_file"=>joinpath(transport_dir, "S2P3_depth80_m2amp04_phys.nc"),
+        "surf_file"=>joinpath(transport_dir, "S2P3_depth80_m2amp04_surf.nc"),
+    )
+)
+
+config_ocean_shelf1D_expts(model, ["baseline"]); tspan=(0,2.0)
+
+
+initial_state, modeldata = PALEOmodel.initialize!(model)
+
+# Check initial derivative:
+# initial_deriv = similar(initial_state)
+# PALEOmodel.SolverFunctions.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0)
+# Check Jacobian:
+# jac, jac_prototype = PALEOmodel.JacobianAD.jac_config_ode(:ForwardDiffSparse, model, initial_state, modeldata, 0.0)
+# J = copy(jac_prototype)
+# jac(jac_prototype, initial_state, nothing, 0.0)
+
+paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
+
+# first run includes JIT time
+sol = PALEOmodel.ODE.integrateForwardDiff(
+    paleorun, initial_state, modeldata, tspan,
+    solvekwargs=(saveat=1/8760, reltol=1e-5, maxiters=1000000),
+)  
+
+# sol = PALEOmodel.ODE.integrateDAEForwardDiff(
+#     paleorun, initial_state, modeldata, tspan,
+#     solvekwargs=(saveat=1/8760, reltol=1e-5, maxiters=200000),
+# )
+
+
+########################################
+# Plot output
+########################################
+colT=collect(range(tspan[1], stop=tspan[end], step=0.5))
+
+# individual plots
+# plotlyjs(size=(750, 565))
+# pager = PALEOmodel.DefaultPlotPager()
+
+# assemble plots onto screens with 6 subplots
+gr(size=(1600, 900))
+pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm)))
+
+plot_shelf_phys(paleorun.output; pager=pager)
+pager(:newpage)
+plot_tracers_conc(
+    paleorun.output; 
+    tracers=["Tfast", "Tslow", "O2"],
+    colT=colT,
+    plot_totals=true,
+    pager=pager
+)
+plot_airsea(paleorun.output; tracers=["O2"], pager=pager)
+
+pager(:newpage) # flush output
+

examples/shelf1D/PALEO_examples_shelf1D_P_O2.jl

@@ -0,0 +1,78 @@
+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_shelf1D_expts.jl")
+include("plot_shelf.jl")
+include("../atmreservoirreaction.jl") # ReactionReservoirAtm
+
+transport_dir = "S2P3_transport_20240614" # folder containing S2P3 physical variables output collated to netcdf files
+
+# P, O2 only population-based phytoplankton model
+model = PB.create_model_from_config(
+    joinpath(@__DIR__, "PALEO_examples_shelf1D_cfg.yaml"),
+    "shelf1D_P_O2";
+    modelpars=Dict(
+        "phys_file"=>joinpath(transport_dir, "S2P3_depth80_m2amp04_phys.nc"),
+        "surf_file"=>joinpath(transport_dir, "S2P3_depth80_m2amp04_surf.nc"),
+    )
+)
+
+# config_ocean_shelf1D_expts(model, ["QE"]); tspan=(0,5.0)
+config_ocean_shelf1D_expts(model, []); tspan=(0,5.0)
+
+initial_state, modeldata = PALEOmodel.initialize!(model)
+
+# Check initial derivative:
+# initial_deriv = similar(initial_state)
+# PALEOmodel.SolverFunctions.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0)
+# Check Jacobian:
+# jac, jac_prototype = PALEOmodel.JacobianAD.jac_config_ode(:ForwardDiffSparse, model, initial_state, modeldata, 0.0)
+# J = copy(jac_prototype)
+# jac(jac_prototype, initial_state, nothing, 0.0)
+
+paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
+
+# first run includes JIT time
+sol = PALEOmodel.ODE.integrateForwardDiff(
+    paleorun, initial_state, modeldata, tspan,
+    solvekwargs=(saveat=1/8760, reltol=1e-5, maxiters=1000000),
+)  
+
+
+########################################
+# Plot output
+########################################
+
+colT=collect(range(tspan[1], stop=tspan[end], step=0.5))
+
+# individual plots
+# plotlyjs(size=(750, 565))
+# pager = PALEOmodel.DefaultPlotPager()
+
+# assemble plots onto screens with 6 subplots
+gr(size=(1600, 900))
+pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm)))
+
+plot_shelf_phys(paleorun.output; pager=pager)
+pager(:newpage)
+plot_tracers_conc(
+    paleorun.output; 
+    tracers=["O2"],
+    colT=colT,
+    plot_totals=true,
+    pager=pager
+)
+plot_airsea(paleorun.output; tracers=["O2"], pager=pager)
+plot_biota(paleorun.output, colT=colT, pager=pager)
+pager(:newpage) # flush output
+
+