SSS and SSS difference from WOA23¶

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

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-test3v2-00532b88/datastore.json"

# What physical field you want, in CF terms (stays the same across models)
variable_standard_name = "sea_surface_salinity"

# Fallback variable name to use if the catalog doesn’t expose standard_name
fallback_variable_names = ["sos", "salinity", "salt", "psalt", "asalt"]

# Reference / observational dataset (we'll use this later when we get to the WOA block)
obs_file_pattern = "/g/data/ik11/inputs/access-om3/woa23/025/2025.08.26/woa23_ts_*"
obs_var_name = "salt"

# Frequency depends on the datastore you’re using:
data_frequency = "1mon"

# 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-test3v2-00532b88/datastore.json" # What physical field you want, in CF terms (stays the same across models) variable_standard_name = "sea_surface_salinity" # Fallback variable name to use if the catalog doesn’t expose standard_name fallback_variable_names = ["sos", "salinity", "salt", "psalt", "asalt"] # Reference / observational dataset (we'll use this later when we get to the WOA block) obs_file_pattern = "/g/data/ik11/inputs/access-om3/woa23/025/2025.08.26/woa23_ts_*" obs_var_name = "salt" # Frequency depends on the datastore you’re using: data_frequency = "1mon" # 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-1.0-beta-5165c0f8/datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/"
nbname = "SSS.ipynb"

# Parameters esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-5165c0f8/datastore.json" papermill = True cwd = "/g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/" nbname = "SSS.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)
from mkfigs_configdoc import MkmdWriter
mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill)
from exptdata_access import get_experiment_info, guess_experiment_from_esm_file
from model_agnostic import get_lon_lat_from_catalog, select_variable
# Infer model information from the datastore path
expt_key, info = guess_experiment_from_esm_file(esm_file)
model_name = info["model"]             # OM3 or CM3

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) from mkfigs_configdoc import MkmdWriter mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill) from exptdata_access import get_experiment_info, guess_experiment_from_esm_file from model_agnostic import get_lon_lat_from_catalog, select_variable # Infer model information from the datastore path expt_key, info = guess_experiment_from_esm_file(esm_file) model_name = info["model"] # OM3 or CM3

In [4]:

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

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

In [5]:

client = Client(threads_per_worker=1)
client

client = Client(threads_per_worker=1) client

Out[5]:

Client

Client-653d4cdf-7389-11f1-adc8-0000018dfe80

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

Launch dashboard in JupyterLab

Cluster Info

LocalCluster

313b4014

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-50d62fd6-68ef-424c-89e2-8cead1127d81

Comm: tcp://127.0.0.1:45635	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:44173	Total threads: 1
Dashboard: http://127.0.0.1:34983/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:39999
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-ir8t7wx0

Worker: 1

Comm: tcp://127.0.0.1:36353	Total threads: 1
Dashboard: http://127.0.0.1:39533/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:33629
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-p_tcpmg1

Worker: 2

Comm: tcp://127.0.0.1:34077	Total threads: 1
Dashboard: http://127.0.0.1:33977/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:38143
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-wazqg_fn

Worker: 3

Comm: tcp://127.0.0.1:35349	Total threads: 1
Dashboard: http://127.0.0.1:35773/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:41319
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-_7ueqolp

Worker: 4

Comm: tcp://127.0.0.1:34277	Total threads: 1
Dashboard: http://127.0.0.1:39423/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:46815
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-achgyfq6

Worker: 5

Comm: tcp://127.0.0.1:37569	Total threads: 1
Dashboard: http://127.0.0.1:39201/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:39067
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-xlblu6wi

Worker: 6

Comm: tcp://127.0.0.1:42049	Total threads: 1
Dashboard: http://127.0.0.1:37213/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:41435
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-1dugpwze

Worker: 7

Comm: tcp://127.0.0.1:44399	Total threads: 1
Dashboard: http://127.0.0.1:42647/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:33307
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-91u71uc7

