| 1 | 
function [data,xax,yax,zax,months,time,dim] = ... | 
| 2 | 
    DiagLoadGradsData(Grads,dad,fln,DiagDebug) | 
| 3 | 
 | 
| 4 | 
 | 
| 5 | 
% Initial assumptions about data. | 
| 6 | 
ShiftData = 0; | 
| 7 | 
MultFiles = 0; | 
| 8 | 
format = 'NONSEQUENTIAL'; | 
| 9 | 
 | 
| 10 | 
 | 
| 11 | 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 12 | 
%                          Parse table file                               % | 
| 13 | 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 14 | 
tablfile = [dad,'/',Grads]; | 
| 15 | 
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) | 
| 22 | 
    rem = file{iline}; tokens = {}; itoken = 0; | 
| 23 | 
    while isstr(rem) | 
| 24 | 
        [token,rem] = strtok(rem); itoken = itoken + 1; | 
| 25 | 
        if ~isempty(token), tokens{itoken} = token; end | 
| 26 | 
    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') | 
| 47 | 
        if ~isequal(tokens{3},'LINEAR') | 
| 48 | 
            error('Unable to deal with nonlinear grads x-axis data!'); | 
| 49 | 
        end | 
| 50 | 
        num = str2num(tokens{2}); | 
| 51 | 
        ini = str2num(tokens{4}); | 
| 52 | 
        inc = str2num(tokens{5}); | 
| 53 | 
        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 | 
| 57 | 
            ShiftData = 1; | 
| 58 | 
        end | 
| 59 | 
        nx = length(xax); | 
| 60 | 
    end | 
| 61 | 
     | 
| 62 | 
    % Read in y-axis information. | 
| 63 | 
    if isequal(tokens{1},'YDEF') | 
| 64 | 
        if ~isequal(tokens{3},'LINEAR') | 
| 65 | 
            error('Unable to deal with nonlinear grads y-axis data!'); | 
| 66 | 
        end | 
| 67 | 
        num = str2num(tokens{2}); | 
| 68 | 
        ini = str2num(tokens{4}); | 
| 69 | 
        inc = str2num(tokens{5}); | 
| 70 | 
        yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax); | 
| 71 | 
    end | 
| 72 | 
     | 
| 73 | 
    % Read in z-axis information. | 
| 74 | 
    if isequal(tokens{1},'ZDEF') | 
| 75 | 
        if isequal(tokens{2},'1') | 
| 76 | 
            zax = [0]; zax_found = 1; dim = 2; | 
| 77 | 
        else | 
| 78 | 
            dim = 3; | 
| 79 | 
            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 | 
| 92 | 
        end | 
| 93 | 
        nz = length(zax); | 
| 94 | 
    end | 
| 95 | 
     | 
| 96 | 
    % Read in t-axis (time) information. | 
| 97 | 
    if isequal(tokens{1},'TDEF'), ini=tokens{4}; | 
| 98 | 
        if ~isequal(tokens{3},'LINEAR') | 
| 99 | 
            error('Currently unable to deal with nonlinear grads t-axis data!'); | 
| 100 | 
        end | 
| 101 | 
        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!'); | 
| 106 | 
        end | 
| 107 | 
        % Determine initial month and initial year.  This is done in an ad | 
| 108 | 
        % hoc manner in that depending on the length start time string, the | 
| 109 | 
        % month is pluck out.  There would be problems is there are start | 
| 110 | 
        % time strings of different lengths have months in different | 
| 111 | 
        % locations.  In this case the month will not be found and an error | 
| 112 | 
        % will be thrown.  The year is only required for data spanning | 
| 113 | 
        % 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; | 
| 123 | 
        elseif isequal(monchar,'FEB'), inimonth = 2; | 
| 124 | 
        elseif isequal(monchar,'MAR'), inimonth = 3; | 
| 125 | 
        elseif isequal(monchar,'APR'), inimonth = 4; | 
| 126 | 
        elseif isequal(monchar,'MAY'), inimonth = 5; | 
