Skip to content
Snippets Groups Projects
Commit fe678fff authored by oboebel's avatar oboebel
Browse files

Added seawater library to calculate depth from float max pressure for single source + topo tracking

parent 3eea3543
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 1883 additions and 0 deletions
lib/vendor/seawater_library/README.txt 0 → 100644
+ 8
− 0
View file @ fe678fff
% README file $Revision: 1.1 $ $Date: 2003/12/12 04:23:22 $
%
% SEAWATER is a toolkit of MATLAB routines for calculating the
% properties of sea water. They are a self contained library and
% are extremely easy to use.
%
% See the file sw_info.m for info on installation and usage
lib/vendor/seawater_library/sw_adtg.m 0 → 100644
+ 112
− 0
View file @ fe678fff
function ADTG = sw_adtg(S,T,P)
% SW_ADTG Adiabatic temperature gradient
%===========================================================================
% SW_ADTG $Id: sw_adtg.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1992.
%
% adtg = sw_adtg(S,T,P)
%
% DESCRIPTION:
% Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78) ]
% T = temperature [degree C (ITS-90)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
% ADTG = adiabatic temperature gradient [degree_C/db]
%
% AUTHOR: Phil Morgan, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Fofonoff, P. and Millard, R.C. Jr
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater. Unesco Tech. Pap. in Mar. Sci., No. 44, 53 pp. Eqn.(31) p.39
%
% Bryden, H. 1973.
% "New Polynomials for thermal expansion, adiabatic temperature gradient
% and potential temperature of sea water."
% DEEP-SEA RES., 1973, Vol20,401-408.
%=========================================================================
% Modifications
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
%-------------
% CHECK INPUTS
%-------------
if nargin ~= 3
error('sw_adtg.m: Must pass 3 parameters ')
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('check_stp: S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('check_stp: P has wrong dimensions')
end %if
%***check_stp
%-------------
% BEGIN
%-------------
T68 = 1.00024 * T;
a0 = 3.5803E-5;
a1 = +8.5258E-6;
a2 = -6.836E-8;
a3 = 6.6228E-10;
b0 = +1.8932E-6;
b1 = -4.2393E-8;
c0 = +1.8741E-8;
c1 = -6.7795E-10;
c2 = +8.733E-12;
c3 = -5.4481E-14;
d0 = -1.1351E-10;
d1 = 2.7759E-12;
e0 = -4.6206E-13;
e1 = +1.8676E-14;
e2 = -2.1687E-16;
ADTG = a0 + (a1 + (a2 + a3.*T68).*T68).*T68 ...
+ (b0 + b1.*T68).*(S-35) ...
+ ( (c0 + (c1 + (c2 + c3.*T68).*T68).*T68) ...
+ (d0 + d1.*T68).*(S-35) ).*P ...
+ ( e0 + (e1 + e2.*T68).*T68 ).*P.*P;
return
%==========================================================================
lib/vendor/seawater_library/sw_alpha.m 0 → 100644
+ 95
− 0
View file @ fe678fff
function [ALPHA] = sw_alpha(S, T, P, keyword)
% SW_ALPHA Thermal expansion coefficient (alpha)
%================================================================
% SW_ALPHA $Id: sw_alpha.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Nathan Bindoff 1993.
%
% USAGE: [ALPHA] = alpha(S, T, P, keyword)
%
% [ALPHA] = alpha(S, T, P, 'temp') %default
% [ALPHA] = alpha(S, PTMP, P, 'ptmp')
%
% DESCRIPTION:
% A function to calculate the thermal expansion coefficient.
%
% INPUT:
% S = salinity [psu (PSS-78) ]
% * PTMP = potential temperature [degree C (ITS-90)]
% * T = temperature [degree C (ITS-90)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% keyword = optional string to identify if temp or ptmp passed.
% = No argument defaults to 'temp'
% = 'temp' assumes (S,T,P) passed. Will execute slower
% as ptmp will be calculated internally.
% = 'ptmp' assumes (S,PTMP,P) passed. Will execute faster.
%
% OUTPUT:
% ALPHA = Thermal expansion coeff (alpha) [degree_C.^-1]
%
% AUTHOR: N.L. Bindoff 1993, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE:
% McDougall, T.J. 1987. "Neutral Surfaces"
% Journal of Physical Oceanography vol 17 pages 1950-1964,
%
% CHECK VALUE:
% See sw_beta.m amd sw_aonb.m
%================================================================
% Modifications
% 93-04-22. Phil Morgan, Help display modified to suit library
% 93-04-23. Phil Morgan, Input argument checking
% 94-10-15. Phil Morgan, Pass S,T,P and keyword for 'ptmp'
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CHECK INPUT ARGUMENTS
if ~(nargin==3 | nargin==4)
error('sw_alpha.m: requires 3 or 4 input arguments')
end %if
if nargin == 3
keyword = 'temp';
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('check_stp: S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('check_stp: P has wrong dimensions')
end %if
%***check_stp
% BEGIN
ALPHA = sw_aonb(S,T,P,keyword).*sw_beta(S,T,P,keyword);
return
%------------------------------------------------------------------------
lib/vendor/seawater_library/sw_aonb.m 0 → 100644
+ 121
− 0
View file @ fe678fff
function [AONB]= aonb(S, T, P, keyword)
% SW_AONB Calculate alpha/beta (a on b)
%================================================================
% SW_AONB $Revision: 1.1 $ $Date: 2003/12/12 04:23:22 $
% Copyright (C) CSIRO, Nathan Bindoff 1993
%
% USAGE: [AONB] = aonb(S, T, P, {keyword} )
%
% [AONB] = aonb(S, T, P, 'temp' ) %default
% [AONB] = aonb(S, PTMP, P, 'ptmp' )
%
% DESCRIPTION
% Calculate alpha/beta. See sw_alpha.m and sw_beta.m
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78) ]
% * PTMP = potential temperature [degree C (ITS-90)]
% * T = temperature [degree C (ITS-90)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% keyword = optional string to identify if temp or ptmp passed.
% = No argument defaults to 'temp'
% = 'temp' assumes (S,T,P) passed. Will execute slower
% as ptmp will be calculated internally.
% = 'ptmp' assumes (S,PTMP,P) passed. Will execute faster.
%
% OUTPUT
% AONB = alpha/beta [psu/degree_C]
%
% AUTHOR: N.L. Bindoff 1993, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE:
% McDougall, T.J. 1987. "Neutral Surfaces"
% Journal of Physical Oceanography vol 17 pages 1950-1964,
%
% CHECK VALUE:
% aonb=0.34763 psu C^-1 at S=40.0 psu, ptmp=10.0 C, p=4000 db
%================================================================
% Modifications
% 93-04-22. Phil Morgan, Help display modified to suit library
% 93-04-23. Phil Morgan, Input argument checking
% 94-10-15. Phil Morgan, Pass S,T,P and keyword for 'ptmp'
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CHECK INPUT ARGUMENTS
if ~(nargin==3 | nargin==4)
error('sw_aonb.m: requires 3 input arguments')
end %if
if nargin == 3
keyword = 'temp';
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('check_stp: S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('check_stp: P has wrong dimensions')
end %if
%***check_stp
% ENSURE WE USE PTMP IN CALCULATIONS
if ~strcmp(lower(keyword),'ptmp')
T = sw_ptmp(S,T,P,0); % now have ptmp
end %if
T = T * 1.00024;
% BEGIN
c1=fliplr([ 0.665157e-1, 0.170907e-1, ...
-0.203814e-3, 0.298357e-5, ...
-0.255019e-7]);
c2=fliplr([ 0.378110e-2, ...
-0.846960e-4]);
c2a=fliplr([0.0 -0.164759e-6, ...
-0.251520e-11]);
c3=[-0.678662e-5];
c4=fliplr([+0.380374e-4, -0.933746e-6, ...
+0.791325e-8]);
c5=[0.512857e-12];
c6=[-0.302285e-13];
%
% Now calaculate the thermal expansion saline contraction ratio adb
%
[m,n] = size(S);
sm35 = S-35.0*ones(m,n);
AONB = polyval(c1,T) + sm35.*(polyval(c2,T)...
+ polyval(c2a,P)) ...
+ sm35.^2*c3 + P.*polyval(c4,T) ...
+ c5*(P.^2).*(T.^2) + c6*P.^3;
return
%----------------------------------------------------------------------
lib/vendor/seawater_library/sw_beta.m 0 → 100644
+ 117
− 0
View file @ fe678fff
function [BETA] = sw_beta(S, T, P, keyword)
% SW_BETA Saline contraction coefficient (beta)
%========================================================================
% SW_BETA $Id: sw_beta.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% % Copyright (C) CSIRO, Nathan Bindoff 1993.
%
% USAGE: [BETA] = sw_beta(S, T, P, {keyword} )
%
% [BETA] = sw_beta(S, T, P, 'temp') %default
% [BETA] = sw_beta(S, PTMP, P, 'ptmp')
%
% DESCRIPTION
% The saline contraction coefficient as defined by T.J. McDougall.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78) ]
% * PTMP = potential temperature [degree C (ITS-90)]
% * T = temperature [degree C (ITS-90)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% keyword = optional string to identify if temp or ptmp passed.
% = No argument defaults to 'temp'
% = 'temp' assumes (S,T,P) passed. Will execute slower
% as ptmp will be calculated internally.
% = 'ptmp' assumes (S,PTMP,P) passed. Will execute faster.
%
% OUTPUT
% BETA = Saline Contraction Coefficient [psu.^-1]
%
% AUTHOR: N.L. Bindoff 1993, Lindsay Pender (Lindsay.pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE:
% McDougall, T.J. 1987. "Neutral Surfaces"
% Journal of Physical Oceanography vol 17 pages 1950-1964,
%
% CHECK VALUE:
% beta=0.72088e-3 psu.^-1 at S=40.0 psu, ptmp = 10.0 C (ITS-68), p=4000 db
%========================================================================
% Modifications
% 93-04-22. Phil Morgan, Help display modified to suit library
% 93-04-23. Phil Morgan, Input argument checking
% 94-10-15. Phil Morgan, Pass S,T,P and keyword for 'ptmp'
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CHECK INPUT ARGUMENTS
if ~(nargin==3 | nargin==4)
error('sw_beta.m: requires 3 or 4 input arguments')
end %if
if nargin == 3
keyword = 'temp';
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('check_stp: S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('check_stp: P has wrong dimensions')
end %if
%***check_stp
% ENSURE WE USE PTMP IN CALCULATIONS
if ~strcmp(lower(keyword),'ptmp')
T = sw_ptmp(S,T,P,0); % now have ptmp
end %if
T = T * 1.00024;
% BEGIN
c1=fliplr([ 0.785567e-3, -0.301985e-5 ...
0.555579e-7, -0.415613e-9]);
c2=fliplr([ -0.356603e-6, 0.788212e-8]);
c3=fliplr([0.0 0.408195e-10, -0.602281e-15]);
c4=[0.515032e-8];
c5=fliplr([-0.121555e-7, 0.192867e-9, -0.213127e-11]);
c6=fliplr([0.176621e-12 -0.175379e-14]);
c7=[0.121551e-17];
%
% Now calaculate the thermal expansion saline contraction ratio adb
%
[m,n] = size(S);
sm35 = S-35*ones(m,n);
BETA = polyval(c1,T) + sm35.*(polyval(c2,T) + ...
polyval(c3,P)) + c4*(sm35.^2) + ...
P.*polyval(c5,T) + (P.^2).*polyval(c6,T) ...
+c7*( P.^3);
return
%------------------------------------------------------------------------
lib/vendor/seawater_library/sw_bfrq.m 0 → 100644
+ 156
− 0
View file @ fe678fff
function [n2,q,p_ave] = sw_bfrq(S,T,P,LAT)
% SW_BFRQ Brunt-Vaisala Frequency Squared (N^2)
%===========================================================================
% SW_BFRQ $Id: sw_bfrq.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1993.
%
% USAGE: [bfrq,vort,p_ave] = sw_bfrq(S,T,P,{LAT})
%
% DESCRIPTION:
% Calculates Brunt-Vaisala Frequency squared (N^2) at the mid depths
% from the equation,
%
% -g d(pdens)
% N2 = ----- x --------
% pdens d(z)
%
% Also returns Potential Vorticity from q = f*N2/g.
%
% INPUT: (all must have same dimensions MxN)
% S = salinity [psu (PSS-78) ]
% T = temperature [degree C (ITS-90)]
% P = pressure [db]
%
% OPTIONAL:
% LAT = Latitude in decimal degrees north [-90..+90]
% May have dimensions 1x1 or 1xN where S(MxN).
% (Will use sw_g instead of the default g=9.8 m^2/s)
% (Will also calc d(z) instead of d(p) in numerator)
% OUTPUT:
% bfrq = Brunt-Vaisala Frequency squared (M-1xN) [s^-2]
% vort = Planetary Potential Vorticity (M-1xN) [(ms)^-1]
% (if isempty(LAT) vort=NaN )
% p_ave = Mid pressure between P grid (M-1xN) [db]
%
% AUTHOR: Phil Morgan 93-06-24, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% A.E. Gill 1982. p.54 eqn 3.7.15
% "Atmosphere-Ocean Dynamics"
% Academic Press: New York. ISBN: 0-12-283522-0
%
% Jackett, D.R. and McDougall, T.J. 1994.
% Minimal adjustment of hydrographic properties to achieve static
% stability. submitted J.Atmos.Ocean.Tech.
%
% Greg Johnson (gjohnson@pmel.noaa.gov)
% added potential vorticity calcuation
%=========================================================================
% Modifications
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CALLER: general purpose
% CALLEE: sw_dens.m sw_pden.m
%$Id: sw_bfrq.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
%-------------
% Check Inputs
%-------------
error(nargchk(3, 4, nargin));
if nargin == 3
LAT = [];
end
% Get S,T,P dimensions and check consistency
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% Check S, T P and have the same shape
if (ms~=mt) | (ns~=nt) | (np~=np)
error('S, T & P must have same dimensions')
end
% Check S and T have length at least of 2
if (ms * ns == 1)
error('Length of T, S and P must be at least 2')
end
% If S, T and P are row vectors - transpose
if ms == 1
S = S';
T = T';
P = P';
transpose = 1;
else
transpose = 0;
end
% If lat passed then verify dimensions
if ~isempty(LAT)
[mL,nL] = size(LAT);
if mL==1 & nL==1
LAT = LAT*ones(size(S));
else
if (ms~=mL) | (ns~=nL) % S & LAT are not the same shape
if (ns==nL) & (mL==1) % copy LATS down each column
LAT = LAT( ones(1,ms), : ); % s.t. dim(S)==dim(LAT)
else
error('Inputs arguments have wrong dimensions')
end
end
end
end
%------
% Begin
%------
if ~isempty(LAT)
% note that sw_g expects height as argument
Z = sw_dpth(P,LAT);
g = sw_g(LAT,-Z);
f = sw_f(LAT);
else
Z = P;
g = 9.8*ones(size(P));
f = NaN*ones(size(P));
end %if
[m,n] = size(P);
iup = 1:m-1;
ilo = 2:m;
p_ave = (P(iup,:)+P(ilo,:) )/2;
pden_up = sw_pden(S(iup,:),T(iup,:),P(iup,:),p_ave);
pden_lo = sw_pden(S(ilo,:),T(ilo,:),P(ilo,:),p_ave);
mid_pden = (pden_up + pden_lo )/2;
dif_pden = pden_up - pden_lo;
mid_g = (g(iup,:)+g(ilo,:))/2;
dif_z = diff(Z);
n2 = -mid_g .* dif_pden ./ (dif_z .* mid_pden);
mid_f = f(iup,:);
q = mid_f .* dif_pden ./ (dif_z .* mid_pden);
if transpose
n2 = n2';
q = q';
p_ave = p_ave';
end
lib/vendor/seawater_library/sw_c3515.m 0 → 100644
+ 40
− 0
View file @ fe678fff
function c3515 = sw_c3515()
% SW_C3515 Conductivity at (35,15,0)
%=========================================================================
% SW_c3515 $Id: sw_c3515.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% % Copyright (C) CSIRO, Phil Morgan 1993.
%
% USAGE: c3515 = sw_c3515
%
% DESCRIPTION:
% Returns conductivity at S=35 psu , T=15 C [ITPS 68] and P=0 db).
%
% INPUT: (none)
%
% OUTPUT:
% c3515 = Conductivity [mmho/cm == mS/cm]
%
% AUTHOR: Phil Morgan 93-04-17 (morgan@ml.csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% R.C. Millard and K. Yang 1992.
% "CTD Calibration and Processing Methods used by Woods Hole
% Oceanographic Institution" Draft April 14, 1992
% (Personal communication)
%=========================================================================
% CALLER: none
% CALLEE: none
%
c3515 = 42.914;
return
%-------------------------------------------------------------------------
lib/vendor/seawater_library/sw_cndr.m 0 → 100644
+ 135
− 0
View file @ fe678fff
function R = sw_cndr(S,T,P)
% SW_CNDR Conductivity ratio R = C(S,T,P)/C(35,15(IPTS-68),0)
%=========================================================================
% SW_CNDR $Id: sw_cndr.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1993.
%
% USAGE: cndr = sw_cndr(S,T,P)
%
% DESCRIPTION:
% Calculates conductivity ratio from S,T,P.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78) ]
% T = temperature [degree C (ITS-90)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
% cndr = Conductivity ratio R = C(S,T,P)/C(35,15(IPTS-68),0) [no units]
%
% AUTHOR: Phil Morgan 93-04-21, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Fofonoff, P. and Millard, R.C. Jr
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
%=========================================================================
% Modifications
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CALLER: general purpose
% CALLEE: sw_salds.m sw_sals.m sw_salrt.m
%--------------
% check inputs
%-------------
if nargin~=3
error('sw_cndr.m: must have 3 input arguments')
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('check_stp: S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('check_stp: P has wrong dimensions')
end %if
%***check_stp
%-------
% BEGIN
%-------
T68 = T * 1.00024;
for i = 1:ms
for j = 1:ns
%---------------------------------------------------------------------
% DO A NEWTON-RAPHSON ITERATION FOR INVERSE INTERPOLATION OF Rt FROM S.
%---------------------------------------------------------------------
S_loop = S(i,j); % S in the loop
T_loop = T(i,j); % T in the loop
Rx_loop = sqrt(S_loop/35.0); % first guess at Rx = sqrt(Rt)
SInc = sw_sals(Rx_loop.*Rx_loop,T_loop); % S INCrement (guess) from Rx
iloop = 0;
end_loop = 0;
while ~end_loop
Rx_loop = Rx_loop + (S_loop - SInc)./sw_salds(Rx_loop,T_loop - 15);
SInc = sw_sals(Rx_loop.*Rx_loop,T_loop);
iloop = iloop + 1;
dels = abs(SInc-S_loop);
if (dels>1.0e-4 & iloop<10)
end_loop = 0;
else
end_loop = 1;
end %if
end %while
Rx(i,j) = Rx_loop;
end %for j
end %for i
%------------------------------------------------------
% ONCE Rt FOUND, CORRESPONDING TO EACH (S,T) EVALUATE R
%------------------------------------------------------
% eqn(4) p.8 Unesco 1983
d1 = 3.426e-2;
d2 = 4.464e-4;
d3 = 4.215e-1;
d4 = -3.107e-3;
e1 = 2.070e-5;
e2 = -6.370e-10;
e3 = 3.989e-15;
A = (d3 + d4.*T68);
B = 1 + d1.*T68 + d2.*T68.^2;
C = P.*(e1 + e2.*P + e3.*P.^2);
% eqn(6) p.9 UNESCO 1983.
Rt = Rx.*Rx;
rt = sw_salrt(T);
Rtrt = rt.*Rt;
D = B - A.*rt.*Rt;
E = rt.*Rt.*A.*(B+C);
R = sqrt(abs(D.^2+4*E)) - D;
R = 0.5*R./A;
lib/vendor/seawater_library/sw_copy.m 0 → 100644
+ 164
− 0
View file @ fe678fff
% SW_COPY Copyright and licence information on SEAWATER library.
% =================================================================
% SW_COPY $Id: sw_copy.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) Phil Morgan 1993
%
% SOFTWARE LICENCE AGREEMENT
%
% 1.0 Grant of Licence
%
% 1.1 The CSIRO Division of Oceanography (herein referred to as
% "CSIRO") hereby grants you (hereinafter referred to as
% the "Licencee"), subject to the Licencee agreeing to
% comply with the terms and conditions of this Agreement, a
% non-transferable, non-exclusive licence to use the
% computer programs described in this document (hereinafter
% referred to as the "Software") for the purpose of
% the Licencee's computing activity.
%
% 1.2 CSIRO hereby grants the Licencee the right to make copies
% of the Software for the purpose of the Licencee's
% computing activity only.
%
% 1.3 The benefit of the rights granted to the Licencee by the
% Licence and this Agreement generally shall be personal to
% the Licencee and the Licencee shall not mortgage, charge,
% assign, rent, lease, sell or otherwise dispose of or
% transfer the same or any part to any third party.
%
% 1.4 Unless otherwise agreed in writing or provided for in
% this Agreement, CSIRO shall be under no obligation or
% responsibility to provide the Licencee with any training,
% maintenance services, enhancements or updates of the
% Software or any services whatsoever.
%
% 2.0 Acknowledgment by the Licencee
%
% 2.1 The Licencee acknowledges and agrees that it shall not:
%
% (i) sell, let for hire or by way of trade, offer or
% exhibit or expose for sale or hire or otherwise
% distribute the Software for the purposes of trade or
% any other purpose;
%
% (ii) authorise or assist any third person to do any
% of the acts set out in (i) above;
%
% (iii) modify the Software source code without advising
% CSIRO.
%
% 2.2 The Licencee agrees that:
%
% (a) CSIRO is the owner of all copyright and other
% Intellectual Property Rights subsisting in the
% Software;
%
% (b) this document must be properly cited in any
% publication reporting results derived from this
% document or obtained from application and use of this
% software. Any of the Licencee's documentation
% describing results generated by the Licencee's
% use of the Software will contain an acknowledgement
% of CSIRO's ownership of the Software;
%
% (c) CSIRO reserves all rights in the Software other than
% the rights granted to the Licencee by this
% Agreement;
%
% (d) each item of the Software will display a banner
% summarising the terms of this Agreement and
% acknowledging the source of the Software, and the
% contents of a banner will not be modified and its
% display will not be inactivated by the Licencee
% without the approval of CSIRO.
%
% 3.0 Indemnity
%
% 3.1 To the full extent permitted by law, CSIRO excludes any
% and all liability in respect of any loss or damage,
% whether personal (includes death or illness) or of
% property and whether direct, consequential or special
% (including consequential financial loss or damage) of
% the Licencee, its officers, agents and employees or any
% third party howsoever caused, which may be suffered or
% incurred or which may arise directly or indirectly in
% respect of or arising out of the Licencee's use or
% inability to use the Software or the failure or omission
% on the part of CSIRO to comply with the conditions and
% warranties under this Licence Agreement. Insofar as
% liability for loss or damages under or pursuant to such
% legislation cannot be excluded, CSIRO's liability for
% loss or damages shall be limited to the amount of One
% Dollar ($1.00).
%
% 3.2 CSIRO make no warranties, expressed or implied, and
% excludes all other warranties representations, terms or
% conditions, whether express or implied, oral or written,
% statutory or otherwise, relating in any way to the
% Software, or to this Agreement, including any implied
% warranty of merchantability or of fitness for particular
% purpose. To the full extent permitted by the law of the
% Commonwealth of Australia or the laws of any State or
% Territory of Australia, any conditions or warranties
% imposed by such legislation are hereby excluded. In so
% far as liability under or pursuant to such legislation
% may not be excluded, CSIRO's liability to the Licencee
% pursuant to this Agreement shall be limited as set out in
% clause 3.1 hereof.
%
% 3.3 The Licencee acknowledges and agrees that the Software
% was developed for CSIRO research purposes and may have
% inherent defects, errors or deficiencies, and that it is
% the responsibility of the Licencee to make its own
% assessment of the suitability of the Software for the
% purpose of the Licencee's computing activity. The
% Licencee will use the Software, and advice, opinions or
% information supplied by CSIRO, its officers, employees or
% agents concerning the Software at the Licencee's own
% risk.
%
% 3.4 The Licencee hereby releases and indemnifies and shall
% continue to release and indemnify CSIRO, its officers,
% employees and agents from and against all actions,
% claims, proceedings or demands (including those brought
% by third parties) which may be bought against it or them,
% whether on their own or jointly with the Licencee and
% whether at common law, in equity or pursuant to statute or
% otherwise, in respect of any loss, death, injury, illness
% or damage (whether personal or property, and whether
% direct or consequential, including consequential
% financial loss) and any infringement of copyright,
% patents, trade marks, designs or other Intellectual
% Property Rights, howsoever arising out of the Licencee's
% exercise of its rights under this Agreement and from and
% against all damages, costs and expenses incurred in
% defending or settling any such claim, proceeding or
% demand.
%
% 3.5 The Licencee's obligation to indemnify CSIRO and its
% officers, employees and agents set out in clause 3.4
% hereof is a continuing obligation separate from
% and independent of the Licencee's other obligations under
% this Agreement, and shall survive all expiration or
% termination of this Agreement.
%
% 4.0 Termination
%
% 4.1 The Licence shall terminate immediately upon the Licencee
% breaching any term or condition of this Agreement whether
% or not CSIRO is aware of the occurrence of the breach at
% the time that it happens.
%
% 4.2 CSIRO may terminate the Licence on reasonable grounds by
% notice in writing to the Licencee, and such notice of
% termination shall be effective immediately upon receipt
% by the Licencee.
%
%=========================================================================
more on
help sw_copy
more off
return
%--------------------------------------------------------------------
lib/vendor/seawater_library/sw_cp.m 0 → 100644
+ 162
− 0
View file @ fe678fff
function cp = sw_cp(S,T,P)
% SW_CP Heat Capacity (Cp) of sea water
%=========================================================================
% SW_CP $Id: sw_cp.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1993.
%
% USAGE: cp = sw_cp(S,T,P)
%
% DESCRIPTION:
% Heat Capacity of Sea Water using UNESCO 1983 polynomial.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78)]
% T = temperature [degree C (ITS-90)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
% cp = Specific Heat Capacity [J kg^-1 C^-1]
%
% AUTHOR: Phil Morgan, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Fofonff, P. and Millard, R.C. Jr
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
%=========================================================================
% Modifications
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CALLER: general purpose
% CALLEE: none
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin ~=3
error('Must pass 3 parameters')
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('P has wrong dimensions')
end %if
%***check_stp
%------
% BEGIN
%------
P = P/10; % to convert db to Bar as used in Unesco routines
T68 = T * 1.00024;
%------------
% eqn 26 p.32
%------------
c0 = 4217.4;
c1 = -3.720283;
c2 = 0.1412855;
c3 = -2.654387e-3;
c4 = 2.093236e-5;
a0 = -7.64357;
a1 = 0.1072763;
a2 = -1.38385e-3;
b0 = 0.1770383;
b1 = -4.07718e-3;
b2 = 5.148e-5;
Cpst0 = (((c4.*T68 + c3).*T68 + c2).*T68 + c1).*T68 + c0 + ...
(a0 + a1.*T68 + a2.*T68.^2).*S + ...
(b0 + b1.*T68 + b2.*T68.^2).*S.*sqrt(S);
%------------
% eqn 28 p.33
%------------
a0 = -4.9592e-1;
a1 = 1.45747e-2;
a2 = -3.13885e-4;
a3 = 2.0357e-6;
a4 = 1.7168e-8;
b0 = 2.4931e-4;
b1 = -1.08645e-5;
b2 = 2.87533e-7;
b3 = -4.0027e-9;
b4 = 2.2956e-11;
c0 = -5.422e-8;
c1 = 2.6380e-9;
c2 = -6.5637e-11;
c3 = 6.136e-13;
del_Cp0t0 = (((((c3.*T68 + c2).*T68 + c1).*T68 + c0).*P + ...
((((b4.*T68 + b3).*T68 + b2).*T68 + b1).*T68 + b0)).*P + ...
((((a4.*T68 + a3).*T68 + a2).*T68 + a1).*T68 + a0)).*P;
%------------
% eqn 29 p.34
%------------
d0 = 4.9247e-3;
d1 = -1.28315e-4;
d2 = 9.802e-7;
d3 = 2.5941e-8;
d4 = -2.9179e-10;
e0 = -1.2331e-4;
e1 = -1.517e-6;
e2 = 3.122e-8;
f0 = -2.9558e-6;
f1 = 1.17054e-7;
f2 = -2.3905e-9;
f3 = 1.8448e-11;
g0 = 9.971e-8;
h0 = 5.540e-10;
h1 = -1.7682e-11;
h2 = 3.513e-13;
j1 = -1.4300e-12;
S3_2 = S.*sqrt(S);
del_Cpstp = [((((d4.*T68 + d3).*T68 + d2).*T68 + d1).*T68 + d0).*S + ...
((e2.*T68 + e1).*T68 + e0).*S3_2].*P + ...
[(((f3.*T68 + f2).*T68 + f1).*T68 + f0).*S + ...
g0.*S3_2].*P.^2 + ...
[((h2.*T68 + h1).*T68 + h0).*S + ...
j1.*T68.*S3_2].*P.^3;
cp = Cpst0 + del_Cp0t0 + del_Cpstp;
lib/vendor/seawater_library/sw_dens.m 0 → 100644
+ 91
− 0
View file @ fe678fff
function dens = sw_dens(S,T,P)
% SW_DENS Density of sea water
%=========================================================================
% SW_DENS $Id: sw_dens.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1992.
%
% USAGE: dens = sw_dens(S,T,P)
%
% DESCRIPTION:
% Density of Sea Water using UNESCO 1983 (EOS 80) polynomial.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78)]
% T = temperature [degree C (ITS-90)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
% dens = density [kg/m^3]
%
% AUTHOR: Phil Morgan 92-11-05, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Fofonoff, P. and Millard, R.C. Jr
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
%
% Millero, F.J., Chen, C.T., Bradshaw, A., and Schleicher, K.
% " A new high pressure equation of state for seawater"
% Deap-Sea Research., 1980, Vol27A, pp255-264.
%=========================================================================
% Modifications
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CALLER: general purpose
% CALLEE: sw_dens0.m sw_seck.m
% UNESCO 1983. eqn.7 p.15
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin ~=3
error('sw_dens.m: Must pass 3 parameters')
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('check_stp: S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('check_stp: P has wrong dimensions')
end %if
%***check_stp
%------
% BEGIN
%------
densP0 = sw_dens0(S,T);
K = sw_seck(S,T,P);
P = P/10; % convert from db to atm pressure units
dens = densP0./(1-P./K);
return
%--------------------------------------------------------------------
lib/vendor/seawater_library/sw_dens0.m 0 → 100644
+ 82
− 0
View file @ fe678fff
function dens = sw_dens0(S,T)
% SW_DENS0 Denisty of sea water at atmospheric pressure
%=========================================================================
% SW_DENS0 $Id: sw_dens0.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1992
%
% USAGE: dens0 = sw_dens0(S,T)
%
% DESCRIPTION:
% Density of Sea Water at atmospheric pressure using
% UNESCO 1983 (EOS 1980) polynomial.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78)]
% T = temperature [degree C (ITS-90)]
%
% OUTPUT:
% dens0 = density [kg/m^3] of salt water with properties S,T,
% P=0 (0 db gauge pressure)
%
% AUTHOR: Phil Morgan 92-11-05, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
%
% Millero, F.J. and Poisson, A.
% International one-atmosphere equation of state of seawater.
% Deep-Sea Res. 1981. Vol28A(6) pp625-629.
%=========================================================================
% Modifications
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CALLER: general purpose, sw_dens.m
% CALLEE: sw_smow.m
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin ~=2
error('sw_dens0.m: Must pass 2 parameters')
end %if
[mS,nS] = size(S);
[mT,nT] = size(T);
if (mS~=mT) | (nS~=nT)
error('sw_dens0.m: S,T inputs must have the same dimensions')
end %if
%----------------------
% DEFINE CONSTANTS
%----------------------
T68 = T * 1.00024;
% UNESCO 1983 eqn(13) p17.
b0 = 8.24493e-1;
b1 = -4.0899e-3;
b2 = 7.6438e-5;
b3 = -8.2467e-7;
b4 = 5.3875e-9;
c0 = -5.72466e-3;
c1 = +1.0227e-4;
c2 = -1.6546e-6;
d0 = 4.8314e-4;
dens = sw_smow(T) + (b0 + (b1 + (b2 + (b3 + b4*T68).*T68).*T68).*T68).*S ...
+ (c0 + (c1 + c2*T68).*T68).*S.*sqrt(S) + d0*S.^2;
return
%--------------------------------------------------------------------
lib/vendor/seawater_library/sw_dist.m 0 → 100644
+ 104
− 0
View file @ fe678fff
function [dist,phaseangle] = sw_dist(lat,lon,units)
% SW_DIST Distance between two lat,lon coordinates
%===================================================================
% SW_DIST $Id: sw_dist.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan & Steve Rintoul 1992.
%
% USAGE: [dist,phaseangle] = sw_dist(lat,lon {,units} )
%
% DESCRIPTION:
% Calculate distance between two positions on glode using the "Plane
% Sailing" method. Also uses simple geometry to calculate the bearing of
% the path between position pairs.
%
% INPUT:
% lat = decimal degrees (+ve N, -ve S) [- 90.. +90]
% lon = decimal degrees (+ve E, -ve W) [-180..+180]
% units = optional string specifing units of distance
% 'nm' = nautical miles (default)
% 'km' = kilometres
%
% OUTPUT:
% dist = distance between positions in units
% phaseangle = angle of line between stations with x axis (East).
% Range of values are -180..+180. (E=0, N=90, S=-90)
%
% AUTHOR: Phil Morgan and Steve Rintoul 92-02-10
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE:
% The PLANE SAILING method as descriibed in "CELESTIAL NAVIGATION" 1989 by
% Dr. P. Gormley. The Australian Antartic Division.
%==================================================================
% Modifications
% 99-06-25. Lindsay Pender, Function name change from distance to sw_dist.
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% CALLER: general purpose
% CALLEE: none
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin > 3
error('sw_dist.m: No more than 3 arguments allowed')
elseif nargin==3
if ~isstr(units)
error('sw_dist.m: units argument must be string')
end %if
elseif nargin==2
units = 'nm'; % default units
else
error('sw_dist.m: wrong number of arguments')
end%if
[mlat,nlat] = size(lat);
if mlat~=1 & nlat~=1
error('sw_dist.m: lat, lon must be vectors. No matrices allowed')
end%if
if length(lat)~=length(lon)
error('sw_dist.m: lat and lon must have same number of elements')
end%if
%-----------------
% DEFINE CONSTANTS
%-----------------
DEG2RAD = (2*pi/360);
RAD2DEG = 1/DEG2RAD;
DEG2MIN = 60;
DEG2NM = 60;
NM2KM = 1.8520; % Defined in Pond & Pickard p303.
% BEGIN
npositions = length(lat);
ind=1:npositions-1; % index to first of position pairs
dlon = diff(lon);
if any(abs(dlon)>180)
flag = find(abs(dlon)>180);
for ii=1:length(flag)
dlon(flag(ii))= -sign(dlon(flag(ii))) * (360 - abs(dlon(flag(ii))) );
end %for
end %if
latrad = abs(lat*DEG2RAD);
dep = cos( (latrad(ind+1)+latrad(ind))./2 ) .* dlon;
dlat = diff(lat);
dist = DEG2NM*sqrt(dlat.^2 + dep.^2); % in n.miles
if strcmp(units,'km') % defaults to n.miles
dist = dist * NM2KM;
end %if
% CALCUALTE ANGLE TO X AXIS
phaseangle = angle(dep+dlat*sqrt(-1))*RAD2DEG;
return
%--------------------------------------------------------------------
lib/vendor/seawater_library/sw_dpth.m 0 → 100644
+ 81
− 0
View file @ fe678fff
function DEPTHM = sw_dpth(P,LAT)
% SW_DPTH Depth from pressure
%===========================================================================
% SW_DPTH $Id: sw_dpth.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1992.
%
% USAGE: dpth = sw_dpth(P,LAT)
%
% DESCRIPTION:
% Calculates depth in metres from pressure in dbars.
%
% INPUT: (all must have same dimensions)
% P = Pressure [db]
% LAT = Latitude in decimal degress north [-90..+90]
% (lat may have dimensions 1x1 or 1xn where P(mxn).
%
% OUTPUT:
% dpth = depth [metres]
%
% AUTHOR: Phil Morgan 92-04-06 (morgan@ml.csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
%=========================================================================
% Modifications
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% CALLER: general purpose
% CALLEE: none
%-------------
% CHECK INPUTS
%-------------
[mP,nP] = size(P);
[mL,nL] = size(LAT);
if mL==1 & nL==1 % LAT scalar - fill to size of P
LAT = LAT*ones(size(P));
elseif nP == nL & mL == 1 % LAT is row vector
LAT = LAT(ones(1, mP), :); % Coppy down each column
elseif mP == mL & nL == 1 % LAT is column vector
LAT = LAT(:, ones(1, nP)); % Copy across each row
elseif mP == mL & nP == nL
% Ok
else
error('sw_depth.m: Inputs arguments have wrong dimensions')
end %if
%-------------
% BEGIN
%-------------
% Eqn 25, p26. Unesco 1983.
DEG2RAD = pi/180;
c1 = +9.72659;
c2 = -2.2512E-5;
c3 = +2.279E-10;
c4 = -1.82E-15;
gam_dash = 2.184e-6;
LAT = abs(LAT);
X = sin(LAT*DEG2RAD); % convert to radians
X = X.*X;
bot_line = 9.780318*(1.0+(5.2788E-3+2.36E-5*X).*X) + gam_dash*0.5*P;
top_line = (((c4*P+c3).*P+c2).*P+c1).*P;
DEPTHM = top_line./bot_line;
return
%===========================================================================
%
lib/vendor/seawater_library/sw_f.m 0 → 100644
+ 58
− 0
View file @ fe678fff
function f = sw_f(lat)
% SW_F Coriolis factor "f"
%===========================================================================
% SW_F $Id: sw_f.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1993.
%
% USAGE: f = sw_f(lat)
%
% DESCRIPTION:
% Calculates the Coriolis factor "f" defined by
% f = 2*Omega*Sin(lat) where Omega = 7.292e-5 radians/sec
%
% INPUT:
% lat = Latitude in decimal degress north [-90..+90]
%
% OUTPUT:
% f = Coriolis Factor "f" [s-1]
%
% AUTHOR: Phil Morgan 93-04-20 (morgan@ml.csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE:
% S. Pond & G.Pickard 2nd Edition 1986
% Introductory Dynamical Oceanogrpahy
% Pergamon Press Sydney. ISBN 0-08-028728-X
%
% A.E. Gill 1982. p.597
% "Atmosphere-Ocean Dynamics"
% Academic Press: New York. ISBN: 0-12-283522-0
%=========================================================================
% CALLER: general purpose
% CALLEE: none
%-------------
% CHECK INPUTS
%-------------
if nargin ~= 1
error('sw_f.m: Requires one input argument')
end %if
%-------------
% BEGIN
%-------------
% Eqn p27. Unesco 1983.
DEG2RAD = pi/180;
OMEGA = 7.292e-5; %s-1 A.E.Gill p.597
f = 2*OMEGA*sin(lat*DEG2RAD);
return
%===========================================================================
lib/vendor/seawater_library/sw_fp.m 0 → 100644
+ 78
− 0
View file @ fe678fff
function fp = sw_fp(S,P)
% SW_FP Freezing point of sea water
%=========================================================================
% SW_FP % $Id: sw_fp.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1993.
%
% USAGE: fp = sw_fp(S,P)
%
% DESCRIPTION:
% Heat Capacity of Sea Water using UNESCO 1983 polynomial.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78)]
% P = pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
% fp = Freezing Point temperature [degree C (ITS-90)]
%
% AUTHOR: Phil Morgan 93-04-20, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Fofonff, P. and Millard, R.C. Jr
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
%=========================================================================
% Modifications
% 99-06-25. Lindsay Pender, Fixed transpose of row vectors.
% 03-12-12. Lindsay Pender, Converted to ITS-90.
% CALLER: general purpose
% CALLEE: none
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin ~=2
error('Must pass 3 parameters')
end %if
[ms,ns] = size(S);
[mp,np] = size(P);
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('P has wrong dimensions')
end %if
%------
% BEGIN
%------
%P = P/10; % to convert db to Bar as used in Unesco routines
%------------
% eqn p.29
%------------
a0 = -0.0575;
a1 = 1.710523e-3;
a2 = -2.154996e-4;
b = -7.53e-4;
fp = (a0.*S + a1.*S.*sqrt(S) + a2.*S.^2 + b.*P) / 1.00024;
lib/vendor/seawater_library/sw_g.m 0 → 100644
+ 71
− 0
View file @ fe678fff
function g = sw_g(LAT,z)
% SW_G Gravitational acceleration
%===========================================================================
% SW_G $Id: sw_g.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1993.
%
% USAGE: g = sw_g(lat,z)
%
% DESCRIPTION:
% Calculates acceleration due to gravity as function of latitude.
%
% INPUT: (all must have same dimensions)
% lat = Latitude in decimal degress north [-90..+90]
% z = height in metres (+ve above sea surface, -ve below)
%
% OUTPUT:
% g = gravity [m/s^2]
%
% AUTHOR: Phil Morgan 93-04-20 (morgan@ml.csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCES:
% Unesco 1983. Algorithms for computation of fundamental properties of
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
%
% A.E. Gill 1982. p.597
% "Atmosphere-Ocean Dynamics"
% Academic Press: New York. ISBN: 0-12-283522-0
%=========================================================================
% CALLER: general purpose
% CALLEE: none
%-------------
% CHECK INPUTS
%-------------
if ~(nargin==1 | nargin==2)
error('sw_g.m: Requires one or two input arguments')
end %if
if nargin == 1
z = zeros(size(LAT));
end %if
[mL,nL] = size(LAT);
[mz,nz] = size(z);
if ~(mL==mz | nL==nz)
error('sw_g.m: Input arguments should have same dimensions')
end %if
%-------------
% BEGIN
%-------------
% Eqn p27. Unesco 1983.
a = 6371000; % mean radius of earth A.E.Gill
DEG2RAD = pi/180;
LAT = abs(LAT);
X = sin(LAT*DEG2RAD); % convert to radians
sin2 = X.*X;
g = 9.780318*(1.0+(5.2788E-3+2.36E-5*sin2).*sin2);
if any(any(z))
g = g./((1+z/a).^2); % from A.E.Gill p.597
end %if
return
%===========================================================================
lib/vendor/seawater_library/sw_gpan.m 0 → 100644
+ 114
− 0
View file @ fe678fff
function [ga] = sw_gpan(S,T,P)
% SW_GPAN Geopotential anomaly
%=========================================================================
% SW_GPAN $Id: sw_gpan.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1992.
%
% USAGE: [gpan]= sw_gpan(S,T,P)
%
% DESCRIPTION:
% Geopotential Anomaly calculated as the integral of svan from the
% the sea surface to the bottom. Thus RELATIVE TO SEA SURFACE.
%
% INPUT: (all must have same dimensions)
% S = salinity [psu (PSS-78)]
% T = temperature [degree C (ITS-90)]
% P = Pressure [db]
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
% gpan = Geopotential Anomaly [m^3 kg^-1 Pa == m^2 s^-2 == J kg^-1]
%
% AUTHOR: Phil Morgan 92-11-05, Lindsay Pender (Lindsay.Pender@csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE: S. Pond & G.Pickard 2nd Edition 1986
% Introductory Dynamical Oceanogrpahy
% Pergamon Press Sydney. ISBN 0-08-028728-X
%
% Note that older literature may use units of "dynamic decimeter' for above.
%
% Adapted method from Pond and Pickard (p76) to calc gpan rel to sea
% surface whereas P&P calculated relative to the deepest common depth.
%=========================================================================
% Modifications
% 03-12-12. Lindsay Pender, Converted to ITS-90.
%
% CALLER: general purpose
% CALLEE: sw_svan.m
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin ~=3
error('Must pass 3 parameters')
end %if
% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
error('S & T must have same dimensions')
end %if
% CHECK OPTIONAL SHAPES FOR P
if mp==1 & np==1 % P is a scalar. Fill to size of S
P = P(1)*ones(ms,ns);
elseif np==ns & mp==1 % P is row vector with same cols as S
P = P( ones(1,ms), : ); % Copy down each column.
elseif mp==ms & np==1 % P is column vector
P = P( :, ones(1,ns) ); % Copy across each row
elseif mp==ms & np==ns % PR is a matrix size(S)
% shape ok
else
error('P has wrong dimensions')
end %if
[mp,np] = size(P);
% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
Transpose = 0;
if mp == 1 % row vector
P = P(:);
T = T(:);
S = S(:);
Transpose = 1;
end %if
%***check_stp
%------
% BEGIN
%------
db2Pascal = 1e4;
[m,n] = size(P);
svan = sw_svan(S,T,P);
mean_svan = 0.5*(svan(2:m,:) + svan(1:m-1,:) );
if n==1
top = svan(1,1).*P(1,1)*db2Pascal;
else
top = svan(1,:).*P(1,:)*db2Pascal;
end %if
%press_diff = diff(P);
delta_ga = (mean_svan.*diff(P))*db2Pascal;
ga = cumsum([ top; delta_ga]);
if Transpose
ga = ga';
end %if
lib/vendor/seawater_library/sw_gvel.m 0 → 100644
+ 63
− 0
View file @ fe678fff
function vel = sw_gvel(ga,lat,lon)
% SW_GVEL Geostrophic velocity
%===================================================================
% GEOVEL $Id: sw_gvel.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1992
%
% USAGE: vel = sw_gvel(ga,lat,lon)
%
% DESCRIPTION:
% Calculates geostrophic velocity given the geopotential anomaly
% and position of each station.
%
% INPUT:
% ga = geopotential anomoly relative to the sea surface.
% dim(mxnstations)
% lat = latitude of each station (+ve = N, -ve = S) [ -90.. +90]
% lon = longitude of each station (+ve = E, -ve = W) [-180..+180]
%
% OUTPUT:
% vel = geostrophic velocity RELATIVE to the sea surface.
% dim(m,nstations-1)
%
% AUTHOR: Phil Morgan 1992/03/26 (morgan@ml.csiro.au)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE: S. Pond & G.Pickard 2nd Edition 1986
% Introductory Dynamical Oceanogrpahy
% Pergamon Press Sydney. ISBN 0-08-028728-X
% Equation 8.9A p73 Pond & Pickard
%
% NOTE: This calls sw_dist.m. You can replace the call to this
% routine if you have a more appropraite distance routine.
%==================================================================
% CALLER: general purpose
% CALLEE: sw_dist.m
%
DEG2RAD = pi/180;
RAD2DEG = 180/pi;
OMEGA = 7.292e-5; % Angular velocity of Earth [radians/sec]
% You may replace the call to sw_dist if you have
% a more appropriate distance routine.
distm = 1000*sw_dist(lat,lon,'km');
[m,n] = size(ga);
f = 2*OMEGA*sin( (lat(1:n-1)+lat(2:n))*DEG2RAD/2 );
lf = f.*distm;
LF = lf(ones(m,1),:);
vel = -( ga(:,2:n)-ga(:,1:n-1) ) ./ LF;
return
%--------------------------------------------------------------------
lib/vendor/seawater_library/sw_info.m 0 → 100644
+ 31
− 0
View file @ fe678fff
% SW_INFO Computational routines for the properties of sea water
%
% SEAWATER - developed by Phil Morgan, CSIRO
%
% DESCRIPTION:
% SEAWATER is a toolkit of MATLAB routines for calculating the
% properties of sea water. They are a self contained library and
% are extremely easy to use and will run on all computers that
% support MATLAB.
%
% MATLAB:
% For information on MATLAB contact info@mathworks.com
%
% DISCLAIMER
% This software is provided "as is" without warranty of any kind.
% See the file sw_copy.m for conditions of use and licence.
%
% MORE INFORMATION:
% http://www.marine.csiro.au/~morgan/seawater
%
% $Id: sw_info.m,v 1.1 2003/12/12 04:23:22 pen078 Exp $
% Copyright (C) CSIRO, Phil Morgan 1993.
%=========================================================================
more on
help sw_info
more off
return
%--------------------------------------------------------------------
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment