update all plots for submission

Author: briochemc

src/average_monthly_seqeff_and_adjprop.jl

@@ -1,4 +1,4 @@
-# qsub -I -P xv83 -l mem=64GB -l storage=scratch/gh0+scratch/xv83 -l walltime=01:00:00 -l ncpus=12
+# qsub -I -P xv83 -l mem=47GB -l storage=scratch/gh0+scratch/xv83 -l walltime=01:00:00 -l ncpus=12
 
 using Pkg
 Pkg.activate(".")
@@ -37,18 +37,13 @@ using NaNStatistics
 # script options
 @show model = "ACCESS-ESM1-5"
 if isempty(ARGS)
-    member = "r1i1p1f1"
-    # experiment = "historical"
-    # time_window = "Jan1850-Dec1859"
-    # time_window = "Jan1990-Dec1999"
     experiment = "ssp370"
+    member = "r1i1p1f1"
     time_window = "Jan2030-Dec2039"
     # time_window = "Jan2090-Dec2099"
-    WRITEDATA = "true"
 else
-    experiment, member, time_window, finalmonth, WRITEDATA = ARGS
+    experiment, member, time_window = ARGS
 end
-WRITEDATA = parse(Bool, WRITEDATA)
 @show experiment
 @show member
 @show time_window
@@ -70,9 +65,9 @@ volcello_ds = open_dataset(joinpath(fixedvarsinputdir, "volcello.nc"))
 areacello_ds = open_dataset(joinpath(fixedvarsinputdir, "areacello.nc"))
 
 # members = ["r$(i)i1p1f1" for i in 1:3]
-members = ["r$(i)i1p1f1" for i in 1:4]
+# members = ["r$(i)i1p1f1" for i in 1:1]
 
-for member in members
+# for member in members
 
     # Gadi directory for input files
     lumpby = "month"
@@ -154,7 +149,8 @@ for member in members
     # Save to netCDF file
     outputfile = joinpath(inputdir, "calE$(upwind_str)_$(κVdeep_str)_$(κH_str)_$(κVML_str).nc")
     @info "Saving mean sequestration efficiency as netCDF file:\n  $(outputfile)"
-    savedataset(ds, path = outputfile, driver = :netcdf, overwrite = true)
+    ds_chunked = setchunks(ds, (x = 60, y = 60, Ti = length(ds.Ti)))
+    savedataset(ds_chunked, path = outputfile, driver = :netcdf, overwrite = true)
 
 
     # Do the same for 𝒢̃
@@ -190,8 +186,9 @@ for member in members
     # Save to netCDF file
     outputfile = joinpath(inputdir, "calgtilde$(upwind_str)_$(κVdeep_str)_$(κH_str)_$(κVML_str).nc")
     @info "Saving mean adjoint propagator as netCDF file:\n  $(outputfile)"
-    savedataset(ds, path = outputfile, driver = :netcdf, overwrite = true)
+    ds_chunked = setchunks(ds, (x = 60, y = 60, Ti = length(ds.Ti)))
+    savedataset(ds_chunked, path = outputfile, driver = :netcdf, overwrite = true)
 
 
-end
+# end
 

src/average_monthly_seqeff_and_adjprop.sh

@@ -0,0 +1,23 @@
+#!/bin/bash
+
+#PBS -P xv83
+#PBS -N avg_member_placeholder_time_window_placeholder
+#PBS -l ncpus=16
+#PBS -q express
+#PBS -l mem=63GB
+#PBS -l jobfs=4GB
+#PBS -l walltime=0:10:00
+#PBS -l storage=scratch/gh0+scratch/xv83
+#PBS -l wd
+#PBS -o output/PBS/
+#PBS -j oe
+
+experiment=ssp370
+time_window=time_window_placeholder
+member=member_placeholder
+
+echo "Going into ACCESS-TMIP"
+cd ~/Projects/TMIP/ACCESS-TMIP
+
+echo "Running script"
+julia src/average_monthly_seqeff_and_adjprop.jl $experiment $member $time_window &> output/average_monthly_seqeff_and_adjprop.$experiment.$member.$time_window.$PBS_JOBID.out

