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: This function is called by:

SUBFUNCTIONS ^

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