Skip to content
Snippets Groups Projects
Commit bd824018 authored by Yves Gouriou's avatar Yves Gouriou
Browse files

Développement de la focntion d'étalonnage.

Ajout de la bilbiothèque CSIRO à SVN (pas nécessaire)
parent 410b30a4
No related branches found
No related tags found
No related merge requests found
Showing with 1896 additions and 0 deletions
No preview for this file type
% SEAWATER Library
% Version 2.0.1 22-Apr-1998
%
% *******************************
% * SEAWATER Library *
% * *
% * Version 2.0.1 *
% * (for Matlab 5.x) *
% * *
% * *
% * Phillip P. Morgan *
% * CSIRO *
% * *
% * Phil.Morgan@marine.csiro.au *
% *******************************
%
% LIST OF ROUTINES:
%
% SW_NEW What's new in this version of seawater.
%
% SW_ADTG Adiabatic temperature gradient
% SW_ALPHA Thermal expansion coefficient (alpha)
% SW_AONB Calculate alpha/beta (a on b)
% SW_BETA Saline contraction coefficient (beta)
% SW_BFRQ Brunt-Vaisala Frequency Squared (N^2)
% SW_COPY Copyright and Licence file
% SW_CP Heat Capacity (Cp) of Sea Water
% SW_DENS Density of sea water
% SW_DENS0 Denisty of sea water at atmospheric pressure
% SW_DIST Distance between two lat, lon coordinates
% SW_DPTH Depth from pressure
% SW_F Coriolis factor "f"
% SW_FP Freezing Point of sea water
% SW_G Gravitational acceleration
% SW_GPAN Geopotential anomaly
% SW_GVEL Geostrophic velocity
% SW_INFO Information on the SEAWATER library.
% SW_PDEN Potential Density
% SW_PRES Pressure from depth
% SW_PTMP Potential temperature
% SW_SALS Salinity of sea water
% SW_SALT Salinity from cndr, T, P
% SW_SATAr Solubility (saturation) of Ar in seawater
% SW_SATN2 Solubility (saturation) of N2 in seawater
% SW_SATO2 Solubility (saturation) of O2 in seawater
% SW_SVAN Specific volume anomaly
% SW_SVEL Sound velocity of sea water
% SW_SMOW Denisty of standard mean ocean water (pure water)
% SW_TEMP Temperature from potential temperature
% SW_TEST Run test suite on library
% SW_VER Version number of SEAWATER library
%
% LOW LEVEL ROUTINES CALLED BY ABOVE: (also available for you to use)
%
% SW_C3515 Conductivity at (35,15,0)
% SW_CNDR Conductivity ratio R = C(S,T,P)/C(35,15,0)
% SW_SALDS Differiential dS/d(sqrt(Rt)) at constant T.
% SW_SALRP Conductivity ratio Rp(S,T,P) = C(S,T,P)/C(S,T,0)
% SW_SALRT Conductivity ratio rt(T) = C(35,T,0)/C(35,15,0)
% SW_SECK Secant bulk modulus (K) of sea water
%=======================================================================
% Contents.m $Revision: 1.6 $ $Date: 1998/04/22 02:12:17 $
% README file $Revision: 1.1 $ $Date: 1994/10/11 02:54:21 $
%
% 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
function [n2, e] = bm_bvfrq( S, T, P, LAT, DP)
%*************************************************************
% ***** brunt-vaisala freq *****
%
% uses 1980 equation of state
%
% in :
% P pressure decibars
% T temperature deg celsius (ipts-68)
% S salinity psu (pss-78)
% LAT latitude
% DP intervalle de pression (db) sur lequel est calcule n2
%
% out :
% n2 bouyancy freq cph
% e stability 1/m
%
% r. millard feb. 1991
%
%*************************************************************
%------
% BEGIN
%------
% change call list xlat to glat & pass gravity m/s^2
% after formulation of n. p. fofonoff & breck owen's
%
% Calcul de la gravit
% note that sw_g expects height as argument
Z = sw_dpth(P,LAT);
gp = sw_g(LAT,-Z);
g2 = gp .* gp;
g2 = g2 * 1.e-4;
% compute least squares estimate of specific volume anamoly gradient
[m,n] = size(P);
n2 = zeros(size(P));
% Calcul du nombre d'observation sur lequel est calcule n2
dp = P(2)-P(1);
nobs = DP / dp;
% Debut du profil
for i = 1:(nobs/2)
j = 1:i+nobs/2;
[n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end
% Fin du profil
for i = m-(nobs/2):m
j = i-nobs/2:m;
[n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end
% Profil
for i = (nobs/2+1):m-(nobs/2+1)
j = i-nobs/2:i+nobs/2;
[n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end
% Corps du calcul
function [n2, e] = cal_bvfrq( S, T, P, g2, gp, nobs, j)
% default gravity & rad/sec to cph conversion
%
radsec = 572.9578;
pav = mean(P(j));
data = sw_svan( S(j), sw_ptmp( S(j), T(j), P(j), pav), pav);
cxy = sum( data .* (P(j)-pav) );
cy = sum( data );
cxx = sum( (P(j)-pav).*(P(j)-pav) );
sigma = sw_pden(S(j), sw_ptmp( S(j), T(j), P(j), pav), P(j), pav);
a0 = cxy/cxx;
v350p = 1./sigma - data;
vbar = v350p(end) + cy/(nobs+1);
dvdp = a0;
e = -g2*dvdp/(vbar)^2;
n2 = radsec*sqrt(abs(e))*sign(e);
% define stability parameter units (1/m)
e = e/gp;
return
\ No newline at end of file
0 25.698 35.221
10 26.673 36.106
20 26.678 36.106
30 26.676 36.107
50 24.528 36.561
75 22.753 36.614
100 21.427 36.637
125 20.633 36.627
150 19.522 36.558
200 18.798 36.555
250 18.431 36.537
300 18.189 36.526
400 17.726 36.477
500 17.165 36.381
600 15.592 36.105
700 13.458 35.766
800 11.109 35.437
900 8.798 35.178
1000 6.292 35.044
1100 5.249 35.004
1200 4.813 34.995
1300 4.554 34.986
1400 4.357 34.977
1500 4.245 34.975
1750 4.028 34.973
2000 3.852 34.975
2500 3.424 34.968
3000 2.963 34.946
3500 2.462 34.920
4000 2.259 34.904
4327 2.221 34.896
\ No newline at end of file
File added
function ADTG = sw_adtg(S,T,P)
% SW_ADTG Adiabatic temperature gradient
%===========================================================================
% SW_ADTG $Revision: 1.4 $ $Date: 1994/10/10 04:16:37 $
% 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 (IPTS-68)]
% 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 92-04-03 (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:
% 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.
%=========================================================================
%-------------
% 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
[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
%-------------
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.*T).*T).*T ...
+ (b0 + b1.*T).*(S-35) ...
+ ( (c0 + (c1 + (c2 + c3.*T).*T).*T) + (d0 + d1.*T).*(S-35) ).*P ...
+ ( e0 + (e1 + e2.*T).*T ).*P.*P;
if Transpose
ADTG = ADTG';
end %if
return
%==========================================================================
function [ALPHA] = sw_alpha(S, T, P, keyword)
% SW_ALPHA Thermal expansion coefficient (alpha)
%================================================================
% SW_ALPHA $Revision: 1.6 $ $Date: 1998/04/21 05:42:10 $
% 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 (IPTS-68)]
% * T = temperature [degree C (IPTS-68)]
% 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
%
% 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'
% 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
[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
ALPHA = sw_aonb(S,T,P,keyword).*sw_beta(S,T,P,keyword);
return
%------------------------------------------------------------------------
function [AONB]= aonb(S, T, P, keyword)
% SW_AONB Calculate alpha/beta (a on b)
%================================================================
% SW_AONB $Revision: 1.5 $ $Date: 1994/11/15 04:10:45 $
% 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 (IPTS-68)]
% * T = temperature [degree C (IPTS-68)]
% 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
%
% 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'
% 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
[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
% ENSURE WE USE PTMP IN CALCULATIONS
if strcmp(lower(keyword),'ptmp')
% already have ptmp
else
T = sw_ptmp(S,T,P,0); % now have ptmp
end %if
% 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
%----------------------------------------------------------------------
function [BETA] = sw_beta(S, T, P, keyword)
% SW_BETA Saline contraction coefficient (beta)
%========================================================================
% SW_BETA $Revision: 1.4 $ $Date: 1994/11/15 04:10:05 $
% % 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 (IPTS-68)]
% * T = temperature [degree C (IPTS-68)]
% 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
%
% 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, 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'
% 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
[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
% ENSURE WE USE PTMP IN CALCULATIONS
if strcmp(lower(keyword),'ptmp')
% already have ptmp
else
T = sw_ptmp(S,T,P,0); % now have ptmp
end %if
% 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
%------------------------------------------------------------------------
function [n2,q,p_ave] = sw_bfrq(S,T,P,LAT)
% SW_BFRQ Brunt-Vaisala Frequency Squared (N^2)
%===========================================================================
% SW_BFRQ $Revision: 1.12 $ $Date: 1994/11/15 04:13:34 $
% 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 (IPTS-68)]
% 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 (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:
% 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
%=========================================================================
% CALLER: general purpose
% CALLEE: sw_dens.m sw_pden.m
%$Id: sw_bfrq.M,v 1.12 1994/11/15 04:13:34 morgan Exp $
%-------------
% CHECK INPUTS
%-------------
if ~(nargin==3 | nargin==4)
error('sw_bfrq.m: Must pass 3 or 4 parameters ')
end %if
if nargin == 3
LAT = [];
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
[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
% IF LAT PASSED THEN VERIFY DIMENSIONS
if ~isempty(LAT)
[mL,nL] = size(LAT);
if mL==1 & nL==1
LAT = LAT*ones(size(S));
%end % Je commente le end et remplace le if suivant par elseif
elseif (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('sw_bfrq.m: Inputs arguments have wrong dimensions')
end %if
end %if
end %if
%------
% 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 %if
return
%-------------------------------------------------------------------
function [n2,q,p_ave] = sw_bfrq_mean(S,T,P,LAT)
% SW_BFRQ Brunt-Vaisala Frequency Squared (N^2)
%===========================================================================
% SW_BFRQ $Revision: 1.12 $ $Date: 1994/11/15 04:13:34 $
% 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 (IPTS-68)]
% 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 (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:
% 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
%=========================================================================
% CALLER: general purpose
% CALLEE: sw_dens.m sw_pden.m
%$Id: sw_bfrq.M,v 1.12 1994/11/15 04:13:34 morgan Exp $
%-------------
% CHECK INPUTS
%-------------
if ~(nargin==3 | nargin==4)
error('sw_bfrq.m: Must pass 3 or 4 parameters ')
end %if
if nargin == 3
LAT = [];
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
[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
% IF LAT PASSED THEN VERIFY DIMENSIONS
if ~isempty(LAT)
[mL,nL] = size(LAT);
if mL==1 & nL==1
LAT = LAT*ones(size(S));
%end % Je commente le end et remplace le if suivant par elseif
elseif (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('sw_bfrq.m: Inputs arguments have wrong dimensions')
end %if
end %if
end %if
%------
% 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);
% Calcul des moyennes
% -------------------
Mmid_pden = mean(mid_pden,2);
Mdif_pden = mean(dif_pden,2);
Mmid_g = mid_g(:,1);
Mdif_z = dif_z(:,1);
n2 = -Mmid_g .* Mdif_pden ./ (Mdif_z .* Mmid_pden);
mid_f = f(iup,:);
Mmid_f = mid_f(:,1);
q = Mmid_f .* Mdif_pden ./ (Mdif_z .* Mmid_pden);
if Transpose
n2 = n2';
q = q';
p_ave = p_ave';
end %if
return
%-------------------------------------------------------------------
function [n2,q,p_ave] = sw_bfrq(S,T,P,LAT)
% SW_BFRQ Brunt-Vaisala Frequency Squared (N^2)
%===========================================================================
% SW_BFRQ $Revision: 1.12 $ $Date: 1994/11/15 04:13:34 $
% 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 (IPTS-68)]
% 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 (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:
% 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
%=========================================================================
% CALLER: general purpose
% CALLEE: sw_dens.m sw_pden.m
%$Id: sw_bfrq.M,v 1.12 1994/11/15 04:13:34 morgan Exp $
%-------------
% CHECK INPUTS
%-------------
if ~(nargin==3 | nargin==4)
error('sw_bfrq.m: Must pass 3 or 4 parameters ')
end %if
if nargin == 3
LAT = [];
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
[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
% IF LAT PASSED THEN VERIFY DIMENSIONS
if ~isempty(LAT)
[mL,nL] = size(LAT);
if mL==1 & nL==1
LAT = LAT*ones(size(S));
end %if
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('sw_bfrq.m: Inputs arguments have wrong dimensions')
end %if
end %if
end %if
%------
% 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 %if
return
%-------------------------------------------------------------------
function c3515 = sw_c3515()
% SW_C3515 Conductivity at (35,15,0)
%=========================================================================
% SW_c3515 $Revision: 1.4 $ $Date: 1994/10/10 04:35:22 $
% % 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
%-------------------------------------------------------------------------
function R = sw_cndr(S,T,P)
% SW_CNDR Conductivity ratio R = C(S,T,P)/C(35,15,0)
%=========================================================================
% SW_CNDR $Revision: 1.3 $ $Date: 1994/10/10 04:36:58 $
% 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 (IPTS-68)]
% 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,0) [no units]
%
% AUTHOR: Phil Morgan 93-04-21 (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:
% 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.
%=========================================================================
% 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
[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
%-------
del_T = T - 15;
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,del_T(i,j));
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.*T);
B = 1 + d1.*T + d2.*T.^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;
return
%-----------------------------------------------------------------------
% SW_COPY Copyright and licence information on SEAWATER library.
% =================================================================
% SW_COPY $Revision: 1.6 $ $Date: 1994/12/06 03:05:54 $
% 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.
%
%=========================================================================
%$Id: sw_copy.M,v 1.6 1994/12/06 03:05:54 morgan Exp $
more on
help sw_copy
more off
return
%--------------------------------------------------------------------
function cp = sw_cp(S,T,P)
% SW_CP Heat Capacity (Cp) of sea water
%=========================================================================
% SW_CP $Revision: 1.3 $ $Date: 1994/10/10 04:38:05 $
% 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 (IPTS-68)]
% 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 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:
% 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.
%=========================================================================
% CALLER: general purpose
% CALLEE: none
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin ~=3
error('sw_cp.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
[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
%------
P = P/10; % to convert db to Bar as used in Unesco routines
%------------
% 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 = c0 + c1.*T + c2.*T.^2 + c3.*T.^3 + c4.*T.^4 + ...
(a0 + a1.*T + a2.*T.^2).*S + ...
(b0 + b1.*T + b2.*T.^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 = (a0 + a1.*T + a2.*T.^2 + a3.*T.^3 + a4.*T.^4).*P + ...
(b0 + b1.*T + b2.*T.^2 + b3.*T.^3 + b4.*T.^4).*P.^2 + ...
(c0 + c1.*T + c2.*T.^2 + c3.*T.^3).*P.^3;
%------------
% 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 = [(d0 + d1.*T + d2.*T.^2 + d3.*T.^3 + d4.*T.^4).*S + ...
(e0 + e1.*T + e2.*T.^2).*S3_2].*P + ...
[(f0 + f1.*T + f2.*T.^2 + f3.*T.^3).*S + ...
g0.*S3_2].*P.^2 + ...
[(h0 + h1.*T + h2.*T.^2).*S + ...
j1.*T.*S3_2].*P.^3;
cp = Cpst0 + del_Cp0t0 + del_Cpstp;
if Transpose
cp = cp';
end %if
return
%--------------------------------------------------------------------
File added
function dens = sw_dens(S,T,P)
% SW_DENS Density of sea water
%=========================================================================
% SW_DENS $Revision: 1.3 $ $Date: 1994/10/10 04:39:16 $
% 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 (IPTS-68)]
% 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 (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:
% 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.
%=========================================================================
% 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
[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
%------
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);
if Transpose
dens = dens';
end %if
return
%--------------------------------------------------------------------
function dens = sw_dens0(S,T)
% SW_DENS0 Denisty of sea water at atmospheric pressure
%=========================================================================
% SW_DENS0 $Revision: 1.3 $ $Date: 1994/10/10 04:54:09 $
% 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 (IPTS-68)]
%
% 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 (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.
%
% Millero, F.J. and Poisson, A.
% International one-atmosphere equation of state of seawater.
% Deep-Sea Res. 1981. Vol28A(6) pp625-629.
%=========================================================================
% 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
Transpose = 0;
if mS == 1 % a row vector
S = S(:);
T = T(:);
Transpose = 1;
end %if
%----------------------
% DEFINE CONSTANTS
%----------------------
% 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*T + b2*T.^2 + b3*T.^3 + b4*T.^4).*S ...
%$$$ + (c0 + c1*T + c2*T.^2).*S.*sqrt(S) + d0*S.^2;
dens = sw_smow(T) + (b0 + (b1 + (b2 + (b3 + b4*T).*T).*T).*T).*S ...
+ (c0 + (c1 + c2*T).*T).*S.*sqrt(S) + d0*S.^2;
if Transpose
dens = dens';
end %if
return
%--------------------------------------------------------------------
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment