| 1 |
mlosch |
1.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