/[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.6 - (hide annotations) (download)
Tue Apr 2 20:51:52 2013 UTC (12 years, 3 months ago) by gforget
Branch: MAIN
Changes since 1.5: +5 -3 lines
- mask obs and model consistently.

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 gforget 1.6 fld_dat=fld_dat.*mygrid.mskC(:,:,1);%land mask
70     fld_dat(find(fld_dat<-99))=NaN;%missing data
71     msk=1+0*fld_dat;%combined mask
72 heimbach 1.1
73 gforget 1.6 fld_mod=v4_read_bin([dirModel fileModel],mm,0);
74     fld_mod=fld_mod.*msk;%mask consistent with fld_dat
75 heimbach 1.1
76     fld_dif(:,:,mm)=fld_mod-fld_dat;
77    
78     fld=fld_mod.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
79     fld=fld_dat.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);
80 gforget 1.5
81     fld=fld_mod.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthMod(1,mm)=nansum(fld);
82     fld=fld_dat.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthObs(1,mm)=nansum(fld);
83 heimbach 1.3
84     for kk=1:5;
85     tmpmsk=0.*mygrid.XC;
86 gforget 1.5 if lonLims(kk+1) > lonLims(kk)
87     tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC < lonLims(kk+1)))=1.;
88 heimbach 1.3 else
89 gforget 1.5 tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC <= 180.))=1.;
90     tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < lonLims(kk+1)))=1.;
91 heimbach 1.3 end
92     tmpmsk=tmpmsk.*(mygrid.YC<0);
93     %
94     fld=fld_mod.*mygrid.RAC.*tmpmsk;
95 gforget 1.5 IceAreaSouthMod(kk+1,mm)=nansum(fld);
96 heimbach 1.3 %
97     fld=fld_dat.*mygrid.RAC.*tmpmsk;
98 gforget 1.5 IceAreaSouthObs(kk+1,mm)=nansum(fld);
99 heimbach 1.3 end
100 heimbach 1.1
101     end;
102     toc;
103     end;
104    
105     fld_rms=sqrt(nanmean(fld_dif.^2,3));
106     fld_cost=fld_rms.^2.*fld_w;
107    
108 gforget 1.5 eval(['save ' dirMat '/cost_seaicearea.mat fld_err fld_rms fld_cost IceArea*;']);
109 heimbach 1.1
110     else;%display previously computed results
111    
112     global mygrid;
113    
114     eval(['load ' dirMat '/cost_seaicearea.mat;']);
115    
116     %figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
117     %figure; m_map_gcmfaces(fld_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
118    
119     figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.1:1.]});
120     myCaption={'modeled-observed rms -- sea ice area (K)'};
121    
122     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
123    
124     end;
125    
126    

  ViewVC Help
Powered by ViewVC 1.1.22