1 |
% |
2 |
% function [XZ,YZ,ZETAr,ZETAp] = calc_vort(U,V,DX,DY,Ymin) |
3 |
% |
4 |
% Computes `vertical' component of |
5 |
% the planetary and relative vorticity |
6 |
% ZETAp(XZ,YZ) and ZETAr(XZ,YZ) from |
7 |
% (U,V) on a C-grid (XU,YU) (XV,YV) |
8 |
% |
9 |
% Ymin is the southern latitude (negative, in degree) |
10 |
% DX is the longitudinal resolution (in degree) |
11 |
% DY is the latitudinal resolution (in degree) |
12 |
% |
13 |
% YZ has NY-1 component. i.e. ZETAr and ZETAp |
14 |
% are not computed on the southern boundary |
15 |
% (at Ymin where V=0) |
16 |
% |
17 |
% (c) acz, Nov. 2002 |
18 |
|
19 |
|
20 |
function [XZ,YZ,ZETAr,ZETAp] = calc_vort(U,V,DX,DY,Ymin) |
21 |
|
22 |
% C-grid |
23 |
% |
24 |
[NX NY] = size(U); %or V |
25 |
XU = [0:DX:(DX*NX-DX)]; |
26 |
XV = XU + DX/2; |
27 |
YU = [(Ymin+DY/2):DY:(-Ymin-DY/2)]; |
28 |
YV = [Ymin:DY:-Ymin-DY]; |
29 |
|
30 |
% Calculate Vorticity |
31 |
% |
32 |
ZETAr = NaN * ones(NX,NY-1); |
33 |
ZETAp = NaN * ones(NX,NY-1); |
34 |
XZ = XU; |
35 |
YZ = YV(2:end); |
36 |
|
37 |
RADIUS = 6371 * 1000; |
38 |
OMEGA = 2 * pi / (24 * 3600); |
39 |
|
40 |
for j = 1:NY-1 |
41 |
|
42 |
for i = 2:NX |
43 |
dy = RADIUS * DY * pi/180; |
44 |
dxN = RADIUS * cos(YU(j+1)*pi/180) * DX * pi/180; |
45 |
dxS = RADIUS * cos(YU(j)*pi/180) * DX * pi/180; |
46 |
h = sqrt( dy^2 - 0.25*(dxS-dxN)^2 ); |
47 |
area = 0.5 * h * (dxS + dxN); %Formule du Trapeze |
48 |
ZETAr(i,j) = - (dy*V(i-1,j+1) + dxN*U(i,j+1) - dy*V(i,j+1) - dxS*U(i,j)) / area; |
49 |
ZETAp(i,j) = 2 * OMEGA * sin(YZ(j)*pi/180); |
50 |
end |
51 |
|
52 |
ZETAr(1,j) = - (dy*V(NX,j+1) + dxN*U(1,j+1) - dy*V(1,j+1) - dxS*U(1,j)) / area; |
53 |
ZETAp(1,j) = 2 * OMEGA * sin(YZ(j)*pi/180); |
54 |
|
55 |
end |
56 |
|