matlab_class/gcmfaces_smooth/diffsmooth2D_extrap_fwd.m

function [FLD]=diffsmooth2D(fld,mskOut,eps);

%object: implementation (gforget@mit.edu) of a diffusive smoother (Weaver and Courtier 2001)

%input: fld     field to be smoothed (masked with NaN)
%       dxCsm,dyCsm     scale in first/second direction
%output:FLD     smoothed field

%asumption: dxCsm/dyCsm are assumed to be given at the positions of U/V points

global mygrid;

dxC=mygrid.DXC; dyC=mygrid.DYC; 
dxG=mygrid.DXG; dyG=mygrid.DYG; 
rA=mygrid.RAC;

dxCsm=dxC; dyCsm=dyC;
mskFreeze=fld; mskFreeze(find(~isnan(mskFreeze)))=0; mskFreeze(find(isnan(mskFreeze)))=1;

%put first guess:
x=convert2array(mygrid.XC);
y=convert2array(mygrid.YC);
z=convert2array(fld);
m=convert2array(mskOut);
tmp1=find(~isnan(z));
tmp2=find(~isnan(m));
zz=z;
zz(tmp2) = griddata(x(tmp1),y(tmp1),z(tmp1),x(tmp2),y(tmp2),'nearest');
fld=convert2array(zz);

%put 0 first guess if needed and switch land mask:
fld(find(isnan(fld)))=0; fld=fld.*mskOut;

%scale the diffusive operator:
tmp0=dxCsm./dxC; tmp0(isnan(mskOut))=NaN; tmp00=nanmax(tmp0);
tmp0=dyCsm./dyC; tmp0(isnan(mskOut))=NaN; tmp00=max([tmp00 nanmax(tmp0)]);
smooth2D_nbt=tmp00;
smooth2D_nbt=ceil(1.1*2*smooth2D_nbt^2);

smooth2D_dt=1;
smooth2D_T=smooth2D_nbt*smooth2D_dt;
smooth2D_Kux=dxCsm.*dxCsm/smooth2D_T/2;
smooth2D_Kvy=dyCsm.*dyCsm/smooth2D_T/2;

%time-stepping loop:
FLD=fld;

test1=0; it=0;
while ~test1;

    it=it+1;

    [dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0);
    tmpU=dTdxAtU.*smooth2D_Kux;
    tmpV=dTdyAtV.*smooth2D_Kvy;
    [fldDIV]=calc_UV_div(tmpU,tmpV);
    dFLDdt=smooth2D_dt*fldDIV./rA;
    dFLDdt=dFLDdt.*mskFreeze;
    FLD=FLD-dFLDdt;

    tmp1=max(abs(dFLDdt));
    if mod(it,10)==0; fprintf([num2str(it) ' del=' num2str(tmp1) ' eps=' num2str(eps) ' \n']); end;
    test1=tmp1<eps;

    FLDstore{it}=dFLDdt;
end;

if 1;
   jj=1; FLDstore2=zeros([size(FLD{jj}) length(FLDstore)]);
   for ii=1:length(FLDstore); FLDstore2(:,:,ii)=FLDstore{ii}{jj}; end;
   tmp1=abs(FLDstore2(:,:,end)); [ii,jj]=find(tmp1==max(tmp1(:)));
   figure; plot(cumsum(squeeze(FLDstore2(ii,jj,:))));
end;


1	gforget	1.1	function [FLD]=diffsmooth2D(fld,mskOut,eps);
2
3			%object: implementation (gforget@mit.edu) of a diffusive smoother (Weaver and Courtier 2001)
4
5			%input: fld field to be smoothed (masked with NaN)
6			% dxCsm,dyCsm scale in first/second direction
7			%output:FLD smoothed field
8
9			%asumption: dxCsm/dyCsm are assumed to be given at the positions of U/V points
10
11			global mygrid;
12
13			dxC=mygrid.DXC; dyC=mygrid.DYC;
14			dxG=mygrid.DXG; dyG=mygrid.DYG;
15			rA=mygrid.RAC;
16
17			dxCsm=dxC; dyCsm=dyC;
18			mskFreeze=fld; mskFreeze(find(~isnan(mskFreeze)))=0; mskFreeze(find(isnan(mskFreeze)))=1;
19
20			%put first guess:
21			x=convert2array(mygrid.XC);
22			y=convert2array(mygrid.YC);
23			z=convert2array(fld);
24			m=convert2array(mskOut);
25			tmp1=find(~isnan(z));
26			tmp2=find(~isnan(m));
27			zz=z;
28			zz(tmp2) = griddata(x(tmp1),y(tmp1),z(tmp1),x(tmp2),y(tmp2),'nearest');
29	gforget	1.2	fld=convert2array(zz);
30	gforget	1.1
31			%put 0 first guess if needed and switch land mask:
32			fld(find(isnan(fld)))=0; fld=fld.*mskOut;
33
34			%scale the diffusive operator:
35			tmp0=dxCsm./dxC; tmp0(isnan(mskOut))=NaN; tmp00=nanmax(tmp0);
36			tmp0=dyCsm./dyC; tmp0(isnan(mskOut))=NaN; tmp00=max([tmp00 nanmax(tmp0)]);
37			smooth2D_nbt=tmp00;
38			smooth2D_nbt=ceil(1.12smooth2D_nbt^2);
39
40			smooth2D_dt=1;
41			smooth2D_T=smooth2D_nbt*smooth2D_dt;
42			smooth2D_Kux=dxCsm.*dxCsm/smooth2D_T/2;
43			smooth2D_Kvy=dyCsm.*dyCsm/smooth2D_T/2;
44
45			%time-stepping loop:
46			FLD=fld;
47
48			test1=0; it=0;
49			while ~test1;
50
51			it=it+1;
52
53			[dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0);
54			tmpU=dTdxAtU.*smooth2D_Kux;
55			tmpV=dTdyAtV.*smooth2D_Kvy;
56			[fldDIV]=calc_UV_div(tmpU,tmpV);
57			dFLDdt=smooth2D_dt*fldDIV./rA;
58			dFLDdt=dFLDdt.*mskFreeze;
59			FLD=FLD-dFLDdt;
60
61			tmp1=max(abs(dFLDdt));
62			if mod(it,10)==0; fprintf([num2str(it) ' del=' num2str(tmp1) ' eps=' num2str(eps) ' \n']); end;
63			test1=tmp1<eps;
64
65			FLDstore{it}=dFLDdt;
66			end;
67
68			if 1;
69			jj=1; FLDstore2=zeros([size(FLD{jj}) length(FLDstore)]);
70			for ii=1:length(FLDstore); FLDstore2(:,:,ii)=FLDstore{ii}{jj}; end;
71			tmp1=abs(FLDstore2(:,:,end)); [ii,jj]=find(tmp1==max(tmp1(:)));
72			figure; plot(cumsum(squeeze(FLDstore2(ii,jj,:))));
73			end;
74
75
76