%................Run kwaj models based on cluster analysis.................
%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 two files (broken up to conserve memory). 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
outFname1 = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File1.mat');
outFname2 = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File2.mat');

% 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 = 1; %Time interval of model output, in s

% 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,'xbeach_input','temp.grd'); % Temporary file created with spatially-varying friction coefficient
dlmwrite(cf_fname,bathy.cf,'\t');

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', {'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, 'nglobalvar',{'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(outFname1,'modelOut','-v7.3');
end

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', {'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, 'nglobalvar',{'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(outFname1,'modelOut','-v7.3');
end

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', {'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, 'nglobalvar',{'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(outFname1,'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', {'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, 'nglobalvar',{'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(outFname2,'modelOut','-v7.3');
end

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', {'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, 'nglobalvar',{'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(outFname2,'modelOut','-v7.3');
end

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', {'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, 'nglobalvar',{'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(outFname2,'modelOut','-v7.3');
end