ACCESS-OM3 depth-latitude temperature and temperature difference from WOA23¶

Related issue: https://github.com/ACCESS-Community-Hub/access-om3-25km-paper-1/issues/11

OM2 paper figure 12 https://github.com/COSIMA/ACCESS-OM2-1-025-010deg-report/blob/master/figures/temperature_biases/zonal_temp_salt_bias.ipynb

In [1]:

parameters

# These first two cells must be in all notebooks!
# It allows us to run all the notebooks at once, this cell has a tag "parameters" which allows us to pass in 
# arguments externally using papermill (see mkfigs.sh for details)

# Set esm_file to the datastore for the main experiment of interest
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json"

# papermill settings. *No need to modify these if running interactively.* 
papermill = False                      # `cwd` and `nbname` will be populated by papermill.
cwd = None                             # current working directory 
nbname = None                          # notebook name

# These first two cells must be in all notebooks! # It allows us to run all the notebooks at once, this cell has a tag "parameters" which allows us to pass in # arguments externally using papermill (see mkfigs.sh for details) # Set esm_file to the datastore for the main experiment of interest esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json" # papermill settings. *No need to modify these if running interactively.* papermill = False # `cwd` and `nbname` will be populated by papermill. cwd = None # current working directory nbname = None # notebook name

In [2]:

injected-parameters

# Parameters
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/"
nbname = "temp-salt-vs-depth-latitude.ipynb"

# Parameters esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/datastore.json" papermill = True cwd = "/g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/" nbname = "temp-salt-vs-depth-latitude.ipynb"

In [3]:

if not papermill: 
    import nci_ipynb, os  # requires conda/analysis3-26.03 or later
    cwd = nci_ipynb.dir()
    nbname = nci_ipynb.name()
    os.chdir(cwd)
import mkfigs_bootstrap  # noqa: adds external/access-model-mkfigs/src to sys.path (stop-gap)
from mkfigs import MkmdWriter
mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill)

if not papermill: import nci_ipynb, os # requires conda/analysis3-26.03 or later cwd = nci_ipynb.dir() nbname = nci_ipynb.name() os.chdir(cwd) import mkfigs_bootstrap # noqa: adds external/access-model-mkfigs/src to sys.path (stop-gap) from mkfigs import MkmdWriter mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill)

In [4]:

IAF = esm_file.find('iaf') > 0
IAF

IAF = esm_file.find('iaf') > 0 IAF

Out[4]:

True

In [5]:

import os
import numpy as np
import xarray as xr
import cftime
import cf_xarray as cfxr
import cf_xarray.units
import pint_xarray
from pint import application_registry as ureg
import gsw
import intake
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import matplotlib as mpl
from distributed import Client
import cmocean
import cartopy.crs as ccrs
import cartopy.feature as cft
from shapely import geometry
from textwrap import wrap
from tqdm.notebook import tqdm
xr.set_options(keep_attrs=True);  # cf_xarray works best when xarray keeps attributes by default

import os import numpy as np import xarray as xr import cftime import cf_xarray as cfxr import cf_xarray.units import pint_xarray from pint import application_registry as ureg import gsw import intake import matplotlib.pyplot as plt from matplotlib.patches import Patch import matplotlib as mpl from distributed import Client import cmocean import cartopy.crs as ccrs import cartopy.feature as cft from shapely import geometry from textwrap import wrap from tqdm.notebook import tqdm xr.set_options(keep_attrs=True); # cf_xarray works best when xarray keeps attributes by default

In [6]:

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

import warnings warnings.simplefilter(action='ignore', category=FutureWarning)

In [7]:

client = Client(threads_per_worker=1)
client

client = Client(threads_per_worker=1) client

Out[7]:

Client

Client-c8dac90b-8670-11f1-89b5-000001f7fe80

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Launch dashboard in JupyterLab

Cluster Info

LocalCluster

bb8afcc8

Dashboard: http://127.0.0.1:8787/status	Workers: 16
Total threads: 16	Total memory: 188.56 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-36afb532-5509-4588-9736-36b09f42f922

Comm: tcp://127.0.0.1:38349	Workers: 0
Dashboard: http://127.0.0.1:8787/status	Total threads: 0
Started: Just now	Total memory: 0 B

Workers

Worker: 0

Comm: tcp://127.0.0.1:35905	Total threads: 1
Dashboard: http://127.0.0.1:38293/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:44705
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-u52j7rkg

Worker: 1

Comm: tcp://127.0.0.1:34959	Total threads: 1
Dashboard: http://127.0.0.1:42743/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:44205
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-uzk06owk

Worker: 2

Comm: tcp://127.0.0.1:43639	Total threads: 1
Dashboard: http://127.0.0.1:43687/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:35697
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-7uxmj1pj

Worker: 3

Comm: tcp://127.0.0.1:38931	Total threads: 1
Dashboard: http://127.0.0.1:44535/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:40807
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-sgngb9qv

Worker: 4

Comm: tcp://127.0.0.1:35969	Total threads: 1
Dashboard: http://127.0.0.1:39877/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:36417
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-x771my_1

Worker: 5

Comm: tcp://127.0.0.1:42315	Total threads: 1
Dashboard: http://127.0.0.1:42205/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:33571
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-_wktxl_3

Worker: 6

Comm: tcp://127.0.0.1:34077	Total threads: 1
Dashboard: http://127.0.0.1:34181/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:46571
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-qn_akcsm

Worker: 7

Comm: tcp://127.0.0.1:44179	Total threads: 1
Dashboard: http://127.0.0.1:42843/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:45477
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-uy7uj6ky

Worker: 8

Comm: tcp://127.0.0.1:41397	Total threads: 1
Dashboard: http://127.0.0.1:39193/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:38613
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-5q23i0id

Worker: 9

Comm: tcp://127.0.0.1:33205	Total threads: 1
Dashboard: http://127.0.0.1:45499/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:44105
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-yrm5dsd2

Worker: 10

Comm: tcp://127.0.0.1:44565	Total threads: 1
Dashboard: http://127.0.0.1:37405/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:44865
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-63iq1uv9

Worker: 11

Comm: tcp://127.0.0.1:44141	Total threads: 1
Dashboard: http://127.0.0.1:42371/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:38331
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-g9usjepf

Worker: 12

Comm: tcp://127.0.0.1:36587	Total threads: 1
Dashboard: http://127.0.0.1:40541/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:44569
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-7lk5znax

Worker: 13

Comm: tcp://127.0.0.1:36603	Total threads: 1
Dashboard: http://127.0.0.1:37829/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:44549
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-clwkv45v

Worker: 14

Comm: tcp://127.0.0.1:36401	Total threads: 1
Dashboard: http://127.0.0.1:42477/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:39917
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-rev6yk4n

Worker: 15

