function csg_obj=load_cs_grid(nr,nb,ng,gdir)
% 
% Plot cube sphere grid lines
% e.g. load_cs_grid( 32, 32, 32,'float64');
% e.g. load_cs_grid(510,510,510,'float64');
%

% Set cartesian/unstructure list toplogy size
csg_obj.nx=0;
csg_obj.ny=0;

% Set cube topology size.
csg_obj.nr=nr;
csg_obj.nb=nb;
csg_obj.ng=ng;

% Set format used to read grid variables
% Input if 64-bit floats. The store in memory as 32-bit floats option below
% can be used to save memory when working with large grids.
gvfmt='float64';
gvprec='double';
gvfmt='float64=>float32';
gvprec='single';

% Calculate cube face sizes
csg_obj.fx=zeros(6,1);
csg_obj.fy=zeros(6,1);
for i=1:6
 if i == 1
  csg_obj.fx(i)=nb; csg_obj.fy(i)=ng;
 end
 if i == 2
  csg_obj.fx(i)=nr; csg_obj.fy(i)=ng;
 end
 if i == 3
  csg_obj.fx(i)=nr; csg_obj.fy(i)=nb;
 end
 if i == 4
  csg_obj.fx(i)=ng; csg_obj.fy(i)=nb;
 end
 if i == 5
  csg_obj.fx(i)=ng; csg_obj.fy(i)=nr;
 end
 if i == 6
  csg_obj.fx(i)=nb; csg_obj.fy(i)=nr;
 end
end

% Read in physical grid information
% This comes from the 6 tile files which hold
% sixteen terms
% XC , YC , DXF, DYF, RA , XG , YG , DXV, 
% DYU, RAZ, DXC, DYC, RAW, RAS, DXG, DYG
% that define the grid as follows:
% XC  - Cell center longitude
% YC  - Cell center latitude
% DXF - Cell face spacing in local X direction in meters and
%       passing through the cell center.
% DYF - 
  xcpos=1;   ycpos=2;  dxfpos=3;  dyfpos=4;
  rapos=5;   xgpos=6;   ygpos=7;  dxvpos=8;  
 dyupos=9; razpos=10; dxcpos=11; dycpos=12; 
rawpos=13; raspos=14; dxgpos=15; dygpos=16;
  csg_obj.xcpos  = xcpos;
  csg_obj.ycpos  = ycpos;
  csg_obj.dxfpos = dxfpos;
  csg_obj.dyfpos = dyfpos;
  csg_obj.rapos  = rapos;
  csg_obj.xgpos  = xgpos;
  csg_obj.ygpos  = ygpos;
  csg_obj.dxvpos = dxvpos;
  csg_obj.dyupos = dyupos;
  csg_obj.razpos = razpos;
  csg_obj.dxcpos = dxcpos;
  csg_obj.dycpos = dycpos;
  csg_obj.rawpos = rawpos;
  csg_obj.raspos = raspos;
  csg_obj.dxgpos = dxgpos;
  csg_obj.dygpos = dygpos;
 
% Each term is a set grid terms for each cube face for a 
% total of
% (nb+1)*(ng+1)*2+(nr+1)*(ng+1)*2+(nr+1)*(nb+1)*2 
% points
% 1 - create a dummy 2d array to hold the terms
nxd=nb+1+nr+1+ng+1+nb+1;
nyd=ng+1+nb+1+nr+1;
csg_obj.gridarr=cast(ones(nxd,nyd,16)*12345.6789,gvprec);
mask_val=csg_obj.gridarr(1,1,1);
% 1- Read tile*.mitgrid
%    holds terms for cube face 1, size (nb+1)*(ng+1)
for i= 1:6
 fx=csg_obj.fx(i);
 fy=csg_obj.fy(i);
 if i == 1
  csg_obj.ilog(1)=1;csg_obj.jlog(1)=1;
 end
 if i == 2
  csg_obj.ilog(2)=csg_obj.ilog(1)+fx+1;
  csg_obj.jlog(2)=1;
 end
 if i == 3
  csg_obj.ilog(3)=csg_obj.ilog(1)+fx+1;
  csg_obj.jlog(3)=csg_obj.jlog(1)+fy+1;
 end
 if i == 4
  csg_obj.ilog(4)=csg_obj.ilog(3)+fx+1;
  csg_obj.jlog(4)=csg_obj.jlog(1)+fy+1;
 end
 if i == 5
  csg_obj.ilog(5)=csg_obj.ilog(3)+fx+1;
  csg_obj.jlog(5)=csg_obj.jlog(4)+fy+1;
 end
 if i == 6
  csg_obj.ilog(6)=csg_obj.ilog(5)+fx+1;
  csg_obj.jlog(6)=csg_obj.jlog(4)+fx+1;
 end
 ilog=csg_obj.ilog(i);
 jlog=csg_obj.jlog(i);
 ihi=ilog+fx;jhi=jlog+fy;
 fnam=sprintf('%s/tile%3.3d.mitgrid',gdir,i);
 fid=fopen(fnam,'r','ieee-be');
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,xcpos )=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,ycpos )=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dxfpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dyfpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,rapos )=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,xgpos )=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,ygpos )=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dxvpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dyupos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,razpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dxcpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dycpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,rawpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,raspos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dxgpos)=reshape(tmp,[(fx+1) (fy+1)]);
 tmp=fread(fid,(fx+1)*(fy+1),gvfmt);
 csg_obj.gridarr(ilog:ihi,jlog:jhi,dygpos)=reshape(tmp,[(fx+1) (fy+1)]);
 fclose(fid);
 csg_obj.index_i_center(ilog:ihi,jlog:jhi)=cast(repmat([ilog:ihi]',1,jhi-jlog+1),'int32');
 csg_obj.index_j_center(ilog:ihi,jlog:jhi)=cast(repmat([jlog:jhi] ,ihi-ilog+1,1),'int32');
 csg_obj.face_center(ilog:ihi,jlog:jhi)=cast(i,'int32');
end
csg_obj.gridarr(csg_obj.gridarr==mask_val)=NaN;
csg_obj.active_mask=csg_obj.gridarr(:,:,1)*0.+1;

return