/[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.6 by gforget, Tue Apr 2 13:55:53 2013 UTC revision 1.13 by gforget, Tue May 10 21:38:51 2016 UTC
# Line 15  function []=insitu_cost(doComp,varargin) Line 15  function []=insitu_cost(doComp,varargin)
15  %  %
16  %example:         insitu_cost('./',1,{'listData','argo_in*'});  %example:         insitu_cost('./',1,{'listData','argo_in*'});
17    
18  %load paths and grid:  gcmfaces_global;
19  %--------------------  global myparms;
20  global mygrid; if isempty(mygrid); grid_load('GRID/',5,'compact'); end;  
   
 %set more optional paramaters to default values  
 global myparms; year0=myparms.yearInAve(1); year1=myparms.yearInAve(2);  
 listData={'argo_*_1992_to_2007*','argo_*_2008_to_2010*',...  
           'argo_*_2011_to_2012*','climode_*',...  
           'ctd_*','itp_*','seals_*','xbt_*'};  
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  dirData='/profiles/output/';  dirData='/profiles/output/';
27  dirMat='mat/'; suffMat='_all';  dirMat='mat/'; suffMat='_all';
# Line 46  for ii=1:nargin-1; Line 43  for ii=1:nargin-1;
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  date0=datenum(year0,1,1); date1=datenum( Line 64  date0=datenum(year0,1,1); date1=datenum(
64    
65  for vv=1:length(listVar);  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],listData,varCur);      [MITprof]=MITprof_stats_load([dirData],listData,varCur);
70            if ~isfield(MITprof,'prof_point');
71          loc_tile=gcmfaces_loc_tile(90,90,MITprof.prof_lon,MITprof.prof_lat);
72          MITprof.prof_point=loc_tile.point;
73        end;
74      %mask out values that are not in year range:      %mask out values that are not in year range:
75      ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);      ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);
76      MITprof.prof(ii,:)=NaN;      MITprof.prof(ii,:)=NaN;
# Line 82  for vv=1:length(listVar); Line 96  for vv=1:length(listVar);
96      eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);      eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);
97      eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);      eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);
98    
99        %upper ocean temporal distribution
100        for bb=1:length(listBas);
101        for ll=1:3;
102          bbb=listBas{bb};
103          if ll==1;     Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
104          elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
105          elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
106          end;
107          msk=v4_basin(bbb);
108          msk(mygrid.YC<Lmin|mygrid.YC>Lmax)=0;
109          %
110          msk=convert2array(msk);
111          MITprof.prof_msk=msk(MITprof.prof_point);    
112          MITprofSub=MITprof_subset(MITprof,'msk',1,'depth',[0 700]);      
113          %
114          tmp_date=(MITprofSub.prof_date-date0)/365;
115          [x,y,z,n]=MITprof_stats(tmp_date,[0:1/4:1+(year1-year0)],MITprofSub.prof,[-5:0.25:5]);
116          z(n==0)=NaN; n(n==0)=NaN;
117          eval([txt '.x=x;']); eval([txt '.y=y;']);
118          eval([txt '.z=z;']); eval([txt '.n=n;']);
119        end;
120        end;
121    
122  end;  end;
123    
124  eval(['save ' dirMat 'insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* listData listVar;']);  if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
125    eval(['save ' dirMat 'cost/insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* ' ...
126          'listData listVar listBas ' listBasTxt ';']);
127    
128  else;%display result  else;%display result
129    
130    if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
131    
132  eval(['load ' dirMat 'insitu_cost' suffMat '.mat;']);  eval(['load ' dirMat 'insitu_cost' suffMat '.mat;']);
133    
134  figureL;  figureL;
# Line 99  for vv=1:length(listVar); Line 140  for vv=1:length(listVar);
140    
141      eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);      eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);
142      eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);      eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);
143      z(n<1e2)=NaN;      z(n<1e2)=NaN; y=y+year0-1;
144      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]);
145      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) ')']);  
146      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));
147    
148      %latitudinal distribution:      %latitudinal distribution:
# Line 118  myCaption={'Cost function (top) for in s Line 158  myCaption={'Cost function (top) for in s
158          '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).'};
159  if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;  if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
160    
161    if ~isempty(listBas);
162    
163    ii=0;
164    ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=3;
165    ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=3;
166    ii=ii+1; listPanels(ii).bas='arct'; listPanels(ii).ll=3;
167    ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=2;
168    ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=2;
169    ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=2;
170    ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=1;
171    ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=1;
172    ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=1;
173    
174    for vv=1:length(listVar);
175      varCur=listVar{vv};
176      figureL;
177      for ii=1:9;
178          bbb=listPanels(ii).bas; ll=listPanels(ii).ll;
179          if ll==1;     Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
180          elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
181          elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
182          end;
183    
184        eval(['x=' txt '.x;']); eval(['y=' txt '.y;']);
185        eval(['z=' txt '.z;']); eval(['n=' txt '.n;']);
186        x=x+year0-1;
187        tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
188        subplot(3,3,ii); contourf(x,y,z,[0:0.05:0.5]); grid on;
189        caxis([0 0.5]); colorbar; title(txt,'Interpreter','none');
190        %title([txt ', normalized misfit pdf'],'Interpreter','none');
191        if mod(ii,3)==1; ylabel('normalized misfit'); end;
192        %if ii<=6; set(gca,'XTick',[]); end;
193      end;
194    
195    myCaption={'Distribution of normalized misfits per basin (panel) as ',...
196            'a function of latitude, for ',varCur};
197    if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
198    
199    end;
200    
201    end;
202    
203  end;  end;
204    
205    
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22