Comm: tcp://127.0.0.1:40837	Total threads: 1
Dashboard: http://127.0.0.1:43771/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:41871
Local directory: /jobfs/174557274.gadi-pbs/dask-scratch-space/worker-1kzohvt9

Define regions¶

In [8]:

regions = { # [minx, maxx, miny, maxy], using model longitude range (-280 to 80)
    "Global": [-280, 80, -90, 90],
    # "Arctic": [-280, 80, 65, 90],
    # "Southern Ocean": [-280, 80, -82, -63],
    # "ACC": [-280, 80, -63, -45],
    # "Southern Pacific": [-210, -70, -45, -20],
    # "Tropical Pacific": [-240, -100, -20, 20],
    # "North Pacific": [-240, -100, 20, 65],
    # "South Atlantic": [-60, 20, -45, -20],
    # "Tropical Atlantic": [-70, 20, -20, 20],
    # "North Atlantic": [-100, 0, 20, 65],
    # "Indian": [30, 120, -45, 20],
    # "Aegean Sea": [18, 27.5, 34, 44],
    # "Black Sea": [27.5, 43, 40.5, 48],
    # "Baltic Sea": [13, 30, 53, 58],
    # "Mediterranean Sea": [0, 35, 31, 41],
    # "Red Sea": [33, 44, 12, 29],
    # "Persian Gulf": [47, 56, 24, 31],
    # "White Sea": [31, 41, 63, 68],
}
regions = {k: dict(zip(["minx", "maxx", "miny", "maxy"], v)) for k, v in regions.items()}  # convert to dicts

for r, d in regions.items():
    for k, x in d.items():
        if k in ["minx", "maxx"] and x != max(-280, min(x, 80)):
            raise ValueError(f"{r} {k} = {x} is outside the range -280 to 80")
        if k in ["miny", "maxy"] and x > 65:
            print(f"{r} {k} changed from {x} to 65 to omit tripolar region")
            d[k] = 65
            

fig = plt.figure(figsize=(10, 4))

colors = mpl.color_sequences['tab20']
try:  # set Global color to white 
    colors[list(regions.keys()).index("Global")] = (1, 1, 1)
except ValueError:
    pass

ax = plt.axes(position=[0.05,0.05,0.7,0.9], projection=ccrs.Robinson())
ax.coastlines(resolution="110m")
ax.gridlines(draw_labels=False)

legend_elements = []
for i, (region, limits) in enumerate(regions.items()):
    ax.add_geometries([geometry.box(**limits)], crs=ccrs.PlateCarree(), color=colors[i], alpha=0.3)
    legend_elements.append(Patch(color=colors[i], alpha=0.3, label=region))

ax.legend(handles=legend_elements, bbox_to_anchor=(1.2, 1.00))
plt.title("Regions (BUG: incorrect shapes north of 65°N)")

regions = { # [minx, maxx, miny, maxy], using model longitude range (-280 to 80) "Global": [-280, 80, -90, 90], # "Arctic": [-280, 80, 65, 90], # "Southern Ocean": [-280, 80, -82, -63], # "ACC": [-280, 80, -63, -45], # "Southern Pacific": [-210, -70, -45, -20], # "Tropical Pacific": [-240, -100, -20, 20], # "North Pacific": [-240, -100, 20, 65], # "South Atlantic": [-60, 20, -45, -20], # "Tropical Atlantic": [-70, 20, -20, 20], # "North Atlantic": [-100, 0, 20, 65], # "Indian": [30, 120, -45, 20], # "Aegean Sea": [18, 27.5, 34, 44], # "Black Sea": [27.5, 43, 40.5, 48], # "Baltic Sea": [13, 30, 53, 58], # "Mediterranean Sea": [0, 35, 31, 41], # "Red Sea": [33, 44, 12, 29], # "Persian Gulf": [47, 56, 24, 31], # "White Sea": [31, 41, 63, 68], } regions = {k: dict(zip(["minx", "maxx", "miny", "maxy"], v)) for k, v in regions.items()} # convert to dicts for r, d in regions.items(): for k, x in d.items(): if k in ["minx", "maxx"] and x != max(-280, min(x, 80)): raise ValueError(f"{r} {k} = {x} is outside the range -280 to 80") if k in ["miny", "maxy"] and x > 65: print(f"{r} {k} changed from {x} to 65 to omit tripolar region") d[k] = 65 fig = plt.figure(figsize=(10, 4)) colors = mpl.color_sequences['tab20'] try: # set Global color to white colors[list(regions.keys()).index("Global")] = (1, 1, 1) except ValueError: pass ax = plt.axes(position=[0.05,0.05,0.7,0.9], projection=ccrs.Robinson()) ax.coastlines(resolution="110m") ax.gridlines(draw_labels=False) legend_elements = [] for i, (region, limits) in enumerate(regions.items()): ax.add_geometries([geometry.box(**limits)], crs=ccrs.PlateCarree(), color=colors[i], alpha=0.3) legend_elements.append(Patch(color=colors[i], alpha=0.3, label=region)) ax.legend(handles=legend_elements, bbox_to_anchor=(1.2, 1.00)) plt.title("Regions (BUG: incorrect shapes north of 65°N)")

Global maxy changed from 90 to 65 to omit tripolar region

Out[8]:

Text(0.5, 1.0, 'Regions (BUG: incorrect shapes north of 65°N)')

Load data from ACCESS-OM3¶

In [9]:

om3exptname = os.path.basename(os.path.dirname(esm_file))
catalogs = [
                esm_file,
                # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm1-d968c801/datastore.json',
                # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm2-5dc49da6/datastore.json',
                # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm3-da330542/datastore.json',
                # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm4-9fd08880/datastore.json',
                # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm5-9b5dbfa9/datastore.json',
                # '/g/data/ol01/access-om3-output/access-om3-025/25km-iaf-test-for-AK-expt-7df5ef4c/datastore.json',
            ]
catalogs

om3exptname = os.path.basename(os.path.dirname(esm_file)) catalogs = [ esm_file, # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm1-d968c801/datastore.json', # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm2-5dc49da6/datastore.json', # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm3-da330542/datastore.json', # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm4-9fd08880/datastore.json', # '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm5-9b5dbfa9/datastore.json', # '/g/data/ol01/access-om3-output/access-om3-025/25km-iaf-test-for-AK-expt-7df5ef4c/datastore.json', ] catalogs

Out[9]:

['/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/datastore.json']

In [10]:

om3datastores = { os.path.normpath(c).split(os.sep)[-2]:
               intake.open_esm_datastore(c,
                                         columns_with_iterables=[
                                            "variable",
                                            "variable_long_name",
                                            "variable_standard_name",
                                            "variable_cell_methods",
                                            "variable_units"]
                                        )
              for c in catalogs }
om3datastores

om3datastores = { os.path.normpath(c).split(os.sep)[-2]: intake.open_esm_datastore(c, columns_with_iterables=[ "variable", "variable_long_name", "variable_standard_name", "variable_cell_methods", "variable_units"] ) for c in catalogs } om3datastores

Out[10]:

{'MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36': <datastore catalog with 16 dataset(s) from 2190 asset(s)>}

In [11]:

om3varnames = [ 'thetao', 'so' ]  # potential temperature and practical salinity

om3vars = {
    vname: {expt: ds.search(variable=vname).to_dask(
                xarray_open_kwargs = dict(
                    chunks={"xh": -1, "yh": "auto", "z_l": "auto"},
                    decode_timedelta=True
                )
            )[vname]
            for expt, ds in om3datastores.items()
          }
    for vname in om3varnames
    }

# omit latitudes with grid bug in this run https://github.com/ACCESS-NRI/ocean_model_grid_generator/issues/7
for vname, d in om3vars.items():
    try:
        d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yh=slice(10, None))
    except ValueError:
        d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yq=slice(10, None))
    except KeyError:
        pass

