Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added more comments to the subpolar.jl example #21

Merged
merged 1 commit into from
Aug 15, 2022
Merged

Added more comments to the subpolar.jl example #21

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
182 changes: 71 additions & 111 deletions example/subpolar.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,18 @@
"""
This script illustrates how to run OceanBioME as a LOBSTER model
This script illustrates how to run OceanBioME as a 1D column model using the LOBSTER biogeochemical model using a PAR timeseries from the North Atlantic subpolar gyre
This script also shows how to read in data files (in netCDF format) for forcing the model
In this case, the script reads in the mixed layer depth from the Mercator Ocean state esimate and uses this to construct an idealized diffusivity profile
The script also reads in a timeseries of the photosynthetically available radiation for forcing the LOBSTER model
References
(1)2005 A four-dimensional mesoscale map of the spring bloom in the northeast Atlantic (POMME experiment): Results of a prognostic model
(2)2001 Impact of sub-mesoscale physics on production and subduction of phytoplankton in an oligotrophic regime
(3)2012 How does dynamical spatial variability impact 234Th-derived estimates of organic export
(4)2001 Bio-optical properties of oceanic waters: A reappraisal
(1) Lévy, M., Gavart, M., Mémery, L., Caniaux, G. and Paci, A., 2005. A four‐dimensional mesoscale map of the spring bloom in the northeast Atlantic (POMME experiment): Results of a prognostic model. Journal of Geophysical Research: Oceans, 110(C7).
(2) Lévy, M., Klein, P. and Treguier, A.M., 2001. Impact of sub-mesoscale physics on production and subduction of phytoplankton in an oligotrophic regime. Journal of marine research, 59(4), pp.535-565.
(3) Resplandy, L., Martin, A.P., Le Moigne, F., Martin, P., Aquilina, A., Mémery, L., Lévy, M. and Sanders, R., 2012. How does dynamical spatial variability impact 234Th-derived estimates of organic export?. Deep Sea Research Part I: Oceanographic Research Papers, 68, pp.24-45.
(4) Morel, A. and Maritorena, S., 2001. Bio‐optical properties of oceanic waters: A reappraisal. Journal of Geophysical Research: Oceans, 106(C4), pp.7163-7180.
(5) Resplandy, L., Lévy, M., d'Ovidio, F. and Merlivat, L., 2009. Impact of submesoscale variability in estimating the air‐sea CO2 exchange: Results from a model study of the POMME experiment. Global Biogeochemical Cycles, 23(1).
"""

using Random # load required modules
# First, load the needed modules
using Random
using Printf
using Plots
using JLD2
Expand All @@ -16,44 +21,40 @@ using HDF5
using Interpolations
using Statistics

using Oceananigans#9e8cae18-63c1-5223-a75c-80ca9d6e9a09
using Oceananigans
using Oceananigans.Units: second,minute, minutes, hour, hours, day, days, year, years

using BGC

params = LOBSTER.default #load parameters in src/parameters/lobster.jl

