profilesMatlabProcessing/profiles_devel/profiles_subgrid_stats_assemble.m

function [myFld,myFld1]=profiles_subgrid_stats_assemble(VV,KK,choiceFld,COORD);
%KK level choice
%VV variable choice
%choiceFld stat choice
%COORD is the coordinate type

doSave=1;
dirOut='./';

%========= PART 1 : load grid & atlases ========

gcmfaces_global;

listSGN=[0 45 30 18 15 10 9 6 5 3 2 1]

global atlas;
if isempty(atlas);
    dirAtlases=[myenv.gcmfaces_dir 'sample_input/OCCAetcONv4GRID/'];
    atlas.coord='sig0';
    eval(['load ' myenv.gcmfaces_dir 'gcmfaces_devel/RCsig0.mat RCsig0;']);
    atlas.RC=RCsig0;
    nr=length(atlas.RC);
    %
    tmp1=read2memory([dirAtlases 'T_OWPv1_M_eccollc_90x50_' atlas.coord '.bin'],[90 1170 nr 12]);
    tmp1(tmp1==0)=NaN;
    atlas.T=convert2gcmfaces(tmp1);
    %
    tmp1=1*(sum(~isnan(atlas.T),4)>2);
    tmp1(tmp1==0)=NaN;
    atlas.mskC=tmp1;
    %
    tmp1=read2memory([dirAtlases 'S_OWPv1_M_eccollc_90x50_' atlas.coord '.bin'],[90 1170 nr 12]);
    tmp1(tmp1==0)=NaN;
    atlas.S=convert2gcmfaces(tmp1);
    %
    tmp1=read2memory([dirAtlases 'D_OWPv1_M_eccollc_90x50_' atlas.coord '.bin'],[90 1170 nr 12]);
    tmp1(tmp1==0)=NaN;
    atlas.D=convert2gcmfaces(tmp1);
    %
    tmp1=read2memory([dirAtlases atlas.coord '_OWPv1_M_eccollc_90x50.bin'],[90 1170 50 12]);
    tmp1(tmp1==0)=NaN;
    atlas.sig=convert2gcmfaces(tmp1);
end;

%========= PART 2 : load and average estimates of std ========

%if result was not completed, then skip:
test0=dir([dirOut VV '_k' num2str(KK) '_' num2str(1) '.mat']);
if isempty(test0); myFld=NaN*mygrid.RAC; myFld1=myFld; return; end;

listStats={'prc10','mea','med','prc90','std','iqr','mad'};

myWeightPower=4

for sgn=listSGN;
    eval(['load ' dirOut VV '_k' num2str(KK) '_' num2str(sgn) '.mat myStat;']);
    %"bootstrap"
    kk=find(listSGN==sgn);
    if kk==1;
        myFld=myStat;
        w=myStat.nb/(sqrt(90*1170)^myWeightPower);
        myFld.nb=w;
        for ff=1:length(listStats);
          eval(['myFld.' listStats{ff} '=myStat.' listStats{ff} '.*w;']);
        end;
    else;
        w=myStat.nb/(sgn^myWeightPower);
        myFld.nb=myFld.nb+w;
        for ff=1:length(listStats);
          eval(['myFld.' listStats{ff} '=myFld.' listStats{ff} '+myStat.' listStats{ff} '.*w;']);
        end;
    end;
end;

%THIS WAS A BUG : myFld.nb=myFld.nb+w;
for ff=1:length(listStats);
  eval(['myFld.' listStats{ff} '=(myFld.' listStats{ff} './myFld.nb);']);
end;
myFld.msk=atlas.mskC(:,:,KK);

%original value & "local" value forcing:
suff=listStats{choiceFld};
eval(['myFld1=myFld.msk.*myFld.' suff ';']);

%========= PART 3 : smoothing setup ========

