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