%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DWI2nc
%
% Generates DWI.nc which contains the following coastal indicators: 
%  volumes (in m2) of different depth zones: 
%   * -1/-5  m
%   * -1/-8  m
%   * -1/-12 m
%   * -5/-8  m
%   * -5/-12 m
%   * -8/-12 m
% For each zone parameters are calculated in two diffent ways:
% 1) regarding the volume under the profile between the upper and lower
% boundaries of the depthzone
% 2) regarding the volume under the profile between the landward and
% seaward boundaries from the reference year relative to the lower boundary
% of the depthzone
%
% A more elaborate (technical) description is found in DWI2nc.docx. 
%
% Tommer Vermaas, oktober 2011
% tommer.vermaas@wur.nl | 0317-481654
% 
% See also: JARKUS, snctools, jarkus_horizontalVolume,
% jarkus_verticalVolume

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2010 Alterra
%       Tommer Vermaas
%       tommer.vermaas@wur.nl | 0317-481654
%
%       Postbus 47
%       6700 AA Wageningen
%
%   This library is free software: you can redistribute it and/or
%   modify it under the terms of the GNU Lesser General Public
%   License as published by the Free Software Foundation, either
%   version 2.1 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
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library. If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

%% settings

% jarkus database for transect info
% url          = jarkus_url
% url = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc';
url='D:\JARKUS\transect_28022013_opendap.nc';

% general settings
OPT.fillvalue    = NaN;
OPT.refdatenum   = datenum(1970,1,1);
OPT.outputfile   = 'DWI.nc';
OPT.disp         = true;
OPT.dump         = true;

%% Read data form jarkus transect.nc
id         = nc_varget(url,'id');
MHW        = nc_varget(url,'mean_high_water');
MLW        = nc_varget(url,'mean_low_water');
lat        = nc_varget(url,'lat');
lon        = nc_varget(url,'lon');
cross      = nc_varget(url,'cross_shore');
along      = nc_varget(url,'alongshore');
areaname   = nc_varget(url,'areaname');
areacode   = nc_varget(url,'areacode');
time       = nc_varget(url,'time');
jarkusyear = str2num(datestr(nc_cf_time(url),'yyyy')); 
x          = nc_varget(url,'x');
y          = nc_varget(url,'y');
rsp_x      = nc_varget(url,'rsp_x');
rsp_y      = nc_varget(url,'rsp_y');
angle      = nc_varget(url,'angle');

%% Define depthzone boundaries and create empty variables
upperB = [-1 -1 -1 -5 -5 -8];
lowerB = [-5 -8 -12 -8 -12 -12];

landwB      = NaN(length(time),length(id), length(upperB));
seawB       = landwB;
hVolume     = landwB;
x_hVolume   = landwB;
hVolume_x   = landwB;
hVolume_y   = landwB;
hVolume_lon = landwB;
hVolume_lat = landwB;
vVolume     = landwB;
x_vVolume   = landwB;
vVolume_x   = landwB;
vVolume_y   = landwB;
vVolume_lon = landwB;
vVolume_lat = landwB;
uB_vVolume  = landwB;

%% Calculate variables for referenceyear (1990)
jref      = find(jarkusyear==1990);

disp('calculating referenceyear');
for d=1:length(upperB)
    for i=1:length(id)
        z = nc_varget(url,'altitude',[jref-1,i-1,0],[1,1,-1]);
        [hVolume(jref,i,d),landwB(jref,i,d),seawB(jref,i,d)] = ...
                jarkus_horizontalVolume(cross(~isnan(z)), z(~isnan(z)), upperB(d), lowerB(d));
    end
end
vVolume(jref,:,:)    = hVolume(jref,:,:);
uB_vVolume(jref,:,:) = repmat(upperB,length(id),1);

%% Loop through other years and calculate variables
disp('calculating other years');
for t=1:length(time)
    disp(['year=' num2str(t)]);
    if t~=jref %skip reference year, is already calculated
        for d=1:length(upperB)
            disp(['depthzone=' num2str(d)]);
            for i=1:length(id)
                z = nc_varget(url,'altitude',[t-1,i-1,0],[1,1,-1]);

                [hVolume(t,i,d),landwB(t,i,d),seawB(t,i,d)] = ...
                        jarkus_horizontalVolume(cross(~isnan(z)), z(~isnan(z)), upperB(d), lowerB(d));
                [vVolume(t,i,d), uB_vVolume(t,i,d)] = ... 
                        jarkus_verticalVolume(cross(~isnan(z)), z(~isnan(z)), landwB(jref,i,d), seawB(jref,i,d), lowerB(d));

            end
        end
    end %endif
