Merge pull request #44 from efiring/salinometer

Add new salinometer functions from GSW-C.

README.rst

@@ -30,8 +30,8 @@ module to wrap the C functions as numpy
 `ufuncs <https://docs.scipy.org/doc/numpy/reference/ufuncs.html>`__,
 and then use an
 autogenerated Python module to add docstrings and handle masked
-arrays.  161 scalar C functions with only double-precision
-arguments and return values are wrapped as ufuncs, and 155 of
+arrays.  165 scalar C functions with only double-precision
+arguments and return values are wrapped as ufuncs, and 158 of
 these are exposed in the ``gsw`` namespace with an additional
 wrapper in Python.
 

gsw/_wrapped_ufuncs.py

@@ -2222,6 +2222,62 @@ def melting_seaice_SA_CT_ratio_poly(SA, CT, p, SA_seaice, t_seaice):
     """
     return _gsw_ufuncs.melting_seaice_sa_ct_ratio_poly(SA, CT, p, SA_seaice, t_seaice)
 
+@match_args_return
+def O2sol(SA, CT, p, lon, lat):
+    """
+    Calculates the oxygen concentration expected at equilibrium with air at
+    an Absolute Pressure of 101325 Pa (sea pressure of 0 dbar) including
+    saturated water vapor.  This function uses the solubility coefficients
+    derived from the data of Benson and Krause (1984), as fitted by Garcia
+    and Gordon (1992, 1993).
+
+    Parameters
+    ----------
+    SA : array-like
+        Absolute Salinity, g/kg
+    CT : array-like
+        Conservative Temperature (ITS-90), degrees C
+    p : array-like
+        Sea pressure (absolute pressure minus 10.1325 dbar), dbar
+    lon : array-like
+        Longitude, -360 to 360 degrees
+    lat : array-like
+        Latitude, -90 to 90 degrees
+
+    Returns
+    -------
+    O2sol : array-like, umol/kg
+        solubility of oxygen in micro-moles per kg
+
+
+    """
+    return _gsw_ufuncs.o2sol(SA, CT, p, lon, lat)
+
+@match_args_return
+def O2sol_SP_pt(SP, pt):
+    """
+    Calculates the oxygen concentration expected at equilibrium with air at
+    an Absolute Pressure of 101325 Pa (sea pressure of 0 dbar) including
+    saturated water vapor.  This function uses the solubility coefficients
+    derived from the data of Benson and Krause (1984), as fitted by Garcia
+    and Gordon (1992, 1993).
+
+    Parameters
+    ----------
+    SP : array-like
+        Practical Salinity (PSS-78), unitless
+    pt : array-like
+        Potential temperature referenced to a sea pressure, degrees C
+
+    Returns
+    -------
+    O2sol : array-like, umol/kg
+        solubility of oxygen in micro-moles per kg
+
+
+    """
+    return _gsw_ufuncs.o2sol_sp_pt(SP, pt)
+
 @match_args_return
 def p_from_z(z, lat):
     """
@@ -3655,6 +3711,34 @@ def SP_from_Sstar(Sstar, p, lon, lat):
     """
     return _gsw_ufuncs.sp_from_sstar(Sstar, p, lon, lat)
 
+@match_args_return
+def SP_salinometer(Rt, t):
+    """
+    Calculates Practical Salinity SP from a salinometer, primarily using the
+    PSS-78 algorithm.  Note that the PSS-78 algorithm for Practical Salinity
+    is only valid in the range 2 < SP < 42.  If the PSS-78 algorithm
+    produces a Practical Salinity that is less than 2 then the Practical
+    Salinity is recalculated with a modified form of the Hill et al. (1986)
+    formula.  The modification of the Hill et al. (1986) expression is to
+    ensure that it is exactly consistent with PSS-78 at SP = 2.
+
+    Parameters
+    ----------
+    Rt : array-like
+        C(SP,t_68,0)/C(SP=35,t_68,0), unitless
+    t : array-like
+        In-situ temperature (ITS-90), degrees C
+
+    Returns
+    -------
+    SP : array-like, unitless
+        Practical Salinity on the PSS-78 scale
+        t may have dimensions 1x1 or Mx1 or 1xN or MxN, where Rt is MxN.
+
+
+    """
+    return _gsw_ufuncs.sp_salinometer(Rt, t)
+
 @match_args_return
 def specvol(SA, CT, p):
     """

