1 |
gforget |
1.1 |
function PT = sw_ptmp(S,T,P,PR) |
2 |
|
|
|
3 |
|
|
% SW_PTMP Potential temperature |
4 |
|
|
%=========================================================================== |
5 |
|
|
% SW_PTMP $Revision: 1.3 $ $Date: 1994/10/10 05:45:13 $ |
6 |
|
|
% Copyright (C) CSIRO, Phil Morgan 1992. |
7 |
|
|
% |
8 |
|
|
% USAGE: ptmp = sw_ptmp(S,T,P,PR) |
9 |
|
|
% |
10 |
|
|
% DESCRIPTION: |
11 |
|
|
% Calculates potential temperature as per UNESCO 1983 report. |
12 |
|
|
% |
13 |
|
|
% INPUT: (all must have same dimensions) |
14 |
|
|
% S = salinity [psu (PSS-78) ] |
15 |
|
|
% T = temperature [degree C (IPTS-68)] |
16 |
|
|
% P = pressure [db] |
17 |
|
|
% PR = Reference pressure [db] |
18 |
|
|
% (P & PR may have dims 1x1, mx1, 1xn or mxn for S(mxn) ) |
19 |
|
|
% |
20 |
|
|
% OUTPUT: |
21 |
|
|
% ptmp = Potential temperature relative to PR [degree C (IPTS-68)] |
22 |
|
|
% |
23 |
|
|
% AUTHOR: Phil Morgan 92-04-06 (morgan@ml.csiro.au) |
24 |
|
|
% |
25 |
|
|
% DISCLAIMER: |
26 |
|
|
% This software is provided "as is" without warranty of any kind. |
27 |
|
|
% See the file sw_copy.m for conditions of use and licence. |
28 |
|
|
% |
29 |
|
|
% REFERENCES: |
30 |
|
|
% Fofonoff, P. and Millard, R.C. Jr |
31 |
|
|
% Unesco 1983. Algorithms for computation of fundamental properties of |
32 |
|
|
% seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp. |
33 |
|
|
% Eqn.(31) p.39 |
34 |
|
|
% |
35 |
|
|
% Bryden, H. 1973. |
36 |
|
|
% "New Polynomials for thermal expansion, adiabatic temperature gradient |
37 |
|
|
% and potential temperature of sea water." |
38 |
|
|
% DEEP-SEA RES., 1973, Vol20,401-408. |
39 |
|
|
%========================================================================= |
40 |
|
|
|
41 |
|
|
% CALLER: general purpose |
42 |
|
|
% CALLEE: sw_adtg.m |
43 |
|
|
|
44 |
|
|
%------------- |
45 |
|
|
% CHECK INPUTS |
46 |
|
|
%------------- |
47 |
|
|
if nargin ~= 4 |
48 |
|
|
error('sw_ptmp.m: Must pass 4 parameters ') |
49 |
|
|
end %if |
50 |
|
|
|
51 |
|
|
% CHECK S,T,P dimensions and verify consistent |
52 |
|
|
[ms,ns] = size(S); |
53 |
|
|
[mt,nt] = size(T); |
54 |
|
|
[mp,np] = size(P); |
55 |
|
|
|
56 |
|
|
[mpr,npr] = size(PR); |
57 |
|
|
|
58 |
|
|
|
59 |
|
|
% CHECK THAT S & T HAVE SAME SHAPE |
60 |
|
|
if (ms~=mt) | (ns~=nt) |
61 |
|
|
error('check_stp: S & T must have same dimensions') |
62 |
|
|
end %if |
63 |
|
|
|
64 |
|
|
% CHECK OPTIONAL SHAPES FOR P |
65 |
|
|
if mp==1 & np==1 % P is a scalar. Fill to size of S |
66 |
|
|
P = P(1)*ones(ms,ns); |
67 |
|
|
elseif np==ns & mp==1 % P is row vector with same cols as S |
68 |
|
|
P = P( ones(1,ms), : ); % Copy down each column. |
69 |
|
|
elseif mp==ms & np==1 % P is column vector |
70 |
|
|
P = P( :, ones(1,ns) ); % Copy across each row |
71 |
|
|
elseif mp==ms & np==ns % PR is a matrix size(S) |
72 |
|
|
% shape ok |
73 |
|
|
else |
74 |
|
|
error('check_stp: P has wrong dimensions') |
75 |
|
|
end %if |
76 |
|
|
[mp,np] = size(P); |
77 |
|
|
|
78 |
|
|
|
79 |
|
|
% CHECK OPTIONAL SHAPES FOR PR |
80 |
|
|
if mpr==1 & npr==1 % PR is a scalar. Fill to size of S |
81 |
|
|
PR = PR(1)*ones(ms,ns); |
82 |
|
|
elseif npr==ns & mpr==1 % PR is row vector with same cols as S |
83 |
|
|
PR = PR( ones(1,ms), : ); % Copy down each column. |
84 |
|
|
elseif mpr==ms & npr==1 % P is column vector |
85 |
|
|
PR = PR( :, ones(1,ns) ); % Copy across each row |
86 |
|
|
elseif mpr==ms & npr==ns % PR is a matrix size(S) |
87 |
|
|
% shape ok |
88 |
|
|
else |
89 |
|
|
error('check_stp: PR has wrong dimensions') |
90 |
|
|
end %if |
91 |
|
|
[mpr,npr] = size(PR); |
92 |
|
|
|
93 |
|
|
|
94 |
|
|
% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN. |
95 |
|
|
Transpose = 0; |
96 |
|
|
if mp == 1 % row vector |
97 |
|
|
P = P(:); |
98 |
|
|
T = T(:); |
99 |
|
|
S = S(:); |
100 |
|
|
|
101 |
|
|
PR = PR(:); |
102 |
|
|
|
103 |
|
|
Transpose = 1; |
104 |
|
|
end %if |
105 |
|
|
%***check_stp |
106 |
|
|
|
107 |
|
|
%------ |
108 |
|
|
% BEGIN |
109 |
|
|
%------ |
110 |
|
|
|
111 |
|
|
% theta1 |
112 |
|
|
del_P = PR - P; |
113 |
|
|
del_th = del_P.*sw_adtg(S,T,P); |
114 |
|
|
th = T + 0.5*del_th; |
115 |
|
|
q = del_th; |
116 |
|
|
|
117 |
|
|
% theta2 |
118 |
|
|
del_th = del_P.*sw_adtg(S,th,P+0.5*del_P); |
119 |
|
|
th = th + (1 - 1/sqrt(2))*(del_th - q); |
120 |
|
|
q = (2-sqrt(2))*del_th + (-2+3/sqrt(2))*q; |
121 |
|
|
|
122 |
|
|
% theta3 |
123 |
|
|
del_th = del_P.*sw_adtg(S,th,P+0.5*del_P); |
124 |
|
|
th = th + (1 + 1/sqrt(2))*(del_th - q); |
125 |
|
|
q = (2 + sqrt(2))*del_th + (-2-3/sqrt(2))*q; |
126 |
|
|
|
127 |
|
|
% theta4 |
128 |
|
|
del_th = del_P.*sw_adtg(S,th,P+del_P); |
129 |
|
|
PT = th + (del_th - 2*q)/6; |
130 |
|
|
|
131 |
|
|
if Transpose |
132 |
|
|
PT = PT'; |
133 |
|
|
end %if |
134 |
|
|
|
135 |
|
|
return |
136 |
|
|
%========================================================================= |