Skip to content

Commit

Permalink
Merge pull request #31 from PALEOtoolkit/napc_tutorial
Browse files Browse the repository at this point in the history 
Check-in of examples/NAPC_tutorial
  • Loading branch information
@ssjos
ssjos authored Jun 17, 2024
2 parents aa02926 + af3ee2e commit 637deaf
Show file tree
Hide file tree
Showing 5 changed files with 558 additions and 0 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,5 @@ LocalPreferences.toml
/docs/build/
/docs/src/collated_examples/
.vscode
.nc
examples/NAPC_tutorial/S2P3_transport_20240614
296 changes: 296 additions & 0 deletions examples/NAPC_tutorial/NAPC_cfg.yaml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,296 @@
########################################################
# NP model from Armstrong 1994 in 1D shelf setting
########################################################

NP_shelf:
parameters:
CIsotope: ScalarData
SIsotope: ScalarData
phys_file: S2P3_depth80_m2amp04_phys.nc
surf_file: S2P3_depth80_m2amp04_surf.nc
domains:
global:
# scalar domain

reactions:
total_O2:
class: ReactionSum
parameters:
vars_to_add: [atm.O2, ocean.O2_total, -106*ocean.P0_total, -106*ocean.POP_total]
variable_links:
sum: total_O2

total_P:
class: ReactionSum
parameters:
vars_to_add: [ocean.P_total, ocean.POP_total, ocean.P0_total]
variable_links:
sum: total_P

fluxAtmtoOceansurface:
reactions:
fluxtarget:
class: ReactionFluxTarget
parameters:
flux_totals: true
fluxlist: ["O2"]

fluxOceanfloor:
reactions:
particulatefluxtarget:
class: ReactionFluxTarget
parameters:
flux_totals: true
target_prefix: particulateflux_
fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] # fluxlist_BioParticulate

solutefluxtarget:
class: ReactionFluxTarget
parameters:
flux_totals: true
target_prefix: soluteflux_
fluxlist: ["DIC::CIsotope", "TAlk", "P", "O2"] #, "SO4::SIsotope", "H2S::SIsotope", "CH4::CIsotope"] # fluxlist_Solute

transferPOP:
class: ReactionFluxToComponents
parameters:
outputflux_prefix: particulateflux_
outputflux_names: ["Corg", "N", "P"]
outputflux_stoich: [106.0, 0.0, 1.0] # must match bioprod stoich
variable_links:
inputflux: sinkflux_POP

atm:
reactions:
reservoir_O2:
# NB: rescaled to 1e-14 * global atmosphere, to workaround loss of numerical precision when tracking
# exchanges with a 1 m^2 ocean column
class: ReactionReservoirAtm
parameters:
moles1atm: 1.77e6 # mol 1e-14*1.77e20
variable_links:
R*: O2*
pRatm: pO2atm
pRnorm: pO2PAL
variable_attributes:
R:norm_value: 3.7e5 # mol 1e-14 * present-day global atmospheric level
R:initial_value: 3.7e5 # mol for 1e-14 * global atmosphere

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
piston_fixed: false

insol:
class: ReactionForceInsolationModernEarth
parameters:
latitude: [50.0] # degrees N
variable_links:
insolation: surface_insol

open_area_fraction:
class: ReactionConst
parameters:
constnames: ["open_area_fraction"]
variable_attributes:
open_area_fraction%initial_value: 1.0

wind_speed:
class: ReactionForceGrid
parameters:
netcdf_file: external%surf_file
data_var: wspeedms
time_var: time
tidx_start: 1
tidx_end: 365
cycle_time: 1.0 # yr (periodic annual forcing)
variable_links:
F: wind_speed

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:
transport1D:
class: ReactionOceanTransportColumn
parameters:
grid_file: external%phys_file

sal:
class: ReactionConst
parameters:
constnames: ["sal"]
variable_attributes:
sal%initial_value: 35.0

temp:
class: ReactionForceGrid
parameters:
netcdf_file: external%phys_file
data_var: temp
constant_offset: 273.15 # convert deg C in netcdf file to K
time_var: time
tidx_start: 1
tidx_end: 8760
cycle_time: 1.0 # yr (periodic annual forcing)
variable_links:
F: temp # Kelvin, Ocean temperature

rho:
class: ReactionForceGrid
parameters:
netcdf_file: external%phys_file
data_var: sigmat
constant_offset: 1000.0 # convert to kg m-3 ocean density
time_var: time
tidx_start: 1
tidx_end: 8760
cycle_time: 1.0 # yr (periodic annual forcing)
variable_links:
F: rho # kg m-3, ocean density

Kz:
class: ReactionForceGrid
parameters:
netcdf_file: external%phys_file
data_var: Kz
time_var: time
tidx_start: 1
tidx_end: 8760
cycle_time: 1.0 # yr (periodic annual forcing)
variable_links:
F: Kz # m^2 s-1 eddy diffusivity on upper cell surfaces

