%................Run kwaj models based on cluster analysis - 10 Hz.........
%Author: James Bramante
%Date: November 8, 2019
%
% Produce models for a range of clustered wave climate parameters
% This scripts runs XBeach for 6 wave climate/sea-level scenarios and saves
% the results in six files. This script produces results similar to 
% Kwajalein_xbeach_runs_publication, but with 10 Hz resolution at several
% locations across the reef transect, instead of 1 Hz resolution at all
% points across the transect. Currently this
% script is set to model each of the wave climate clusters in each scenario
% (~30 clusters in each scenario) for a single sea-level associated with
% that scenario. User can input different tidal elevations to offset the
% sea level in those scenarios to simulate e.g. high or low tide.
%
% This script requires the presence of the XBeach Toolbox for MATLAB. This
% toolbox should be available at: https://oss.deltares.nl/web/xbeach/tools
%
% The three data files required for this script are:
% 1 Kwajalein_transect_model_bathymetry.mat - The file containing a struct
%                                             outlining the bathymetry of 
%                                             the modeled reef transect
% 2 Storlazzi_historical_all_GCM_clusters.mat - These two files contain the
%                                               clustered wave data from 
%                                               the WaveWatch III forecasts 
%                                               of Shope et al. (2016).
%                                               Only Historical and RCP8.5 
%                                               results for the 2080-2100 
%                                               period were used
% 3 Storlazzi_2100_85_all_GCM_clusters.mat

close all; clear all; clc;

%% Finalish Variables
% User input directories
baseDir = ''; % Directory containing all input files
BINARYDIR = ''; % Directory containing the XBeachX executable
outputDir = ''; % Directory to save XBeach model setup and results before collation

% Files for input and output
HIST_CLUSTERS = fullfile(baseDir,'Storlazzi_historical_all_GCM_clusters.mat'); % File containing Kmeans clusters of Historical scenario modeled wave climate
GCM2100_CLUSTERS = fullfile(baseDir,'Storlazzi_2100_85_all_GCM_clusters.mat'); %File containing Kmeans clusters of RCP8.5 2100 scenario modeled wave climate
bathyFile = fullfile(baseDir,'Kwajalein_transect_model_bathymetry.mat'); % Bathymetric file
outFnamePrefix = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_10Hz_File');

% Tidal offsets, in m, for each model scenario
histTideOffset = 0;
waveClimateTideOffset = 0;
slr0p5TideOffset = 0.5;
slr1TideOffset = 1;
slr1p5TideOffset = 1.5;
slr2TideOffset = 2;

ALFA = 180; % Bathymetric transect rotation from W-E
MHW = 1.978 - 1.457; %Mean high water at Kwajalein tidal station
MLW = 0.929 - 1.457; %Mean low water at Kwajalein tidal station
SPINUP = 6*3600; %Spinup time for each model run in seconds
MODEL_LENGTH = 4096; %Length of each model run in seconds
TINT = 0.1; %Time interval of model output, in s
outPoints = [66,100,120,143,148.25]*4-3; %Grid indices at which to report output variables

% Read in the cluster values
histClusters = load(HIST_CLUSTERS);
histClusters = histClusters.clusters;
gcmClusters = load(GCM2100_CLUSTERS);
gcmClusters = gcmClusters.clusters;


% Compile model inputs
% tides = [MLW,0,MHW];
tides = [0];
duration = SPINUP + MODEL_LENGTH; % Total model duration
modelOut = struct();
bathy = load(bathyFile);
bathy = bathy.bathy;

% Write out cf file
cf_fname = fullfile(baseDir,'temp.grd'); % Temporary file created with spatially-varying friction coefficient
dlmwrite(cf_fname,bathy.cf,'\t');
outPoints2 = {};
for ii = 1 : length(outPoints)
    outPoints2{ii} = sprintf('%5.1f 0.0',bathy.x(outPoints(ii)));
end