if IAF:
    for vname, d in om3vars.items():
        for expt in d:
            try:
                d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True)
            except KeyError:
                pass

om3varnames = [ 'thetao', 'so' ] # potential temperature and practical salinity om3vars = { vname: {expt: ds.search(variable=vname).to_dask( xarray_open_kwargs = dict( chunks={"xh": -1, "yh": "auto", "z_l": "auto"}, decode_timedelta=True ) )[vname] for expt, ds in om3datastores.items() } for vname in om3varnames } # omit latitudes with grid bug in this run https://github.com/ACCESS-NRI/ocean_model_grid_generator/issues/7 for vname, d in om3vars.items(): try: d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yh=slice(10, None)) except ValueError: d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yq=slice(10, None)) except KeyError: pass if IAF: for vname, d in om3vars.items(): for expt in d: try: d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True) except KeyError: pass

Load data from ACCESS-OM2¶

use control - see https://forum.access-hive.org.au/t/access-om2-control-experiments/258#p-747-ryf-7

In [12]:

if IAF:
    om2exptname = '025deg_jra55_iaf_omip2_cycle1'
else:
    om2exptname = '025deg_jra55_ryf9091_gadi'  # monthly temp only for 1900-1903, and 2300-, with annual data in the gap
om2datastores = {om2exptname: intake.cat.access_nri[om2exptname]}

if IAF: om2exptname = '025deg_jra55_iaf_omip2_cycle1' else: om2exptname = '025deg_jra55_ryf9091_gadi' # monthly temp only for 1900-1903, and 2300-, with annual data in the gap om2datastores = {om2exptname: intake.cat.access_nri[om2exptname]}

In [13]:

om2varnames = [ 'pot_temp', 'salt' ]  # potential temperature and practical salinity

om2vars = {
    vname: {expt: ds.search(variable=vname, frequency="1mon").to_dask(
                xarray_open_kwargs = dict(
                    chunks={"xt_ocean": -1, "yt_ocean": "auto", "st_ocean" : "auto"},
                    decode_timedelta=True
                )
            )[vname]
            for expt, ds in om2datastores.items()
          }
    for vname in om2varnames
    }

if IAF:
    for vname, d in om2vars.items():
        for expt in d:
            try:
                d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True)
            except KeyError:
                pass

for expt, da in om2vars['pot_temp'].items():
    try:
        om2vars['pot_temp'][expt] = om2vars['pot_temp'][expt] - 273.15  # convert to degC
        om2vars['pot_temp'][expt].attrs['units'] = 'degC'
    except KeyError:
        pass

om2varnames = [ 'pot_temp', 'salt' ] # potential temperature and practical salinity om2vars = { vname: {expt: ds.search(variable=vname, frequency="1mon").to_dask( xarray_open_kwargs = dict( chunks={"xt_ocean": -1, "yt_ocean": "auto", "st_ocean" : "auto"}, decode_timedelta=True ) )[vname] for expt, ds in om2datastores.items() } for vname in om2varnames } if IAF: for vname, d in om2vars.items(): for expt in d: try: d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True) except KeyError: pass for expt, da in om2vars['pot_temp'].items(): try: om2vars['pot_temp'][expt] = om2vars['pot_temp'][expt] - 273.15 # convert to degC om2vars['pot_temp'][expt].attrs['units'] = 'degC' except KeyError: pass

Set time range¶

In [14]:

lasttime = []
for vname, d in om3vars.items():
    for expt, da in d.items():
        lasttime.append(da.time.values[-1])
for vname, d in om2vars.items():
    for expt, da in d.items():
        lasttime.append(da.time.values[-1])

if IAF:
    # datestop = cftime.DatetimeProlepticGregorian(2018, 12, 31, 12, 0, 0, 0, has_year_zero=True)  # to match OM2 IAF run
    datestop = min(lasttime) # earliest final date in data
    datelist = list(cftime.to_tuple(datestop))
    datelist[0] -= 25  # last 25 years
    datestart = cftime.datetime(*datelist, calendar=datestop.calendar)
    # bracket to start/end of year
    datestart = cftime.datetime(cftime.to_tuple(datestart)[0],  1, 1, calendar=datestop.calendar)
    datestop  = cftime.datetime(cftime.to_tuple(datestop)[0], 12, 31, 23, 59, 59, calendar=datestop.calendar)
    timerange = slice(datestart, datestop)
else:
    timerange = slice(cftime.DatetimeNoLeap(1932, 1, 1, 0, 0, 0, 0),
                      cftime.DatetimeNoLeap(1952, 1, 1, 0, 0, 0, 0))

# timerange = slice(None, None)

timerange

lasttime = [] for vname, d in om3vars.items(): for expt, da in d.items(): lasttime.append(da.time.values[-1]) for vname, d in om2vars.items(): for expt, da in d.items(): lasttime.append(da.time.values[-1]) if IAF: # datestop = cftime.DatetimeProlepticGregorian(2018, 12, 31, 12, 0, 0, 0, has_year_zero=True) # to match OM2 IAF run datestop = min(lasttime) # earliest final date in data datelist = list(cftime.to_tuple(datestop)) datelist[0] -= 25 # last 25 years datestart = cftime.datetime(*datelist, calendar=datestop.calendar) # bracket to start/end of year datestart = cftime.datetime(cftime.to_tuple(datestart)[0], 1, 1, calendar=datestop.calendar) datestop = cftime.datetime(cftime.to_tuple(datestop)[0], 12, 31, 23, 59, 59, calendar=datestop.calendar) timerange = slice(datestart, datestop) else: timerange = slice(cftime.DatetimeNoLeap(1932, 1, 1, 0, 0, 0, 0), cftime.DatetimeNoLeap(1952, 1, 1, 0, 0, 0, 0)) # timerange = slice(None, None) timerange

Out[14]:

