MITgcm_contrib/arnaud_matlab/calc_vort.m

%
%  function [XZ,YZ,ZETAr,ZETAp] = calc_vort(U,V,DX,DY,Ymin)
%
% Computes `vertical' component of 
% the planetary and relative vorticity
% ZETAp(XZ,YZ) and ZETAr(XZ,YZ) from 
% (U,V) on a C-grid (XU,YU) (XV,YV)
% 
% Ymin is the southern latitude (negative, in degree)
% DX is the longitudinal resolution (in degree)
% DY is the latitudinal resolution (in degree)
%
% YZ has NY-1 component. i.e. ZETAr and ZETAp 
% are not computed on the southern boundary
% (at Ymin where V=0)
%
% (c) acz, Nov. 2002


  function [XZ,YZ,ZETAr,ZETAp] = calc_vort(U,V,DX,DY,Ymin)

% C-grid
%
  [NX NY] = size(U); %or V
  XU = [0:DX:(DX*NX-DX)];
  XV = XU + DX/2;
  YU = [(Ymin+DY/2):DY:(-Ymin-DY/2)]; 
  YV = [Ymin:DY:-Ymin-DY];

% Calculate Vorticity
%
  ZETAr = NaN * ones(NX,NY-1);
  ZETAp = NaN * ones(NX,NY-1);
  XZ = XU;
  YZ = YV(2:end);

  RADIUS = 6371 * 1000;
  OMEGA = 2 * pi / (24 * 3600);

  for j = 1:NY-1

    for i = 2:NX
      dy = RADIUS * DY * pi/180;
      dxN = RADIUS * cos(YU(j+1)*pi/180) * DX * pi/180;
      dxS = RADIUS * cos(YU(j)*pi/180) * DX * pi/180;
      h = sqrt( dy^2 - 0.25*(dxS-dxN)^2 );
      area = 0.5 * h * (dxS + dxN); %Formule du Trapeze 
      ZETAr(i,j) = - (dy*V(i-1,j+1) + dxN*U(i,j+1) - dy*V(i,j+1) - dxS*U(i,j)) / area;
      ZETAp(i,j) = 2 * OMEGA * sin(YZ(j)*pi/180);
    end

    ZETAr(1,j) = - (dy*V(NX,j+1) + dxN*U(1,j+1) - dy*V(1,j+1) - dxS*U(1,j)) / area;
    ZETAp(1,j) = 2 * OMEGA * sin(YZ(j)*pi/180);

  end

1	edhill	1.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) = - (dyV(i-1,j+1) + dxNU(i,j+1) - dyV(i,j+1) - dxSU(i,j)) / area;
49			ZETAp(i,j) = 2 * OMEGA * sin(YZ(j)*pi/180);
50			end
51
52			ZETAr(1,j) = - (dyV(NX,j+1) + dxNU(1,j+1) - dyV(1,j+1) - dxSU(1,j)) / area;
53			ZETAp(1,j) = 2 * OMEGA * sin(YZ(j)*pi/180);
54
55			end
56