matlab_class/ecco_v4/cost_seaicearea.m

function []=cost_seaicearea(dirModel,dirMat,doComp,dirTex,nameTex);
%object:        compute cost function term for sea ice data
%inputs:        dimodel is the model directory
%               dirMat is the directory where diagnozed .mat files will be saved
%                     -> set it to '' to use the default [dirModel 'mat/']
%               doComp is a switch (1->compute; 0->display)
%optional:      dirTex is the directory where tex and figures files are created
%                 (if not specified then display all results to screen instead)
%               nameTex is the tex file name (default : 'myPlots')

if isempty(dirMat); dirMat=[dirModel 'mat' filesep]; else; dirMat=[dirMat filesep]; end;
if isempty(dir(dirMat));     mkdir([dirMat]); end;

%determine if and where to create tex and figures files
dirMat=[dirMat filesep];
if isempty(who('dirTex'));
  addToTex=0;
else;
  if ~ischar(dirTex); error('mis-specified dirTex'); end;
  addToTex=1;
  if isempty(who('nameTex')); nameTex='myPlots'; end;
  fileTex=[dirTex filesep nameTex '.tex'];
end;

%grid, params and inputs
gcmfaces_global; global myparms;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
if isfield(myparms,'yearsta'); yearsta=myparms.yearsta; yearend=myparms.yearend;
else; yearsta=1992; yearend=2011;
end;

[lonPairs,latPairs,names] = line_greatC_TUV_MASKS_core2_antarctic;
lonLims=[lonPairs(1:5,1);lonPairs(1,1)];

if doComp;
%grid, params and inputs
fld_err=ones(90,1170);
fld_err=convert2gcmfaces(fld_err);
fld_w=fld_err.^-2;

if isempty(dir([dirModel 'barfiles']));
  dirEcco=dirModel;
else;
  dirEcco=[dirModel 'barfiles' filesep];
end;

nameData='nsidc79_monthly_';
if ~isempty(dir([dirModel nameData '*']));
  dirData=dirModel;
elseif ~isempty(dir([dirModel 'inputfiles/' nameData '*']));
  dirData=[dirModel 'inputfiles/'];
else;
  error(['could not find ' nameData '*']);
end;

fileModel=dir([dirEcco 'm_siv4_area*data']); 
fileModel=fileModel.name;

nyears=yearend-yearsta+1;
nmonths=12*nyears;

%misfits :
fld_dif=convert2gcmfaces(NaN*ones(90,90*13,nmonths));
fld_nsidc=convert2gcmfaces(NaN*ones(90,90*13,nmonths));

%monthly mean climatology :
climMod=convert2gcmfaces(zeros(90,90*13,12));
climObs=convert2gcmfaces(zeros(90,90*13,12));
climMsk=convert2gcmfaces(zeros(90,90*13,12));
climNb=convert2gcmfaces(zeros(90,90*13,12));

%monthly integrals :
IceAreaNorthMod=NaN*zeros(1,nmonths);
IceAreaNorthObs=NaN*zeros(1,nmonths);

IceAreaSouthMod=NaN*zeros(6,nmonths);
IceAreaSouthObs=NaN*zeros(6,nmonths);

%computational loop :
for ycur=yearsta:yearend;
 tic;
 for mcur=1:12;
  mm=(ycur-yearsta)*12+mcur;

  fld_dat=read_bin([dirData nameData num2str(ycur)],mcur,0);

  fld_dat=fld_dat.*mygrid.mskC(:,:,1);%land mask
  fld_dat(find(fld_dat<-99))=NaN;%missing data
  msk=1+0*fld_dat;%combined mask

  fld_mod=read_bin([dirEcco fileModel],mm,0);
  fld_mod=fld_mod.*msk;%mask consistent with fld_dat

  %misfits : 
  fld_dif(:,:,mm)=fld_mod-fld_dat;
  fld_nsidc(:,:,mm)=fld_dat;

  %climatology :
  tmp1=msk; tmp1(isnan(tmp1))=0; climMsk(:,:,mcur)=climMsk(:,:,mcur)+tmp1;
  tmp1=fld_mod; tmp1(isnan(tmp1))=0; climMod(:,:,mcur)=climMod(:,:,mcur)+tmp1;
  tmp1=fld_dat; tmp1(isnan(tmp1))=0; climObs(:,:,mcur)=climObs(:,:,mcur)+tmp1;

  %integrals :
  fld=fld_mod.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
  fld=fld_dat.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);

  fld=fld_mod.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthMod(1,mm)=nansum(fld);
  fld=fld_dat.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthObs(1,mm)=nansum(fld);

    for kk=1:5; 
        tmpmsk=0.*mygrid.XC;
        if lonLims(kk+1) > lonLims(kk)
          tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC < lonLims(kk+1)))=1.;
        else 
          tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC <= 180.))=1.;
          tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < lonLims(kk+1)))=1.;
        end
        tmpmsk=tmpmsk.*(mygrid.YC<0);
        %
        fld=fld_mod.*mygrid.RAC.*tmpmsk;
        IceAreaSouthMod(kk+1,mm)=nansum(fld);
        %
        fld=fld_dat.*mygrid.RAC.*tmpmsk;
        IceAreaSouthObs(kk+1,mm)=nansum(fld);
    end

 end;
 toc;