slice(cftime.datetime(1993, 1, 1, 0, 0, 0, 0, calendar='proleptic_gregorian', has_year_zero=True), cftime.datetime(2018, 12, 31, 23, 59, 59, 0, calendar='proleptic_gregorian', has_year_zero=True), None)

Calculate zonal mean for each region¶

In [15]:

%%time
om3vars_mean = {
    region:
    {
        vname: {
                expt: da.sel(yh=slice(b['miny'], b['maxy']))\
                        .sel(xh=slice(b['minx'], b['maxx']))\
                        .sel(time=timerange)\
                        .mean('xh').mean('time').load()
                for expt, da in tqdm(vdatadict.items(), desc=f'      {vname} experiments', )
              }
        for vname, vdatadict in tqdm(om3vars.items(), desc=f'   {region} variables')
    }
    for region, b in tqdm(regions.items(), desc='regions')
}

%%time om3vars_mean = { region: { vname: { expt: da.sel(yh=slice(b['miny'], b['maxy']))\ .sel(xh=slice(b['minx'], b['maxx']))\ .sel(time=timerange)\ .mean('xh').mean('time').load() for expt, da in tqdm(vdatadict.items(), desc=f' {vname} experiments', ) } for vname, vdatadict in tqdm(om3vars.items(), desc=f' {region} variables') } for region, b in tqdm(regions.items(), desc='regions') }

CPU times: user 45.1 s, sys: 15.4 s, total: 1min
Wall time: 2min 54s

In [16]:

%%time
om2vars_mean = {
    region:
    {
        vname: {
                expt: da.sel(yt_ocean=slice(b['miny'], b['maxy']))\
                        .sel(xt_ocean=slice(b['minx'], b['maxx']))\
                        .sel(time=timerange)\
                        .mean('xt_ocean').mean('time').load()
                for expt, da in tqdm(vdatadict.items(), desc=f'      {vname} experiments', )
              }
        for vname, vdatadict in tqdm(om2vars.items(), desc=f'   {region} variables')
    }
    for region, b in tqdm(regions.items(), desc='regions')
}

%%time om2vars_mean = { region: { vname: { expt: da.sel(yt_ocean=slice(b['miny'], b['maxy']))\ .sel(xt_ocean=slice(b['minx'], b['maxx']))\ .sel(time=timerange)\ .mean('xt_ocean').mean('time').load() for expt, da in tqdm(vdatadict.items(), desc=f' {vname} experiments', ) } for vname, vdatadict in tqdm(om2vars.items(), desc=f' {region} variables') } for region, b in tqdm(regions.items(), desc='regions') }

CPU times: user 32.8 s, sys: 11.8 s, total: 44.6 s
Wall time: 2min 15s

Load data from WOA¶

NB: we're taking a mean of the monthly means. This reduces seasonal bias relative to the annual mean WOA (I think!), but will also involve more interpolated values.

The alternative would be to use the WOA annual mean /g/data/ik11/observations/woa23/woa23_decav_[ts]00_04.nc but this needs to be converted to the appropriate temperature and salinity, and regridded to model grid.

In [17]:

# fragile! obs_file_pattern depends on model run!
# obs_file_pattern = "/g/data/ik11/inputs/access-om3/woa23/025/2025.10.24/woa23_ts_*"
obs_file_pattern = "/g/data/ik11/inputs/access-om3/woa23/025/2026.06.11/woa23_ts_*"  # for MC_25km_jra_iaf+wombatlite-test4-d28e0359

om3mask = np.isfinite(om3vars['thetao'][om3exptname].isel(time=0)).load()

ds_obs = xr.open_mfdataset(
        obs_file_pattern,
        chunks={"time": -1})\
    .cf.mean("time")\
    .rename({'lon': 'xh', 'lat': 'yh', 'depth': 'z_l'})\
    .where(om3mask)  # mask out interpolated values at dry points in model

woa23vars = { k: ds_obs[k] for k in ["ptemp", "salt"] } # for OM3: potential temperature and practical salinity

# calculate zonal mean for each region
woa23vars_mean = {
    region:
    {
        vname: da.cf.sel( latitude=slice(b['miny'], b['maxy']))\
                 .cf.sel(longitude=slice(b['minx'], b['maxx']))\
                 .cf.mean('longitude')\
                 .load()
        for vname, da in tqdm(woa23vars.items(), desc=f'   {region} variables')
    }
    for region, b in tqdm(regions.items(), desc='regions')
}

# fragile! obs_file_pattern depends on model run! # obs_file_pattern = "/g/data/ik11/inputs/access-om3/woa23/025/2025.10.24/woa23_ts_*" obs_file_pattern = "/g/data/ik11/inputs/access-om3/woa23/025/2026.06.11/woa23_ts_*" # for MC_25km_jra_iaf+wombatlite-test4-d28e0359 om3mask = np.isfinite(om3vars['thetao'][om3exptname].isel(time=0)).load() ds_obs = xr.open_mfdataset( obs_file_pattern, chunks={"time": -1})\ .cf.mean("time")\ .rename({'lon': 'xh', 'lat': 'yh', 'depth': 'z_l'})\ .where(om3mask) # mask out interpolated values at dry points in model woa23vars = { k: ds_obs[k] for k in ["ptemp", "salt"] } # for OM3: potential temperature and practical salinity # calculate zonal mean for each region woa23vars_mean = { region: { vname: da.cf.sel( latitude=slice(b['miny'], b['maxy']))\ .cf.sel(longitude=slice(b['minx'], b['maxx']))\ .cf.mean('longitude')\ .load() for vname, da in tqdm(woa23vars.items(), desc=f' {region} variables') } for region, b in tqdm(regions.items(), desc='regions') }

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:38: UserWarning: An error occurred while loading the ERDDAP data fetcher, it will not be available !
<class 'ImportError'>
cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py)
  warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:50: UserWarning: An error occurred while loading the ArgoVis data fetcher, it will not be available !
<class 'ImportError'>
cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py)
  warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:66: UserWarning: An error occurred while loading the GDAC data fetcher, it will not be available !
<class 'ImportError'>
cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py)
  warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:91: UserWarning: An error occurred while loading the ERDDAP index fetcher, it will not be available !
<class 'ImportError'>
cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py)
  warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:107: UserWarning: An error occurred while loading the GDAC index fetcher, it will not be available !
<class 'ImportError'>
cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py)
  warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'argo' loading failed:
cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py)
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)

WARNING:pint.util:Redefining 'degrees_north' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degrees_N' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degreesN' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degree_north' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degree_N' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degreeN' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degrees_east' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degrees_E' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degreesE' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degree_east' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degree_E' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'degreeE' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

WARNING:pint.util:Redefining 'gpm' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/distributed/client.py:3387: UserWarning: Sending large graph of size 93.97 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/distributed/client.py:3387: UserWarning: Sending large graph of size 93.97 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

In [18]:

