close all; clear all;

% goal is to set-up fast 2DH overwash simulattion for XBeach course

% location xb_testbed Santa Rosa island
loc = 'D:\subversion_xbeach\testbed\input\2DH_Santa_Rosa_Island\';

%% define models to be setup
name     = {'inundation_base','overwash_base','collision_base','inundation_broad','inundation_narrow','inundation_nourishment_t=0','inundation_nourishment_t=30','inundation_harbor_t=0','inundation_harbor_t=1'};
s        = [20 20 20 2 100 20 20 20 20];                         % coefficient for directional spreading
gammajsp = [3.3 3.3 3.3 1 10 3.3 3.3 3.3 3.3];                   % peak enhancement factor
bedconf  = [1 1 1 1 1 3 3 2 2];                                  % base, port, nourishment
tr       = [15 9 3 15 15 15 15 15 15]*3600;   % get surge height, wave height and wave period at tr hours from original surge hydrograph 
option   = [3 3 3 3 3 3 3 3 3];

%% boundary conditions
wt = xb_read_waves([loc,'jonswaplist.txt']);
wt = xb_peel(wt);
wt.mainang = 270;

surge = xb_read_tide([loc,'surge.txt']);
surge = xb_peel(surge);

% compute zs0 for each sallenger regime
for i = 1:length(name)
    indwt(i) = min(find(tr(i)<cumsum(wt.duration)));
    zs0(i,1) = interp1(surge.time,surge.tide(:,1),tr(i));
    zs0(i,2) = interp1(surge.time,surge.tide(:,2),max(tr(i)-7200,1)); %2 hour delay off surge in bay; idem shape
end

%% bathy
% OPTION I: from XB testbed
% start with Santa Rosa island:
x = load([loc,'x.grd']);
y = load([loc,'y.grd']);
z = load([loc,'z.dep']);

% figure; pcolor(xt,yt,zt); shading flat;% 
% figure; plot(surge.time,surge.tide(:,1),'b'); hold on; plot(surge.time,surge.tide(:,2),'r');
% figure; subplot(211); plot(cumsum(wt.duration),wt.Hm0,'r'); subplot(212); plot(cumsum(wt.duration),wt.Tp,'r');

% 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));
xt = x(indymin:indymax,indxmin:indxmax);
yt = y(indymin:indymax,indxmin:indxmax);
zt = z(indymin:indymax,indxmin:indxmax);

% add three columns:
ztt = zeros(size(zt,1)+4,size(zt,2));
ztt(3:end-2,:) = zt;
ztt(1:2,:) = repmat(zt(1,:),2,1);
ztt(end-1:end,:) = repmat(zt(end,:),2,1);

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

% OPTION II: read grid and bathy from Delft3D model
% own model
% [xt2,yt2] = wlgrid('read',[name{3},'\xy.grd']);
% zt2 = wldep('read',[name{3},'\bed.dep'],size(xt2));

