enderton/Diagnostics/DiagLoadGradsData.m

function [data,xax,yax,zax,months,time,dim] = ...
    DiagLoadGradsData(Grads,dad,fln,DiagDebug)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          Parse table file                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ShiftData = 0;
format = 'NONSEQUENTIAL';

tablfile = [dad,'/',Grads,'.tabl'];
file = textread(tablfile,'%s','delimiter','\n','whitespace','');

for iline = 1:length(file)
    rem = file{iline}; tokens = {}; itoken = 0;
    while isstr(rem)
        [token,rem] = strtok(rem); itoken = itoken + 1;
        if ~isempty(token), tokens{itoken} = token; end
    end
    
    if isequal(tokens{1},'XDEF')
        if ~isequal(tokens{3},'LINEAR')
            error('Unable to deal with nonlinear grads x-axis data!');
        end
        num = str2num(tokens{2});
        ini = str2num(tokens{4});
        inc = str2num(tokens{5});
        xax = [ini:inc:ini+(num-1)*inc];
        if min(xax) >= 0 && max(xax) > 180
            ShiftData = 1;
        end
        nx = length(xax);
    end
    
    if isequal(tokens{1},'YDEF')
        if ~isequal(tokens{3},'LINEAR')
            error('Unable to deal with nonlinear grads y-axis data!');
        end
        num = str2num(tokens{2});
        ini = str2num(tokens{4});
        inc = str2num(tokens{5});
        yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);
    end
    
    if isequal(tokens{1},'ZDEF')
        if isequal(tokens{2},'1')
            zax = [0]; zax_found = 1; dim = 2;
        else
            if ~isequal(tokens{3},'LINEAR')
                error('Unable to deal with nonlinear grads z-axis data!');
            end
            num = str2num(tokens{2});
            ini = str2num(tokens{4});
            inc = str2num(tokens{5});
            zax = [ini:inc:ini+(num-1)*inc]; dim = 3;
        end
        nz = length(zax);
    end
    
    if isequal(tokens{1},'TDEF'), ini=tokens{4};
        if ~isequal(tokens{3},'LINEAR')
            error('Currently unable to deal with nonlinear grads z-axis data!');
        end
        if ~isequal(tokens{5},'1mo')
            error('Currently unable to deal with non monthly mean grads data!');
        end
        if length(ini) == 13
            monchar=ini(9:11);
        end
        if     isequal(monchar,'JAN'), inimonth = 1;
        elseif isequal(monchar,'FEB'), inimonth = 2;
        elseif isequal(monchar,'MAR'), inimonth = 3;
        elseif isequal(monchar,'APR'), inimonth = 4;
        elseif isequal(monchar,'MAY'), inimonth = 5;
        elseif isequal(monchar,'JUN'), inimonth = 6;
        elseif isequal(monchar,'JUL'), inimonth = 7;
        elseif isequal(monchar,'AUG'), inimonth = 8;
        elseif isequal(monchar,'SEP'), inimonth = 9;
        elseif isequal(monchar,'OCT'), inimonth = 10;
        elseif isequal(monchar,'NOV'), inimonth = 11;
        elseif isequal(monchar,'DEC'), inimonth = 12; end
        num = str2num(tokens{2});
        months=[inimonth:num+inimonth-1];
        time=months/12; nt = length(months);
        
    end
    
    if isequal(tokens{1},'VARS')
        nv = str2num(tokens{2});
        VARSline = iline;
    end
    if isequal(tokens{1},fln)
        FIELDline = iline;
        ivar = iline - VARSline;
    end
    
    if isequal(tokens{1},'UNDEF')
        undef = str2num(tokens{2});
    end
    
    if isequal(tokens{1},'FORMAT')
        if isequal(tokens{2},'SEQUENTIAL')
            format = 'SEQUENTIAL';
        else
            disp(['Unrecognized grads FORMAT:  ',tokens{2}]);
        end
    end 
    
    if isequal(tokens{1},'DSET')
        datafile = [dad,tokens{2}(2:end)];
    end  
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Verification and debugging                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if DiagDebug, disp(['  Debug -- grads data format:  ',format]); end

% Verify that things are there.
parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','ivar'};
GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln   };
for ii = 1:length(parameters)
    try
        eval([parameters{ii},';']);
    catch
        error(['GRADS information not found:  ',GRADSnames{ii}]);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              Read data file                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tol = 1e-5;
fid=fopen(datafile,'r','b');
data = fread(fid,'real*4');
fclose(fid);
if isequal(format,'SEQUENTIAL')
    index=true([nx*ny*nz*nv*nt+2*nv*nt,1]);
    index([1:nx*ny*nz+2:end])=false;
    index([2:nx*ny*nz+2:end])=false;
    data = data(index);
end
data = reshape(data,[nx,ny,nz,nv,nt]);
data = squeeze(data(:,:,:,ivar,:));
data( abs((data-undef)/undef) < tol ) = NaN;
if ShiftData
    indexWestHemi = xax>=180;
    indexEastHemi = xax<180;
    data = cat(1,data(indexWestHemi,:,:),data(indexEastHemi,:,:));
    xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