# fragile! obs_file_pattern depends on model run!
obs_file_pattern = "/g/data/ik11/observations/woa13/025_KDS50/woa13_ts_??_mom025.nc"  # for 025deg_jra55_iaf_omip2_cycle1

om2mask = np.isfinite(om2vars['pot_temp'][om2exptname].isel(time=0)).load()

ds_obs = xr.open_mfdataset(
        obs_file_pattern,
        chunks={"time": -1},
        decode_times=xr.coders.CFDatetimeCoder(use_cftime=True))\
    .cf.mean("time")\
    .rename({'GRID_X_T': 'xt_ocean', 'GRID_Y_T': 'yt_ocean', 'ZT': 'st_ocean'})

# fix to handle slightly different st_ocean precision
ds_obs.coords['st_ocean'] = om2mask.coords['st_ocean'].values
       
ds_obs = ds_obs.where(om2mask)  # mask out interpolated values at dry points in model

woa13vars = { k: ds_obs[k] for k in ["temp", "salt"] } # for OM2: potential temperature and practical salinity

# calculate zonal mean for each region
woa13vars_mean = {
    region:
    {
        vname: da.cf.sel( latitude=slice(b['miny'], b['maxy']))\
                 .cf.sel(longitude=slice(b['minx'], b['maxx']))\
                 .cf.mean('longitude')\
                 .load()
        for vname, da in tqdm(woa13vars.items(), desc=f'   {region} variables')
    }
    for region, b in tqdm(regions.items(), desc='regions')
}

# fragile! obs_file_pattern depends on model run! obs_file_pattern = "/g/data/ik11/observations/woa13/025_KDS50/woa13_ts_??_mom025.nc" # for 025deg_jra55_iaf_omip2_cycle1 om2mask = np.isfinite(om2vars['pot_temp'][om2exptname].isel(time=0)).load() ds_obs = xr.open_mfdataset( obs_file_pattern, chunks={"time": -1}, decode_times=xr.coders.CFDatetimeCoder(use_cftime=True))\ .cf.mean("time")\ .rename({'GRID_X_T': 'xt_ocean', 'GRID_Y_T': 'yt_ocean', 'ZT': 'st_ocean'}) # fix to handle slightly different st_ocean precision ds_obs.coords['st_ocean'] = om2mask.coords['st_ocean'].values ds_obs = ds_obs.where(om2mask) # mask out interpolated values at dry points in model woa13vars = { k: ds_obs[k] for k in ["temp", "salt"] } # for OM2: potential temperature and practical salinity # calculate zonal mean for each region woa13vars_mean = { region: { vname: da.cf.sel( latitude=slice(b['miny'], b['maxy']))\ .cf.sel(longitude=slice(b['minx'], b['maxx']))\ .cf.mean('longitude')\ .load() for vname, da in tqdm(woa13vars.items(), desc=f' {region} variables') } for region, b in tqdm(regions.items(), desc='regions') }

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/distributed/client.py:3387: UserWarning: Sending large graph of size 58.09 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/distributed/client.py:3387: UserWarning: Sending large graph of size 58.09 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

alternative approach for WOA data - not used¶

Would need to interpolate WOA data onto model grid

# use annual decadal average product (not initial condition, which is January/JFM) WOA_t_insitu = xr.open_dataset('/g/data/ik11/observations/woa23/woa23_decav_t00_04.nc', decode_times=False)['t_an'].isel(time=0, drop=True) WOA_s_psu = xr.open_dataset('/g/data/ik11/observations/woa23/woa23_decav_s00_04.nc', decode_times=False)['s_an'].isel(time=0, drop=True)# convert WOA in situ temperature to conservative temperature, # based on the conversion code for the model initial conditions # https://github.com/ACCESS-NRI/initial_conditions_access-om3/blob/8d0b700/setup_WOA_initial_conditions.py#L153-L175 lat = WOA_t_insitu.lat[:] lon = WOA_t_insitu.lon[:] longitude, latitude = np.meshgrid(lon, lat) pressure = gsw.p_from_z(-WOA_t_insitu.depth, lat) pressure_tile = np.tile(pressure, (len(lon), 1, 1)).swapaxes(0, 2).swapaxes(0, 1) WOA_s_absolute = gsw.SA_from_SP(WOA_s_psu, pressure_tile, longitude, latitude) WOA_s_absolute.attrs.pop('standard_name', None) WOA_s_absolute.attrs['long_name'] = "Absolute Salinity" WOA_s_absolute.attrs['units'] = "g kg-1" WOA_s_absolute.name = "WOA_s_absolute" WOA_t_conservative = gsw.CT_from_t(WOA_s_absolute, WOA_t_insitu, pressure_tile) WOA_t_conservative.attrs.pop('standard_name', None) WOA_t_conservative.attrs['long_name'] = "Conservative Temperature" WOA_t_conservative.attrs['units'] = "degC" WOA_t_conservative.name = "WOA_t_conservative" WOA_t_potential = gsw.pt0_from_t(WOA_s_absolute, WOA_t_insitu, pressure_tile) WOA_t_potential.attrs.pop('standard_name', None) WOA_t_potential.attrs['long_name'] = "Potential Temperature referenced to 0 dbar" WOA_t_potential.attrs['units'] = "degC" WOA_t_potential.name = "WOA_t_potential" del pressure del pressure_tile# Calculate WOA zonal means WOA_t_conservative_xmean = WOA_t_conservative.mean('lon') WOA_t_potential_xmean = WOA_t_potential.mean('lon')

Regridders¶

See https://github.com/aekiss/ice_analysis/blob/main/ice_maps.ipynb

import xesmfdef WOA_regridder(ds_in, ds_out, method='bilinear', **kwargs): '''Return a function that takes one dataarray argument on the same grid as ds_in and returns that dataarray regridded onto the grid from ds_out. ds_in -- WOA dataset with coords "lon" and "depth" ds_out -- ''' weightfn = '_'.join(['regrid', 'weights', os.path.splitext(os.path.basename(grid_in))[0], 'to', os.path.splitext(os.path.basename(grid_out))[0], method])+'.nc' ds_in = xr.open_dataset(grid_in)\ .rename({'xgrid': 'x', 'ygrid': 'y', 'longitude': 'lon', 'latitude': 'lat'}) try: # for /g/data/ik11/grids/ocean_grid_10.nc and /g/data/ik11/grids/ocean_grid_025.nc ds_out = xr.open_dataset(grid_out).rename({'xt_ocean': 'x', 'yt_ocean': 'y', 'geolon_t': 'lon', 'geolat_t': 'lat'}) except ValueError: # for /g/data/ik11/grids/ocean_grid_01.nc - see https://github.com/COSIMA/access-om2/issues/236 ds_out = xr.open_dataset(grid_out).rename({'grid_x_T': 'x', 'grid_y_T': 'y', 'geolon_t': 'lon', 'geolat_t': 'lat'}) ds_out_fix = xr.open_dataset('/g/data/ik11/outputs/access-om2-01/01deg_jra55v140_iaf_cycle3/output488/ocean/ocean-3d-temp-1-daily-mean-ym_1958_01.nc')\ .rename({'xt_ocean': 'x', 'yt_ocean': 'y'}) ds_out = ds_out.assign_coords({'x': ds_out_fix['x'], 'y': ds_out_fix['y']}) # make regridder rg = xesmf.Regridder(ds_in, ds_out, method=method, reuse_weights=os.path.exists(weightfn), filename=weightfn, **kwargs) def outf(da): # define regridding function to output outds = rg(da) outds.coords['x'] = ds_out['x'] outds.coords['y'] = ds_out['y'] return outds.rename({'x': 'xt_ocean', 'y': 'yt_ocean', 'lon': 'longitude', 'lat': 'latitude'}) return outf

