-
Notifications
You must be signed in to change notification settings - Fork 21
/
Copy pathcolumn.jl
133 lines (93 loc) · 6.74 KB
/
column.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
# This script is not continuously tested so may be out of date
# It is compatible with the version of OceanBioME available at JOSS paper release (v0.X.Y)
using OceanBioME, Oceananigans, Printf, JLD2
using OceanBioME.SLatissimaModel: SLatissima
using Oceananigans.Units
const year = years = 365days
@inline PAR⁰(x, y, t) = 60 * (1 - cos((t + 15days) * 2π / year))*(1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days) ^ 2))) + 2
@inline H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0)
@inline fmld1(t) = H(t, 50days, year) * (1 / (1 +exp(-(t - 100days) / 5days))) * (1 / (1 + exp((t - 330days) / 25days)))
@inline MLD(t) = - (10 + 340 * (1 - fmld1(year - eps(year)) * exp(- mod(t, year) / 25days) - fmld1(mod(t, year))))
@inline κₜ(x, y, z, t) = 1e-2 * (1 + tanh((z - MLD(t)) / 10)) / 2 + 1e-4
@inline t_function(x, y, z, t) = 2.4 * cos(t * 2π / year + 50day) + 10
# ## Grid and PAR field
# Define the grid and an extra Oceananigans field for the PAR to be stored in
Lx, Ly = 20, 20
grid = RectilinearGrid(size=(1, 1, 50), extent=(Lx, Ly, 200))
#####
##### Initial warmup
#####
CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = (args...) -> 35)
N_fixation(x, y, t) = - 10 / year
N_upwelling = Relaxation(rate = 1/10days, mask = (x, y, z) -> ifelse(z < -150, 1.0, 0.0), target = 1.0)
model = NonhydrostaticModel(; grid,
closure = ScalarDiffusivity(ν = κₜ, κ = κₜ),
biogeochemistry = LOBSTER(; grid,
surface_photosynthetically_active_radiation = PAR⁰,
carbonates = true,
variable_redfield = true),
boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),
NO₃ = FieldBoundaryConditions(top = FluxBoundaryCondition(N_fixation))),
forcing = (NO₃ = N_upwelling, ),
advection = nothing)
set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2409.0)
simulation = Simulation(model, Δt = 10minutes, stop_time = 3year)
simulation.callbacks[:nan_checker] = Callback(Oceananigans.Simulations.NaNChecker(; fields = model.tracers, erroring = true))
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n",
iteration(sim),
prettytime(sim),
prettytime(sim.Δt),
prettytime(sim.run_wall_time))
simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))
scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON))
simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite())
wizard = TimeStepWizard(cfl = 0.15, diffusive_cfl = 0.15, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))
simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([3years, 5years - 30days]),
prefix = "warmup_checkpoint", overwrite_existing = true)
run!(simulation) # initial warmup
simulation.stop_time = 5years - 30days
filename = "column"
simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.tracers, model.auxiliary_fields), filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true)
run!(simulation, pickup = false)
# ## Kelp Particle setup
@info "Setting up kelp particles"
n = 5 # number of kelp fronds
z₀ = [-21:5:-1;] * 1.0 # depth of kelp fronds
particles = SLatissima(; x = ones(n) * Lx / 2, y = ones(n) * Ly / 2, z = z₀,
A = ones(n) * 10.0, N = ones(n) * 0.014, C = ones(n) * 0.4,
latitude = 57.5,
scalefactor = 500.0,
pescribed_temperature = t_function)
CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = (args...) -> 35)
N_fixation(x, y, t) = - 10 / year
N_upwelling = Relaxation(rate = 1/10days, mask = (x, y, z) -> ifelse(z < -150, 1.0, 0.0), target = 1.0)
model = NonhydrostaticModel(; grid,
closure = ScalarDiffusivity(ν = κₜ, κ = κₜ),
biogeochemistry = LOBSTER(; grid,
surface_photosynthetically_active_radiation = PAR⁰,
carbonates = true,
variable_redfield = true,
particles),
boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),
NO₃ = FieldBoundaryConditions(top = FluxBoundaryCondition(N_fixation))),
forcing = (NO₃ = N_upwelling, ),
advection = nothing)
model.clock.time = 5year - 30days
set!(model, "warmup_checkpoint_iteration634986.jld2")
simulation = Simulation(model, Δt=1minutes, stop_time=7years)
simulation.callbacks[:nan_checker] = Callback(Oceananigans.Simulations.NaNChecker(; fields = model.tracers, erroring = true))
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n",
iteration(sim),
prettytime(sim),
prettytime(sim.Δt),
prettytime(sim.run_wall_time))
simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))
filename = "kelp"
simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.tracers, model.auxiliary_fields), filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true)
simulation.output_writers[:particles] = JLD2OutputWriter(model, (; particles), filename = "$(filename)_particles.jld2", schedule = TimeInterval(1day), overwrite_existing=true)
simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite())
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))
# ## Run!
# Finally we run the simulation
run!(simulation)