| 127 | 
        elseif isequal(monchar,'JUN'), inimonth = 6; | 
| 128 | 
        elseif isequal(monchar,'JUL'), inimonth = 7; | 
| 129 | 
        elseif isequal(monchar,'AUG'), inimonth = 8; | 
| 130 | 
        elseif isequal(monchar,'SEP'), inimonth = 9; | 
| 131 | 
        elseif isequal(monchar,'OCT'), inimonth = 10; | 
| 132 | 
        elseif isequal(monchar,'NOV'), inimonth = 11; | 
| 133 | 
        elseif isequal(monchar,'DEC'), inimonth = 12; | 
| 134 | 
        else, error('Can''t determine initial month in GRADS data!'); end | 
| 135 | 
        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]; | 
| 150 | 
        time=months/12; nt = length(months); | 
| 151 | 
    end | 
| 152 | 
     | 
| 153 | 
    % Find where the desired variable is located. | 
| 154 | 
    if isequal(tokens{1},'VARS') | 
| 155 | 
        nv = str2num(tokens{2}); | 
| 156 | 
        VARSline = iline; | 
| 157 | 
    end | 
| 158 | 
    if isequal(tokens{1},fln) | 
| 159 | 
        FIELDline = iline; | 
| 160 | 
        ivar = iline - VARSline; | 
| 161 | 
    end | 
| 162 | 
     | 
| 163 | 
    % Find the values defonting undefined data. | 
| 164 | 
    if isequal(tokens{1},'UNDEF') | 
| 165 | 
        undef = str2num(tokens{2}); | 
| 166 | 
    end | 
| 167 | 
     | 
| 168 | 
    % Determine data format type. | 
| 169 | 
    if isequal(tokens{1},'FORMAT') | 
| 170 | 
        if isequal(tokens{2},'SEQUENTIAL') | isequal(tokens{2},'sequential') | 
| 171 | 
            format = 'SEQUENTIAL'; | 
| 172 | 
        else | 
| 173 | 
            disp(['Unrecognized grads FORMAT:  ',tokens{2}]); | 
| 174 | 
        end | 
| 175 | 
    end  | 
| 176 | 
end | 
| 177 | 
 | 
| 178 | 
 | 
| 179 | 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 180 | 
%                       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 | 
| 185 | 
parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','ivar'}; | 
| 186 | 
GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln   }; | 
| 187 | 
for ii = 1:length(parameters) | 
| 188 | 
    try | 
| 189 | 
        eval([parameters{ii},';']); | 
| 190 | 
    catch | 
| 191 | 
        error(['GRADS information not found:  ',GRADSnames{ii}]); | 
| 192 | 
    end | 
| 193 | 
end | 
| 194 | 
 | 
| 195 | 
 | 
| 196 | 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 197 | 
%                              Read data file                             % | 
| 198 | 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 199 | 
 | 
| 200 | 
% Read in data from a single file. | 
| 201 | 
if ~MultFiles | 
| 202 | 
    fid=fopen(datafile,'r','b'); | 
| 203 | 
    data = fread(fid,'real*4'); | 
| 204 | 
    fclose(fid); | 
| 205 | 
        if isequal(format,'SEQUENTIAL') | isequal(format,'sequential') | 
| 206 | 
        index=true([nx*ny*nz*nv*nt+2*nv*nt,1]); | 
| 207 | 
        index([1:nx*ny*nz+2:end])=false; | 
| 208 | 
        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 | 
| 235 | 
data = squeeze(data); | 
| 236 | 
 | 
| 237 | 
% Turn undefined values into NaN. | 
| 238 | 
tol = 1e-5; | 
| 239 | 
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 | 
| 243 | 
    indexWestHemi = xax>=180; | 
| 244 | 
    indexEastHemi = xax<180; | 
| 245 | 
    data = cat(1,data(indexWestHemi,:,:,:),data(indexEastHemi,:,:,:)); | 
| 246 | 
    xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi)); | 
| 247 | 
end |