Worker: 8

Comm: tcp://127.0.0.1:38111	Total threads: 1
Dashboard: http://127.0.0.1:38965/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:45829
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-nrcxk885

Worker: 9

Comm: tcp://127.0.0.1:45157	Total threads: 1
Dashboard: http://127.0.0.1:42169/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:32933
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-tqn4u54h

Worker: 10

Comm: tcp://127.0.0.1:45113	Total threads: 1
Dashboard: http://127.0.0.1:34391/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:35601
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-d7i9u2p9

Worker: 11

Comm: tcp://127.0.0.1:44341	Total threads: 1
Dashboard: http://127.0.0.1:39309/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:42525
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-ldcsg3nm

Worker: 12

Comm: tcp://127.0.0.1:41465	Total threads: 1
Dashboard: http://127.0.0.1:35755/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:39995
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-l9pyzyw_

Worker: 13

Comm: tcp://127.0.0.1:35647	Total threads: 1
Dashboard: http://127.0.0.1:39291/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:41299
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-3m6evy8r

Worker: 14

Comm: tcp://127.0.0.1:33981	Total threads: 1
Dashboard: http://127.0.0.1:35613/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:46611
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-ru4v7njl

Worker: 15

Comm: tcp://127.0.0.1:39281	Total threads: 1
Dashboard: http://127.0.0.1:42091/status	Memory: 11.78 GiB
Nanny: tcp://127.0.0.1:35293
Local directory: /jobfs/172626275.gadi-pbs/dask-scratch-space/worker-2b7u4c1g

Define plot function¶

In [6]:

blue_marble = plt.imread('/g/data/ik11/grids/BlueMarble.tiff')
blue_marble_extent = (-180, 180, -90, 90)

blue_marble = plt.imread('/g/data/ik11/grids/BlueMarble.tiff') blue_marble_extent = (-180, 180, -90, 90)

In [7]:

def plot(dat, title=None, projection=None, add_blue_marble=True, **kwargs):
    """
    Generic 2D map plot helper for model / obs / bias fields.

    - Works for any variable (SST, SSH, salinity, etc.).
    - Assumes the last two dimensions are (y, x) horizontal dims.
    - Uses contourf with a cartopy projection.
    """

    import matplotlib.pyplot as plt
    import cartopy.crs as ccrs
    from textwrap import wrap

    # Make sure we're working with a DataArray
    da = dat
    if not isinstance(da, xr.DataArray):
        raise TypeError("plot() expects an xarray.DataArray")

    # Handle simple time dimension cases (e.g. time-mean or single time step)
    if da.ndim == 3 and "time" in da.dims and da.sizes["time"] == 1:
        da = da.isel(time=0, drop=True)

    if da.ndim != 2:
        raise ValueError(
            f"plot() expects a 2D field or a 3D field with singleton time; got dims {da.dims}"
        )

    # Infer horizontal dims as the last two dims
    ydim, xdim = da.dims[-2], da.dims[-1]

    # Title and colourbar label from attributes
    long_name = da.attrs.get("long_name", da.name or "Field")
    units = da.attrs.get("units", "")

    if title is None:
        title = long_name

    cbar_label = long_name if units == "" else f"{long_name} [{units}]"
    cbar_label = "\n".join(wrap(cbar_label, 45))

    # Projection: default to Robinson
    if projection is None:
        projection = ccrs.Robinson(central_longitude=-100)

    fig = plt.figure(figsize=(12, 6))
    ax = plt.axes(projection=projection)

    # Main filled contour plot
    da.plot.contourf(
        ax=ax,
        x=xdim,
        y=ydim,
        transform=ccrs.PlateCarree(),
        cbar_kwargs={
            "label": cbar_label,
            "fraction": 0.03,
            "aspect": 15,
            "shrink": 0.7,
        },
        **kwargs,
    )

    # Optional blue marble overlay if available
    if add_blue_marble and "blue_marble" in globals() and "blue_marble_extent" in globals():
        ax.imshow(
            blue_marble,
            extent=blue_marble_extent,
            transform=ccrs.PlateCarree(),
            origin="upper",
        )

    ax.coastlines()
    plt.title(title)
    plt.tight_layout()

