| 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),2*li,2*lj) )=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),2*li,2*lj-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),2*li-1,2*lj-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),2*li-1,2*lj) )=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 ); |
Parent Directory
| 
Revision Log
| 
Revision Graph