profilesMatlabProcessing/profiles_devel/profiles_subgrid_stats_assemble.m

function [myFld,myFld1]=profiles_subgrid_stats_assemble(choiceVar,lev,choiceDir,choiceFld);
%... one level at a time

doSave=1
choiceGrid='v4';
%choiceDir=2;
%choiceFld=2;

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

gcmfaces_global;

%listSGN=[0 8 4 2 1];%for cs24
if strcmp(choiceGrid,'cs96');
    %full suite for cs96
    listSGN=[0 48 32 24 16 12 8 6 4 3 2 1];
    dirGrid='/net/weddell/raid3/gforget/grids/gridCompleted/cube_FM/';
    %dirGrid='/Users/gforget/mywork/projects_inprogress/2012mayInputs/insitu/processed/';
    grid_load_native([dirGrid 'cube_96/'],6);
else;
    ;%for ecco v4
    listSGN=[0 45 30 18 15 10 9 6 5 3 2 1]
    dirGrid='/net/nares/raid10/gforget/2012julyIters/GRID/';
    %dirGrid='/net/weddell/raid3/gforget/ecco_v4/GRID/';
    %dirGrid='GRID/';
    %grid_load(dirGrid,5,'compact');
end;

global atlasT atlasS;
if isempty(atlasT);
    dirAtlas='/net/weddell/raid3/gforget/ecco_v4/input_files/';
    %dirAtlas='/Users/gforget/mywork/projects_inprogress/2012mayInputs/atlas/bin/';
    %atlasT=v4_read_bin([dirAtlas 'T_OWPv1_Y_eccollc_90x50.bin']);
    %atlasS=v4_read_bin([dirAtlas 'S_OWPv1_Y_eccollc_90x50.bin']);
end;

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

if choiceDir==1;
dirOut='/net/nares/raid10/gforget/2013mayInputs/insitu/stats_depth/';
elseif choiceDir==2;
dirOut='/net/nares/raid10/gforget/2012julyIters/ecco_it0003/mode_water_analyses/profiles_subgrid_r3it3/'
elseif choiceDir==3;
dirOut='/net/nares/raid10/gforget/2012julyIters/ecco_it0003/mode_water_analyses/profiles_subgrid_argo/'
end;

%dirOut='~/mywork/projects_inprogress/2012mayInputs/insitu/tmp/';

%if result was not completed, then skip:
test0=dir([dirOut choiceGrid '_' choiceVar '_k' num2str(lev) '_' 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 choiceGrid '_' choiceVar '_k' num2str(lev) '_' 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=mygrid.mskC(:,:,lev);

%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' choiceVar '(:,:,lev);']);
    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' choiceVar '(:,:,lev);']);
    [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 choiceGrid '_' choiceVar '_k' num2str(lev) '_' suff '.mat myFld myFld1;']);
end;

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