matlab_class/gcmfaces_smooth/diffsmooth2D_extrap_inv.m

function [FLD]=diffsmooth2D_extrap_inv(fld,mskOut,varargin);
%object:        extrapolate an incomplete field to create a full field, by 
%                       solving for a diffusion equation equilibrium state.
%inputs:        fld     incomplete field of interest (masked with NaN)
%               mskOut  land mask (1s and NaNs) for the full field (output)
%output:        FLD     full field
%optional:      doFormMatrix (1 by default)
%                       if set to 1 then compute the dFLDdt_op matrix
%                       if set to 0 then use precomputed one (global dFLDdt_op)
%               doMSKOUT (1 by default)
%                       if set to 1 then first detect closed regions lacking
%                         data points to extrapolate from (e.g. abyssal canions)
%                         and restrict mskOut to MSKOUT accordingly
%                       if set to 0 then omit this step
%principle : 
%in points where mskFreeze is 1 solve diffusion equation
%in points where mskFreeze is 0 solve the trivial FLD=fld

gcmfaces_global;

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

if doMSKOUT;
  %problematic regions will show MSKOUT==0;
  tmp1=1+0*fld; [MSKOUT]=diffsmooth2D_extrap_inv(tmp1,mskOut,doFormMatrix,0);
  %if no problematic region then no need to recompute dFLDdt_op 
  if sum(MSKOUT==0)==0; doFormMatrix=0; end;
  %finalize MSKOUT
  MSKOUT(MSKOUT==0)=NaN; 
  MSKOUT(~isnan(MSKOUT))=1;
end;

dxC=mygrid.DXC; dyC=mygrid.DYC; 
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);
nn=length(kk);
NN=tmp1; NN(kk)=[1:nn]; %NN=convert2array(NN); 

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_conv(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);%right hand side of FLD=fld 
mskFreeze_vec=convert2array(mskFreeze);
FLD_vec(find(mskFreeze_vec==1))=0;%right hand side of diffusion equation
FLD_vec=FLD_vec(kk);

INV_op=1-mskFreeze_vec(kk);
INV_op=sparse([1:nn],[1:nn],INV_op,nn,nn);%identity where mskFreeze is 0 
INV_op=INV_op+dFLDdt_op;%add diffusion operator where mskFreeze is 1

warning('off','MATLAB:nearlySingularMatrix');
warning('off','MATLAB:singularMatrix');
INV_vec=INV_op\FLD_vec;%solve
warning('on','MATLAB:nearlySingularMatrix');
warning('on','MATLAB:singularMatrix');

INV_fld=convert2array(mskOut);
INV_fld(find(~isnan(INV_fld)))=INV_vec;
FLD=convert2array(INV_fld);%reformat

if doMSKOUT; FLD=FLD.*MSKOUT; end;%mask problematic regions