sinkfloat:
class: ReactionSinkFloat
parameters:
transportfloor: true

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.1 # m-1

reservoir_P0:
class: ReactionReservoirTotal
parameters:
variable_links:
R*: P0*
variable_attributes:
R:norm_value: 0.001
R:initial_value: 0.001
R_conc:vertical_movement: 0.0 # m d-1
R_conc:advect: true
R_conc:specific_light_extinction: 0.0 # m^2 mol-1

P0_growth:
class: ReactionPhytoplankton
parameters:
mumax: 1.4 #units="d-1" max growth rate
Ks: 1e-3 #units="mol m-3 nutrient limitation factor
insol_scale: 0.005 #units="" insolation scale factor
O2_thresh: 1e-6 #units="mol m-3" oxygen threshold constant
stoich_PtoO2: -106.0 #units="" stoichiometry, negative by convention
variable_links:
PX*: P0*

P0_decay:
class: ReactionParticleDecay
parameters:
decay_timescale: 0.17 # yr
variable_links:
Particle*: P0*
decayflux: POP_sms

reservoir_POP:
class: ReactionReservoirTotal
variable_links:
R*: POP*
variable_attributes:
R:initial_value: 0.0 # concentration m-3
R_conc:vertical_movement: -100.0 # m d-1
R_conc:advect: true

POP_decay:
class: ReactionParticleDecay
parameters:
decay_timescale: 0.05 # yr
variable_links:
Particle*: POP*
decayflux: POP_decay

POP_decaycomponents:
class: ReactionFluxToComponents
parameters:
outputflux_prefix: remin_
outputflux_names: ["Corg", "N", "P"]
outputflux_stoich: [106.0, 0.0, 1.0] # must match bioprod stoich
variable_links:
inputflux: POP_decay

reminocean:
class: ReactionReminO2
parameters:
variable_links:
soluteflux_P: P_sms
soluteflux_O2: O2_sms


oceanfloor:
reactions:
reminoceanfloor:
class: ReactionReminO2
parameters:
variable_links:
remin*: particulateflux*
soluteflux_P: fluxOceanfloor.soluteflux_P
soluteflux_O2: fluxOceanfloor.soluteflux_O2

transfer_particulatefluxOceanfloor:
class: ReactionFluxTransfer
parameters:
transfer_matrix: Identity
input_fluxes: fluxOceanfloor.particulateflux_$fluxname$
output_fluxes: particulateflux_$fluxname$
variable_links:

transfer_solutefluxOceanfloor:
class: ReactionFluxTransfer
parameters:
transfer_matrix: Identity
input_fluxes: fluxOceanfloor.soluteflux_$fluxname$
output_fluxes: ocean.oceanfloor.$fluxname$_sms
variable_links:
69 changes: 69 additions & 0 deletions examples/NAPC_tutorial/NAPC_examples.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
using Logging
using OrdinaryDiffEq
using Sundials
using Plots

@info "Start $(@__FILE__)"

@info "importing modules ... (may take a few seconds)"
import PALEOboxes as PB
import PALEOocean
using PALEOmodel
@info " ... done"

global_logger(ConsoleLogger(stderr,Logging.Info))

# include("NAPC_expts.jl")
include("NAPC_reactions.jl")
include("atmreservoirreaction.jl") # for 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__, "NAPC_cfg.yaml"),
"NP_shelf";
modelpars=Dict(
"phys_file"=>joinpath(transport_dir, "S2P3_depth80_m2amp04_phys.nc"),
"surf_file"=>joinpath(transport_dir, "S2P3_depth80_m2amp04_surf.nc"),
)
)

tspan=(0,5.0)

initial_state, modeldata = PALEOmodel.initialize!(model)

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(start=tspan[1], stop=tspan[end], step=0.5))

# individual plots
#plotlyjs(size=(500, 300))
#pager = PALEOmodel.DefaultPlotPager()

# assemble plots onto screens with 6 subplots
gr(size=(1600, 800)) # gr backend (plotly merges captions with multiple panels)
pager=PALEOmodel.PlotPager((2, 3), (legend_background_color=nothing, margin=(5, :mm)))

pager(plot(title="Total P", paleorun.output, "global.total_P"))
pager(plot(title="Total O2", paleorun.output, "global.total_O2"))
pager(:newpage) # flush output

pager(heatmap(title="O2", paleorun.output, "ocean.O2", (column=1, ), clim=(0,Inf)))
pager(heatmap(title="P", paleorun.output, "ocean.P", (column=1, ), clim=(0,Inf)))

pager(heatmap(title="P0", paleorun.output, "ocean.P0", (column=1, ), clim=(0,Inf)))
pager(heatmap(title="P0_growth", paleorun.output, "ocean.P0_growth_rate", (column=1, ), clim=(0,Inf)))
pager(:newpage) # flush output

@info "End $(@__FILE__)"
Loading

0 comments on commit 637deaf

Please sign in to comment.