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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          Parse table file                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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
            xax = xax - 180;
        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;