Plot zonal means for each region¶

In [19]:

for region, d in om3vars_mean.items():
    om3variable = "thetao"  # potential temperature
    om2variable = "pot_temp"  # potential temperature
    woa23variable = "ptemp"  # potential temperature
    woa13variable = "temp"  # potential temperature

    npanels = 2 * (len(d[om3variable]) + len(om2vars_mean[region][om2variable]) + 1)
    nrows = (npanels+1)//2
    fig, axs = plt.subplots(nrows, 2,
                            figsize=(9, 3*nrows),
                            sharex=True, sharey=True)
    fig.suptitle(f'{region} zonal mean, {timerange.start.strftime("%Y-%m-%d")} - {timerange.stop.strftime("%Y-%m-%d")}', y=0.98)

    # Column 1: Temperature
    cmap = cmocean.cm.balance

    vlim = -1
    for da in d[om3variable].values():
        vlim = max(vlim, np.abs(da - woa23vars_mean[region][woa23variable]).max().data)
    for da in om2vars_mean[region][om2variable].values():
        vlim = max(vlim, np.abs(da - woa13vars_mean[region][woa13variable]).max().data)

    vlim = 1.8 # to match Kiss et al 2020 fig 12

    # ACCESS-OM3
    for i, (expt, da) in enumerate(d[om3variable].items()):
        # print(da)
        plt.subplot(nrows, 2, 2*i+1).set_facecolor('gray')  # gray land points
        p = (da - woa23vars_mean[region][woa23variable]).plot.contourf(
            levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False,
        )
        (da - woa23vars_mean[region][woa23variable]).plot.contour(
            levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"]
        )
        plt.gca().set_ylim([0, None])
        plt.gca().invert_yaxis()
        plt.title(f"({chr(ord('a')+2*i)}) ACCESS-OM3 {expt}\n{da.attrs['long_name']} ({om3variable}) minus WOA23", fontsize=9)
        plt.xlabel('')
        plt.ylabel('Depth [m]')

    # ACCESS-OM2
    for i, (expt, da) in enumerate(om2vars_mean[region][om2variable].items(), start=i+1):
        plt.subplot(nrows, 2, 2*i+1).set_facecolor('gray')  # gray land points
        p1 = (da - woa13vars_mean[region][woa13variable]).plot.contourf(
            levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False,
        )
        (da - woa13vars_mean[region][woa13variable]).plot.contour(
            levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"]
        )
        plt.gca().set_ylim([0, None])
        plt.gca().invert_yaxis()
        plt.title(f"({chr(ord('a')+2*i)}) ACCESS-OM2 {expt}\n{da.attrs['long_name']} ({om2variable}) minus WOA13", fontsize=9)
        plt.xlabel('')
        plt.ylabel('Depth [m]')

    # WOA23
    i += 1
    cmap = cmocean.cm.thermal
    da = woa23vars_mean[region][woa23variable]
    plt.subplot(nrows, 2, 2*i+1).set_facecolor('gray')  # gray land points
    p4 = da.plot.contourf(
        levels=101, cmap=cmap, add_colorbar=False,
        vmin=0, vmax=25  # to match Kiss et al 2020 fig 12
    )
    da.plot.contour(
        levels=25, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"]
    )
    plt.gca().set_ylim([0, None])
    plt.gca().invert_yaxis()
    plt.title(f"({chr(ord('a')+2*i)}) WOA23 {da.attrs['standard_name'].replace('_', ' ')}", fontsize=9)
    plt.xlabel('Latitude (°N)')
    plt.ylabel('Depth [m]')

    # Column 2: Salt
    om3variable = "so"
    om2variable = "salt"
    woa23variable = "salt"
    woa13variable = "salt"
    cmap = cmocean.cm.delta

    vlim = -1
    for da in d[om3variable].values():
        vlim = max(vlim, np.abs(da - woa23vars_mean[region][woa23variable]).max().data)
    for da in om2vars_mean[region][om2variable].values():
        vlim = max(vlim, np.abs(da - woa13vars_mean[region][woa13variable]).max().data)

    vlim = 0.8 # to match Kiss et al 2020 fig 12

    # ACCESS-OM3
    for i, (expt, da) in enumerate(d[om3variable].items()):
        # print(da)
        plt.subplot(nrows, 2, 2*i+2).set_facecolor('gray')  # gray land points
        (da - woa23vars_mean[region][woa23variable]).plot.contourf(
            levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False,
        )
        (da - woa23vars_mean[region][woa23variable]).plot.contour(
            levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"]
        )
        plt.gca().set_ylim([0, None])
        plt.gca().invert_yaxis()
        plt.title(f"({chr(ord('a')+2*i+1)}) ACCESS-OM3 {expt}\n{da.attrs['long_name']} ({om3variable}) minus WOA23", fontsize=9)
        plt.xlabel('')

    # ACCESS-OM2
    for i, (expt, da) in enumerate(om2vars_mean[region][om2variable].items(), start=i+1):
        plt.subplot(nrows, 2, 2*i+2).set_facecolor('gray')  # gray land points
        s1 = (da - woa13vars_mean[region][woa13variable]).plot.contourf(
            levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False,
        )
        (da - woa13vars_mean[region][woa13variable]).plot.contour(
            levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"]
        )
        plt.gca().set_ylim([0, None])
        plt.gca().invert_yaxis()
        plt.title(f"({chr(ord('a')+2*i+1)}) ACCESS-OM2 {expt}\n{da.attrs['long_name']} ({om2variable}) minus WOA13", fontsize=9)
        plt.xlabel('')

    # WOA23
    cmap = cmocean.cm.haline
    i += 1
    da = woa23vars_mean[region][woa23variable]
    plt.subplot(nrows, 2, 2*i+2).set_facecolor('gray')  # gray land points
    s4 = da.plot.contourf(
        levels=101, cmap=cmap, add_colorbar=False,
        vmin=33, vmax=35.8  # to match Kiss et al 2020 fig 12
    )
    da.plot.contour(
        levels=25, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"]
    )
    plt.gca().set_ylim([0, None])
    plt.gca().invert_yaxis()
    plt.title(f"({chr(ord('a')+2*i+1)}) WOA23 {da.attrs['standard_name'].replace('_', ' ')}", fontsize=9)
    plt.xlabel('Latitude (°N)')

    for ax in fig.get_axes():
        ax.label_outer()
    plt.tight_layout()

