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;

%load grid
gcmfaces_global;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;

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

%%%dirData='/net/weddell/raid3/gforget/ecco_v4/input_files/';
dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';

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

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

fld_dif=convert2gcmfaces(NaN*ones(90,90*13,nmonths));

IceAreaNorthMod=zeros(1,nmonths);
IceAreaNorthObs=zeros(1,nmonths);
IceAreaNorthModMin=zeros(1,nyears);
IceAreaNorthObsMin=zeros(1,nyears);

IceAreaSouthMod=zeros(6,nmonths);
IceAreaSouthObs=zeros(6,nmonths);
IceAreaSouthModMax=zeros(6,nyears);
IceAreaSouthObsMax=zeros(6,nyears);

[core2lon,core2lat,core2names] = line_greatC_TUV_MASKS_core2_antarctic;
core2lon(6)=core2lon(1); 
core2lat(6)=core2lat(1);

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

%%%    fld_dat=v4_read_bin([dirData 'gael_quick_interp/g_SST_monthly_r2_' num2str(ycur)],mcur,0);
%%%    fld_dat=v4_read_bin([dirData 'input_nsidc_postproc/NSIDC0079_v4_gael_' num2str(ycur)],mcur,0);
    fld_dat=v4_read_bin([dirData 'input_nsidc_all/nsidc79_monthly_' num2str(ycur)],mcur,0);

  fld_dat(find(fld_dat==0))=NaN;
  fld_dat(find(fld_dat<-99))=NaN;

  fld_mod=v4_read_bin([dirModel fileModel],mm,1).*mygrid.mskC(:,:,1);
 
  fld_dif(:,:,mm)=fld_mod-fld_dat;

  fld=fld_mod.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
  fld=fld_dat.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);
  if mcur==9;
   IceAreaNorthModMin(ycur-yearsta+1) = IceAreaNorthMod(mm);
   IceAreaNorthObsMin(ycur-yearsta+1) = IceAreaNorthObs(mm);
  end

    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);
    if mcur==9;
       IceAreaSouthModMax(1,ycur-yearsta+1) = IceAreaSouthMod(1,mm);
       IceAreaSouthObsMax(1,ycur-yearsta+1) = IceAreaSouthObs(1,mm);
    end

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

    end

 end;
 toc;
end;

fld_rms=sqrt(nanmean(fld_dif.^2,3));
fld_cost=fld_rms.^2.*fld_w;

