update to TEOS-10

etc/config.move.g.e20.G.TL319_t13.control.001

@@ -60,6 +60,7 @@ time_avg=monthly
 
 [options]
 array = MOVE
+eos = teos10
 td_geo = False
 endpoint = False
 georef_level = 4625.

etc/config.rapid.g.e20.G.TL319_t13.control.001

@@ -75,6 +75,7 @@ time_avg=monthly
 
 [options]
 array = RAPID
+eos = teos10
 td_geo = False
 endpoint = False
 #georef_level = 4625.

etc/config.samba.g.e20.G.TL319_t13.control.001

@@ -72,6 +72,7 @@ time_avg=monthly
 
 [options]
 array = SAMBA
+eos = teos10
 td_geo = False
 endpoint = False
 georef_level = 4700.

metric/geostrophy.py

@@ -6,17 +6,18 @@
 import numpy as np
 import copy
 
+import gsw
+
 from . import constants
 from . import utils
 
 
 
-def eos_insitu(t, s, p):
+def eos80_insitu_dens(t, s, p):
     """
-    Returns in situ density of seawater as calculated by the NEMO
-    routine eos_insitu.f90. Computes the density referenced to
-    a specified depth/pressure from potential temperature and salinity
-    using the Jackett and McDougall (1994) equation of state.
+    Computes insitu density from potential temperature and salinity using the
+    1980 Equation of State of Seawater (EOS80; Jackett and McDougall 1994).
+    https://doi.org/10.1175/1520-0426(1995)012<0381:MAOHPT>2.0.CO;2
 
     """
     # Convert to double precision
@@ -63,7 +64,27 @@ def eos_insitu(t, s, p):
     return rho
 
 