def plot(dat, title=None, projection=None, add_blue_marble=True, **kwargs): """ Generic 2D map plot helper for model / obs / bias fields. - Works for any variable (SST, SSH, salinity, etc.). - Assumes the last two dimensions are (y, x) horizontal dims. - Uses contourf with a cartopy projection. """ import matplotlib.pyplot as plt import cartopy.crs as ccrs from textwrap import wrap # Make sure we're working with a DataArray da = dat if not isinstance(da, xr.DataArray): raise TypeError("plot() expects an xarray.DataArray") # Handle simple time dimension cases (e.g. time-mean or single time step) if da.ndim == 3 and "time" in da.dims and da.sizes["time"] == 1: da = da.isel(time=0, drop=True) if da.ndim != 2: raise ValueError( f"plot() expects a 2D field or a 3D field with singleton time; got dims {da.dims}" ) # Infer horizontal dims as the last two dims ydim, xdim = da.dims[-2], da.dims[-1] # Title and colourbar label from attributes long_name = da.attrs.get("long_name", da.name or "Field") units = da.attrs.get("units", "") if title is None: title = long_name cbar_label = long_name if units == "" else f"{long_name} [{units}]" cbar_label = "\n".join(wrap(cbar_label, 45)) # Projection: default to Robinson if projection is None: projection = ccrs.Robinson(central_longitude=-100) fig = plt.figure(figsize=(12, 6)) ax = plt.axes(projection=projection) # Main filled contour plot da.plot.contourf( ax=ax, x=xdim, y=ydim, transform=ccrs.PlateCarree(), cbar_kwargs={ "label": cbar_label, "fraction": 0.03, "aspect": 15, "shrink": 0.7, }, **kwargs, ) # Optional blue marble overlay if available if add_blue_marble and "blue_marble" in globals() and "blue_marble_extent" in globals(): ax.imshow( blue_marble, extent=blue_marble_extent, transform=ccrs.PlateCarree(), origin="upper", ) ax.coastlines() plt.title(title) plt.tight_layout()

Load and plot data from ESM datastore¶

In [8]:

exptname=os.path.basename(os.path.dirname(esm_file))
print("Experiment name:", exptname)

datastore = intake.open_esm_datastore(
    esm_file,
    columns_with_iterables=[
        "variable",
        "variable_long_name",
        "variable_standard_name",
        "variable_cell_methods",
        "variable_units"
    ]
)

exptname=os.path.basename(os.path.dirname(esm_file)) print("Experiment name:", exptname) datastore = intake.open_esm_datastore( esm_file, columns_with_iterables=[ "variable", "variable_long_name", "variable_standard_name", "variable_cell_methods", "variable_units" ] )

Experiment name: MC_25km_jra_iaf-1.0-beta-5165c0f8

In [9]:

da_model = select_variable(
    datastore,
    variable_standard_name,
    fallback_variable_names,
    data_frequency=data_frequency,
)
print("Selected variable:", da_model.name)
print("dims:", da_model.dims)

da_model = select_variable( datastore, variable_standard_name, fallback_variable_names, data_frequency=data_frequency, ) print("Selected variable:", da_model.name) print("dims:", da_model.dims)

CF-based lookup failed: KeyError("Receive multiple variables for key 'sea_surface_salinity': {'sos_max', 'sos_min'}. Expected only one. Please pass a list ['sea_surface_salinity'] instead to get all variables matching 'sea_surface_salinity'.")
No files matched frequency='1mon'; using any available frequency.

