/[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.4 by molod, Tue Jun 28 21:33:51 2005 UTC revision 1.5 by enderton, Mon Sep 5 18:48:27 2005 UTC
# Line 2  function [data,xax,yax,zax,months,time,d Line 2  function [data,xax,yax,zax,months,time,d
2      DiagLoadGradsData(Grads,dad,fln,DiagDebug)      DiagLoadGradsData(Grads,dad,fln,DiagDebug)
3    
4    
5  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  % Initial assumptions about data.
 %                          Parse table file                               %  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
6  ShiftData = 0;  ShiftData = 0;
7    MultFiles = 0;
8  format = 'NONSEQUENTIAL';  format = 'NONSEQUENTIAL';
9    
10    
11    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
12    %                          Parse table file                               %
13    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14  tablfile = [dad,'/',Grads];  tablfile = [dad,'/',Grads];
15  file = textread(tablfile,'%s','delimiter','\n','whitespace','');  file = textread(tablfile,'%s','delimiter','\n','whitespace','');
16    
17    % Cycle though lines of the GRADS table file, pasring the file line-by-line
18    % to properly read in and format the data.  The first line of every file
19    % has an identifier denoting the information on that line.  The following
20    % blocks of code parse the information according to different identifiers.
21  for iline = 1:length(file)  for iline = 1:length(file)
22      rem = file{iline}; tokens = {}; itoken = 0;      rem = file{iline}; tokens = {}; itoken = 0;
23      while isstr(rem)      while isstr(rem)
# Line 18  for iline = 1:length(file) Line 25  for iline = 1:length(file)
25          if ~isempty(token), tokens{itoken} = token; end          if ~isempty(token), tokens{itoken} = token; end
26      end      end
27            
28        % Determine file(s) to read in.
29        if isequal(tokens{1},'DSET')
30            if isempty(findstr(tokens{2},'%'))
31                datafile = [dad,tokens{2}(2:end)];
32            else
33                str = tokens{2}(2:end);
34                alldatafile = strrep(str,'%y2','*');
35                alldatafile = strrep(alldatafile,'%m2','*');
36                files = ls([dad,alldatafile]);
37                breaks = [0,find(isspace(files))];
38                for ibreak = 1:length(breaks)-1
39                    datafile{ibreak} = files([breaks(ibreak)+1:breaks(ibreak+1)-1]);
40                end
41                MultFiles = 1;
42            end
43        end  
44        
45        % Read in x-axis information.
46      if isequal(tokens{1},'XDEF')      if isequal(tokens{1},'XDEF')
47          if ~isequal(tokens{3},'LINEAR')          if ~isequal(tokens{3},'LINEAR')
48              error('Unable to deal with nonlinear grads x-axis data!');              error('Unable to deal with nonlinear grads x-axis data!');
# Line 26  for iline = 1:length(file) Line 51  for iline = 1:length(file)
51          ini = str2num(tokens{4});          ini = str2num(tokens{4});
52          inc = str2num(tokens{5});          inc = str2num(tokens{5});
53          xax = [ini:inc:ini+(num-1)*inc];          xax = [ini:inc:ini+(num-1)*inc];
54            % If the x-axis is from 0 to 360, the axis information and the data
55            % will later have to be shifted to a -180 to 180 axis.
56          if min(xax) >= 0 && max(xax) > 180          if min(xax) >= 0 && max(xax) > 180
57              ShiftData = 1;              ShiftData = 1;
58          end          end
59          nx = length(xax);          nx = length(xax);
60      end      end
61            
62        % Read in y-axis information.
63      if isequal(tokens{1},'YDEF')      if isequal(tokens{1},'YDEF')
64          if ~isequal(tokens{3},'LINEAR')          if ~isequal(tokens{3},'LINEAR')
65              error('Unable to deal with nonlinear grads y-axis data!');              error('Unable to deal with nonlinear grads y-axis data!');
# Line 42  for iline = 1:length(file) Line 70  for iline = 1:length(file)
70          yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);          yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);
71      end      end
72            
73        % Read in z-axis information.
74      if isequal(tokens{1},'ZDEF')      if isequal(tokens{1},'ZDEF')
75          if isequal(tokens{2},'1')          if isequal(tokens{2},'1')
76              zax = [0]; zax_found = 1; dim = 2;              zax = [0]; zax_found = 1; dim = 2;
77          else          else
78              if ~isequal(tokens{3},'LINEAR')              dim = 3;
79                  error('Unable to deal with nonlinear grads z-axis data!');              if isequal(tokens{3},'LINEAR')
80                    num = str2num(tokens{2});
81                    ini = str2num(tokens{4});
82                    inc = str2num(tokens{5});
83                    zax = [ini:inc:ini+(num-1)*inc];
84                elseif isequal(tokens{3},'LEVELS')
85                    num = str2num(tokens{2});
86                    for inum = 1:num
87                        zax(inum) = 100.*str2num(tokens{inum+3});
88                    end
89                else
90                    error('Unknown information in z-axis data!');
91              end              end
             num = str2num(tokens{2});  
             ini = str2num(tokens{4});  
             inc = str2num(tokens{5});  
             zax = [ini:inc:ini+(num-1)*inc]; dim = 3;  