# from https://github.com/COSIMA/ACCESS-OM2-1-025-010deg-report/blob/master/figures/temperature_biases/zonal_temp_salt_bias.ipynb
    plt.subplots_adjust(bottom=0.14, hspace=0.24)
    
    ax1 = plt.axes([0.1,0.06,0.19,0.01])
    cb1 = plt.colorbar(p1, cax=ax1, orientation='horizontal', ticks=[-1.5, -0.75, 0, 0.75, 1.5])
    cb1.ax.xaxis.minorticks_off()
    cb1.ax.tick_params(labelsize=8)
    cb1.set_label('Temperature bias (°C)')
    
    ax2 = plt.axes([0.32,0.06,0.19,0.01])
    cb2 = plt.colorbar(p4, cax=ax2, orientation='horizontal', ticks=[0, 3, 6, 9, 12, 15, 18, 21, 24])
    cb2.ax.xaxis.minorticks_off()
    cb2.ax.tick_params(labelsize=8)
    cb2.set_label('Temperature (°C)')
    
    ax3 = plt.axes([0.57,0.06,0.19,0.01])
    cb3 = plt.colorbar(s1, cax=ax3, orientation='horizontal', ticks=[-0.8, -0.4, 0, 0.4, 0.8])
    cb3.ax.xaxis.minorticks_off()
    cb3.ax.tick_params(labelsize=8)
    cb3.set_label('Salinity bias (psu)')
    
    ax4 = plt.axes([0.8,0.06,0.19,0.01])
    cb4 = plt.colorbar(s4, cax=ax4, orientation='horizontal', ticks=[33, 33.5, 34, 34.5, 35, 35.5])
    cb4.ax.xaxis.minorticks_off()
    cb4.ax.tick_params(labelsize=8)
    cb4.set_label('Salinity (psu)')
    mkmd.savefig(fig,f"{region} zonal mean temperature and salinity transects", f"{region} zonal mean temperature and salinity bias of ACCESS-OM3 and ACCESS-OM2 relative to WOA. GitHub issues: [#11](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/11) Zonally averaged T/S bias, [#13](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/13) WOCE/GO-SHIP meridional transects", dpi=150)

for region, d in om3vars_mean.items(): om3variable = "thetao" # potential temperature om2variable = "pot_temp" # potential temperature woa23variable = "ptemp" # potential temperature woa13variable = "temp" # potential temperature npanels = 2 * (len(d[om3variable]) + len(om2vars_mean[region][om2variable]) + 1) nrows = (npanels+1)//2 fig, axs = plt.subplots(nrows, 2, figsize=(9, 3*nrows), sharex=True, sharey=True) fig.suptitle(f'{region} zonal mean, {timerange.start.strftime("%Y-%m-%d")} - {timerange.stop.strftime("%Y-%m-%d")}', y=0.98) # Column 1: Temperature cmap = cmocean.cm.balance vlim = -1 for da in d[om3variable].values(): vlim = max(vlim, np.abs(da - woa23vars_mean[region][woa23variable]).max().data) for da in om2vars_mean[region][om2variable].values(): vlim = max(vlim, np.abs(da - woa13vars_mean[region][woa13variable]).max().data) vlim = 1.8 # to match Kiss et al 2020 fig 12 # ACCESS-OM3 for i, (expt, da) in enumerate(d[om3variable].items()): # print(da) plt.subplot(nrows, 2, 2*i+1).set_facecolor('gray') # gray land points p = (da - woa23vars_mean[region][woa23variable]).plot.contourf( levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False, ) (da - woa23vars_mean[region][woa23variable]).plot.contour( levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"({chr(ord('a')+2*i)}) ACCESS-OM3 {expt}\n{da.attrs['long_name']} ({om3variable}) minus WOA23", fontsize=9) plt.xlabel('') plt.ylabel('Depth [m]') # ACCESS-OM2 for i, (expt, da) in enumerate(om2vars_mean[region][om2variable].items(), start=i+1): plt.subplot(nrows, 2, 2*i+1).set_facecolor('gray') # gray land points p1 = (da - woa13vars_mean[region][woa13variable]).plot.contourf( levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False, ) (da - woa13vars_mean[region][woa13variable]).plot.contour( levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"({chr(ord('a')+2*i)}) ACCESS-OM2 {expt}\n{da.attrs['long_name']} ({om2variable}) minus WOA13", fontsize=9) plt.xlabel('') plt.ylabel('Depth [m]') # WOA23 i += 1 cmap = cmocean.cm.thermal da = woa23vars_mean[region][woa23variable] plt.subplot(nrows, 2, 2*i+1).set_facecolor('gray') # gray land points p4 = da.plot.contourf( levels=101, cmap=cmap, add_colorbar=False, vmin=0, vmax=25 # to match Kiss et al 2020 fig 12 ) da.plot.contour( levels=25, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"({chr(ord('a')+2*i)}) WOA23 {da.attrs['standard_name'].replace('_', ' ')}", fontsize=9) plt.xlabel('Latitude (°N)') plt.ylabel('Depth [m]') # Column 2: Salt om3variable = "so" om2variable = "salt" woa23variable = "salt" woa13variable = "salt" cmap = cmocean.cm.delta vlim = -1 for da in d[om3variable].values(): vlim = max(vlim, np.abs(da - woa23vars_mean[region][woa23variable]).max().data) for da in om2vars_mean[region][om2variable].values(): vlim = max(vlim, np.abs(da - woa13vars_mean[region][woa13variable]).max().data) vlim = 0.8 # to match Kiss et al 2020 fig 12 # ACCESS-OM3 for i, (expt, da) in enumerate(d[om3variable].items()): # print(da) plt.subplot(nrows, 2, 2*i+2).set_facecolor('gray') # gray land points (da - woa23vars_mean[region][woa23variable]).plot.contourf( levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False, ) (da - woa23vars_mean[region][woa23variable]).plot.contour( levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"({chr(ord('a')+2*i+1)}) ACCESS-OM3 {expt}\n{da.attrs['long_name']} ({om3variable}) minus WOA23", fontsize=9) plt.xlabel('') # ACCESS-OM2 for i, (expt, da) in enumerate(om2vars_mean[region][om2variable].items(), start=i+1): plt.subplot(nrows, 2, 2*i+2).set_facecolor('gray') # gray land points s1 = (da - woa13vars_mean[region][woa13variable]).plot.contourf( levels=101, vmin=-vlim, vmax=vlim, cmap=cmap, add_colorbar=False, ) (da - woa13vars_mean[region][woa13variable]).plot.contour( levels=25, vmin=-vlim, vmax=vlim, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"({chr(ord('a')+2*i+1)}) ACCESS-OM2 {expt}\n{da.attrs['long_name']} ({om2variable}) minus WOA13", fontsize=9) plt.xlabel('') # WOA23 cmap = cmocean.cm.haline i += 1 da = woa23vars_mean[region][woa23variable] plt.subplot(nrows, 2, 2*i+2).set_facecolor('gray') # gray land points s4 = da.plot.contourf( levels=101, cmap=cmap, add_colorbar=False, vmin=33, vmax=35.8 # to match Kiss et al 2020 fig 12 ) da.plot.contour( levels=25, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"({chr(ord('a')+2*i+1)}) WOA23 {da.attrs['standard_name'].replace('_', ' ')}", fontsize=9) plt.xlabel('Latitude (°N)') for ax in fig.get_axes(): ax.label_outer() plt.tight_layout() # from https://github.com/COSIMA/ACCESS-OM2-1-025-010deg-report/blob/master/figures/temperature_biases/zonal_temp_salt_bias.ipynb plt.subplots_adjust(bottom=0.14, hspace=0.24) ax1 = plt.axes([0.1,0.06,0.19,0.01]) cb1 = plt.colorbar(p1, cax=ax1, orientation='horizontal', ticks=[-1.5, -0.75, 0, 0.75, 1.5]) cb1.ax.xaxis.minorticks_off() cb1.ax.tick_params(labelsize=8) cb1.set_label('Temperature bias (°C)') ax2 = plt.axes([0.32,0.06,0.19,0.01]) cb2 = plt.colorbar(p4, cax=ax2, orientation='horizontal', ticks=[0, 3, 6, 9, 12, 15, 18, 21, 24]) cb2.ax.xaxis.minorticks_off() cb2.ax.tick_params(labelsize=8) cb2.set_label('Temperature (°C)') ax3 = plt.axes([0.57,0.06,0.19,0.01]) cb3 = plt.colorbar(s1, cax=ax3, orientation='horizontal', ticks=[-0.8, -0.4, 0, 0.4, 0.8]) cb3.ax.xaxis.minorticks_off() cb3.ax.tick_params(labelsize=8) cb3.set_label('Salinity bias (psu)') ax4 = plt.axes([0.8,0.06,0.19,0.01]) cb4 = plt.colorbar(s4, cax=ax4, orientation='horizontal', ticks=[33, 33.5, 34, 34.5, 35, 35.5]) cb4.ax.xaxis.minorticks_off() cb4.ax.tick_params(labelsize=8) cb4.set_label('Salinity (psu)') mkmd.savefig(fig,f"{region} zonal mean temperature and salinity transects", f"{region} zonal mean temperature and salinity bias of ACCESS-OM3 and ACCESS-OM2 relative to WOA. GitHub issues: [#11](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/11) Zonally averaged T/S bias, [#13](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/13) WOCE/GO-SHIP meridional transects", dpi=150)

Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/mkmd/temp-salt-vs-depth-latitude_01.png
Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/mkmd/temp-salt-vs-depth-latitude.md

Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite_bs2-mpudig-backscat2-0578cc36/mkmd/temp-salt-vs-depth-latitude.md successfully.