Selected variable via fallback name: sos
Selected variable: sos
dims: ('time', 'yh', 'xh')

In [10]:

# Load model grid and attach true lon/lat to the model field
static = xr.open_dataset(os.path.join(os.path.dirname(esm_file), "output000", "access-om3.mom6.static.nc"))

geolon = static["geolon"]
geolat = static["geolat"]

model_all = da_model.cf.assign_coords(
    {
        "longitude": geolon,
        "latitude": geolat,
    }
)
print("Attached geolon/geolat from static grid.")
print("model_all dims:", model_all.dims)

# Load model grid and attach true lon/lat to the model field static = xr.open_dataset(os.path.join(os.path.dirname(esm_file), "output000", "access-om3.mom6.static.nc")) geolon = static["geolon"] geolat = static["geolat"] model_all = da_model.cf.assign_coords( { "longitude": geolon, "latitude": geolat, } ) print("Attached geolon/geolat from static grid.") print("model_all dims:", model_all.dims)

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'>)

Attached geolon/geolat from static grid.
model_all dims: ('time', 'yh', 'xh')

In [11]:

averaging_mode = "last_n_years"   # or "full_period" / "fixed_period"
averaging_last_n_years = 10
averaging_start_date = None
averaging_end_date   = None

averaging_mode = "last_n_years" # or "full_period" / "fixed_period" averaging_last_n_years = 10 averaging_start_date = None averaging_end_date = None

In [12]:

# Put model time axis on a common calendar
model_all = model_all.convert_calendar("proleptic_gregorian", use_cftime=True)

# Inspect full time coverage
t0 = model_all.time.values[0]
t1 = model_all.time.values[-1]
print("Full model time range:", t0, "→", t1)

# Decide the averaging window based on configuration
if averaging_mode == "full_period":
    datestart = t0
    datestop = t1

elif averaging_mode == "last_n_years":
    datestop = t1
    # Subtract N years in a calendar-aware way
    datelist = list(cftime.to_tuple(datestop))
    datelist[0] -= averaging_last_n_years
    datestart = cftime.datetime(*datelist, calendar=datestop.calendar)

elif averaging_mode == "fixed_period":
    if averaging_start_date is None or averaging_end_date is None:
        raise ValueError(
            "averaging_mode='fixed_period' requires averaging_start_date and averaging_end_date"
        )
    # Use cftime to create start/end on the same calendar
    datestart = xr.cftime_range(
        start=averaging_start_date,
        periods=1,
        calendar="proleptic_gregorian",
    )[0]
    datestop = xr.cftime_range(
        start=averaging_end_date,
        periods=1,
        calendar="proleptic_gregorian",
    )[0]
else:
    raise ValueError(f"Unknown averaging_mode: {averaging_mode!r}")

timerange = slice(datestart, datestop)
print("Averaging window:", timerange)

# Restrict model to this configured window
model_window = model_all.cf.sel(time=timerange)
print("Windowed dims:", model_window.dims)

# Put model time axis on a common calendar model_all = model_all.convert_calendar("proleptic_gregorian", use_cftime=True) # Inspect full time coverage t0 = model_all.time.values[0] t1 = model_all.time.values[-1] print("Full model time range:", t0, "→", t1) # Decide the averaging window based on configuration if averaging_mode == "full_period": datestart = t0 datestop = t1 elif averaging_mode == "last_n_years": datestop = t1 # Subtract N years in a calendar-aware way datelist = list(cftime.to_tuple(datestop)) datelist[0] -= averaging_last_n_years datestart = cftime.datetime(*datelist, calendar=datestop.calendar) elif averaging_mode == "fixed_period": if averaging_start_date is None or averaging_end_date is None: raise ValueError( "averaging_mode='fixed_period' requires averaging_start_date and averaging_end_date" ) # Use cftime to create start/end on the same calendar datestart = xr.cftime_range( start=averaging_start_date, periods=1, calendar="proleptic_gregorian", )[0] datestop = xr.cftime_range( start=averaging_end_date, periods=1, calendar="proleptic_gregorian", )[0] else: raise ValueError(f"Unknown averaging_mode: {averaging_mode!r}") timerange = slice(datestart, datestop) print("Averaging window:", timerange) # Restrict model to this configured window model_window = model_all.cf.sel(time=timerange) print("Windowed dims:", model_window.dims)

