lon_east = -153 ;
lon_west = -117 ;
lat_south = 24 ;
lat_north = 48 ;

dx = 0.125 ;
dy = 0.125 ;

lon1d = (lon_east:dx:lon_west) ;
lat1d = (lat_south:dx:lat_north) ;

[lon,lat] =  meshgrid(lon1d,lat1d) ;

CUT_OFF = 0.0 ;
LIM_BATHYMETRY = 0.4 ;
DRY_VAL = 999999 ;

set(groot,'DefaultFigureColormap',jet)
bin_dir = '/mnt/datahome/tutos/TUTORIAL_GRIDGEN/bin';
ref_dir = '/mnt/datahome/tutos/TUTORIAL_GRIDGEN/reference';
nml_dir = '/mnt/datahome/tutos/TUTORIAL_GRIDGEN/namelist';
data_dir = '/mnt/datahome/tutos/TUTORIAL_GRIDGEN/data';

addpath(bin_dir,'-END');

xvar = 'lon' ;
yvar = 'lat' ;
zvar = 'elevation' ;

depth = generate_grid('rect',lon,lat,ref_dir,'gebco',LIM_BATHYMETRY,CUT_OFF,DRY_VAL,xvar,yvar,zvar) ;

% plot bathy
figure(1);clf;
d=depth;d(d==DRY_VAL)=nan; pcolor(lon,lat,d); shading flat; colorbar

%land-sea mask
m1 = ones(size(depth));
m1(depth == DRY_VAL) = 0;

% coastal boundaries
lon_start = min(min(lon))-dx+360;
lon_end = max(max(lon))+dx+360;
lat_start = min(min(lat))-dy;
lat_end = max(max(lat))+dy;
coord = [lat_start lon_start lat_end lon_end];
MIN_DIST = 4; % minimum distance between edge of polygon and boundary

load([ref_dir,'/coastal_bound_full.mat']);
[b,N] = compute_boundary(coord,bound);

% plot coastal boundaries
figure(3);clf;
for i = 1:N
plot(b(i).x,b(i).y);
hold on;
end

SPLIT_LIM = 0.5;
b_split = split_boundary(b,SPLIT_LIM,MIN_DIST);

Nb = length(b_split);
figure(1);clf;
for i = 1:Nb
plot(b_split(i).x,b_split(i).y);
hold on;
end;