def doplots(om3variable = "thetao", om2variable = "pot_temp", vmin = None, vmax = None, cmap = cmocean.cm.thermal): for region, d in om3vars_mean.items(): npanels = len(d[om3variable]) + len(om2vars_mean[region][om2variable]) nrows = (npanels+1)//2 fig, axs = plt.subplots(nrows, 2, figsize=(9, 3*nrows), sharex=True, sharey=True) fig.suptitle(f'{region} zonal mean, {timerange.start.strftime("%Y-%m-%d")} - {timerange.stop.strftime("%Y-%m-%d")}') if vmin == None: # calculate min over all panels (different for each region) vmin = 1e20 for i, (expt, da) in enumerate(d[om3variable].items(), start=1): vmin = min(vmin, da.min().data) for i2, (expt, da) in enumerate(om2vars_mean[region][om2variable].items(), start=i+1): vmin = min(vmin, da.min().data) if vmax == None: # calculate max over all panels (different for each region) vmax = -1e20 for i, (expt, da) in enumerate(d[om3variable].items(), start=1): vmax = max(vmax, da.max().data) for i2, (expt, da) in enumerate(om2vars_mean[region][om2variable].items(), start=i+1): vmax = max(vmax, da.max().data) # ACCESS-OM3 for i, (expt, da) in enumerate(d[om3variable].items(), start=1): # print(da) plt.subplot(nrows, 2, i).set_facecolor('gray') # gray land points p = da.plot.contourf( levels=101, vmin=vmin, vmax=vmax, cmap=cmap, add_colorbar=False, ) da.plot.contour( levels=25, vmin=vmin, vmax=vmax, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"ACCESS-OM3 {expt} {om3variable}", fontsize=8) plt.xlabel('') plt.ylabel('Depth [m]') label = f"ACCESS-OM3: {da.attrs['long_name']} ({om3variable}) [{da.attrs['units']}]" # ACCESS-OM2 for i2, (expt, da) in enumerate(om2vars_mean[region][om2variable].items(), start=i+1): plt.subplot(nrows, 2, i2).set_facecolor('gray') # gray land points da.plot.contourf( levels=101, vmin=vmin, vmax=vmax, cmap=cmap, add_colorbar=False, ) da.plot.contour( levels=25, vmin=vmin, vmax=vmax, add_colorbar=False, colors=["k"], linewidths=[0.2], linestyles=["-"] ) plt.gca().set_ylim([0, None]) plt.gca().invert_yaxis() plt.title(f"ACCESS-OM2 {expt} {om2variable}", fontsize=8) plt.xlabel('') plt.ylabel('Depth [m]') label += f"\nACCESS-OM2: {da.attrs['long_name']} ({om2variable}) [{da.attrs['units']}]" for ax in fig.get_axes(): ax.label_outer() plt.tight_layout() fig.colorbar(p, ax=axs, orientation='horizontal', fraction=.1, shrink=0.7, aspect=40, label=label) mkmd.savefig(fig,f"Zonal mean temperature and salinity transects: {region} ({om3variable})", "ACCESS-OM3 and ACCESS-OM2 zonal mean "+om3variable+" GitHub issues: [#11](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/11) Zonally averaged T/S bias, [#13](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/13) WOCE/GO-SHIP meridional transects", dpi=150) doplots(om3variable="thetao", om2variable="pot_temp")doplots(om3variable="so", om2variable="salt", cmap=cmocean.cm.haline)

In [20]:

client.close()

client.close()