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/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat)); eval(['!mkdir ' dirMat ';']); end;

%determine if and where to create tex and figures files
dirMat=[dirMat '/'];
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 nameTex '.tex'];
end;

if doComp;

%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;

fld_err=ones(90,1170);
fld_err=convert2gcmfaces(fld_err);
fld_w=fld_err.^-2;

dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';

fileModel=dir([dirModel 'barfiles/gbar_area*data']); 
fileModel=['barfiles/' fileModel.name];

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

%misfits :
fld_dif=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);

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

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

  fld_dat=v4_read_bin([dirData 'input_nsidc_all/nsidc79_monthly_' 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=v4_read_bin([dirModel fileModel],mm,0);
  fld_mod=fld_mod.*msk;%mask consistent with fld_dat

  %misfits : 
  fld_dif(:,:,mm)=fld_mod-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 :
fld_rms=sqrt(nanmean(fld_dif.^2,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;

eval(['save ' dirMat '/cost_seaicearea.mat fld_err fld_rms IceArea* clim*;']);

else;%display previously computed results

global mygrid;

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

%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
%figure; m_map_gcmfaces(fld_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});

figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.1:1.]});
myCaption={'modeled-observed rms -- sea ice area (K)'};

if addToTex; write2tex(fileTex,2,myCaption,gcf); 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			if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
12			if isempty(dir(dirMat)); eval(['!mkdir ' dirMat ';']); end;
13
14			%determine if and where to create tex and figures files
15			dirMat=[dirMat '/'];
16			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			fileTex=[dirTex nameTex '.tex'];
23	heimbach	1.1	end;
24
25			if doComp;
26
27	gforget	1.5	%grid, params and inputs
28
29			gcmfaces_global; global myparms;
30	heimbach	1.1	if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
31			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
32	gforget	1.5	if isfield(myparms,'yearsta'); yearsta=myparms.yearsta; yearend=myparms.yearend;
33			else; yearsta=1992; yearend=2011;
34			end;
35	heimbach	1.1
36			fld_err=ones(90,1170);
37			fld_err=convert2gcmfaces(fld_err);
38			fld_w=fld_err.^-2;
39
40	heimbach	1.3	dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';
41	heimbach	1.1
42			fileModel=dir([dirModel 'barfiles/gbar_area*data']);
43			fileModel=['barfiles/' fileModel.name];
44
45	heimbach	1.3	nyears=yearend-yearsta+1;
46			nmonths=12*nyears;
47
48	gforget	1.7	%misfits :
49	heimbach	1.3	fld_dif=convert2gcmfaces(NaNones(90,9013,nmonths));
50
51	gforget	1.7	%monthly mean climatology :
52			climMod=convert2gcmfaces(zeros(90,90*13,12));
53			climObs=convert2gcmfaces(zeros(90,90*13,12));
54			climMsk=convert2gcmfaces(zeros(90,90*13,12));
55			climNb=convert2gcmfaces(zeros(90,90*13,12));
56
57			%monthly integrals :
58	gforget	1.4	IceAreaNorthMod=NaN*zeros(1,nmonths);
59			IceAreaNorthObs=NaN*zeros(1,nmonths);
60
61			IceAreaSouthMod=NaN*zeros(6,nmonths);
62			IceAreaSouthObs=NaN*zeros(6,nmonths);
63	heimbach	1.3
64	gforget	1.5	[lonPairs,latPairs,names] = line_greatC_TUV_MASKS_core2_antarctic;
65			lonLims=[lonPairs(1:5,1);lonPairs(1,1)];
66	heimbach	1.1
67	gforget	1.7	%computational loop :
68	heimbach	1.3	for ycur=yearsta:yearend;
69	heimbach	1.1	tic;
70			for mcur=1:12;
71	heimbach	1.3	mm=(ycur-yearsta)*12+mcur;
72
73	gforget	1.7	fld_dat=v4_read_bin([dirData 'input_nsidc_all/nsidc79_monthly_' num2str(ycur)],mcur,0);
74	heimbach	1.1
75	gforget	1.6	fld_dat=fld_dat.*mygrid.mskC(:,:,1);%land mask
76			fld_dat(find(fld_dat<-99))=NaN;%missing data
77			msk=1+0*fld_dat;%combined mask
78	heimbach	1.1
79	gforget	1.6	fld_mod=v4_read_bin([dirModel fileModel],mm,0);
80			fld_mod=fld_mod.*msk;%mask consistent with fld_dat
81	gforget	1.7
82			%misfits :
83	heimbach	1.1	fld_dif(:,:,mm)=fld_mod-fld_dat;
84
85	gforget	1.7	%climatology :
86			tmp1=msk; tmp1(isnan(tmp1))=0; climMsk(:,:,mcur)=climMsk(:,:,mcur)+tmp1;
87			tmp1=fld_mod; tmp1(isnan(tmp1))=0; climMod(:,:,mcur)=climMod(:,:,mcur)+tmp1;
88			tmp1=fld_dat; tmp1(isnan(tmp1))=0; climObs(:,:,mcur)=climObs(:,:,mcur)+tmp1;
89
90			%integrals :
91	heimbach	1.1	fld=fld_mod.mygrid.RAC.(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
92			fld=fld_dat.mygrid.RAC.(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);
93	gforget	1.5
94			fld=fld_mod.mygrid.RAC.(mygrid.YC<0); IceAreaSouthMod(1,mm)=nansum(fld);
95			fld=fld_dat.mygrid.RAC.(mygrid.YC<0); IceAreaSouthObs(1,mm)=nansum(fld);
96	heimbach	1.3
97			for kk=1:5;
98			tmpmsk=0.*mygrid.XC;
99	gforget	1.5	if lonLims(kk+1) > lonLims(kk)
100			tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC < lonLims(kk+1)))=1.;
101	heimbach	1.3	else
102	gforget	1.5	tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC <= 180.))=1.;
103			tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < lonLims(kk+1)))=1.;
104	heimbach	1.3	end
105			tmpmsk=tmpmsk.*(mygrid.YC<0);
106			%
107			fld=fld_mod.mygrid.RAC.tmpmsk;
108	gforget	1.5	IceAreaSouthMod(kk+1,mm)=nansum(fld);
109	heimbach	1.3	%
110			fld=fld_dat.mygrid.RAC.tmpmsk;
111	gforget	1.5	IceAreaSouthObs(kk+1,mm)=nansum(fld);
112	heimbach	1.3	end
113	heimbach	1.1
114			end;
115			toc;
116			end;
117
118	gforget	1.7	%misfits :
119	heimbach	1.1	fld_rms=sqrt(nanmean(fld_dif.^2,3));
120
121	gforget	1.7	%climatology :
122			for mcur=1:12;
123			tmp1=climMsk(:,:,mcur); tmp1(tmp1==0)=NaN;
124			climNb(:,:,mcur)=tmp1;
125			climMod(:,:,mcur)=climMod(:,:,mcur)./tmp1;
126			climObs(:,:,mcur)=climObs(:,:,mcur)./tmp1;
127			end;
128			clear climMsk;
129
130			eval(['save ' dirMat '/cost_seaicearea.mat fld_err fld_rms IceArea* clim*;']);
131	heimbach	1.1
132			else;%display previously computed results
133
134			global mygrid;
135
136			eval(['load ' dirMat '/cost_seaicearea.mat;']);
137
138			%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
139			%figure; m_map_gcmfaces(fld_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
140
141			figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.1:1.]});
142			myCaption={'modeled-observed rms -- sea ice area (K)'};
143
144			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
145
146			end;
147
148