Skip to content

Commit

Permalink
Version 3.0
Browse files Browse the repository at this point in the history
Major improvements

48-term density equation introduced (intended for use by observational,
theoretical oceanographers and ocean modellers). The use of a single
density equation will ensure a unified field.

the Absolute Salinity Anomaly takes into account evaporation and
dilution.
  • Loading branch information
PaulMBarker committed May 29, 2015
1 parent f0c5be3 commit ce64332
Show file tree
Hide file tree
Showing 209 changed files with 17,422 additions and 6,121 deletions.
238 changes: 238 additions & 0 deletions Toolbox/Contents.m

Large diffs are not rendered by default.

43 changes: 43 additions & 0 deletions Toolbox/gsw_Abs_Pressure_from_p.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
function Absolute_Pressure = gsw_Abs_Pressure_from_p(p)

% gsw_Abs_Pressure_from_p Absolute Pressure
%==========================================================================
%
% USAGE:
% Absolute_Pressure = gsw_Abs_Pressure_from_p(p)
%
% DESCRIPTION:
% Calculates Absolute Pressure from sea pressure. Note that Absolute
% Pressure is in Pa NOT dbar.
%
% INPUT:
% p = sea pressure [ dbar ]
%
% OUTPUT:
% Absolute_Pressure = Absolute Pressure [ Pa ]
%
% AUTHOR:
% Trevor McDougall and Paul Barker [ help_gsw@csiro.au ]
%
% VERSION NUMBER: 3.0 (29th March, 2011)
%
% REFERENCES:
% IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of
% seawater - 2010: Calculation and use of thermodynamic properties.
% Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
% UNESCO (English), 196 pp. Available from http://www.TEOS-10.org.
% See Eqn. (2.2.1) of this TEOS-10 Manual.
%
% The software is available from http://www.TEOS-10.org
%
%==========================================================================

if ~(nargin == 1)
error('gsw_Abs_Pressure_from_p: Requires one input')
end

db2Pa = 1e4;

Absolute_Pressure = p*db2Pa + 101325;

end
37 changes: 37 additions & 0 deletions Toolbox/gsw_C3515.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
function C3515 = gsw_C3515

% gsw_C3515 Conductivity of SSW at SP=35,t_68=15,p=0
%==========================================================================
%
% USAGE:
% C3515 = gsw_C3515
%
% DESCRIPTION:
% This function provides the present estimate of Conductivity, C, of
% Standard Seawater (SSW) at (SP=35, t_68=15, p=0) which is
% 42.9140 mS/cm (=4.29140 S/m) (Culkin and Smith, 1980; UNESCO, 1983).
%
% OUTPUT:
% C3515 = Conductivity at (SP=35, t_68=15, p=0) [ mS/cm ]
%
% AUTHOR:
% Trevor McDougall and Paul Barker [ help_gsw@csiro.au ]
%
% VERSION NUMBER: 3.0 (29th March, 2011)
%
% REFERENCES:
% Culkin and Smith, 1980: Determination of the Concentration of Potassium
% Chloride Solution Having the Same Electrical Conductivity, at 15C and
% Infinite Frequency, as Standard Seawater of Salinity 35.0000
% (Chlorinity 19.37394), IEEE J. Oceanic Eng, 5, 22-23.
%
% Unesco, 1983: Algorithms for computation of fundamental properties of
% seawater. Unesco Technical Papers in Marine Science, 44, 53 pp.
%
% The software is available from http://www.TEOS-10.org
%
%==========================================================================

C3515 = 42.9140;

end
59 changes: 40 additions & 19 deletions Toolbox/gsw_CT_first_derivatives.m
Original file line number Diff line number Diff line change
Expand Up @@ -25,29 +25,23 @@
% CT_SA = The derivative of Conservative Temperature with respect to
% Absolute Salinity at constant potential temperature
% (the regular potential temperature which has reference
% sea pressure of 0 dbar). The CT_SA output has units of:
% [ K/(g/kg)]
% sea pressure of 0 dbar).
% The CT_SA output has units of: [ K/(g/kg)]
% CT_pt = The derivative of Conservative Temperature with respect to
% potential temperature (the regular one with pr = 0 dbar)
% at constant SA. CT_pt is dimensionless. [ unitless ]
%
% AUTHOR:
% Trevor McDougall and Paul Barker [ help_gsw@csiro.au ]
% Trevor McDougall and Paul Barker [ help_gsw@csiro.au ]
%
% VERSION NUMBER: 2.0 (26th August, 2010)
% VERSION NUMBER: 3.0 (11th April 2011)
%
% REFERENCES:
% IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of
% seawater - 2010: Calculation and use of thermodynamic properties.
% Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
% UNESCO (English), 196 pp. Available from http://www.TEOS-10.org.
% See Eqns. (A.12.3) and (A.12.9a,b) of this TEOS-10 Manual.
%
% McDougall T. J., D. R. Jackett, P. M. Barker, C. Roberts-Thomson, R.
% Feistel and R. W. Hallberg, 2010: A computationally efficient 25-term
% expression for the density of seawater in terms of Conservative
% Temperature, and related properties of seawater. To be submitted
% to Ocean Science Discussions.
% See Eqns. (A.12.3a,b) and (A.15.8) of this TEOS-10 Manual.
%
% This software is available from http://www.TEOS-10.org
%
Expand All @@ -73,8 +67,8 @@
end

