function varargout = desb_prepare_xbeach(varargin)
%DESB_PREPARE_XBEACH  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = desb_prepare_xbeach(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   desb_prepare_xbeach
%
%   See also

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2011 Delft University of Technology
%       Kees den Heijer
%
%       C.denHeijer@TUDelft.nl
%
%       Faculty of Civil Engineering and Geosciences
%       P.O. Box 5048
%       2600 GA Delft
%       The Netherlands
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://www.OpenEarth.eu">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute
% your own tools.

%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 08 Sep 2011
% Created with Matlab version: 7.7.0.471 (R2008b)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

%%
OPT = struct(...
    'id', 7003925,...
    'irun', 2,...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'WL', 5,...
    'Hsig', 9,...
    'Tp', 12,...
    'D50', 225e-6,...
    'duration', 5*3600,...
    'remote', false);

OPT = setproperty(OPT, varargin{:});

%% read data
tr = jarkus_transects(...
    'id', OPT.id,...
    'url', fullfile(fileparts(which(mfilename)), '..', 'data', 'jarkus.nc'));

%% derive relevant info
T = year(tr.time+datenum(1970,1,1));
if OPT.irun == 1
    T = NaN;
end

if OPT.irun > 1
    ncfile = fullfile(fileparts(which(mfilename)), '..', 'data', 'jarkus_volumes.nc');
    id = nc_varget(ncfile, 'id');
    csb = nc_varget(ncfile, 'cross_shore_boundary');
    xswitch = csb(ismember(id, OPT.id), 3);
end

%% prepare XBeach in loop
ntransects = size(tr.altitude,2);
nT = length(T);
for itransect = 1:ntransects
    basepath = fullfile(OPT.maindir, sprintf('%08i', OPT.id(itransect)), 'xbeach', sprintf('run%02i', OPT.irun));
    
    % create maximum transect as basis for the grid
    nnid = ~isnan(sum(squeeze(tr.altitude(:,itransect,:))));
    x = tr.cross_shore(nnid)';
    z = max(squeeze(tr.altitude(:,itransect,nnid)))';
    
    [x ix] = sort(-x);
    z = z(ix);
    
    morfac = 5;
    zsea = min(floor([min(OPT.WL(:)) - max(OPT.Tp(:))/.55 -20]));
    xb_setpref('grid_finalise', {'zmin', zsea})
    xb = xb_generate_model(...
        'bathy', {'x', x, 'z', z, 'xgrid', {'dxmin', 2}, 'finalise', {'seaward_extend'}, 'superfast', true},...
        'tide', {'front' OPT.WL},...
        'waves', {...
        'type', 'jonswap',...
        'Hm0', OPT.Hsig,...
        'Tp', OPT.Tp,...
        'duration', 1800*morfac},...
        'wavegrid', {'nbins' 1},...
        'settings', {...
        'morfac', morfac,...
        'posdwn', -1,...
        'tintg', 100,...
        'tstop', ceil(OPT.duration/morfac)*morfac,...
        'D50', OPT.D50,...
        'D90', OPT.D50*1.5,...
        'random' 0,...
        'outputformat', 'netcdf',...
        'meanvars', {'H' 'zs'},...
        });
    % remove swtable from structure
    xb = xb_del(xb, 'swtable');
    
    xori = xb_get(xb, 'xori');
    xfile = xb_get(xb, 'xfile.xfile');
    depfile = xb_get(xb, 'depfile.depfile');
    
    
    for iT = 1:nT
        
        
        if OPT.irun ~= 1
            % actual transect
            nnid = ~isnan(squeeze(tr.altitude(iT,itransect,:)));
            z = squeeze(tr.altitude(iT,itransect,nnid));
            casename = sprintf('%i', T(iT));
        end
        x = tr.cross_shore(nnid)';
        if OPT.irun > 0
            % mean transect
            nnid = ~isnan(sum(squeeze(tr.altitude(:,itransect,:))));
            zm = mean(squeeze(tr.altitude(:,itransect,nnid)))';
            x = tr.cross_shore(nnid)';
            if OPT.irun == 1
                % use mean transect
                z = zm;
                casename = 'mean';
            elseif OPT.irun == 2
                % use actual seaward tail and mean for the rest of the
                % transect
                z(x<xswitch(itransect)) = zm(x<xswitch(itransect));
            elseif OPT.irun == 3
                % use mean seaward tail and actual data for the rest of the
                % profile
                z(x>xswitch(itransect)) = zm(x>xswitch(itransect));
            end
        end
        
        depfilea = interp1(-xori-x, z, xfile);
        depfilea(isnan(depfilea)) = depfile(isnan(depfilea));
        
        xb = xb_set(xb, 'depfile.depfile', depfilea);
        
        calcdir = fullfile(basepath, casename);
        mkpath(calcdir)
        xb_write_input(calcdir, xb);
        if OPT.remote
            fid = fopen(fullfile(calcdir, sprintf('%s.sh', casename)), 'w');
            fprintf(fid, '#!/bin/sh\n. /opt/sge/InitSGE\n. ~/.bashrc\n/u/heijer/bin/xbeach_netcdf\n');
            fclose(fid);
        end
    end
end