/[MITgcm]/MITgcm_contrib/gael/matlab_class/ecco_v4/cost_seaicearea.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/ecco_v4/cost_seaicearea.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.5 - (hide annotations) (download)
Tue Apr 2 20:37:53 2013 UTC (12 years, 3 months ago) by gforget
Branch: MAIN
Changes since 1.4: +18 -36 lines
- cosmetic changes.

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     %%%dirData='/net/weddell/raid3/gforget/ecco_v4/input_files/';
41 heimbach 1.3 dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';
42 heimbach 1.1
43     fileModel=dir([dirModel 'barfiles/gbar_area*data']);
44     fileModel=['barfiles/' fileModel.name];
45    
46 heimbach 1.3 nyears=yearend-yearsta+1;
47     nmonths=12*nyears;
48    
49     fld_dif=convert2gcmfaces(NaN*ones(90,90*13,nmonths));
50    
51 gforget 1.4 IceAreaNorthMod=NaN*zeros(1,nmonths);
52     IceAreaNorthObs=NaN*zeros(1,nmonths);
53    
54     IceAreaSouthMod=NaN*zeros(6,nmonths);
55     IceAreaSouthObs=NaN*zeros(6,nmonths);
56 heimbach 1.3
57 gforget 1.5 [lonPairs,latPairs,names] = line_greatC_TUV_MASKS_core2_antarctic;
58     lonLims=[lonPairs(1:5,1);lonPairs(1,1)];
59 heimbach 1.1
60 heimbach 1.3 for ycur=yearsta:yearend;
61 heimbach 1.1 tic;
62     for mcur=1:12;
63 heimbach 1.3 mm=(ycur-yearsta)*12+mcur;
64    
65 heimbach 1.1 %%% fld_dat=v4_read_bin([dirData 'gael_quick_interp/g_SST_monthly_r2_' num2str(ycur)],mcur,0);
66     %%% fld_dat=v4_read_bin([dirData 'input_nsidc_postproc/NSIDC0079_v4_gael_' num2str(ycur)],mcur,0);
67     fld_dat=v4_read_bin([dirData 'input_nsidc_all/nsidc79_monthly_' num2str(ycur)],mcur,0);
68    
69     fld_dat(find(fld_dat==0))=NaN;
70     fld_dat(find(fld_dat<-99))=NaN;
71    
72 gforget 1.4 fld_mod=v4_read_bin([dirModel fileModel],mm,0).*mygrid.mskC(:,:,1);
73 heimbach 1.1
74     fld_dif(:,:,mm)=fld_mod-fld_dat;
75    
76     fld=fld_mod.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
77     fld=fld_dat.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);
78 gforget 1.5
79     fld=fld_mod.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthMod(1,mm)=nansum(fld);
80     fld=fld_dat.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthObs(1,mm)=nansum(fld);
81 heimbach 1.3
82     for kk=1:5;
83     tmpmsk=0.*mygrid.XC;
84 gforget 1.5 if lonLims(kk+1) > lonLims(kk)
85     tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC < lonLims(kk+1)))=1.;
86 heimbach 1.3 else
87 gforget 1.5 tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC <= 180.))=1.;
88     tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < lonLims(kk+1)))=1.;
89 heimbach 1.3 end
90     tmpmsk=tmpmsk.*(mygrid.YC<0);
91     %
92     fld=fld_mod.*mygrid.RAC.*tmpmsk;
93 gforget 1.5 IceAreaSouthMod(kk+1,mm)=nansum(fld);
94 heimbach 1.3 %
95     fld=fld_dat.*mygrid.RAC.*tmpmsk;
96 gforget 1.5 IceAreaSouthObs(kk+1,mm)=nansum(fld);
97 heimbach 1.3 end
98 heimbach 1.1
99     end;
100     toc;
101     end;
102    
103     fld_rms=sqrt(nanmean(fld_dif.^2,3));
104     fld_cost=fld_rms.^2.*fld_w;
105    
106 gforget 1.5 eval(['save ' dirMat '/cost_seaicearea.mat fld_err fld_rms fld_cost IceArea*;']);
107 heimbach 1.1
108     else;%display previously computed results
109    
110     global mygrid;
111    
112     eval(['load ' dirMat '/cost_seaicearea.mat;']);
113    
114     %figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
115     %figure; m_map_gcmfaces(fld_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
116    
117     figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.1:1.]});
118     myCaption={'modeled-observed rms -- sea ice area (K)'};
119    
120     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
121    
122     end;
123    
124    

  ViewVC Help
Powered by ViewVC 1.1.22