for x = 1 : size(histClusters,1)
    for y = 1 : length(tides)
        % Run the Historical clusters, no SLR scenario first
        xbm = xb_generate_model( ...
            'bathy',    {'x', bathy.x, 'z', bathy.z, 'optimize', 0, 'ne', zeros(1,length(bathy.x)), 'posdwn', 0},...
            'tide',     {'front', tides(y)+histTideOffset, 'back', tides(y)+histTideOffset}, ...
            'waves',    {'type','jonswap','Hm0', histClusters(x,2), 'Tp', histClusters(x,3), 'duration', duration, 'mainang', histClusters(x,4)+ALFA}, ...
            'settings', {'nglobalvar',0,'nonh',1,'sprdthr',0.08,'dtheta', 10,'tstart',duration-MODEL_LENGTH,'tstop',duration,...
            'thetanaut' , 1, 'dy', bathy.dy, 'alfa',0,'maxbrsteep',0.5, 'morphology', 0,'swave',0,'sedtrans',0,...
            'wind',0,'lat',8.5,'xori', 0, 'fw', 0, 'bedfriction', 'cf', 'bedfricfile', cf_fname,...
            'tint', TINT, 'npoints',outPoints2,'npointvar',{'hh','u','taubx'}} ...
        );

        %% Now run the model! Note that the 'binary' option points to the FOLDER
        %containing xbeach.exe, not xbeach.exe itself
        xbr = xb_run(xbm,'path',outputDir,'binary',BINARYDIR);
        while xb_check_run(xbr)%Wait for model to stop running
        end

        modelOut.gcmhist.(['c',num2str(x)]).(regexprep(regexprep(['t',num2str(tides(y))],'-','n'),'\.','')) = xb_read_output(xbr);
    end
    save(strcat(outFnamePrefix,'1.mat'),'modelOut','-v7.3');
end

clear modelOut;
modelOut = struct();
for x = 1 :size(gcmClusters,1)
    for y = 1 : length(tides)
        % Run the GCM clusters, no SLR scenario
        xbm = xb_generate_model( ...
            'bathy',    {'x', bathy.x, 'z', bathy.z, 'optimize', 0, 'ne', zeros(1,length(bathy.x)), 'posdwn', 0},...
            'tide',     {'front', tides(y)+waveClimateTideOffset, 'back', tides(y)+waveClimateTideOffset}, ...
            'waves',    {'type','jonswap','Hm0', gcmClusters(x,2), 'Tp', gcmClusters(x,3), 'duration', duration, 'mainang', gcmClusters(x,4)+ALFA}, ...
            'settings', {'nglobalvar',0,'nonh',1,'sprdthr',0.08,'dtheta', 10,'tstart',duration-MODEL_LENGTH,'tstop',duration,...
            'thetanaut' , 1, 'dy', bathy.dy, 'alfa',0,'maxbrsteep',0.5, 'sedtrans', 0, 'morphology', 0,'swave',0,...
            'wind',0,'lat',8.5,'xori', 0, 'fw', 0, 'bedfriction', 'cf', 'bedfricfile', cf_fname,...
            'tint', TINT, 'npoints',outPoints2,'npointvar',{'hh','u','taubx'}} ...
        );

        %% Now run the model! Note that the 'binary' option points to the FOLDER
        %containing xbeach.exe, not xbeach.exe itself
        xbr = xb_run(xbm,'path',outputDir,'binary',BINARYDIR);
        while xb_check_run(xbr)%Wait for model to stop running
        end

        modelOut.gcm2100_withoutWL.(['c',num2str(x)]).(regexprep(regexprep(['t',num2str(tides(y))],'-','n'),'\.','')) = xb_read_output(xbr);
    end
    save(strcat(outFnamePrefix,'2.mat'),'modelOut','-v7.3');
end

clear modelOut;
modelOut = struct();
for x = 1 : size(histClusters,1)
    for y = 1 : length(tides)
        % Run the 2 m SLR scenario next
        xbm = xb_generate_model( ...
            'bathy',    {'x', bathy.x, 'z', bathy.z, 'optimize', 0, 'ne', zeros(1,length(bathy.x)), 'posdwn', 0},...
            'tide',     {'front', tides(y)+slr2TideOffset, 'back', tides(y)+slr2TideOffset}, ...
            'waves',    {'type','jonswap','Hm0', histClusters(x,2), 'Tp', histClusters(x,3), 'duration', duration, 'mainang', histClusters(x,4)+ALFA}, ...
            'settings', {'nglobalvar',0,'nonh',1,'sprdthr',0.08,'dtheta', 10,'tstart',duration-MODEL_LENGTH,'tstop',duration,...
            'thetanaut' , 1, 'dy', bathy.dy, 'alfa',0,'maxbrsteep',0.5, 'morphology', 0,'swave',0,'sedtrans',0,...
            'wind',0,'lat',8.5,'xori', 0, 'fw', 0, 'bedfriction', 'cf', 'bedfricfile', cf_fname,...
            'tint', TINT, 'npoints',outPoints2,'npointvar',{'hh','u','taubx'}} ...
        );

        %% Now run the model! Note that the 'binary' option points to the FOLDER
        %containing xbeach.exe, not xbeach.exe itself
        xbr = xb_run(xbm,'path',outputDir,'binary',BINARYDIR);
        while xb_check_run(xbr)%Wait for model to stop running
        end

        modelOut.gcmhist_2.(['c',num2str(x)]).(regexprep(regexprep(['t',num2str(tides(y))],'-','n'),'\.','')) = xb_read_output(xbr);
    end
    save(strcat(outFnamePrefix,'3.mat'),'modelOut','-v7.3');