######## 1: import annual cycle physics data #Global Ocean 1/12° Physics Analysis and Forecast updated Daily
filename1 = "subpolar_physics.nc" # One may need to modify path if necessary.
#ncinfo(filename1) # ncinfo can be used to check out info in .nc files.
time_series_second = (0:364)days # Start from zero if we don't use extrapolation. We cannot use extrapolation if we want to use annual cycle.
so = ncread(filename1, "so"); #salinity Practical Salinity Unit
# Step 1: import data from the Mercator Ocean state estimate: Global Ocean 1/12° Physics Analysis and Forecast updated Daily
# The temperature and salinity are needed to calculate the air-sea CO2 flux. The mixed layer depth is used to construct an idealized diffusivity profile
filename1 = "subpolar_physics.nc" # Specify the file containing the Mercator model output. One may need to modify path if necessary.
time_series_second = (0:364)days # Create an array of times, one per day in units of seconds. Start from zero if we don't use extrapolation. We cannot use extrapolation if we want to use annual cycle.
so = ncread(filename1, "so"); # Read the salinity in Practical Salinity Units
so_scale_factor = ncgetatt(filename1, "so", "scale_factor") #Real_Value = (Display_Value X scale_factor) + add_offset
so_add_offset = ncgetatt(filename1, "so", "add_offset")
salinity = mean(so, dims=(1,2))[1:365]*so_scale_factor.+so_add_offset # use [1:365] because sometimes one year has 366 days.
salinity_itp = LinearInterpolation(time_series_second, salinity) #converted to interpolations to access them at arbitary time, how to use it: salinity_itp(mod(timeinseconds,364days))
#plot(salinity)
thetao = ncread(filename1, "thetao"); #temperature, Degrees Celsius
salinity = mean(so, dims=(1,2))[1:365]*so_scale_factor.+so_add_offset # NEEDS EXPLANATION
salinity_itp = LinearInterpolation(time_series_second, salinity) # Create a function to interpolate the salinity to a specified time. To use this function, call: salinity_itp(mod(timeinseconds,364days))
thetao = ncread(filename1, "thetao"); # Read the temperature in units of Degrees Celsius
thetao_scale_factor = ncgetatt(filename1, "thetao", "scale_factor")
thetao_add_offset = ncgetatt(filename1, "thetao", "add_offset")
temperature = mean(thetao, dims=(1,2))[1:365]*thetao_scale_factor.+thetao_add_offset
temperature_itp = LinearInterpolation(time_series_second, temperature)
#plot(temperature)
mlotst = ncread(filename1, "mlotst"); #mixed_layer_depth,Meters
temperature = mean(thetao, dims=(1,2))[1:365]*thetao_scale_factor.+thetao_add_offset # NEEDS EXPLANATION
temperature_itp = LinearInterpolation(time_series_second, temperature) # Create a function to interpolate the temperature to a specified time
mlotst = ncread(filename1, "mlotst"); # Read the mixed layer depth in units of eters
mlotst_scale_factor = ncgetatt(filename1, "mlotst", "scale_factor")
mlotst_add_offset = ncgetatt(filename1, "mlotst", "add_offset")
mixed_layer_depth = mean(mlotst, dims=(1,2))[1:365]*mlotst_scale_factor.+mlotst_add_offset
mld_itp = LinearInterpolation(time_series_second, mixed_layer_depth)
#plot(mixed_layer_depth)
mixed_layer_depth = mean(mlotst, dims=(1,2))[1:365]*mlotst_scale_factor.+mlotst_add_offset # NEEDS EXPLANATION
mld_itp = LinearInterpolation(time_series_second, mixed_layer_depth) # Create a function to interpolate the mixed layer depth to a specified time

######## 2: import annual cycle chl data (to calculate PAR_func) #Global Ocean Biogeochemistry Analysis and Forecast
# Step 2: import annual cycle chl data (to calculate PAR_func) #Global Ocean Biogeochemistry Analysis and Forecast
filename2 = "subpolar_chl.nc" #subpolar_chl.nc
chl = ncread(filename2, "chl"); #chl scale_factor=1, add_offset=0
chl_mean = mean(chl, dims=(1,2))[1,1,:,1:365] # mg m-3, unit no need to change.
depth_chl = ncread(filename2, "depth");
#heatmap(1:365, -depth_chl[end:-1:1], chl_mean[end:-1:1,:])

######## 3: import annual cycle surface photosynthetic available radiation (PAR) data #Ocean Color VIIRS-SNPP PAR daily 9km
# Step 3: import annual cycle surface photosynthetic available radiation (PAR) data #Ocean Color VIIRS-SNPP PAR daily 9km
path="./subpolar/" #subtropical #./subpolar/
par_mean_timeseries=zeros(365)
for i in 1:365 #https://discourse.julialang.org/t/leading-zeros/30450
Expand All @@ -69,9 +70,11 @@ surface_PAR_itp = LinearInterpolation((0:364)day, par_mean_timeseries)
surface_PAR(t) = surface_PAR_itp(mod(t, 364days)) # will be used as part of Option 2 of PAR_field

#=
The light (PAR) can be specified in 1 of 2 ways:

