%
%  %Id:  $
%
%  Ed Hill
%
%  Create a netCDF file defining a global Lat-Lon grid for SCRIP.

clear all
close all

debug_lev = 0;

%
%  JMC wants: 128x64
%
nlat = 64;
nlon = 128;
latcell = 180/nlat;
loncell = 360/nlon;
lat = linspace(-90+latcell/2, 90-latcell/2, nlat);
%
%  Note: the APE grid needs to be [0, MAX_EAST]
%   lon = linspace(-180+loncell/2, 180-loncell/2, nlon);
lon = linspace(0+loncell/2, 360-loncell/2, nlon);

dlat_t = diff(lat);
dlon_t = diff(lon);
dlat = 0.5*( [dlat_t(1) ; dlat_t(:)] + ...
             [dlat_t(:) ; dlat_t(length(dlat_t))] );
dlon = 0.5*( [dlon_t(1) ; dlon_t(:)] + ...
             [dlon_t(:) ; dlon_t(length(dlon_t))] );
clear dlat_t dlon_t

grid_size = nlat*nlon;
grid_rank = 2;
grid_corners = 4;

! rm -f remap_grid_LL128x64.cdf

nc = netcdf(['remap_grid_LL128x64.cdf'], 'clobber');
nc.title = 'Global 2-degree Lat-Lon 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'}(:) = [ length(lon) length(lat) ];

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 ilon = 1:length(lon)
  for ilat = 1:length(lat)
    j = j + 1;
    latc = lat(ilat);
    lonc = lon(ilon);
    cen_lat(j) = latc;
    cen_lon(j) = lonc;
    %  Note:  grid points *MUST* be ordered CCW and here
    %    we start from the lower-left corner
    da = dlat(ilat)/2.0;
    do = dlon(ilon)/2.0;
    cor_lat(j,:) = [ latc-da latc-da latc+da latc+da ];
    cor_lon(j,:) = [ lonc-do lonc+do lonc+do lonc-do ];
%     if debug_lev > 1
%       pause(1)
%       hold on
%       plot([cor_lon(j,:) cor_lon(j,1)],...
%            [cor_lat(j,:) cor_lat(j,1)],'o-g')
%       plot(cen_lon(j),cen_lat(j),'xr')
%       hold off
%     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)