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];
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);
        elseif length(ini)==7
            monchar=ini(1:3);
        else
         error('Cannot parse TDEF correctly');
        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') | 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') | 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
datatmp(:,:,inimonth:num+inimonth-1) = data;
data=datatmp;
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	molod	1.4	tablfile = [dad,'/',Grads];
12	enderton	1.1	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	molod	1.4	elseif length(ini)==7
70			monchar=ini(1:3);
71			else
72			error('Cannot parse TDEF correctly');
73	enderton	1.1	end
74			if isequal(monchar,'JAN'), inimonth = 1;
75			elseif isequal(monchar,'FEB'), inimonth = 2;
76			elseif isequal(monchar,'MAR'), inimonth = 3;
77			elseif isequal(monchar,'APR'), inimonth = 4;
78			elseif isequal(monchar,'MAY'), inimonth = 5;
79			elseif isequal(monchar,'JUN'), inimonth = 6;
80			elseif isequal(monchar,'JUL'), inimonth = 7;
81			elseif isequal(monchar,'AUG'), inimonth = 8;
82			elseif isequal(monchar,'SEP'), inimonth = 9;
83			elseif isequal(monchar,'OCT'), inimonth = 10;
84			elseif isequal(monchar,'NOV'), inimonth = 11;
85			elseif isequal(monchar,'DEC'), inimonth = 12; end
86			num = str2num(tokens{2});
87	enderton	1.2	months=[inimonth:num+inimonth-1];
88			time=months/12; nt = length(months);
89	enderton	1.1
90			end
91
92			if isequal(tokens{1},'VARS')
93			nv = str2num(tokens{2});
94			VARSline = iline;
95			end
96			if isequal(tokens{1},fln)
97			FIELDline = iline;
98	enderton	1.2	ivar = iline - VARSline;
99	enderton	1.1	end
100
101			if isequal(tokens{1},'UNDEF')
102			undef = str2num(tokens{2});
103			end
104
105	enderton	1.2	if isequal(tokens{1},'FORMAT')
106	molod	1.4	if isequal(tokens{2},'SEQUENTIAL') \| isequal(tokens{2},'sequential')
107	enderton	1.2	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	enderton	1.1	end
118
119	enderton	1.2
120			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121			% Verification and debugging %
122			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123
124			if DiagDebug, disp([' Debug -- grads data format: ',format]); end
125
126	enderton	1.1	% Verify that things are there.
127	enderton	1.2	parameters = {'xax' ,'yax' ,'zax' ,'time','nv' ,'undef','ivar'};
128			GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln };
129	enderton	1.1	for ii = 1:length(parameters)
130			try
131			eval([parameters{ii},';']);
132			catch
133			error(['GRADS information not found: ',GRADSnames{ii}]);
134			end
135			end
136
137	enderton	1.2
138			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139			% Read data file %
140			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141
142	enderton	1.1	tol = 1e-5;
143			fid=fopen(datafile,'r','b');
144			data = fread(fid,'real*4');
145			fclose(fid);
146	molod	1.4	if isequal(format,'SEQUENTIAL') \| isequal(format,'sequential')
147	enderton	1.2	index=true([nxnynznvnt+2nvnt,1]);
148			index([1:nxnynz+2:end])=false;
149			index([2:nxnynz+2:end])=false;
150			data = data(index);
151			end
152	enderton	1.1	data = reshape(data,[nx,ny,nz,nv,nt]);
153	enderton	1.2	data = squeeze(data(:,:,:,ivar,:));
154	molod	1.3	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	molod	1.4	end
161			datatmp(:,:,inimonth:num+inimonth-1) = data;
162			data=datatmp;