src/average_monthly_seqeff_and_adjprop_multijob.sh

@@ -0,0 +1,10 @@
+for member in r{14..14}i1p1f1; do
+    for time_window in Jan2030-Dec2039 Jan2090-Dec2099; do
+        # change member_placeholders in timestep_adjoint_propagator.sh and create file and remove it
+        sed "s/member_placeholder/$member/g" src/average_monthly_seqeff_and_adjprop.sh > src/tmp.sh
+        sed -i "s/time_window_placeholder/$time_window/g" src/tmp.sh
+        qsub src/tmp.sh
+        rm src/tmp.sh
+    done
+done
+

src/plot_AA_age.jl

@@ -1,4 +1,4 @@
-# qsub -I -P xv83 -l mem=64GB -l storage=scratch/gh0+scratch/xv83 -l walltime=00:20:00 -l ncpus=12
+# qsub -I -P xv83 -q express -l mem=47GB -l storage=scratch/gh0+scratch/xv83 -l walltime=00:30:00 -l ncpus=12
 
 using Pkg
 Pkg.activate(".")

src/plot_ACCESS_quantilesandmean.jl

@@ -1,4 +1,4 @@
-# qsub -I -P xv83 -l mem=47GB -q express -l storage=scratch/gh0+scratch/xv83 -l walltime=01:00:00 -l ncpus=12
+# qsub -I -P xv83 -q express -l mem=47GB -l storage=scratch/gh0+scratch/xv83 -l walltime=01:00:00 -l ncpus=12
 
 using Pkg
 Pkg.activate(".")
@@ -72,8 +72,6 @@ lev = zt
 indices = makeindices(gridmetrics.v3D)
 (; wet3D, N) = indices
 
-
-
 # Preferred diffusivities
 κVdeep = 3.0e-5 # m^2/s
 κVML = 1.0      # m^2/s
@@ -86,7 +84,7 @@ upwind_str = upwind ? "" : "_centered"
 upwind_str2 = upwind ? "upwind" : "centered"
 
 # Use yearly time-stepped simulations or monthly ones?
-yearly = true
+yearly = false
 yearly_str = yearly ? "_yearly" : ""
 yearly_str2 = yearly ? "(yearly)" : ""
 
@@ -99,34 +97,29 @@ end
 
 # To avoid loading and carrying 100s of GB of data around,
 # preprocess each member before and only save the 2D data needed for plots.
