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