mlosch/interp_llc/xyrecur_fv.m

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 );
1	mlosch	1.1	function [H2,Hx2,Hy2]=xyrecur_fv(He,Hx,Hy,niter)
2			%
3			% Returns the gridded data on a grid of half the size
4			% with the data "averaged" to the cell-centers using
5			% a "Finite Volume" algorithm (see Adcroft, JPO 2000).
6			%
7			% e.g.
8			% [H2,Hx2,Hy2]=xyrecur_fv(He,Hx,Hy,3);
9
10			% Check that Hx and Hy don't violate basic rules
11
12			debuglevels=1;
13			useavg=0;
14			%%% keyboard
15
16			if ~isempty(Hx)
17			disp('Hx,Hy defined on entry so checking consistency...')
18			check_HxHy(He,Hx,Hy)
19			end
20
21			% Recursively iterate
22			for nit=1:niter,
23
24			hxhyempty=isempty(Hx);
25
26			[nx,ny,nt]=size(He);
27			disp(sprintf(' Recursion level %i (%i,%i) -> (%i,%i) ... ',nit,nx,ny,nx/2,ny/2));
28
29			ii=1:2:nx-1;
30			jj=1:2:ny-1;
31			ip=2:2:nx;
32			im=[nx 2:2:nx-2];
33			iM=ii([end 1:end-1]);
34			jp=2:2:ny;
35			jm=[2 2:2:ny-2];
36			jM=jj([end 1:end-1]);
37
38			% Adjust He for angled barriers in coarse cell (H2)
39			disp(' Adjust fine grid for angled barriers (NE)')
40			if hxhyempty
41			Hcrnr=min(max(He(ii,jp,:),He(ip,jp,:)),max(He(ip,jj,:),He(ip,jp,:)));
42			else
43			Hcrnr=min(Hx(ip,jp,:),Hy(ip,jp,:));
44			end
45			[li,lj]=find( Hcrnr>He(ip,jp) );
46			if ~isempty(li) & ~isempty(lj)
47			He( sub2ind(size(He),2li,2lj) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
48			end
49			clear Hcrnr li lj
50			disp(' Adjust fine grid for angled barriers (SE)')
51			if hxhyempty
52			Hcrnr=min(max(He(ii,jj,:),He(ip,jj,:)),max(He(ip,jj,:),He(ip,jp,:)));
53			else
54			Hcrnr=min(Hx(ip,jj,:),Hy(ip,jp,:));
55			end
56			[li,lj]=find( Hcrnr>He(ip,jj) );
57			if ~isempty(li) & ~isempty(lj)
58			He( sub2ind(size(He),2li,2lj-1) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
59			end
60			clear Hcrnr li lj
61			disp(' Adjust fine grid for angled barriers (SW)')
62			if hxhyempty
63			Hcrnr=min(max(He(ii,jj,:),He(ip,jj,:)),max(He(ii,jj,:),He(ii,jp,:)));
64			else
65			Hcrnr=min(Hx(ip,jj,:),Hy(ii,jp,:));
66			end
67			[li,lj]=find( Hcrnr>He(ii,jj) );
68			if ~isempty(li) & ~isempty(lj)
69			He( sub2ind(size(He),2li-1,2lj-1) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
70			end
71			clear Hcrnr li lj
72			disp(' Adjust fine grid for angled barriers (NW)')
73			if hxhyempty
74			Hcrnr=min(max(He(ii,jp,:),He(ip,jp,:)),max(He(ii,jj,:),He(ii,jp,:)));
75			else
76			Hcrnr=min(Hx(ip,jp,:),Hy(ii,jp,:));
77			end
78			[li,lj]=find( Hcrnr>He(ii,jp) );
79			if ~isempty(li) & ~isempty(lj)
80			He( sub2ind(size(He),2li-1,2lj) )=Hcrnr( sub2ind(size(Hcrnr),li,lj) );
81			end
82			clear Hcrnr li lj
83
84			if ~hxhyempty
85			disp(' Regenerating Hx,Hy after adjusting He')
86			[Nx,Ny]=gen_hxhy(He);
87			Hx=max(Hx,Nx);
88			Hy=max(Hy,Ny);
89			clear Nx Ny
90			end
91
92			% Four-point average in X-Y
93			disp(' Four-point average for H2')
94			if useavg
95			H2=(He(ii,jj,:)+He(ip,jj,:)+He(ii,jp,:)+He(ip,jp,:))/4;
96			else
97			H2=min(min(He(ii,jj,:),He(ip,jj,:)),min(He(ii,jp,:),He(ip,jp,:)));
98			end
99
100			% Two-point average in X
101			disp(' Line average for Hx2')
102			if hxhyempty
103			Hx2=(max(He(ii,jj,:),He(im,jj,:))+max(He(ii,jp,:),He(im,jp,:)))/2;
104			else
105			if useavg
106			Hx2=(Hx(ii,jj,:)+Hx(ii,jp,:))/2;
107			else
108			Hx2=min(Hx(ii,jj,:),Hx(ii,jp,:));
109			end
110			end
111
112			% Two-point average in Y
113			disp(' Line average for Hy2')
114			if hxhyempty
115			Hy2=(max(He(ii,jj,:),He(ii,jm,:))+max(He(ip,jj,:),He(ip,jm,:)))/2;
116			else
117			if useavg
118			Hy2=(Hy(ii,jj,:)+Hy(ip,jj,:))/2;
119			else
120			Hy2=min(Hy(ii,jj,:),Hy(ip,jj,:));
121			end
122			end
123
124			% Straight barriers intersecting coarse cell
125			disp(' Filling cell due to staight zonal barriers')
126			if hxhyempty
127			Hx2mid=min( max(He(ii,jj,:),He(ip,jj,:)) , max(He(ii,jp,:),He(ip,jp,:)) );
128			else
129			Hx2mid=min(Hx(ip,jj,:),Hx(ip,jp,:));
130			end
131			Hx2sides=max( Hx2, Hx2([2:end 1],:,:) );
132			ll=find(Hx2mid>Hx2sides & Hx2mid>H2);
133			H2(ll)=Hx2mid(ll);
134			disp(' Filling cell due to staight meridional barriers')
135			if hxhyempty
136			Hy2mid=min( max(He(ii,jj,:),He(ii,jp,:)) , max(He(ip,jj,:),He(ip,jp,:)) );
137			else
138			Hy2mid=min(Hy(ii,jp,:),Hy(ip,jp,:));
139			end
140			Hy2sides=max( Hy2, Hy2(:,[2:end 1],:) );
141			ll=find(Hy2mid>Hy2sides & Hy2mid>H2);
142			H2(ll)=Hy2mid(ll);
143
144			% Make Hx2,Hy2 consistent with coarse H2
145			disp(' Limiting Hx2 by neighbouring H2')
146			Hx2=max( max(Hx2, H2), H2([end 1:end-1],:,:));
147			disp(' Limiting Hy2 by neighbouring H2')
148			Hy2=max( max(Hy2, H2), H2(:,[end 1:end-1],:));
149
150			% Check that Hx2 and Hy2 don't violate basic rules
151			%check_HxHy(H2,Hx2,Hy2)
152
153			% Re-cycle data for another iteration
154			He=H2;
155			Hx=Hx2;
156			Hy=Hy2;
157
158			if debuglevels
159			eval(sprintf('save LEVEL%2.2i He Hx Hy',nit));
160			end
161
162			end
163
164			% Fill-up cell-centers where all sides are higher
165			disp('Filling holes in H2 (bowls)')
166			H2sides=min( Hx2, Hx2([2:end 1],:,:) );
167			H2sides=min( H2sides, Hy2 );
168			H2sides=min( H2sides, Hy2(:,[2:end end],:) );
169			%H2=max( H2, H2sides );