The literature on sound speed in water (freshwater and seawater) is rather confusing. There are 6 primary references in the literature (refs. [1]-[6] below). Treating them in chronological order, their content can be summarized as follows:

1. Greenspan and Tschiegg (1959): Tables of c vs T for distilled water, 0<T<100 deg C; no pressure is given and it is presumably standard atmospheric P. Also given is a ploynomial fit with standard deviation 0.026 m/sec or 17 ppm.
2. Wilson (1960a): Tables of c for seawater, -3 < T < 30 deg C; 1.033 (=1 atm)< P < 1000 kg/cm2; 33 < S < 37 ppt. Note that P is absolute pressure. A polynomial fit is reported with standard deviation 0.22 m/sec.
3. Wilson (1960b): Tables of sound speed in distilled water, for 0 < T < 100 degC, 14.7 < P < 14,000 psia (ie, 1~1000 atm). The accuracy is 0.1 ppt.
4. Wilson (1960c): A brief note giving an updated fit equation covering a wider range of parameters (especially salinity): -4 < T < 30 deg C, 1.0 < P < 1000 kg/cm2, 0 < S <37 ppt.
5. Del Grosso (1974): Provides a fit formula and tables of differences relative to results of previous formulas. Ranges are not given explicitly but tables in the paper cover 0 < T < 35 degC, 0 < S < 43 ppt, 0 < P < 1000 kg/cm2. Note that he uses gauge pressure, ie P=0 corresponds to 1 atm = 1.033 kg/cm2 absolute.
6. Chen and Millero (1977): Tables and fit for 0 < T < 40 deg C, 0 < P < 1000 bars, 5 < S < 40 ppt. Note pressure is gauge pressure (P=0 corresponds to 1 atm). Standard deviation of fit is 0.19 m/sec.

Attached is source code for Fortran-77 functions to compute the sound velocity for specified S,T,P using the fits described in refs. [1], [4] (which covers [2] and [3]), [5] and [6] within their valid ranges. The functions have been arranged for uniform input (T in degC, P in kg/cm2 absolute, S in ppt) and return c=-1 if input parameters are out of the range of validity for the fit coded. Naturally there is no warranty express or implied for the use of these functions!

Following the source code is a table of c vs T according to the various authors, for fresh water, and one point for seawater at T=1 degC, P=500 kg/cm2 abs, S=35 ppt.

[1] M. Greenspan and C. Tschiegg, (1959), JASA 31:75.

[2] W. Wilson, (1959), JASA 31:1067.

[3] W. Wilson, (1960a), JASA 32:641.

[4] W. Wilson, (1960b), JASA 32:1357.

[5] V. Del Grosso, (1974), JASA 56:1084.

[6] C. Chen and F. Millero, (1977), JASA 62:1129.

	double precision function greensp(T,P,S)
	implicit double precision (a-z)
c	calculate c(m/sec) in fresh water at 1 atm given t(degC)
c	according to m. greenspan and c. tschiegg, JASA 31:75 (1959)
	if (s.gt.0.or.p.gt.1.033) then
		greensp=-1.0
		return
		endif
	c=1402.736
	c=c + (5.03358)*t +(-0.0579506)*t**2
     #     + (3.31636e-04)*t**3 + (-1.45262e-06)*T**4
     #     + (3.0449e-09)*t**5
	greensp=c
	return
	end

	double precision function wilson(t,p,s)
c	find c(m/sec) in water with (T(degC), P(kg/cm2 abs), S(ppt))
c	according to wilson, JASA 1960
	implicit double precision (a-z)
	logical notok
c	test for in-range
	wilson=-1.0
	notok=.false.
	if ((t.lt.(-4.0)).or.(t.gt.30.0)) notok=.true.
	if ((s.lt.0).or.(s.gt.37.0)) notok=.true.
	if ((p.lt.0).or.(p.gt.1000.0)) notok=.true.
	if (notok) return
