A Phosphorus Cycling Model

A Phosphorus Cycling Model

Tip

This example is also available as a Jupyter notebook: p_cycle_DIP_DOP_POP.ipynb

Tracer equations

We consider a simple model for the cycling of phosphorus with 3 state variables consisting of phosphate (PO₄) AKA dissolved inorganic phosphorus (DIP), dissolved organic phosphorus (DOP), and particulate organic phosphorus (POP).

The dissolved phases are transported by advection and diffusion whereas the particulate phase sinks rapidly down the water column without any appreciable transport by the circulation.

The governing equations are:

\[\frac{\partial}{\partial t} DIP + \nabla \cdot \left[\boldsymbol{u} + \mathbf{K}\cdot\nabla \right] DIP = -\gamma(DIP) + \kappa_\mathsf{D} \, DOP,\]
\[\frac{\partial}{\partial t} DOP + \nabla \cdot \left[\boldsymbol{u} + \mathbf{K}\cdot \nabla \right] DOP = \sigma \, \gamma(DIP) + \kappa_\mathsf{P} \, POP - \kappa_\mathsf{D} \, DOP,\]

and

\[\frac{\partial}{\partial t} POP + \frac{\partial}{\partial z} \left[w_\mathsf{P} \, POP\right] = (1-\sigma) \, \gamma(DIP) - \kappa_\mathsf{P} \, POP,\]

where $\boldsymbol{u}$ is the fluid velocity and $\mathbf{K}$ is the eddy-diffusion tensor. Thus, $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right]$ is a differential operator that represents the transport by the ocean circulation. The function $\gamma(DIP)$ represents the biological uptake of DIP by phytoplankton.

Oxygen participates to this cycle too and satisfies its own tracer equation

\[\frac{\partial}{\partial t} DO_2 + \nabla \cdot \left[\boldsymbol{u} + \mathbf{K}\cdot\nabla \right] DO_2 = -r_{\mathsf{O}_2:\mathsf{P}} \, \kappa_\mathsf{D} \, DOP + \boldsymbol{\Lambda}(DO_2 - [O_2]_{\mathsf{sat}})\]

where $\Lambda$ is the air-sea gas exchange operator.

These tracer equations depend on a number of scalars, that we list below

SymbolDefinition
$w_\mathsf{P}$depth dependent particle sinking speed
$\sigma$fraction of the organic matter production allocated to the dissolved phase
$\kappa_\mathsf{D}$respiration rate for dissolved organic matter (DOP → DIP)
$\kappa_\mathsf{P}$dissolution rate for particulate organic matter (POP → DOP)
$r_{\mathsf{O}_2:\mathsf{P}}$number of moles of O₂ needed to respire 1 mole of DOP
using AIBECS

Load the circulation and grid

