%
%  %Id:  $
%
%  Ed Hill
%
%  Create a netCDF file defining a CS32 grid for SCRIP.

clear all
close all

debug_lev = 0;

%==================================================================
%  Read the tile00?.mitgrid files
gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
          'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
ne = 32;
nep1 = ne + 1;
iface = 1;
for iface = 1:6
  fname = sprintf('CS32/tile%03d.mitgrid', iface);
  gid = fopen(fname, 'r', 'ieee-be');
  tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
  fclose(gid);
  %  surf(tmp(:,:,1)), view(2), shading interp
  %  for jj = 1:length(gvars)
  for jj = 1:7
    comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
                   [gvars{jj}], iface, jj);
    eval(comm);
  end
end
%  surf(XC(:,:,1)), view(2), shading interp
%  subplot(2,1,1),  a = [1:10];     surf(XC(a,a,1)), view(2)
%  subplot(2,1,2),  a = [(nep1-10):nep1];  surf(XC(a,a,1)), view(2)
%  surf(YC(:,:,1)), view(2), shading interp
%  surf(XG(:,:,1)), view(2), shading interp
%  surf(YG(:,:,1)), view(2), shading interp
%==================================================================


grid_size = ne * ne * 6;
grid_rank = 2;
grid_corners = 4;

! rm -f remap_grid_CS32.cdf

nc = netcdf(['remap_grid_CS32.cdf'], 'clobber');
nc.title = 'Global CS32 grid';
nc.author = 'Ed Hill';
nc('grid_size') = grid_size;
nc('grid_rank') = grid_rank;
nc('grid_corners') = grid_corners;

nc{'grid_dims'} = ncint('grid_rank');
nc{'grid_dims'}(:) = [ ne 6*ne ];

nc{'grid_imask'} = ncint('grid_size');
nc{'grid_imask'}(:) = ones(length(grid_size),1);

nc{'grid_center_lat'} = ncdouble('grid_size');
nc{'grid_center_lon'} = ncdouble('grid_size');
nc{'grid_center_lat'}.units = 'degrees';
nc{'grid_center_lon'}.units = 'degrees';
nc{'grid_corner_lat'} = ncdouble('grid_size', 'grid_corners');
nc{'grid_corner_lon'} = ncdouble('grid_size', 'grid_corners');
nc{'grid_corner_lat'}.units = 'degrees';
nc{'grid_corner_lon'}.units = 'degrees';

cen_lat = zeros(grid_size,1);
cen_lon = zeros(grid_size,1);
cor_lat = zeros(grid_size,4);
cor_lon = zeros(grid_size,4);
j = 0;
for k = 1:6
  for iy = 1:ne
    for ix = 1:ne
      j = j + 1;
      cen_lat(j) = YC(ix,iy,k);
      cen_lon(j) = XC(ix,iy,k);
      %  Note:  grid points *MUST* be ordered CCW and here
      %    we start from the lower-left corner
      a = [ YG(ix,iy,k) YG(ix+1,iy,k) YG(ix+1,iy+1,k) YG(ix,iy+1,k) ];
      b = [ XG(ix,iy,k) XG(ix+1,iy,k) XG(ix+1,iy+1,k) XG(ix,iy+1,k) ];
      %  fixup the poles
      for ii = 1:4
        if (abs(b(ii)-XC(ix,iy,k)) > 150)
          b(ii) = b(ii) + sign(XC(ix,iy,k))*360;
        end
        if (90 - abs(a(ii))) < 0.01
          b(ii) = XC(ix,iy,k);
        end
      end
      cor_lat(j,:) = a;
      cor_lon(j,:) = b;
%       if debug_lev > 1
%         pause(1)
%         hold on
%         plot([cor_lon(j,:) cor_lon(j,1)],...
%              [cor_lat(j,:) cor_lat(j,1)],'o-b')
%         plot(cen_lon(j),cen_lat(j),'xr')
%         hold off
%       end
    end
  end
end
nc{'grid_center_lat'}(:) = cen_lat;
nc{'grid_center_lon'}(:) = cen_lon;
nc{'grid_corner_lat'}(:) = cor_lat;
nc{'grid_corner_lon'}(:) = cor_lon;

close(nc)