/[MITgcm]/MITgcm_contrib/enderton/Diagnostics/DiagLoadGradsData.m
ViewVC logotype

Diff of /MITgcm_contrib/enderton/Diagnostics/DiagLoadGradsData.m

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

revision 1.1 by enderton, Mon Jan 31 15:43:26 2005 UTC revision 1.4 by molod, Tue Jun 28 21:33:51 2005 UTC
# Line 1  Line 1 
1  function [data,xax,yax,zax,months,time,dim] = ...  function [data,xax,yax,zax,months,time,dim] = ...
2      DiagLoadGradsData(Grads,dad,fln)      DiagLoadGradsData(Grads,dad,fln,DiagDebug)
3    
4  % Read table file.  
5  tablfile = [dad,'/',Grads,'.tabl'];  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6  datafile = [dad,'/',Grads,'.data'];  %                          Parse table file                               %
7    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8    ShiftData = 0;
9    format = 'NONSEQUENTIAL';
10    
11    tablfile = [dad,'/',Grads];
12  file = textread(tablfile,'%s','delimiter','\n','whitespace','');  file = textread(tablfile,'%s','delimiter','\n','whitespace','');
13    
14  for iline = 1:length(file)  for iline = 1:length(file)
15      rem = file{iline}; tokens = {}; itoken = 0;      rem = file{iline}; tokens = {}; itoken = 0;
16      while isstr(rem)      while isstr(rem)
# Line 21  for iline = 1:length(file) Line 27  for iline = 1:length(file)
27          inc = str2num(tokens{5});          inc = str2num(tokens{5});
28          xax = [ini:inc:ini+(num-1)*inc];          xax = [ini:inc:ini+(num-1)*inc];
29          if min(xax) >= 0 && max(xax) > 180          if min(xax) >= 0 && max(xax) > 180
30              xax = xax - 180;              ShiftData = 1;
31          end          end
32          nx = length(xax);          nx = length(xax);
33      end      end
# Line 60  for iline = 1:length(file) Line 66  for iline = 1:length(file)
66          end          end
67          if length(ini) == 13          if length(ini) == 13
68              monchar=ini(9:11);              monchar=ini(9:11);
69            elseif length(ini)==7
70                monchar=ini(1:3);
71            else
72             error('Cannot parse TDEF correctly');
73          end          end
74          if     isequal(monchar,'JAN'), inimonth = 1;          if     isequal(monchar,'JAN'), inimonth = 1;
75          elseif isequal(monchar,'FEB'), inimonth = 2;          elseif isequal(monchar,'FEB'), inimonth = 2;
# Line 74  for iline = 1:length(file) Line 84  for iline = 1:length(file)
84          elseif isequal(monchar,'NOV'), inimonth = 11;          elseif isequal(monchar,'NOV'), inimonth = 11;
85          elseif isequal(monchar,'DEC'), inimonth = 12; end          elseif isequal(monchar,'DEC'), inimonth = 12; end
86          num = str2num(tokens{2});          num = str2num(tokens{2});
87          months=[inimonth:num]; time=months/12; nt = length(months);          months=[inimonth:num+inimonth-1];
88            time=months/12; nt = length(months);
89                    
90      end      end
91            
# Line 82  for iline = 1:length(file) Line 93  for iline = 1:length(file)
93          nv = str2num(tokens{2});          nv = str2num(tokens{2});
94          VARSline = iline;          VARSline = iline;
95      end      end
       
96      if isequal(tokens{1},fln)      if isequal(tokens{1},fln)
97          FIELDline = iline;          FIELDline = iline;
98          index = iline - VARSline;          ivar = iline - VARSline;
99      end      end
100            
101      if isequal(tokens{1},'UNDEF')      if isequal(tokens{1},'UNDEF')
102          undef = str2num(tokens{2});          undef = str2num(tokens{2});
103      end      end
104            
105        if isequal(tokens{1},'FORMAT')
106            if isequal(tokens{2},'SEQUENTIAL') | isequal(tokens{2},'sequential')
107                format = 'SEQUENTIAL';
108            else
109                disp(['Unrecognized grads FORMAT:  ',tokens{2}]);
110            end
111        end
112        
113        if isequal(tokens{1},'DSET')
114            datafile = [dad,tokens{2}(2:end)];
115        end  
116        
117  end  end
118    
119    
120    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121    %                       Verification and debugging                        %
122    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123    
124    if DiagDebug, disp(['  Debug -- grads data format:  ',format]); end
125    
126  % Verify that things are there.  % Verify that things are there.
127  parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','index'};  parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','ivar'};
128  GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln    };  GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln   };
129  for ii = 1:length(parameters)  for ii = 1:length(parameters)
130      try      try
131          eval([parameters{ii},';']);          eval([parameters{ii},';']);
# Line 105  for ii = 1:length(parameters) Line 134  for ii = 1:length(parameters)
134      end      end
135  end  end
136    
137  % Read in grads data using axis infomation.  
138    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139    %                              Read data file                             %
140    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141    
142  tol = 1e-5;  tol = 1e-5;
143  fid=fopen(datafile,'r','b');  fid=fopen(datafile,'r','b');
144  data = fread(fid,'real*4');  data = fread(fid,'real*4');
145  fclose(fid);  fclose(fid);
146    if isequal(format,'SEQUENTIAL') | isequal(format,'sequential')
147        index=true([nx*ny*nz*nv*nt+2*nv*nt,1]);
148        index([1:nx*ny*nz+2:end])=false;
149        index([2:nx*ny*nz+2:end])=false;
150        data = data(index);
151    end
152  data = reshape(data,[nx,ny,nz,nv,nt]);  data = reshape(data,[nx,ny,nz,nv,nt]);
 data = squeeze(data(:,:,:,index,:));  
 data( abs((data-undef)/undef) < tol ) = NaN;  
153    data = squeeze(data(:,:,:,ivar,:));
154    data( abs((data-undef)/undef) < tol ) = NaN;
155    if ShiftData
156        indexWestHemi = xax>=180;
157        indexEastHemi = xax<180;
158        data = cat(1,data(indexWestHemi,:,:),data(indexEastHemi,:,:));
159        xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
160    end
161    datatmp(:,:,inimonth:num+inimonth-1) = data;
162    data=datatmp;
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22