if ms == 1
SA = SA';
pt = pt';
SA = SA.';
pt = pt.';
transposed = 1;
else
transposed = 0;
Expand All @@ -86,19 +80,46 @@

cp0 = 3991.86795711963; % from Eqn. 3.3.3 of IOC et al. (2010).
n0 = 0;
n1 = 1;
n2 = 2;
pr0 = zeros(size(SA));
abs_pt = 273.15 + pt;

CT_SA = (gsw_gibbs(n1,n0,n0,SA,pt,pr0) -...
abs_pt.*gsw_gibbs(n1,n1,n0,SA,pt,pr0))./cp0;

CT_pt = - (abs_pt.*gsw_gibbs(n0,n2,n0,SA,pt,pr0))./cp0;

%--------------------------------------------------------------------------

sfac = 0.0248826675584615; % sfac = 1/(40*(35.16504/35)).
x2 = sfac.*SA;
x = sqrt(x2);
y_pt = 0.025*pt;

g_SA_T_mod = 1187.3715515697959 + ...
x.*(-1480.222530425046 + x.*(2175.341332000392 + x.*(-980.14153344888 + 220.542973797483.*x) + ...
y_pt.*(-548.4580073635929 + y_pt.*(592.4012338275047 + y_pt.*(-274.2361238716608 + 49.9394019139016.*y_pt)))) + ...
y_pt.*(-258.3988055868252 + y_pt.*(-90.2046337756875 + y_pt.*10.50720794170734))) + ...
y_pt.*(3520.125411988816 + y_pt.*(-1351.605895580406 + ...
y_pt.*(731.4083582010072 + y_pt.*(-216.60324087531103 + 25.56203650166196.*y_pt))));
g_SA_T_mod = 0.5*sfac*0.025*g_SA_T_mod;

g_SA_mod = 8645.36753595126 + ...
x.*(-7296.43987145382 + x.*(8103.20462414788 + ...
y_pt.*(2175.341332000392 + y_pt.*(-274.2290036817964 + ...
y_pt.*(197.4670779425016 + y_pt.*(-68.5590309679152 + 9.98788038278032.*y_pt)))) + ...
x.*(-5458.34205214835 - 980.14153344888.*y_pt + ...
x.*(2247.60742726704 - 340.1237483177863.*x + 220.542973797483.*y_pt))) + ...
y_pt.*(-1480.222530425046 + ...
y_pt.*(-129.1994027934126 + ...
y_pt.*(-30.0682112585625 + y_pt.*(2.626801985426835 ))))) + ...
y_pt.*(1187.3715515697959 + ...
y_pt.*(1760.062705994408 + y_pt.*(-450.535298526802 + ...
y_pt.*(182.8520895502518 + y_pt.*(-43.3206481750622 + 4.26033941694366.*y_pt)))));
g_SA_mod = 0.5*sfac*g_SA_mod;

CT_SA = (g_SA_mod - abs_pt.*g_SA_T_mod)./cp0;

if transposed
CT_SA = CT_SA';
CT_pt = CT_pt';
CT_SA = CT_SA.';
CT_pt = CT_pt.';
end

end
177 changes: 177 additions & 0 deletions Toolbox/gsw_CT_freezing.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,177 @@
function CT_freezing = gsw_CT_freezing(SA,p,saturation_fraction)

% gsw_CT_freezing Conservative Temperature at which seawater freezes
%==========================================================================
%
% USAGE:
% CT_freezing = gsw_CT_freezing(SA,p,saturation_fraction)
%
% DESCRIPTION:
% Calculates the Conservative Temperature at which seawater freezes.
%
% INPUT:
% SA = Absolute Salinity [ g/kg ]
% p = sea pressure [ dbar ]
% ( i.e. absolute pressure - 10.1325 dbar )
%
% OPTIONAL:
% saturation_fraction = the saturation fraction of dissolved air in
% seawater
% (i.e., saturation_fraction must be between 0 and 1, and the default
% is 1, completely saturated)
%
% p & saturation_fraction (if provided) may have dimensions 1x1 or Mx1 or
% 1xN or MxN, where SA is MxN.
%
% OUTPUT:
% CT_freezing = Conservative Temperature at freezing of seawater [ deg C ]
% That is, the freezing temperature expressed in
% terms of Conservative Temperature (ITS-90).
%
% AUTHOR:
% Trevor McDougall, Paul Barker and Rainer Feistal [ help@teos-10.org ]
%
% VERSION NUMBER: 3.0 (4th November, 2011)
%
% REFERENCES:
% IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of
% seawater - 2010: Calculation and use of thermodynamic properties.
% Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
% UNESCO (English), 196 pp. Available from http://www.TEOS-10.org.
% See sections 3.33 and 3.34 of this TEOS-10 Manual.
%
% The software is available from http://www.TEOS-10.org
%
%==========================================================================

