Home > wafo > misc > wdensity.m

# wdensity

## PURPOSE

Returns the water density

## SYNOPSIS

rho = wdensity(S,T,P)

## DESCRIPTION

``` WDENSITY  Returns the water density

CALL:  rho = wdensity(S,T,P);

rho =  water density [kg/m^3]
S   = salinity       [psu]       (Default 35)
T   = temperature    [degrees C] (default 4)
P   = pressure       [db]        (default 0)

Example:
S = linspace(20,40)'; T = linspace(4,20)';
[S1 T1]=meshgrid(S,T);
sc = contour(S,T,wdensity(S1,T1)); clabel(sc)
xlabel('Salinity'),ylabel('Temperature')```

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. clear Clear variables and functions from memory. error Display message and abort function.
This function is called by:
 specoptset Create or alter SPECTRUM OPTIONS structure. tran Computes transfer functions based on linear wave theory

## SOURCE CODE

```001 function rho = wdensity(S,T,P)
002 %WDENSITY  Returns the water density
003 %
004 %  CALL:  rho = wdensity(S,T,P);
005 %
006 %   rho =  water density [kg/m^3]
007 %   S   = salinity       [psu]       (Default 35)
008 %   T   = temperature    [degrees C] (default 4)
009 %   P   = pressure       [db]        (default 0)
010 %
011 % Example:
012 %  S = linspace(20,40)'; T = linspace(4,20)';
013 %  [S1 T1]=meshgrid(S,T);
014 %  sc = contour(S,T,wdensity(S1,T1)); clabel(sc)
015 %  xlabel('Salinity'),ylabel('Temperature')
016
017 % REFERENCES:
018 %  Fofonoff, P. and Millard, R.C. Jr (1983)
019 %  Algorithms for computation of fundamental properties of
020 %  seawater, Unesco Tech. Pap. in Mar. Sci., No. 44, 53, pp 17-18
021 %
022 %  Millero, F.J & Poisson, A. (1981)
023 %  International one-atmosphere equation of state for seawater.
024 %  Deep-Sea Research, Vol. 28A, No.6, pp 625-629.
025
026 %"An overview of the 1995 SWARM shallow water internal
027 %    acoustic scattering experiment", Apel, et el. IEEE J.
028 %    Ocean. Eng., Vol 22, No. 3, July 1997.
029
030 %    "A Surface-Trapped Intrusion of Slope Water onto the Continental
031 %    Shelf in the Middle Atlantic Bight", Gawarkiewicz, et. el. ,
032 %    Geophys. Res. Let. 23(25), 3763-3766, 1996.
033
034 %rho = 1026; Old value
035
036 error(nargchk(0,3,nargin))
037 if nargin<1|isempty(S), S = 35; end
038 if nargin<2|isempty(T), T = 4; end
039 if nargin<3|isempty(P), P = 0; end
040 [ec, S,T,P] = comnsize(S,T,P);
041 if ec>0,
042   error('S, T and P must be scalar or of common size')
043 end
044
045 % Density of fresh water at atmospheric pressure
046 %-------------------------------------------------
047 rho0 = 999.842594;
048 rho = rho0(ones(size(S)));
049
050 k = find(T~=0);
051 if any(k),
052   Tk = T(k);
053   a  = [6.536332e-9, -1.120083e-6, 1.001685e-4, -9.095290e-3, 6.793952e-2];
054   rho(k) = rho(k)+(a(5)+(a(4)+(a(3)+(a(2)+a(1)*Tk).*Tk).*Tk).*Tk).*Tk;
055   %rho(k) = rho(k) + polyval(a,Tk);
056 end
057
058 % Density of salt water at atmospheric pressure:
059 %------------------------------------------------
060 k1 = find(S~=0);
061 if any(k1),
062   c3 = -5.72466e-3;
063   b5 = 8.24493e-1;
064   tmp  = b5(ones(size(k1)));
065   tmp1 = c3(ones(size(k1)));
066   Tk   = T(k1);
067   k12  = find(Tk~=0);
068   if any(k12),
069     Tk = Tk(k12);
070     b = [ 5.3875e-9, -8.2467e-7,  7.6438e-5, -4.0899e-3, 0];
071     c = [ -1.6546e-6, 1.0227e-4, 0];
072     tmp(k12)  = tmp(k12)+(b(4) + (b(3) + (b(2) + b(1)*Tk).*Tk).*Tk).*Tk;
073     tmp1(k12) = tmp1(k12)+(c(2) + c(1)*Tk).*Tk;
074     clear Tk
075   end
076   d = 4.8314e-4;
077   Sk = S(k1);
078   rho(k1) = rho(k1) + (tmp + tmp1.*sqrt(Sk) + d*Sk).*Sk;
079 end
080
081 k2 = find(P~=0);
082 if any(k2),
083   Pk  = P(k2);
084   K   = seck(S(k2),T(k2),Pk);
085   Pk  = Pk/10;  % convert from db to atm pressure units
086   rho(k2) = rho(k2)./(1-Pk./K);
087 end
088 return
089
090
091 function K = seck(S,T,P)
092 %SECK  Secant bulk modulus (K) of sea water
093 %
094 % CALL:  K = seck(S,T,P);
095 %
096 %   S = salinity    [psu]
097 %   T = temperature [degrees C]
098 %   P = pressure    [db]
099 %
100 %   K = Secant Bulk Modulus  [bars]
101 %
102 %    Secant Bulk Modulus (K) of Sea Water using Equation of state 1980.
103 %    UNESCO polynomial implementation.
104
105
106 error(nargchk(3,3,nargin))
107
108
109 % Fresh water terms of K at atmosphere pressure.
110 %------------------------------------------------
111
112 K0 = 19652.21;
113 K  = K0(ones(size(S)));
114
115 k1 = find(T~=0);
116 if any(k1),
117   Tk = T(k1);
118   e = [-5.155288E-5, 1.360477E-2, -2.327105, 148.4206,0];
119   K(k1)  = K(k1) + (e(4) + (e(3) + (e(2) + e(1)*Tk).*Tk).*Tk).*Tk;   % eqn 19
120 end
121
122
123
124 % Salt water terms of K at atmosphere pressure
125 %---------------------------------------------
126
127 SR = sqrt(S);
128 k1 = find(S~=0);
129 if any(k1),
130   Tk = T(k1);
131   f = [-6.1670E-5, 1.09987E-2,-0.603459,54.6746];
132   g = [-5.3009E-4, 1.6483E-2,7.944E-2];
133
134   K(k1) = K(k1) + (  f(4) + (f(3) + (f(2) + f(1)*Tk).*Tk).*Tk ...
135       + (g(3) + (g(2) + g(1)*Tk).*Tk).*SR(k1)).*S(k1);      % eqn 16
136 end
137
138
139 k2 = find(P~=0);
140 if any(k2),
141   A0 = 3.239908; B0 = 8.50935E-5;
142   A  = A0(ones(size(k2)));
143   B  = B0(ones(size(k2)));
144   Tk = T(k2);
145
146   k12 = find(Tk~=0);
147   if any(k12)
148     Tk1 = Tk(k12);
149     k12 = k(k12);
150     h  = [-5.77905E-7, 1.16092E-4, 1.43713E-3, 0];
151     A(k12)  = A(k12) + (h(3) + (h(2) + h(1).*Tk1).*Tk1).*Tk1;
152
153     b = [ 5.2787E-8, -6.12293E-6,0];
154     B(k12)  = B(k12) + (b(2) + b(1)*Tk1).*Tk1;
155   end
156   J = 1.91075E-4;
157   I =[ -1.6078E-6, -1.0981E-5  2.2838E-3];
158   A  = A + (I(3) + (I(2) + I(1)*Tk).*Tk + J*SR(k2)).*S(k2);
159
160   m = [ 9.1697E-10, 2.0816E-8, -9.9348E-7];
161   B = B + (m(3) + (m(2) + m(1)*Tk).*Tk).*S(k2);
162
163   Pk = P(k2)/10;  %convert from db to atmospheric pressure units
164   K(k2) = K(k2) + (A + B.*Pk).*Pk;  % eqn 15
165 end
166 return
167
168
169
170
171
172
173
174
175
176```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003