Merge pull request #4 from OceanBioME/jsw/add-fig-to-readme

Update README.md

README.md

@@ -5,6 +5,19 @@ Currently, it contains a script for a repeat year forcing/free runing near globa
 
 **Please note** this is a work in progress and the scripts were written as demonstrations - the physics and biogeochemistry are not complete or calibrated.
 
+https://github.com/OceanBioME/GlobalOceanBioME/assets/26657828/ca2dd956-5fa5-4bb3-b564-743cbd379bd7
+
+<details>
+  <summary>More detailed description of proof of concept above</summary>
+  This proof of concept is a 1°x1° near-global run between 75°N and 75°S with 48 vertical levels.
+
+  The physics includes isopycnal skew-symmetric diffusivity and a convective adjustment boundary layer closure, and the biogeochemistry is a simple 4-variable nutrient, phytoplankton, zooplankton, detritus model (<a href="https://oceanbiome.github.io/OceanBioME.jl/stable/model_components/biogeochemical/NPZ/">detailed here</a>). 
+
+  The physics is initialised with temperature and salinity, and the top boundary is forced with data, from ECCO v4 (as per <a href="https://github.com/CliMA/ClimaOcean.jl">ClimaOcean NearGlobalSimulation</a>), but the biogeochemistry in this simple example is initialised to be uniform across the globe (hence the abrupt change at the start). The surface photosynthetically available radiation for the biogeochemical model is interpolated from NASA ocean colour data.
+
+  This model runs in around 45 minutes per year to run on an Nvidia A100 GPU.
+</details>
+
 ## To run the proof of concept script
 
 ```

scripts/one_degree_proof_of_concept.jl

@@ -53,13 +53,13 @@ simulation = one_degree_near_global_simulation(architecture;
     interior_background_vertical_viscosity = 1e-4,
     surface_background_vertical_viscosity = 1e-4,
     biogeochemistry,
-    biogeochemistry_kwargs = (surface_photosynthetically_active_radiation = OneDegreeSurfacePAR(architecture), scale_negatives = true),
+    biogeochemistry_kwargs = (surface_photosynthetically_active_radiation = OneDegreeSurfacePAR(architecture), scale_negatives = true, invalid_fill_value = 0),
     tracers = (:N, :P, :Z, :D, :T, :S) # have to specify since NPZD adds T but not S, and buoyancy requires both
 )
 
 # Define output
 slices_save_interval = 2day
-fields_save_interval = 1days
+fields_save_interval = 3hour
 Nx, Ny, Nz = size(simulation.model.grid)
 
 dir = "bgc" 
@@ -69,7 +69,7 @@ with_isopycnal_skew_symmetric_diffusivity || (output_prefix *= "_no_gm")
 
 simulation.output_writers[:checkpointer] = Checkpointer(simulation.model; dir,
                                                         prefix = output_prefix * "_checkpointer",
-                                                        schedule = TimeInterval(365days),
+                                                        schedule = TimeInterval(5days),
                                                         cleanup = true,
                                                         overwrite_existing = true)
 
@@ -79,26 +79,39 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocit
                                                       schedule = TimeInterval(fields_save_interval),
                                                       filename = output_prefix * "_fields",
                                                       with_halos = true,
-                                                      overwrite_existing = true)
+                                                      overwrite_existing = true,
+indices = (:, :, model.grid.Nz))
+#=
+""" Load initial conditions from Copernicus models.
 
-set!(model, N = 10.0, P = 0.1, Z = 0.01)
+P and N are a direct downsampling and unit conversion from https://doi.org/10.48670/moi-00015
+Z is downsampled and unit converted from https://doi.org/10.48670/moi-00020 and then divided by 
+the mixed layer depth from https://doi.org/10.48670/moi-00024, and then applied uniformly over
+the mixed region
 
-simulation.Δt = 20minute
-simulation.stop_time = time(simulation) + 1days
+I think the units are wrong still. For now assume zooplankton are 138:106:16:1 O₂:C:N:P redfield ratio.
+"""
+file = jldopen(datadep"2010_near_global_bgc/initial_conditions.jld2")
 
-@info "Running a simulation with Δt = $(prettytime(simulation.Δt)) from $(prettytime(simulation.model.clock.time)) until $(prettytime(simulation.stop_time))"
+N_init = on_architecture(architecture, file["N"])
+P_init = on_architecture(architecture, file["P"].*16/106)# mmolC -> mmolN
+Z_init = on_architecture(architecture, file["Z"].*16/(138*16*2+106*12+16*14+1*31)) # mg -> mmolN
+
+close(file)
+
+set!(model, N = N_init, P = P_init, Z = Z_init)
+=#
+set!(model, N = 10.0, P = 0.1, Z = 0.01)
 
 prog(sim) = @info "$(prettytime(time(sim))) in $(prettytime(sim.run_wall_time)) with Δt = $(prettytime(sim.Δt))"
 
 add_callback!(simulation, prog, IterationInterval(10))
 
-run!(simulation)
-
 simulation.callbacks[:nan_checker] = Callback(Oceananigans.Simulations.NaNChecker(; fields = merge(model.tracers, model.velocities), erroring = true), IterationInterval(10))
 
-simulation.Δt = 10minutes
-simulation.stop_time = start_time + 3 * 365days
+simulation.Δt = 20minute
+simulation.stop_time = 345days + 365days*3
 
 @info "Running a simulation with Δt = $(prettytime(simulation.Δt)) from $(prettytime(simulation.model.clock.time)) until $(prettytime(simulation.stop_time))"
 
-run!(simulation; pickup = false)
+run!(simulation; pickup = false)

scripts/plot_one_degree.jl

@@ -0,0 +1,22 @@
+using Oceananigans, CairoMakie, GeoMakie
+
+P = FieldTimeSeries("bgc/near_global_360_150_48_ConvectiveAdjustmentVerticalDiffusivity_fields.jld2", "P");
+
+n = Observable(1)
+
+P_plt = @lift interior(P[$n], :, :, 1)
+
+fig = Figure()
+
+ax = GeoAxis(fig[1, 1];title, dest = "+proj=natearth2")
+
+xc, yc, zf = nodes(P.grid, Center(), Center(), Face())
+
+hm = heatmap!(ax, xc, yc, P_plt, colorrange = (0.0001, 5), colorscale = log10)
+
+Colorbar(fig[2, 1], hm, vertical= false, label = "Phytoplankton (mmolN/m³)")
+
+CairoMakie.record(fig, "P.mp4", 1:10:length(P.times)) do i; 
+    n[] = i; 
+    @info "$(n.val)"; 
+end