end;

%misfits :
mis_rms=sqrt(nanmean(fld_dif.^2,3));
obs_std=nanstd(fld_nsidc,[],3);
mod_std=nanstd(fld_nsidc+fld_dif,[],3);

%climatology :
for mcur=1:12;
  tmp1=climMsk(:,:,mcur); tmp1(tmp1==0)=NaN;
  climNb(:,:,mcur)=tmp1;
  climMod(:,:,mcur)=climMod(:,:,mcur)./tmp1;
  climObs(:,:,mcur)=climObs(:,:,mcur)./tmp1;
end;
clear climMsk;

if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
eval(['save ' dirMat '/cost/cost_seaicearea.mat fld_err mis_rms obs_std mod_std IceArea* clim*;']);

else;%display previously computed results

if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;

eval(['load ' dirMat '/cost_seaicearea.mat;']);

%variance maps:
figure; m_map_gcmfaces(mis_rms,0,{'myCaxis',[0:0.1:1.]});
myCaption={'modeled-observed rms -- sea ice concentration'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure; m_map_gcmfaces(obs_std,0,{'myCaxis',[0:0.1:1.]});
myCaption={'observed std -- sea ice concentration'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure; m_map_gcmfaces(mod_std,0,{'myCaxis',[0:0.1:1.]});
myCaption={'modelled std -- sea ice concentration'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

%arctic maps
figureL; 
subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),2,{'myCaxis',[0:0.1:1.]});
subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),2,{'myCaxis',[0:0.1:1.]});
subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),2,{'myCaxis',[0:0.1:1.]});
subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),2,{'myCaxis',[0:0.1:1.]});
myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

%southern ocean maps
figureL;
subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),3,{'myCaxis',[0:0.1:1.]});
subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),3,{'myCaxis',[0:0.1:1.]});
subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),3,{'myCaxis',[0:0.1:1.]});
subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),3,{'myCaxis',[0:0.1:1.]});
myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

%northern/southern integrals
ny=length(IceAreaNorthMod)/12;
yy=yearsta+[0:ny-1];

figureL;
subplot(1,2,1);
plot(yy,IceAreaNorthMod(3:12:end),'LineWidth',2); hold on;
plot(yy,IceAreaNorthObs(3:12:end),'r','LineWidth',2);
plot(yy,IceAreaNorthMod(9:12:end),'LineWidth',2);
plot(yy,IceAreaNorthObs(9:12:end),'r','LineWidth',2);
axis([yearsta yearsta+ny-1 0 20e12]);
ylabel('m^2'); title('Northern Hemisphere');
subplot(1,2,2);
plot(yy,IceAreaSouthMod(1,3:12:end),'LineWidth',2); hold on;
plot(yy,IceAreaSouthObs(1,3:12:end),'r','LineWidth',2);
plot(yy,IceAreaSouthMod(1,9:12:end),'LineWidth',2);
plot(yy,IceAreaSouthObs(1,9:12:end),'r','LineWidth',2);
axis([yearsta yearsta+ny-1 0 20e12]);
ylabel('m^2'); title('Southern Hemisphere');

myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in March and September',...
           'in Northern Hemisphere (left) and Southern Hemisphere (right)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

%southern basin integrals

