close all; clear all;

% addpath c:\Delft3D\w32\matlab\
% addpath c:\Delft3D\w32\matlab\toolbox\
% wlsettings

% boundary conditions
name     = {'base','oblique','breakw','river'} 
tstart   = 5400;
duration = 3600*3+tstart;
zs0      = [0 0];
Hm0      = [6 6 6 2];
Tp       = [15 15 15 8];
dir      = [270 250 250 250];
s        = [20 20 20 20];
gammajsp = [3.3 3.3 3.3 3.3];
nbin     = [1 5 5 5];
struc    = [0 0 1 0];
backb    = {'abs_2d','abs_2d','abs_2d','wall'};

%% bathymetry and topography
% for all cases we use same bathy and topo
x = load('SantaRosaBase\x.grd');
y = load('SantaRosaBase\y.grd');
z = load('SantaRosaBase\z.dep');

% figure; surf(x,y,z); colorbar; shading flat;

% select bathy-topo we want to use
indxmax = min(find(x(1,:)>2100));
indxmin = min(find(x(1,:)>1000));
indymin = min(find(y(:,1)>100));
indymax = min(find(y(:,1)>1900));
x = x(indymin:indymax,indxmin:indxmax);
y = y(indymin:indymax,indxmin:indxmax);
z = z(indymin:indymax,indxmin:indxmax);

% fix topo at lateral boundary
indzt1 = min(find(z(1,:)>zs0(1,1))):1:max(find(z(1,:)>zs0(1,2)));
indztend = min(find(z(end,:)>zs0(1,1))):1:max(find(z(1,:)>zs0(1,2)));
z(1:3,indzt1) =  max(z(1:3,indzt1),max(zs0)+2);
z(end-2:end,indztend) = max(z(end-2:end,indztend),max(zs0)+2);
    
% smooth offshore and bay-side bathy-topo
z(:,1) = z(1,1);
indxintp = min(find(x(1,:)>1800));
zbay = -5;
z(:,end) = zbay;
z(:,indxintp:end) = interp2([x(:,indxintp) x(:,end)],[y(:,indxintp) y(:,end)], [z(:,indxintp) z(:,end)], x(:,indxintp:end), y(:,indxintp:end));  

%% generate model
xb_setpref('grid_finalise', {'zmin',-20,'z0',zbay,'nlat',25,'noff',3,'nbay',20,'slope',1/70});