end
1	enderton	1.1	function [data,xax,yax,zax,months,time,dim] = ...
2	enderton	1.2	DiagLoadGradsData(Grads,dad,fln,DiagDebug)
3
4
5			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6			% Parse table file %
7			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8	molod	1.3	ShiftData = 0;
9	enderton	1.2	format = 'NONSEQUENTIAL';
10	enderton	1.1
11			tablfile = [dad,'/',Grads,'.tabl'];
12			file = textread(tablfile,'%s','delimiter','\n','whitespace','');
13	enderton	1.2
14	enderton	1.1	for iline = 1:length(file)
15			rem = file{iline}; tokens = {}; itoken = 0;
16			while isstr(rem)
17			[token,rem] = strtok(rem); itoken = itoken + 1;
18			if ~isempty(token), tokens{itoken} = token; end
19			end
20
21			if isequal(tokens{1},'XDEF')
22			if ~isequal(tokens{3},'LINEAR')
23			error('Unable to deal with nonlinear grads x-axis data!');
24			end
25			num = str2num(tokens{2});
26			ini = str2num(tokens{4});
27			inc = str2num(tokens{5});
28			xax = [ini:inc:ini+(num-1)*inc];
29			if min(xax) >= 0 && max(xax) > 180
30	molod	1.3	ShiftData = 1;
31	enderton	1.1	end
32			nx = length(xax);
33			end
34
35			if isequal(tokens{1},'YDEF')
36			if ~isequal(tokens{3},'LINEAR')
37			error('Unable to deal with nonlinear grads y-axis data!');
38			end
39			num = str2num(tokens{2});
40			ini = str2num(tokens{4});
41			inc = str2num(tokens{5});
42			yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);
43			end
44
45			if isequal(tokens{1},'ZDEF')
46			if isequal(tokens{2},'1')
47			zax = [0]; zax_found = 1; dim = 2;
48			else
49			if ~isequal(tokens{3},'LINEAR')
50			error('Unable to deal with nonlinear grads z-axis data!');
51			end
52			num = str2num(tokens{2});
53			ini = str2num(tokens{4});
54			inc = str2num(tokens{5});
55			zax = [ini:inc:ini+(num-1)*inc]; dim = 3;
56			end
57			nz = length(zax);
58			end
59
60			if isequal(tokens{1},'TDEF'), ini=tokens{4};
61			if ~isequal(tokens{3},'LINEAR')
62			error('Currently unable to deal with nonlinear grads z-axis data!');
63			end
64			if ~isequal(tokens{5},'1mo')
65			error('Currently unable to deal with non monthly mean grads data!');
66			end
67			if length(ini) == 13
68			monchar=ini(9:11);
69			end
70			if isequal(monchar,'JAN'), inimonth = 1;
71			elseif isequal(monchar,'FEB'), inimonth = 2;
72			elseif isequal(monchar,'MAR'), inimonth = 3;
73			elseif isequal(monchar,'APR'), inimonth = 4;
74			elseif isequal(monchar,'MAY'), inimonth = 5;
75			elseif isequal(monchar,'JUN'), inimonth = 6;
76			elseif isequal(monchar,'JUL'), inimonth = 7;
77			elseif isequal(monchar,'AUG'), inimonth = 8;
78			elseif isequal(monchar,'SEP'), inimonth = 9;
79			elseif isequal(monchar,'OCT'), inimonth = 10;
80			elseif isequal(monchar,'NOV'), inimonth = 11;
81			elseif isequal(monchar,'DEC'), inimonth = 12; end
82			num = str2num(tokens{2});
83	enderton	1.2	months=[inimonth:num+inimonth-1];
84			time=months/12; nt = length(months);
85	enderton	1.1
86			end
87
88			if isequal(tokens{1},'VARS')
89			nv = str2num(tokens{2});
90			VARSline = iline;
91			end
92			if isequal(tokens{1},fln)
93			FIELDline = iline;
94	enderton	1.2	ivar = iline - VARSline;
95	enderton	1.1	end
96
97			if isequal(tokens{1},'UNDEF')
98			undef = str2num(tokens{2});
99			end
100
101	enderton	1.2	if isequal(tokens{1},'FORMAT')
102			if isequal(tokens{2},'SEQUENTIAL')
103			format = 'SEQUENTIAL';
104			else
105			disp(['Unrecognized grads FORMAT: ',tokens{2}]);
106			end
107			end
108
109			if isequal(tokens{1},'DSET')
110			datafile = [dad,tokens{2}(2:end)];
111			end
112
113	enderton	1.1	end
114
115	enderton	1.2
116			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117			% Verification and debugging %
118			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119
120			if DiagDebug, disp([' Debug -- grads data format: ',format]); end
121
122	enderton	1.1	% Verify that things are there.
123	enderton	1.2	parameters = {'xax' ,'yax' ,'zax' ,'time','nv' ,'undef','ivar'};
124			GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln };
125	enderton	1.1	for ii = 1:length(parameters)
126			try
127			eval([parameters{ii},';']);
128			catch
129			error(['GRADS information not found: ',GRADSnames{ii}]);
130			end
131			end
132
133	enderton	1.2
134			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135			% Read data file %
136			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137
138	enderton	1.1	tol = 1e-5;
139			fid=fopen(datafile,'r','b');
140			data = fread(fid,'real*4');
141			fclose(fid);
142	enderton	1.2	if isequal(format,'SEQUENTIAL')
143			index=true([nxnynznvnt+2nvnt,1]);
144			index([1:nxnynz+2:end])=false;
145			index([2:nxnynz+2:end])=false;
146			data = data(index);
147			end
148	enderton	1.1	data = reshape(data,[nx,ny,nz,nv,nt]);
149	enderton	1.2	data = squeeze(data(:,:,:,ivar,:));
150	molod	1.3	data( abs((data-undef)/undef) < tol ) = NaN;
151			if ShiftData
152			indexWestHemi = xax>=180;
153			indexEastHemi = xax<180;
154			data = cat(1,data(indexWestHemi,:,:),data(indexEastHemi,:,:));
155			xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
156			end