for mm=[3 9];
  figureL;
  for ii=1:6;
    subplot(3,2,ii);
    if ii==1; iiTxt='Entire Southern Ocean';
    else; iiTxt=sprintf('%dE to %dE',lonLims(ii-1),lonLims(ii));
    end;
    plot(yy,IceAreaSouthMod(ii,mm:12:end),'LineWidth',2); hold on;
    plot(yy,IceAreaSouthObs(ii,mm:12:end),'r','LineWidth',2);
    aa=axis; aa(1:2)=[yearsta yearsta+ny-1]; axis(aa);
    ylabel('m^2'); title(iiTxt);
  end;
  if mm==3; mmTxt='March'; elseif mm==9; mmTxt='September'; else; '???'; end;
  myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in ',mmTxt,' per Southern Ocean sector'};
  if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
end;

end;


1	gforget	1.2	function []=cost_seaicearea(dirModel,dirMat,doComp,dirTex,nameTex);
2	heimbach	1.1	%object: compute cost function term for sea ice data
3			%inputs: dimodel is the model directory
4			% dirMat is the directory where diagnozed .mat files will be saved
5			% -> set it to '' to use the default [dirModel 'mat/']
6			% doComp is a switch (1->compute; 0->display)
7			%optional: dirTex is the directory where tex and figures files are created
8			% (if not specified then display all results to screen instead)
9	gforget	1.2	% nameTex is the tex file name (default : 'myPlots')
10	heimbach	1.1
11	gforget	1.15	if isempty(dirMat); dirMat=[dirModel 'mat' filesep]; else; dirMat=[dirMat filesep]; end;
12	gforget	1.13	if isempty(dir(dirMat)); mkdir([dirMat]); end;
13	heimbach	1.1
14			%determine if and where to create tex and figures files
15	gforget	1.15	dirMat=[dirMat filesep];
16	heimbach	1.1	if isempty(who('dirTex'));
17			addToTex=0;
18			else;
19			if ~ischar(dirTex); error('mis-specified dirTex'); end;
20	gforget	1.2	addToTex=1;
21			if isempty(who('nameTex')); nameTex='myPlots'; end;
22	gforget	1.15	fileTex=[dirTex filesep nameTex '.tex'];
23	heimbach	1.1	end;
24
25	gforget	1.5	%grid, params and inputs
26			gcmfaces_global; global myparms;
27	heimbach	1.1	if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
28			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
29	gforget	1.5	if isfield(myparms,'yearsta'); yearsta=myparms.yearsta; yearend=myparms.yearend;
30			else; yearsta=1992; yearend=2011;
31			end;
32	heimbach	1.1
33	gforget	1.8	[lonPairs,latPairs,names] = line_greatC_TUV_MASKS_core2_antarctic;
34			lonLims=[lonPairs(1:5,1);lonPairs(1,1)];
35
36			if doComp;
37			%grid, params and inputs
38	heimbach	1.1	fld_err=ones(90,1170);
39			fld_err=convert2gcmfaces(fld_err);
40			fld_w=fld_err.^-2;
41
42	gforget	1.14	if isempty(dir([dirModel 'barfiles']));
43			dirEcco=dirModel;
44			else;
45			dirEcco=[dirModel 'barfiles' filesep];
46			end;
47
48	gforget	1.16	nameData='nsidc79_monthly_';
49			if ~isempty(dir([dirModel nameData '*']));
50	gforget	1.14	dirData=dirModel;
51	gforget	1.17	elseif ~isempty(dir([dirModel 'inputfiles/' nameData '*']));
52			dirData=[dirModel 'inputfiles/'];
53	gforget	1.14	else;
54	gforget	1.16	error(['could not find ' nameData '*']);
55	gforget	1.14	end;
56	heimbach	1.1
57	gforget	1.14	fileModel=dir([dirEcco 'm_siv4_area*data']);
58			fileModel=fileModel.name;
59	heimbach	1.1
60	heimbach	1.3	nyears=yearend-yearsta+1;
61			nmonths=12*nyears;
62
63	gforget	1.7	%misfits :
64	heimbach	1.3	fld_dif=convert2gcmfaces(NaNones(90,9013,nmonths));
65	gforget	1.8	fld_nsidc=convert2gcmfaces(NaNones(90,9013,nmonths));
66	heimbach	1.3
67	gforget	1.7	%monthly mean climatology :
68			climMod=convert2gcmfaces(zeros(90,90*13,12));
69			climObs=convert2gcmfaces(zeros(90,90*13,12));
70			climMsk=convert2gcmfaces(zeros(90,90*13,12));
71			climNb=convert2gcmfaces(zeros(90,90*13,12));
72
73			%monthly integrals :
74	gforget	1.4	IceAreaNorthMod=NaN*zeros(1,nmonths);
75			IceAreaNorthObs=NaN*zeros(1,nmonths);
76
77			IceAreaSouthMod=NaN*zeros(6,nmonths);
78			IceAreaSouthObs=NaN*zeros(6,nmonths);
79	heimbach	1.3
80	gforget	1.7	%computational loop :
81	heimbach	1.3	for ycur=yearsta:yearend;
82	heimbach	1.1	tic;
83			for mcur=1:12;
84	heimbach	1.3	mm=(ycur-yearsta)*12+mcur;
85
86	gforget	1.16	fld_dat=read_bin([dirData nameData num2str(ycur)],mcur,0);
87	heimbach	1.1
88	gforget	1.6	fld_dat=fld_dat.*mygrid.mskC(:,:,1);%land mask
89			fld_dat(find(fld_dat<-99))=NaN;%missing data
90			msk=1+0*fld_dat;%combined mask
91	heimbach	1.1
92	gforget	1.16	fld_mod=read_bin([dirEcco fileModel],mm,0);
93	gforget	1.6	fld_mod=fld_mod.*msk;%mask consistent with fld_dat
94	gforget	1.7
95			%misfits :
96	heimbach	1.1	fld_dif(:,:,mm)=fld_mod-fld_dat;
97	gforget	1.8	fld_nsidc(:,:,mm)=fld_dat;
98	heimbach	1.1
99	gforget	1.7	%climatology :
100			tmp1=msk; tmp1(isnan(tmp1))=0; climMsk(:,:,mcur)=climMsk(:,:,mcur)+tmp1;
101			tmp1=fld_mod; tmp1(isnan(tmp1))=0; climMod(:,:,mcur)=climMod(:,:,mcur)+tmp1;
102			tmp1=fld_dat; tmp1(isnan(tmp1))=0; climObs(:,:,mcur)=climObs(:,:,mcur)+tmp1;
103
104			%integrals :
105	heimbach	1.1	fld=fld_mod.mygrid.RAC.(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
106			fld=fld_dat.mygrid.RAC.(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);
107	gforget	1.5
108			fld=fld_mod.mygrid.RAC.(mygrid.YC<0); IceAreaSouthMod(1,mm)=nansum(fld);
109			fld=fld_dat.mygrid.RAC.(mygrid.YC<0); IceAreaSouthObs(1,mm)=nansum(fld);
110	heimbach	1.3
111			for kk=1:5;
112			tmpmsk=0.*mygrid.XC;
113	gforget	1.5	if lonLims(kk+1) > lonLims(kk)
114			tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC < lonLims(kk+1)))=1.;
115	heimbach	1.3	else
116	gforget	1.5	tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC <= 180.))=1.;
117			tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < lonLims(kk+1)))=1.;
118	heimbach	1.3	end
119			tmpmsk=tmpmsk.*(mygrid.YC<0);
120			%
121			fld=fld_mod.mygrid.RAC.tmpmsk;
122	gforget	1.5	IceAreaSouthMod(kk+1,mm)=nansum(fld);
123	heimbach	1.3	%
124			fld=fld_dat.mygrid.RAC.tmpmsk;
125	gforget	1.5	IceAreaSouthObs(kk+1,mm)=nansum(fld);
126	heimbach	1.3	end
127	heimbach	1.1
128			end;
129			toc;
130			end;
131
132	gforget	1.7	%misfits :
133	gforget	1.8	mis_rms=sqrt(nanmean(fld_dif.^2,3));
134			obs_std=nanstd(fld_nsidc,[],3);
135			mod_std=nanstd(fld_nsidc+fld_dif,[],3);
136	heimbach	1.1
137	gforget	1.7	%climatology :
138			for mcur=1:12;
139			tmp1=climMsk(:,:,mcur); tmp1(tmp1==0)=NaN;
140			climNb(:,:,mcur)=tmp1;
141			climMod(:,:,mcur)=climMod(:,:,mcur)./tmp1;
142			climObs(:,:,mcur)=climObs(:,:,mcur)./tmp1;
143			end;
144			clear climMsk;
145
146	gforget	1.10	if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
147			eval(['save ' dirMat '/cost/cost_seaicearea.mat fld_err mis_rms obs_std mod_std IceArea* clim*;']);
148	heimbach	1.1
149			else;%display previously computed results
150
151	gforget	1.9	if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
152
153	gforget	1.8	eval(['load ' dirMat '/cost_seaicearea.mat;']);
154
155			%variance maps:
156			figure; m_map_gcmfaces(mis_rms,0,{'myCaxis',[0:0.1:1.]});
157			myCaption={'modeled-observed rms -- sea ice concentration'};
158			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
159
160			figure; m_map_gcmfaces(obs_std,0,{'myCaxis',[0:0.1:1.]});
161			myCaption={'observed std -- sea ice concentration'};
162			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
163
164			figure; m_map_gcmfaces(mod_std,0,{'myCaxis',[0:0.1:1.]});
165			myCaption={'modelled std -- sea ice concentration'};
166			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
167	heimbach	1.1
168	gforget	1.8	%arctic maps
169			figureL;
170			subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),2,{'myCaxis',[0:0.1:1.]});
171			subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),2,{'myCaxis',[0:0.1:1.]});
172			subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),2,{'myCaxis',[0:0.1:1.]});
173			subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),2,{'myCaxis',[0:0.1:1.]});
174	gforget	1.12	myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'};
175	gforget	1.8	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
176	heimbach	1.1
177	gforget	1.8	%southern ocean maps
178			figureL;
179			subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),3,{'myCaxis',[0:0.1:1.]});
180			subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),3,{'myCaxis',[0:0.1:1.]});
181			subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),3,{'myCaxis',[0:0.1:1.]});
182			subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),3,{'myCaxis',[0:0.1:1.]});
183	gforget	1.12	myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'};
184	gforget	1.8	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
185	heimbach	1.1
186	gforget	1.8	%northern/southern integrals
187			ny=length(IceAreaNorthMod)/12;
188			yy=yearsta+[0:ny-1];
189
190			figureL;
191			subplot(1,2,1);
192			plot(yy,IceAreaNorthMod(3:12:end),'LineWidth',2); hold on;
193			plot(yy,IceAreaNorthObs(3:12:end),'r','LineWidth',2);
194			plot(yy,IceAreaNorthMod(9:12:end),'LineWidth',2);
195			plot(yy,IceAreaNorthObs(9:12:end),'r','LineWidth',2);
196			axis([yearsta yearsta+ny-1 0 20e12]);
197			ylabel('m^2'); title('Northern Hemisphere');
198			subplot(1,2,2);
199			plot(yy,IceAreaSouthMod(1,3:12:end),'LineWidth',2); hold on;
200			plot(yy,IceAreaSouthObs(1,3:12:end),'r','LineWidth',2);
201			plot(yy,IceAreaSouthMod(1,9:12:end),'LineWidth',2);
202			plot(yy,IceAreaSouthObs(1,9:12:end),'r','LineWidth',2);
203			axis([yearsta yearsta+ny-1 0 20e12]);
204			ylabel('m^2'); title('Southern Hemisphere');
205	heimbach	1.1
206	gforget	1.12	myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in March and September',...
207	gforget	1.8	'in Northern Hemisphere (left) and Southern Hemisphere (right)'};
208	heimbach	1.1	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
209
210	gforget	1.8	%southern basin integrals
211
212			for mm=[3 9];
213			figureL;
214			for ii=1:6;
215			subplot(3,2,ii);
216			if ii==1; iiTxt='Entire Southern Ocean';
217			else; iiTxt=sprintf('%dE to %dE',lonLims(ii-1),lonLims(ii));
218			end;
219			plot(yy,IceAreaSouthMod(ii,mm:12:end),'LineWidth',2); hold on;
220			plot(yy,IceAreaSouthObs(ii,mm:12:end),'r','LineWidth',2);
221			aa=axis; aa(1:2)=[yearsta yearsta+ny-1]; axis(aa);
222			ylabel('m^2'); title(iiTxt);
223			end;
224			if mm==3; mmTxt='March'; elseif mm==9; mmTxt='September'; else; '???'; end;
225	gforget	1.12	myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in ',mmTxt,' per Southern Ocean sector'};
226	gforget	1.8	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
227			end;
228
229	heimbach	1.1	end;
230
231