Option 1: use PAR as a function, see below PAR_func(x, y, z, t) = PAR_extrap(z, mod(t, 364days))
Obtain PAR_extrap(z, mod(t, 364days)) with annual cycle data chl_mean and surface PAR data par_mean_timeseries as input
Take the following configs:
Use the following steps:
* comment out auxiliary_fields = (PAR=PAR_field, ) in Model instantiation model = NonhydrostaticModel(..., #auxiliary_fields = (PAR=PAR_field, ))
* comment out simulation.callbacks[:update_par]
* indicate PAR_func as an input in bgc = Setup.Oceananigans()
Expand All @@ -80,7 +83,7 @@ Option 2: calculate PAR as a field
Initialize PAR_field = Oceananigans.Fields.Field{Center, Center, Center}(grid)
Update PAR_field through simulation.callbacks[:update_par] and src/Light.jl. Instead of using annual cycle data chl_mean as input, one uses chl as
a function of phytoplankton to calculate PAR. See (A22) in reference (2).
Take the following configs:
Use the following steps:
* keep auxiliary_fields = (PAR=PAR_field, ) in Model instantiation model = NonhydrostaticModel(..., #auxiliary_fields = (PAR=PAR_field, ))
* keep simulation.callbacks[:update_par]
* indicate PAR_field as an input in bgc = Setup.Oceananigans()
Expand All @@ -89,33 +92,27 @@ Take the following configs:
PAR = zeros(length(depth_chl),365) # initialize PAR
PAR_r = zeros(length(depth_chl),365) # break it down into two wavebands: red and blue.
PAR_b = zeros(length(depth_chl),365)
#PAR_g = zeros(length(depth_chl),365)
PAR[1,:] = par_mean_timeseries # surface PAR was obtained from Ocean Color shown above
PAR_r[1,:] = par_mean_timeseries/2 # visible light is split equally between two bands.
PAR_b[1,:] = par_mean_timeseries/2
#PAR_g[1,:] = par_mean_timeseries/3
for i =2:length(depth_chl) # equations refer to the APPENDIX of reference (2) above
PAR_r[i,:]=PAR_r[i-1,:].*exp.(-(params.k_r0.+params.Χ_rp*(chl_mean[i-1,:]).^params.e_r)*(depth_chl[i]-depth_chl[i-1]))
PAR_b[i,:]=PAR_b[i-1,:].*exp.(-(params.k_b0.+params.Χ_bp*(chl_mean[i-1,:]).^params.e_b)*(depth_chl[i]-depth_chl[i-1]))
#PAR_g[i,:]=PAR_g[i-1,:].*exp.(-(params.k_g0.+params.Χ_gp*(chl_mean[i-1,:]).^params.e_g)*(depth_chl[i]-depth_chl[i-1]))
#PAR[i,:]=params.β₁*PAR_b[i,:]+params.β₂*PAR_g[i,:]+params.β₃*PAR_r[i,:]
PAR[i,:]=PAR_b[i,:]+PAR_r[i,:]
end
#heatmap(1:365, -depth_chl[end:-1:1], PAR[end:-1:1,:])
PAR_itp = Interpolations.interpolate((-depth_chl[end:-1:1], (0:364)day), PAR[end:-1:1,:], Gridded(Linear())) #interpolate to access to them at arbitary time and depth
PAR_itp = Interpolations.interpolate((-depth_chl[end:-1:1], (0:364)day), PAR[end:-1:1,:], Gridded(Linear())) # create an interpolation function to allow access to the PAR at arbitary time and depth
PAR_extrap = extrapolate(PAR_itp, (Line(),Throw())) # PAR_extrap(z, mod(t,364days)) Interpolations.extrapolate Method

# Simulation duration
duration=1days #2years
duration=2years

# Define the grid
Lx = 1 #500
Ly = 500
Lx = 1
Ly = 1
Nx = 1
Ny = 1
Nz = 150#150 # number of points in the vertical direction
Lz = 600 # domain depth # subpolar mixed layer depth max 427m

Nz = 150 # number of points in the vertical direction
Lz = 600 # domain depth, for reference, the maximum mixed layer depth is 427m

refinement = 5 # controls spacing near surface (higher means finer spaced) #1.2 5 -4.9118e9
stretching = 2.5 # controls rate of stretching at bottom #24 2.5 5.754
Expand All @@ -132,39 +129,23 @@ grid = RectilinearGrid(size = (Nx, Ny, Nz),
x = (0, Lx),
y = (0, Ly),
z = z_faces)
#=
grid = RectilinearGrid(
size=(Nx, Ny, Nz),
extent=(Lx, Ly, Lz))
=#

B₀ = 0e-8 #m²s⁻³ #buoyancy tracer
N² = 9e-6 #dbdz=N^2, s⁻²
buoyancy_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(B₀),
bottom = GradientBoundaryCondition(N²))

########## u boundary condition
u₁₀ = 0 # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2.5e-3 # dimensionless drag coefficient
ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean
ρₐ = 1.225 # kg m⁻³, average density of air at sea-level
Qᵘ = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

t_function(x, y, z, t) = temperature_itp(mod(t, 364days)) .+ 273.15 # will be used to calculate air-sea-flux
s_function(x, y, z, t) = salinity_itp(mod(t, 364days)) # will be used to calculate air-sea-flux
t_function(x, y, z, t) = temperature_itp(mod(t, 364days)) .+ 273.15 # will be used to calculate air-sea CO2 flux
s_function(x, y, z, t) = salinity_itp(mod(t, 364days)) # will be used to calculate air-sea CO2 flux

PAR_field = Oceananigans.Fields.Field{Center, Center, Center}(grid) #initialize a PAR field
PAR_func(x, y, z, t) = PAR_extrap(z, mod(t, 364days)) #LoadError: cannot define function PAR; it already has a value. So I renamed it as PAR_func. z goes first, then t.
PAR_func(x, y, z, t) = PAR_extrap(z, mod(t, 364days)) # Define the PAR as a function. z goes first, then t.

# Specify the boundary condition for DIC based on the air-sea CO2 flux
dic_bc = Boundaries.setupdicflux(params; forcings=(T=t_function, S=s_function)) # calculate air-sea CO2 flux using DIC, ALK, temperature and salinity.

dic_bc = Boundaries.setupdicflux(params; forcings=(T=t_function, S=s_function)) #calculate air-sea-flux using DIC, ALK, temperature and salinity.
bgc = Setup.Oceananigans(:LOBSTER, grid, params, PAR_field, topboundaries=(DIC=dic_bc, ), optional_sets=(:carbonates, )) #input either PAR_func or PAR_field
# set up the OceanBioME model with the specified biogeochemical model, light, and boundary conditions
bgc = Setup.Oceananigans(:LOBSTER, grid, params, PAR_field, topboundaries=(DIC=dic_bc, ), optional_sets=(:carbonates, )) # input either PAR_func or PAR_field

#npz = Setup.Oceananigans(:NPZ, grid, NPZ.defaults)
#κₜ(x, y, z, t) = 1e-2*max(1-(z+50)^2/50^2,0)+1e-5;
# create a function with the vertical turbulent vertical diffusivity. This is an idealized functional form, but the depth of mixing is based on an interpolation to the mixed layer depth from the Mercator Ocean state estimate
κₜ(x, y, z, t) = 8e-2*max(1-(z+mld_itp(mod(t,364days))/2)^2/(mld_itp(mod(t,364days))/2)^2,0)+1e-4; #setup viscosity and diffusivity in the following Model instantiation

###Model instantiation
# Now, create a 'model' to run in Oceananignas
model = NonhydrostaticModel(advection = UpwindBiasedFifthOrder(),
timestepper = :RungeKutta3,
grid = grid,
Expand All @@ -173,36 +154,28 @@ model = NonhydrostaticModel(advection = UpwindBiasedFifthOrder(),
buoyancy = BuoyancyTracer(),
closure = ScalarDiffusivity(ν=κₜ, κ=κₜ),
forcing = bgc.forcing,
boundary_conditions = merge((u=u_bcs, b=buoyancy_bcs), bgc.boundary_conditions),
boundary_conditions = bgc.boundary_conditions,
auxiliary_fields = (PAR=PAR_field, ) #comment out this line if using functional form of PAR
)

## Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise


#set initial conditions
initial_mixed_layer_depth = -100 # m
stratification(z) = z < initial_mixed_layer_depth ? N² * z : N² * (initial_mixed_layer_depth)
bᵢ(x, y, z) = stratification(z) #

#could initiate P from chl data?
Pᵢ(x,y,z)= (tanh((z+250)/100)+1)/2*(0.038)+0.002 # ((tanh((z+100)/50)-1)/2*0.23+0.23)*16/106
Zᵢ(x,y,z)= (tanh((z+250)/100)+1)/2*(0.038)+0.008 # ((tanh((z+100)/50)-1)/2*0.3+0.3)*16/106
# Initialize the biogeochemical variables
# These initial conditions are set rather arbitrarily in the hope that the model will converge to a repeatable annual cycle if run long enough
Pᵢ(x,y,z)= (tanh((z+250)/100)+1)/2*(0.038)+0.002
Zᵢ(x,y,z)= (tanh((z+250)/100)+1)/2*(0.038)+0.008
Dᵢ(x,y,z)=0
DDᵢ(x,y,z)=0
NO₃ᵢ(x,y,z)= (1-tanh((z+300)/150))/2*6+11.4 # # 17.5*(1-tanh((z+100)/10))/2
NH₄ᵢ(x,y,z)= (1-tanh((z+300)/150))/2*0.05+0.05 #1e-1*(1-tanh((z+100)/10))/2
NO₃ᵢ(x,y,z)= (1-tanh((z+300)/150))/2*6+11.4
NH₄ᵢ(x,y,z)= (1-tanh((z+300)/150))/2*0.05+0.05
DOMᵢ(x,y,z)= 0
DICᵢ(x,y,z)= 2200 # mmol/m^-3
ALKᵢ(x,y,z)= 2400 # mmol/m^-3
DICᵢ(x,y,z)= 2200 # mmol/m^-3
ALKᵢ(x,y,z)= 2400 # mmol/m^-3

## `set!` the `model` fields using functions or constants:
set!(model, b=bᵢ, P=Pᵢ, Z=Zᵢ, D=Dᵢ, DD=DDᵢ, NO₃=NO₃ᵢ, NH₄=NH₄ᵢ, DOM=DOMᵢ,DIC=DICᵢ,ALK=ALKᵢ, u=0, v=0, w=0)
## set the initial conditions using functions or constants:
set!(model, P=Pᵢ, Z=Zᵢ, D=Dᵢ, DD=DDᵢ, NO₃=NO₃ᵢ, NH₄=NH₄ᵢ, DOM=DOMᵢ,DIC=DICᵢ,ALK=ALKᵢ, u=0, v=0, w=0, b=0)

## Setting up a simulation

simulation = Simulation(model, Δt=50, stop_time=duration) #Δt=0.5*(Lz/Nz)^2/1e-2,
## Set up th simulation
simulation = Simulation(model, Δt=50, stop_time=duration) # use a timestep of 50 seconds
# create a model 'callback' to update the light (PAR) profile every 1 timestep
simulation.callbacks[:update_par] = Callback(Light.update_2λ!, IterationInterval(1), merge(params, (surface_PAR=surface_PAR,)))#comment out if using PAR as a function, PAR_func

## Print a progress message
Expand All @@ -212,9 +185,11 @@ progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time:
prettytime(sim.Δt),
prettytime(sim.run_wall_time))

# Create a message to display every 100 timesteps
simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))

