function [A, lat, lon] = readBgrid(fname)
% function [A, lat, lon] = readBgrid(fname)
% INPUT
% fname     name of *.b grid file
%
% OUTPUT
% A         [lat x lon] grid of data; missing or 'no transformation'
%               values are now NaN
% lat       latitudes (decimal degrees) corresponding to 
%               rows of A
% lon       longitudes (decimal degrees) corresponding to 
%               columns of A


% open file
fid=fopen(fname,'r');
if fid<1
    error(sprintf('cannot open file %s',fname));
end

% read the header to extract information on how to read
% the remainder of the file
status = fseek(fid,4,'bof'); % skip first 4 bytes (Fortran header)
gridcorner = fread(fid,2,'double');
dlatlon = fread(fid,2,'double');
nlat = fread(fid,1,'long');
nlon = fread(fid,1,'long');
ikind = fread(fid,1,'int');
fseek(fid,4,'cof'); % end of header

% 'ikind' tells the type of data in remainder of file
switch ikind
    case -1
        fprintf('2 byte int\n');
        itype = 'uint16';
    case 0
        fprintf('4 byte int\n');
        itype = 'uint32';
    case 1
        fprintf('4 byte real\n');
        itype = 'single';
    case 2
        fprintf('2 byte int\n');
        itype = 'uint16';
    otherwise
        error('invalid ikind: %d',ikind);
end

% process records sequentially to fill data matrix
tic;
for i = 1:nlat
    fseek(fid,4,'cof'); % skip Fortran buffer
    A(i,:) = fread(fid,nlon,itype); % read 'nlon' values
    fseek(fid,4,'cof'); % skip Fortran buffer for end of record
end
toc;
fclose(fid);

% replace -999.00 values with NaN
q = A==-999.0;
A(q) = NaN;

% create lat/lon vectors
lat = gridcorner(1) + [0:dlatlon(1):dlatlon(1)*(nlat-1)];
lon = gridcorner(2) + [0:dlatlon(2):dlatlon(2)*(nlon-1)];

return
