% latlongrid.m creates the model grid and stores the grid information
% various arrays. Nominally set are:
%    xc, yc, xg, yg, [2d]
%    surface area of grid-cell on sphere (zA) [2d]
%
% It also stores the above information in a MATLAB file GRID.mat
% usage:
%        load GRID.mat                   (loads everything above)
%        load GRID.mat nxc nyc xc yc     (loads just these variables/arrays)
%  
% Uses data files: n/a
% Creates data files: GRID.mat
% Creates arrays: xc,yc,xg,yg,xcf,ycf,zA

load DATAPATH 

periodic=1;

% read grid data
vnms = { 'XC' 'YC' 'XG' 'YG' 'rA'};
d2r = pi/180.0;
ginfo = {};
for k = 1:5
  nc = netcdf(fullfile(gridpath,sprintf(fnamepattern,k)),'nowrite');
  for j = 1:length(vnms)
    eval(sprintf('ginfo(k).%s = nc{''%s''}(:);',vnms{j},vnms{j}));
  end
  nc = close(nc);
end

% assemble grid
xc = [ginfo(1).XC;ginfo(2).XC;ginfo(4).XC(end:-1:1,:)'; ...
      ginfo(5).XC(end:-1:1,:)'];
yc = [ginfo(1).YC;ginfo(2).YC;ginfo(4).YC(end:-1:1,:)'; ...
      ginfo(5).YC(end:-1:1,:)'];
xg = [ginfo(1).XG(1:end-1,:);ginfo(2).XG(1:end-1,:); ...
      ginfo(4).XG(end:-1:1,1:end-1)';ginfo(5).XG(end:-1:1,:)'];
yg = [ginfo(1).YG(1:end-1,:);ginfo(2).YG(1:end-1,:); ...
      ginfo(4).YG(end:-1:1,1:end-1)';ginfo(5).YG(end:-1:1,:)'];
zA = [ginfo(1).rA;ginfo(2).rA;ginfo(4).rA(end:-1:1,:)'; ...
      ginfo(5).rA(end:-1:1,:)'];

% transform longitudes:
i1 = find(xc<0); xc(i1) = xc(i1) + 360;
i1 = find(xg<0); xg(i1) = xg(i1) + 360;
warning('remember the hack for xg');
xg(end,:) = 360; xg(1,:) = 0;

% do the rotation now for the side faces
xc = xc+s_lon;
xg = xg+s_lon;
yc = yc+s_lat;
yg = yg+s_lat;

disp( sprintf('# of recursive levels = %i',nfinepow) );

% Create fine grid for interpolation computations
nfine=2^nfinepow;

% first for the cap (remember the rotation)
xc3 = ginfo(3).XC+s_lon;
yc3 = ginfo(3).YC+s_lat;
xg3 = ginfo(3).XG+s_lon;
yg3 = ginfo(3).YG+s_lat;

clear xcf3 ycf3
xcf = zeros(size(xc3)*nfine);
ycf = zeros(size(yc3)*nfine);
xfine = [0.5:1:nfine]/nfine;
disp('generating fine grid')
% we need special treatment of the polar caps!
% $$$ xg1 = xg3; i1 = find(xg1<0); xg1(i1) = xg1(i1) + 360;
% $$$ i1 = find(xg1>=270); xg1(i1) = NaN;
% $$$ xg2 = xg3; i2 = find(xg2>=90); xg2(i2) = NaN;
for ix = 1:size(xc3,1)
  for iy = 1:size(yc3,2)
    xxg1 = xg3(ix:ix+1,iy:iy+1);
    if ~isempty(find(xxg1(:) >= 100)) & ~isempty(find(xxg1(:) < 100))
      i1 = find(xxg1(:) < 0); xxg1(i1) = xxg1(i1)+360;
    end
    xx = interp2([0 1],[0 1],xxg1',xfine,xfine')';
    xcf3((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = xx;
    yy = interp2([0 1],[0 1],yg3(ix:ix+1,iy:iy+1)',xfine,xfine')';
    ycf3((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = yy;
  end
end

xcf30{1}=xcf3;
for k=2:4
  xcf30{k}=xcf30{k-1};
  ii = find(xcf30{k}< (k-3)*90);
  xcf30{k}(ii) = xcf30{k}(ii) + 360;
end

xcf3 = [rot90(xcf30{1});xcf30{2};rot90(xcf30{3},3);rot90(xcf30{4},2)];
ycf3 = [rot90(ycf3);ycf3;rot90(rot90(rot90(ycf3)));rot90(rot90(ycf3))];
%i1 = find(xcf3<0); xcf3(i1) = xcf3(i1) + 360;


xcf = zeros(size(xc)*nfine);
ycf = zeros(size(yc)*nfine);
xfine = [0.5:1:nfine]/nfine;
for ix = 1:size(xc,1)
  for iy = 1:size(yc,2)
    xx = interp2([0 1],[0 1],xg(ix:ix+1,iy:iy+1)',xfine,xfine')';
    xcf((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = xx;
    yy = interp2([0 1],[0 1],yg(ix:ix+1,iy:iy+1)',xfine,xfine')';
    ycf((ix-1)*nfine+1:ix*nfine,(iy-1)*nfine+1:iy*nfine) = yy;
  end
end

xcf = [xcf xcf3];
ycf = [ycf ycf3];

save HFGRID.mat xcf ycf nfinepow

disp( sprintf('Fine grid: %ix%i',size(xcf,1),size(ycf,2) ));
%disp( sprintf('Fine grid res.: %gx%g',[min(dlonf(:)) min(dlatf(:))]) );

% now create the "cap" for the coarse grid
xc3 = ginfo(3).XC;
xc3 = [rot90(xc3);xc3;rot90(rot90(rot90(xc3)));rot90(rot90(xc3))];
yc3 = ginfo(3).YC;
yc3 = [rot90(yc3);yc3;rot90(rot90(rot90(yc3)));rot90(rot90(yc3))];
xg3 = ginfo(3).XG;
xg3 = [rot90(xg3(:,1:end-1));xg3(1:end-1,:);rot90(rot90(rot90(xg3(:,2:end))));rot90(rot90(xg3))];
yg3 = ginfo(3).YG;
yg3 = [rot90(yg3(:,1:end-1));yg3(1:end-1,:);rot90(rot90(rot90(yg3(:,2:end))));rot90(rot90(yg3))];
zA3 = ginfo(3).rA;
zA3 = [rot90(zA3);zA3;rot90(rot90(rot90(zA3)));rot90(rot90(zA3))];

i1 = find(xc3<0); xc3(i1) = xc3(i1) + 360;
i1 = find(xg3<0); xg3(i1) = xg3(i1) + 360;

% assemble
xc = [xc xc3+s_lon];
yc = [yc yc3+s_lat];
xg = [xg xg3+s_lon];
yg = [yg yg3+s_lat];
zA = [zA zA3];

nxc=size(xc,1);
disp([sprintf('nxc = %i factors(nxc)=',nxc) sprintf(' %i',factor(nxc))]);
nyc=size(yc,2);
disp([sprintf('nyc = %i factors(nyc)=',nyc) sprintf(' %i',factor(nyc))]);

lon_lo=round(min(xg(:)));
lon_L=round(max(xg(:)))-lon_lo;
lat_lo=round(min(yg(:)));
lat_L=round(max(yg(:)))-lat_lo;
lon_hi=lon_lo+lon_L;
lat_hi=lat_lo+lat_L;
  
xu=(xg(:,1:end-1)+xg(:,2:end))/2;
yu=(yg(:,1:end-1)+yg(:,2:end))/2;
xv=(xg(1:end-1,:)+xg(2:end,:))/2;
yv=(yg(1:end-1,:)+yg(2:end,:))/2;

disp( sprintf('[lon_lo lon hi lat_lo lat_hi] = %f %f %f %f',[lon_lo lon_hi lat_lo lat_hi]));
disp( sprintf('Periodic = %i',periodic) );

% Structure
GRID.xc=xc;
GRID.yc=yc;
GRID.xg=xg;
GRID.yg=yg;
GRID.rac=zA;

% Store data in HGRID.mat
nt=1;
save HGRID.mat nxc nyc nt xc yc xu yu xv yv xg yg zA ...
               lon_lo lon_hi lat_lo lat_hi periodic GRID