const wet3d, grd, T_Circulation = OCIM1.load()
([0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 1.0 1.0 … 1.0 1.0; 0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 1.0 1.0 … 1.0 1.0; 0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 1.0 1.0 … 1.0 1.0; 0.0 0.0 … 0.0 0.0]

...

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], Dict{String,Any}("XT" => [1.0 3.0 … 357.0 359.0; 1.0 3.0 … 357.0 359.0; … ; 1.0 3.0 … 357.0 359.0; 1.0 3.0 … 357.0 359.0],"Areat" => [8.442828573400822e8 8.442828573400822e8 … 8.442828573400822e8 8.442828573400822e8; 2.53184242418399e9 2.53184242418399e9 … 2.53184242418399e9 2.53184242418399e9; … ; 2.531842424184019e9 2.531842424184019e9 … 2.531842424184019e9 2.531842424184019e9; 8.442828573400851e8 8.442828573400851e8 … 8.442828573400851e8 8.442828573400851e8],"YC" => [-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0],"dxu" => [2.0 2.0 … 2.0 2.0],"YT3d" => [-89.01098901098901 -89.01098901098901 … -89.01098901098901 -89.01098901098901; -87.03296703296704 -87.03296703296704 … -87.03296703296704 -87.03296703296704; … ; 87.03296703296702 87.03296703296702 … 87.03296703296702 87.03296703296702; 89.01098901098901 89.01098901098901 … 89.01098901098901 89.01098901098901]

[-89.01098901098901 -89.01098901098901 … -89.01098901098901 -89.01098901098901; -87.03296703296704 -87.03296703296704 … -87.03296703296704 -87.03296703296704; … ; 87.03296703296702 87.03296703296702 … 87.03296703296702 87.03296703296702; 89.01098901098901 89.01098901098901 … 89.01098901098901 89.01098901098901]

[-89.01098901098901 -89.01098901098901 … -89.01098901098901 -89.01098901098901; -87.03296703296704 -87.03296703296704 … -87.03296703296704 -87.03296703296704; … ; 87.03296703296702 87.03296703296702 … 87.03296703296702 87.03296703296702; 89.01098901098901 89.01098901098901 … 89.01098901098901 89.01098901098901]

...

[-89.01098901098901 -89.01098901098901 … -89.01098901098901 -89.01098901098901; -87.03296703296704 -87.03296703296704 … -87.03296703296704 -87.03296703296704; … ; 87.03296703296702 87.03296703296702 … 87.03296703296702 87.03296703296702; 89.01098901098901 89.01098901098901 … 89.01098901098901 89.01098901098901]

[-89.01098901098901 -89.01098901098901 … -89.01098901098901 -89.01098901098901; -87.03296703296704 -87.03296703296704 … -87.03296703296704 -87.03296703296704; … ; 87.03296703296702 87.03296703296702 … 87.03296703296702 87.03296703296702; 89.01098901098901 89.01098901098901 … 89.01098901098901 89.01098901098901]

[-89.01098901098901 -89.01098901098901 … -89.01098901098901 -89.01098901098901; -87.03296703296704 -87.03296703296704 … -87.03296703296704 -87.03296703296704; … ; 87.03296703296702 87.03296703296702 … 87.03296703296702 87.03296703296702; 89.01098901098901 89.01098901098901 … 89.01098901098901 89.01098901098901],"YV3d" => [-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0]

[-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0]

[-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0]

...

[-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0]

[-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0]

[-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0],"DYT3d" => [219946.0087474788 219946.0087474788 … 219946.0087474788 219946.0087474788; 219946.00874747802 219946.00874747802 … 219946.00874747802 219946.00874747802; … ; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796]

[219946.0087474788 219946.0087474788 … 219946.0087474788 219946.0087474788; 219946.00874747802 219946.00874747802 … 219946.00874747802 219946.00874747802; … ; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796]

[219946.0087474788 219946.0087474788 … 219946.0087474788 219946.0087474788; 219946.00874747802 219946.00874747802 … 219946.00874747802 219946.00874747802; … ; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796]

...

[219946.0087474788 219946.0087474788 … 219946.0087474788 219946.0087474788; 219946.00874747802 219946.00874747802 … 219946.00874747802 219946.00874747802; … ; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796]

[219946.0087474788 219946.0087474788 … 219946.0087474788 219946.0087474788; 219946.00874747802 219946.00874747802 … 219946.00874747802 219946.00874747802; … ; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796]

[219946.0087474788 219946.0087474788 … 219946.0087474788 219946.0087474788; 219946.00874747802 219946.00874747802 … 219946.00874747802 219946.00874747802; … ; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796],"dxv" => [2.0 2.0 … 2.0 2.0],"YV" => [-88.02197802197801 -88.02197802197801 … -88.02197802197801 -88.02197802197801; -86.04395604395604 -86.04395604395604 … -86.04395604395604 -86.04395604395604; … ; 88.02197802197801 88.02197802197801 … 88.02197802197801 88.02197802197801; 90.0 90.0 … 90.0 90.0],"DYV" => [219946.00874747802 219946.00874747802 … 219946.00874747802 219946.00874747802; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; … ; 219946.0087474796 219946.0087474796 … 219946.0087474796 219946.0087474796; 219946.0087474788 219946.0087474788 … 219946.0087474788 219946.0087474788]…), 
  [1     ,      1]  =  0.000197784
  [2     ,      1]  =  1.05541e-8
  [10384 ,      1]  =  -2.09257e-7
  [10442 ,      1]  =  -0.000191611
  [10443 ,      1]  =  4.80961e-9
  [20825 ,      1]  =  -1.83059e-9
  [20883 ,      1]  =  7.967e-9
  [1     ,      2]  =  -5.98052e-8
  [2     ,      2]  =  0.000187532
  ⋮
  [200160, 200159]  =  -2.04829e-8
  [197886, 200160]  =  2.41231e-9
  [199766, 200160]  =  6.70985e-9
  [199777, 200160]  =  -1.26357e-9
  [199778, 200160]  =  -7.22216e-9
  [199779, 200160]  =  7.59325e-9
  [199790, 200160]  =  -7.41023e-9
  [200156, 200160]  =  -3.07748e-8
  [200159, 200160]  =  -2.14741e-8
  [200160, 200160]  =  5.22074e-8)