-if yearly
-    ℰ_file0 = "/scratch/xv83/TMIP/data/$model/$experiment/$(first(members))/$(time_window)/seqeff$(upwind_str)_$(κVdeep_str)_$(κH_str)_$(κVML_str)$(yearly_str).nc"
-    ℰ_ds0 = open_dataset(ℰ_file0)
-    years = ℰ_ds0.Ti |> Array
-    ℰ1050_ensemble = reduce((x,y) -> cat(x, y, dims = 4), map(members) do member
-        @info "loading $member ℰ"
-        ℰ_file = "/scratch/xv83/TMIP/data/$model/$experiment/$member/$(time_window)/seqeff$(upwind_str)_$(κVdeep_str)_$(κH_str)_$(κVML_str)$(yearly_str).nc"
-        ℰ_ds = open_dataset(ℰ_file)
-        ℰ = readcubedata(ℰ_ds.seqeff)
-        ℰ10 = map(
-            ts -> yearatquantile(ts, 0.9),
-            view(ℰ, i, j, :) for i in 1:size(ℰ,1), j in 1:size(ℰ,2)
-        )
-        ℰ50 = map(
-            ts -> yearatquantile(ts, 0.5),
-            view(ℰ, i, j, :) for i in 1:size(ℰ,1), j in 1:size(ℰ,2)
-        )
-        [ℰ10;;; ℰ50]
-    end)
-    ℰ1050_ensemblemean = dropdims(mean(ℰ1050_ensemble, dims = 4), dims = 4)
-    ℰ1050_ensemblemax = dropdims(maximum(ℰ1050_ensemble, dims = 4), dims = 4)
-    ℰ1050_ensemblemin = dropdims(minimum(ℰ1050_ensemble, dims = 4), dims = 4)
-    ℰ1050_ensemblerange = ℰ1050_ensemblemax - ℰ1050_ensemblemin
-
-
-else
-    # TODO (don't forget to deal with months)
-end
+varname = yearly ? "seqeff" : "calE"
+ℰ_file0 = "/scratch/xv83/TMIP/data/$model/$experiment/$(first(members))/$(time_window)/$(varname)$(upwind_str)_$(κVdeep_str)_$(κH_str)_$(κVML_str)$(yearly_str).nc"
+ℰ_ds0 = open_dataset(ℰ_file0)
+years = ℰ_ds0.Ti |> Array
+τℰ1050_ensemble = reduce((x,y) -> cat(x, y, dims = 4), map(members) do member
+    @info "loading $member ℰ"
+    ℰ_file = "/scratch/xv83/TMIP/data/$model/$experiment/$member/$(time_window)/$(varname)$(upwind_str)_$(κVdeep_str)_$(κH_str)_$(κVML_str)$(yearly_str).nc"
+    ℰ_ds = open_dataset(ℰ_file)
+    ℰ = readcubedata(ℰ_ds[varname])
+    τℰ10 = map(
+        ts -> yearatquantile(ts, 0.9),
+        view(ℰ, i, j, :) for i in 1:size(ℰ,1), j in 1:size(ℰ,2)
+    )
+    τℰ50 = map(
+        ts -> yearatquantile(ts, 0.5),
+        view(ℰ, i, j, :) for i in 1:size(ℰ,1), j in 1:size(ℰ,2)
+    )
+    [τℰ10;;; τℰ50]
+end)
+τℰ1050_ensemblemean = dropdims(mean(τℰ1050_ensemble, dims = 4), dims = 4)
+τℰ1050_ensemblemax = dropdims(maximum(τℰ1050_ensemble, dims = 4), dims = 4)
+τℰ1050_ensemblemin = dropdims(minimum(τℰ1050_ensemble, dims = 4), dims = 4)
+τℰ1050_ensemblerange = τℰ1050_ensemblemax - τℰ1050_ensemblemin
 
 include("plotting_functions.jl") # load seafloorvalue function
 
@@ -157,16 +150,18 @@ fig = Figure(size = (ncols * 500, nrows * 250 + 100), fontsize = 18)
 yticks = -60:30:60
 xticks = -120:60:120 + 360
 
-datamean = (Γout_ensemblemean, ℰ1050_ensemblemean[:,:,2], ℰ1050_ensemblemean[:,:,1])
-datarange = (Γout_ensemblerange, ℰ1050_ensemblerange[:,:,2], ℰ1050_ensemblerange[:,:,1])
+datamean = (Γout_ensemblemean, τℰ1050_ensemblemean[:,:,2], τℰ1050_ensemblemean[:,:,1])
+datarange = (Γout_ensemblerange, τℰ1050_ensemblerange[:,:,2], τℰ1050_ensemblerange[:,:,1])
 𝒓 = rich("r", font = :bold_italic)
 Γstr = rich("Γ", superscript("†"), rich("‾", offset = (-0.55, 0.25)), rich("‾", offset = (-0.85, 0.25)))
 Γfun = rich(Γstr, rich("(", 𝒓, ")", offset = (0.4, 0)))
 # Qstr = rich("Q", font = :italic)
 ℰstr = rich("ℰ", rich("‾", offset = (-0.5, 0.15)))
 # Q10 = rich(Qstr, "(0.1)")
 # Q50 = rich(Qstr, "(0.5)")
-rowlabels = (rich("Mean time, ", Γfun), rich("Median, ", ℰstr, " = 50 %"), rich("10th percentile, ", ℰstr, " = 90 %"))
+# rowlabels = (rich("Mean time, ", Γfun), rich("Median time (", ℰstr, " = 50 %)"), rich("10th percentile time (", ℰstr, " = 90 %)"))
+rowlabels = (rich("Mean time, ", Γfun), rich("Median time (", ℰstr, " = 50 %)"), rich("10th %ile time (", ℰstr, " = 90 %)"))
+# rowlabels = (rich("Mean time, ", Γfun), rich("Median time (", ℰstr, " = 50 %)"), rich("10 % time (", ℰstr, " = 90 %)"))
 
 
 for (irow, (x2Dmean, x2Drange, text)) in enumerate(zip(datamean, datarange, rowlabels))