end

% Create a second output variable/file to avoid files that are too large
% and taking up too much memory
clear modelOut;
modelOut = struct();

for x = 1 : size(histClusters,1)
    for y = 1 : length(tides)
        % Run the 0.5 m SLR scenario
        xbm = xb_generate_model( ...
            'bathy',    {'x', bathy.x, 'z', bathy.z, 'optimize', 0, 'ne', zeros(1,length(bathy.x)), 'posdwn', 0},...
            'tide',     {'front', tides(y)+slr0p5TideOffset, 'back', tides(y)+slr0p5TideOffset}, ...
            'waves',    {'type','jonswap','Hm0', histClusters(x,2), 'Tp', histClusters(x,3), 'duration', duration, 'mainang', histClusters(x,4)+ALFA}, ...
            'settings', {'nglobalvar',0,'nonh',1,'sprdthr',0.08,'dtheta', 10,'tstart',duration-MODEL_LENGTH,'tstop',duration,...
            'thetanaut' , 1, 'dy', bathy.dy, 'alfa',0,'maxbrsteep',0.5, 'morphology', 0,'swave',0,'sedtrans',0,...
            'wind',0,'lat',8.5,'xori', 0, 'fw', 0, 'bedfriction', 'cf', 'bedfricfile', cf_fname,...
            'tint', TINT, 'npoints',outPoints2,'npointvar',{'hh','u','taubx'}} ...
        );

        %% Now run the model! Note that the 'binary' option points to the FOLDER
        %containing xbeach.exe, not xbeach.exe itself
        xbr = xb_run(xbm,'path',outputDir,'binary',BINARYDIR);
        while xb_check_run(xbr)%Wait for model to stop running
        end

        modelOut.gcmhist_0p5.(['c',num2str(x)]).(regexprep(regexprep(['t',num2str(tides(y))],'-','n'),'\.','')) = xb_read_output(xbr);
    end
    save(strcat(outFnamePrefix,'4.mat'),'modelOut','-v7.3');
end

clear modelOut;
modelOut = struct();
for x = 1 : size(histClusters,1)
    for y = 1 : length(tides)
        % Run the 1 m SLR scenario next
        xbm = xb_generate_model( ...
            'bathy',    {'x', bathy.x, 'z', bathy.z, 'optimize', 0, 'ne', zeros(1,length(bathy.x)), 'posdwn', 0},...
            'tide',     {'front', tides(y)+slr1TideOffset, 'back', tides(y)+slr1TideOffset}, ...
            'waves',    {'type','jonswap','Hm0', histClusters(x,2), 'Tp', histClusters(x,3), 'duration', duration, 'mainang', histClusters(x,4)+ALFA}, ...
            'settings', {'nglobalvar',0,'nonh',1,'sprdthr',0.08,'dtheta', 10,'tstart',duration-MODEL_LENGTH,'tstop',duration,...
            'thetanaut' , 1, 'dy', bathy.dy, 'alfa',0,'maxbrsteep',0.5, 'morphology', 0,'swave',0,'sedtrans',0,...
            'wind',0,'lat',8.5,'xori', 0, 'fw', 0, 'bedfriction', 'cf', 'bedfricfile', cf_fname,...
            'tint', TINT, 'npoints',outPoints2,'npointvar',{'hh','u','taubx'}} ...
        );

        %% Now run the model! Note that the 'binary' option points to the FOLDER
        %containing xbeach.exe, not xbeach.exe itself
        xbr = xb_run(xbm,'path',outputDir,'binary',BINARYDIR);
        while xb_check_run(xbr)%Wait for model to stop running
        end

        modelOut.gcmhist_1.(['c',num2str(x)]).(regexprep(regexprep(['t',num2str(tides(y))],'-','n'),'\.','')) = xb_read_output(xbr);
    end
    save(strcat(outFnamePrefix,'5.mat'),'modelOut','-v7.3');
end

