-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgsw_geostrophic_velocity.m
188 lines (168 loc) · 8.1 KB
/
gsw_geostrophic_velocity.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
function [geostrophic_velocity, mid_lat, mid_long] = gsw_geostrophic_velocity(geo_strf,long,lat,p)
% gsw_geostrophic_velocity geostrophic velocity
%==========================================================================
% USAGE:
% [geostrophic_velocity, mid_lat, mid_long] = ...
% gsw_geostrophic_velocity(geo_strf,long,lat,p)
%
% DESCRIPTION:
% Calculates geostrophic velocity relative to the sea surface, given a
% geostrophic streamfunction and the position (longitude, latitude and
% pressure (long, lat & p)) of each station in sequence along an ocean
% section. The data can be from a single isobaric or "density" surface,
% or from a series of such surfaces.
%
% INPUT:
% geo_strf = geostrophic streamfunction. This geostrophic streamfunction
% can be any of, for example,
% (1) geo_strf_dyn_height_CT (in an isobaric surface)
% (2) geo_strf_Montgomery_CT (in a specific volume anomaly surface)
% (3) geo_strf_Cunninhgam_CT (in an approximately neutral surface
% e.g. a potential denisty surface).
% (4) geo_strf_isopycnal_CT (in an approximately neutral surface
% e.g. a potential denisty surface, a Neutral Density
% surface or an omega surface (Klocker et al., 2009)).
%
% long = longitude in decimal degrees [ 0 ... +360 ]
% or [ -180 ... +180 ]
% lat = latitude in decimal degrees north [ -90 ... +90 ]
%
% There needs to be more than one station. The input geo_strf has
% dimensions (M("bottles") x N(stations)); that is, geo_strf has
% dimensions (M(surfaces) x N(stations)).
%
% Note. The ith "bottle" of each station (i.e. the ith row of geo_strf)
% must be on the same ith surface, whether that surface be,
% (1) an isobaric surface,
% (2) a specific volume anomaly surface,
% or some type of approximately neutral surface (cases (3) & (4)).
%
% lat & long need to have dimensions 1xN or MxN, where geo_strf is MxN.
% p may have dimensions 1x1 or Mx1 or 1xN or MxN, where geo_strf is MxN.
%
% OPTIONAL:
% p = sea pressure ( default is 0 ) [ dbar ]
% ( i.e. absolute pressure - 10.1325 dbar )
% Note. This optional input is used to obtain an accurate distance,
% "dist", taking into account that the radius from the centre of the
% Earth depends on the depth below the sea surface.
%
% OUTPUT:
% geostrophic_velocity = geostrophic velocity RELATIVE to the sea surface.
% It has dimensions (Mx(N-1)) and the relative
% geostrophic velocity is located at mid_long,
% mid_lat (and at the mid-point pressure).
% mid_long = mid point longitude
% (the range corresponds to that entered)
% mid_lat = mid point latitude, [ -90 ... +90 ]
% (in decimal degrees north)
%
%
% AUTHOR:
% Paul Barker, Trevor McDougall and Phil Morgan [ help@teos-10.org ]
%
% VERSION NUMBER: 3.05 (27th January 2015)
%
% REFERENCES:
% Cunningham, S. A., 2000: Circulation and volume flux of the North
% Atlantic using syoptic hydrographic data in a Bernoulli inverse.
% J. Marine Res., 58, 1-35.
%
% 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.27 - 3.3.30 of this TEOS-10 Manual.
%
% Jackett, D. R. and T. J. McDougall, 1997: A neutral density variable
% for the world’s oceans. Journal of Physical Oceanography, 27, 237-263.
%
% Klocker, A., T. J. McDougall and D. R. Jackett, 2009: A new method for
% forming approximately neutral surfaces. Ocean Sci., 5, 155-172.
%
% McDougall, T. J. and A. Klocker, 2010: An approximate geostrophic
% streamfunction for use in density surfaces. Ocean Modelling, 32,
% 105-117.
% See Eqn. (62), of this paper, for definition of the McDougall-Klocker
% geostrophic streamfunction.
%
% Montgomery, R. B., 1937: A suggested method for representing gradient
% flow in isentropic surfaces. Bull. Amer. Meteor. Soc. 18, 210-212.
%
% The software is available from http://www.TEOS-10.org
%
%==========================================================================
%--------------------------------------------------------------------------
% Check variables and resize if necessary
%--------------------------------------------------------------------------
if ~(nargin == 3 | nargin == 4)
error('gsw_geostrophic_velocity: Requires three or four inputs')
end %if
[mg,ng] = size(geo_strf);
[mlo,nlo] = size(long);
[mla,nla] = size(lat);
if (nlo == 1) & (mlo == 1)
error('gsw_geostrophic_velocity: need more than one station')
elseif (ng == nlo) & (mlo == 1) % long is row vector,
long = long(ones(1,mg), :); % copy down each column.
elseif (ng == mlo) & (nlo == 1) % long is a transposed row vector,
long = long.'; % transposed then
long = long(ones(1,mg), :); % copy down each column.
elseif (mg == mlo) & (ng == nlo)
% ok
else
error('gsw_geostrophic_velocity: Inputs array dimensions arguments do not agree')
end %if
long(long < 0) = long(long < 0) + 360;
if (nla == 1) & (mla == 1)
error('gsw_geostrophic_velocity: need more than one station')
elseif (ng == nla) & (mla == 1) % lat is row vector,
lat = lat(ones(1,mg), :); % copy down each column.
elseif (ng == mla) & (nla == 1) % long is a transposed row vector,
lat = lat.'; % transposed then
lat = lat(ones(1,mg), :); % copy down each column.
elseif (mg == mla) & (ng == nla)
% ok
else
error('gsw_geostrophic_velocity: Inputs array dimensions arguments do not agree')
end %if
if exist('p','var')
[mp,np] = size(p);
if (mp == 1) & (np == 1) % p scalar - fill to size of geo_str
p = p*ones(size(geo_strf));
elseif (ng == np) & (mp == 1) % p is row vector,
p = p(ones(1,mg), :); % copy down each column.
elseif (mg == mp) & (np == 1) % p is column vector,
p = p(:,ones(1,ng)); % copy across each row.
elseif (ng == mp) & (np == 1) % p is a transposed row vector,
p = p.'; % transposed then
p = p(ones(1,mg), :); % copy down each column.
elseif (mg == mp) & (ng == np)
% ok
else
error('gsw_geostrophic_velocity: Inputs array dimensions arguments do not agree')
end %if
else
p = zeros(size(geo_strf));
end
%--------------------------------------------------------------------------
% Start of the calculation
%--------------------------------------------------------------------------
dist = gsw_distance(long,lat,p); % Note that dist is in m (not km)
mid_lat = 0.5*(lat(:,1:ng-1) + lat(:,2:ng));
f = gsw_f(mid_lat); % This gets the Coriolis parameter
geostrophic_velocity = (geo_strf(:,2:ng) - geo_strf(:,1:ng-1))./(dist.*f);
mid_long = 0.5*(long(:,1:ng-1) + long(:,2:ng));
diff_long = (long(1,1:ng-1) - long(1,2:ng));
long2 = long;
long2(long2 > 180) = long2(long2 > 180) - 360;
mid_long2 = 0.5*(long2(:,1:ng-1) + long2(:,2:ng));
diff_long2 = (long2(1,1:ng-1) - long2(1,2:ng));
[dummy,gmc] = min([abs(diff_long); abs(diff_long2)]);
mid_long(gmc == 2) = mid_long2(gmc == 2);
mid_long(mid_long < 0) = mid_long(mid_long < 0) + 360;
mid_long(mid_long > 180) = mid_long(mid_long > 180) - 360;
% Note. This geostrophic velocity difference, v - v0, when positive, is
% directed to the left of the horizontal vector which points from one
% cast to the next.
end