Radiocarbon with OCIM

using AIBECS
t = empty_parameter_table()               # initialize table of parameters
add_parameter!(t, :λ, 1 / (5730*log(2))u"yr") # add the radioactive decay e-folding timescale
add_parameter!(t, :κ, 50u"m" / 10u"yr")
initialize_Parameters_type(t, "C14_OCIM_Parameters")
t

We generate the parameters via

p = C14_OCIM_Parameters()

     λ = 2.52e-04 [yr⁻¹] (fixed)
     κ = 5.00e+00 [m yr⁻¹] (fixed)

We load the OCIM matrix via

const mask, grd, T_OCIM = OCIM1.load()
T¹⁴C(p) = T_OCIM
iwet = indices_of_wet_boxes(mask);  # index to wet gridboxes

200160-element Array{Int64,1}:
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
      ⋮
 390450
 390451
 390539
 390540
 390541
 390542
 390630
 390631
 390632

and some useful constants

dz1 = grd["dzt"][1]               # thickness of the top layer
z = vector_of_depths(mask, grd)
nb = number_of_wet_boxes(mask)

200160

source minus sinks

function SMS¹⁴C(R, p) # source minus sink of age
    λ = p.λ
    κ = p.κ
    return κ * (z .< 20) .* (1.0 .- R) / dz1 - λ.*R
end

SMS¹⁴C (generic function with 1 method)

Generate state function and Jacobian

x = ones(nb)
Ts = (T¹⁴C,)           # bundles all the transport matrices in a tuple
SMSs = (SMS¹⁴C,) # bundles all the source-sink functions in a tuple
F, ∇ₓF = state_function_and_Jacobian(Ts, SMSs, nb) # generates the state function (and its Jacobian!)
F(x,p) # to check it works

200160-element Array{Float64,1}:
 -7.978366759081164e-12
 -7.978683586612089e-12
 -7.97824989349255e-12 
 -7.978615770633085e-12
 -7.978460079492284e-12
 -7.978297176901567e-12
 -7.978408730498e-12   
 -7.978270699884135e-12
 -7.978355798768889e-12
 -7.97849788386643e-12 
  ⋮                    
 -7.978402643025237e-12
 -7.978402643038472e-12
 -7.97840264304509e-12 
 -7.978402643015311e-12
 -7.978402643041781e-12
 -7.978402643041781e-12
 -7.978402643031855e-12
 -7.978402643041781e-12
 -7.978402643041781e-12

Set up the problem in AIBECS

prob = SteadyStateProblem(F, ∇ₓF, x, p)
R = solve(prob, CTKAlg())
C14_age = -log.(R) / p.λ;

200160-element Array{Float64,1}:
 3.942273096560254e10 
 3.980081330050604e10 
 4.039395435254564e10 
 3.998797942292548e10 
 3.742935022296953e10 
 3.669643544760427e10 
 3.753248046209904e10 
 3.703244437629167e10 
 3.577737704586382e10 
 3.402738749624555e10 
 ⋮                    
 3.2269848925595623e10
 3.2059184113358418e10
 3.2539013376892643e10
 3.240752129878833e10 
 3.2187009397144234e10
 3.1753835571645058e10
 3.237799887583362e10 
 3.225366161852383e10 
 3.21011466642212e10  

Now we plot

lat, lon = vec(grd["yt"]), vec(grd["xt"]);
C14_age_3d = NaN * mask     # creates a 3D array of NaNs
C14_age_3d[iwet] = C14_age   # Fills the wet grid boxes with the age values
size(C14_age)            # Just to check the size

(200160,)

Chose the depth index of the slice we want to plot

depth = vec(grd["zt"])
iz = findfirst(depth .> 700)
iz, depth[iz]

(11, 757.7081696779262)

get slice and convert to years

C14_age_3d_1000m_yr = C14_age_3d[:,:,iz] * ustrip(1.0u"s" |> u"yr")

91×180 Array{Float64,2}:
  NaN      NaN      NaN      NaN      …  NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN         NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN         NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN         NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN         NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN      …  NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN         NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN         NaN      NaN      NaN       NaN    
  NaN      NaN      NaN      NaN         NaN      NaN      NaN       NaN    
 1342.82   NaN      NaN      NaN         NaN      NaN      NaN      1342.79 
    ⋮                                 ⋱                                     
  551.565  551.434  552.414  553.257     552.514  551.17   551.323   552.009
  548.607  550.837  552.478  553.964     545.618  544.095  545.933   547.237
  555.314  556.118  557.729  556.933     552.949  551.708  553.125   554.407
  561.17   561.348  559.361  555.091  …  563.126  560.959  558.9     560.849
  566.044  568.603  566.295  563.76      565.157  564.015  562.776   564.358
  557.612  560.431  562.178  563.778     552.645  552.367  553.436   555.091
  552.427  554.31   556.009  557.609     548.683  549.213  549.925   550.787
  549.492  550.181  550.84   551.239     547.947  548.061  548.057   548.836
  NaN      NaN      NaN      NaN      …  NaN      NaN      NaN       NaN    

and plot

ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
clf()
ccrs = pyimport("cartopy.crs")
ax = subplot(projection=ccrs.Robinson(central_longitude=-155.0))
ax.coastlines()
lon_cyc = [lon; 360+lon[1]] # making it cyclic for Cartopy
age_cyc = hcat(C14_age_3d_1000m_yr, C14_age_3d_1000m_yr[:,1])
p = contourf(lon_cyc, lat, age_cyc, levels=0:100:3600, transform=ccrs.PlateCarree(), zorder=-1)
colorbar(p, orientation="horizontal")
gcf()

This page was generated using Literate.jl.