Define useful constants and arrays

const iwet = indices_of_wet_boxes(wet3d)
const nb = number_of_wet_boxes(wet3d)
const v = vector_of_volumes(wet3d, grd)
const z = vector_of_depths(wet3d, grd)
const ztop = vector_of_top_depths(wet3d, grd)
200160-element Array{Float64,1}:
    0.0           
    0.0           
    0.0           
    0.0           
    0.0           
    0.0           
    0.0           
    0.0           
    0.0           
    0.0           
    ⋮             
 5116.506284367635
 5116.506284367635
 5116.506284367635
 5116.506284367635
 5116.506284367635
 5116.506284367635
 5116.506284367635
 5116.506284367635
 5116.506284367635

And matrices

const DIV = buildDIV(wet3d, iwet, grd)
const Iabove = buildIabove(wet3d, iwet) ;
200160×200160 SparseMatrixCSC{Bool,Int64} with 189719 stored entries:
  [10442 ,      1]  =  1
  [10443 ,      2]  =  1
  [10444 ,      3]  =  1
  [10445 ,      4]  =  1
  [10446 ,      5]  =  1
  [10447 ,      6]  =  1
  [10448 ,      7]  =  1
  [10449 ,      8]  =  1
  [10450 ,      9]  =  1
  ⋮
  [200151, 199755]  =  1
  [200152, 199756]  =  1
  [200153, 199757]  =  1
  [200154, 199764]  =  1
  [200155, 199765]  =  1
  [200156, 199766]  =  1
  [200157, 199767]  =  1
  [200158, 199776]  =  1
  [200159, 199777]  =  1
  [200160, 199778]  =  1

Transport matrices

T_DIP(p) = T_Circulation
T_DOP(p) = T_Circulation
T_DO2(p) = T_Circulation
const S₀ = buildPFD(ones(nb), DIV, Iabove)
const S′ = buildPFD(ztop, DIV, Iabove)
function T_POP(p)
    w₀, w′ = p.w₀, p.w′
    return w₀ * S₀ + w′ * S′
end
T_all = (T_DIP, T_DOP, T_POP, T_DO2)
(Main.ex-p_cycle_DIP_DOP_POP.T_DIP, Main.ex-p_cycle_DIP_DOP_POP.T_DOP, Main.ex-p_cycle_DIP_DOP_POP.T_POP, Main.ex-p_cycle_DIP_DOP_POP.T_DO2)

Because AIBECS will solve for the steady state solution directly without time-stepping the goverining equations to equilibrium, we don't have any opportunity to specify any intial conditions. Initial conditions are how the total amount of conserved elements get specified in most global biogeochemical modelels. Thus to specify the total inventory of P in AIBECS we add a very weak resporing term to the DIP equation. The time-scale for this restoring term is chosen to be very long compared to the timescale with which the ocean circulation homogenizes a tracer. Because of this long timescale we call it the geological restoring term, but geochemists who work on geological processes don't like that name! In any event the long timescale allows us to prescribe the total inventory of P without having any appreciable impact on the 3d distribution of P.

Sources minus sinks

Geological Restoring
function geores(x, p)
    τg, xgeo = p.τg, p.xgeo
    return (xgeo .- x) / τg
end
geores (generic function with 1 method)
Uptake of phosphate (DIP)
relu(x) = (x .≥ 0) .* x
function uptake(DIP, p)
    τu, ku, z₀ = p.τu, p.ku, p.z₀
    DIP⁺ = relu(DIP)
    return 1/τu * DIP⁺.^2 ./ (DIP⁺ .+ ku) .* (z .≤ z₀)
