| 1 |
function [Tm,xc,yc,zc] = interp3ddata_int(xl,yl,zl,TL,xc,yc,zc,msk); |
| 2 |
% |
| 3 |
% [T,x,y,z] = interp3ddata_int(xl,yl,zl,TL,xc,yc,zc,msk); |
| 4 |
|
| 5 |
nxc=size(xc,1); |
| 6 |
nyc=size(yc,2); |
| 7 |
[nzc]=prod(size(zc)); |
| 8 |
[nxl,nyl,nzl]=size(TL); |
| 9 |
|
| 10 |
% Back out the grid spacing for model |
| 11 |
dz(1)=-zc(1)*2; |
| 12 |
zf(1)=0; |
| 13 |
for k=1:nzc-1; |
| 14 |
zf(k+1)=zf(k)-dz(k); |
| 15 |
dz(k+1)=2*(zf(k+1)-zc(k+1)); |
| 16 |
end |
| 17 |
zf(nzc+1)=zf(nzc)-dz(nzc); |
| 18 |
|
| 19 |
% Assume some grid spacing for Levitus |
| 20 |
zlf=[0 (zl+zl([2:end end]))/2]; |
| 21 |
dzl=zlf(1:nzl)-zlf(2:nzl+1); |
| 22 |
|
| 23 |
% Use NaN's to represent no data |
| 24 |
TL=sq(TL); |
| 25 |
|
| 26 |
% Create halo in X direction |
| 27 |
disp('Adding wrap around data') |
| 28 |
xl=[xl(end-2:end)-360 xl xl(1:3)+360]; |
| 29 |
for k=1:size(TL,3); |
| 30 |
TLL(:,:,k)=[TL(end-2:end,:,k)' TL(:,:,k)' TL(1:3,:,k)']'; |
| 31 |
end |
| 32 |
TL=TLL; clear TLL |
| 33 |
|
| 34 |
% Interpolate horiztonally |
| 35 |
for k=1:size(TL,3); |
| 36 |
disp( sprintf('Interpolating (x-y) level %i',k) ); |
| 37 |
Tl=interp2(yl,xl,TL(:,:,k),yc,xc,'cubic'); |
| 38 |
Te=interp2(yl,xl,TL(:,:,k),yc,xc,'linear'); |
| 39 |
jj=find( isnan(Tl) ) ; Tl(jj)=Te(jj); |
| 40 |
Te=interp2(yl,xl,TL(:,:,k),yc,xc,'nearest'); |
| 41 |
jj=find( isnan(Tl) ) ; Tl(jj)=Te(jj); |
| 42 |
%Tl(find(isnan(Tl)))=0; |
| 43 |
Tl=xyexpand(Tl.*sq(msk(:,:,1)),11); |
| 44 |
Tm(:,:,k)=Tl.*sq(msk(:,:,1)); |
| 45 |
end |
| 46 |
clear Tl Te jj TL |
| 47 |
|
| 48 |
|
| 49 |
disp('Vertical interpolation') |
| 50 |
Tm( find(isnan(Tm)) )=0; |
| 51 |
hm=reshape( -sum(reshape(msk,[nxc*nyc nzc])*dz',3) ,[nxc nyc]); |
| 52 |
hl=reshape( -sum(reshape(mask(Tm),[nxc*nyc nzl])*dzl',3) ,[nxc nyc]); |
| 53 |
for kk=1:nzc; |
| 54 |
z1m=max(zf(kk),hm); |
| 55 |
z2m=max(zf(kk+1),hm); |
| 56 |
dm=max(z1m-z2m,0); |
| 57 |
dw=0; |
| 58 |
ti=0; |
| 59 |
for k=1:nzl; |
| 60 |
z1l=max(zlf(k),hl); |
| 61 |
z2l=max(zlf(k+1),hl); |
| 62 |
z1=min(z1m,z1l); |
| 63 |
% z2=min(z2m,z2l); |
| 64 |
z2=max(z2m,z2l); |
| 65 |
dd=max(z1-z2,0); |
| 66 |
dw=dw+dd; |
| 67 |
ti=ti+dd.*Tm(:,:,k); |
| 68 |
end |
| 69 |
Tl(:,:,kk)=ti./sq(dw); |
| 70 |
end |
| 71 |
Tm=Tl.*msk; |
| 72 |
clear Tl |
| 73 |
|
| 74 |
Tmsk=mask(Tm); |
| 75 |
it=0; |
| 76 |
while( sum(Tmsk(:)-msk(:)) & it<15) |
| 77 |
it=it+1; |
| 78 |
disp( sprintf('Number of points needing vertical extrapolation = %i', ... |
| 79 |
sum(msk(:)-Tmsk(:)) ) ); |
| 80 |
Te=2*Tm(:,:,[1 1:end-1])-Tm(:,:,[1 1 1:end-2]); |
| 81 |
jj=find( Tmsk==0 & msk==1 ); |
| 82 |
Tm(jj)=Te(jj); |
| 83 |
Tmsk=mask(Tm); |
| 84 |
end |
| 85 |
|
| 86 |
% Remove NaN's |
| 87 |
Tm( find(isnan(Tm)) )=0; |
Parent Directory
| 
Revision Log
| 
Revision Graph