Meridional Overturning
mom6_tools.moc collection of functions for computing and plotting meridional overturning circulation.
The goal of this notebook is the following:
server as an example on to compute a meridional overturning streamfunction (global and Atalntic) from CESM/MOM output;
evaluate model experiments by comparing transports against observed estimates and other model results.
[1]:
%load_ext autoreload
%autoreload 2
[2]:
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")
import matplotlib
import numpy as np
import xarray as xr
# mom6_tools
from mom6_tools.MOM6grid import MOM6grid
from mom6_tools.moc import *
from mom6_tools.jobqueue import get_cluster
from mom6_tools.m6toolbox import genBasinMasks, add_global_attrs
from mom6_tools.m6toolbox import cime_xmlquery
import matplotlib.pyplot as plt
import yaml, os, intake
Basemap module not found. Some regional plots may not function properly
[3]:
# Read in the yaml file
diag_config_yml_path = "diag_config.yml"
diag_config_yml = yaml.load(open(diag_config_yml_path,'r'), Loader=yaml.Loader)
[4]:
caseroot = diag_config_yml['Case']['CASEROOT']
casename = cime_xmlquery(caseroot, 'CASE')
DOUT_S = cime_xmlquery(caseroot, 'DOUT_S')
if DOUT_S:
OUTDIR = cime_xmlquery(caseroot, 'DOUT_S_ROOT')+'/ocn/hist/'
else:
OUTDIR = cime_xmlquery(caseroot, 'RUNDIR')
print('Output directory is:', OUTDIR)
print('Casename is:', casename)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[4], line 2
1 caseroot = diag_config_yml['Case']['CASEROOT']
----> 2 casename = cime_xmlquery(caseroot, 'CASE')
3 DOUT_S = cime_xmlquery(caseroot, 'DOUT_S')
4 if DOUT_S:
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/mom6_tools/m6toolbox.py:47, in cime_xmlquery(caseroot, varname)
45 """run CIME's xmlquery for varname in the directory caseroot, return the value"""
46 try:
---> 47 value = subprocess.check_output(
48 ["./xmlquery", "-N", "--value", varname],
49 stderr=subprocess.STDOUT,
50 cwd=caseroot,
51 )
52 except subprocess.CalledProcessError:
53 value = subprocess.check_output(
54 ["./xmlquery", "--value", varname], stderr=subprocess.STDOUT, cwd=caseroot
55 )
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:421, in check_output(timeout, *popenargs, **kwargs)
418 empty = b''
419 kwargs['input'] = empty
--> 421 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
422 **kwargs).stdout
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:503, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
500 kwargs['stdout'] = PIPE
501 kwargs['stderr'] = PIPE
--> 503 with Popen(*popenargs, **kwargs) as process:
504 try:
505 stdout, stderr = process.communicate(input, timeout=timeout)
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:971, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize)
967 if self.text_mode:
968 self.stderr = io.TextIOWrapper(self.stderr,
969 encoding=encoding, errors=errors)
--> 971 self._execute_child(args, executable, preexec_fn, close_fds,
972 pass_fds, cwd, env,
973 startupinfo, creationflags, shell,
974 p2cread, p2cwrite,
975 c2pread, c2pwrite,
976 errread, errwrite,
977 restore_signals,
978 gid, gids, uid, umask,
979 start_new_session)
980 except:
981 # Cleanup if the child failed starting.
982 for f in filter(None, (self.stdin, self.stdout, self.stderr)):
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:1863, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session)
1861 if errno_num != 0:
1862 err_msg = os.strerror(errno_num)
-> 1863 raise child_exception_type(errno_num, err_msg, err_filename)
1864 raise child_exception_type(err_msg)
FileNotFoundError: [Errno 2] No such file or directory: '/glade/work/gmarques/cesm.cases/G/g.e30_a07c_cesm.GJRAv4.TL319_t232_wgx3_hycom1_N75.2025.130/'
[5]:
# The following parameters must be set accordingly
######################################################
# add your name and email address below
author = 'Gustavo Marques (gmarques@ucar.edu)'
######################################################
# create an empty class object
class args:
pass
# load avg dates
avg = diag_config_yml['Avg']
args.infile = OUTDIR
args.monthly = casename+diag_config_yml['Fnames']['z']
args.sigma2 = casename+diag_config_yml['Fnames']['rho2']
args.static = casename+diag_config_yml['Fnames']['static']
args.geom = casename+diag_config_yml['Fnames']['geom']
args.start_date = avg['start_date']
args.end_date = avg['end_date']
args.case_name = casename
args.label = ''
args.savefigs = False
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 13
10 # load avg dates
11 avg = diag_config_yml['Avg']
---> 13 args.infile = OUTDIR
14 args.monthly = casename+diag_config_yml['Fnames']['z']
15 args.sigma2 = casename+diag_config_yml['Fnames']['rho2']
NameError: name 'OUTDIR' is not defined
[6]:
# read grid info
geom_file = OUTDIR+'/'+args.geom
if os.path.exists(geom_file):
grd = MOM6grid(OUTDIR+'/'+args.static, geom_file)
else:
grd = MOM6grid(OUTDIR+'/'+args.static)
try:
depth = grd.depth_ocean
except:
depth = grd.deptho
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 2
1 # read grid info
----> 2 geom_file = OUTDIR+'/'+args.geom
3 if os.path.exists(geom_file):
4 grd = MOM6grid(OUTDIR+'/'+args.static, geom_file)
NameError: name 'OUTDIR' is not defined
[7]:
# remove Nan's, otherwise genBasinMasks won't work
depth[np.isnan(depth)] = 0.0
basin_code = genBasinMasks(grd.geolon, grd.geolat, depth, verbose=False)
basin_code_xr = genBasinMasks(grd.geolon, grd.geolat, depth, verbose=False, xda=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 # remove Nan's, otherwise genBasinMasks won't work
----> 2 depth[np.isnan(depth)] = 0.0
3 basin_code = genBasinMasks(grd.geolon, grd.geolat, depth, verbose=False)
4 basin_code_xr = genBasinMasks(grd.geolon, grd.geolat, depth, verbose=False, xda=True)
NameError: name 'depth' is not defined
[8]:
parallel, cluster, client = get_cluster(6, cluster_class='PBSCluster')
client
Starting a dask cluster: PBSCluster
Requesting 6 workers...
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[8], line 1
----> 1 parallel, cluster, client = get_cluster(6, cluster_class='PBSCluster')
2 client
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/mom6_tools/jobqueue.py:300, in get_cluster(nw, cluster_class, args, config, **kwargs)
297 print('Requesting {} workers... \n'.format(nw))
298 if is_jobqueue_cluster:
299 # dask.config.set({'distributed.dashboard.link': '/proxy/{port}/status'})
--> 300 cluster = cluster_class(**kwargs)
301 cluster.scale(nw)
302 else:
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/dask_jobqueue/core.py:656, in JobQueueCluster.__init__(self, n_workers, job_cls, loop, security, shared_temp_directory, silence_logs, name, asynchronous, dashboard_address, host, scheduler_options, scheduler_cls, interface, protocol, config_name, **job_kwargs)
653 self._job_kwargs = job_kwargs
655 worker = {"cls": self.job_cls, "options": self._job_kwargs}
--> 656 if "processes" in self._job_kwargs and self._job_kwargs["processes"] > 1:
657 worker["group"] = [
658 "-" + str(i) for i in range(self._job_kwargs["processes"])
659 ]
661 self._dummy_job # trigger property to ensure that the job is valid
TypeError: '>' not supported between instances of 'NoneType' and 'int'
[9]:
def preprocess(ds):
variables = ['vmo','vhml','vhGM']
for v in variables:
if v not in ds.variables:
ds[v] = xr.zeros_like(ds.vo)
return ds[variables]
[10]:
print('\n Reading dataset...')
# load data
%time ds = xr.open_mfdataset(OUTDIR+'/'+args.monthly, parallel=True, \
combine="nested", concat_dim="time", \
preprocess=preprocess).chunk({"time": 12})
Reading dataset...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
File <timed exec>:1
NameError: name 'OUTDIR' is not defined
[11]:
attrs = {
'description': 'Annual mean meridional thickness flux by components ',
'reduction_method': 'annual mean weighted by days in each month',
'casename': casename
}
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 4
1 attrs = {
2 'description': 'Annual mean meridional thickness flux by components ',
3 'reduction_method': 'annual mean weighted by days in each month',
----> 4 'casename': casename
5 }
NameError: name 'casename' is not defined
[12]:
ds_ann = m6toolbox.weighted_temporal_mean_vars(ds,attrs=attrs)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 ds_ann = m6toolbox.weighted_temporal_mean_vars(ds,attrs=attrs)
NameError: name 'ds' is not defined
[13]:
print('\n Selecting data between {} and {}...'.format(args.start_date, args.end_date))
%time ds_sel = ds_ann.sel(time=slice(args.start_date, args.end_date))
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[13], line 1
----> 1 print('\n Selecting data between {} and {}...'.format(args.start_date, args.end_date))
2 get_ipython().run_line_magic('time', 'ds_sel = ds_ann.sel(time=slice(args.start_date, args.end_date))')
AttributeError: type object 'args' has no attribute 'start_date'
Compute temporal mean for each term
[14]:
stream = True
# create a ndarray subclass
class C(np.ndarray): pass
[15]:
print('\n Computing time mean...')
%time ds_mean = ds_sel.mean('time').compute()
Computing time mean...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
File <timed exec>:1
NameError: name 'ds_sel' is not defined
[16]:
# create a ndarray subclass
class C(np.ndarray): pass
if 'vmo' in ds.variables:
varName = 'vmo'; conversion_factor = 1.e-9
elif 'vh' in ds.variables:
varName = 'vh'; conversion_factor = 1.e-6
if 'zw' in ds.variables: conversion_factor = 1.e-9 # Backwards compatible for when we had wrong units for 'vh'
else: raise Exception('Could not find "vh" or "vmo" in file "%s"'%(args.infile+args.monthly))
tmp = np.ma.masked_invalid(ds_sel[varName].mean('time').values)
tmp = tmp[:].filled(0.)
VHmod = tmp.view(C)
VHmod.units = ds_sel[varName].units
Zmod = m6toolbox.get_z(ds, depth, varName)
if args.case_name != '': case_name = args.case_name + ' ' + args.label
else: case_name = rootGroup.title + ' ' + args.label
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 4
1 # create a ndarray subclass
2 class C(np.ndarray): pass
----> 4 if 'vmo' in ds.variables:
5 varName = 'vmo'; conversion_factor = 1.e-9
6 elif 'vh' in ds.variables:
NameError: name 'ds' is not defined
Global MOC
[17]:
%matplotlib inline
# Global MOC
m6plot.setFigureSize([16,9],576,debug=False)
axis = plt.gca()
cmap = plt.get_cmap('dunnePM')
zg = Zmod.min(axis=-1); psiPlot = MOCpsi(VHmod)*conversion_factor
psiPlot = 0.5 * (psiPlot[0:-1,:]+psiPlot[1::,:])
yyg = grd.geolat_c[:,:].max(axis=-1)+0*zg
ci=m6plot.pmCI(0.,40.,5.)
plotPsi(yyg, zg, psiPlot, ci, 'Global MOC [Sv]')
plt.xlabel(r'Latitude [$\degree$N]')
plt.suptitle(case_name)
findExtrema(yyg, zg, psiPlot, max_lat=-30.)
findExtrema(yyg, zg, psiPlot, min_lat=25., min_depth=250.)
findExtrema(yyg, zg, psiPlot, min_depth=2000., mult=-1.)
plt.gca().invert_yaxis()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 7
5 axis = plt.gca()
6 cmap = plt.get_cmap('dunnePM')
----> 7 zg = Zmod.min(axis=-1); psiPlot = MOCpsi(VHmod)*conversion_factor
8 psiPlot = 0.5 * (psiPlot[0:-1,:]+psiPlot[1::,:])
9 yyg = grd.geolat_c[:,:].max(axis=-1)+0*zg
NameError: name 'Zmod' is not defined
[18]:
# create dataset to store results
moc = xr.Dataset(data_vars={ 'moc' : (('z_l','yq'), psiPlot),
'amoc' : (('z_l','yq'), np.zeros((psiPlot.shape))),
'moc_FFM' : (('z_l','yq'), np.zeros((psiPlot.shape))),
'moc_GM' : (('z_l','yq'), np.zeros((psiPlot.shape))),
'amoc_45' : (('time'), np.zeros((ds_ann.time.shape))),
'moc_GM_ACC' : (('time'), np.zeros((ds_ann.time.shape))),
'amoc_26' : (('time'), np.zeros((ds_ann.time.shape))) },
coords={'z_l': ds.z_l, 'yq':ds.yq, 'time':ds_ann.time})
attrs = {'description': 'MOC time-mean sections and time-series', 'unit': 'Sv', 'start_date': avg['start_date'],
'end_date': avg['end_date']}
add_global_attrs(moc,attrs)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 2
1 # create dataset to store results
----> 2 moc = xr.Dataset(data_vars={ 'moc' : (('z_l','yq'), psiPlot),
3 'amoc' : (('z_l','yq'), np.zeros((psiPlot.shape))),
4 'moc_FFM' : (('z_l','yq'), np.zeros((psiPlot.shape))),
5 'moc_GM' : (('z_l','yq'), np.zeros((psiPlot.shape))),
6 'amoc_45' : (('time'), np.zeros((ds_ann.time.shape))),
7 'moc_GM_ACC' : (('time'), np.zeros((ds_ann.time.shape))),
8 'amoc_26' : (('time'), np.zeros((ds_ann.time.shape))) },
9 coords={'z_l': ds.z_l, 'yq':ds.yq, 'time':ds_ann.time})
10 attrs = {'description': 'MOC time-mean sections and time-series', 'unit': 'Sv', 'start_date': avg['start_date'],
11 'end_date': avg['end_date']}
12 add_global_attrs(moc,attrs)
NameError: name 'psiPlot' is not defined
[19]:
print('Saving netCDF files...')
os.makedirs('ncfiles', exist_ok=True)
moc.to_netcdf('ncfiles/'+str(casename)+'_MOC.nc')
Saving netCDF files...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 3
1 print('Saving netCDF files...')
2 os.makedirs('ncfiles', exist_ok=True)
----> 3 moc.to_netcdf('ncfiles/'+str(casename)+'_MOC.nc')
NameError: name 'moc' is not defined
Atlantic MOC
[20]:
m6plot.setFigureSize([16,9],576,debug=False)
cmap = plt.get_cmap('dunnePM')
m = 0*basin_code; m[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)]=1
ci=m6plot.pmCI(0.,22.,2.)
z = (m*Zmod).min(axis=-1); psiPlot = MOCpsi(VHmod, vmsk=m*np.roll(m,-1,axis=-2))*conversion_factor
psiPlot = 0.5 * (psiPlot[0:-1,:]+psiPlot[1::,:])
yy = grd.geolat_c[:,:].max(axis=-1)+0*z
plotPsi(yy, z, psiPlot, ci, 'Atlantic MOC [Sv]')
plt.xlabel(r'Latitude [$\degree$N]')
plt.suptitle(case_name)
findExtrema(yy, z, psiPlot, min_lat=26.5, max_lat=27., min_depth=250.) # RAPID
findExtrema(yy, z, psiPlot, min_lat=44, max_lat=46., min_depth=250.) # RAPID
findExtrema(yy, z, psiPlot, max_lat=-33.)
findExtrema(yy, z, psiPlot)
findExtrema(yy, z, psiPlot, min_lat=5.)
plt.gca().invert_yaxis()
moc['amoc'].data = psiPlot
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 3
1 m6plot.setFigureSize([16,9],576,debug=False)
2 cmap = plt.get_cmap('dunnePM')
----> 3 m = 0*basin_code; m[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)]=1
4 ci=m6plot.pmCI(0.,22.,2.)
5 z = (m*Zmod).min(axis=-1); psiPlot = MOCpsi(VHmod, vmsk=m*np.roll(m,-1,axis=-2))*conversion_factor
NameError: name 'basin_code' is not defined
<Figure size 1024x576 with 0 Axes>
AMOC profile at 26N
[21]:
catalog = intake.open_catalog(diag_config_yml['oce_cat'])
rapid_vertical = catalog["moc-rapid"].to_dask()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[21], line 1
----> 1 catalog = intake.open_catalog(diag_config_yml['oce_cat'])
2 rapid_vertical = catalog["moc-rapid"].to_dask()
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/__init__.py:186, in open_catalog(uri, **kwargs)
179 raise ValueError(
180 f"Unknown catalog driver '{driver}'. "
181 "Do you need to install a new driver from the plugin directory? "
182 "https://intake.readthedocs.io/en/latest/plugin-directory.html\n"
183 f"Current registry: {list(sorted(registry))}"
184 )
185 try:
--> 186 return registry[driver](uri, **kwargs)
187 except VersionError:
188 # warn that we are switching to V2? The file will be read twice
189 return from_yaml_file(uri, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:617, in YAMLFileCatalog.__init__(self, path, text, autoreload, **kwargs)
615 self.filesystem = kwargs.pop("fs", None)
616 self.access = "name" not in kwargs
--> 617 super(YAMLFileCatalog, self).__init__(**kwargs)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/base.py:128, in Catalog.__init__(self, entries, name, description, metadata, ttl, getenv, getshell, persist_mode, storage_options, user_parameters)
126 self.updated = time.time()
127 self._entries = entries if entries is not None else self._make_entries_container()
--> 128 self.force_reload()
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/base.py:186, in Catalog.force_reload(self)
184 """Imperative reload data now"""
185 self.updated = time.time()
--> 186 self._load()
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:652, in YAMLFileCatalog._load(self, reload)
650 logger.warning("Use of '!template' deprecated - fixing")
651 text = text.replace("!template ", "")
--> 652 self.parse(text)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:730, in YAMLFileCatalog.parse(self, text)
728 # Second, we validate the schema and semantics
729 context = dict(root=self._dir)
--> 730 result = CatalogParser(data, context=context, getenv=self.getenv, getshell=self.getshell)
731 if result.errors:
732 raise exceptions.ValidationError(
733 "Catalog '{}' has validation errors:\n\n{}"
734 "".format(self.path, "\n".join(result.errors)),
735 result.errors,
736 )
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:342, in CatalogParser.__init__(self, data, getenv, getshell, context)
340 self.getenv = getenv
341 self.getshell = getshell
--> 342 self._data = self._parse(data)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:565, in CatalogParser._parse(self, data)
561 if (data.get("version", None) or data.get("metadata", {}).get("version", None) or 1) > 1:
562 raise VersionError("Not a V1 Catalog; perhaps use intake.open_catalog")
564 return dict(
--> 565 plugin_sources=self._parse_plugins(data),
566 data_sources=self._parse_data_sources(data),
567 metadata=data.get("metadata", {}),
568 name=data.get("name"),
569 description=data.get("description"),
570 )
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:400, in CatalogParser._parse_plugins(self, data)
397 elif "module" in plugin_source:
398 import intake
--> 400 intake.import_name(plugin_source["module"])
401 elif "dir" in plugin_source:
402 self.error(
403 "The key 'dir', and in general the feature of registering "
404 "plugins from a directory of Python scripts outside of "
405 "sys.path, is no longer supported. Use 'module'.",
406 plugin_source,
407 )
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/utils.py:27, in import_name(name)
25 modname = name.split(":", 1)[0]
26 logger.debug("Importing: '%s'" % modname)
---> 27 mod = importlib.import_module(modname)
28 if ":" in name:
29 end = name.split(":")[1]
File ~/.asdf/installs/python/3.10.20/lib/python3.10/importlib/__init__.py:126, in import_module(name, package)
124 break
125 level += 1
--> 126 return _bootstrap._gcd_import(name[level:], package, level)
File <frozen importlib._bootstrap>:1050, in _gcd_import(name, package, level)
File <frozen importlib._bootstrap>:1027, in _find_and_load(name, import_)
File <frozen importlib._bootstrap>:1004, in _find_and_load_unlocked(name, import_)
ModuleNotFoundError: No module named 'intake_xarray'
[22]:
if 'zl' in ds:
zl=ds.zl.values
elif 'z_l' in ds:
zl=ds.z_l.values
else:
raise ValueError("Dataset does not have vertical coordinate zl or z_l")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 1
----> 1 if 'zl' in ds:
2 zl=ds.zl.values
3 elif 'z_l' in ds:
NameError: name 'ds' is not defined
[23]:
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
ax.plot(rapid_vertical.stream_function_mar.mean('time'), rapid_vertical.depth, 'k', label='RAPID')
ax.plot(moc['amoc'].sel(yq=26, method='nearest'), moc.z_l, label=case_name)
ax.legend()
plt.gca().invert_yaxis()
plt.grid()
ax.set_xlabel('AMOC @ 26N [Sv]')
ax.set_ylabel('Depth [m]');
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 2
1 fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
----> 2 ax.plot(rapid_vertical.stream_function_mar.mean('time'), rapid_vertical.depth, 'k', label='RAPID')
3 ax.plot(moc['amoc'].sel(yq=26, method='nearest'), moc.z_l, label=case_name)
4 ax.legend()
NameError: name 'rapid_vertical' is not defined
AMOC time series
[24]:
dtime = ds_ann.time.values
amoc_26 = np.zeros(len(dtime))
amoc_45 = np.zeros(len(dtime))
moc_GM_ACC = np.zeros(len(dtime))
# loop in time
for t in range(len(dtime)):
tmp = np.ma.masked_invalid(ds_ann[varName].sel(time=dtime[t]).values)
tmp = tmp[:].filled(0.)
psi = MOCpsi(tmp, vmsk=m*np.roll(m,-1,axis=-2))*conversion_factor
psi = 0.5 * (psi[0:-1,:]+psi[1::,:])
amoc_26[t] = findExtrema(yy, z, psi, min_lat=26.5, max_lat=27., plot=False)
amoc_45[t] = findExtrema(yy, z, psi, min_lat=44., max_lat=46., plot=False)
tmp_GM = np.ma.masked_invalid(ds_ann['vhGM'][t,:].values)
tmp_GM = tmp_GM[:].filled(0.)
psiGM = MOCpsi(tmp_GM)*conversion_factor
psiGM = 0.5 * (psiGM[0:-1,:]+psiGM[1::,:])
moc_GM_ACC[t] = findExtrema(yyg, zg, psiGM, min_lat=-65., max_lat=-30, mult=-1., plot=False)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 1
----> 1 dtime = ds_ann.time.values
3 amoc_26 = np.zeros(len(dtime))
4 amoc_45 = np.zeros(len(dtime))
NameError: name 'ds_ann' is not defined
[25]:
# add dataarays to the moc dataset
moc['amoc_26'].data = amoc_26
moc['amoc_45'].data = amoc_45
moc['moc_GM_ACC'].data = moc_GM_ACC
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[25], line 2
1 # add dataarays to the moc dataset
----> 2 moc['amoc_26'].data = amoc_26
3 moc['amoc_45'].data = amoc_45
4 moc['moc_GM_ACC'].data = moc_GM_ACC
NameError: name 'amoc_26' is not defined
[26]:
# load datasets from oce catalog
amoc_core_26 = catalog["moc-core2-26p5"].to_dask()
amoc_pop_26 = catalog["moc-pop-jra-26"].to_dask()
rapid = m6toolbox.weighted_temporal_mean_vars(catalog["transports-rapid"].to_dask())
amoc_core_45 = catalog["moc-core2-45"].to_dask()
amoc_pop_45 = catalog["moc-pop-jra-45"].to_dask()
#list(catalog)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[26], line 2
1 # load datasets from oce catalog
----> 2 amoc_core_26 = catalog["moc-core2-26p5"].to_dask()
3 amoc_pop_26 = catalog["moc-pop-jra-26"].to_dask()
4 rapid = m6toolbox.weighted_temporal_mean_vars(catalog["transports-rapid"].to_dask())
NameError: name 'catalog' is not defined
AMOC @ 26 \(^o\) N
[27]:
# plot
fig = plt.figure(figsize=(12, 6))
plt.plot(np.arange(len(moc.time))+1958.5 ,moc['amoc_26'].values, color='k', label=case_name, lw=2)
# core data
core_mean = amoc_core_26['MOC'].mean(axis=0).data
core_std = amoc_core_26['MOC'].std(axis=0).data
plt.plot(amoc_core_26.time,core_mean, 'k', label='CORE II (group mean)', color='#1B2ACC', lw=1)
plt.fill_between(amoc_core_26.time, core_mean-core_std, core_mean+core_std,
alpha=0.25, edgecolor='#1B2ACC', facecolor='#089FFF')
# pop data
plt.plot(np.arange(len(amoc_pop_26.time))+1958.5 ,amoc_pop_26.AMOC_26n.values, color='r', label='POP', lw=1)
# rapid
plt.plot(np.arange(len(rapid.time))+2004.5 ,rapid.moc_mar_hc10.values, color='green', label='RAPID', lw=1)
#plt.plot(np.arange(len(rapid_filtered.time))+2004.5 ,rapid_filtered.values, color='green', label='RAPID', lw=1)
plt.title('AMOC @ 26 $^o$ N', fontsize=16)
plt.ylim(5,20)
plt.xlim(1948,1958.5+len(moc.time))
plt.xlabel('Time [years]', fontsize=16); plt.ylabel('Sv', fontsize=16)
plt.legend(fontsize=13, ncol=2)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[27], line 3
1 # plot
2 fig = plt.figure(figsize=(12, 6))
----> 3 plt.plot(np.arange(len(moc.time))+1958.5 ,moc['amoc_26'].values, color='k', label=case_name, lw=2)
4 # core data
5 core_mean = amoc_core_26['MOC'].mean(axis=0).data
NameError: name 'moc' is not defined
<Figure size 1200x600 with 0 Axes>
AMOC @ 45 \(^o\) N
[28]:
# plot
fig = plt.figure(figsize=(12, 6))
plt.plot(np.arange(len(moc.time))+1958.5 ,moc['amoc_45'].values, color='k', label=case_name, lw=2)
# core data
core_mean = amoc_core_45['MOC'].mean(axis=0).data
core_std = amoc_core_45['MOC'].std(axis=0).data
plt.plot(amoc_core_45.time,core_mean, 'k', label='CORE II (group mean)', color='#1B2ACC', lw=1)
plt.fill_between(amoc_core_45.time, core_mean-core_std, core_mean+core_std,
alpha=0.25, edgecolor='#1B2ACC', facecolor='#089FFF')
# pop data
plt.plot(np.arange(len(amoc_pop_45.time))+1958. ,amoc_pop_45.AMOC_45n.values, color='r', label='POP', lw=1)
plt.title('AMOC @ 45 $^o$ N', fontsize=16)
plt.ylim(5,20)
plt.xlim(1948,1958+len(moc.time))
plt.xlabel('Time [years]', fontsize=16); plt.ylabel('Sv', fontsize=16)
plt.legend(fontsize=13, ncol=3)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[28], line 3
1 # plot
2 fig = plt.figure(figsize=(12, 6))
----> 3 plt.plot(np.arange(len(moc.time))+1958.5 ,moc['amoc_45'].values, color='k', label=case_name, lw=2)
4 # core data
5 core_mean = amoc_core_45['MOC'].mean(axis=0).data
NameError: name 'moc' is not defined
<Figure size 1200x600 with 0 Axes>
Submesoscale-induced Global MOC
[29]:
# create a ndarray subclass
class C(np.ndarray): pass
if 'vhml' in ds.variables:
varName = 'vhml'; conversion_factor = 1.e-9
else: raise Exception('Could not find "vhml" in file "%s"'%(args.infile+args.monthly))
tmp = np.ma.masked_invalid(ds_mean[varName].values)
tmp = tmp[:].filled(0.)
VHmod = tmp.view(C)
VHmod.units = ds[varName].units
# Global MOC
m6plot.setFigureSize([16,9],576,debug=False)
axis = plt.gca()
cmap = plt.get_cmap('dunnePM')
z = Zmod.min(axis=-1); psiPlot = MOCpsi(VHmod)*conversion_factor
psiPlot = 0.5 * (psiPlot[0:-1,:]+psiPlot[1::,:])
#yy = y[1:,:].max(axis=-1)+0*z
yy = grd.geolat_c[:,:].max(axis=-1)+0*z
ci=m6plot.pmCI(0.,20.,2.)
plotPsi(yy, z, psiPlot, ci, 'Global MOC [Sv] due to vhML', zval=[0.,-400.,-6500.])
plt.xlabel(r'Latitude [$\degree$N]')
plt.suptitle(case_name)
plt.gca().invert_yaxis()
moc['moc_FFM'].data = psiPlot
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[29], line 4
1 # create a ndarray subclass
2 class C(np.ndarray): pass
----> 4 if 'vhml' in ds.variables:
5 varName = 'vhml'; conversion_factor = 1.e-9
6 else: raise Exception('Could not find "vhml" in file "%s"'%(args.infile+args.monthly))
NameError: name 'ds' is not defined
Eddy(GM)-induced Global MOC
[30]:
# create a ndarray subclass
class C(np.ndarray): pass
if 'vhGM' in ds.variables:
varName = 'vhGM'; conversion_factor = 1.e-9
else: raise Exception('Could not find "vhGM" in file "%s"'%(args.infile+args.monthly))
tmp = np.ma.masked_invalid(ds_mean[varName].values)
tmp = tmp[:].filled(0.)
VHmod = tmp.view(C)
VHmod.units = ds[varName].units
# Global MOC
m6plot.setFigureSize([16,9],576,debug=False)
axis = plt.gca()
cmap = plt.get_cmap('dunnePM')
z = Zmod.min(axis=-1); psiPlot = MOCpsi(VHmod)*conversion_factor
psiPlot = 0.5 * (psiPlot[0:-1,:]+psiPlot[1::,:])
yy = grd.geolat_c[:,:].max(axis=-1)+0*z
ci=m6plot.pmCI(0.,20.,1.)
plotPsi(yy, z, psiPlot, ci, 'Global MOC [Sv] due to GM')
plt.xlabel(r'Latitude [$\degree$N]')
plt.suptitle(case_name)
findExtrema(yy, z, psiPlot, min_lat=-65., max_lat=-30, mult=-1.)
plt.gca().invert_yaxis()
moc['moc_GM'].data = psiPlot
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[30], line 4
1 # create a ndarray subclass
2 class C(np.ndarray): pass
----> 4 if 'vhGM' in ds.variables:
5 varName = 'vhGM'; conversion_factor = 1.e-9
6 else: raise Exception('Could not find "vhGM" in file "%s"'%(args.infile+args.monthly))
NameError: name 'ds' is not defined
[31]:
print('Saving netCDF files...')
moc.to_netcdf('ncfiles/'+str(casename)+'_MOC.nc')
Saving netCDF files...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[31], line 2
1 print('Saving netCDF files...')
----> 2 moc.to_netcdf('ncfiles/'+str(casename)+'_MOC.nc')
NameError: name 'moc' is not defined
Sigma-2 space (to be implemented…)
[32]:
def calc_moc_rho(vmo):
# Sum over the zonal direction and integrate along density
integ_layers = (
vmo.sum("xh").cumsum("rho2_l") - vmo.sum("xh").sum("rho2_l")
) / rho0 / 1.0e6 + 0.1
# The result of the integration over layers is evaluated at the interfaces
# with psi = 0 as the bottom boundary condition for the integration
bottom_condition = xr.zeros_like(integ_layers.isel({"rho2_l": 0}))
# combine bottom condition with data array
# psi_raw = xr.concat([integ_layers, bottom_condition], dim='rho2_l')
psi_raw = xr.concat([bottom_condition, integ_layers], dim="rho2_l")
# rename to correct dimension and add correct vertical coordinate
psi = psi_raw.rename({"rho2_l": "rho2_i"}).transpose("rho2_i", "yq")
psi["rho2_i"] = xr.concat([vmo.rho2_l[0]*0, vmo.rho2_l], dim="rho2_l").rename({"rho2_l": "rho2_i"})
#psi = psi.assign_coords(rho2_i=rho2_i)
psi.name = "psi"
return psi.load()
[33]:
print('\n Reading dataset...')
# load data
%time ds_sigma2 = xr.open_mfdataset(OUTDIR+'/'+args.sigma2, parallel=True, \
combine="nested", concat_dim="time", \
preprocess=preprocess).chunk({"time": 12})
Reading dataset...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
File <timed exec>:1
NameError: name 'OUTDIR' is not defined
[34]:
ds_ann_sigma2 = m6toolbox.weighted_temporal_mean_vars(ds_sigma2,attrs=attrs)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[34], line 1
----> 1 ds_ann_sigma2 = m6toolbox.weighted_temporal_mean_vars(ds_sigma2,attrs=attrs)
NameError: name 'ds_sigma2' is not defined
[35]:
print('\n Selecting data between {} and {}...'.format(args.start_date, args.end_date))
%time ds_sel_sigma2 = ds_ann_sigma2.sel(time=slice(args.start_date, args.end_date))
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[35], line 1
----> 1 print('\n Selecting data between {} and {}...'.format(args.start_date, args.end_date))
2 get_ipython().run_line_magic('time', 'ds_sel_sigma2 = ds_ann_sigma2.sel(time=slice(args.start_date, args.end_date))')
AttributeError: type object 'args' has no attribute 'start_date'
[36]:
print('\n Computing time mean...')
%time ds_mean_sigma2 = ds_sel_sigma2.mean('time').compute()
Computing time mean...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
File <timed exec>:1
NameError: name 'ds_sel_sigma2' is not defined
[37]:
# The following is from John Krasting's notebook. Thanks, John!
# https://github.com/jkrasting/mar/blob/main/src/gfdlnb/notebooks/ocean/Global_Meridional_Overturning.ipynb
rho0 = 1035.
ax = plt.subplot(1, 1, 1)
vmo = ds_mean_sigma2.vmo.where(ds_mean_sigma2.vmo < 1e14)
psi = calc_moc_rho(vmo)
levels = np.arange(-40, 42, 2)
cb = ax.contourf(psi.yq, psi.rho2_i, psi, levels=levels, cmap="RdBu_r")
cs = ax.contour(psi.yq, psi.rho2_i, psi, levels=levels, colors="k", linewidths=0.3)
stats = {}
y = psi.yq
r = psi.rho2_i
Y, R = np.meshgrid(y, r)
mask = (Y >= -75) & (Y <= -50) & (R >= 1036.8) & (R <= 1037.4)
_psi = psi.values
min_psi = np.min(_psi[mask])
min_index = np.unravel_index(np.argmin(_psi[mask]), _psi[mask].shape)
min_y = Y[mask][min_index]
min_r = R[mask][min_index]
stats = {
"min_psi": round(float(min_psi),2),
"min_lat": round(float(min_y),2),
"min_rho": round(float(min_r),2),
}
square_size = 10
ax.plot(
min_y,
min_r,
marker="s",
color="magenta",
markersize=square_size + 5,
markeredgewidth=2,
markeredgecolor="magenta",
markerfacecolor="none",
)
ax.annotate(
f"{min_psi:.1f} Sv",
xy=(min_y, min_r),
xytext=(10, 10),
textcoords="offset points",
arrowprops=dict(arrowstyle="->", color="black"),
bbox=dict(
facecolor="white",
edgecolor="black",
boxstyle="round,pad=0.2",
linewidth=0.5,
alpha=0.7,
),
)
ax.set_yscale("splitscale", zval=[1037.25, 1036.5, 1028.8])
plt.colorbar(cb)
date_range = 'Years '+args.start_date+' to '+ args.end_date
ax.text(.99,1.03,date_range,ha="right", transform=ax.transAxes)
ax.text(.01,1.03,f"Global",ha="left", transform=ax.transAxes)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[37], line 6
4 rho0 = 1035.
5 ax = plt.subplot(1, 1, 1)
----> 6 vmo = ds_mean_sigma2.vmo.where(ds_mean_sigma2.vmo < 1e14)
7 psi = calc_moc_rho(vmo)
9 levels = np.arange(-40, 42, 2)
NameError: name 'ds_mean_sigma2' is not defined
[38]:
# release workers
client.close(); cluster.close()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[38], line 2
1 # release workers
----> 2 client.close(); cluster.close()
NameError: name 'client' is not defined