matlab_class/gcmfaces_smooth/diffsmooth2D_extrap_inv.m

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

%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;

%check for domain edge points where no exchange is possible:
tmp1=mskOut; tmp1(:)=1; tmp2=exch_T_N(tmp1);
for iF=1:mskOut.nFaces;
tmp3=mskOut{iF}; tmp4=tmp2{iF}; 
tmp4=tmp4(2:end-1,1:end-2)+tmp4(2:end-1,3:end)+tmp4(1:end-2,2:end-1)+tmp4(3:end,2:end-1);
if ~isempty(find(isnan(tmp4)&~isnan(tmp3))); fprintf('warning: mask was modified\n'); end;
tmp3(isnan(tmp4))=NaN; mskOut{iF}=tmp3; 
end;

%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;

%form matrix problem:
tmp1=convert2array(mskOut);
kk=find(~isnan(tmp1));
KK=tmp1; KK(kk)=kk; KK=convert2array(KK,fld);
nn=length(kk);
NN=tmp1; NN(kk)=[1:nn]; %NN=convert2array(NN,fld); 

if nargin==3; doFormMatrix=varargin{1}; else; doFormMatrix=1; end;

global dFLDdt_op; 
if doFormMatrix==1;

dFLDdt_op=sparse([],[],[],nn,nn,nn*5);

for iF=1:fld.nFaces; for ii=1:3; for jj=1:3;
    FLDones=fld; FLDones(find(~isnan(fld)))=0; 
    FLDones{iF}(ii:3:end,jj:3:end)=1;
    FLDones(find(isnan(fld)))=NaN;

    FLDkkFROMtmp=fld; FLDkkFROMtmp(find(~isnan(fld)))=0;
    FLDkkFROMtmp{iF}(ii:3:end,jj:3:end)=KK{iF}(ii:3:end,jj:3:end);
    FLDkkFROMtmp(find(isnan(fld)))=0;

    FLDkkFROM=exch_T_N(FLDkkFROMtmp); 
    FLDkkFROM(find(isnan(FLDkkFROM)))=0;
    for iF2=1:fld.nFaces; 
       tmp1=FLDkkFROM{iF2}; tmp2=zeros(size(tmp1)-2);
       for ii2=1:3; for jj2=1:3; tmp2=tmp2+tmp1(ii2:end-3+ii2,jj2:end-3+jj2); end; end;
       FLDkkFROM{iF2}=tmp2; 
    end;
    %clear FLDkkFROMtmp;

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

    dFLDdt=convert2array(dFLDdt); FLDkkFROM=convert2array(FLDkkFROM); FLDkkTO=convert2array(KK);
    tmp1=find(dFLDdt~=0&~isnan(dFLDdt)); 
    dFLDdt=dFLDdt(tmp1); FLDkkFROM=FLDkkFROM(tmp1); FLDkkTO=FLDkkTO(tmp1);
    dFLDdt_op=dFLDdt_op+sparse(NN(FLDkkTO),NN(FLDkkFROM),dFLDdt,nn,nn);

end; end; end;

end;%if doFormMatrix==1;

%figure; spy(dFLDdt_op);

FLD_vec=convert2array(fld);
mskFreeze_vec=convert2array(mskFreeze);
FLD_vec(find(mskFreeze_vec==1))=0;
FLD_vec=FLD_vec(kk);

INV_op=1-mskFreeze_vec(kk);
INV_op=sparse([1:nn],[1:nn],INV_op,nn,nn);
INV_op=INV_op+dFLDdt_op;
INV_vec=INV_op\FLD_vec;
INV_fld=convert2array(mskOut);
INV_fld(find(~isnan(INV_fld)))=INV_vec;

FLD=convert2array(INV_fld,fld);

return;

%time step using matrix:
FLD_vec=convert2array(fld);
FLD_vec=FLD_vec(find(~isnan(FLD_vec)));
dFLDdt_vec=dFLDdt_op*FLD_vec;
dFLDdt_fld=convert2array(FLD); 
dFLDdt_fld(find(~isnan(dFLDdt_fld)))=dFLDdt_vec;
dFLDdt_fld(find(isnan(dFLDdt_fld)))=0;

