MITgcm_contrib/arnaud_matlab/calc_hadv.m

%
%  function [XT,YT,TADV] = calc_hadv(T,U,V,DX,DY,Ymin)
%
% Computes horizontal advection of scalar T 
% on the XT, YT (physics) grid at a given vertical level.
% NB: (T,U,V) are 2-D fields
%
% (U,V) on a C-grid (XU,YU) (XV,YV)
% NB: The way it is computed is consistent with
% the flux form used by MIT-GCM (i.e. if
% T DIV is added one recovers the flux form)
% 
% Ymin is the southern latitude (negative, in degree)
% DX is the longitudinal resolution (in degree)
% DY is the latitudinal resolution (in degree)
%
% (c) acz, Jul. 2003


  function [XT,YT,TADV] = calc_hadv(T,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];
  XT = XV; YT = YU;

% Constants
  RADIUS = 6371 * 1000;
  DYG = RADIUS * DY * pi/180;
  DXG = RADIUS * DX * pi/180;

% Calculate zonal advection on U-grid
  advu = zeros(NX,NY);
  for i = 1:NX-1
    advu(i+1,:) = U(i+1,:) .* (T(i+1,:)-T(i,:));
  end
  advu(1,:) = U(1,:) .* (T(1,:)-T(NX,:));

% Average advu on T-grid
  advuTG = zeros(NX,NY);
  AG = cos(YT*pi/180) * DYG * DXG;
  for i = 1:NX-1
    advuTG(i,:) = DYG*( advu(i,:)+advu(i+1,:) ) ./ (2*AG);
  end
  advuTG(NX,:) = DYG*( advu(NX,:)+advu(1,:) ) ./ (2*AG);

% Calculate meridional advection on V-grid
  advv = zeros(NX,NY); %note advv(:,1) = 0 because v(:,1)=0
  for j = 2:NY
    advv(:,j) = V(:,j) .* (T(:,j)-T(:,j-1));
  end

% Average advv on T-grid
  advvTG = zeros(NX,NY);
  for j = 1:NY-1
    AG(j) = DYG * DXG * ( cos(YV(j+1)*pi/180)+cos(YV(j)*pi/180) )/2;
  end
  AG(NY) = DYG * DXG * ( cos((YV(NY)+DY)*pi/180)+cos(YV(NY)*pi/180) )/2;
  DDXG = DXG * cos((YV+DY)*pi/180);
  for j = 1:NY-1
    advvTG(:,j) = DDXG(j)*( advv(:,j)+advv(:,j+1) ) ./ (2*AG(j));
  end
  advvTG(:,NY) = DDXG(NY)*( advv(:,NY)+0 ) ./ (2*AG(NY));

% Horizontal advection
  TADV = advuTG + advvTG;

1	edhill	1.1	%
2			% function [XT,YT,TADV] = calc_hadv(T,U,V,DX,DY,Ymin)
3			%
4			% Computes horizontal advection of scalar T
5			% on the XT, YT (physics) grid at a given vertical level.
6			% NB: (T,U,V) are 2-D fields
7			%
8			% (U,V) on a C-grid (XU,YU) (XV,YV)
9			% NB: The way it is computed is consistent with
10			% the flux form used by MIT-GCM (i.e. if
11			% T DIV is added one recovers the flux form)
12			%
13			% Ymin is the southern latitude (negative, in degree)
14			% DX is the longitudinal resolution (in degree)
15			% DY is the latitudinal resolution (in degree)
16			%
17			% (c) acz, Jul. 2003
18
19
20			function [XT,YT,TADV] = calc_hadv(T,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			XT = XV; YT = YU;
30
31			% Constants
32			RADIUS = 6371 * 1000;
33			DYG = RADIUS * DY * pi/180;
34			DXG = RADIUS * DX * pi/180;
35
36			% Calculate zonal advection on U-grid
37			advu = zeros(NX,NY);
38			for i = 1:NX-1
39			advu(i+1,:) = U(i+1,:) .* (T(i+1,:)-T(i,:));
40			end
41			advu(1,:) = U(1,:) .* (T(1,:)-T(NX,:));
42
43			% Average advu on T-grid
44			advuTG = zeros(NX,NY);
45			AG = cos(YTpi/180) DYG * DXG;
46			for i = 1:NX-1
47			advuTG(i,:) = DYG( advu(i,:)+advu(i+1,:) ) ./ (2AG);
48			end
49			advuTG(NX,:) = DYG( advu(NX,:)+advu(1,:) ) ./ (2AG);
50
51			% Calculate meridional advection on V-grid
52			advv = zeros(NX,NY); %note advv(:,1) = 0 because v(:,1)=0
53			for j = 2:NY
54			advv(:,j) = V(:,j) .* (T(:,j)-T(:,j-1));
55			end
56
57			% Average advv on T-grid
58			advvTG = zeros(NX,NY);
59			for j = 1:NY-1
60			AG(j) = DYG * DXG * ( cos(YV(j+1)pi/180)+cos(YV(j)pi/180) )/2;
61			end
62			AG(NY) = DYG * DXG * ( cos((YV(NY)+DY)pi/180)+cos(YV(NY)pi/180) )/2;
63			DDXG = DXG * cos((YV+DY)*pi/180);
64			for j = 1:NY-1
65			advvTG(:,j) = DDXG(j)( advv(:,j)+advv(:,j+1) ) ./ (2AG(j));
66			end
67			advvTG(:,NY) = DDXG(NY)( advv(:,NY)+0 ) ./ (2AG(NY));
68
69			% Horizontal advection
70			TADV = advuTG + advvTG;
71