-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgsw_specvol_first_derivatives_wrt_enthalpy_CT_exact.m
123 lines (107 loc) · 4.35 KB
/
gsw_specvol_first_derivatives_wrt_enthalpy_CT_exact.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
function [v_SA, v_h] = gsw_specvol_first_derivatives_wrt_enthalpy_CT_exact(SA,CT,p)
% gsw_specvol_first_derivatives_wrt_enthalpy_CT_exact first derivatives
% of specific volume with respect to enthalpy
% =========================================================================
%
% USAGE:
% [v_SA, v_h] = ...
% gsw_specvol_first_derivatives_wrt_enthalpy_CT_exact(SA,CT,p)
%
% DESCRIPTION:
% Calculates the following two first-order derivatives of specific
% volume (v),
% (1) v_SA, first-order derivative with respect to Absolute Salinity
% at constant h & p.
% (2) v_h, first-order derivative with respect to h at constant SA & p.
%
% Note that this function uses the full Gibbs function. There is an
% alternative to calling this function, namely
% gsw_specvol_first_derivatives_wrt_enthalpy(SA,CT,p) which uses the
% computationally efficient 75 term expression for specific volume in
% terms of SA, CT and p (Roquet et al., 2015).
%
% INPUT:
% SA = Absolute Salinity [ g/kg ]
% CT = Conservative Temperature (ITS-90) [ deg C ]
% p = sea pressure [ dbar ]
% ( i.e. absolute pressure - 10.1325 dbar )
%
% SA & CT need to have the same dimensions.
% p may have dimensions 1x1 or Mx1 or 1xN or MxN, where SA & CT are MxN.
%
% OUTPUT:
% v_SA = The first derivative of specific volume with respect to
% Absolute Salinity at constant h & p. [ J/(kg (g/kg)^2) ]
% v_h = The first derivative of specific volume with respect to
% h at constant SA & p. [ J/(kg K(g/kg)) ]
%
% AUTHOR:
% Trevor McDougall and Paul Barker. [ help@teos-10.org ]
%
% VERSION NUMBER: 3.05 (27th January 2015)
%
% 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.
%
% Roquet, F., G. Madec, T.J. McDougall, P.M. Barker, 2015: Accurate
% polynomial expressions for the density and specifc volume of seawater
% using the TEOS-10 standard. Ocean Modelling.
%
% This software is available from http://www.TEOS-10.org
%
%==========================================================================
%--------------------------------------------------------------------------
% Check variables and resize if necessary
%--------------------------------------------------------------------------
if ~(nargin == 3)
error('gsw_specvol_first_derivatives_wrt_enthalpy_CT_exact: Requires three inputs')
end %if
if ~(nargout == 2)
error('gsw_specvol_first_derivatives_wrt_enthalpy_CT_exact: Requires two outputs')
end %if
[ms,ns] = size(SA);
[mt,nt] = size(CT);
[mp,np] = size(p);
if (ms ~= mt | ns ~= nt )
error('gsw_specvol_first_derivatives_wrt_enthalpy_CT_exact: SA and CT do not have the same dimensions')
end %if
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_specvol_first_derivatives_wrt_enthalpy_CT_exact: The dimensions of p do not agree')
end %if
if ms == 1
SA = SA.';
CT = CT.';
p = p.';
transposed = 1;
else
transposed = 0;
end
%--------------------------------------------------------------------------
% Start of the calculation
%--------------------------------------------------------------------------
% This line ensures that SA is non-negative.
SA(SA < 0) = 0;
[vCT_SA, vCT_CT, dummy] = gsw_specvol_first_derivatives_CT_exact(SA,CT,p);
[h_SA, h_CT] = gsw_enthalpy_first_derivatives_CT_exact(SA,CT,p);
rec_h_CT = 1./h_CT;
v_SA = vCT_SA - (vCT_CT.*h_SA).*rec_h_CT;
v_h = vCT_CT.*rec_h_CT;
if transposed
v_SA = v_SA.';
v_h = v_h.';
end
end