function grd = cms_read_grid(filename)
% CMS_READ_GRID Reads a CMS nonuniform Cartesian grid
%
% Description:
%   Reads a CMS grid file and creates
%   a structure variable containing the grid variables
%
% Usage:
%   grd = cms_read_grid(filename);
%
% Input:
%   filename - Input file name including full path (e.g. 'Flow_grd.h5')
%
% Output:
%   grd - structure variable containing the following variables:
%     nx - Grid size in x direction [-]
%     ny - Grid size in y direction [-]
%     x0 - Grid origin in x direction [m]
%     y0 - Grid origin in y direction [m]
%     angle - Grid angle [deg]
%     x(1:nx) - Cell-center coordinates in x direction [m]
%     y(1:ny) - Cell-center coordinates in y direction [m]
%     depth(1:ny,1:nx) - Cell-centered depth [m]
%     active(1:ny,1:nx) - Cell activity (1-active and 0-inactive) [logical]
%     Plus other optional datasets such as mannings, d50, etc.
%
% Author: Alex Sanchez, USACE-ERDC-CHL

info = h5info(filename);
s = info.Groups;
simlabel = s.Name(2:end);
a = s.Attributes;

for i=1:length(a)
    switch a(i).Name
        case 'Bearing'
            grd.angle = 360 - double(a(i).Value);
            if(grd.angle==360) 
                grd.angle = 0;
            end
        case 'Origin'
            grd.x0 = double(a(i).Value(1));
            grd.y0 = double(a(i).Value(2));
        case 'NumI'
            grd.nx = double(a(i).Value);
        case 'NumJ'
            grd.ny = double(a(i).Value);
    end
    
end
dpath = ['/',simlabel,'/'];
coordi = h5read(filename,[dpath,'CoordsI']);
coordj = h5read(filename,[dpath,'CoordsJ']);
grd.dx = coordi;
grd.dy = coordj;
grd.dx(2:end) = coordi(2:end) - coordi(1:end-1);
grd.dy(2:end) = coordj(2:end) - coordj(1:end-1);
grd.x = grd.x0 + coordi - grd.dx/2;
grd.y = grd.y0 + coordj - grd.dy/2;

% grd.depth = h5read(filename,[dpath,'Datasets/Depth/Values']);
for i=1:length(s.Groups)
    if strcmpi(s.Groups(i).Name,[dpath,'Datasets'])
        for j=1:length(s.Groups(i).Groups)
            datasetname = s.Groups(i).Groups(j).Name;
            ind = strfind(datasetname,'/');
            name = datasetname(ind(end)+1:end);
            %Remove/replace invalid characters for field names
            name = regexprep(name,' ','_');
            name = regexprep(name,'(','');
            name = regexprep(name,')','');
            if strcmpi(name,'ManningsN') %remove extra N
                name = 'mannings';
            end
            name = lower(name);
            grd.(name) = double(h5read(filename,[datasetname,'/Values']));
            grd.(name) = reshape(grd.(name),[grd.nx,grd.ny])';
        end
    end
end

grd.active = h5read(filename,[dpath,'PROPERTIES/CellTypes']);
grd.active = (grd.active==1);
grd.active = reshape(grd.active,[grd.nx,grd.ny])';
grd.depth(~grd.active) = NaN;
grd.depth(grd.depth==-999) = NaN;

return
end
