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