Skip to content

Commit 916bbca

Browse files
briochemcbriochemc
committed
Add MOC plot for Reviewer 2
1 parent 6b9f93d commit 916bbca

File tree

2 files changed

+296
-10
lines changed

2 files changed

+296
-10
lines changed

‎src/plot_MOCs_for_Reviewer2.jl‎

Lines changed: 285 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,285 @@
1+
# qsub -I -P xv83 -q express -l mem=47GB -l storage=scratch/gh0+scratch/xv83 -l walltime=01:00:00 -l ncpus=12
2+
# This is Fig. S7? in Pasquier et al. (GRL, 2025) # TODO: check figure number
3+
4+
using Pkg
5+
Pkg.activate(".")
6+
Pkg.instantiate()
7+
8+
using OceanTransportMatrixBuilder
9+
using NetCDF
10+
using YAXArrays
11+
using DataFrames
12+
using DimensionalData
13+
# using SparseArrays
14+
# using LinearAlgebra
15+
using Unitful
16+
using Unitful: s, yr, m, kg
17+
using CairoMakie
18+
using GeoMakie
19+
using Interpolations
20+
using OceanBasins
21+
using Statistics
22+
using StatsBase
23+
using NaNStatistics
24+
using Format
25+
using GeometryBasics
26+
using LibGEOS
27+
using GeometryOps
28+
29+
include("plotting_functions.jl")
30+
31+
model = "ACCESS-ESM1-5"
32+
33+
time_windows = [
34+
"Jan1850-Dec1859",
35+
"Jan2030-Dec2039",
36+
"Jan2090-Dec2099",
37+
]
38+
39+
experiments = [
40+
"historical",
41+
"ssp370",
42+
"ssp370",
43+
]
44+
experiments2 = [
45+
"historical",
46+
"SSP-3.70",
47+
"SSP-3.70",
48+
]
49+
50+
members = ["r$(r)i1p1f1" for r in 1:40]
51+
52+
lumpby = "month"
53+
months = 1:12
54+
55+
56+
# Gadi directory for input files
57+
fixedvarsinputdir = "/scratch/xv83/TMIP/data/$model"
58+
# Load areacello and volcello for grid geometry
59+
volcello_ds = open_dataset(joinpath(fixedvarsinputdir, "volcello.nc"))
60+
areacello_ds = open_dataset(joinpath(fixedvarsinputdir, "areacello.nc"))
61+
62+
63+
# Load fixed variables in memory
64+
areacello = readcubedata(areacello_ds.areacello)
65+
volcello = readcubedata(volcello_ds.volcello)
66+
lon = readcubedata(volcello_ds.lon)
67+
lat = readcubedata(volcello_ds.lat)
68+
lev = volcello_ds.lev
69+
# Identify the vertices keys (vary across CMIPs / models)
70+
volcello_keys = propertynames(volcello_ds)
71+
lon_vertices_key = volcello_keys[findfirst(x -> occursin("lon", x) & occursin("vert", x), string.(volcello_keys))]
72+
lat_vertices_key = volcello_keys[findfirst(x -> occursin("lat", x) & occursin("vert", x), string.(volcello_keys))]
73+
lon_vertices = readcubedata(getproperty(volcello_ds, lon_vertices_key))
74+
lat_vertices = readcubedata(getproperty(volcello_ds, lat_vertices_key))
75+
76+
# Make makegridmetrics
77+
gridmetrics = makegridmetrics(; areacello, volcello, lon, lat, lev, lon_vertices, lat_vertices)
78+
(; lon_vertices, lat_vertices) = gridmetrics
79+
80+
# Make indices
81+
indices = makeindices(gridmetrics.v3D)
82+
83+
84+
ϕsnorth = map(time_windows, experiments) do time_window, experiment
85+
86+
ϕsnorth_ensemblemean = mean(map(members) do member
87+
88+
inputdir = "/scratch/xv83/TMIP/data/$model/$experiment/$member/$(time_window)"
89+
cycloinputdir = joinpath(inputdir, "cyclomonth")
90+
umo_ds = open_dataset(joinpath(cycloinputdir, "umo.nc"))
91+
vmo_ds = open_dataset(joinpath(cycloinputdir, "vmo.nc"))
92+
ψᵢGM_ds = open_dataset(joinpath(cycloinputdir, "tx_trans_gm.nc"))
93+
ψⱼGM_ds = open_dataset(joinpath(cycloinputdir, "ty_trans_gm.nc"))
94+
ψᵢsubmeso_ds = open_dataset(joinpath(cycloinputdir, "tx_trans_submeso.nc"))
95+
ψⱼsubmeso_ds = open_dataset(joinpath(cycloinputdir, "ty_trans_submeso.nc"))
96+
mlotst_ds = open_dataset(joinpath(cycloinputdir, "mlotst.nc"))
97+
98+
mean_days_in_month = umo_ds.mean_days_in_month |> Array
99+
w = Weights(mean_days_in_month)
100+
101+
umo = dropdims(mean(readcubedata(umo_ds.umo), w; dims=:month), dims=:month)
102+
vmo = dropdims(mean(readcubedata(vmo_ds.vmo), w; dims=:month), dims=:month)
103+
104+
ψᵢGM = dropdims(mean(readcubedata(ψᵢGM_ds.tx_trans_gm), w; dims=:month), dims=:month)
105+
ψⱼGM = dropdims(mean(readcubedata(ψⱼGM_ds.ty_trans_gm), w; dims=:month), dims=:month)
106+
ψᵢsubmeso = dropdims(mean(readcubedata(ψᵢsubmeso_ds.tx_trans_submeso), w; dims=:month), dims=:month)
107+
ψⱼsubmeso = dropdims(mean(readcubedata(ψⱼsubmeso_ds.ty_trans_submeso), w; dims=:month), dims=:month)
108+
109+
# Replace missing values and convert to arrays
110+
# I think latest YAXArrays converts _FillValues to missing
111+
ψᵢGM = replace(ψᵢGM |> Array, missing => 0) .|> Float64
112+
ψⱼGM = replace(ψⱼGM |> Array, missing => 0) .|> Float64
113+
ψᵢsubmeso = replace(ψᵢsubmeso |> Array, missing => 0) .|> Float64
114+
ψⱼsubmeso = replace(ψⱼsubmeso |> Array, missing => 0) .|> Float64
115+
116+
# Also remove missings in umo and vmo
117+
umo = replace(umo, missing => 0)
118+
vmo = replace(vmo, missing => 0)
119+
120+
# Take the vertical diff of zonal/meridional transport diagnostics to get their mass transport
121+
(nx, ny, _) = size(ψᵢGM)
122+
ϕᵢGM = diff([fill(0.0, nx, ny, 1);;; ψᵢGM |> Array], dims=3)
123+
ϕⱼGM = diff([fill(0.0, nx, ny, 1);;; ψⱼGM |> Array], dims=3)
124+
ϕᵢsubmeso = diff([fill(0.0, nx, ny, 1);;; ψᵢsubmeso |> Array], dims=3)
125+
ϕⱼsubmeso = diff([fill(0.0, nx, ny, 1);;; ψⱼsubmeso |> Array], dims=3)
126+
127+
# TODO fix incompatible dimensions betwewen umo and ϕᵢGM/ϕᵢsubmeso Dim{:i} and Dim{:xu_ocean}
128+
ϕ = let umo = umo + ϕᵢGM + ϕᵢsubmeso, vmo = vmo + ϕⱼGM + ϕⱼsubmeso
129+
facefluxesfrommasstransport(; umo, vmo, gridmetrics, indices)
130+
end
131+
132+
return ϕ.north
133+
end)
134+
135+
end
136+
137+
# unpack model grid
138+
(; lon, lat, zt, v3D,) = gridmetrics
139+
lev = zt
140+
# unpack indices
141+
(; wet3D, N) = indices
142+
143+
# basins
144+
basin_keys = (:ATL, :INDOPAC, :GBL)
145+
basin_strs = ("Atlantic", "Indo-Pacific", "Global")
146+
OCEANS = OceanBasins.oceanpolygons()
147+
isatlnoSO(lat, lon, o) = isatlantic(lat, lon, o) .& (lat .≥ -30)
148+
isindopacific(lat, lon, o) = (isindian(lat, lon, o) .| ispacific(lat, lon, o)) .& (lat .≥ -30)
149+
isglobal(lat, lon, o) = trues(size(lat))
150+
basin_functions = (isatlnoSO, isindopacific, isglobal)
151+
basin_values = (reshape(f(lat[:], lon[:], OCEANS), size(lat)) for f in basin_functions)
152+
basins = (; (basin_keys .=> basin_values)...)
153+
basin_latlims_values = [clamp.(0 .* (-5, +5) .+ extrema(lat[.!isnan.(v3D[:,:,1]) .& basin[:,:,1]]), -80, 80) for basin in basins]
154+
basin_latlims = (; (basin_keys .=> basin_latlims_values)...)
155+
156+
157+
158+
# Plot meridional overturning circulation for each basin
159+
160+
ρ = 1035.0 # density (kg/m^3)
161+
162+
levels = -24:2:24
163+
colormap = cgrad(:curl, length(levels) + 1; categorical=true, rev=true)
164+
extendlow = colormap[1]
165+
extendhigh = colormap[end]
166+
colormap = cgrad(colormap[2:end-1]; categorical=true)
167+
168+
fig = Figure(size = (1000, 250length(ϕsnorth)), fontsize = 18)
169+
axs = Array{Any, 2}(undef, (length(ϕsnorth), length(basin_keys)))
170+
contours = Array{Any, 2}(undef, (length(ϕsnorth), length(basin_keys)))
171+
for (icol, (basin_key, basin)) in enumerate(pairs(basins))
172+
173+
for (irow, (x3D, time_window, experiment)) in enumerate(zip(ϕsnorth, time_windows, experiments2))
174+
175+
x2D = dropdims(reverse(nancumsum(reverse(nansum(basin .* x3D, dims = 1), dims=3), dims = 3), dims=3), dims = 1) # kg/s
176+
x2Dmask = zonalaverage(1, gridmetrics; mask = basin) .> 0
177+
x2D[.!x2Dmask] .= NaN
178+
179+
# convert to Sv
180+
x2D = x2D / 1e6 / ρ # Sv
181+
182+
local ax = Axis(fig[irow, icol],
183+
backgroundcolor = :lightgray,
184+
xgridvisible = true, ygridvisible = true,
185+
xgridcolor = (:black, 0.05), ygridcolor = (:black, 0.05),
186+
ylabel = "depth (m)")
187+
188+
X = dropdims(maximum(lat, dims=1), dims=1)
189+
Y = zt
190+
Z = x2D
191+
co = contourf!(ax, X, Y, Z;
192+
levels,
193+
colormap,
194+
nan_color = :lightgray,
195+
extendlow,
196+
extendhigh,
197+
)
198+
translate!(co, 0, 0, -100)
199+
contours[irow, icol] = co
200+
201+
xlim = basin_latlims[basin_key]
202+
# basin2 = LONGTEXT[basin]
203+
204+
ax.yticks = (ztick, zticklabel)
205+
xticks = -90:30:90
206+
ax.xticks = (xticks, latticklabel.(xticks))
207+
ylims!(ax, zlim)
208+
# xlims!(ax, (-90, 90))
209+
xlims!(ax, xlim)
210+
211+
212+
hidexdecorations!(ax,
213+
label = irow < size(axs, 1), ticklabels = irow < size(axs, 1),
214+
ticks = irow < size(axs, 1), grid = false)
215+
hideydecorations!(ax,
216+
label = icol > 1, ticklabels = icol > 1,
217+
ticks = icol > 1, grid = false)
218+
219+
(icol == 1) && Label(fig[irow, 0]; text = "$experiment $(time_window[4:7])s", rotation = π/2, tellheight = false)
220+
221+
axs[irow, icol] = ax
222+
end
223+
224+
Label(fig[0, icol], "$(basin_strs[icol]) MOC", tellwidth = false)
225+
226+
227+
end
228+
229+
cb = Colorbar(fig[1:size(axs, 1), length(basin_keys) + 1], contours[1, 1];
230+
vertical = true, flipaxis = true,
231+
# ticks = (, cbarticklabelformat.(levels)),
232+
tickformat = x -> map(t -> replace(format("{:+d}", t), "-" => "−"), x),
233+
label = "MOC (Sv)",
234+
)
235+
cb.height = Relative(0.666)
236+
237+
for (ibasin, xlims) in enumerate(basin_latlims)
238+
# Label(f[0, ibasin], LONGTEXT[basin], fontsize=20, tellwidth=false)
239+
colsize!(fig.layout, ibasin, Auto(xlims[2] - xlims[1]))
240+
end
241+
242+
rowgap!(fig.layout, 15)
243+
# rowgap!(fig.layout, 4, 15)
244+
# # # rowgap!(fig.layout, 5, 10)
245+
colgap!(fig.layout, 15)
246+
# title = "$model ensemble-mean MOC"
247+
# Label(fig[-1, 1:3]; text = title, fontsize = 20, tellwidth = false)
248+
249+
labels = [
250+
"a" "b" "c"
251+
"d" "e" "f"
252+
"g" "h" "i"
253+
]
254+
255+
labeloptions = (
256+
font = :bold,
257+
align = (:left, :bottom),
258+
offset = (5, 2),
259+
space = :relative,
260+
fontsize = 24
261+
)
262+
263+
for (ax, label) in zip(axs, labels)
264+
txt = text!(ax, 0, 0; text = label, labeloptions..., strokecolor = :white, strokewidth = 3)
265+
translate!(txt, 0, 0, 100)
266+
txt= text!(ax, 0, 0; text = label, labeloptions...)
267+
translate!(txt, 0, 0, 100)
268+
end
269+
270+
fig
271+
272+
outputdir = "output/plots"
273+
274+
# save plot
275+
outputfile = joinpath(outputdir, "$(model)_MOC_change_for_Reviewer2.png")
276+
@info "Saving MOC as image file:\n $(outputfile)"
277+
save(outputfile, fig)
278+
outputfile = joinpath(outputdir, "$(model)_MOC_change_for_Reviewer2.pdf")
279+
@info "Saving MOC as image file:\n $(outputfile)"
280+
save(outputfile, fig)
281+
282+
283+
284+
285+

‎src/plot_MOCs_withGM.jl‎

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,9 @@ using Interpolations
1717
using OceanBasins
1818
using NaNStatistics
1919
using Format
20+
using GeometryBasics
21+
using LibGEOS
22+
using GeometryOps
2023

2124
include("plotting_functions.jl")
2225

@@ -26,10 +29,12 @@ model = "ACCESS-ESM1-5"
2629

2730
CMIP_version = model == "ACCESS1-3" ? "CMIP5" : "CMIP6"
2831

29-
experiment = "historical"
32+
experiment = "ssp370"
33+
# experiment = "historical"
3034
# experiment = "piControl"
3135

32-
time_window = "Jan1990-Dec1999"
36+
# time_window = "Jan1990-Dec1999"
37+
time_window = "Jan2030-Dec2090"
3338
# time_window = "Jan1071-Dec1100" # <- last 30 years of ACCESS-ESM1-5 piControl
3439
# time_window = "Jan1420-Dec1449" # <- last 30 years of ACCESS-CM2 piControl
3540

@@ -48,7 +53,7 @@ requiredvariables = ["umo", "vmo", "mlotst", "volcello", "areacello", "agessc"]
4853
hasrequireddata(member, variable_name) = isfile(joinpath(inputdirfun(member), "$variable_name.nc"))
4954
hasrequireddata(member) = all(variable_name -> hasrequireddata(member, variable_name), requiredvariables)
5055
members = readdir(joinpath(DATADIR, "$model/$experiment"))
51-
members = [m for m in members if m ≠ ".DS_Store"]
56+
members = [m for m in members if m ≠ ".DS_Store" && startswith(m, "r")]
5257