92          end          end
93          nz = length(zax);          nz = length(zax);
94      end      end
95            
96        % Read in t-axis (time) information.
97      if isequal(tokens{1},'TDEF'), ini=tokens{4};      if isequal(tokens{1},'TDEF'), ini=tokens{4};
98          if ~isequal(tokens{3},'LINEAR')          if ~isequal(tokens{3},'LINEAR')
99              error('Currently unable to deal with nonlinear grads z-axis data!');              error('Currently unable to deal with nonlinear grads t-axis data!');
100          end          end
101          if ~isequal(tokens{5},'1mo')          if ~isequal(tokens{5},'1mo') & ~isequal(tokens{5},'01mo')
102                % Note that monthly mean assumptions are also made when data
103                % spanning multiple files are read in later on and that this
104                % part must be updated if non-monthly mean data is allowed.
105              error('Currently unable to deal with non monthly mean grads data!');              error('Currently unable to deal with non monthly mean grads data!');
106          end          end
107          if length(ini) == 13          % Determine initial month and initial year.  This is done in an ad
108              monchar=ini(9:11);          % hoc manner in that depending on the length start time string, the
109          elseif length(ini)==7          % month is pluck out.  There would be problems is there are start
110              monchar=ini(1:3);          % time strings of different lengths have months in different
111          else          % locations.  In this case the month will not be found and an error
112           error('Cannot parse TDEF correctly');          % will be thrown.  The year is only required for data spanning
113          end          % multiple files, which thus fas has a start time string length of
114            % 12.  If it is not found, and error of undefined yrchar will be
115            % cast and the code will need to be appropriately updated.
116            if     length(ini) == 13, monchar=ini(9:11);
117            elseif length(ini) == 12, monchar=ini(6:8);  yrchar = ini(9:12);
118            elseif length(ini) == 10, monchar=ini(6:8);
119            elseif length(ini) == 7,  monchar=ini(1:3);
120            else, error('Cannot parse TDEF correctly'); end
121            monchar = upper(monchar);
122          if     isequal(monchar,'JAN'), inimonth = 1;          if     isequal(monchar,'JAN'), inimonth = 1;
123          elseif isequal(monchar,'FEB'), inimonth = 2;          elseif isequal(monchar,'FEB'), inimonth = 2;
124          elseif isequal(monchar,'MAR'), inimonth = 3;          elseif isequal(monchar,'MAR'), inimonth = 3;
# Line 82  for iline = 1:length(file) Line 130  for iline = 1:length(file)
130          elseif isequal(monchar,'SEP'), inimonth = 9;          elseif isequal(monchar,'SEP'), inimonth = 9;
131          elseif isequal(monchar,'OCT'), inimonth = 10;          elseif isequal(monchar,'OCT'), inimonth = 10;
132          elseif isequal(monchar,'NOV'), inimonth = 11;          elseif isequal(monchar,'NOV'), inimonth = 11;
133          elseif isequal(monchar,'DEC'), inimonth = 12; end          elseif isequal(monchar,'DEC'), inimonth = 12;
134            else, error('Can''t determine initial month in GRADS data!'); end
135          num = str2num(tokens{2});          num = str2num(tokens{2});
136            % Determine num and inimonth parameters accounting for the unknown
137            % number of data files present when dealing with multiple data
138            % files.  Assumption of monthly mean data is made.
139            if MultFiles;
140                initFileFound = 0;
141                for ifile = 1:length(datafile)
142                    initFileFound = ~isempty(findstr(yrchar,datafile{ifile})) ...
143                                    | initFileFound;
144                end
145                if ~initFileFound, inimonth = 1; end
146                num = 12.*length(datafile);
147            end
148            % Assumption of monthly mean data is made.
149          months=[inimonth:num+inimonth-1];          months=[inimonth:num+inimonth-1];
150          time=months/12; nt = length(months);          time=months/12; nt = length(months);
           
