Temperature and Salinity at Depth Levels
Goals of this notebook:
serve as an example of how to post-process CESM/MOM6 output;
create time averages of T/S fields at depth levels and compared agains observations (WOA18).
Temprature and salinity comparisons (model vs obs) at selected depth levels are grouped into the following regions: Global, Antarctic, and Arctic.
[1]:
%load_ext autoreload
%autoreload 2
[2]:
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")
from mom6_tools.MOM6grid import MOM6grid
from mom6_tools.m6plot import xycompare, polarcomparison
from mom6_tools.m6toolbox import cime_xmlquery, weighted_temporal_mean_vars
from mom6_tools import m6toolbox
from mom6_tools.jobqueue import get_cluster
import yaml, intake, os
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from IPython.display import display, Markdown, Latex
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
######################################################
# create an empty class object
class args:
pass
# load avg dates
avg = diag_config_yml['Avg']
args.start_date = avg['start_date']
args.end_date = avg['end_date']
args.casename = casename
args.obs = "woa-2018-tx2_3v2-annual-all"
args.monthly = casename+diag_config_yml['Fnames']['z']
args.static = casename+diag_config_yml['Fnames']['static']
args.geom = casename+diag_config_yml['Fnames']['geom']
args.savefigs = False
args.nw = 6 # requesting 6 workers
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 13
11 args.start_date = avg['start_date']
12 args.end_date = avg['end_date']
---> 13 args.casename = casename
14 args.obs = "woa-2018-tx2_3v2-annual-all"
15 args.monthly = casename+diag_config_yml['Fnames']['z']
NameError: name 'casename' is not defined
[6]:
# load mom6 grid
# read grid info
geom_file = OUTDIR+'/'+args.geom
if os.path.exists(geom_file):
grd_xr = MOM6grid(OUTDIR+'/'+args.static, geom_file, xrformat=True)
grd = MOM6grid(OUTDIR+'/'+args.static, geom_file)
else:
grd_xr = MOM6grid(OUTDIR+'/'+args.static, xrformat=True)
grd = MOM6grid(OUTDIR+'/'+args.static)
try:
area = grd_xr.area_t.where(grd_xr.wet > 0).values
except:
area = grd_xr.areacello.where(grd_xr.wet > 0).values
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 3
1 # load mom6 grid
2 # read grid info
----> 3 geom_file = OUTDIR+'/'+args.geom
4 if os.path.exists(geom_file):
5 grd_xr = MOM6grid(OUTDIR+'/'+args.static, geom_file, xrformat=True)
NameError: name 'OUTDIR' is not defined
[7]:
parallel, cluster, client = get_cluster(args.nw, cluster_class='PBSCluster')
client
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[7], line 1
----> 1 parallel, cluster, client = get_cluster(args.nw, cluster_class='PBSCluster')
2 client
AttributeError: type object 'args' has no attribute 'nw'
[8]:
client
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 1
----> 1 client
NameError: name 'client' is not defined
[9]:
# load history files
def preprocess(ds):
''' Return a dataset desired variables'''
variables = ['thetao', 'so']
return ds[variables]
ds = xr.open_mfdataset(OUTDIR+args.monthly, \
parallel=True, data_vars='minimal', \
coords='minimal', compat='override', preprocess=preprocess)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 8
5 variables = ['thetao', 'so']
6 return ds[variables]
----> 8 ds = xr.open_mfdataset(OUTDIR+args.monthly, \
9 parallel=True, data_vars='minimal', \
10 coords='minimal', compat='override', preprocess=preprocess)
NameError: name 'OUTDIR' is not defined
[10]:
# Select data
%time ds_sel = ds.sel(time=slice(args.start_date, args.end_date))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
File <timed exec>:1
NameError: name 'ds' is not defined
[11]:
# compute annual mean and then average in time
ds_ann = weighted_temporal_mean_vars(ds_sel)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 2
1 # compute annual mean and then average in time
----> 2 ds_ann = weighted_temporal_mean_vars(ds_sel)
NameError: name 'ds_sel' is not defined
[12]:
thetao_mean = ds_ann.thetao.mean('time')
temp = np.ma.masked_invalid(thetao_mean.values)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 thetao_mean = ds_ann.thetao.mean('time')
2 temp = np.ma.masked_invalid(thetao_mean.values)
NameError: name 'ds_ann' is not defined
[13]:
so_mean = ds_ann.so.mean('time')
salt = np.ma.masked_invalid(so_mean.values)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 so_mean = ds_ann.so.mean('time')
2 salt = np.ma.masked_invalid(so_mean.values)
NameError: name 'ds_ann' is not defined
[14]:
# load WOA18 data
catalog = intake.open_catalog(diag_config_yml['oce_cat'])
woa18 = catalog[args.obs].to_dask()
woa18 = woa18.rename({'z_l' : 'depth'});
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[14], line 2
1 # load WOA18 data
----> 2 catalog = intake.open_catalog(diag_config_yml['oce_cat'])
3 woa18 = catalog[args.obs].to_dask()
4 woa18 = woa18.rename({'z_l' : 'depth'});
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'
[15]:
print('Saving netCDF files...')
os.makedirs('ncfiles', exist_ok=True)
attrs = {'description': 'model - obs at depth levels',
'start_date': args.start_date,
'end_date': args.end_date,
'casename': casename,
'obs': args.obs,
}
Saving netCDF files...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 8
1 print('Saving netCDF files...')
3 os.makedirs('ncfiles', exist_ok=True)
5 attrs = {'description': 'model - obs at depth levels',
6 'start_date': args.start_date,
7 'end_date': args.end_date,
----> 8 'casename': casename,
9 'obs': args.obs,
10 }
NameError: name 'casename' is not defined
[16]:
temp_bias = np.ma.masked_invalid(thetao_mean.values - woa18.thetao.values)
ds_thetao = xr.Dataset(data_vars={ 'thetao' : (('z_l','yh','xh'), thetao_mean.values),
'thetao_bias' : (('z_l','yh','xh'), temp_bias)},
coords={'z_l' : ds.z_l, 'yh' : grd.yh, 'xh' : grd.xh})
m6toolbox.add_global_attrs(ds_thetao,attrs)
ds_thetao.to_netcdf('ncfiles/'+str(casename)+'_thetao_time_mean.nc')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 temp_bias = np.ma.masked_invalid(thetao_mean.values - woa18.thetao.values)
2 ds_thetao = xr.Dataset(data_vars={ 'thetao' : (('z_l','yh','xh'), thetao_mean.values),
3 'thetao_bias' : (('z_l','yh','xh'), temp_bias)},
4 coords={'z_l' : ds.z_l, 'yh' : grd.yh, 'xh' : grd.xh})
5 m6toolbox.add_global_attrs(ds_thetao,attrs)
NameError: name 'thetao_mean' is not defined
[17]:
so_bias = np.ma.masked_invalid(so_mean.values - woa18.so.values)
ds_so = xr.Dataset(data_vars={ 'so' : (('z_l','yh','xh'), so_mean.values),
'so_bias' : (('z_l','yh','xh'), so_bias)},
coords={'z_l' : ds.z_l, 'yh' : grd.yh, 'xh' : grd.xh})
m6toolbox.add_global_attrs(ds_so,attrs)
ds_so.to_netcdf('ncfiles/'+str(casename)+'_so_time_mean.nc')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 1
----> 1 so_bias = np.ma.masked_invalid(so_mean.values - woa18.so.values)
2 ds_so = xr.Dataset(data_vars={ 'so' : (('z_l','yh','xh'), so_mean.values),
3 'so_bias' : (('z_l','yh','xh'), so_bias)},
4 coords={'z_l' : ds.z_l, 'yh' : grd.yh, 'xh' : grd.xh})
5 m6toolbox.add_global_attrs(ds_so,attrs)
NameError: name 'so_mean' is not defined
[18]:
client.close(); cluster.close()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 client.close(); cluster.close()
NameError: name 'client' is not defined
Global
[19]:
%matplotlib inline
km = len(woa18['depth'])
for k in range(km):
if ds['z_l'][k].values < 1000.0:
temp_obs = np.ma.masked_invalid(woa18['thetao'][k,:].values)
xycompare(temp[k,:] , temp_obs, grd.geolon, grd.geolat, area=area,
title1 = 'model temperature, depth ='+str(ds['z_l'][k].values)+ 'm',
title2 = 'observed temperature, depth ='+str(woa18['depth'][k].values)+ 'm',
suptitle=casename + ', averaged '+str(args.start_date)+ ' and ' +str(args.end_date),
clim=(-1.9,30.), dcolormap=plt.cm.bwr,
extend='both', dextend='neither', dlim=(-2,2),
show= True)
salt_obs = np.ma.masked_invalid(woa18['so'][k,:].values)
xycompare( salt[k,:] , salt_obs, grd.geolon, grd.geolat, area=area,
title1 = 'model salinity, depth ='+str(ds['z_l'][k].values)+ 'm',
title2 = 'observed salinity, depth ='+str(woa18['depth'][k].values)+ 'm',
suptitle=casename, clim=(30,39.), dcolormap=plt.cm.bwr,
extend='both', dextend='neither', dlim=(-2,2),
show= True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 2
1 get_ipython().run_line_magic('matplotlib', 'inline')
----> 2 km = len(woa18['depth'])
3 for k in range(km):
4 if ds['z_l'][k].values < 1000.0:
NameError: name 'woa18' is not defined
Antarctic
[20]:
# loop over depths and compare TS fields
km = len(woa18['depth'])
for k in range(km):
if (ds['z_l'][k].values < 100.):
temp_obs = np.ma.masked_invalid(woa18['thetao'][k,:].values)
polarcomparison(temp[k,:] , temp_obs, grd,
title1 = 'model temperature, depth ='+str(ds['z_l'][k].values)+ 'm',
title2 = 'observed temperature, depth ='+str(woa18['depth'][k].values)+ 'm',
extend='both', dextend='neither', clim=(-1.9,10.5), dlim=(-2,2), dcolormap=plt.cm.bwr,
suptitle=casename + ', averaged '+str(args.start_date)+ ' to ' +str(args.end_date),
proj='SP', show= True)
salt_obs = np.ma.masked_invalid(woa18['so'][k,:].values)
polarcomparison( salt[k,:] , salt_obs, grd,
title1 = 'model salinity, depth ='+str(ds['z_l'][k].values)+ 'm',
title2 = 'observed salinity, depth ='+str(woa18['depth'][k].values)+ 'm',
extend='both', dextend='neither', clim=(33.,35.), dlim=(-2,2), dcolormap=plt.cm.bwr,
suptitle=casename + ', averaged '+str(args.start_date)+ ' to ' +str(args.end_date),
proj='SP', show= True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 2
1 # loop over depths and compare TS fields
----> 2 km = len(woa18['depth'])
3 for k in range(km):
4 if (ds['z_l'][k].values < 100.):
NameError: name 'woa18' is not defined
Arctic
[21]:
# loop over depths and compare TS fields
km = len(woa18['depth'])
for k in range(km):
if (ds['z_l'][k].values < 100.):
temp_obs = np.ma.masked_invalid(woa18['thetao'][k,:].values)
polarcomparison(temp[k,:] , temp_obs, grd,
title1 = 'model temperature, depth ='+str(ds['z_l'][k].values)+ 'm',
title2 = 'observed temperature, depth ='+str(woa18['depth'][k].values)+ 'm',
extend='both', dextend='neither', clim=(-1.9,10.5), dlim=(-2,2), dcolormap=plt.cm.bwr,
suptitle=casename + ', averaged '+str(args.start_date)+ ' to ' +str(args.end_date),
proj='NP', show= True)
salt_obs = np.ma.masked_invalid(woa18['so'][k,:].values)
polarcomparison( salt[k,:] , salt_obs, grd,
title1 = 'model salinity, depth ='+str(ds['z_l'][k].values)+ 'm',
title2 = 'observed salinity, depth ='+str(woa18['depth'][k].values)+ 'm',
extend='both', dextend='neither', clim=(32.,34.5), dlim=(-2,2), dcolormap=plt.cm.bwr,
suptitle=casename + ', averaged '+str(args.start_date)+ ' to ' +str(args.end_date),
proj='NP', show= True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 2
1 # loop over depths and compare TS fields
----> 2 km = len(woa18['depth'])
3 for k in range(km):
4 if (ds['z_l'][k].values < 100.):
NameError: name 'woa18' is not defined