-def calc_dh(t_on_v, s_on_v, get_rho=True):
+def teos10_insitu_dens(t, s, z, lat, lon):
+    """
+    Computes the insitu density from potential temperature and salinity using the
+    Thermodynamic Equation of Seawater 2010 (TEOS-10; IOC, SCOR and IAPSO, 2010).
+    http://www.teos-10.org/pubs/TEOS-10_Manual.pdf
+
+    """
+
+    depth = np.ones_like(t) * z[None,:,None]
+    lat = np.ones_like(t) * lat[None,None,:]
+    lon = np.ones_like(t) * lon[None,None,:]
+
+    p = gsw.p_from_z(-depth, lat)
+    SA = gsw.SA_from_SP(s, p, lon, lat)
+    CT = gsw.CT_from_pt(SA, t)
+    rho = gsw.rho(SA,CT,p)
+
+    return rho
+
+
+def calc_dh(t_on_v, s_on_v, eos, get_rho=True):
     """
     Return ZonalSections containing dynamic heights calculated from
     from temperature and salinity interpolated onto velocity boundaries.
@@ -72,8 +93,14 @@ def calc_dh(t_on_v, s_on_v, get_rho=True):
     # Calculate in situ density at bounds
     rho = copy.deepcopy(t_on_v)
     rho.data = None # Density not needed at v mid-points
-    rho.bounds_data = eos_insitu(t_on_v.bounds_data, s_on_v.bounds_data,
-                                 t_on_v.z_as_bounds_data)
+    if eos == 'eos80':
+      rho.bounds_data = eos80_insitu_dens(t_on_v.bounds_data, s_on_v.bounds_data,
+                                          t_on_v.z_as_bounds_data)
+    elif eos == 'teos10':
+      rho.bounds_data = teos10_insitu_dens(t_on_v.bounds_data, s_on_v.bounds_data,
+                                           t_on_v.z, t_on_v.ybounds, t_on_v.xbounds)
+    else:
+      raise RuntimeError('EOS based on either the 1980 equation of state (eos="eos80") or the international thermodynamic equation of seawater - 2010 (eos="teos10") are implemented')
 
     # Calculate dynamic height relative to a reference level
     dh = copy.deepcopy(rho)

metric/metric.py

@@ -29,6 +29,10 @@ def compute_amoc_transport(line_args):
     if args.shift_date is not None:
         config.set('output', 'shift_date', args.shift_date)
 
+    # Update shift_date in config file
+    if utils.get_config_opt(config, 'options', 'eos') is None:
+        config.set('options', 'eos', 'teos10')
+
     # Check files and options 
     if not utils.path.exists(args.config_file):
       raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), args.config_file)

metric/move/output.py

@@ -52,6 +52,7 @@ def create_netcdf(config, move_trans, model_trans, bdry_trans, int_trans, int_mo
     bdry_maxlon = config.getfloat('options','bdry_maxlon')
     int_maxlon = config.getfloat('options','int_maxlon')
     georef = config.getfloat('options','georef_level')
+    eos = config.get('options','eos')
 
     try:
         vref_level =  config.getfloat('options','vref_level')
@@ -111,6 +112,7 @@ def create_netcdf(config, move_trans, model_trans, bdry_trans, int_trans, int_mo
       dataset.reference_to_model_velocity = vref_level
     if config.getboolean('options', 'endpoint'):
       dataset.geostrophic_computation = 'endpoint'
+    dataset.eos = eos
     dataset.rhocp = move_trans.rhocp
     dataset.contact = 'fredc.ucar.edu'
     dataset.code_reference = 'https://github.com/NCAR/metric'

metric/move/transports.py

@@ -273,7 +273,8 @@ def calc_transports_from_sections(config, v, t_on_v, s_on_v, ssh_on_v):
     intmin,intmax = utils.get_indrange(v.x, bdry_maxlon, int_maxlon)         # boundary component
 
     # Calculate dynamic heights
-    dh, rho = calc_dh(t_on_v, s_on_v, get_rho=True)
+    eos = config.get('options','eos')
+    dh, rho = calc_dh(t_on_v, s_on_v, eos, get_rho=True)
 
     # Calculate geostrophic transports
     if config.getboolean('options', 'td_geo'):

metric/rapid/output.py

@@ -54,6 +54,7 @@ def create_netcdf(config, rapid_trans, model_trans, fc_trans,
     int_maxlon = config.getfloat('options','int_maxlon')
     georef = config.getfloat('options','georef_level')
     ek_level = config.getfloat('options','ekman_depth')
+    eos = config.get('options','eos')
 
     try:
         vref_level =  config.getfloat('options','vref_level')
@@ -113,6 +114,7 @@ def create_netcdf(config, rapid_trans, model_trans, fc_trans,
       dataset.reference_to_model_velocity = vref_level
     if config.getboolean('options', 'endpoint'):
       dataset.geostrophic_computation = 'endpoint'
+    dataset.eos = eos
     dataset.rhocp = rapid_trans.rhocp
     dataset.ekman_level = ek_level
     dataset.contact = 'fredc.ucar.edu'

metric/rapid/transports.py

@@ -271,7 +271,8 @@ def calc_transports_from_sections(config, v, tau, t_on_v, s_on_v, ssh_on_v):
     intmin, intmax = utils.get_indrange(v.x, wbw_maxlon, int_maxlon) # Gyre interior
 
     # Calculate dynamic heights
-    dh, rho = calc_dh(t_on_v, s_on_v, get_rho=True)
+    eos = config.get('options','eos')
+    dh, rho = calc_dh(t_on_v, s_on_v, eos, get_rho=True)
 
     # Calculate geostrophic transports
     if config.getboolean('options', 'td_geo'):

metric/samba/output.py

@@ -54,6 +54,7 @@ def create_netcdf(config, samba_trans, model_trans, wbw_trans,
     ebw_maxlon = config.getfloat('options','ebw_maxlon')
     georef = config.getfloat('options','georef_level')
     ek_level = config.getfloat('options','ekman_depth')
+    eos = config.get('options','eos')
 
     try:
         vref_level =  config.getfloat('options','vref_level')
@@ -113,6 +114,7 @@ def create_netcdf(config, samba_trans, model_trans, wbw_trans,
       dataset.reference_to_model_velocity = vref_level
     if config.getboolean('options', 'endpoint'):
       dataset.geostrophic_computation = 'endpoint'
+    dataset.eos = eos
     dataset.rhocp = samba_trans.rhocp
     dataset.ekman_level = ek_level
     dataset.contact = 'fredc.ucar.edu'

metric/samba/transports.py

@@ -271,7 +271,8 @@ def calc_transports_from_sections(config, v, tau, t_on_v, s_on_v, ssh_on_v):
     ebwmin, ebwmax = utils.get_indrange(v.x, int_maxlon, ebw_maxlon) # EBW
 
     # Calculate dynamic heights
-    dh, rho = calc_dh(t_on_v, s_on_v, get_rho=True)
+    eos = config.get('options','eos')
+    dh, rho = calc_dh(t_on_v, s_on_v, eos, get_rho=True)
 
     # Calculate geostrophic transports
     if config.getboolean('options', 'td_geo'):

requirements.txt

@@ -4,3 +4,4 @@ netCDF4
 ConfigParser
 argparse
 matplotlib
+gsw