end

for d=1:length(upperB)
    x_hVolume(:,:,d) = landwB(:,:,d) + (hVolume(:,:,d) ./ (upperB(d)-lowerB(d)));
    x_vVolume(:,:,d) = landwB(:,:,d) + (vVolume(:,:,d) ./ (uB_vVolume(:,:,d)-lowerB(d)));
end


%% convert to RD coordinates and WGS84
%to RD
disp('converting coordinates');
for t=1:length(time)
    for d=1:length(upperB)
        for i=1:length(id)
            hVolume_x(t,i,d) = interp1(cross', x(i,:), x_hVolume(t,i,d));
            hVolume_y(t,i,d) = interp1(cross', y(i,:), x_hVolume(t,i,d));
            vVolume_x(t,i,d) = interp1(cross', x(i,:), x_vVolume(t,i,d));
            vVolume_y(t,i,d) = interp1(cross', y(i,:), x_vVolume(t,i,d));
        end
    end
end

%%
%to WGS84
for d=1:length(upperB)
    [hVolume_lon(:,:,d), hVolume_lat(:,:,d)] = ...
        convertCoordinates(hVolume_x(:,:,d),hVolume_y(:,:,d),'CS1.code',28992,'CS2.code',4326);
    [vVolume_lon(:,:,d), vVolume_lat(:,:,d)] = ...
        convertCoordinates(vVolume_x(:,:,d),vVolume_y(:,:,d),'CS1.code',28992,'CS2.code',4326);
end

disp('finished with calculations, build netCDF-file');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Create empty file and add overall meta info
nc_create_empty(OPT.outputfile)

attval =  sprintf('created on %s', datestr(now));

nc_attput(OPT.outputfile, nc_global, 'title'         , 'deep water indicators');
nc_attput(OPT.outputfile, nc_global, 'institution'   , 'Rijkswaterstaat'); 
nc_attput(OPT.outputfile, nc_global, 'source'        , 'on shore measurements');
nc_attput(OPT.outputfile, nc_global, 'history'       , ['Data processed from Jarkus database, ', attval]);
nc_attput(OPT.outputfile, nc_global, 'references'    , 'original source: http://www.rijkswaterstaat.nl , Deltares storage: http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/ , converted with script DWI2nc.m');
nc_attput(OPT.outputfile, nc_global, 'Conventions'   , 'CF-1.5');
nc_attput(OPT.outputfile, nc_global, 'spatial_ref', 'COMPD_CS["Amersfoort / RD New + NAP",PROJCS["Amersfoort / RD New",GEOGCS["Amersfoort",DATUM["Amersfoort",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[565.04,49.91,465.84,-0.40939438743923684,-0.35970519561431136,1.868491000350572,0.8409828680306614],AUTHORITY["EPSG","6289"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","4289"]],PROJECTION["Oblique Stereographic",AUTHORITY["EPSG","9809"]],PARAMETER["central_meridian",5.387638888888891],PARAMETER["latitude_of_origin",52.15616055555556],PARAMETER["scale_factor",0.9999079],PARAMETER["false_easting",155000.0],PARAMETER["false_northing",463000.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","28992"]],VERT_CS["Normaal Amsterdams Peil",VERT_DATUM["Normaal Amsterdams Peil",2005,AUTHORITY["EPSG","5109"]],UNIT["m",1.0],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","5709"]],AUTHORITY["EPSG","7415"]]');
nc_attput(OPT.outputfile, nc_global, 'reference_year', jarkusyear(jref));

%% create dimensions
nc_adddim(OPT.outputfile, 'time', length(time)); 
nc_adddim(OPT.outputfile, 'alongshore', length(along)); 
nc_adddim(OPT.outputfile, 'depthzone', length(upperB)); %new dimension
nc_adddim(OPT.outputfile, 'stringsize', 100);

%% add variables
i = 0;

%% add epsg
[lon_dummy,lat_dummy,OPT2] = convertCoordinates(0,0,'CS1.code',28992,'CS2.code',4326);
epsg = 28992;

i = i + 1;
nc(i).Name      = 'epsg';
nc(i).Nctype    = 'int';
nc(i).Dimension = {};
nc(i).Attribute = struct('Name', ...
   {'grid_mapping_name', ...
    'semi_major_axis', ...
    'semi_minor_axis', ...
    'inverse_flattening', ...
    'latitude_of_projection_origin', ...
    'longitude_of_projection_origin', ...
    'false_easting', ...
    'false_northing', ...
    'scale_factor_at_projection_origin',...
    'comment'}, ...
    'Value', ...
    {OPT2.proj_conv1.method.name,    ...
    OPT2.CS1.ellips.semi_major_axis, ...
    OPT2.CS1.ellips.semi_minor_axis, ...
    OPT2.CS1.ellips.inv_flattening,  ...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'Latitude of natural origin'    )), ...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'Longitude of natural origin'   )),...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'False easting'                 )),...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'False northing'                )),...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'Scale factor at natural origin')),...
    'value is equal to EPSG code'});