#=
SHOULD THIS BE REMOVED?
#this can move into LOBSTER.jl at some point
pco2_bc = zeros(2,round(Int,duration/1day)+3) #Int(duration/simulation.Δt)
function pco2(sim) #https://clima.github.io/OceananigansDocumentation/stable/generated/baroclinic_adjustment/
Expand All @@ -223,43 +198,28 @@ function pco2(sim) #https://clima.github.io/OceananigansDocumentation/stable/ge
pco2_bc[1,round(Int,sim.model.clock.time/1day)+1] = sim.model.clock.time/1day
#sim.model.clock.iteration
end

simulation.callbacks[:pco2] = Callback(pco2, IterationInterval(Int(1day/simulation.Δt))) #callback every 1 day
=#

# We then set up the simulation:

# Vertical slice
# Specify which data to save
# Extract vertical profiles
simulation.output_writers[:profiles] =
JLD2OutputWriter(model, merge(model.velocities, model.tracers, model.auxiliary_fields),
filename = "profile_subpolar_test_PARfield.jld2",
indices = (1, 1, :),
schedule = TimeInterval(1days), #TimeInterval(1days),
overwrite_existing = true)

#simulation.output_writers[:particles] = JLD2OutputWriter(model, (particles=model.particles,),
# prefix = "particles",
# schedule = TimeInterval(1minute),
# force = true)

# We're ready:

# Run the simulation
run!(simulation)

# REMOVE THESE LINES?
#jldsave("pco2_water_subpolar.jld2"; pco2_bc)
#jldopen("pco2_water_subpolar.jld2", "r")


# ## Turbulence visualization
#
# We animate the data saved in `ocean_wind_mixing_and_convection.jld2`.
# We prepare for animating the flow by creating coordinate arrays,
# opening the file, building a vector of the iterations that we saved
# data at, and defining functions for computing colorbar limits:

## Coordinate arrays
xw, yw, zw = nodes(model.velocities.w)
xb, yb, zb = nodes(model.tracers.b)

# Load and plot the results
results = BGC.Plot.load_tracers(simulation)
BGC.Plot.profiles(results)

# Save the plot to a PDF file
savefig("subpolar_test_PARfield.pdf")