end
uptake (generic function with 1 method)
Remineralization DOP into DIP
function remineralization(DOP, p)
    κDOP = p.κDOP
    return κDOP * DOP
end
remineralization (generic function with 1 method)
Dissolution of POP into DOP
function dissolution(POP, p)
    κPOP = p.κPOP
    return κPOP * POP
end
dissolution (generic function with 1 method)
Air-sea gas exchange
dz1 = grd["dzt"][1]               # thickness of the top layer
z = vec(grd["ZT3d"])[iwet]        # depth of the gridbox centers
using WorldOceanAtlasTools
WOA = WorldOceanAtlasTools
μDO2 , σ²DO2 = WOA.fit_to_grid(grd,2018,"O2sat","annual","1°","an")
([0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.9986298060417176 0.9981694602966309 … 1.000671501159668 0.9994396591186524; 1.0175252294540404 1.0174463748931886 … 1.0179427027702332 1.0176871752738952]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.9550732585362026 0.9562181799752372 … 0.9535434395926339 0.9541520363943917; 0.9528527014596122 0.9531717082432338 … 0.9524201202392578 0.9525935690743583]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.8482290013631184 0.8513909212748211 … 0.8417107359568279 0.8448687076568604; 0.836611909866333 0.8373468748728435 … 0.8351414712270101 0.8358594926198323]

...

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.8468556912740072 0.8483416112263997 … 0.8462580108642578 0.8461976095346304; 0.8393070729573567 0.8395611953735351 … 0.8417618124825614 0.8413471494402204]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [Inf Inf … Inf Inf; Inf Inf … Inf Inf; … ; 4.483668471052279e-5 4.483668471052279e-5 … 4.483668471052279e-5 4.483668471052279e-5; 6.510639911157342e-5 6.60700174907106e-5 … 5.990546374687256e-5 6.299638850869087e-5]

[Inf Inf … Inf Inf; Inf Inf … Inf Inf; … ; 0.0014753435689954131 0.0014058725722332383 … 0.0016401173756970821 0.0015562222898595792; 0.0020962392960074796 0.002079495734697957 … 0.002135013684531441 0.0021153090527070164]

[Inf Inf … Inf Inf; Inf Inf … Inf Inf; … ; 6.02427243913553e-5 6.371971117605426e-5 … 5.616889233288021e-5 5.7314006471460744e-5; 6.847410248901724e-5 7.391922987526414e-5 … 6.143014941535512e-5 6.413035838413634e-5]

...

[Inf Inf … Inf Inf; Inf Inf … Inf Inf; … ; 4.483668471052279e-5 4.483668471052279e-5 … 4.483668471052279e-5 4.483668471052279e-5; 4.483668471052279e-5 4.483668471052279e-5 … 4.483668471052279e-5 4.483668471052279e-5]

[Inf Inf … Inf Inf; Inf Inf … Inf Inf; … ; Inf Inf … Inf Inf; Inf Inf … Inf Inf]

[Inf Inf … Inf Inf; Inf Inf … Inf Inf; … ; Inf Inf … Inf Inf; Inf Inf … Inf Inf])

This below was commented out?

function airsea(DO2, p)
    κDO2 = p.κDO2
    return κDO2 * (z .< 20) .* (1.0 .- DO2) / dz1
end
airsea (generic function with 1 method)

Add them up into sms functions (Sources Minus Sinks)

function sms_DIP(DIP, DOP, POP, p)
    return -uptake(DIP, p) + remineralization(DOP, p) + geores(DIP, p)
end
function sms_DOP(DIP, DOP, POP, p)
    σ = p.σ
    return σ * uptake(DIP, p) - remineralization(DOP, p) + dissolution(POP, p)
end
function sms_POP(DIP, DOP, POP, p)
    σ = p.σ
    return (1 - σ) * uptake(DIP, p) - dissolution(POP, p)
end
sms_all = (sms_DIP, sms_DOP, sms_POP) # bundles all the source-sink functions in a tuple
(Main.ex-p_cycle_DIP_DOP_POP.sms_DIP, Main.ex-p_cycle_DIP_DOP_POP.sms_DOP, Main.ex-p_cycle_DIP_DOP_POP.sms_POP)

