| 1 |
gforget |
1.1 |
function [fldOut]=gcmfaces_remap_3d(lon,lat,depIn,fldIn); |
| 2 |
|
|
%object: use bin average to remap a lat-lon-depth field to mygrid.mskC |
| 3 |
|
|
%input: lon,lat,dep,fld are the gridded product arrays (2D,2D,1D,3D) |
| 4 |
|
|
%assumption:fld should show NaN for missing values |
| 5 |
|
|
|
| 6 |
|
|
gcmfaces_global; |
| 7 |
|
|
|
| 8 |
|
|
%initiate triangulation: |
| 9 |
|
|
gcmfaces_bindata; |
| 10 |
|
|
|
| 11 |
|
|
%vertical levels: |
| 12 |
|
|
nrIn=length(depIn); |
| 13 |
|
|
depOut=-mygrid.RC; |
| 14 |
|
|
nrOut=length(depOut); |
| 15 |
|
|
|
| 16 |
|
|
%corresponding model mask: chosen to avoid vertical extrapolation later |
| 17 |
|
|
% 1) from model just above |
| 18 |
|
|
mskC=repmat(0*mygrid.XC,[1 1 nrIn]); |
| 19 |
|
|
KK2=zeros(1,nrIn); |
| 20 |
|
|
for kk=1:nrIn; |
| 21 |
|
|
kk2=max(find(depOut<=depIn(kk)));%model level just above |
| 22 |
|
|
if isempty(kk2); kk2=1; end;%pathological case |
| 23 |
|
|
KK2(kk)=kk2; |
| 24 |
|
|
mskC(:,:,kk)=mygrid.mskC(:,:,kk2); |
| 25 |
|
|
end; |
| 26 |
|
|
% 2) shift to above atlas level mask |
| 27 |
|
|
for kk=nrIn:-1:2; |
| 28 |
|
|
mskC(:,:,kk)=mskC(:,:,kk-1);%atlas level just above |
| 29 |
|
|
end; |
| 30 |
|
|
|
| 31 |
|
|
%test vertical interpolation to mygrid.RC of the mask: |
| 32 |
|
|
% => the printed result should be 0 |
| 33 |
|
|
if 0; |
| 34 |
|
|
mskCout=atlas_interp_vert(mskC,depIn,depOut); |
| 35 |
|
|
aa=convert2gcmfaces(isnan(mskCout)&~isnan(mygrid.mskC)); |
| 36 |
|
|
nansum(aa(:)) |
| 37 |
|
|
end; |
| 38 |
|
|
|
| 39 |
|
|
%horizontal mapping to mygrid.XC, mygrid.YC, mskC: |
| 40 |
|
|
fld1=mskC; |
| 41 |
|
|
tocs=[]; |
| 42 |
|
|
for kk=1:nrIn; |
| 43 |
|
|
tic; |
| 44 |
|
|
fld1(:,:,kk)=gcmfaces_remap_2d(lon,lat,fldIn(:,:,kk),3,mskC(:,:,kk)); |
| 45 |
|
|
tocs=[tocs;toc]; |
| 46 |
|
|
[kk nrIn mean(tocs)] |
| 47 |
|
|
end; |
| 48 |
|
|
|
| 49 |
|
|
%vertical interpolation to mygrid.RC |
| 50 |
|
|
fld2=atlas_interp_vert(fld1,depIn,depOut); |
| 51 |
|
|
|
| 52 |
|
|
|
| 53 |
|
|
%final mask handling: |
| 54 |
|
|
%1) remask with mygrid.mskC (mskC was wider by design) |
| 55 |
|
|
msk1=mygrid.mskC; tmp1=msk1.*fld2; |
| 56 |
|
|
%2) extrapolate vertically (some deep canions may still miss data) |
| 57 |
|
|
fldOut=atlas_extrap_vert(tmp1,msk1,depOut); |
| 58 |
|
|
|
| 59 |
|
|
|
| 60 |
|
|
function [fldOut]=atlas_interp_vert(fldIn,depthIn,depthOut); |
| 61 |
|
|
%vertical interpolation |
| 62 |
|
|
|
| 63 |
|
|
gcmfaces_global; |
| 64 |
|
|
|
| 65 |
|
|
nrIn=length(depthIn); |
| 66 |
|
|
nrOut=length(depthOut); |
| 67 |
|
|
|
| 68 |
|
|
%map depthIn to depthOut: |
| 69 |
|
|
coeff=interp1(depthIn,[1:nrIn],depthOut); |
| 70 |
|
|
tmp1=find(isnan(coeff)); coeff(tmp1)=nrIn-0.01; %quick fix |
| 71 |
|
|
|
| 72 |
|
|
%vertical interpolation to mygrid.RC: |
| 73 |
|
|
fldOut=repmat(0*mygrid.XC,[1 1 nrOut]); |
| 74 |
|
|
for kk=1:nrOut; |
| 75 |
|
|
tmp1=coeff(kk); tmp2=floor(tmp1); tmp1=tmp1-tmp2; |
| 76 |
|
|
if tmp2==nrIn; |
| 77 |
|
|
tmp3=fldIn(:,:,tmp2); |
| 78 |
|
|
else; |
| 79 |
|
|
tmp3=(1-tmp1)*fldIn(:,:,tmp2)+tmp1*fldIn(:,:,tmp2+1); |
| 80 |
|
|
end; |
| 81 |
|
|
fldOut(:,:,kk)=tmp3; |
| 82 |
|
|
end; |
| 83 |
|
|
|
| 84 |
|
|
function [fldOut]=atlas_extrap_vert(fldIn,mskOut,depthIn); |
| 85 |
|
|
%vertical extrapolation when needed |
| 86 |
|
|
|
| 87 |
|
|
fldOut=convert2array(fldIn); |
| 88 |
|
|
mskOut=convert2array(mskOut); |
| 89 |
|
|
nrIn=length(depthIn); |
| 90 |
|
|
|
| 91 |
|
|
for kk=2:nrIn; |
| 92 |
|
|
tmp_kk=fldOut(:,:,kk); |
| 93 |
|
|
msk_kk=mskOut(:,:,kk); |
| 94 |
|
|
[II,JJ]=find(isnan(tmp_kk)&~isnan(msk_kk)); |
| 95 |
|
|
if ~isempty(II)&kk==2; |
| 96 |
|
|
tmp_kkM1=fldOut(:,:,kk-1); |
| 97 |
|
|
tmp1=find(isnan(tmp_kk)&~isnan(msk_kk)&~isnan(tmp_kkM1)); |
| 98 |
|
|
tmp_kk(tmp1)=tmp_kkM1(tmp1); |
| 99 |
|
|
fldOut(:,:,kk)=tmp_kk; |
| 100 |
|
|
elseif ~isempty(II); |
| 101 |
|
|
tmp_kkM1=fldOut(:,:,kk-1); |
| 102 |
|
|
tmp_kkM2=fldOut(:,:,kk-2); |
| 103 |
|
|
coeff=(depthIn(kk)-depthIn(kk-1))/(depthIn(kk-1)-depthIn(kk-2)); |
| 104 |
|
|
tmp_kk_x=tmp_kkM1+coeff*(tmp_kkM1-tmp_kkM2); |
| 105 |
|
|
tmp1=find(isnan(tmp_kk)&~isnan(msk_kk)&~isnan(tmp_kk_x)); |
| 106 |
|
|
tmp_kk(tmp1)=tmp_kk_x(tmp1); |
| 107 |
|
|
fldOut(:,:,kk)=tmp_kk; |
| 108 |
|
|
end; |
| 109 |
|
|
end; |
| 110 |
|
|
|
| 111 |
|
|
fldOut=convert2array(fldOut); |
| 112 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph