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

% Check that Hx and Hy don't violate basic rules

debuglevels=1;
useavg=0;
%%% keyboard

if ~isempty(Hx)
 disp('Hx,Hy defined on entry so checking consistency...')
 check_HxHy(He,Hx,Hy)
end

% Recursively iterate
for nit=1:niter,

hxhyempty=isempty(Hx);

[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];
iM=ii([end 1:end-1]);
jp=2:2:ny;
jm=[2 2:2:ny-2];
jM=jj([end 1:end-1]);

% Adjust He for angled barriers in coarse cell (H2)
disp('   Adjust fine grid for angled barriers (NE)')
if hxhyempty
 Hcrnr=min(max(He(ii,jp,:),He(ip,jp,:)),max(He(ip,jj,:),He(ip,jp,:)));
else
 Hcrnr=min(Hx(ip,jp,:),Hy(ip,jp,:));
end
[li,lj]=find( Hcrnr>He(ip,jp) );
if ~isempty(li) & ~isempty(lj)
He( sub2ind(size(He),2*li,2*lj) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
end
clear Hcrnr li lj
disp('   Adjust fine grid for angled barriers (SE)')
if hxhyempty
 Hcrnr=min(max(He(ii,jj,:),He(ip,jj,:)),max(He(ip,jj,:),He(ip,jp,:)));
else
 Hcrnr=min(Hx(ip,jj,:),Hy(ip,jp,:));
end
[li,lj]=find( Hcrnr>He(ip,jj) );
if ~isempty(li) & ~isempty(lj)
He( sub2ind(size(He),2*li,2*lj-1) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
end
clear Hcrnr li lj
disp('   Adjust fine grid for angled barriers (SW)')
if hxhyempty
 Hcrnr=min(max(He(ii,jj,:),He(ip,jj,:)),max(He(ii,jj,:),He(ii,jp,:)));
else
 Hcrnr=min(Hx(ip,jj,:),Hy(ii,jp,:));
end
[li,lj]=find( Hcrnr>He(ii,jj) );
if ~isempty(li) & ~isempty(lj)
He( sub2ind(size(He),2*li-1,2*lj-1) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
end
clear Hcrnr li lj
disp('   Adjust fine grid for angled barriers (NW)')
if hxhyempty
 Hcrnr=min(max(He(ii,jp,:),He(ip,jp,:)),max(He(ii,jj,:),He(ii,jp,:)));
else
 Hcrnr=min(Hx(ip,jp,:),Hy(ii,jp,:));
end
[li,lj]=find( Hcrnr>He(ii,jp) );
if ~isempty(li) & ~isempty(lj)
He( sub2ind(size(He),2*li-1,2*lj) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
end
clear Hcrnr li lj

if ~hxhyempty
 disp('   Regenerating Hx,Hy after adjusting He')
 [Nx,Ny]=gen_hxhy(He);
 Hx=max(Hx,Nx);
 Hy=max(Hy,Ny);
 clear Nx Ny
end

% Four-point average in X-Y
disp('   Four-point average for H2')
if useavg
 H2=(He(ii,jj,:)+He(ip,jj,:)+He(ii,jp,:)+He(ip,jp,:))/4;
else
 H2=min(min(He(ii,jj,:),He(ip,jj,:)),min(He(ii,jp,:),He(ip,jp,:)));
end

% Two-point average in X
disp('   Line average for Hx2')
if hxhyempty
 Hx2=(max(He(ii,jj,:),He(im,jj,:))+max(He(ii,jp,:),He(im,jp,:)))/2;
else
 if useavg
  Hx2=(Hx(ii,jj,:)+Hx(ii,jp,:))/2;
 else
  Hx2=min(Hx(ii,jj,:),Hx(ii,jp,:));
 end
end

% Two-point average in Y
disp('   Line average for Hy2')
if hxhyempty
 Hy2=(max(He(ii,jj,:),He(ii,jm,:))+max(He(ip,jj,:),He(ip,jm,:)))/2;
else
 if useavg
  Hy2=(Hy(ii,jj,:)+Hy(ip,jj,:))/2;
 else
  Hy2=min(Hy(ii,jj,:),Hy(ip,jj,:));
 end
end

% Straight barriers intersecting coarse cell
disp('   Filling cell due to staight zonal barriers')
if hxhyempty
 Hx2mid=min( max(He(ii,jj,:),He(ip,jj,:)) , max(He(ii,jp,:),He(ip,jp,:)) );
else
 Hx2mid=min(Hx(ip,jj,:),Hx(ip,jp,:));
end
Hx2sides=max( Hx2, Hx2([2:end 1],:,:) );
ll=find(Hx2mid>Hx2sides & Hx2mid>H2);
H2(ll)=Hx2mid(ll);
disp('   Filling cell due to staight meridional barriers')
if hxhyempty
 Hy2mid=min( max(He(ii,jj,:),He(ii,jp,:)) , max(He(ip,jj,:),He(ip,jp,:)) );
else
 Hy2mid=min(Hy(ii,jp,:),Hy(ip,jp,:));
end
Hy2sides=max( Hy2, Hy2(:,[2:end 1],:) );
ll=find(Hy2mid>Hy2sides & Hy2mid>H2);
H2(ll)=Hy2mid(ll);

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

% 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
disp('Filling holes in H2 (bowls)')
H2sides=min( Hx2, Hx2([2:end 1],:,:) );
H2sides=min( H2sides, Hy2 );
H2sides=min( H2sides, Hy2(:,[2:end end],:) );
%H2=max( H2, H2sides );