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