if 0;%simple smoothing, which does not account for no. of obs
    eval(['myFld.mean=myFld.msk.*atlas' VV '(:,:,KK);']);
    myFld.sm0=diffsmooth2D(myFld1,mygrid.DXC*3,mygrid.DYC*3);
    dxy=3*sqrt(mygrid.RAC);
    myFld.sm1=diffsmooth2D(myFld1,dxy,dxy);
    myFld.sm2=diffsmooth2Drotated(myFld1,dxy,dxy/10,myFld.mean);
end;


%scale the diffusive operator:
dxLarge=3*sqrt(mygrid.RAC);
dxSmall=0.1*dxLarge;

%time scale:
tmp0=dxLarge./mygrid.DXC; tmp0(isnan(myFld1))=NaN; tmp00=nanmax(tmp0);
tmp0=dxLarge./mygrid.DYC; tmp0(isnan(myFld1))=NaN; tmp00=max([tmp00 nanmax(tmp0)]);
nbt=tmp00;
nbt=ceil(1.1*2*nbt^2);

dt=1;
T=nbt*dt;

%build diffusion operator:
kLarge=dxLarge.*dxLarge/T/2;
kSmall=dxSmall.*dxSmall/T/2;

if 1;%isotropic diffusion, rather than slanted diffusion
    Kux=dxLarge.*dxLarge/T/2;
    Kvy=dxLarge.*dxLarge/T/2;
    Kuy=[]; Kvx=[];
else;%slanted diffusion
    eval(['myFld.mean=myFld.msk.*atlas' VV '(:,:,KK);']);
    [Kux,Kuy,Kvx,Kvy]=diffrotated(kLarge,kSmall,myFld.mean);
end;

%finalize diffusion/smoothing problem set-up:
myOp.dt=1;
% myOp.nbt=nbt;
myOp.eps=1e-3;
myOp.Kux=Kux;
myOp.Kuy=Kuy;
myOp.Kvx=Kvx;
myOp.Kvy=Kvy;

%========= PART 4 : relaxation term setup ========

%1) set relaxation strength: (local <-> smoother)
%---------------------------

%use the myFld.nb index, modified as follows
w=myFld.nb;
%I do a linear transiton in log10 
w=log10(w);
%by mapping [-2 2] to [2 -1]
w=(-1-3*(w-2)/(2+2));
%go back to original units (~nb obs) and scale by nbt (nbt = 1 smoother)
w=nbt*exp( w*log(10) );
%ensure stability
w(w<1)=1;
%enforce minimum forcing
w(w>1e3*nbt)=1e3*nbt;
if 0;
    %figureL; m_map_gcmfaces(log10(myFld.nb),0,{'myCaxis',[-4 3]});
    figureL; m_map_gcmfaces(log10(w/nbt),0,{'myCaxis',[-2 2]}); return;
end;
myOp.tau=w*myOp.dt;

% myOp.tau=0.5*myOp.dt;
% myOp.tau=nbt*myOp.dt;
% myOp.tau=nbt;

%2) set relaxation field: ("local" value)
%------------------------

fldRelax=myFld1;


%========= PART 5 : resolve smoothing/relaxation problem ========
%
%   here we integrate to a balance between 
%       "local" value (relaxation term) 
%       vs  smoothing (diffusion)

myFld=gcmfaces_timestep(myOp,myFld1,fldRelax);

%plot / save result:
%===================

if 0; 
  figureL; m_map_gcmfaces(log10(myFld),0,{'myCaxis',[-1.5 0.5]});
end;

if doSave;
  eval(['save ' dirOut VV '_k' num2str(KK) '_' suff '.mat myFld myFld1;']);
end;

