/[MITgcm]/MITgcm_contrib/arnaud_matlab/calc_vort.m
ViewVC logotype

Contents of /MITgcm_contrib/arnaud_matlab/calc_vort.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Mon Aug 25 16:00:59 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
Initial checkin of Arnaud's MatLAB diagnostics

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

  ViewVC Help
Powered by ViewVC 1.1.22