1 |
% |
2 |
% [] = B_COMPUTE_RELATIVE_VORTICITY(SNAPSHOT) |
3 |
% |
4 |
% For a time snapshot, this program computes the |
5 |
% 3D relative vorticity field from 3D |
6 |
% horizontal speed fields U,V (x,y,z) as: |
7 |
% OMEGA = ( -dVdz ; dUdz ; dVdx - dUdy ) |
8 |
% = ( Ox ; Oy ; ZETA ) |
9 |
% 3 output files are created. |
10 |
% |
11 |
% 06/07/2006 |
12 |
% gmaze@mit.edu |
13 |
% |
14 |
|
15 |
function [] = B_compute_relative_vorticity(snapshot) |
16 |
|
17 |
|
18 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
19 |
% Setup |
20 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
21 |
global sla netcdf_UVEL netcdf_VVEL netcdf_domain netcdf_suff |
22 |
pv_checkpath |
23 |
|
24 |
|
25 |
%% U,V files name: |
26 |
filU = strcat(netcdf_UVEL,'.',netcdf_domain); |
27 |
filV = strcat(netcdf_VVEL,'.',netcdf_domain); |
28 |
|
29 |
|
30 |
%% Path and extension to find them: |
31 |
pathname = strcat('netcdf-files',sla,snapshot,sla); |
32 |
ext = strcat('.',netcdf_suff); |
33 |
|
34 |
|
35 |
%% Load files: |
36 |
ferfile = strcat(pathname,sla,filU,ext); |
37 |
ncU = netcdf(ferfile,'nowrite'); |
38 |
[Ulon Ulat Udpt] = coordfromnc(ncU); |
39 |
|
40 |
ferfile = strcat(pathname,sla,filV,ext); |
41 |
ncV = netcdf(ferfile,'nowrite'); |
42 |
[Vlon Vlat Vdpt] = coordfromnc(ncV); |
43 |
|
44 |
clear ext ferfile |
45 |
|
46 |
%% Optionnal flags |
47 |
computeZETA = 1; % Compute ZETA or not ? |
48 |
global toshow % Turn to 1 to follow the computing process |
49 |
|
50 |
|
51 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
52 |
% VERTICAL COMPONENT: ZETA % |
53 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
54 |
|
55 |
% U field is on the zonal side of the c-grid and |
56 |
% V field on the meridional one. |
57 |
% So computing meridional gradient for U and |
58 |
% zonal gradient for V makes the relative vorticity |
59 |
% zeta defined on the corner of the c-grid. |
60 |
|
61 |
%%%%%%%%%%%%%% |
62 |
%% Dimensions of ZETA field: |
63 |
if toshow,disp('Dim'),end |
64 |
ny = length(Ulat)-1; |
65 |
nx = length(Vlon)-1; |
66 |
nz = length(Udpt); % Note that Udpt=Vdpt |
67 |
ZETA_lon = Ulon(1:nx); |
68 |
ZETA_lat = Vlat(1:ny); |
69 |
|
70 |
%%%%%%%%%%%%%% |
71 |
%% Pre-allocation: |
72 |
if toshow,disp('Pre-allocate'),end |
73 |
ZETA = zeros(nz,ny-1,nx-1).*NaN; |
74 |
dx = zeros(ny-1,nx-1); |
75 |
dy = zeros(ny-1,nx-1); |
76 |
|
77 |
%%%%%%%%%%%%%% |
78 |
%% Compute relative vorticity for each z-level: |
79 |
if computeZETA |
80 |
for iz=1:nz |
81 |
if toshow |
82 |
disp(strcat('Computing \zeta at depth : ',num2str(Udpt(iz)),... |
83 |
'm (',num2str(iz),'/',num2str(nz),')' )); |
84 |
end |
85 |
|
86 |
% Get velocities: |
87 |
U = ncU{4}(iz,:,:); |
88 |
V = ncV{4}(iz,:,:); |
89 |
|
90 |
% And now compute the vertical component of relative vorticity: |
91 |
% (TO DO: m_lldist accepts tables as input, so this part may be |
92 |
% done without x,y loop ...) |
93 |
for iy = 1 : ny |
94 |
for ix = 1 : nx |
95 |
if iz==1 % It's more efficient to make this test each time than |
96 |
% recomputing distance each time. m_lldist is a slow routine. |
97 |
% ZETA axis and grid distance: |
98 |
dx(iy,ix) = m_lldist([Vlon(ix+1) Vlon(ix)],[1 1]*Vlat(iy)); |
99 |
dy(iy,ix) = m_lldist([1 1]*Vlon(ix),[Ulat(iy+1) Ulat(iy)]); |
100 |
end %if |
101 |
% Horizontal gradients and ZETA: |
102 |
dVdx = ( V(iy,ix+1)-V(iy,ix) ) / dx(iy,ix) ; |
103 |
dUdy = ( U(iy+1,ix)-U(iy,ix) ) / dy(iy,ix) ; |
104 |
ZETA(iz,iy,ix) = dVdx - dUdy; |
105 |
end %for ix |
106 |
end %for iy |
107 |
|
108 |
end %for iz |
109 |
|
110 |
%%%%%%%%%%%%%% |
111 |
%% Netcdf record: |
112 |
|
113 |
% General informations: |
114 |
netfil = strcat('ZETA','.',netcdf_domain,'.',netcdf_suff); |
115 |
units = '1/s'; |
116 |
ncid = 'ZETA'; |
117 |
longname = 'Vertical Component of the Relative Vorticity'; |
118 |
uniquename = 'vertical_relative_vorticity'; |
119 |
|
120 |
% Open output file: |
121 |
nc = netcdf(strcat(pathname,sla,netfil),'clobber'); |
122 |
|
123 |
% Define axis: |
124 |
nc('X') = nx; |
125 |
nc('Y') = ny; |
126 |
nc('Z') = nz; |
127 |
|
128 |
nc{'X'} = 'X'; |
129 |
nc{'Y'} = 'Y'; |
130 |
nc{'Z'} = 'Z'; |
131 |
|
132 |
nc{'X'} = ncfloat('X'); |
133 |
nc{'X'}.uniquename = ncchar('X'); |
134 |
nc{'X'}.long_name = ncchar('longitude'); |
135 |
nc{'X'}.gridtype = nclong(0); |
136 |
nc{'X'}.units = ncchar('degrees_east'); |
137 |
nc{'X'}(:) = ZETA_lon; |
138 |
|
139 |
nc{'Y'} = ncfloat('Y'); |
140 |
nc{'Y'}.uniquename = ncchar('Y'); |
141 |
nc{'Y'}.long_name = ncchar('latitude'); |
142 |
nc{'Y'}.gridtype = nclong(0); |
143 |
nc{'Y'}.units = ncchar('degrees_north'); |
144 |
nc{'Y'}(:) = ZETA_lat; |
145 |
|
146 |
nc{'Z'} = ncfloat('Z'); |
147 |
nc{'Z'}.uniquename = ncchar('Z'); |
148 |
nc{'Z'}.long_name = ncchar('depth'); |
149 |
nc{'Z'}.gridtype = nclong(0); |
150 |
nc{'Z'}.units = ncchar('m'); |
151 |
nc{'Z'}(:) = Udpt; |
152 |
|
153 |
% And main field: |
154 |
nc{ncid} = ncfloat('Z', 'Y', 'X'); |
155 |
nc{ncid}.units = ncchar(units); |
156 |
nc{ncid}.missing_value = ncfloat(NaN); |
157 |
nc{ncid}.FillValue_ = ncfloat(NaN); |
158 |
nc{ncid}.longname = ncchar(longname); |
159 |
nc{ncid}.uniquename = ncchar(uniquename); |
160 |
nc{ncid}(:,:,:) = ZETA; |
161 |
|
162 |
nc=close(nc); |
163 |
|
164 |
clear x y z U V dx dy nx ny nz DVdx dUdy |
165 |
|
166 |
end %if compute ZETA |
167 |
|
168 |
|
169 |
%%%%%%%%%%%%%%%%%%%%%%%%% |
170 |
% HORIZONTAL COMPONENTS % |
171 |
%%%%%%%%%%%%%%%%%%%%%%%%% |
172 |
if toshow, disp('') |
173 |
disp('Now compute horizontal components of relative vorticity ...'); end |
174 |
|
175 |
% U and V are defined on the same Z grid. |
176 |
|
177 |
%%%%%%%%%%%%%% |
178 |
%% Dimensions of OMEGA x and y fields: |
179 |
if toshow,disp('Dim'),end |
180 |
O_nx = [length(Vlon) length(Ulon)]; |
181 |
O_ny = [length(Vlat) length(Ulat)]; |
182 |
O_nz = length(Udpt) - 1; % Idem Vdpt |
183 |
|
184 |
%%%%%%%%%%%%%% |
185 |
%% Pre-allocations: |
186 |
if toshow,disp('Pre-allocate'),end |
187 |
Ox = zeros(O_nz,O_ny(1),O_nx(1)).*NaN; |
188 |
Oy = zeros(O_nz,O_ny(2),O_nx(2)).*NaN; |
189 |
|
190 |
%%%%%%%%%%%%%% |
191 |
%% Horizontal components: |
192 |
|
193 |
%% Vertical grid differences: |
194 |
dZ = diff(Udpt); |
195 |
Odpt = Udpt(1:O_nz) + dZ/2; |
196 |
|
197 |
%% Zonal component of OMEGA: |
198 |
if toshow,disp('Zonal direction ...'); end |
199 |
[a dZ_3D c] = meshgrid(Vlat,dZ,Vlon); clear a c |
200 |
V = ncV{4}(:,:,:); |
201 |
Ox = - ( V(2:O_nz+1,:,:) - V(1:O_nz,:,:) ) ./ dZ_3D; |
202 |
clear V dZ_3D % For memory use |
203 |
|
204 |
%% Meridional component of OMEGA: |
205 |
if toshow,disp('Meridional direction ...'); end |
206 |
[a dZ_3D c] = meshgrid(Ulat,dZ,Ulon); clear a c |
207 |
U = ncU{4}(:,:,:); |
208 |
Oy = ( U(2:O_nz+1,:,:) - U(1:O_nz,:,:) ) ./ dZ_3D; |
209 |
clear U dZ_3D % For memory use |
210 |
|
211 |
clear dZ |
212 |
|
213 |
%%%%%%%%%%%%%% |
214 |
%% Record Zonal component: |
215 |
if toshow,disp('Records ...'); end |
216 |
|
217 |
% General informations: |
218 |
netfil = strcat('OMEGAX','.',netcdf_domain,'.',netcdf_suff); |
219 |
units = '1/s'; |
220 |
ncid = 'OMEGAX'; |
221 |
longname = 'Zonal Component of the Relative Vorticity'; |
222 |
uniquename = 'zonal_relative_vorticity'; |
223 |
|
224 |
% Open output file: |
225 |
nc = netcdf(strcat(pathname,sla,netfil),'clobber'); |
226 |
|
227 |
% Define axis: |
228 |
nc('X') = O_nx(1); |
229 |
nc('Y') = O_ny(1); |
230 |
nc('Z') = O_nz; |
231 |
|
232 |
nc{'X'} = 'X'; |
233 |
nc{'Y'} = 'Y'; |
234 |
nc{'Z'} = 'Z'; |
235 |
|
236 |
nc{'X'} = ncfloat('X'); |
237 |
nc{'X'}.uniquename = ncchar('X'); |
238 |
nc{'X'}.long_name = ncchar('longitude'); |
239 |
nc{'X'}.gridtype = nclong(0); |
240 |
nc{'X'}.units = ncchar('degrees_east'); |
241 |
nc{'X'}(:) = Vlon; |
242 |
|
243 |
nc{'Y'} = ncfloat('Y'); |
244 |
nc{'Y'}.uniquename = ncchar('Y'); |
245 |
nc{'Y'}.long_name = ncchar('latitude'); |
246 |
nc{'Y'}.gridtype = nclong(0); |
247 |
nc{'Y'}.units = ncchar('degrees_north'); |
248 |
nc{'Y'}(:) = Vlat; |
249 |
|
250 |
nc{'Z'} = ncfloat('Z'); |
251 |
nc{'Z'}.uniquename = ncchar('Z'); |
252 |
nc{'Z'}.long_name = ncchar('depth'); |
253 |
nc{'Z'}.gridtype = nclong(0); |
254 |
nc{'Z'}.units = ncchar('m'); |
255 |
nc{'Z'}(:) = Odpt; |
256 |
|
257 |
% And main field: |
258 |
nc{ncid} = ncfloat('Z', 'Y', 'X'); |
259 |
nc{ncid}.units = ncchar(units); |
260 |
nc{ncid}.missing_value = ncfloat(NaN); |
261 |
nc{ncid}.FillValue_ = ncfloat(NaN); |
262 |
nc{ncid}.longname = ncchar(longname); |
263 |
nc{ncid}.uniquename = ncchar(uniquename); |
264 |
nc{ncid}(:,:,:) = Ox; |
265 |
|
266 |
nc=close(nc); |
267 |
|
268 |
%%%%%%%%%%%%%% |
269 |
%% Record Meridional component: |
270 |
% General informations: |
271 |
netfil = strcat('OMEGAY','.',netcdf_domain,'.',netcdf_suff); |
272 |
units = '1/s'; |
273 |
ncid = 'OMEGAY'; |
274 |
longname = 'Meridional Component of the Relative Vorticity'; |
275 |
uniquename = 'meridional_relative_vorticity'; |
276 |
|
277 |
% Open output file: |
278 |
nc = netcdf(strcat(pathname,sla,netfil),'clobber'); |
279 |
|
280 |
% Define axis: |
281 |
nc('X') = O_nx(2); |
282 |
nc('Y') = O_ny(2); |
283 |
nc('Z') = O_nz; |
284 |
|
285 |
nc{'X'} = 'X'; |
286 |
nc{'Y'} = 'Y'; |
287 |
nc{'Z'} = 'Z'; |
288 |
|
289 |
nc{'X'} = ncfloat('X'); |
290 |
nc{'X'}.uniquename = ncchar('X'); |
291 |
nc{'X'}.long_name = ncchar('longitude'); |
292 |
nc{'X'}.gridtype = nclong(0); |
293 |
nc{'X'}.units = ncchar('degrees_east'); |
294 |
nc{'X'}(:) = Ulon; |
295 |
|
296 |
nc{'Y'} = ncfloat('Y'); |
297 |
nc{'Y'}.uniquename = ncchar('Y'); |
298 |
nc{'Y'}.long_name = ncchar('latitude'); |
299 |
nc{'Y'}.gridtype = nclong(0); |
300 |
nc{'Y'}.units = ncchar('degrees_north'); |
301 |
nc{'Y'}(:) = Ulat; |
302 |
|
303 |
nc{'Z'} = ncfloat('Z'); |
304 |
nc{'Z'}.uniquename = ncchar('Z'); |
305 |
nc{'Z'}.long_name = ncchar('depth'); |
306 |
nc{'Z'}.gridtype = nclong(0); |
307 |
nc{'Z'}.units = ncchar('m'); |
308 |
nc{'Z'}(:) = Odpt; |
309 |
|
310 |
% And main field: |
311 |
nc{ncid} = ncfloat('Z', 'Y', 'X'); |
312 |
nc{ncid}.units = ncchar(units); |
313 |
nc{ncid}.missing_value = ncfloat(NaN); |
314 |
nc{ncid}.FillValue_ = ncfloat(NaN); |
315 |
nc{ncid}.longname = ncchar(longname); |
316 |
nc{ncid}.uniquename = ncchar(uniquename); |
317 |
nc{ncid}(:,:,:) = Oy; |
318 |
|
319 |
nc=close(nc); |
320 |
|