
function [Zsmax,Hsmax,Tpmax,DirMmax] = hazardbc_example(inpath,zs0file,wavefile,outpath,outfile,wtype,Dori)
% Matlab script generated by netCDF kickstarter
% at 2015-12-14T12:20Z
% http://publicwiki.deltares.nl/display/OET/netCDF%20kickstarter
%% READ FILES FROM MODEL RESULTS (XBEACH)

% inpath  = 'E:\RISCKIT_FEWS3\Risc-kit\Modules\XBeach\PF\swan_d3d_hindcast\model'
% outpath = 'E:\RISC-KIT_project\WP3\BayesianNetwork\model\nc_out_test'
% zs0file = 'zs0file.txt' 
% wavefile= 'bcfile_B0.txt' % output is for instat=5 or 'jonswap.txt'
% outfile = 'hazardbc_test.nc'
% wtype   = 1 for SWAN2D files and 2 for jonshap table file
% Dori    = Wave direction perpendicular to the coastline  Ria Formosa 220

zs0input  = fullfile(inpath,zs0file);
waveinput  = fullfile(inpath,wavefile);
ncfile   = fullfile(outpath,outfile);
%% LOAD INPUT FILES 
zs0=load(zs0input);
switch wtype
    case 1 % for SWAN2D instat =5
        T=readtable(waveinput);
        C=cell2mat(table2cell(T));
        inds=regexp(str2mat(C(1,:)),' ');
        
        wtime=str2num(C(:,1:inds(1)-1));
        wfile= str2mat(C(:,inds(end)+1:end));
        
        for ii=1:length(wtime)
            spec = swan_io_spectrum(fullfile(inpath,wfile(ii,:)))
            spec1D=abs(trapz(spec.directions,squeeze(spec.VaDens(1,:,:))'));
            Hs(ii) = swan_hs(spec.frequency,spec1D)
            Tm01(ii) = swan_tm01(spec.frequency,spec1D)
            Tm02(ii) = swan_tm02(spec.frequency,spec1D)
            Tp(ii)   =1/spec.frequency(find(spec1D==max(spec1D)))
            spec1DD=abs(trapz(spec.frequency,squeeze(spec.VaDens(1,:,:))));
            DirP(ii)   =spec.directions(find(spec1DD==max(spec1DD)))
            DirM(ii) = sum(spec.directions.*spec1DD/sum(spec1DD))
        end
        
        %% INTERPOLATE
        tstopT= max(zs0(end,1));
        tstopW= max(cumsum(wtime)-3600)
        if tstopT == tstopW
            Hs_int=interp1(cumsum(wtime)-3600,Hs,zs0(:,1));
            Tm01_int=interp1(cumsum(wtime)-3600,Tm01,zs0(:,1));
            Tm02_int=interp1(cumsum(wtime)-3600,Tm02,zs0(:,1));
            Tp_int=interp1(cumsum(wtime)-3600,Tp,zs0(:,1));
            DirP_int=interp1(cumsum(wtime)-3600,DirP,zs0(:,1));
            DirM_int=interp1(cumsum(wtime)-3600,DirM,zs0(:,1));
        else
            disp('END TIME OF TIMESIRIES IS NOT THE SAME')
        end
    case 2 % for jonswap.txt instat= jons_table
        wave=load(wavefile);
 
        wtime= cumsum(wave(:,6));
        Hs   = wave(:,1);
        Tp   = wave(:,2);
        DirM = wave(:,3);
        
        tstopT= max(zs0(end,1));
        tstopW= max(wtime);
        if tstopT == tstopW
            Hs_int=interp1(wtime-3600,Hs,zs0(:,1));
            Tp_int=interp1(wtime-3600,Tp,zs0(:,1));  
            DirM_int=interp1(wtime-3600,DirM,zs0(:,1));
        else
            disp('END TIME OF TIMESIRIES IS NOT THE SAME')
        end
end



%% FIND MAX FOR EACH PARAMETER (or any other Statistic)

indZs0= find(zs0(:,2)==max(zs0(:,2)));
indHs= find(Hs_int==max(Hs_int));
% indTm01= find(Tm01_int==max(Tm01_int));
% indTm02= find(Tm02_int==max(Tm02_int));
indTp= find(Tp_int==max(Tp_int));
%angP=DirP_int-Dori;
angM=DirM_int-Dori;
%indDirP= find(angP==min(angP));
indDirM= find(angM==min(angM));
%%  PARAMETERS SET
tstep= mean(diff(zs0(:,1))); % dt in second 
time_span= 6*60*60/tstep; % time window in index (to find max in a 12h window)
indTime=unique([indZs0 indHs indTp]);
time_dim=size(indTime);
timemax = zs0(indTime,1);

indTime_min = max(indTime-time_span,1);
indTime_max = min(indTime+time_span,tstopT/tstep);
for ii= 1:time_dim(2)
    Zsmax(ii)=max(zs0(indTime_min(ii):indTime_max(ii),2)); % find max tide level 
    Hsmax(ii)=mean(Hs_int(indTime_min(ii):indTime_max(ii))); % find mean Hs during tide
    %Tm01max=mean(Tm01_int(indTime_min(ii):indTime_max(ii)));
    Tpmax(ii)=mean(Tp_int(indTime_min(ii):indTime_max(ii)));% find mean Tp during tide
    DirMmax(ii)=mean(DirM_int(indTime_min(ii):indTime_max(ii))); %find mean Dir during tide
end



time_units ='minutes since 1970-01-01 00:00:00.0 +0000';
time_now =now;
time_max= (timemax-datenum(time_units(15:end-6))).*(60*24);
analysis_time= (time_now-datenum(time_units(15:end-6))).*(60*24);

%% CREATE EMPTY NETCDF FILE
%ncfile = string with the name of the output file  e.g.

nc_create_empty(ncfile);

%% ADD DIMENSIONS%% ADD DIMENSIONS
nc_adddim(ncfile, 'time', time_dim(2));
nc_adddim(ncfile, 'analysis_time', 1);
nc_adddim(ncfile, 'stations', 1);
nc_adddim(ncfile, 'char_leng_id', 64);
nc_adddim(ncfile, 'char_leng_name', 255);
%% ADD GLOBAL ATTRIBUTES
% see http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/formats/DataDiscoveryAttConvention.html
nc_attput(ncfile, nc_global, 'Conventions', 'CF-1.6');
nc_attput(ncfile, nc_global, 'title', 'Data');
nc_attput(ncfile, nc_global, 'institution', 'UAlg');
nc_attput(ncfile, nc_global, 'source', 'Export NETCDF-CF_TIMESEIRES that mimic FEWS output');
nc_attput(ncfile, nc_global, 'history', [datestr(now,31) ' exported']);
nc_attput(ncfile, nc_global, 'references', '');
nc_attput(ncfile, nc_global, 'Metadata_Conventions', 'Unidata Dataset Discovery v1.0');
nc_attput(ncfile, nc_global, 'summary', 'Data exported WITHOUT using FEWS');
nc_attput(ncfile, nc_global, 'date_created', datestr(now,31));

%% ADD VARIABLES
var = struct(...
    'Name', 'time',...
    'Datatype', 'double',...
    'Dimension', {{'time'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'time', 'long_name', '');
nc_attput(ncfile, 'time', 'standard_name', 'time');
nc_attput(ncfile, 'time', 'units', 'minutes since 1970-01-01 00:00:00.0 +0000');
nc_attput(ncfile, 'time', 'axis', 'T');

var = struct(...
    'Name', 'analysis_time',...
    'Datatype', 'double',...
    'Dimension', {{'analysis_time'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'analysis_time', 'long_name', 'forecast_reference_time');
nc_attput(ncfile, 'analysis_time', 'standard_name', 'forecasr_reference_time');
nc_attput(ncfile, 'analysis_time', 'units', 'minutes since 1970-01-01 00:00:00.0 +0000');


var = struct(...
    'Name', 'lat',...
    'Datatype', 'double',...
    'Dimension', {{'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'lat', 'long_name', 'Station coordinates, latitude');
nc_attput(ncfile, 'lat', 'standard_name', 'latitude');
nc_attput(ncfile, 'lat', 'units', 'degrees_north');
nc_attput(ncfile, 'lat', 'axis', 'Y');
nc_attput(ncfile, 'lat', '_FillValue', 9.96921e+03);


var = struct(...
    'Name', 'lon',...
    'Datatype', 'double',...
    'Dimension', {{'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'lon', 'long_name', 'Station coordinates, longitude');
nc_attput(ncfile, 'lon', 'standard_name', 'longitude');
nc_attput(ncfile, 'lon', 'units', 'degrees_east');
nc_attput(ncfile, 'lon', 'axis', 'X');
nc_attput(ncfile, 'lon', '_FillValue', 9.96920996838687e+036);


var = struct(...
    'Name', 'x',...
    'Datatype', 'double',...
    'Dimension', {{'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'x', 'long_name', 'x coordinate according to WGS 1984');
nc_attput(ncfile, 'x', 'standard_name', 'projection_x_coordinate');
nc_attput(ncfile, 'x', 'units', 'degrees');
nc_attput(ncfile, 'x', 'axis', 'X');
nc_attput(ncfile, 'x', '_FillValue', 9.96920996838687e+036);

var = struct(...
    'Name', 'y',...
    'Datatype', 'double',...
    'Dimension', {{'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'y', 'long_name', 'y coordinate according to WGS 1984');
nc_attput(ncfile, 'y', 'standard_name', 'projection_y_coordinate');
nc_attput(ncfile, 'y', 'units', 'degrees');
nc_attput(ncfile, 'y', 'axis', 'Y');
nc_attput(ncfile, 'y', '_FillValue', 9.96920996838687e+036);

var = struct(...
    'Name', 'station_id',...
    'Datatype', 'char',...
    'Dimension', {{'char_leng_id'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'station_id', 'long_name', 'station identification code');
nc_attput(ncfile, 'station_id', 'cf_role', 'timeseries_id');

var = struct(...
    'Name', 'station_name',...
    'Datatype', 'char',...
    'Dimension', {{'char_leng_name'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'station_name', 'long_name', 'station name');

var = struct(...
    'Name', 'Zs',...
    'Datatype', 'double',...
    'Dimension', {{'time' 'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'Zs', 'long_name', 'Zs');
nc_attput(ncfile, 'Zs', 'standard_name', 'sea_surface_height_above_sea_level');
nc_attput(ncfile, 'Zs', 'units', 'm');
nc_attput(ncfile, 'Zs', 'coordinates', 'lat lon time'); % analysis_time in FEWS nc but (bug???) why lat lot?
nc_attput(ncfile, 'Zs', '_FillValue', 9.96920996838687e+036);

var = struct(...
    'Name', 'Hs',...
    'Datatype', 'double',...
    'Dimension', {{'time' 'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'Hs', 'long_name', 'Hs');
nc_attput(ncfile, 'Hs', 'standard_name', 'sea_surface_wave_significant_height');
nc_attput(ncfile, 'Hs', 'units', 'm');
nc_attput(ncfile, 'Hs', 'coordinates', 'lat lon time'); % analysis_time in FEWS nc but (bug???) why lat lot?
nc_attput(ncfile, 'Hs', '_FillValue', 9.96920996838687e+036);

var = struct(...
    'Name', 'th',...
    'Datatype', 'double',...
    'Dimension', {{'time' 'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'th', 'long_name', 'th');
nc_attput(ncfile, 'th', 'standard_name', 'sea_surface_wave_from_direction');
nc_attput(ncfile, 'th', 'units', 'degree');
nc_attput(ncfile, 'th', 'coordinates', 'lat lon time'); % analysis_time in FEWS nc but (bug???) why lat lot?
nc_attput(ncfile, 'th', '_FillValue', 9.96920996838687e+036);

var = struct(...
    'Name', 'Tp',...
    'Datatype', 'double',...
    'Dimension', {{'time' 'stations'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'Tp', 'long_name', 'Tp');
nc_attput(ncfile, 'Tp', 'standard_name', 'sea_surface_wave_period_at_variance_spectral_density_maximum');
nc_attput(ncfile, 'Tp', 'nodc_name', '');
nc_attput(ncfile, 'Tp', 'units', 's');
nc_attput(ncfile, 'Tp', 'coordinates', 'lat lon time'); % analysis_time in FEWS nc but (bug???) why lat lot?
nc_attput(ncfile, 'Tp', '_FillValue', 9.96920996838687e+036);

%% Store Variables
nc_varput(ncfile, 'time', time_max)
nc_varput(ncfile, 'analysis_time', analysis_time)
nc_varput(ncfile, 'Zs', Zsmax')
nc_varput(ncfile, 'Hs', Hsmax')
nc_varput(ncfile, 'Tp', Tpmax')
nc_varput(ncfile, 'th', DirMmax')