c	ok
	c0 = 1449.14
	ct1= 4.5721e00
	ct2=-4.4532e-02
	ct3=-2.6045e-04
	ct4= 7.9851e-06
	cp1=1.60272e-01
	cp2=1.0268e-05
	cp3=3.5216e-09
	cp4=-3.3603e-12
	dcs=1.39799*(S - 35.0)+1.69202e-03*(S - 35.0)**2
	cstp1=-1.1244e-02
	cstp2=7.7711e-07
	cstp3=7.7016e-05
	cstp4=-1.2943e-07
	cstp5=3.1580e-08
	cstp6= 1.5790e-09
	cstp7=-1.8607e-04
	cstp8=7.4812e-06
	cstp9=4.5283e-08
	cstp10=-2.5294e-07
	cstp11=1.8563e-09
	cstp12=-1.9646e-10
	dct=t*(ct1+t*(ct2+t*(ct3+t*ct4)))
	dcp=p*(cp1+p*(cp2+p*(cp3+p*cp4)))
	dcstp=(s-35.0)*(cstp1*t + cstp2*t**2 + cstp3*p
     #        + cstp4*p**2 + cstp5*p*t + cstp6*p*t**2)
     #        + p*(cstp7*t + cstp8*t**2 + cstp9*t**3)
     #        + p*p*(cstp10*t + cstp11*t**2)
     #        + p**3*cstp12*t
	c=c0+dct+dcp+dcs+dcstp
	wilson=c
	return
	end

	double precision function delgros(t,pa,s)
	implicit double precision (a-z)
	logical notok
c	returns c(m/s) in water with (T(deg C), P(kg/cm2 abs), S(ppt))
c	according to V. Del Grosso, JASA 56:1084 (1974)
c	convert kg/cm2 absolute to gauge pressure
	P=Pa-1.033
c	test for in-range
	delgros=-1.0
	notok=.false.
	if ((t.lt.0).or.(t.gt.35.0)) notok=.true.
	if ((s.lt.0).or.(s.gt.43.0)) notok=.true.
	if ((p.lt.0).or.(p.gt.1000.0)) notok=.true.
	if (notok) return
c	ok
	c0=1402.392
	ct1= 0.501109398873E+01
	ct2= -0.550946843172E-01
	ct3=+0.221535969240E-03
	cs1= 0.132952290781E+01
	cs2= +0.128955756844e-03
	cp1= 0.156059257041e00
	cp2= +0.244998688441e-04
	cp3= -0.883392332513e-08
	c1= -0.127562783426e-01
	c2= +0.635191613389e-02
	c3= +0.265484716608e-07
	c4= -0.159349479045e-05
	c5= +0.522116437235e-09
	c6= -0.438031096213e-06
	c7= -0.161674495909e-08
	c8= +0.968403156410e-04
	c9= +0.485639620015e-05
	c10= -0.340597039004e-03
	dct=t*(ct1 +t*(ct2 + t*ct3))
	dcs=s*(cs1 + s*cs2)
	dcp=p*(cp1 + p*(cp2 + p*cp3))
	dcstp=c1*t*s + c2*t*p + c3*(t*p)**2 + c4*t*p**2
     #       + c5*t*p**3 + c6*p*t**3 + c7*(s*p)**2
     #       + c8*s*t**2 + c9*t*p*s**2 +c10*t*s*p
	c=c0+dct+dcs+dcp+dcstp
	delgros=c
	return
	end


double precision chenmil(t,p0,s)	
c	* sound speed according to Chen and Millero (1977) JASA,62,1129
implicit none
double precision s,t,p0
double precision a,a0,a1,a2,a3,b,b0,b1
double c,c0,c1,c2,c3,p,sr,d,sv
c	convert from absolute to bars
	p = p0 - 1.033
c	test for in-range
	wilson=-1.0
	notok=.false.
	if ((t.lt.0.0).or.(t.gt.40.0)) notok=.true.
	if ((s.lt.5.0).or.(s.gt.40.0)) notok=.true.
	if ((p.lt.0).or.(p.gt.1000.0)) notok=.true.
	if (notok) return
c	ok
	sr = sqrt(s)
	d = 1.727e-3 - 7.9836e-6 * p
	b1 = 7.3637e-5 + 1.7945e-7 * t
	b0 = -1.922e-2 - 4.42e-5 * t
	b = b0 + b1 * p
	a3 = (-3.389e-13 * t + 6.649e-12) * t + 1.100e-10
	a2 = ((7.988e-12 * t - 1.6002e-10) * t 
     #    + 9.1041e-9) * t - 3.9064e-7
	a1 = (((-2.0122e-10 * t + 1.0507e-8) * t 
     #    - 6.4885e-8) * t - 1.2580e-5) * t + 9.4742e-5
	a0 = (((-3.21e-8 * t + 2.006e-6) * t 
     #     + 7.164e-5) * t -1.262e-2) * t + 1.389
	a = ((a3 * p + a2) * p + a1) * p + a0
	c3 = (-2.3643e-12 * t + 3.8504e-10) * t - 9.7729e-9
	c2 = (((1.0405e-12 * t -2.5335e-10) * t 
     #    + 2.5974e-8) * t - 1.7107e-6) * t + 3.1260e-5
	c1 = (((-6.1185e-10 * t + 1.3621e-7) * t 
     #   - 8.1788e-6) * t + 6.8982e-4) * t
		 + 0.153563
	c0 = ((((3.1464e-9 * t - 1.47800e-6) * t 
     #   + 3.3420e-4) * t - 5.80852e-2) * t
     #	 + 5.03711) * t + 1402.388
	c = ((c3 * p + c2) * p + c1) * p + c0
	chenmil = c + (a + b * sr + d * s) * s
	return
	end

Fresh water, surface
  T,degC        Greenspan       Wilson          DelGrosso
       0         0.14027E+04     0.14024E+04     0.14024E+04
       1         0.14077E+04     0.14074E+04     0.14073E+04
       2         0.14126E+04     0.14122E+04     0.14122E+04
       3         0.14173E+04     0.14169E+04     0.14169E+04
       4         0.14220E+04     0.14216E+04     0.14216E+04
       5         0.14265E+04     0.14261E+04     0.14261E+04
       6         0.14309E+04     0.14306E+04     0.14305E+04
       7         0.14352E+04     0.14350E+04     0.14348E+04
       8         0.14395E+04     0.14392E+04     0.14391E+04
       9         0.14436E+04     0.14434E+04     0.14432E+04
      10         0.14476E+04     0.14475E+04     0.14472E+04
      11         0.14515E+04     0.14514E+04     0.14511E+04
      12         0.14553E+04     0.14553E+04     0.14550E+04
      13         0.14591E+04     0.14591E+04     0.14587E+04
      14         0.14627E+04     0.14628E+04     0.14624E+04
      15         0.14662E+04     0.14664E+04     0.14659E+04
      16         0.14697E+04     0.14699E+04     0.14694E+04
      17         0.14731E+04     0.14734E+04     0.14727E+04
      18         0.14764E+04     0.14767E+04     0.14760E+04
      19         0.14795E+04     0.14800E+04     0.14792E+04
      20         0.14827E+04     0.14831E+04     0.14823E+04
Seawater, p=500, s=35
  T,degC        Greenspan       Wilson          DelGrosso
       1        -0.10000E+01     0.15364E+04     0.15358E+04