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;
dirData='./';

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

gcmfaces_global;

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

global atlas;
test0=isempty(atlas);%no previously loaded atlas?
if ~test0; test0=~strcmp(atlas.coord,COORD); end;%same coordinate as before?
if test0;
    atlas=[];
    atlas.coord=COORD;
    %
    dirAtlases=[myenv.gcmfaces_dir 'sample_input/OCCAetcONv4GRID/'];
    if ~strcmp(atlas.coord,'depth');
      eval(['load ' myenv.gcmfaces_dir 'gcmfaces_devel/RCsig0.mat RCsig0;']);
      atlas.RC=RCsig0;
      txt=['_' atlas.coord];
    else;
      atlas.RC=-mygrid.RC;
      txt='';
    end;
    nr=length(atlas.RC);
    %
    if ~strcmp(atlas.coord,'depth');
      tmp1=read2memory([dirAtlases 'D_OWPv1_M_eccollc_90x50' txt '.bin'],[90 1170 nr 12]);
      tmp1(tmp1==0)=NaN;
      atlas.D=convert2gcmfaces(tmp1);
      %
      tmp1=1*(sum(~isnan(atlas.D),4)>2);
      tmp1(tmp1==0)=NaN;
      %
      tmp2=nanmean(atlas.D,4);                  
      tmp2=1*(tmp2<2000); 
      tmp2(tmp2==0)=NaN;
      %
      atlas.mskC=tmp1.*tmp2;
    else;
      atlas.mskC=mygrid.mskC;
    end;
    %
    if 0;%only needed for slanted diffusion
      tmp1=read2memory([dirAtlases 'T_OWPv1_M_eccollc_90x50' txt '.bin'],[90 1170 nr 12]);
      tmp1(tmp1==0)=NaN;
      atlas.T=convert2gcmfaces(tmp1);
      %
      tmp1=read2memory([dirAtlases 'S_OWPv1_M_eccollc_90x50' txt '.bin'],[90 1170 nr 12]);
      tmp1(tmp1==0)=NaN;
      atlas.S=convert2gcmfaces(tmp1);
    end;
end;

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

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