function [H2,Hx2,Hy2]=xyrecur_fv(He,Hx,Hy,niter)
%
% Returns the gridded data on a grid of half the size
% with the data "averaged" to the cell-centers using
% a "Finite Volume" algorithm (see Adcroft, JPO 2000).
%
% e.g.
%      [H2,Hx2,Hy2]=xyrecur_fv(He,Hx,Hy,3);

debuglevels=0;

%H2=He;
%Hx2=Hx;
%Hy2=Hy;

% Check that Hx and Hy don't violate basic rules
check_HxHy(He,Hx,Hy)

% Recursively iterate
for nit=1:niter,

clear H2 Hx2 Hy2
[nx,ny,nt]=size(He);
disp(sprintf(' Recursion level %i (%i,%i) -> (%i,%i) ... ',nit,nx,ny,nx/2,ny/2));

%ii=1:2:nx-1;
%jj=1:2:ny-1;
%ip=2:2:nx;
%im=[nx 2:2:nx-2];
%jp=2:2:ny;
%jm=[2 2:2:ny-2];
%I=1:nx/2; IM=I([end 1:end-1]);
%J=1:ny/2; JM=J([end 1:end-1]);

for J=1:ny/2; jj=1+2*(J-1);
  jp=mymod(jj+1,ny);jm=mymod(jj-1,ny);JM=mymod(J-1,ny/2);
for I=1:nx/2; ii=1+2*(I-1);
  ip=mymod(ii+1,nx);im=mymod(ii-1,nx);IM=mymod(I-1,nx/2);

% Four-point average in X-Y
%new? H2=(He(ii,jj,:)+He(ii+1,jj,:)+He(ii,jj+1,:)+He(ii+1,jj+1,:))/4;
H2(I,J,:)=(He(ii,jj,:)+He(ip,jj,:)+He(ii,jp,:)+He(ip,jp,:))/4;

% Two-point average in X
%new? Hx2=(Hx(ii,jj,:)+Hx(ii,jj+1,:))/2;
Hx2(I,J,:)=(Hx(ii,jj,:)+Hx(ii,jp,:))/2;

% Two-point average in Y
%new? Hy2=(Hy(ii,jj,:)+Hy(ii+1,jj,:))/2;
Hy2(I,J,:)=(Hy(ii,jj,:)+Hy(ip,jj,:))/2;

end;end;
for J=1:ny/2; jj=1+2*(J-1);
  jp=mymod(jj+1,ny),jm=mymod(jj-1,ny),JM=mymod(J-1,ny/2),
for I=1:nx/2; ii=1+2*(I-1);
  ip=mymod(ii+1,nx),im=mymod(ii-1,nx),IM=mymod(I-1,nx/2),

% Blocking from either side?
Hx2e=(Hx(ip,jj,:)+Hx(ip,jp,:))/2;
Hx2w=(Hx(im,jj,:)+Hx(im,jp,:))/2;
Hx2(I,J,:)=max( max( Hx2(I,J,:), Hx2e ) , Hx2w);

% Blocking from either side?
Hy2n=(Hy(ii,jp,:)+Hy(ip,jp,:))/2;
Hy2s=(Hy(ii,jm,:)+Hy(ip,jm,:))/2;
Hy2(I,J,:)=max( max( Hy2(I,J,:), Hy2n ) , Hy2s);

% Make Hx2,Hy2 consistent with coarse H2
%Hx2=max( max(Hx2, H2), H2([end 1:end-1],:,:));
%Hy2=max( max(Hy2, H2), H2(:,[end 1:end-1],:));

% Fill-up cell-centers where all sides are higher
%H2sides=min( Hx2, Hx2([2:end 1],:,:) );
%H2sides=min( H2sides, Hy2 );
%H2sides=min( H2sides, Hy2(:,[2:end end],:) );
%H2=max( H2, H2sides );

% Fill up if barriers intersects coarse cell
%%Hx2mid=(Hx(ip,jj,:)+Hx(ip,jp,:))/2;
%%Hy2mid=(Hy(ii,jp,:)+Hy(ip,jp,:))/2;
%Hx2mid=min(Hx(ip,jj,:),Hx(ip,jp,:));
%Hy2mid=min(Hy(ii,jp,:),Hy(ip,jp,:));
%H2mid=max(max(H2, Hy2mid),Hx2mid);
%ll=find(H2mid>H2sides);
%H2(ll)=H2mid(ll);

% Make Hx2,Hy2 consistent with coarse H2
%Hx2=max( max(Hx2, H2), H2([end 1:end-1],:,:));
%Hy2=max( max(Hy2, H2), H2(:,[end 1:end-1],:));
Hx2(I,J,:)=max( max(Hx2(I,J,:), H2(I,J,:)), H2(IM,J,:));
Hy2(I,J,:)=max( max(Hy2(I,J,:), H2(I,J,:)), H2(I,JM,:));

end;end;

% Check that Hx2 and Hy2 don't violate basic rules
check_HxHy(H2,Hx2,Hy2)

% Re-cycle data for another iteration
He=H2;
Hx=Hx2;
Hy=Hy2;

if debuglevels
 eval(sprintf('save LEVEL%2.2i He Hx Hy',nit));
end

end

% Fill-up cell-centers where all sides are higher
H2sides=min( Hx2, Hx2([2:end 1],:,:) );
H2sides=min( H2sides, Hy2 );
H2sides=min( H2sides, Hy2(:,[2:end end],:) );
H2=max( H2, H2sides );



function []=QQQQQcheck_HxHy(H2,Hx2,Hy2)
% Check that Hx2 and Hy2 don't violate basic rules
% This should *not* ever happen
return; % spead up process
[Hxcc,Hycc]=gen_hxhy(H2,x2,y2);
ic=find(Hx2<Hxcc-0.001);
if prod(size(ic)) ~= 0
 [Hx2(ic)/1e3 Hxcc(ic)/1e3 Hx2(ic)-Hxcc(ic)]
 error('The X cross-sectional area turned out inconsistently')
end
jc=find(Hy2<Hycc-0.001);
if prod(size(jc)) ~= 0
 [Hy2(jc)/1e3 Hycc(jc)/1e3 Hy2(jc)-Hycc(jc)]
 error('The Y cross-sectional area turned out inconsistently')
end


function [n] = mymod(n,N)
n=mod(n+N-1,N)+1;