v_sa_p returned by gsw_specvol_second_derivatives is incorrect

[Hallberg-NOAA]

Self consistency testing within the MOM6 ocean model strongly suggests that the value of the second derivative of specific volume with respect to CT and P, as calculated in gsw_specvol_second_derivatives() is incorrect.

Specifically, I believe that line 109 of gsw_specvol_second_derivatives.F90, v_sa_p_part = ... + 4.0_r8*b004*z)) should be changed to v_sa_p_part = (... + 4.0_r8*b004*z)))/xs. The extra factor in the denominator of the derivative arises because xs is the square root of salinity with an offset. With this change the MOM6 self-consistency testing is indicating consistency between the various second derivatives of density returned from gsw_rho_second_derivatives(). For reference, compare this second derivative expression with the corresponding expression for the first derivative of specific volume with salinity, v_sa, on line 64, of gsw_specvol_first_derivatives.F90, which our testing indicates is correct, and which does include this extra xs factor in its denominator.

[efiring]

The corresponding code in GSW-C includes the division by xs (last line):

        if (v_sa_p != NULL) {

            v_sa_p_part = b001 + xs*(b101 + xs*(b201 + xs*(b301
                + b401*xs))) + ys*(b011 + xs*(b111 + xs*(b211
                + b311*xs)) + ys*(b021 + xs*(b121 + b221*xs)
                + ys*(b031 + b131*xs + b041*ys))) + z*(2.0*(b002 + xs*(b102
                + xs*(b202 + b302*xs)) + ys*(b012 + xs*(b112
                + b212*xs) + ys*(b022
                + b122*xs + b032*ys))) + z*(3.0*(b003 + b103*xs + b013*ys)
                + 4.0*b004*z));

            *v_sa_p = 1e-8*0.5*gsw_sfac*v_sa_p_part/xs;

        }

[efiring]

The fix to the C was ported from the Matlab, line 322 is this diff: https://github.com/TEOS-10/GSW-Matlab/commit/38c9635d6fd93e74c2648e4ee23cec49c1f58530

[ocefpaf]

Fixed in #27