for i=1:4
    
    i
    
    % non erodible layer
    ne = ones(size(z))*1; 
    
    xb2dh = xb_generate_model( 'bathy',{'x',x,'y',y,'z',z,'ne',ne...
                               'xgrid',{'dxmin',5,'wl',5,'ppwl',16}, ...
                               'xrange',[30:60],...
                               'ygrid',{'dymin',40,'dymax',60,'area_size',0.7},...
                               'finalise',{'seaward_extend','landward_extend','lateral_extend'}}, ...
                               'waves',{'Hm0',Hm0(i),'Tp',Tp(i),'s',s(i),'gammajsp',gammajsp(i),'mainang',dir(i),'duration',duration},...
                               'wavegrid',{'nbins',nbin(i)},...
                               'tide',{'time',[0 duration],'front',[zs0(1) zs0(1)],'back',[zs0(2) zs0(2)]},...
                               'settings',{'front','abs_2d','back',backb{i},...
                                           'random',0,...
                                           'morfac',10,'mortstart',tstart,...'
                                           'wetslp',0.15,'dryslp',1,...
                                           'D50',0.0002,'D90',0.0003,...
                                           'struct',struc(i),...
                                           'C',55,...
                                           'tstart',0,'tstop',duration',...
                                           'gridform','delft3d',...
                                           'xyfile','xy.grd',...
                                           'outputformat','netcdf',...
                                           'tintg',120,'globalvars',{'zb','zs','u','v','H','Fx','Fy','ceqsg','ccg','Susg','Svsg','sedero','dzav'}, ...
                                           'tintm',duration-tstart,'meanvars',{'zs','u','v','H','ceqsg','ccg','Susg','Svsg'}},...
                               'write', true, ...
                               'path', name{i} ...
                               );0;

   

    % transform to Delft3D grid
    xb_peel(xb2dh);                
    xgr = xfile.data.value;
    ygr = yfile.data.value;
    dep = depfile.data.value;
    ne  = ne_layer.data.value;
    dims = size(xgr)-1;
    [xori,yori] = xb_get(xb2dh,'xori','yori');
    
    % make bay walls at lateral boundary
    indz1 = find(dep(1,:)==max(dep(1,:)));
    indzend = find(dep(end,:)==max(dep(end,:)));
    dep(1:3,indz1:end) =  max(dep(1:3,indz1:end),max(zs0)+2);
    dep(end-2:end,indzend:end) = max(dep(end-2:end,indzend:end),max(zs0)+2);
    
    % make offshore breakwater
    if i==3
        % specify breakwater
        xstart   = 1200;
        Wbreak   = 50;
        Lbreak   = 800;
        ystart   = mean(ygr(:,1));
        %
        indx1 = min(find(xgr(1,:)>xstart));
        indx2 = min(find(xgr(1,:)>xstart+Wbreak));
        indy1 = min(find(ygr(:,1)>ystart-Lbreak/2));
        indy2 = min(find(ygr(:,1)>ystart+Lbreak/2));
        dep(indy1:indy2,indx1:indx2) = 2;
        ne(indy1:indy2,indx1:indx2) = 0;
    end
    
    
    xgri = xgr+xori;
    ygri = ygr+yori;

    ok1=wlgrid('write','xy.grd',xgri',ygri');
    
    % add -999 row and column for quickin
    depw = [dep -999*ones(size(dep,1),1); ...
                -999*ones(1,size(dep,2)+1)];
    new  = [ne  -999*ones(size(dep,1),1); ...
                -999*ones(1,size(dep,2)+1)];

    % and fix params
    par = xb_read_params([name{i},'\params.txt'])
    par = xb_del(par,'xfile','yfile','posdown','swtable','xori','yori');
    
    % river input
    if i==4
        Wriver = 100;
        ystart   = mean(ygri(:,1));
        par = xb_set(par,'ndischarge',1,'ntdischarge',2,'disch_loc_file','discharge_locations.txt','disch_timeseries_file','discharge_timeseries.txt');
        disloc = [xgri(1,end) ystart-Wriver/2 xgri(1,end) ystart+Wriver/2];
        dist = [0 450; duration 450;];
        save('discharge_locations.txt','-ascii','disloc');
        save('discharge_timeseries.txt','-ascii','dist');
        movefile('discharge_locations.txt',name{i});
        movefile('discharge_timeseries.txt',name{i});
        indy1 = min(find(ygri(:,1)>ystart-Wriver/2));
        indy2 = min(find(ygri(:,1)>ystart+Wriver/2));
        depw(:,indy1:indy2) = min(depw(:,indy1:indy2),-3);
    
    end
    xb_write_params([name{i},'\params.txt'],par);
    
    ok2=wldep('write','bed.dep',depw');
    ok3=wldep('write','nebed.dep',new');
    movefile('xy.grd',name{i});
    movefile('bed.dep',name{i});
    movefile('nebed.dep',name{i});
    deletefile([name{i},'\x.grd']);
    deletefile([name{i},'\y.grd']);
    deletefile([name{i},'\swtable.txt']);

end
%% check what we've got

figure; surf(xgr,ygr,dep); colorbar; axis equal;

% indp = size(x,1);
% figure; plot(x(indp,:),z(indp,:),'r'); hold on; 
% plot([x(indp,1) x(indp,end)],[zs0(1) zs0(1)],'b'); plot(x(indp,:),z(indp,:),'k');
% plot(xgr(1,:),dep(1,:),'g');