clear modelOut;
modelOut = struct();
for x = 1 : size(histClusters,1)
    for y = 1 : length(tides)
        % Run the SLR 1.5 m scenario next
        xbm = xb_generate_model( ...
            'bathy',    {'x', bathy.x, 'z', bathy.z, 'optimize', 0, 'ne', zeros(1,length(bathy.x)), 'posdwn', 0},...
            'tide',     {'front', tides(y)+slr1p5TideOffset, 'back', tides(y)+slr1p5TideOffset}, ...
            'waves',    {'type','jonswap','Hm0', histClusters(x,2), 'Tp', histClusters(x,3), 'duration', duration, 'mainang', histClusters(x,4)+ALFA}, ...
            'settings', {'nglobalvar',0,'nonh',1,'sprdthr',0.08,'dtheta', 10,'tstart',duration-MODEL_LENGTH,'tstop',duration,...
            'thetanaut' , 1, 'dy', bathy.dy, 'alfa',0,'maxbrsteep',0.5, 'morphology', 0,'swave',0,'sedtrans',0,...
            'wind',0,'lat',8.5,'xori', 0, 'fw', 0, 'bedfriction', 'cf', 'bedfricfile', cf_fname,...
            'tint', TINT, 'npoints',outPoints2,'npointvar',{'hh','u','taubx'}} ...
        );

        %% Now run the model! Note that the 'binary' option points to the FOLDER
        %containing xbeach.exe, not xbeach.exe itself
        xbr = xb_run(xbm,'path',outputDir,'binary',BINARYDIR);
        while xb_check_run(xbr)%Wait for model to stop running
        end

        modelOut.gcmhist_1p5.(['c',num2str(x)]).(regexprep(regexprep(['t',num2str(tides(y))],'-','n'),'\.','')) = xb_read_output(xbr);
    end
    save(strcat(outFnamePrefix,'6.mat'),'modelOut','-v7.3');
end


% Now, restructure the data into a more useful format for plotting
fnames1 = {'GCM_forecasts_XBEACHX_Results_Revised_File1.mat','GCM_forecasts_XBEACHX_Results_Revised_File2','GCM_forecasts_XBEACHX_Results_Revised_File3'};
fnames2 = {'GCM_forecasts_XBEACHX_Results_Revised_File4.mat','GCM_forecasts_XBEACHX_Results_Revised_File5','GCM_forecasts_XBEACHX_Results_Revised_File6'};
out_fname = 'Modeled_velocity_10hz.mat';
out_loc_names = {'d261','d397','d477','d569','d590'};
velocity_ind = 3;
depth_ind = 2;

modelOutPts = struct();
for ii = 1 : length(fnames1)
    modelOut = load(fullfile(outputDir,fnames1{ii}));
    modelOut = modelOut.modelOut;
    scenarioName = fields(modelOut);
    scenarioName = scenarioName{1};
    modelOut = modelOut.(scenarioName);
    modelOutPts.(scenarioName) = struct();
    clusters = fields(modelOut);
    for jj = 1 : length(clusters)
        modelOutPts.(scenarioName).(clusters{jj}) = struct();
        tide = fields(modelOut.(clusters{jj}));
        tide = tide{1};
        locs_len = length(modelOut.(clusters{jj}).(tide).data);
        for kk = 2 : locs_len
            modelOutPts.(scenarioName).(clusters{jj}).(out_loc_names{kk-1}).raw = ...
                modelOut.(clusters{jj}).(tide).data(kk).value(:,velocity_ind);
            modelOutPts.(scenarioName).(clusters{jj}).(out_loc_names{kk-1}).rawhh = ...
                modelOut.(clusters{jj}).(tide).data(kk).value(:,depth_ind);
        end
    end
end

modelOutPts2 = struct();
for ii = 1 : length(fnames2)
    modelOut = load(fullfile(outputDir,fnames2{ii}));
    modelOut = modelOut.modelOut;
    scenarioName = fields(modelOut);
    scenarioName = scenarioName{1};
    modelOut = modelOut.(scenarioName);
    modelOutPts2.(scenarioName) = struct();
    clusters = fields(modelOut);
    for jj = 1 : length(clusters)
        modelOutPts2.(scenarioName).(clusters{jj}) = struct();
        tide = fields(modelOut.(clusters{jj}));
        tide = tide{1};
        locs_len = length(modelOut.(clusters{jj}).(tide).data);
        for kk = 2 : locs_len
            modelOutPts2.(scenarioName).(clusters{jj}).(out_loc_names{kk-1}).raw = ...
                modelOut.(clusters{jj}).(tide).data(kk).value(:,velocity_ind);
            modelOutPts2.(scenarioName).(clusters{jj}).(out_loc_names{kk-1}).rawhh = ...
                modelOut.(clusters{jj}).(tide).data(kk).value(:,depth_ind);
        end
    end
end

save(fullfile(outputDir,out_fname),'modelOutPts','modelOutPts2'); 