% OPTION III: DJ's delft3d models
modnam = {'base_case','harbor','nourishment'};
for i = 1:3    
    [xt3,yt3] = wlgrid('read',['D:\documents\XBeach\CoastalSediments\USB\Model_Examples\Delft3D\',modnam{i},'\flow.grd']);
    zt3a = wldep('read',['D:\documents\XBeach\CoastalSediments\USB\Model_Examples\Delft3D\',modnam{i},'\flow_',modnam{i},'.dep'],size(xt3)+1);
    
    if i<0
        % fix bay
        inddj = min(find(zt3a(:,1)==-10));
        xt4 = xt3(1:inddj-1,:);
        yt4 = yt3(1:inddj-1,:);
        zt4a = -zt3a(1:inddj-1,1:end-1);
        zt4a(zt4a==10)=1.75;

        % smooth transition from topo to bay bathy
        dxt = 217.5;
        zt4a(end+1,:) = zbay;
        xt4(end+1,:) = xt4(end,:)+dxt;
        yt4(end+1,:) = yt4(end,:);
    else % DJ included bay in model
        xt4 = xt3;
        yt4 = yt3;
        zt4a = -zt3a(1:end-1,1:end-1);
    end
    
    xDJ{i} = xt4';
    yDJ{i} = yt4';
    zDJ{i} = zt4a';
    
    if i==2 %get hard layer for harbor
        zt3b = wldep('read',['D:\documents\XBeach\CoastalSediments\USB\Model_Examples\Delft3D\',modnam{i},'\underlayer_harbor_03.dep'],size(xt3)+1);
        % zt4b = zt3b(1:inddj-1,1:end-1);
        % zt4b(end+1,:) = 10;
        % zt4b = zt4b';
        zt3b = zt3b';
    end
end
% figure; 
% figure; 
% subplot(121); 
% surf(xt4,yt4,zt4a); caxis([-20 10]);
% subplot(122);
% surf(xt4,yt4,zt4b); caxis([0 15]);

%% create model
xb_setpref('grid_finalise', {'zmin',min(min(zt4a)),'z0',zbay,'nlat',18,'noff',3,'nbay',20});

for i = [1 4 5 6 7] %1:length(name)
    
    i
    
    tstart = 5400;
    duration = 3600*3+tstart;
    
    % make XBeach model
    if option(i)==1
    
        xb2dh = xb_generate_model( 'bathy',{'x',xt,'y',yt,'z',zt,...
                                   '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',wt.Hm0(indwt(i)),'Tp',wt.Tp(indwt(i)),'s',s(i),'gammajsp',gammajsp(i),'duration',duration},...
                                   'wavegrid',{'nbins',1},...
                                   'tide',{'time',[0 duration],'front',[zs0(i,1) zs0(i,1)],'back',[zs0(i,2) zs0(i,2)]},...
                                   'settings',{'front','abs_2d','back','abs_2d',...
                                               'morfac',10,'mortstart',tstart,...'
                                               'wetslp',0.15,'dryslp',1,...
                                               'D50',0.0002,'D90',0.0003,...
                                               'C',55,...
                                               'tstart',tstart,'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} ...
                                   );
                               
    elseif option(i)==2 %Delft3D model
        
        xb2dh = xb_generate_model( 'bathy',{'x',xt2','y',yt2','z',zt2',...
                                   'xgrid',{'dxmin',5,'wl',5,'ppwl',16}, ...
                                   'xrange',[30:60],...
                                   'ygrid',{'dymin',40,'dymax',60,'area_size',0.7},...
                                   'finalise',{}},...
                                   'waves',{'Hm0',wt.Hm0(indwt(i)),'Tp',wt.Tp(indwt(i)),'s',s(i),'gammajsp',gammajsp(i),'duration',duration},...
                                   'wavegrid',{'nbins',1},...
                                   'tide',{'time',[0 duration],'front',[zs0(i,1) zs0(i,1)],'back',[zs0(i,2) zs0(i,2)]},...
                                   'settings',{'front','abs_2d','back','abs_2d',...
                                               'morfac',10,'mortstart',tstart,... 
                                               'tstart',tstart,'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} ...
                                   );          
                               
   elseif option(i)==3 %Delft3D model DJ
        
        % npoints = {'0 1800'; '2000 1800'; '3000 1800'; '3500 1800'; '4000 1800'; '4500 1800'; '5000 1800'; '5500 1800';};
%         npoints(1,:) = npoints(1,:)+xori;
%         npoints(2,:) = npoints(2,:)+yori;
        
       
        settings = {'front','abs_2d','back','abs_2d','epsi',0.005,'random',0,...
                    'morfac',10,'morstart',tstart,... 
                    'tstart',0,'tstop',duration,...
                    'gridform','delft3d',...
                    'xyfile','xy.grd',...
                    'outputformat','netcdf',...
                    'tintg',300,'globalvars',{'zb','zs','u','ue','v','ve','H','Fx','Fy','ceqsg','ccg','Susg','Svsg','sedero','dzav'}, ...
                    'tintm',duration-tstart,'meanvars',{'zs','u','ue','v','ve','H'}};
        
        bathy = {'x',xDJ{bedconf(i)},'y',yDJ{bedconf(i)},'z',zDJ{bedconf(i)},...
                 'rotate',false,...
                 'xgrid',{'dxmin',10,'wl',5,'ppwl',16}, ...
                 'xrange',[30:60],...
                 'ygrid',{'dymin',40,'dymax',60,'area_size',0.7},...
                 'finalise',{'seaward_extend','landward_extend'}};

        if bedconf(i)==2
            bathy = {bathy{:},'ne',zt4b};
            settings = {settings{:},'struct',1};
        end   
        
        xb2dh = xb_generate_model( 'bathy',bathy,...
                                   'waves',{'Hm0',wt.Hm0(indwt(i)),'Tp',wt.Tp(indwt(i)),'s',s(i),'gammajsp',gammajsp(i),'duration',duration},...
                                   'wavegrid',{'nbins',1},...
                                   'tide',{'time',[0 duration],'front',[zs0(i,1) zs0(i,1)],'back',[zs0(i,2) zs0(i,2)]},...
                                   'settings',settings,...
                                   'write', true, ...
                                   'path', ['USB\Model_Examples\XBeach\',name{i}] ...
                                   );          
                                                                                                                     
    end
    
    % Translate from xbeach grid to delft3d grid
    if i==1
        xb_peel(xb2dh);                
        xgr = xfile.data.value;
        ygr = yfile.data.value;
        dep = depfile.data.value;
        dims = size(xgr)-1;
        [xori,yori] = xb_get(xb2dh,'xori','yori');
        
        npxt = [0 2000 3000 3500 4000 4500 5000 5500];
        for j = 1:length(npxt)
            indp(j) = find(abs(xgr(45,:)-npxt(j))==min(abs(xgr(45,:)-npxt(j))));
            npoints{j,1} = [num2str(xgr(45,indp(j))+xori),'       ',num2str(ygr(45,indp(j))+yori)];
        end
     
    else
        xb_peel(xb2dh);
        xgrt = xfile.data.value;
        ygrt = yfile.data.value;
        dept = depfile.data.value;
        dep = interp2(xgrt,ygrt,dept,xgr,ygr);
    end
    
    % fix xori and yori
    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)];
            
    ok2=wldep('writeraw','bed.dep',depw');
    movefile('xy.grd',['USB\Model_Examples\XBeach\',name{i}]);
    movefile('bed.dep',['USB\Model_Examples\XBeach\',name{i}]);
    deletefile(['USB\Model_Examples\XBeach\',name{i},'\x.grd']);
    deletefile(['USB\Model_Examples\XBeach\',name{i},'\y.grd']);
    deletefile(['USB\Model_Examples\XBeach\',name{i},'\swtable.txt']);
    
    pt = xb_read_params(['USB\Model_Examples\XBeach\',name{i},'\params.txt'])
    pt = xb_set(pt,'nx',dims(2),'ny',dims(1),'tintp',30,'pointvars',{'zs','u','ue','v','ve','H'},'npoints',npoints);
    pt = xb_del(pt,'xfile','yfile','posdown','swtable','xori','yori');
    xb_write_params(['USB\Model_Examples\XBeach\',name{i},'\params.txt'],pt);
    
    clear xb2dh
end                          

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

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

ind1 = 1; ind2 = 50;
figure; plot(xgr(ind1:ind2,:),dep(ind1:ind2,:),'b-o');