1	gforget	1.5	function [FLD]=diffsmooth2D_extrap_inv(fld,mskOut,varargin);
2			%object: extrapolate an incomplete field to create a full field, by
3			% solving for a diffusion equation equilibrium state.
4			%inputs: fld incomplete field of interest (masked with NaN)
5			% mskOut land mask (1s and NaNs) for the full field (output)
6			%output: FLD full field
7	gforget	1.7	%optional: doFormMatrix (1 by default)
8			% if set to 1 then compute the dFLDdt_op matrix
9			% if set to 0 then use precomputed one (global dFLDdt_op)
10			% doMSKOUT (1 by default)
11			% if set to 1 then first detect closed regions lacking
12			% data points to extrapolate from (e.g. abyssal canions)
13			% and restrict mskOut to MSKOUT accordingly
14			% if set to 0 then omit this step
15			%principle :
16			%in points where mskFreeze is 1 solve diffusion equation
17			%in points where mskFreeze is 0 solve the trivial FLD=fld
18	gforget	1.1
19	gforget	1.7	gcmfaces_global;
20
21			if nargin==3; doFormMatrix=varargin{1}; else; doFormMatrix=1; end;
22			if nargin==4; doMSKOUT=varargin{2}; else; doMSKOUT=1; end;
23
24			if doMSKOUT;
25			%problematic regions will show MSKOUT==0;
26			tmp1=1+0*fld; [MSKOUT]=diffsmooth2D_extrap_inv(tmp1,mskOut,doFormMatrix,0);
27			%if no problematic region then no need to recompute dFLDdt_op
28			if sum(MSKOUT==0)==0; doFormMatrix=0; end;
29			%finalize MSKOUT
30			MSKOUT(MSKOUT==0)=NaN;
31			MSKOUT(~isnan(MSKOUT))=1;
32			end;
33	gforget	1.1
34			dxC=mygrid.DXC; dyC=mygrid.DYC;
35			rA=mygrid.RAC;
36
37			dxCsm=dxC; dyCsm=dyC;
38			mskFreeze=fld; mskFreeze(find(~isnan(mskFreeze)))=0; mskFreeze(find(isnan(mskFreeze)))=1;
39
40			%check for domain edge points where no exchange is possible:
41			tmp1=mskOut; tmp1(:)=1; tmp2=exch_T_N(tmp1);
42			for iF=1:mskOut.nFaces;
43			tmp3=mskOut{iF}; tmp4=tmp2{iF};
44			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);
45			if ~isempty(find(isnan(tmp4)&~isnan(tmp3))); fprintf('warning: mask was modified\n'); end;
46			tmp3(isnan(tmp4))=NaN; mskOut{iF}=tmp3;
47			end;
48
49			%put 0 first guess if needed and switch land mask:
50			fld(find(isnan(fld)))=0; fld=fld.*mskOut;
51
52			%scale the diffusive operator:
53			tmp0=dxCsm./dxC; tmp0(isnan(mskOut))=NaN; tmp00=nanmax(tmp0);
54			tmp0=dyCsm./dyC; tmp0(isnan(mskOut))=NaN; tmp00=max([tmp00 nanmax(tmp0)]);
55			smooth2D_nbt=tmp00;
56			smooth2D_nbt=ceil(1.12smooth2D_nbt^2);
57
58			smooth2D_dt=1;
59			smooth2D_T=smooth2D_nbt*smooth2D_dt;
60			smooth2D_Kux=dxCsm.*dxCsm/smooth2D_T/2;
61			smooth2D_Kvy=dyCsm.*dyCsm/smooth2D_T/2;
62
63			%form matrix problem:
64			tmp1=convert2array(mskOut);
65			kk=find(~isnan(tmp1));
66	gforget	1.4	KK=tmp1; KK(kk)=kk; KK=convert2array(KK);
67	gforget	1.1	nn=length(kk);
68	gforget	1.4	NN=tmp1; NN(kk)=[1:nn]; %NN=convert2array(NN);
69	gforget	1.1
70	gforget	1.2	global dFLDdt_op;
71			if doFormMatrix==1;
72
73	gforget	1.1	dFLDdt_op=sparse([],[],[],nn,nn,nn*5);
74
75			for iF=1:fld.nFaces; for ii=1:3; for jj=1:3;
76			FLDones=fld; FLDones(find(~isnan(fld)))=0;
77			FLDones{iF}(ii:3:end,jj:3:end)=1;
78			FLDones(find(isnan(fld)))=NaN;
79
80			FLDkkFROMtmp=fld; FLDkkFROMtmp(find(~isnan(fld)))=0;
81			FLDkkFROMtmp{iF}(ii:3:end,jj:3:end)=KK{iF}(ii:3:end,jj:3:end);
82			FLDkkFROMtmp(find(isnan(fld)))=0;
83
84			FLDkkFROM=exch_T_N(FLDkkFROMtmp);
85	gforget	1.3	FLDkkFROM(find(isnan(FLDkkFROM)))=0;
86	gforget	1.1	for iF2=1:fld.nFaces;
87			tmp1=FLDkkFROM{iF2}; tmp2=zeros(size(tmp1)-2);
88			for ii2=1:3; for jj2=1:3; tmp2=tmp2+tmp1(ii2:end-3+ii2,jj2:end-3+jj2); end; end;
89			FLDkkFROM{iF2}=tmp2;
90			end;
91			%clear FLDkkFROMtmp;
92
93			[dTdxAtU,dTdyAtV]=calc_T_grad(FLDones,0);
94			tmpU=dTdxAtU.*smooth2D_Kux;
95			tmpV=dTdyAtV.*smooth2D_Kvy;
96	gforget	1.8	[fldDIV]=calc_UV_conv(tmpU,tmpV);
97	gforget	1.1	dFLDdt=smooth2D_dt*fldDIV./rA;
98			dFLDdt=dFLDdt.*mskFreeze;
99
100			dFLDdt=convert2array(dFLDdt); FLDkkFROM=convert2array(FLDkkFROM); FLDkkTO=convert2array(KK);
101			tmp1=find(dFLDdt~=0&~isnan(dFLDdt));
102			dFLDdt=dFLDdt(tmp1); FLDkkFROM=FLDkkFROM(tmp1); FLDkkTO=FLDkkTO(tmp1);
103			dFLDdt_op=dFLDdt_op+sparse(NN(FLDkkTO),NN(FLDkkFROM),dFLDdt,nn,nn);
104
105			end; end; end;
106
107	gforget	1.2	end;%if doFormMatrix==1;
108
109	gforget	1.1	%figure; spy(dFLDdt_op);
110
111	gforget	1.7	FLD_vec=convert2array(fld);%right hand side of FLD=fld
112	gforget	1.1	mskFreeze_vec=convert2array(mskFreeze);
113	gforget	1.7	FLD_vec(find(mskFreeze_vec==1))=0;%right hand side of diffusion equation
114	gforget	1.1	FLD_vec=FLD_vec(kk);
115
116			INV_op=1-mskFreeze_vec(kk);
117	gforget	1.7	INV_op=sparse([1:nn],[1:nn],INV_op,nn,nn);%identity where mskFreeze is 0
118			INV_op=INV_op+dFLDdt_op;%add diffusion operator where mskFreeze is 1
119
120	gforget	1.9	warning('off','MATLAB:nearlySingularMatrix');
121			warning('off','MATLAB:singularMatrix');
122	gforget	1.7	INV_vec=INV_op\FLD_vec;%solve
123	gforget	1.9	warning('on','MATLAB:nearlySingularMatrix');
124			warning('on','MATLAB:singularMatrix');
125
126	gforget	1.1	INV_fld=convert2array(mskOut);
127			INV_fld(find(~isnan(INV_fld)))=INV_vec;
128	gforget	1.7	FLD=convert2array(INV_fld);%reformat
129	gforget	1.1
130	gforget	1.7	if doMSKOUT; FLD=FLD.*MSKOUT; end;%mask problematic regions
131	gforget	1.1
132