5358
# sort members by r, i, p[, f]
5459
member_regex = CMIP_version == "CMIP6" ? r"r(\d+)i(\d+)p(\d+)f(\d+)" : r"r(\d+)i(\d+)p(\d+)"
@@ -68,8 +73,9 @@ println()
6873
# for member in members[dataavailability.has_it_all]
6974
# for member in [last(members)]
7075
# member = first(members)
71-
member = "r5i1p1f1" # <- presumably the ACCESS-ESM1-5 member which has GM transport available ("HI-05")
76+
# member = "r5i1p1f1" # <- presumably the ACCESS-ESM1-5 member which has GM transport available ("HI-05")
7277
# member = "r1i1p1" # <- presumably the ACCESS1-3 member that Matt Chamberlain used for GM transport
78+
member = "r1i1p1f1"
7379

7480
inputdir = inputdirfun(member)
7581

@@ -113,17 +119,12 @@ member = "r5i1p1f1" # <- presumably the ACCESS-ESM1-5 member which has GM transp
113119
# Some parameter values
114120
ρ = 1035.0 # density (kg/m^3)
115121

116-
# unpack model grid
117-
(; lon, lat, zt, v3D,) = gridmetrics
118-
lev = zt
119-
# unpack indices
120-
(; wet3D, N) = indices
121122

122123

123124

124125
# basins
125126
basin_keys = (:ATL, :INDOPAC, :GBL)
126-
basin_strs = ("Atlantic", "IndoPacific", "Global")
127+
basin_strs = ("Atlantic", "Indo-Pacific", "Global")
127128
isatlnoSO(lat, lon, o) = isatlantic(lat, lon, o) .& (lat .≥ -30)
128129
isindopacific(lat, lon, o) = (isindian(lat, lon, o) .| ispacific(lat, lon, o)) .& (lat .≥ -30)
129130
isglobal(lat, lon, o) = trues(size(lat))

0 commit comments

Comments
 (0)