eval(['save ' dirMat '/cost_seaicearea.mat fld_err fld_rms fld_cost IceAreaNorthMod IceAreaSouthMod IceAreaNorthObs IceAreaSouthObs IceAreaNorthModMin IceAreaNorthObsMin IceAreaSouthModMax IceAreaSouthObsMax;']);

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			%load grid
28			gcmfaces_global;
29			if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
30			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
31
32			fld_err=ones(90,1170);
33			fld_err=convert2gcmfaces(fld_err);
34			fld_w=fld_err.^-2;
35
36			%%%dirData='/net/weddell/raid3/gforget/ecco_v4/input_files/';
37	heimbach	1.3	dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';
38	heimbach	1.1
39			fileModel=dir([dirModel 'barfiles/gbar_area*data']);
40			fileModel=['barfiles/' fileModel.name];
41
42	heimbach	1.3	yearsta=1992;
43			yearend=2011;
44			nyears=yearend-yearsta+1;
45			nmonths=12*nyears;
46
47			fld_dif=convert2gcmfaces(NaNones(90,9013,nmonths));
48
49			IceAreaNorthMod=zeros(1,nmonths);
50			IceAreaNorthObs=zeros(1,nmonths);
51			IceAreaNorthModMin=zeros(1,nyears);
52			IceAreaNorthObsMin=zeros(1,nyears);
53
54			IceAreaSouthMod=zeros(6,nmonths);
55			IceAreaSouthObs=zeros(6,nmonths);
56			IceAreaSouthModMax=zeros(6,nyears);
57			IceAreaSouthObsMax=zeros(6,nyears);
58
59			[core2lon,core2lat,core2names] = line_greatC_TUV_MASKS_core2_antarctic;
60			core2lon(6)=core2lon(1);
61			core2lat(6)=core2lat(1);
62	heimbach	1.1
63	heimbach	1.3	for ycur=yearsta:yearend;
64	heimbach	1.1	tic;
65			for mcur=1:12;
66	heimbach	1.3	mm=(ycur-yearsta)*12+mcur;
67
68	heimbach	1.1	%%% fld_dat=v4_read_bin([dirData 'gael_quick_interp/g_SST_monthly_r2_' num2str(ycur)],mcur,0);
69			%%% fld_dat=v4_read_bin([dirData 'input_nsidc_postproc/NSIDC0079_v4_gael_' num2str(ycur)],mcur,0);
70			fld_dat=v4_read_bin([dirData 'input_nsidc_all/nsidc79_monthly_' num2str(ycur)],mcur,0);
71
72			fld_dat(find(fld_dat==0))=NaN;
73			fld_dat(find(fld_dat<-99))=NaN;
74
75			fld_mod=v4_read_bin([dirModel fileModel],mm,1).*mygrid.mskC(:,:,1);
76
77			fld_dif(:,:,mm)=fld_mod-fld_dat;
78
79			fld=fld_mod.mygrid.RAC.(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
80			fld=fld_dat.mygrid.RAC.(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);
81	heimbach	1.3	if mcur==9;
82			IceAreaNorthModMin(ycur-yearsta+1) = IceAreaNorthMod(mm);
83			IceAreaNorthObsMin(ycur-yearsta+1) = IceAreaNorthObs(mm);
84			end
85
86			fld=fld_mod.mygrid.RAC.(mygrid.YC<0); IceAreaSouthMod(1,mm)=nansum(fld);
87			fld=fld_dat.mygrid.RAC.(mygrid.YC<0); IceAreaSouthObs(1,mm)=nansum(fld);
88			if mcur==9;
89			IceAreaSouthModMax(1,ycur-yearsta+1) = IceAreaSouthMod(1,mm);
90			IceAreaSouthObsMax(1,ycur-yearsta+1) = IceAreaSouthObs(1,mm);
91			end
92
93			for kk=1:5;
94			iii=kk+1;
95			tmpmsk=0.*mygrid.XC;
96			if core2lon(kk+1) > core2lon(kk)
97			tmpmsk(find(mygrid.XC >= core2lon(kk) & mygrid.XC < core2lon(kk+1)))=1.;
98			else
99			tmpmsk(find(mygrid.XC >= core2lon(kk) & mygrid.XC <= 180.))=1.;
100			tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < core2lon(kk+1)))=1.;
101			end
102			tmpmsk=tmpmsk.*(mygrid.YC<0);
103			%
104			fld=fld_mod.mygrid.RAC.tmpmsk;
105			IceAreaSouthMod(iii,mm)=nansum(fld);
106			%
107			fld=fld_dat.mygrid.RAC.tmpmsk;
108			IceAreaSouthObs(iii,mm)=nansum(fld);
109			%
110			if mcur==9;
111			IceAreaSouthModMax(iii,ycur-yearsta+1) = IceAreaSouthMod(iii,mm);
112			IceAreaSouthObsMax(iii,ycur-yearsta+1) = IceAreaSouthObs(iii,mm);
113			end
114
115			end
116	heimbach	1.1
117			end;
118			toc;
119			end;
120
121			fld_rms=sqrt(nanmean(fld_dif.^2,3));
122			fld_cost=fld_rms.^2.*fld_w;
123
124	heimbach	1.3	eval(['save ' dirMat '/cost_seaicearea.mat fld_err fld_rms fld_cost IceAreaNorthMod IceAreaSouthMod IceAreaNorthObs IceAreaSouthObs IceAreaNorthModMin IceAreaNorthObsMin IceAreaSouthModMax IceAreaSouthObsMax;']);
125	heimbach	1.1
126			else;%display previously computed results
127
128			global mygrid;
129
130			eval(['load ' dirMat '/cost_seaicearea.mat;']);
131
132			%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
133			%figure; m_map_gcmfaces(fld_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
134
135			figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.1:1.]});
136			myCaption={'modeled-observed rms -- sea ice area (K)'};
137
138			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
139
140			end;
141
142