151      end      end
152            
153        % Find where the desired variable is located.
154      if isequal(tokens{1},'VARS')      if isequal(tokens{1},'VARS')
155          nv = str2num(tokens{2});          nv = str2num(tokens{2});
156          VARSline = iline;          VARSline = iline;
# Line 98  for iline = 1:length(file) Line 160  for iline = 1:length(file)
160          ivar = iline - VARSline;          ivar = iline - VARSline;
161      end      end
162            
163        % Find the values defonting undefined data.
164      if isequal(tokens{1},'UNDEF')      if isequal(tokens{1},'UNDEF')
165          undef = str2num(tokens{2});          undef = str2num(tokens{2});
166      end      end
167            
168        % Determine data format type.
169      if isequal(tokens{1},'FORMAT')      if isequal(tokens{1},'FORMAT')
170          if isequal(tokens{2},'SEQUENTIAL') | isequal(tokens{2},'sequential')          if isequal(tokens{2},'SEQUENTIAL') | isequal(tokens{2},'sequential')
171              format = 'SEQUENTIAL';              format = 'SEQUENTIAL';
# Line 109  for iline = 1:length(file) Line 173  for iline = 1:length(file)
173              disp(['Unrecognized grads FORMAT:  ',tokens{2}]);              disp(['Unrecognized grads FORMAT:  ',tokens{2}]);
174          end          end
175      end      end
       
     if isequal(tokens{1},'DSET')  
         datafile = [dad,tokens{2}(2:end)];  
     end    
       
176  end  end
177    
178    
# Line 121  end Line 180  end
180  %                       Verification and debugging                        %  %                       Verification and debugging                        %
181  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182    
183    % Verify that everything we need from the table file has been read in.
184  if DiagDebug, disp(['  Debug -- grads data format:  ',format]); end  if DiagDebug, disp(['  Debug -- grads data format:  ',format]); end
   
 % Verify that things are there.  
185  parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','ivar'};  parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','ivar'};
186  GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln   };  GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln   };
187  for ii = 1:length(parameters)  for ii = 1:length(parameters)
# Line 139  end Line 197  end
197  %                              Read data file                             %  %                              Read data file                             %
198  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199    
200  tol = 1e-5;  % Read in data from a single file.
201  fid=fopen(datafile,'r','b');  if ~MultFiles
202  data = fread(fid,'real*4');      fid=fopen(datafile,'r','b');
203  fclose(fid);      data = fread(fid,'real*4');
204  if isequal(format,'SEQUENTIAL') | isequal(format,'sequential')      fclose(fid);
205      index=true([nx*ny*nz*nv*nt+2*nv*nt,1]);          if isequal(format,'SEQUENTIAL') | isequal(format,'sequential')
206      index([1:nx*ny*nz+2:end])=false;          index=true([nx*ny*nz*nv*nt+2*nv*nt,1]);
207      index([2:nx*ny*nz+2:end])=false;          index([1:nx*ny*nz+2:end])=false;
208      data = data(index);          index([2:nx*ny*nz+2:end])=false;
209            data = data(index);
210            end
211        data = reshape(data,[nx,ny,nz,nv,nt]);
212        data = data(:,:,:,ivar,:);
213        datatmp(:,:,:,:,inimonth:num+inimonth-1) = data;
214        datatmp(:,:,:,:,1:inimonth-1) = NaN;
215        data=datatmp;
216    % Read in data from multiple files.
217    else
218        for ifile = 1:length(datafile)
219            if DiagDebug, disp(['  Debug -- Reading grads ',...
220                                'data file:  ',datafile{ifile}]); end
221            fid=fopen(datafile{ifile},'r','b');
222            datatemp = fread(fid,'real*4');
223            fclose(fid);
224            % Assumption of monthly mean data is made.
225            if MultFiles, nt = 12; end
226            if initFileFound & (inimonth ~= 1)
227                error(['Currently not able to handle Grads data spanning ',...
228                        'multiple files not starting in January']);
229            end
230            datatemp = reshape(datatemp,[nx,ny,nz,nv,nt]);
231            datatemp = squeeze(datatemp(:,:,:,ivar,:));
232            data(1:nx,1:ny,1:nz,nt.*(ifile-1)+1:nt.*ifile) = datatemp;
233        end
234  end  end
235  data = reshape(data,[nx,ny,nz,nv,nt]);  data = squeeze(data);
236  data = squeeze(data(:,:,:,ivar,:));  
237    % Turn undefined values into NaN.
238    tol = 1e-5;
239  data( abs((data-undef)/undef) < tol ) = NaN;  data( abs((data-undef)/undef) < tol ) = NaN;
240    
241    % Shift data such that it is on a -180 to 180 longitude axis.
242  if ShiftData  if ShiftData
243      indexWestHemi = xax>=180;      indexWestHemi = xax>=180;
244      indexEastHemi = xax<180;      indexEastHemi = xax<180;
245      data = cat(1,data(indexWestHemi,:,:),data(indexEastHemi,:,:));      data = cat(1,data(indexWestHemi,:,:,:),data(indexEastHemi,:,:,:));
246      xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));      xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
247  end  end
 datatmp(:,:,inimonth:num+inimonth-1) = data;  
 data=datatmp;  
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22