sound_speed_sea_unesco.m

function C = sound_speed_sea_unesco(T,S,P)
% Inputs
%   T: temperature \ degree Celsius \ 0 < T < 40 
%   S: salinity \ ppt \ 5 < S < 40
%   P: pressure \ kPa \ 0 < P < 100000 
% Outputs
%   C: speed of sound in seawater \ m/s

    C00=1402.388; C01=5.03830; C02=-5.81090e-2; C03=3.3432e-4; 
    C04=-1.47797e-6; C05=3.1419e-9;
    C10=0.153563; C11=6.8999e-4; C12=-8.1829e-6; C13=1.3632e-7;
    C14=-6.1260e-10; C20=3.1260e-5; C21=-1.7111e-6; 
    C22=2.5986e-8; C23=-2.5353e-10; C24=1.0415e-12; C30=-9.7729e-9;
    C31=3.8513e-10; C32=-2.3654e-12;
    A00=1.389; A01=-1.262e-2; A02=7.166e-5; A03=2.008e-6;
    A04=-3.21e-8;
    A10=9.4742e-5; A11=-1.2583e-5; A12=-6.4928e-8; A13=1.0515e-8;
    A14=-2.0142e-10;
    A20=-3.9064e-7; A21=9.1061e-9; A22=-1.6009e-10; A23=7.994e-12;
    A30=1.100e-10; A31=6.651e-12; A32=-3.391e-13;
    B00=-1.922e-2; B01=-4.42e-5; B10=7.3637e-5; B11=1.7950e-7;
    D00=1.727e-3;
    D10=-7.9836e-6;

    p = P*1e-2;

    Cw = (C00 + C01*T + C02*(T.^2) + C03*(T.^3) + C04*(T.^4) + C05*(T.^5)) ...
        + (C10 + C11*T + C12*(T.^2) + C13*(T.^3) + C14*(T.^4)).*p ...
        + (C20 + C21*T + C22*(T.^2) + C23*(T.^3) + C24*(T.^4)).*(p.^2) ...
        + (C30 + C31*T + C32*(T.^2)).*(p.^3);

    A = (A00 + A01*T + A02*(T.^2) + A03*(T.^3) + A04*(T.^4)) ...
        + (A10 + A11*T + A12*(T.^2) + A13*(T.^3) + A14*(T.^4)).*p ...
        + (A20 + A21*T + A22*(T.^2) + A23*(T.^3)).*(p.^2) ...
        + (A30 + A31*T + A32*(T.^2)).*(p.^3);

    B = B00 + B01*T + (B10 + B11*T).*p;

    D = D00+(D10.*p);

    C = Cw + A.*S + B.*(S.^(3/2)) + D.*(S.^2);
end