%% add original variables
% id
i = i + 1;
nc(i).Name         = 'id';
nc(i).Nctype       = 'int';
nc(i).Dimension    = {'alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'identifier');
nc(i).Attribute(2) = struct('Name', 'comment'   ,'Value', 'sum of area code (x1000000) and alongshore coordinate');

% areacode
i = i + 1;
nc(i).Name         = 'areacode';
nc(i).Nctype       = 'int';
nc(i).Dimension    = {'alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'area code');
nc(i).Attribute(2) = struct('Name', 'comment'   ,'Value', 'codes for the 15 coastal areas as defined by rijkswaterstaat');

% areaname
i = i + 1;
nc(i).Name         = 'areaname';
nc(i).Nctype       = 'char';
nc(i).Dimension    = {'alongshore','stringsize'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'area name');
nc(i).Attribute(2) = struct('Name', 'comment'   ,'Value', 'names for the 15 coastal areas as defined by rijkswaterstaat');

% alongshore
i = i + 1;
nc(i).Name         = 'alongshore';
nc(i).Nctype       = 'double'; 
nc(i).Dimension    = {'alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'alongshore coordinate');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'alongshore coordinate within the 15 coastal areas as defined by rijkswaterstaat');

% time
i = i + 1;
nc(i).Name         = 'time';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time'};
nc(i).Attribute(1) = struct('Name', 'standard_name' ,'Value', 'time');
nc(i).Attribute(2) = struct('Name', 'axis'          ,'Value', 'T');
nc(i).Attribute(3) = struct('Name', 'units'         ,'Value', 'days since 1970-01-01 00:00 +1:00');
nc(i).Attribute(4) = struct('Name', 'comment'       ,'Value', 'measurement year');

%% add new (calculated) variables

% upper boundary (upperB)
i = i + 1;
nc(i).Name         = 'upper_boundary';
nc(i).Nctype       = 'int';
nc(i).Dimension    = {'depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Upper boundary of depth zone');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Depth relative to NAP');

% lower boundary (lowerB)
i = i + 1;
nc(i).Name         = 'lower_boundary';
nc(i).Nctype       = 'int';
nc(i).Dimension    = {'depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Lower boundary of depth zone');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Depth relative to NAP');

% landward boundary (landwB)
i = i + 1;
nc(i).Name         = 'landward_boundary';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Landward boundary of depth zone');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Distance from RSP to upper boundary');

% seaward boundary (seawB)
i = i + 1;
nc(i).Name         = 'seaward_boundary';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Seaward boundary of depth zone');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Distance from RSP to lower boundary');

% volume between horizontal boundaries (hVolume)
i = i + 1;
nc(i).Name         = 'horizontal_volume';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Volume between upper and lower boundaries');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm3/m');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Volume between upper and lower boundaries');

% distance to volume between horizontal boundaries (x_hVolume)
i = i + 1;
nc(i).Name         = 'horizontal_volume_rsp';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Distcance to volume between upper and lower boundaries');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Sum of landwardboundary and normalised volume (divided by upper-lower boundary)');

% position of volume in RD: x (hVolume_x)
i = i + 1;
nc(i).Name         = 'horizontal_volume_x';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in RD: x');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in RD: x');

