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

%  clear all ; close all
dbug = 0;

%  Get 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
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
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
        ginfo(k).bathy(i,j) = NaN;
      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;
  surf( ginfo(k).cor(:,:,1), ...
        ginfo(k).cor(:,:,2), ...
        ginfo(k).cor(:,:,3), ginfo(k).b_oce )
end
axis equal, view(2)