@@ -179,7 +174,7 @@ for (irow, (x2Dmean, x2Drange, text)) in enumerate(zip(datamean, datarange, rowl
     colormap = cgrad(colormap[1:end-1], categorical = true)
     colorrange = extrema(levels)
 
-    axs[irow, icol] = ax = Axis(fig[irow, icol]; yticks, xticks, xtickformat, ytickformat)
+    axs[irow, icol] = ax = Axis(fig[irow, icol]; yticks, xticks, xtickformat, ytickformat, aspect = DataAspect())
 
     contours[irow, icol] = if usecontourf
         plotcontourfmap!(ax, x2Dmean, gridmetrics; levels, colormap)
@@ -198,7 +193,7 @@ for (irow, (x2Dmean, x2Drange, text)) in enumerate(zip(datamean, datarange, rowl
     colormap = cgrad(colormap[1:end-1], categorical = true)
     colorrange = extrema(levels)
 
-    axs[irow, icol] = ax = Axis(fig[irow, icol]; yticks, xticks, xtickformat, ytickformat)
+    axs[irow, icol] = ax = Axis(fig[irow, icol]; yticks, xticks, xtickformat, ytickformat, aspect = DataAspect())
 
     contours[irow, icol] = if usecontourf
         plotcontourfmap!(ax, x2Drange, gridmetrics; levels, colormap)
@@ -214,7 +209,7 @@ for (irow, (x2Dmean, x2Drange, text)) in enumerate(zip(datamean, datarange, rowl
 end
 
 
-label = rich("ensemble mean (years)")
+label = rich("ensemble mean $(time_window[4:7])s characteristic timescales (years)")
 cb = Colorbar(fig[nrows + 1, 1], contours[1, 1]; label, vertical = false, flipaxis = false, ticks = 0:1000:3000)
 cb.width = Relative(2/3)
 
@@ -240,12 +235,18 @@ for (ax, label) in zip(axs, labels)
 end
 
 # Label(fig[0, 1:2]; text = "$(time_window[4:7])s Seafloor Reemergence Time ($(length(members)) members)", fontsize = 24, tellwidth = false)
-Label(fig[0, 1:2]; text = "$(time_window[4:7])s Characteristic Timescales of Reemergence ($(length(members)) members)$(yearly_str2)", fontsize = 24, tellwidth = false)
+# Label(fig[0, 1:2]; text = "$(time_window[4:7])s Characteristic Timescales of Reemergence ($(length(members)) members)$(yearly_str2)", fontsize = 24, tellwidth = false)
 rowgap!(fig.layout, 10)
 colgap!(fig.layout, 10)
 
+colsize!(fig.layout, 1, Aspect(1, 2.0))
+colsize!(fig.layout, 2, Aspect(1, 2.0))
+resize_to_layout!(fig)
+
 # save plot
 suffix = usecontourf ? "_ctrf" : ""
+
+
 outputfile = joinpath(outputdir, "reemergencetime$(upwind_str)_$(κVdeep_str)_$(κH_str)_$(κVML_str)$(yearly_str)_$(time_window)$(suffix).png")
 @info "Saving reemergencetime on sea floor as image file:\n  $(outputfile)"
 save(outputfile, fig)
@@ -255,7 +256,7 @@ save(outputfile, fig)
 
 
 
-# datamean = (Γout_ensemblemean, ℰ50_ensemblemean, ℰ10_ensemblemean)
+# datamean = (Γout_ensemblemean, τℰ50_ensemblemean, τℰ10_ensemblemean)
 ikeep = .!isnan.(datamean[1]) .& (seafloorvalue(Z3D, wet3D) .> 3000)
 data = datamean[2][ikeep] ./ datamean[1][ikeep] .- 1
 weights = Weights(gridmetrics.area2D[ikeep])