% position of depthzone volume in RD: y (hVolume_y)
i = i + 1;
nc(i).Name         = 'horizontal_volume_y';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in RD: y');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in RD: y');

% position of volume in WGS84: lat (hVolume_lat)
i = i + 1;
nc(i).Name         = 'horizontal_volume_lat';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in WGS84: lat');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'degrees_west');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in WGS84: lat');

% position of depthzone volume in WGS84: lon (hVolume_lon)
i = i + 1;
nc(i).Name         = 'horizontal_volume_lon';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in WGS84: lon');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'degrees_east');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in WGS84: lon');

% volume between vertical boundaries (vVolume)
i = i + 1;
nc(i).Name         = 'vertical_volume';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Volume between landward and seaward boundaries');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm3/m');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Reference level is lower boundary');

% distance to volume between vertical boundaries (x_vVolume)
i = i + 1;
nc(i).Name         = 'vertical_volume_rsp';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Distcance to volume between landward and seaward boundaries');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Sum of landwardboundary and normalised volume (divided by uB_vVolume-reference)');

% position of volume in RD: x (vVolume_x)
i = i + 1;
nc(i).Name         = 'vertical_volume_x';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in RD: x');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in RD: x');

% position of volume in RD: y (vVolume_y)
i = i + 1;
nc(i).Name         = 'vertical_volume_y';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in RD: y');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in RD: y');

% position of volume in WGS84: lat (vVolume_lat)
i = i + 1;
nc(i).Name         = 'vertical_volume_lat';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in WGS84: lat');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'degrees_west');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in WGS84: lat');

% position of depthzone volume in WGS84: lon (vVolume_lon)
i = i + 1;
nc(i).Name         = 'vertical_volume_lon';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Position of volume in WGS84: lon');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'degrees_east');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Position of volume in WGS84: lon');

% z at landward boundary (uB_vVolume)
i = i + 1;
nc(i).Name         = 'vertical_volume_upperboundary';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore','depthzone'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Depth at landward boundary');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Depth relative to NAP');

%% create variables with attributes
for i = 1:length(nc)
    if OPT.disp;disp(['adding ',num2str(i),' ',nc(i).Name]);end
    nc_addvar(OPT.outputfile, nc(i));
end

%% store variables
nc_varput(OPT.outputfile, 'id'                      , id);
nc_varput(OPT.outputfile, 'areacode'                , areacode);
nc_varput(OPT.outputfile, 'areaname'                , areaname);
nc_varput(OPT.outputfile, 'alongshore'              , along);
nc_varput(OPT.outputfile, 'time'                    , time);

nc_varput(OPT.outputfile, 'upper_boundary'          , upperB);
nc_varput(OPT.outputfile, 'lower_boundary'          , lowerB);
nc_varput(OPT.outputfile, 'landward_boundary'       , landwB);
nc_varput(OPT.outputfile, 'seaward_boundary'        , seawB);

nc_varput(OPT.outputfile, 'horizontal_volume'       , hVolume);
nc_varput(OPT.outputfile, 'horizontal_volume_rsp'   , x_hVolume);
nc_varput(OPT.outputfile, 'horizontal_volume_x'     , hVolume_x);
nc_varput(OPT.outputfile, 'horizontal_volume_y'     , hVolume_y);
nc_varput(OPT.outputfile, 'horizontal_volume_lon'   , hVolume_lon);
nc_varput(OPT.outputfile, 'horizontal_volume_lat'   , hVolume_lat);

nc_varput(OPT.outputfile, 'vertical_volume'         , vVolume);
nc_varput(OPT.outputfile, 'vertical_volume_rsp'     , x_vVolume);
nc_varput(OPT.outputfile, 'vertical_volume_x'       , vVolume_x);
nc_varput(OPT.outputfile, 'vertical_volume_y'       , vVolume_y);
nc_varput(OPT.outputfile, 'vertical_volume_lon'     , vVolume_lon);
nc_varput(OPT.outputfile, 'vertical_volume_lat'     , vVolume_lat);
nc_varput(OPT.outputfile, 'vertical_volume_upperboundary' , uB_vVolume);

%% check
if OPT.dump
    nc_dump(OPT.outputfile);
end

%% end
