/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/ecco_v4/insitu_cost.m
ViewVC logotype

Diff of /MITgcm_contrib/gael/profilesMatlabProcessing/ecco_v4/insitu_cost.m

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

revision 1.3 by gforget, Fri Jul 8 12:48:59 2011 UTC revision 1.12 by gforget, Mon Feb 15 18:07:31 2016 UTC
# Line 1  Line 1 
1  function []=insitu_cost(dirData,doComp,varargin);  function []=insitu_cost(doComp,varargin);
2  %[]=insitu_cost(dirData,listData,listVar,varargin);  %object:        computes or displays cost function statistics
3  %object:        computes or displays cost function statistics  %input:         doComp states whether to compute & save (1) or load & display (0)
 %input:         dirData is the data main directory name  
 %               doComp states whether to compute & save (1) or load & display (0)  
4  %optional must take the following form {'name',param1,param2,...} Active ones are  %optional must take the following form {'name',param1,param2,...} Active ones are
5  %               {'subdirData',subdirData} is the data subdirectory name ('profiles/output/' by default)  %               {'dirData',dirData} is the data subdirectory name ('profiles/output/' by default)
6  %               {'listData',data1,data2,...} where data1 e.g. is 'argo_indian.nc'  %               {'listData',data1,data2,...} where data1 e.g. is 'argo_indian.nc'
7  %               {'listVar',var1,var2} where var1 and var2 are 'T' and 'S'  %               {'listVar',var1,var2} where var1 and var2 are 'T' and 'S'
8  %               {'years',year0,year1} to limit computation to [year0 year1] interval  %               {'years',year0,year1} to limit computation to [year0 year1] interval
9  %               {'subdirMat',subdirMat} where subdirMat is the output subdirectory ('mat/' by default)  %               {'dirMat',dirMat} where dirMat is the output subdirectory ('mat/' by default)
10  %               {'suffMat',suffMat} is the mat file suffix ('all' by default)  %               {'suffMat',suffMat} is the mat file suffix ('all' by default)
11  %               {'subdirTex',subdirTex} where subdirTex is the tex directory ('tex/' by default)  %               {'dirTex',dirTex} where dirTex is the tex directory ('tex/' by default)
12    %               {'nameTex',nameTex} is the tex file name ('myPlots' by default)
13  %               {'addToTex',addToTex} states whether (1) or not (0 default) to  %               {'addToTex',addToTex} states whether (1) or not (0 default) to
14  %                       augment an ongoing tex file (see write2tex.m, subdirTex)  %                        augment an ongoing tex file (see write2tex.m, dirTex)
15  %  %
16  %example:        insitu_cost('./',1,{'listData','argo_in*'});  %example:         insitu_cost('./',1,{'listData','argo_in*'});
17    
18    gcmfaces_global;
19    global myparms;
20    
 %load paths and grid:  
 %--------------------  
 global mygrid; if isempty(mygrid); grid_load('GRID/',5,'compact'); end;  
   
 %set more optional paramaters to default values  
 year0=1992; year1=2008;  
 listData={'WOD09_XBT*','WOD09_CTD*','argo_in*','argo_at*','argo_pa*','seals*'};  
