mlosch/interp_llc/xyrecur_fv_loopbased.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);

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;
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	debuglevels=0;
11
12	%H2=He;
13	%Hx2=Hx;
14	%Hy2=Hy;
15
16	% Check that Hx and Hy don't violate basic rules
17	check_HxHy(He,Hx,Hy)
18
19	% Recursively iterate
20	for nit=1:niter,
21
22	clear H2 Hx2 Hy2
23	[nx,ny,nt]=size(He);
24	disp(sprintf(' Recursion level %i (%i,%i) -> (%i,%i) ... ',nit,nx,ny,nx/2,ny/2));
25
26	%ii=1:2:nx-1;
27	%jj=1:2:ny-1;
28	%ip=2:2:nx;
29	%im=[nx 2:2:nx-2];
30	%jp=2:2:ny;
31	%jm=[2 2:2:ny-2];
32	%I=1:nx/2; IM=I([end 1:end-1]);
33	%J=1:ny/2; JM=J([end 1:end-1]);
34
35	for J=1:ny/2; jj=1+2*(J-1);
36	jp=mymod(jj+1,ny);jm=mymod(jj-1,ny);JM=mymod(J-1,ny/2);
37	for I=1:nx/2; ii=1+2*(I-1);
38	ip=mymod(ii+1,nx);im=mymod(ii-1,nx);IM=mymod(I-1,nx/2);
39
40	% Four-point average in X-Y
41	%new? H2=(He(ii,jj,:)+He(ii+1,jj,:)+He(ii,jj+1,:)+He(ii+1,jj+1,:))/4;
42	H2(I,J,:)=(He(ii,jj,:)+He(ip,jj,:)+He(ii,jp,:)+He(ip,jp,:))/4;
43
44	% Two-point average in X
45	%new? Hx2=(Hx(ii,jj,:)+Hx(ii,jj+1,:))/2;
46	Hx2(I,J,:)=(Hx(ii,jj,:)+Hx(ii,jp,:))/2;
47
48	% Two-point average in Y
49	%new? Hy2=(Hy(ii,jj,:)+Hy(ii+1,jj,:))/2;
50	Hy2(I,J,:)=(Hy(ii,jj,:)+Hy(ip,jj,:))/2;
51
52	end;end;
53	for J=1:ny/2; jj=1+2*(J-1);
54	jp=mymod(jj+1,ny),jm=mymod(jj-1,ny),JM=mymod(J-1,ny/2),
55	for I=1:nx/2; ii=1+2*(I-1);
56	ip=mymod(ii+1,nx),im=mymod(ii-1,nx),IM=mymod(I-1,nx/2),
57
58	% Blocking from either side?
59	Hx2e=(Hx(ip,jj,:)+Hx(ip,jp,:))/2;
60	Hx2w=(Hx(im,jj,:)+Hx(im,jp,:))/2;
61	Hx2(I,J,:)=max( max( Hx2(I,J,:), Hx2e ) , Hx2w);
62
63	% Blocking from either side?
64	Hy2n=(Hy(ii,jp,:)+Hy(ip,jp,:))/2;
65	Hy2s=(Hy(ii,jm,:)+Hy(ip,jm,:))/2;
66	Hy2(I,J,:)=max( max( Hy2(I,J,:), Hy2n ) , Hy2s);
67
68	% Make Hx2,Hy2 consistent with coarse H2
69	%Hx2=max( max(Hx2, H2), H2([end 1:end-1],:,:));
70	%Hy2=max( max(Hy2, H2), H2(:,[end 1:end-1],:));
71
72	% Fill-up cell-centers where all sides are higher
73	%H2sides=min( Hx2, Hx2([2:end 1],:,:) );
74	%H2sides=min( H2sides, Hy2 );
75	%H2sides=min( H2sides, Hy2(:,[2:end end],:) );
76	%H2=max( H2, H2sides );
77
78	% Fill up if barriers intersects coarse cell
79	%%Hx2mid=(Hx(ip,jj,:)+Hx(ip,jp,:))/2;
80	%%Hy2mid=(Hy(ii,jp,:)+Hy(ip,jp,:))/2;
81	%Hx2mid=min(Hx(ip,jj,:),Hx(ip,jp,:));
82	%Hy2mid=min(Hy(ii,jp,:),Hy(ip,jp,:));
83	%H2mid=max(max(H2, Hy2mid),Hx2mid);
84	%ll=find(H2mid>H2sides);
85	%H2(ll)=H2mid(ll);
86
87	% Make Hx2,Hy2 consistent with coarse H2
88	%Hx2=max( max(Hx2, H2), H2([end 1:end-1],:,:));
89	%Hy2=max( max(Hy2, H2), H2(:,[end 1:end-1],:));
90	Hx2(I,J,:)=max( max(Hx2(I,J,:), H2(I,J,:)), H2(IM,J,:));
91	Hy2(I,J,:)=max( max(Hy2(I,J,:), H2(I,J,:)), H2(I,JM,:));
92
93	end;end;
94
95	% Check that Hx2 and Hy2 don't violate basic rules
96	check_HxHy(H2,Hx2,Hy2)
97
98	% Re-cycle data for another iteration
99	He=H2;
100	Hx=Hx2;
101	Hy=Hy2;
102
103	if debuglevels
104	eval(sprintf('save LEVEL%2.2i He Hx Hy',nit));
105	end
106
107	end
108
109	% Fill-up cell-centers where all sides are higher
110	H2sides=min( Hx2, Hx2([2:end 1],:,:) );
111	H2sides=min( H2sides, Hy2 );
112	H2sides=min( H2sides, Hy2(:,[2:end end],:) );
113	H2=max( H2, H2sides );
114
115
116
117	function []=QQQQQcheck_HxHy(H2,Hx2,Hy2)
118	% Check that Hx2 and Hy2 don't violate basic rules
119	% This should not ever happen
120	return; % spead up process
121	[Hxcc,Hycc]=gen_hxhy(H2,x2,y2);
122	ic=find(Hx2<Hxcc-0.001);
123	if prod(size(ic)) ~= 0
124	[Hx2(ic)/1e3 Hxcc(ic)/1e3 Hx2(ic)-Hxcc(ic)]
125	error('The X cross-sectional area turned out inconsistently')
126	end
127	jc=find(Hy2<Hycc-0.001);
128	if prod(size(jc)) ~= 0
129	[Hy2(jc)/1e3 Hycc(jc)/1e3 Hy2(jc)-Hycc(jc)]
130	error('The Y cross-sectional area turned out inconsistently')
131	end
132
133
134	function [n] = mymod(n,N)
135	n=mod(n+N-1,N)+1;