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