function varargout = desb_calculate_profile_volumes(varargin)
%DESB_CALCULATE_PROFILE_VOLUMES  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = desb_calculate_profile_volumes(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   desb_calculate_profile_volumes
%
%   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: 07 Sep 2011
% Created with Matlab version: 7.12.0.635 (R2011a)

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

%%
%%
OPT = struct(...
    'id', 7003925,...
    'irun', 0,...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'jarkus_url', fullfile(fileparts(which(mfilename)), '..', 'data', 'jarkus.nc'));

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

%% read data
tr = jarkus_transects(...
    'id', OPT.id,...
    'url', OPT.jarkus_url);

nT = size(tr.altitude,1);
ntransects = size(tr.altitude,2);

[xmin xtransition xmax zmin] = deal(NaN(size(OPT.id)));
csb = NaN(ntransects,5);
[Vs Vw] = deal(NaN(nT, ntransects, size(csb,2)-1));

for itransect = 1:ntransects
    basepath = fullfile(OPT.maindir, sprintf('%08i', OPT.id(itransect)), 'durosplus', sprintf('run%02i', OPT.irun));
    ncfile = fullfile(basepath, sprintf('%08i_pp.nc', OPT.id(itransect)));
    xmax = max(nc_varget(ncfile, 'x_max'));
    xtransition = mean(nc_varget(ncfile, 'x_transition'));
    xmin = min(nc_varget(ncfile, 'x_min'));
    nnid = ~isnan(sum(squeeze(tr.altitude(:,itransect,:))));
    x = tr.cross_shore;
    csb(itransect,:) = [min(x(nnid)) xmin xtransition xmax max(x(nnid))];
    zmin(itransect) = min(min(squeeze(tr.altitude(:,itransect,:))));
    zmax(itransect) = max(max(squeeze(tr.altitude(:,itransect,:))));
    for iT = 1:nT
        for izone = 1:size(csb,2)-1
            [Vs(iT,itransect,izone) resv] = jarkus_getVolume(tr.cross_shore, squeeze(tr.altitude(iT,itransect,:)),...
                'LowerBoundary', zmin(itransect),...
                'LandwardBoundary', csb(itransect,izone),...
                'SeawardBoundary', csb(itransect,izone+1));
            [Vw(iT,itransect,izone) resw] = jarkus_getVolume(minmax(tr.cross_shore), zmax(itransect)*ones(1,2),...
                'x2', tr.cross_shore,...
                'z2', squeeze(tr.altitude(iT,itransect,:)),...
                'LandwardBoundary', csb(itransect,izone),...
                'SeawardBoundary', csb(itransect,izone+1));
        end
    end
    
end

%% write data to netcdf file
outputfile = fullfile(fileparts(which(mfilename)), '..', 'data', 'jarkus_volumes.nc');

nc_create_empty(outputfile);

nc_add_dimension(outputfile, 'alongshore', ntransects)
nc_add_dimension(outputfile, 'time', nT)
nc_add_dimension(outputfile, 'zone', size(Vs,3))
nc_add_dimension(outputfile, 'cross_shore_boundary', size(Vs,3)+1)

% read info from jarkus file
info = nc_info(OPT.jarkus_url);

% add variable and its attributes
varname = 'id';
varid = ismember({info.Dataset.Name}, varname);
s.Name = varname;
s.Nctype    = info.Dataset(varid).Nctype;
s.Dimension = info.Dataset(varid).Dimension;
s.Attribute = rmfield(info.Dataset(varid).Attribute, {'Nctype', 'Datatype'});
nc_addvar(outputfile, s);

s.Name      = 'time';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, ...
    'Value', {'time', 'days since 1970-01-01 00:00 +1:00', 'T'});
nc_addvar(outputfile, s);

s.Name      = 'zone';
s.Nctype    = nc_int;
s.Dimension = {'zone'};
s.Attribute = struct('Name', {'long_name', 'comments'}, ...
    'Value', {'cross-shore zone identifier', 'indication of cross-shore zone: 1, landward side; 2, middle part; 3, seaward side'});
nc_addvar(outputfile, s);

s.Name      = 'cross_shore_boundary';
s.Nctype    = nc_double;
s.Dimension = {'alongshore' 'cross_shore_boundary'};
s.Attribute = struct('Name', {'long_name'}, ...
    'Value', {'cross-shore boundary for volume calculation'});
nc_addvar(outputfile, s);

s.Name      = 'upper_boundary';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name' 'units'}, ...
    'Value', {'upper boundary for water/air volume calculation' 'm'});
nc_addvar(outputfile, s);

s.Name      = 'lower_boundary';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name' 'units'}, ...
    'Value', {'lower boundary for sediment volume calculation' 'm'});
nc_addvar(outputfile, s);

s.Name      = 'volume_sediment';
s.Nctype    = nc_double;
s.Dimension = {'time' 'alongshore' 'zone'};
s.Attribute = struct('Name', {'long_name' 'units'}, ...
    'Value', {'profile volume between cross-shore zone boundaries and above lower boundary' 'm^3/m'});
nc_addvar(outputfile, s);

s.Name      = 'volume_rest';
s.Nctype    = nc_double;
s.Dimension = {'time' 'alongshore' 'zone'};
s.Attribute = struct('Name', {'long_name' 'units'}, ...
    'Value', {'volume above profile between cross-shore zone boundaries and below upper boundary' 'm^3/m'});
nc_addvar(outputfile, s);

% put values into variables
nc_varput(outputfile, 'id', tr.id)
nc_varput(outputfile, 'time',  tr.time)
nc_varput(outputfile, 'zone',  1:size(Vs,3))
nc_varput(outputfile, 'cross_shore_boundary', csb)
nc_varput(outputfile, 'upper_boundary', zmax)
nc_varput(outputfile, 'lower_boundary', zmin)
nc_varput(outputfile, 'volume_sediment', Vs)
nc_varput(outputfile, 'volume_rest', Vw)

history = sprintf('created on %s', datestr(now));
nc_attput(outputfile, nc_global, 'history', history);

producer = sprintf('$HeadURL$ $Revision$');
nc_attput(outputfile, nc_global, 'producer', producer);

description = sprintf('Profile volumes in various cross-shore zones');
nc_attput(outputfile, nc_global, 'description', description);