Parameters

Build the parameters type and p₀

t = empty_parameter_table()    # initialize table of parameters
add_parameter!(t, :xgeo, 2.17u"mmol/m^3",
    variance_obs = ustrip(upreferred(0.1 * 2.17u"mmol/m^3"))^2,
    description = "Geological mean P concentration",
    LaTeX = "\\state^\\mathrm{geo}")
add_parameter!(t, :τg, 1.0u"Myr",
    description = "Geological restoring timescale",
    LaTeX = "\\tau_\\mathrm{geo}")
add_parameter!(t, :ku, 10.0u"μmol/m^3",
    optimizable = true,
    description = "Half-saturation constant (Michaelis-Menten)",
    LaTeX = "k_\\vec{u}")
add_parameter!(t, :z₀, 80.0u"m",
    description = "Depth of the euphotic layer base",
    LaTeX = "z_0")
add_parameter!(t, :w₀, 1.0u"m/d",
    optimizable = true,
    description = "Sinking velocity at surface",
    LaTeX = "w_0")
add_parameter!(t, :w′, 1/4.4625u"d",
    optimizable = true,
    description = "Vertical gradient of sinking velocity",
    LaTeX = "w'")
add_parameter!(t, :κDOP, 1/0.25u"yr",
    optimizable = true,
    description = "Remineralization rate constant (DOP to DIP)",
    LaTeX = "\\kappa")
add_parameter!(t, :κPOP, 1/5.25u"d",
    optimizable = true,
    description = "Dissolution rate constant (POP to DOP)",
    LaTeX = "\\kappa")
add_parameter!(t, :σ, 0.3u"1",
    description = "Fraction of quick local uptake recycling",
    LaTeX = "\\sigma")
add_parameter!(t, :τu, 30.0u"d",
    optimizable = true,
    description = "Maximum uptake rate timescale",
    LaTeX = "\\tau_\\vec{u}")
initialize_Parameters_type(t, "Pcycle_Parameters")   # Generate the parameter type

Generate state function and Jacobian

nt = length(T_all)    # number of tracers
n = nt * nb           # total dimension of the state vector
p = Pcycle_Parameters() # parameters
x = p.xgeo * ones(n) # initial iterate
800640-element Array{Float64,1}:
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 ⋮      
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217
 0.00217

F, ∇ₓF = statefunctionandJacobian(Tall, sms_all, nb)

and solve

prob = SteadyStateProblem(F, ∇ₓF, x, p) s = solve(prob, CTKAlg());

unpack state

function unpackPstate(s,mask) DIP = NaNmask; DOP = NaNmask; POP = NaN*mask iwet = findall(x-> x==1, vec(mask)) nwet = length(iwet) idip = 1:nwet; idop = idip.+nwet; ipop = idop.+nwet DIP[iwet] = s[idip]; DOP[iwet] = s[idop]; POP[iwet] = s[ipop] return DIP, DOP, POP end DIP, DOP, POP = unpackPstate(s,wet3d);

We will plot the concentration of DIP at a given depth horizon

depth = vec(grd["zt"])
iz = findfirst(depth .> 200)
iz, depth[iz]
(6, 246.7350746268657)

dip = DIP[:,:,iz] * ustrip(1.0u"mol/m^3"|>u"mmol/m^3"); lat, lon = vec(grd["yt"]), vec(grd["xt"]);

and plot

ENV["MPLBACKEND"]="qt5agg" using PyPlot, PyCall using Conda; Conda.add("Cartopy") clf() ccrs = pyimport("cartopy.crs") ax = subplot(projection=ccrs.Robinson(central_longitude=-155.0)) ax.coastlines()

making it cyclic for Cartopy

loncyc = [lon; 360+lon[1]] dipcyc = hcat(dip, dip[:,1])

And plot

p = contourf(loncyc, lat, dipcyc, levels=0:0.2:3.6, transform=ccrs.PlateCarree(), zorder=-1) colorbar(p, orientation="horizontal"); gcf()

This page was generated using Literate.jl.