*Fix bug in calculate_spec_vol_linear with spv_ref

  Corrected a sign error in calculate_spec_vol_array_linear and
calculate_spec_vol_scalar_linear when a reference specific volume is provided.
This bug will cause any configurations with EQN_OF_STATE="LINEAR" and
BOUSSINESQ=False (neither of which is the default value) to have the wrong sign
of the pressure gradients and other serious problems, like implausible sea
surface and internal interface heights.  This combination of parameters would
never be used in a realistic ocean model.  There are no impacted cases in any of
the MOM6-examples tests cases, nor those used in the ESMG or dev/NCAR test
suites, and it is very unlikely that any such case would work at all.  This bug
was present in the original version of the calculate_spec_vol_linear routines,
but it was only discovered after the implementation of the comprehensive
equation of state unit testing.  This will change answers in configurations that
could not have worked as viable ocean models, but answers are not impacted in
any known configuration, and all solutions in test cases are bitwise identical.

Author: Hallberg-NOAA

src/equation_of_state/MOM_EOS_linear.F90

@@ -119,7 +119,7 @@ subroutine calculate_spec_vol_scalar_linear(T, S, pressure, specvol, &
   real, optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
 
   if (present(spv_ref)) then
-    specvol = ((1.0 - Rho_T0_S0*spv_ref) + spv_ref*(dRho_dT*T + dRho_dS*S)) / &
+    specvol = ((1.0 - Rho_T0_S0*spv_ref) - spv_ref*(dRho_dT*T + dRho_dS*S)) / &
              ( Rho_T0_S0 + (dRho_dT*T + dRho_dS*S))
   else
     specvol = 1.0 / ( Rho_T0_S0 + (dRho_dT*T + dRho_dS*S))
@@ -148,7 +148,7 @@ subroutine calculate_spec_vol_array_linear(T, S, pressure, specvol, start, npts,
   integer :: j
 
   if (present(spv_ref)) then ; do j=start,start+npts-1
-    specvol(j) = ((1.0 - Rho_T0_S0*spv_ref) + spv_ref*(dRho_dT*T(j) + dRho_dS*S(j))) / &
+    specvol(j) = ((1.0 - Rho_T0_S0*spv_ref) - spv_ref*(dRho_dT*T(j) + dRho_dS*S(j))) / &
                  ( Rho_T0_S0 + (dRho_dT*T(j) + dRho_dS*S(j)))
   enddo ; else ; do j=start,start+npts-1
     specvol(j) = 1.0 / ( Rho_T0_S0 + (dRho_dT*T(j) + dRho_dS*S(j)))