lon_east = -15.00 ;
lon_west = 10.00 ;
lat_south = 43.00 ;
lat_north = 60.00 ;

dx = 0.0167 ;
dy = 0.0167 ;

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 ;

fname='GASC-1M';

set(groot,'DefaultFigureColormap',jet)
bin_dir = '/home1/datahome/cpeureux/tools/SKIM/gridgen/bin';
ref_dir = '/home1/datahome/cpeureux/tools/SKIM/gridgen/reference';
nml_dir = '/home1/datahome/cpeureux/tools/SKIM/gridgen/namelist';
data_dir = '/home1/datahome/cpeureux/tools/SKIM/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 = 10; % 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;

LIM_VAL = 0.5;
OFFSET = max([dx dy]);
m2 = clean_mask(lon,lat,m1,b,LIM_VAL,OFFSET);

LAKE_TOL = 100;
IS_GLOBAL = 0;
[m4,mask_map] = remove_lake(m2,LAKE_TOL,IS_GLOBAL);

figure(1);clf;
pcolor(lon,lat,m4);shading flat;caxis([0 3]);colorbar

OBSTR_OFFSET = 1;
[sx1,sy1] = create_obstr(lon,lat,b,m2,OBSTR_OFFSET,OBSTR_OFFSET);

depth_scale = 1000;
obstr_scale = 100;
d = round((depth)*depth_scale);
write_ww3file([data_dir,'/',fname,'.bot'],d);
write_ww3file([data_dir,'/',fname,'.mask_nobound'],m4);
d1 = round((sx1)*obstr_scale);
d2 = round((sy1)*obstr_scale);
write_ww3obstr([data_dir,'/',fname,'.obst'],d1,d2);
write_ww3meta ([data_dir,'/',fname],[nml_dir,'/gridgen.',fname,'.nml'],'RECT',lon,lat,1/depth_scale,1/obstr_scale,1.0);
