%
%  Ed Hill
%
%  Generate approximate bathymetry for the llc grid
%

%  clear all ; close all

do_print = 0;
dbug = 0;

%  Get the ETOPO2 data
disp('Reading the ETOPO2 data...');
nlat =  5400;
nlon = 10800;
lons = linspace(-180, (180 - 2/60), nlon);
lats = linspace(90, -(90 - 2/60), nlat);
fid = fopen('ETOPO2.raw.bin', 'r', 'ieee-be');
et2 =  reshape( fread(fid, nlat*nlon, 'int16'), [ nlon nlat ]);
fid = fclose(fid);

if dbug > 10
  surf( et2(1:500,1:500) ), view(2), shading flat
  ilon = [  1:10:5000 ];
  ilat = [ 1:10:5400 ];
  surf(lons(ilon), lats(ilat), et2(ilon,ilat)'), view(2), ...
      shading flat, axis equal
end

%  Get the llc grid information
disp('Reading the grid information...');
fnc_in = 'llc_p90_%d.nc';
vnms = { 'XC' 'YC' 'XG' 'YG' };
d2r = pi/180.0;
ginfo = {};
for k = 1:5
  nc = netcdf(sprintf(fnc_in,k),'nowrite');
  for j = 1:length(vnms)
    eval(sprintf('ginfo(k).%s = nc{''%s''}(:);',vnms{j},vnms{j}));
  end
  nc = close(nc);
  
  %  compute 3D coords of the C and G points
  cor = zeros( [ size(ginfo(k).XG) 3 ]);
  [cor(:,:,1),cor(:,:,2),cor(:,:,3)] = sph2cart( ginfo(k).XG*d2r, ...
                                                 ginfo(k).YG*d2r,1 );
  cen = zeros( [ size(ginfo(k).XC) 3 ]);
  [cen(:,:,1),cen(:,:,2),cen(:,:,3)] = sph2cart( ginfo(k).XC*d2r, ...
                                                 ginfo(k).YC*d2r,1 );
  ginfo(k).cor = cor;
  ginfo(k).cen = cen;
  
  %  empty bathy
  ginfo(k).bathy = zeros(size( size(ginfo(k).XC) ));
end

%  Do a quick-and-dirty regrid
disp('Performing the regrid...');
s_lon = -37;
s_lat =   0.0;
for k = 1:5
  disp(sprintf('  k = %d',k));
  %
  %  /  i,j+1 i+1,j+1 \
  %  \  i,j   i+i,j   /
  for j = 1:size(ginfo(k).XC,2)
    jr = [ j j+1 ];
    for i = 1:size(ginfo(k).XC,1)
      ir = [ i i+1 ];
      lon_min = min(min( ginfo(k).XG(ir,jr) + s_lon ));
      lon_max = max(max( ginfo(k).XG(ir,jr) + s_lon ));
      lat_min = min(min( ginfo(k).YG(ir,jr) + s_lat ));
      lat_max = max(max( ginfo(k).YG(ir,jr) + s_lat ));
      % [ lon_min lon_max lat_min lat_max ]
      lon_mm = [ mod(lon_min + 180,360)-180 ...
                 mod(lon_max + 180,360)-180 ];
      lon_min = min(lon_mm);
      lon_max = max(lon_mm);

      d_lon = lon_max - lon_min;
      if (abs(d_lon) > 45 && lat_min < 88)
        i_lon = find( lons >= lon_max  ...
                      |  lons <= lon_min );
      else
        i_lon = find(lon_min <= lons  &  lons <= lon_max);
      end
      i_lat = find(lat_min <= lats & lats <= lat_max);
      
      %  average the bathy
      bv = et2(i_lon,i_lat);
      if length(bv(:)) > 1
        ginfo(k).bathy(i,j) = mean(bv(:));
      else
        %  If not averaging, then interpolate
        %  ginfo(k).bathy(i,j) = NaN;
        i_lon = interp1(lons,[1:length(lons)],lon_min,'nearest');
        i_lat = interp1(lats,[1:length(lats)],lat_min,'nearest');
        ginfo(k).bathy(i,j) = et2(i_lon,i_lat);
      end
    end
  end
end

%  Plot the resulting bathymetry
k = 3;
for k = 1:5
  disp(sprintf('  k = %d',k));
  ginfo(k).b_oce = ginfo(k).bathy;
  %  ginfo(k).b_oce(find(ginfo(k).b_oce >= 0.0)) = NaN;
  ginfo(k).b_oce(find(ginfo(k).b_oce >= 0.0)) = 1000;
  if k == 2
    hold on
  end
  %  figure(k)
  surf( ginfo(k).cor(:,:,1), ...
        ginfo(k).cor(:,:,2), ...
        ginfo(k).cor(:,:,3), ginfo(k).b_oce )
  if k == 5
    hold off
  end
end
axis equal, view(2)


%  Plot the bathymetry at 1/2 resolution
k = 3;
for k = 1:5
  disp(sprintf('  k = %d',k));
  ginfo(k).c_half = ginfo(k).cor(1:2:end,1:2:end,:);
  ginfo(k).b_half = ginfo(k).bathy(1:2:end,1:2:end);
  ginfo(k).b_half(find(ginfo(k).b_half >= 0.0)) = 1000;
  if k == 2
    hold on
  end
  %  figure(k)
  surf( ginfo(k).c_half(:,:,1), ...
        ginfo(k).c_half(:,:,2), ...
        ginfo(k).c_half(:,:,3), ginfo(k).b_half )
  if k == 5
    hold off
  end
end
axis equal, view(2)


if do_print > 0
  print -depsc llc_bathy.eps
end


%  Write the bathymetry to netCDF compatible with MNC
id = 'bathy';
units = 'm';
for k = 1:5
  disp(sprintf('  k = %d',k));
  ginfo(k).b_out = ginfo(k).bathy;
  ginfo(k).b_out(find(ginfo(k).bathy >= 0.0)) = 0.0;
  bfile = sprintf('llc_p90_bathy.f%03d.nc',k);
  
  nc = netcdf(bfile, 'clobber');
  nc.reference = 'bathymetry';
  nc.author = 'Ed Hill <eh3@mit.edu>';
  nc.date = eval('date');
  nc('X') = size( ginfo(k).XC, 1 );
  nc('Y') = size( ginfo(k).XC, 2 );
  nc{ id } = { 'Y', 'X' };
  nc{ id }.units = units;
  nc{ id }(:) = ginfo(k).b_out';
  nc = close(nc);
end

%  Write the bathymetry to a set of 90-by-90 per-tile MDSIO-compatible
%  files
%       face  is  js
iti = [   1    1   1  ; ...
          1   91   1  ; ...
          1  181   1  ; ...
          1  271   1  ; ...
          2    1   1  ; ...
          2   91   1  ; ...
          2  181   1  ; ...
          2  271   1  ; ...
          3    1   1  ; ...
          4    1   1  ; ...
          4    1  91  ; ...
          4    1 181  ; ...
          4    1 271  ; ...
          5    1   1  ; ...
          5    1  91  ; ...
          5    1 181  ; ...
          5    1 271  ];
for k = 1:size(iti,1)
  disp(sprintf('  k = %d',k));
  fname = sprintf('bathy.%03d.001.data',k);
  fid = fopen(fname, 'w', 'ieee-be');
  %  be careful of index orientation
  is = iti(k,3);  ie = is + 90 - 1;
  js = iti(k,2);  je = js + 90 - 1;
  tmp = ginfo(iti(k,1)).b_out(is:ie,js:je)';
  surf(tmp), view(2), axis equal, shading flat, pause(1)
  fwrite(fid,tmp,'real*8',0,'ieee-be');
  fclose(fid);
end