%--------------------------------------------------------------------------
% Check variables and resize if necessary
%--------------------------------------------------------------------------

if ~(nargin == 2 | nargin == 3)
error('gsw_CT_freezing: Requires either two or three inputs')
end %if

if ~exist('saturation_fraction','var')
saturation_fraction = 1;
end

if (saturation_fraction < 0 | saturation_fraction > 1)
error('gsw_CT_freezing: saturation_fraction MUST be between zero and one.')
end

[ms,ns] = size(SA);
[mp,np] = size(p);
[map,nap] = size(saturation_fraction);

if (mp == 1) & (np == 1) % p scalar - fill to size of SA
p = p*ones(size(SA));
elseif (ns == np) & (mp == 1) % p is row vector,
p = p(ones(1,ms), :); % copy down each column.
elseif (ms == mp) & (np == 1) % p is column vector,
p = p(:,ones(1,ns)); % copy across each row.
elseif (ns == mp) & (np == 1) % p is a transposed row vector,
p = p.'; % transposed then
p = p(ones(1,ms), :); % copy down each column.
elseif (ms == mp) & (ns == np)
% ok
else
error('gsw_CT_freezing: Inputs array dimensions arguments do not agree')
end %if

if (map == 1) & (nap == 1) % saturation_fraction scalar
saturation_fraction = saturation_fraction*ones(size(SA)); % fill to size of SA
elseif (ns == nap) & (map == 1) % saturation_fraction is row vector,
saturation_fraction = saturation_fraction(ones(1,ms), :); % copy down each column.
elseif (ms == map) & (nap == 1) % saturation_fraction is column vector,
saturation_fraction = saturation_fraction(:,ones(1,ns)); % copy across each row.
elseif (ns == map) & (nap == 1) % saturation_fraction is a transposed row vector,
saturation_fraction = saturation_fraction.'; % transposed then
saturation_fraction = saturation_fraction(ones(1,ms), :); % copy down each column.
elseif (ms == map) & (ns == nap)
% ok
else
error('gsw_CT_freezing: Inputs array dimensions arguments do not agree')
end %if

if ms == 1
SA = SA.';
p = p.';
saturation_fraction = saturation_fraction.';
transposed = 1;
else
transposed = 0;
end

%--------------------------------------------------------------------------
% Start of the calculation
%--------------------------------------------------------------------------

% These few lines ensure that SA is non-negative.
[I_neg_SA] = find(SA < 0);
if ~isempty(I_neg_SA)
error(' gsw_CT_freezing: SA must be non-negative!')
end

c0 = 0.017947064327968736;
%
c1 = -6.076099099929818;
c2 = 4.883198653547851;
c3 = -11.88081601230542;
c4 = 13.34658511480257;
c5 = -8.722761043208607;
c6 = 2.082038908808201;
%
c7 = -7.389420998107497;
c8 = -2.110913185058476;
c9 = 0.2295491578006229;
%
c10 = -0.9891538123307282;
c11 = -0.08987150128406496;
c12 = 0.3831132432071728;
c13 = 1.054318231187074;
c14 = 1.065556599652796;
c15 = -0.7997496801694032;
c16 = 0.3850133554097069;
c17 = -2.078616693017569;
c18 = 0.8756340772729538;
c19 = -2.079022768390933;
c20 = 1.596435439942262;
c21 = 0.1338002171109174;
c22 = 1.242891021876471;

SA_r = SA.*1e-2;
x = sqrt(SA_r);
p_r = p.*1e-4;

CT_freezing = c0 ...
+ SA_r.*(c1 + x.*(c2 + x.*(c3 + x.*(c4 + x.*(c5 + c6.*x))))) ...
+ p_r.*(c7 + p_r.*(c8 + c9.*p_r)) ...
+ SA_r.*p_r.*(c10 + p_r.*(c12 + p_r.*(c15 + c21.*SA_r)) + SA_r.*(c13 + c17.*p_r + c19.*SA_r) ...
+ x.*(c11 + p_r.*(c14 + c18.*p_r) + SA_r.*(c16 + c20.*p_r + c22.*SA_r)));

% The error of this fit ranges between -5e-4 K and 6e-4 K when compared
% with the Conservative Temperature calculated from the exact in-situ
% freezing temperature which is found by a Newton-Raphson iteration of the
% equality of the chemical potentials of water in seawater and in ice.
% (Note that the in-situ freezing temperature can be found by this exact
% method using the function sea_ice_freezingtemperature_si in the SIA
% library).

% Adjust for the effects of dissolved air
a = 0.014289763856964; % Note that a = 0.502500117621/35.16504.
b = 0.057000649899720;
CT_freezing = CT_freezing ...
- saturation_fraction.*(1e-3).*(2.4 - a.*SA).*(1 + b.*(1 - SA./35.16504));

[Iout_of_range] = find(p > 10000 | SA > 120 | ...
p + SA.*71.428571428571402 > 13571.42857142857);
if ~isempty(Iout_of_range)
CT_freezing(Iout_of_range) = NaN;
end

if transposed
CT_freezing = CT_freezing.';
end

end
Loading

0 comments on commit ce64332

Please sign in to comment.