Merge pull request #60 from ocefpaf/add_gsw_infunnel

Add gsw infunnel

.pre-commit-config.yaml

@@ -1,14 +1,14 @@
 
 repos:
 - repo: https://github.com/codespell-project/codespell
-  rev: v2.2.1
+  rev: v2.2.6
   hooks:
     - id: codespell
       args:
         - --ignore-words-list=preformed
 
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.3.0
+  rev: v4.6.0
   hooks:
     - id: trailing-whitespace
     - id: end-of-file-fixer

gsw_check_functions.c

@@ -8,6 +8,7 @@
 #include "gswteos-10.h"
 #include "gsw_check_data.h"
 
+
 #define test_func(name, arglist, value, var) \
         for (i=0; i<count; i++) { \
             value[i] = gsw_ ## name arglist; \
@@ -96,6 +97,33 @@ double          tf_poly[cast_m*cast_n];
 double          h[cast_m*cast_n];
 double          z[cast_m*cast_n];
 
+int assert_equal(int actual, int expected, const char *test_name) {
+    if (actual == expected) {
+        printf("%s: ............................... passed\n", test_name);
+    } else {
+        printf("%s: Failed (Expected: %d, Actual: %d)\n", test_name, expected, actual);
+        return gsw_error_flag=1;
+    }
+}
+
+
+void test_infunnel() {
+    int result;
+    // Inside the funnel.
+    result = gsw_infunnel(35.0, 5.0, 400);
+    assert_equal(result, 1, "gsw_infunnel -> inside funnel");
+
+    result = gsw_infunnel(20.0, 7.0, 2000);
+    assert_equal(result, 1, "gsw_infunnel -> inside funnel");
+
+    // Outside the funnel
+    result = gsw_infunnel(10.0, 0.0, 9000);
+    assert_equal(result, 0, "gsw_infunnel -> invalid pressure");
+
+    result = gsw_infunnel(45.0, 10.0, 400);
+    assert_equal(result, 0, "gsw_infunnel -> invalid salinity");
+}
+
 int
 main(int argc, char **argv)
 {
@@ -124,7 +152,6 @@ main(int argc, char **argv)
 " of TEOS-10 (version 3.05).\n");
 
         check_count = 1;
-
         section_title("Practical Salinity, PSS-78");
 
         test_func(c_from_sp, (sp[i],t[i],p[i]), c,c_from_sp);
@@ -556,6 +583,8 @@ main(int argc, char **argv)
         test_func(o2sol, (sa[i],ct[i],p[i],lon[i],lat[i]), value, o2sol);
         test_func(o2sol_sp_pt, (sp[i],pt[i]), value, o2sol_sp_pt);
 
+        test_infunnel();
+
         if (gsw_error_flag)
         {
             printf("\nYour installation of the Gibbs SeaWater (GSW) "

gsw_oceanographic_toolbox.c

@@ -180,6 +180,39 @@ gsw_add_mean(double *data_in, double *data_out)
 }
 /*
 !==========================================================================
+function gsw_infunnel(sa,ct,p)
+!==========================================================================
+
+! "oceanographic funnel" check for the 75-term equation
+!
+! sa     : Absolute Salinity                                 [g/kg]
+! ct     : Conservative Temperature                          [deg C]
+! p      : sea pressure                                      [dbar]
+!
+! gsw_infunnel : 0, if SA, CT and p are outside the "funnel"
+!                1, if SA, CT and p are inside the "funnel"
+!
+!  Note. The term "funnel" (McDougall et al., 2003) describes the range of
+!    SA, CT and p over which the error in the fit of the computationally
+!    efficient 75-term expression for specific volume in terms of SA, CT
+!    and p was calculated (Roquet et al., 2015).
+*/
+int
+gsw_infunnel(double sa, double ct, double p)
+{
+    return !(p > 8000 ||
+        sa < 0 ||
+        sa > 42 ||
+        (p < 500 && ct < gsw_ct_freezing(sa, p, 0)) ||
+        (p >= 500 && p < 6500 && sa < p * 5e-3 - 2.5) ||
+        (p >= 500 && p < 6500 && ct > (31.66666666666667 - p * 3.333333333333334e-3)) ||
+        (p >= 500 && ct < gsw_ct_freezing(sa, 500, 0)) ||
+        (p >= 6500 && sa < 30) ||
+        (p >= 6500 && ct > 10.0)
+    );
+}
+/*
+!==========================================================================
 function gsw_adiabatic_lapse_rate_from_ct(sa,ct,p)
 !==========================================================================
 
@@ -215,7 +248,7 @@ elemental function gsw_adiabatic_lapse_rate_ice (t, p)
 !         ( i.e. absolute pressure - 10.1325 dbar )
 !
 !    Note.  The output is in unit of degrees Celsius per Pa,
-!      (or equivilently K/Pa) not in units of K/dbar.
+!      (or equivalently K/Pa) not in units of K/dbar.
 !--------------------------------------------------------------------------
 */
 double
@@ -2781,12 +2814,12 @@ elemental subroutine gsw_frazil_properties (sa_bulk, h_bulk, p, &
 !  throughout the code).
 !
 !  When the mass fraction w_Ih_final is calculated as being a positive
-!  value, the seawater-ice mixture is at thermodynamic equlibrium.
+!  value, the seawater-ice mixture is at thermodynamic equilibrium.
 !
 !  This code returns w_Ih_final = 0 when the input bulk enthalpy, h_bulk,
 !  is sufficiently large (i.e. sufficiently "warm") so that there is no ice
 !  present in the final state.  In this case the final state consists of
-!  only seawater rather than being an equlibrium mixture of seawater and
+!  only seawater rather than being an equilibrium mixture of seawater and
 !  ice which occurs when w_Ih_final is positive.  Note that when
 !  w_Ih_final = 0, the final seawater is not at the freezing temperature.
 !
@@ -2978,12 +3011,12 @@ elemental subroutine gsw_frazil_properties_potential (sa_bulk, h_pot_bulk,&
 !  is assumed to be zero throughout the code).
 !
 !  When the mass fraction w_Ih_final is calculated as being a positive
-!  value, the seawater-ice mixture is at thermodynamic equlibrium.
+!  value, the seawater-ice mixture is at thermodynamic equilibrium.
 !
 !  This code returns w_Ih_final = 0 when the input bulk enthalpy, h_bulk,
 !  is sufficiently large (i.e. sufficiently "warm") so that there is no ice
 !  present in the final state.  In this case the final state consists of
-!  only seawater rather than being an equlibrium mixture of seawater and
+!  only seawater rather than being an equilibrium mixture of seawater and
 !  ice which occurs when w_Ih_final is positive.  Note that when
 !  w_Ih_final = 0, the final seawater is not at the freezing temperature.
 !
@@ -3249,12 +3282,12 @@ elemental subroutine gsw_frazil_properties_potential_poly (sa_bulk, &
 !  is assumed to be zero throughout the code).
 !
 !  When the mass fraction w_Ih_final is calculated as being a positive
-!  value, the seawater-ice mixture is at thermodynamic equlibrium.
+!  value, the seawater-ice mixture is at thermodynamic equilibrium.
 !
 !  This code returns w_Ih_final = 0 when the input bulk enthalpy, h_bulk,
 !  is sufficiently large (i.e. sufficiently "warm") so that there is no ice
 !  present in the final state.  In this case the final state consists of
-!  only seawater rather than being an equlibrium mixture of seawater and
+!  only seawater rather than being an equilibrium mixture of seawater and
 !  ice which occurs when w_Ih_final is positive.  Note that when
 !  w_Ih_final = 0, the final seawater is not at the freezing temperature.
 !
@@ -5803,12 +5836,12 @@ elemental subroutine gsw_melting_ice_into_seawater (sa, ct, p, w_ih, t_ih,&
 !  This code takes the seawater to contain no dissolved air.
 !
 !  When the mass fraction w_Ih_final is calculated as being a positive
-!  value, the seawater-ice mixture is at thermodynamic equlibrium.
+!  value, the seawater-ice mixture is at thermodynamic equilibrium.
 !
 !  This code returns w_Ih_final = 0 when the input bulk enthalpy, h_bulk,
 !  is sufficiently large (i.e. sufficiently "warm") so that there is no ice
 !  present in the final state.  In this case the final state consists of
-!  only seawater rather than being an equlibrium mixture of seawater and
+!  only seawater rather than being an equilibrium mixture of seawater and
 !  ice which occurs when w_Ih_final is positive.  Note that when
 !  w_Ih_final = 0, the final seawater is not at the freezing temperature.
 !

gswteos-10.h

@@ -37,6 +37,7 @@ extern void   gsw_add_barrier(double *input_data, double lon, double lat,
                 double long_grid, double lat_grid, double dlong_grid,
                 double dlat_grid, double *output_data);
 extern void   gsw_add_mean(double *data_in, double *data_out);
+extern int gsw_infunnel(double sa, double ct, double p);
 extern double gsw_adiabatic_lapse_rate_from_ct(double sa, double ct, double p);
 extern double gsw_adiabatic_lapse_rate_ice(double t, double p);
 extern double gsw_alpha(double sa, double ct, double p);