%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.2	function [FLD]=diffsmooth2D(fld,mskOut,varargin);
2	gforget	1.1
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			%check for domain edge points where no exchange is possible:
21			tmp1=mskOut; tmp1(:)=1; tmp2=exch_T_N(tmp1);
22			for iF=1:mskOut.nFaces;
23			tmp3=mskOut{iF}; tmp4=tmp2{iF};
24			tmp4=tmp4(2:end-1,1:end-2)+tmp4(2:end-1,3:end)+tmp4(1:end-2,2:end-1)+tmp4(3:end,2:end-1);
25			if ~isempty(find(isnan(tmp4)&~isnan(tmp3))); fprintf('warning: mask was modified\n'); end;
26			tmp3(isnan(tmp4))=NaN; mskOut{iF}=tmp3;
27			end;
28
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.12smooth2D_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			%form matrix problem:
44			tmp1=convert2array(mskOut);
45			kk=find(~isnan(tmp1));
46			KK=tmp1; KK(kk)=kk; KK=convert2array(KK,fld);
47			nn=length(kk);
48			NN=tmp1; NN(kk)=[1:nn]; %NN=convert2array(NN,fld);
49
50	gforget	1.2	if nargin==3; doFormMatrix=varargin{1}; else; doFormMatrix=1; end;
51
52			global dFLDdt_op;
53			if doFormMatrix==1;
54
55	gforget	1.1	dFLDdt_op=sparse([],[],[],nn,nn,nn*5);
56
57			for iF=1:fld.nFaces; for ii=1:3; for jj=1:3;
58			FLDones=fld; FLDones(find(~isnan(fld)))=0;
59			FLDones{iF}(ii:3:end,jj:3:end)=1;
60			FLDones(find(isnan(fld)))=NaN;
61
62			FLDkkFROMtmp=fld; FLDkkFROMtmp(find(~isnan(fld)))=0;
63			FLDkkFROMtmp{iF}(ii:3:end,jj:3:end)=KK{iF}(ii:3:end,jj:3:end);
64			FLDkkFROMtmp(find(isnan(fld)))=0;
65
66			FLDkkFROM=exch_T_N(FLDkkFROMtmp);
67	gforget	1.3	FLDkkFROM(find(isnan(FLDkkFROM)))=0;
68	gforget	1.1	for iF2=1:fld.nFaces;
69			tmp1=FLDkkFROM{iF2}; tmp2=zeros(size(tmp1)-2);
70			for ii2=1:3; for jj2=1:3; tmp2=tmp2+tmp1(ii2:end-3+ii2,jj2:end-3+jj2); end; end;
71			FLDkkFROM{iF2}=tmp2;
72			end;
73			%clear FLDkkFROMtmp;
74
75			[dTdxAtU,dTdyAtV]=calc_T_grad(FLDones,0);
76			tmpU=dTdxAtU.*smooth2D_Kux;
77			tmpV=dTdyAtV.*smooth2D_Kvy;
78			[fldDIV]=calc_UV_div(tmpU,tmpV);
79			dFLDdt=smooth2D_dt*fldDIV./rA;
80			dFLDdt=dFLDdt.*mskFreeze;
81
82			dFLDdt=convert2array(dFLDdt); FLDkkFROM=convert2array(FLDkkFROM); FLDkkTO=convert2array(KK);
83			tmp1=find(dFLDdt~=0&~isnan(dFLDdt));
84			dFLDdt=dFLDdt(tmp1); FLDkkFROM=FLDkkFROM(tmp1); FLDkkTO=FLDkkTO(tmp1);
85			dFLDdt_op=dFLDdt_op+sparse(NN(FLDkkTO),NN(FLDkkFROM),dFLDdt,nn,nn);
86
87			end; end; end;
88
89	gforget	1.2	end;%if doFormMatrix==1;
90
91	gforget	1.1	%figure; spy(dFLDdt_op);
92
93			FLD_vec=convert2array(fld);
94			mskFreeze_vec=convert2array(mskFreeze);
95			FLD_vec(find(mskFreeze_vec==1))=0;
96			FLD_vec=FLD_vec(kk);
97
98			INV_op=1-mskFreeze_vec(kk);
99			INV_op=sparse([1:nn],[1:nn],INV_op,nn,nn);
100			INV_op=INV_op+dFLDdt_op;
101			INV_vec=INV_op\FLD_vec;
102			INV_fld=convert2array(mskOut);
103			INV_fld(find(~isnan(INV_fld)))=INV_vec;
104
105			FLD=convert2array(INV_fld,fld);
106
107			return;
108
109			%time step using matrix:
110			FLD_vec=convert2array(fld);
111			FLD_vec=FLD_vec(find(~isnan(FLD_vec)));
112			dFLDdt_vec=dFLDdt_op*FLD_vec;
113			dFLDdt_fld=convert2array(FLD);
114			dFLDdt_fld(find(~isnan(dFLDdt_fld)))=dFLDdt_vec;
115			dFLDdt_fld(find(isnan(dFLDdt_fld)))=0;
116
117			%time-stepping loop:
118			FLD=fld;
119
120			test1=0; it=0;
121			while ~test1;
122
123			it=it+1;
124
125			[dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0);
126			tmpU=dTdxAtU.*smooth2D_Kux;
127			tmpV=dTdyAtV.*smooth2D_Kvy;
128			[fldDIV]=calc_UV_div(tmpU,tmpV);
129			dFLDdt=smooth2D_dt*fldDIV./rA;
130			dFLDdt=dFLDdt.*mskFreeze;
131			FLD=FLD-dFLDdt;
132
133			tmp1=max(abs(dFLDdt));
134			if mod(it,10)==0; fprintf([num2str(it) ' del=' num2str(tmp1) ' eps=' num2str(eps) ' \n']); end;
135			test1=tmp1<eps;
136
137			FLDstore{it}=dFLDdt;
138			end;
139
140			if 1;
141			jj=1; FLDstore2=zeros([size(FLD{jj}) length(FLDstore)]);
142			for ii=1:length(FLDstore); FLDstore2(:,:,ii)=FLDstore{ii}{jj}; end;
143			tmp1=abs(FLDstore2(:,:,end)); [ii,jj]=find(tmp1==max(tmp1(:)));
144			figure; plot(cumsum(squeeze(FLDstore2(ii,jj,:))));
145			end;
146
147
148