Full model time range: 1958-01-01 12:00:00 → 2014-12-31 12:00:00
Averaging window: slice(cftime.datetime(2004, 12, 31, 12, 0, 0, 0, calendar='proleptic_gregorian', has_year_zero=True), cftime.DatetimeProlepticGregorian(2014, 12, 31, 12, 0, 0, 0, has_year_zero=True), None)
Windowed dims: ('time', 'yh', 'xh')

In [13]:

model = (
    model_window
    .cf.mean("time")
    .pint.quantify()
    .pint.to("psu")           # only works if original units are convertible
    .pint.dequantify()
)

model = ( model_window .cf.mean("time") .pint.quantify() .pint.to("psu") # only works if original units are convertible .pint.dequantify() )

In [14]:

long_name = model.attrs.get("long_name", "SSS")

plot(
    model,
    levels=25,
    vmin=28,
    vmax=40,
    extend="both",
    cmap=cm.cm.haline,
    title=(
        f"{long_name} "
        f"{datestart.strftime('%Y-%m-%d')} – {datestop.strftime('%Y-%m-%d')} "
        f"mean\nin {exptname}"
    ),
)
mkmd.savefig(plt.gcf(), "Sea Surface Salinity", "ACCESS-OM3 sea surface salinity. [GitHub issue: Global SSS bias](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/10)")

long_name = model.attrs.get("long_name", "SSS") plot( model, levels=25, vmin=28, vmax=40, extend="both", cmap=cm.cm.haline, title=( f"{long_name} " f"{datestart.strftime('%Y-%m-%d')} – {datestop.strftime('%Y-%m-%d')} " f"mean\nin {exptname}" ), ) mkmd.savefig(plt.gcf(), "Sea Surface Salinity", "ACCESS-OM3 sea surface salinity. [GitHub issue: Global SSS bias](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/10)")

Saved /g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/mkmd/SSS_01.png
Adding entry to per-notebook markdown: /g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/mkmd/SSS.md
Lines appended to /g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/mkmd/SSS.md successfully.

Load data from WOA23 (annual mean; Jan is ACCESS-OM3 initial condition)¶

In [15]:

# Load observational / reference SST (e.g. WOA23) in a model-agnostic way

# 1) Open the obs dataset using the pattern & var name from the config cell
ds_obs = xr.open_mfdataset(
    obs_file_pattern,
    chunks={"time": -1},
)

da_obs = ds_obs[obs_var_name]

# 2) Take surface level + time mean using CF-aware indexing
#    (this will work as long as depth/time have CF metadata)
try:
    da_obs_sfc = da_obs.cf.isel(Z=0)          # first depth level
except Exception:
    # Fallback: if depth is literally called 'depth'
    da_obs_sfc = da_obs.isel(depth=0)

obs = da_obs_sfc.cf.mean("time")

# 3) Make sure longitude / latitude coords are available in CF terms
if {"yh", "xh"}.issubset(obs.dims):
    obs_grid = obs
elif {"lat", "lon"}.issubset(obs.dims):
    obs_grid = obs.rename({"lat": "yh", "lon": "xh"})
else:
    raise ValueError("Obs grid dims not recognized; expected yh/xh or lat/lon")

obs = obs_grid.cf.assign_coords(
    {
        "longitude": geolon,
        "latitude": geolat,
    }
)

print("Obs dims:", obs.dims)
print("Obs coords:", list(obs.coords))

# Load observational / reference SST (e.g. WOA23) in a model-agnostic way # 1) Open the obs dataset using the pattern & var name from the config cell ds_obs = xr.open_mfdataset( obs_file_pattern, chunks={"time": -1}, ) da_obs = ds_obs[obs_var_name] # 2) Take surface level + time mean using CF-aware indexing # (this will work as long as depth/time have CF metadata) try: da_obs_sfc = da_obs.cf.isel(Z=0) # first depth level except Exception: # Fallback: if depth is literally called 'depth' da_obs_sfc = da_obs.isel(depth=0) obs = da_obs_sfc.cf.mean("time") # 3) Make sure longitude / latitude coords are available in CF terms if {"yh", "xh"}.issubset(obs.dims): obs_grid = obs elif {"lat", "lon"}.issubset(obs.dims): obs_grid = obs.rename({"lat": "yh", "lon": "xh"}) else: raise ValueError("Obs grid dims not recognized; expected yh/xh or lat/lon") obs = obs_grid.cf.assign_coords( { "longitude": geolon, "latitude": geolat, } ) print("Obs dims:", obs.dims) print("Obs coords:", list(obs.coords))

Obs dims: ('yh', 'xh')
Obs coords: ['depth', 'xh', 'yh']

Plot model minus WOA23¶

First plot full range, then a sequence at specified ranges.

In [16]:

# 1) Obs is already on the model grid with geolon/geolat coordinates
obs_on_model = obs

# 2) Compute bias directly (units already consistent)
bias = (model - obs_on_model).load()
bias.attrs = model.attrs

# 3) Plot the bias
long_name = model.attrs.get("long_name", "SSS")

plot(
    bias,
    levels=51,
    extend="both",
    cmap="seismic",
    title=(
        f"{long_name} bias "
        f"{datestart.strftime('%Y-%m-%d')} – {datestop.strftime('%Y-%m-%d')} "
        f"mean\n{exptname} minus WOA23"
    ),
)

# 1) Obs is already on the model grid with geolon/geolat coordinates obs_on_model = obs # 2) Compute bias directly (units already consistent) bias = (model - obs_on_model).load() bias.attrs = model.attrs # 3) Plot the bias long_name = model.attrs.get("long_name", "SSS") plot( bias, levels=51, extend="both", cmap="seismic", title=( f"{long_name} bias " f"{datestart.strftime('%Y-%m-%d')} – {datestop.strftime('%Y-%m-%d')} " f"mean\n{exptname} minus WOA23" ), )

In [17]:

plot(
    bias,
    levels=51,
    vmin=-5,
    vmax=5,
    extend="both",
    cmap="seismic",
    title=(
        f"{long_name} bias "
        f"{datestart.strftime('%Y-%m-%d')} – {datestop.strftime('%Y-%m-%d')} "
        f"mean\n{exptname} minus WOA23"
    ),
)
mkmd.savefig(plt.gcf(), "Sea Surface Salinity Bias", "ACCESS-OM3 sea surface salinity minus WOA2023. [GitHub issue: Global SSS bias](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/10)")

plot( bias, levels=51, vmin=-5, vmax=5, extend="both", cmap="seismic", title=( f"{long_name} bias " f"{datestart.strftime('%Y-%m-%d')} – {datestop.strftime('%Y-%m-%d')} " f"mean\n{exptname} minus WOA23" ), ) mkmd.savefig(plt.gcf(), "Sea Surface Salinity Bias", "ACCESS-OM3 sea surface salinity minus WOA2023. [GitHub issue: Global SSS bias](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/10)")

Saved /g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/mkmd/SSS_02.png
Adding entry to per-notebook markdown: /g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/mkmd/SSS.md
Lines appended to /g/data/tm70/cyb561/repos/access-om3-paper-1/notebooks/mkfigs_output_MC_25km_jra_iaf-1.0-beta-5165c0f8/mkmd/SSS.md successfully.

In [18]:

client.close()

client.close()