1	gforget	1.3	function [myFld,myFld1]=profiles_subgrid_stats_assemble(VV,KK,choiceFld,COORD);
2			%KK level choice
3			%VV variable choice
4			%choiceFld stat choice
5			%COORD is the coordinate type
6	gforget	1.1
7	gforget	1.3	doSave=1;
8			dirOut='./';
9	gforget	1.1
10			%========= PART 1 : load grid & atlases ========
11
12			gcmfaces_global;
13
14	gforget	1.3	listSGN=[0 45 30 18 15 10 9 6 5 3 2 1]
15
16			global atlas;
17			if isempty(atlas);
18			dirAtlases=[myenv.gcmfaces_dir 'sample_input/OCCAetcONv4GRID/'];
19			atlas.coord='sig0';
20			eval(['load ' myenv.gcmfaces_dir 'gcmfaces_devel/RCsig0.mat RCsig0;']);
21			atlas.RC=RCsig0;
22			nr=length(atlas.RC);
23			%
24			tmp1=read2memory([dirAtlases 'T_OWPv1_M_eccollc_90x50_' atlas.coord '.bin'],[90 1170 nr 12]);
25			tmp1(tmp1==0)=NaN;
26			atlas.T=convert2gcmfaces(tmp1);
27			%
28			tmp1=1*(sum(~isnan(atlas.T),4)>2);
29			tmp1(tmp1==0)=NaN;
30			atlas.mskC=tmp1;
31			%
32			tmp1=read2memory([dirAtlases 'S_OWPv1_M_eccollc_90x50_' atlas.coord '.bin'],[90 1170 nr 12]);
33			tmp1(tmp1==0)=NaN;
34			atlas.S=convert2gcmfaces(tmp1);
35			%
36			tmp1=read2memory([dirAtlases 'D_OWPv1_M_eccollc_90x50_' atlas.coord '.bin'],[90 1170 nr 12]);
37			tmp1(tmp1==0)=NaN;
38			atlas.D=convert2gcmfaces(tmp1);
39			%
40			tmp1=read2memory([dirAtlases atlas.coord '_OWPv1_M_eccollc_90x50.bin'],[90 1170 50 12]);
41			tmp1(tmp1==0)=NaN;
42			atlas.sig=convert2gcmfaces(tmp1);
43	gforget	1.1	end;
44
45			%========= PART 2 : load and average estimates of std ========
46
47			%if result was not completed, then skip:
48	gforget	1.3	test0=dir([dirOut VV '_k' num2str(KK) '_' num2str(1) '.mat']);
49	gforget	1.1	if isempty(test0); myFld=NaN*mygrid.RAC; myFld1=myFld; return; end;
50
51	gforget	1.2	listStats={'prc10','mea','med','prc90','std','iqr','mad'};
52	gforget	1.1
53			myWeightPower=4
54
55			for sgn=listSGN;
56	gforget	1.3	eval(['load ' dirOut VV '_k' num2str(KK) '_' num2str(sgn) '.mat myStat;']);
57	gforget	1.1	%"bootstrap"
58			kk=find(listSGN==sgn);
59			if kk==1;
60			myFld=myStat;
61			w=myStat.nb/(sqrt(90*1170)^myWeightPower);
62			myFld.nb=w;
63	gforget	1.2	for ff=1:length(listStats);
64			eval(['myFld.' listStats{ff} '=myStat.' listStats{ff} '.*w;']);
65			end;
66	gforget	1.1	else;
67			w=myStat.nb/(sgn^myWeightPower);
68	gforget	1.2	myFld.nb=myFld.nb+w;
69			for ff=1:length(listStats);
70			eval(['myFld.' listStats{ff} '=myFld.' listStats{ff} '+myStat.' listStats{ff} '.*w;']);
71	gforget	1.1	end;
72			end;
73			end;
74
75			%THIS WAS A BUG : myFld.nb=myFld.nb+w;
76	gforget	1.2	for ff=1:length(listStats);
77			eval(['myFld.' listStats{ff} '=(myFld.' listStats{ff} './myFld.nb);']);
78			end;
79	gforget	1.3	myFld.msk=atlas.mskC(:,:,KK);
80	gforget	1.1
81			%original value & "local" value forcing:
82	gforget	1.2	suff=listStats{choiceFld};
83			eval(['myFld1=myFld.msk.*myFld.' suff ';']);
84	gforget	1.1
85			%========= PART 3 : smoothing setup ========
86
87			if 0;%simple smoothing, which does not account for no. of obs
88	gforget	1.3	eval(['myFld.mean=myFld.msk.*atlas' VV '(:,:,KK);']);
89	gforget	1.1	myFld.sm0=diffsmooth2D(myFld1,mygrid.DXC3,mygrid.DYC3);
90			dxy=3*sqrt(mygrid.RAC);
91			myFld.sm1=diffsmooth2D(myFld1,dxy,dxy);
92			myFld.sm2=diffsmooth2Drotated(myFld1,dxy,dxy/10,myFld.mean);
93			end;
94
95
96			%scale the diffusive operator:
97			dxLarge=3*sqrt(mygrid.RAC);
98			dxSmall=0.1*dxLarge;
99
100			%time scale:
101			tmp0=dxLarge./mygrid.DXC; tmp0(isnan(myFld1))=NaN; tmp00=nanmax(tmp0);
102			tmp0=dxLarge./mygrid.DYC; tmp0(isnan(myFld1))=NaN; tmp00=max([tmp00 nanmax(tmp0)]);
103			nbt=tmp00;
104			nbt=ceil(1.12nbt^2);
105
106			dt=1;
107			T=nbt*dt;
108
109			%build diffusion operator:
110			kLarge=dxLarge.*dxLarge/T/2;
111			kSmall=dxSmall.*dxSmall/T/2;
112
113			if 1;%isotropic diffusion, rather than slanted diffusion
114			Kux=dxLarge.*dxLarge/T/2;
115			Kvy=dxLarge.*dxLarge/T/2;
116			Kuy=[]; Kvx=[];
117			else;%slanted diffusion
118	gforget	1.3	eval(['myFld.mean=myFld.msk.*atlas' VV '(:,:,KK);']);
119	gforget	1.1	[Kux,Kuy,Kvx,Kvy]=diffrotated(kLarge,kSmall,myFld.mean);
120			end;
121
122			%finalize diffusion/smoothing problem set-up:
123			myOp.dt=1;
124			% myOp.nbt=nbt;
125			myOp.eps=1e-3;
126			myOp.Kux=Kux;
127			myOp.Kuy=Kuy;
128			myOp.Kvx=Kvx;
129			myOp.Kvy=Kvy;
130
131			%========= PART 4 : relaxation term setup ========
132
133			%1) set relaxation strength: (local <-> smoother)
134			%---------------------------
135
136			%use the myFld.nb index, modified as follows
137			w=myFld.nb;
138			%I do a linear transiton in log10
139			w=log10(w);
140			%by mapping [-2 2] to [2 -1]
141			w=(-1-3*(w-2)/(2+2));
142			%go back to original units (~nb obs) and scale by nbt (nbt = 1 smoother)
143			w=nbtexp( wlog(10) );
144			%ensure stability
145			w(w<1)=1;
146			%enforce minimum forcing
147			w(w>1e3nbt)=1e3nbt;
148			if 0;
149			%figureL; m_map_gcmfaces(log10(myFld.nb),0,{'myCaxis',[-4 3]});
150			figureL; m_map_gcmfaces(log10(w/nbt),0,{'myCaxis',[-2 2]}); return;
151			end;
152			myOp.tau=w*myOp.dt;
153
154			% myOp.tau=0.5*myOp.dt;
155			% myOp.tau=nbt*myOp.dt;
156			% myOp.tau=nbt;
157
158			%2) set relaxation field: ("local" value)
159			%------------------------
160
161			fldRelax=myFld1;
162
163
164			%========= PART 5 : resolve smoothing/relaxation problem ========
165			%
166			% here we integrate to a balance between
167			% "local" value (relaxation term)
168			% vs smoothing (diffusion)
169
170			myFld=gcmfaces_timestep(myOp,myFld1,fldRelax);
171
172			%plot / save result:
173			%===================
174
175			if 0;
176			figureL; m_map_gcmfaces(log10(myFld),0,{'myCaxis',[-1.5 0.5]});
177			end;
178
179			if doSave;
180	gforget	1.3	eval(['save ' dirOut VV '_k' num2str(KK) '_' suff '.mat myFld myFld1;']);
181	gforget	1.1	end;
182