21  listVar={'T','S'};  listVar={'T','S'};
22    listBas={'atlExt','pacExt','indExt','arct'};
23    if sum([90 1170]~=mygrid.ioSize)>0; listBas={}; end;
24    listBasTxt=''; for bb=1:length(listBas); listBasTxt=[listBasTxt ' ' listBas{bb} '*']; end;
25  %shorter test case: listData={'argo_in*'}; listVar={'T'};  %shorter test case: listData={'argo_in*'}; listVar={'T'};
26  subdirData='/profiles/output/';  dirData='/profiles/output/';
27  subdirMat='mat/'; suffMat='_all';  dirMat='mat/'; suffMat='_all';
28  subdirTex='tex/'; addToTex=0;  dirTex='tex/'; addToTex=0; nameTex='myPlots';
29  %set more optional paramaters to user defined values  %set more optional paramaters to user defined values
30  for ii=1:nargin-2;  for ii=1:nargin-1;
31     if strcmp(varargin{ii}{1},'listData'); listData={varargin{ii}{2:end}};     if strcmp(varargin{ii}{1},'listData'); listData={varargin{ii}{2:end}};
32     elseif strcmp(varargin{ii}{1},'listVar'); listVar={varargin{ii}{2:end}};     elseif strcmp(varargin{ii}{1},'listVar'); listVar={varargin{ii}{2:end}};
33     elseif strcmp(varargin{ii}{1},'years'); year0=varargin{ii}{2}; year1=varargin{ii}{3};     elseif strcmp(varargin{ii}{1},'years'); year0=varargin{ii}{2}; year1=varargin{ii}{3};
34     elseif strcmp(varargin{ii}{1},'subdirMat'); subdirMat=varargin{ii}{2};     elseif strcmp(varargin{ii}{1},'dirMat'); dirMat=varargin{ii}{2};
35     elseif strcmp(varargin{ii}{1},'suffMat'); suffMat=['_' varargin{ii}{2}];     elseif strcmp(varargin{ii}{1},'suffMat'); suffMat=['_' varargin{ii}{2}];
36     elseif strcmp(varargin{ii}{1},'subdirTex'); subdirTex=varargin{ii}{2};     elseif strcmp(varargin{ii}{1},'dirTex'); dirTex=varargin{ii}{2};
37     elseif strcmp(varargin{ii}{1},'addToTex'); addToTex=varargin{ii}{2};     elseif strcmp(varargin{ii}{1},'addToTex'); addToTex=varargin{ii}{2};
38     elseif strcmp(varargin{ii}{1},'subdirData'); subdirData=varargin{ii}{2};     elseif strcmp(varargin{ii}{1},'nameTex'); nameTex=varargin{ii}{2};
39       elseif strcmp(varargin{ii}{1},'dirData'); dirData=varargin{ii}{2};
40     else;     else;
41         warning('inputCheck:insitu_cost',...         warning('inputCheck:insitu_cost',...
42             ['unknown option ''' varargin{ii}{1} ''' was ignored']);             ['unknown option ''' varargin{ii}{1} ''' was ignored']);
43     end;     end;
44  end;  end;
45    
46    if isempty(myparms)|isempty(mygrid);
47      load([dirMat 'diags_grid_parms.mat']);
48    end;
49    
50    if isempty(whos('year0'));
51      year0=myparms.yearInAve(1); year1=myparms.yearInAve(2);
52    end;
53    
54    if isempty(whos('listData'));
55      listData=dir([dirMat 'profiles/output/*.nc'])
56      listData={listData(:).name};
57      for ff=1:length(listData); listData{ff}=[listData{ff}(1:end-3) '*']; end;
58    end;
59    
60  if doComp;  if doComp;
61    
62  %time limits:  %time limits:
# Line 53  for vv=1:length(listVar); Line 66  for vv=1:length(listVar);
66      varCur=listVar{vv};      varCur=listVar{vv};
67            
68      %get the cost square root:      %get the cost square root:
69      [MITprof]=MITprof_stats_load([dirData '/' subdirData],listData,varCur);      [MITprof]=MITprof_stats_load([dirData],listData,varCur);
70            
71      %mask out values that are not in year range:      %mask out values that are not in year range:
72      ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);      ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);
# Line 80  for vv=1:length(listVar); Line 93  for vv=1:length(listVar);
93      eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);      eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);
94      eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);      eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);
95    
96        %upper ocean temporal distribution
97        for bb=1:length(listBas);
98        for ll=1:3;
99          bbb=listBas{bb};
100          if ll==1;     Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
101          elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
102          elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
103          end;
104          msk=v4_basin(bbb);
105          msk(mygrid.YC<Lmin|mygrid.YC>Lmax)=0;
106          %
107          msk=convert2array(msk);
108          MITprof.prof_msk=msk(MITprof.prof_point);    
109          MITprofSub=MITprof_subset(MITprof,'msk',1,'depth',[0 700]);      
110          %
111          tmp_date=(MITprofSub.prof_date-date0)/365;
112          [x,y,z,n]=MITprof_stats(tmp_date,[0:1/4:1+(year1-year0)],MITprofSub.prof,[-5:0.25:5]);
113          z(n==0)=NaN; n(n==0)=NaN;
114          eval([txt '.x=x;']); eval([txt '.y=y;']);
115          eval([txt '.z=z;']); eval([txt '.n=n;']);
116        end;
117        end;
118    
119  end;  end;
120    
121  eval(['save ' dirData '/' subdirMat 'insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* listData listVar;']);  if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
122    eval(['save ' dirMat 'cost/insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* ' ...
123          'listData listVar listBas ' listBasTxt ';']);
124    
125  else;%display result  else;%display result
126    
127  eval(['load ' dirData '/' subdirMat 'insitu_cost' suffMat '.mat;']);  if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
128    
129    eval(['load ' dirMat 'insitu_cost' suffMat '.mat;']);
130    
131  figureL;  figureL;
132    
# Line 97  for vv=1:length(listVar); Line 137  for vv=1:length(listVar);
137    
138      eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);      eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);
139      eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);      eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);
140      z(n<1e2)=NaN;      z(n<1e2)=NaN; y=y+year0-1;
141      subplot(2,length(listVar),vv); depthStretchPlot('pcolor',{y,x,z},[0:200:1000 1500:500:6000],[0 1000 4000]);      subplot(2,length(listVar),vv); depthStretchPlot('pcolor',{y,x,z},[0:200:1000 1500:500:6000],[0 1000 4000]);
142      caxis([0 8]); shading flat; colorbar; ylabel('depth (in m)');      caxis([0 8]); shading flat; colorbar; ylabel('depth (in m)');
     xlabel(['date (in years since ' num2str(year0) ')']);  
143      title(sprintf('mean cost for %s: %0.3g (median: %0.3g)',varCur,costCur));      title(sprintf('mean cost for %s: %0.3g (median: %0.3g)',varCur,costCur));
144    
145      %latitudinal distribution:      %latitudinal distribution:
# Line 113  for vv=1:length(listVar); Line 152  for vv=1:length(listVar);
152  end;  end;
153    
154  myCaption={'Cost function (top) for in situ profiles, as a function of depth and time. ',...  myCaption={'Cost function (top) for in situ profiles, as a function of depth and time. ',...
155          'Distribution of normalized misfits (bottom) as a function of latitude. For T (left) and S (right).'};          'Distribution of normalized misfits (bottom) as a function of latitude. For T (left) and S (right).'};
156  if addToTex; write2tex([dirData '/' subdirTex '/myPlots.tex'],2,myCaption,gcf); end;  if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
157    
158    if ~isempty(listBas);
159    
160    ii=0;
161    ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=3;
162    ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=3;
163    ii=ii+1; listPanels(ii).bas='arct'; listPanels(ii).ll=3;
164    ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=2;
165    ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=2;
166    ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=2;
167    ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=1;
168    ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=1;
169    ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=1;
170    
171    for vv=1:length(listVar);
172      varCur=listVar{vv};
173      figureL;
174      for ii=1:9;
175          bbb=listPanels(ii).bas; ll=listPanels(ii).ll;
176          if ll==1;     Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
177          elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
178          elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
179          end;
180    
181        eval(['x=' txt '.x;']); eval(['y=' txt '.y;']);
182        eval(['z=' txt '.z;']); eval(['n=' txt '.n;']);
183        x=x+year0-1;
184        tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
185        subplot(3,3,ii); contourf(x,y,z,[0:0.05:0.5]); grid on;
186        caxis([0 0.5]); colorbar; title(txt,'Interpreter','none');
187        %title([txt ', normalized misfit pdf'],'Interpreter','none');
188        if mod(ii,3)==1; ylabel('normalized misfit'); end;
189        %if ii<=6; set(gca,'XTick',[]); end;
190      end;
191    
192    myCaption={'Distribution of normalized misfits per basin (panel) as ',...
193            'a function of latitude, for ',varCur};
194    if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
195    
196    end;
197    
198    end;
199    
200  end;  end;
201    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22