| 1 |
gforget |
1.3 |
function [FLD]=diffsmooth2D_extrap_fwd(fld,mskOut,eps); |
| 2 |
|
|
%object: extrapolate an incomplete field to create a full field, by |
| 3 |
|
|
% time-stepping a diffusion equation to near-equilibrium. |
| 4 |
|
|
%inputs: fld incomplete field of interest (masked with NaN) |
| 5 |
|
|
% mskOut land mask (1s and NaNs) for the full field (output) |
| 6 |
|
|
% eps convergence criterium |
| 7 |
|
|
%output: FLD full field |
| 8 |
gforget |
1.1 |
|
| 9 |
|
|
global mygrid; |
| 10 |
|
|
|
| 11 |
|
|
dxC=mygrid.DXC; dyC=mygrid.DYC; |
| 12 |
|
|
dxG=mygrid.DXG; dyG=mygrid.DYG; |
| 13 |
|
|
rA=mygrid.RAC; |
| 14 |
|
|
|
| 15 |
|
|
dxCsm=dxC; dyCsm=dyC; |
| 16 |
|
|
mskFreeze=fld; mskFreeze(find(~isnan(mskFreeze)))=0; mskFreeze(find(isnan(mskFreeze)))=1; |
| 17 |
|
|
|
| 18 |
|
|
%put first guess: |
| 19 |
|
|
x=convert2array(mygrid.XC); |
| 20 |
|
|
y=convert2array(mygrid.YC); |
| 21 |
|
|
z=convert2array(fld); |
| 22 |
|
|
m=convert2array(mskOut); |
| 23 |
|
|
tmp1=find(~isnan(z)); |
| 24 |
|
|
tmp2=find(~isnan(m)); |
| 25 |
|
|
zz=z; |
| 26 |
|
|
zz(tmp2) = griddata(x(tmp1),y(tmp1),z(tmp1),x(tmp2),y(tmp2),'nearest'); |
| 27 |
gforget |
1.2 |
fld=convert2array(zz); |
| 28 |
gforget |
1.1 |
|
| 29 |
|
|
%put 0 first guess if needed and switch land mask: |
| 30 |
|
|
fld(find(isnan(fld)))=0; fld=fld.*mskOut; |
| 31 |
|
|
|
| 32 |
|
|
%scale the diffusive operator: |
| 33 |
|
|
tmp0=dxCsm./dxC; tmp0(isnan(mskOut))=NaN; tmp00=nanmax(tmp0); |
| 34 |
|
|
tmp0=dyCsm./dyC; tmp0(isnan(mskOut))=NaN; tmp00=max([tmp00 nanmax(tmp0)]); |
| 35 |
|
|
smooth2D_nbt=tmp00; |
| 36 |
|
|
smooth2D_nbt=ceil(1.1*2*smooth2D_nbt^2); |
| 37 |
|
|
|
| 38 |
|
|
smooth2D_dt=1; |
| 39 |
|
|
smooth2D_T=smooth2D_nbt*smooth2D_dt; |
| 40 |
|
|
smooth2D_Kux=dxCsm.*dxCsm/smooth2D_T/2; |
| 41 |
|
|
smooth2D_Kvy=dyCsm.*dyCsm/smooth2D_T/2; |
| 42 |
|
|
|
| 43 |
|
|
%time-stepping loop: |
| 44 |
|
|
FLD=fld; |
| 45 |
|
|
|
| 46 |
|
|
test1=0; it=0; |
| 47 |
|
|
while ~test1; |
| 48 |
|
|
|
| 49 |
|
|
it=it+1; |
| 50 |
|
|
|
| 51 |
|
|
[dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0); |
| 52 |
|
|
tmpU=dTdxAtU.*smooth2D_Kux; |
| 53 |
|
|
tmpV=dTdyAtV.*smooth2D_Kvy; |
| 54 |
|
|
[fldDIV]=calc_UV_div(tmpU,tmpV); |
| 55 |
|
|
dFLDdt=smooth2D_dt*fldDIV./rA; |
| 56 |
|
|
dFLDdt=dFLDdt.*mskFreeze; |
| 57 |
|
|
FLD=FLD-dFLDdt; |
| 58 |
|
|
|
| 59 |
|
|
tmp1=max(abs(dFLDdt)); |
| 60 |
|
|
if mod(it,10)==0; fprintf([num2str(it) ' del=' num2str(tmp1) ' eps=' num2str(eps) ' \n']); end; |
| 61 |
|
|
test1=tmp1<eps; |
| 62 |
|
|
|
| 63 |
|
|
FLDstore{it}=dFLDdt; |
| 64 |
|
|
end; |
| 65 |
|
|
|
| 66 |
|
|
if 1; |
| 67 |
|
|
jj=1; FLDstore2=zeros([size(FLD{jj}) length(FLDstore)]); |
| 68 |
|
|
for ii=1:length(FLDstore); FLDstore2(:,:,ii)=FLDstore{ii}{jj}; end; |
| 69 |
|
|
tmp1=abs(FLDstore2(:,:,end)); [ii,jj]=find(tmp1==max(tmp1(:))); |
| 70 |
|
|
figure; plot(cumsum(squeeze(FLDstore2(ii,jj,:)))); |
| 71 |
|
|
end; |
| 72 |
|
|
|
| 73 |
|
|
|
| 74 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph