%
%  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