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