src/_ufuncs.c

@@ -643,6 +643,7 @@ static void *data_melting_ice_equilibrium_sa_ct_ratio[] = {&gsw_melting_ice_equi
 static void *data_melting_ice_equilibrium_sa_ct_ratio_poly[] = {&gsw_melting_ice_equilibrium_sa_ct_ratio_poly};
 static void *data_melting_seaice_equilibrium_sa_ct_ratio[] = {&gsw_melting_seaice_equilibrium_sa_ct_ratio};
 static void *data_melting_seaice_equilibrium_sa_ct_ratio_poly[] = {&gsw_melting_seaice_equilibrium_sa_ct_ratio_poly};
+static void *data_o2sol_sp_pt[] = {&gsw_o2sol_sp_pt};
 static void *data_p_from_z[] = {&gsw_p_from_z};
 static void *data_pot_enthalpy_ice_freezing[] = {&gsw_pot_enthalpy_ice_freezing};
 static void *data_pot_enthalpy_ice_freezing_poly[] = {&gsw_pot_enthalpy_ice_freezing_poly};
@@ -657,6 +658,7 @@ static void *data_sigma2[] = {&gsw_sigma2};
 static void *data_sigma3[] = {&gsw_sigma3};
 static void *data_sigma4[] = {&gsw_sigma4};
 static void *data_sound_speed_ice[] = {&gsw_sound_speed_ice};
+static void *data_sp_salinometer[] = {&gsw_sp_salinometer};
 static void *data_specvol_ice[] = {&gsw_specvol_ice};
 static void *data_spiciness0[] = {&gsw_spiciness0};
 static void *data_spiciness1[] = {&gsw_spiciness1};
@@ -728,6 +730,7 @@ static void *data_sstar_from_sa[] = {&gsw_sstar_from_sa};
 static void *data_sstar_from_sp[] = {&gsw_sstar_from_sp};
 static void *data_melting_seaice_sa_ct_ratio[] = {&gsw_melting_seaice_sa_ct_ratio};
 static void *data_melting_seaice_sa_ct_ratio_poly[] = {&gsw_melting_seaice_sa_ct_ratio_poly};
+static void *data_o2sol[] = {&gsw_o2sol};
 static void *data_ct_first_derivatives[] = {&gsw_ct_first_derivatives};
 static void *data_entropy_first_derivatives[] = {&gsw_entropy_first_derivatives};
 static void *data_pot_enthalpy_ice_freezing_first_derivatives[] = {&gsw_pot_enthalpy_ice_freezing_first_derivatives};
@@ -1283,6 +1286,18 @@ PyMODINIT_FUNC PyInit__gsw_ufuncs(void)
     PyDict_SetItemString(d, "melting_seaice_equilibrium_sa_ct_ratio_poly", ufunc_ptr);
     Py_DECREF(ufunc_ptr);
 
+    ufunc_ptr = PyUFunc_FromFuncAndData(funcs_dd_d,
+                                    data_o2sol_sp_pt,
+                                    types_dd_d,
+                                    1, 2, 1,  // ndatatypes, nin, nout
+                                    PyUFunc_None,
+                                    "o2sol_sp_pt",
+                                    "o2sol_sp_pt_docstring",
+                                    0);
+
+    PyDict_SetItemString(d, "o2sol_sp_pt", ufunc_ptr);
+    Py_DECREF(ufunc_ptr);
+
     ufunc_ptr = PyUFunc_FromFuncAndData(funcs_dd_d,
                                     data_p_from_z,
                                     types_dd_d,
@@ -1451,6 +1466,18 @@ PyMODINIT_FUNC PyInit__gsw_ufuncs(void)
     PyDict_SetItemString(d, "sound_speed_ice", ufunc_ptr);
     Py_DECREF(ufunc_ptr);
 
+    ufunc_ptr = PyUFunc_FromFuncAndData(funcs_dd_d,
+                                    data_sp_salinometer,
+                                    types_dd_d,
+                                    1, 2, 1,  // ndatatypes, nin, nout
+                                    PyUFunc_None,
+                                    "sp_salinometer",
+                                    "sp_salinometer_docstring",
+                                    0);
+
+    PyDict_SetItemString(d, "sp_salinometer", ufunc_ptr);
+    Py_DECREF(ufunc_ptr);
+
     ufunc_ptr = PyUFunc_FromFuncAndData(funcs_dd_d,
                                     data_specvol_ice,
                                     types_dd_d,
@@ -2303,6 +2330,18 @@ PyMODINIT_FUNC PyInit__gsw_ufuncs(void)
     PyDict_SetItemString(d, "melting_seaice_sa_ct_ratio_poly", ufunc_ptr);
     Py_DECREF(ufunc_ptr);
 
+    ufunc_ptr = PyUFunc_FromFuncAndData(funcs_ddddd_d,
+                                    data_o2sol,
+                                    types_ddddd_d,
+                                    1, 5, 1,  // ndatatypes, nin, nout
+                                    PyUFunc_None,
+                                    "o2sol",
+                                    "o2sol_docstring",
+                                    0);
+
+    PyDict_SetItemString(d, "o2sol", ufunc_ptr);
+    Py_DECREF(ufunc_ptr);
+
     ufunc_ptr = PyUFunc_FromFuncAndData(funcs_dd_dd,
                                     data_ct_first_derivatives,
                                     types_dd_dd,

src/c_gsw/gsw_oceanographic_toolbox.c

@@ -6391,6 +6391,115 @@ gsw_nsquared(double *sa, double *ct, double *p, double *lat, int nz,
 }
 /*
 !==========================================================================
+function gsw_o2sol(sa, ct, p, lon, lat)
+!==========================================================================
+!
+! Calculates the oxygen concentration expected at equilibrium with air at
+! an Absolute Pressure of 101325 Pa (sea pressure of 0 dbar) including
+! saturated water vapor.  This function uses the solubility coefficients
+! derived from the data of Benson and Krause (1984), as fitted by Garcia
+! and Gordon (1992, 1993).
+!
+! Note that this algorithm has not been approved by IOC and is not work
+! from SCOR/IAPSO Working Group 127. It is included in the GSW
+! Oceanographic Toolbox as it seems to be oceanographic best practice.
+!
+! SA  :  Absolute Salinity of seawater                           [ g/kg ]
+! CT  :  Conservative Temperature of seawater (ITS-90)           [ deg C ]
+! p   :  sea pressure at which the melting occurs                [ dbar ]
+!         ( i.e. absolute pressure - 10.1325 dbar )
+! lat : latitude                                                 [deg]
+! lon : longitude                                                [deg]
+!
+! gsw_o2sol : olubility of oxygen in micro-moles per kg          [umol/kg]
+*/
+double
+gsw_o2sol(double sa, double ct, double p, double lon, double lat)
+{
+    GSW_TEOS10_CONSTANTS;
+    double sp, pt, pt68, x, y, o2sol,
+           a0, a1, a2, a3, a4, a5, b0, b1, b2, b3, c0;
+
+    sp = gsw_sp_from_sa(sa, p, lon, lat);
+    x = sp;
+    pt = gsw_pt_from_ct(sa, ct);
+
+    pt68 = pt*1.00024;
+
+    y = log((298.15 - pt68)/(gsw_t0 + pt68));
+
+    a0 =  5.80871;
+    a1 =  3.20291;
+    a2 =  4.17887;
+    a3 =  5.10006;
+    a4 = -9.86643e-2;
+    a5 =  3.80369;
+    b0 = -7.01577e-3;
+    b1 = -7.70028e-3;
+    b2 = -1.13864e-2;
+    b3 = -9.51519e-3;
+    c0 = -2.75915e-7;
+
+    o2sol = exp(a0 + y*(a1 + y*(a2 + y*(a3 + y*(a4 + a5*y))))
+                  + x*(b0 + y*(b1 + y*(b2 + b3*y)) + c0*x));
+
+    return o2sol;
+
+}
+/*
+!==========================================================================
+function gsw_o2sol_sp_pt(sp, pt)
+!==========================================================================
+!
+! Calculates the oxygen concentration expected at equilibrium with air at
+! an Absolute Pressure of 101325 Pa (sea pressure of 0 dbar) including
+! saturated water vapor.  This function uses the solubility coefficients
+! derived from the data of Benson and Krause (1984), as fitted by Garcia
+! and Gordon (1992, 1993).
+!
+! Note that this algorithm has not been approved by IOC and is not work
+! from SCOR/IAPSO Working Group 127. It is included in the GSW
+! Oceanographic Toolbox as it seems to be oceanographic best practice.
+!
+! SP  :  Practical Salinity  (PSS-78)                         [ unitless ]
+! pt  :  potential temperature (ITS-90) referenced               [ dbar ]
+!         to one standard atmosphere (0 dbar).
+!
+! gsw_o2sol_sp_pt : olubility of oxygen in micro-moles per kg     [umol/kg]
+*/
+double
+gsw_o2sol_sp_pt(double sp, double pt)
+{
+    GSW_TEOS10_CONSTANTS;
+    double pt68, x, y, o2sol,
+           a0, a1, a2, a3, a4, a5, b0, b1, b2, b3, c0;
+
+    x = sp;
+
+    pt68 = pt*1.00024;
+
+    y = log((298.15 - pt68)/(gsw_t0 + pt68));
+
+    a0 =  5.80871;
+    a1 =  3.20291;
+    a2 =  4.17887;
+    a3 =  5.10006;
+    a4 = -9.86643e-2;
+    a5 =  3.80369;
+    b0 = -7.01577e-3;
+    b1 = -7.70028e-3;
+    b2 = -1.13864e-2;
+    b3 = -9.51519e-3;
+    c0 = -2.75915e-7;
+
+    o2sol = exp(a0 + y*(a1 + y*(a2 + y*(a3 + y*(a4 + a5*y))))
+                  + x*(b0 + y*(b1 + y*(b2 + b3*y)) + c0*x));
+
+    return o2sol;
+
+}
+/*
+!==========================================================================
 function gsw_p_from_z(z,lat)
 !==========================================================================
 
@@ -9513,6 +9622,66 @@ gsw_sp_from_sstar(double sstar, double p, double lon, double lat)
 }
 /*
 !==========================================================================
+function gsw_sp_salinometer(rt,t)
+!==========================================================================
+!  Calculates Practical Salinity SP from a salinometer, primarily using the
+!  PSS-78 algorithm.  Note that the PSS-78 algorithm for Practical Salinity
+!  is only valid in the range 2 < SP < 42.  If the PSS-78 algorithm
+!  produces a Practical Salinity that is less than 2 then the Practical
+!  Salinity is recalculated with a modified form of the Hill et al. (1986)
+!  formula.  The modification of the Hill et al. (1986) expression is to
+!  ensure that it is exactly consistent with PSS-78 at SP = 2.
+!
+!  A laboratory salinometer has the ratio of conductivities, Rt, as an
+!  output, and the present function uses this conductivity ratio and the
+!  temperature t of the salinometer bath as the two input variables.
+!
+!  rt  = C(SP,t_68,0)/C(SP=35,t_68,0)                          [ unitless ]
+!  t   = temperature of the bath of the salinometer,
+!        measured on the ITS-90 scale (ITS-90)                 [ deg C ]
+!
+!  gsw_sp_salinometer = Practical Salinity on the PSS-78 scale [ unitless ]
+*/
+double
+gsw_sp_salinometer(double rt, double t)
+{
+  GSW_SP_COEFFICIENTS;
+  double t68, ft68, rtx, sp, hill_ratio,
+         x, sqrty, part1, part2, sp_hill_raw;
+
+  if (rt < 0){
+    return NAN;
+  }
+
+  t68 = t*1.00024;
+  ft68 = (t68 - 15)/(1 + k*(t68 - 15));
+
+  rtx = sqrt(rt);
+
+  sp = a0 + (a1 + (a2 + (a3 + (a4 + a5 * rtx) * rtx) * rtx) * rtx) * rtx
+    + ft68 * (b0 + (b1 + (b2+ (b3 + (b4 + b5 * rtx) * rtx) * rtx) * rtx) * rtx);
+
+    /*
+     ! The following section of the code is designed for SP < 2 based on the
+     ! Hill et al. (1986) algorithm.  This algorithm is adjusted so that it is
+     ! exactly equal to the PSS-78 algorithm at SP = 2.
+    */
+
+   if (sp < 2) {
+       hill_ratio  = gsw_hill_ratio_at_sp2(t);
+       x           = 400e0*rt;
+       sqrty       = 10e0*rtx;
+       part1       = 1e0 + x*(1.5e0 + x);
+       part2       = 1e0 + sqrty*(1e0 + sqrty*(1e0 + sqrty));
+       sp_hill_raw = sp - a0/part1 - b0*ft68/part2;
+       sp          = hill_ratio*sp_hill_raw;
+   }
+
+  return sp;
+
+}
+/*
+!==========================================================================
 function gsw_specvol(sa,ct,p)
 !==========================================================================
 !

src/c_gsw/gswteos-10.h

@@ -174,6 +174,8 @@ extern double gsw_melting_seaice_sa_ct_ratio_poly(double sa, double ct,
                 double p, double sa_seaice, double t_seaice);
 extern void   gsw_nsquared(double *sa, double *ct, double *p, double *lat,
                 int nz, double *n2, double *p_mid);
+extern double gsw_o2sol(double sa, double ct, double p, double lon, double lat);
+extern double gsw_o2sol_sp_pt(double sp, double pt);
 extern double gsw_pot_enthalpy_from_pt_ice(double pt0_ice);
 extern double gsw_pot_enthalpy_from_pt_ice_poly(double pt0_ice);
 extern double gsw_pot_enthalpy_ice_freezing(double sa, double p);
@@ -267,6 +269,7 @@ extern double gsw_sp_from_sa(double sa, double p, double lon, double lat);
 extern double gsw_sp_from_sk(double sk);
 extern double gsw_sp_from_sr(double sr);
 extern double gsw_sp_from_sstar(double sstar, double p,double lon,double lat);
+extern double gsw_sp_salinometer(double rt, double t);
 extern double gsw_spiciness0(double sa, double ct);
 extern double gsw_spiciness1(double sa, double ct);
 extern double gsw_spiciness2(double sa, double ct);

tools/docstring_parts.py

@@ -14,6 +14,7 @@
 t_seaice =
 "In-situ temperature of the sea ice at pressure p (ITS-90), degrees C",
 t = "In-situ temperature (ITS-90), degrees C",
+Rt = "C(SP,t_68,0)/C(SP=35,t_68,0), unitless",
 CT = "Conservative Temperature (ITS-90), degrees C",
 C = "Conductivity, mS/cm",
 p = "Sea pressure (absolute pressure minus 10.1325 dbar), dbar",