%................Publication figures for wave climate paper................
%Author: James Bramante
%Date: November 5, 2019

% This script is made to interact with 6 Matlab data files and plot most of
% the figures in Bramante et al (submitted) JGR: Earth Surface. To generate
% the figures, this script must be run to completion once to generate a
% couple additional necessary files (N.B. this script has not been
% optimized for efficiency, and running the whole script could take days). 
% After files are generated, individual figures can be generated by first
% running all of the lines before the first section (%%) break and then
% running all of the code within a section denoted with (%%). To print the
% figures to vector files for editing outside of MATLAB, uncomment the
% print commands located at the end of each block of code producing a
% figure.
% The six files required for this script are:
% 1 GCM_forecasts_XBEACHX_Results_File1.mat - Files containing all model
%                                             output (warning: they are >5
%                                             GB large)
% 2 GCM_forecasts_XBEACHX_Results_File2.mat
% 3 Modeled_velocity_10hz.mat               - File containing model results
%                                             at 10Hz for select points
%                                             along the reef transect
% 3 Kwajalein_transect_model_bathymetry.mat - The file containing a struct
%                                             outlining the bathymetry of 
%                                             the modeled reef transect
% 4 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
% 5 Storlazzi_2100_85_all_GCM_clusters.mat
% 6 Wavewatch_mean_wave_parameters.mat - Mean wave characteristics
%                                        calculated for the closest
%                                        open-ocean cell to Kwajalein's SE
%                                        coast from the NCEP WAVEWATCH III
%                                        reanalysis Phase 1
% REFERENCES
% Ruessink, B.G., Ramaekers, G., & van Rijn, L.C. (2012). On the parameterization of the free-stream non-linear wave orbital motion in nearshore morphodynamic models. Coastal Engineering, 65, 56-63.
% 
% Shope, J. B., Storlazzi, C. D., Erikson, L. H., & Hegermiller, C. A. (2016). Changes to extreme wave climates of islands within the western tropical pacific throughout the 21st century under RCP 4.5 and RCP 8.5, with implications for island vulnerability and sustainability. Global and Planetary Change, 141, 25-38.
% 
% Soulsby, R.L., & Damgaard, J.S. (2005). Bedload sediment transport in coastal waters. Coastal Engineering, 52, 673-689, doi:10.1016/j.coastaleng.2005.04.003
% 
% van Thiel de Vries, J. (2009). Dune erosion during storm surges. Phd thesis. Delft University of Technology., Delft, Netherlands.
% 
% Wu, Z. & Huang, N.E. (2009). Ensemble empirical mode decomposition: A noise-assisted data analysis method. Advances in Adaptive Data Analysis, 1, 1-41.
%

clear all; close all; clc;

% User variables
baseDir = '';
figureOutDir = '';

%% Final variables controlling figure size, etc
% Figure style variables
figureFontSize = 8;
maxWidth = 19.0; % Max and min dimensions of figures, in cm
minWidth = 9.5;
maxLength = 23.0;
minLength = 11.5;
pltColors = {[0.5,0,0.5],[0,0,1],[0,0,1],[0,0.3922,1],[0,0.7843,1],[0,0.5,0.5]};
pltDashes = {'-','--','-','--','-','--'};
pltSymbols = {'^','s','d','o','*'};

% Physical and numerical constants
rho = 1027; % Water density in kg/m^3
rhoS = 1850; % Carbonate sediment density in kg/m^3
phiS = 32*pi/180; % Angle of repose of sand, in radians
G = 9.8; % Gravitational acceleration, in m/s^2
TS_LENGTH = 4097; % Length of output model timeseries
MAXITER = 50000; % Maximum number of iterations for Monte Carlo simulation
modelSampleFreq = 1; % Frequency, in Hz, of output model timeseries
McCallCoeff = 1; % Coefficient applied to inertial contribution to shear stress
phiTau = 20*pi/180; % Phase angle separation between velocity and shear stress

% Indices of model output for water depth (HH), bottom shear stress (TAU), 
% and depth-averaged velocity in x-direction (U)
HHIND = 2;
TAUIND = 3;
UIND = 4;

% Model transect indices from which to pull shear stress distributions
distInds = [261,397,477,569,590];
distcf = [0.15,0.015,0.015,0.015,0.015]; % Friction coefficients for each shear stress distribution
distNames = {'Fore reef','Outer reef flat','Mid reef flat','Inner reef flat','Beach'};

% Linear indices of reef regions in bathy vector
reefCrestInd = 353; 
shallowReefInd = 25;
innerFlatInd = 441;
beachInd = 561;
backBeachInd = 641;

% Linear index of bathy vector and depths (m) of wave gauges
e1Loc = 89; e1El = -8.22;
e2Loc = 413; e2El = -0.362;
e3Loc = 477; e3El = -0.438;
e4Loc = 573; e4El = -0.048;

% Load clusters
cluster2100 = fullfile(baseDir,'Storlazzi_2100_85_all_GCM_clusters.mat');
histcluster = fullfile(baseDir,'Storlazzi_historical_all_GCM_clusters.mat');
% Clusters with oblique waves > 180 degrees from shore normal are excluded
gcmHistIncidentClusters = ones(30,1);
gcmHistIncidentClusters([1,15,18,30]) = 0;
gcm2100IncidentClusters = ones(30,1);
gcm2100IncidentClusters([4,17,26,27]) = 0;
histcluster = load(histcluster);
histcluster = histcluster.clusters;
histcluster(:,end) = histcluster(:,end) .* gcmHistIncidentClusters;
histcluster(:,end) = histcluster(:,end) ./ sum(histcluster(:,end)); %Re-normalize annual proportion for each wave cluster, ignoring oblique clusters
cluster2100 = load(cluster2100);
cluster2100 = cluster2100.clusters;
cluster2100(:,end) = cluster2100(:,end) .* gcm2100IncidentClusters;
cluster2100(:,end) = cluster2100(:,end) ./ sum(cluster2100(:,end)); %Re-normalize annual proportion for each wave cluster, ignoring oblique clusters

%..........................................................................
%% 1.) Figure 2. Bathymetry for XBeach model
bathyFname = fullfile(baseDir,'Kwajalein_transect_model_bathymetry.mat');
load(bathyFname);
xBounds = [-150,150]; % Distance from reef crest over which to plot
MHW = 1.088 + 0.0117 - 0.567; %Station MHW datum adjusted for local MSL during deployment, in m
MLW = 0.039 + 0.0117 - 0.567; %Station MLW datum adjusted for local MSL during deployment, in m

% Initialize figure
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,minWidth,minLength*0.5],'Name','Model bathymetry');
fig.Renderer = 'painters';

bathyx = bathy.x - bathy.x(reefCrestInd); % Make reef crest the origin
plot(bathyx,bathy.z,'k-'); % Plot bathymetry
hold on;
plot([min(xBounds),bathyx(find(bathy.z>=0,1,'first'))],[0,0],'b--'); % Plot MSL line
scatter(bathyx([e1Loc,e2Loc,e3Loc,e4Loc]),[e1El,e2El,e3El,e4El],50,'g^','filled'); % Plot locations of wave sensors
% Indicate location of fringing reef regions
plot(bathyx([shallowReefInd,shallowReefInd]),[bathy.z(shallowReefInd)-1,bathy.z(shallowReefInd)+1],'k-');
plot(bathyx([reefCrestInd,reefCrestInd]),[bathy.z(reefCrestInd)-1,bathy.z(reefCrestInd)+1],'k-');
plot(bathyx([innerFlatInd,innerFlatInd]),[bathy.z(innerFlatInd)-1,bathy.z(innerFlatInd)+1],'k-');
plot(bathyx([beachInd,beachInd]),[bathy.z(beachInd)-1,bathy.z(beachInd)+1],'k-');
plot(bathyx([backBeachInd,backBeachInd]),[bathy.z(backBeachInd)-1,bathy.z(backBeachInd)+1],'k-');
% Plot lines indicating low and high tide
[~,MHWind] = min(abs(bathy.z(bathyx<200) - MHW));
[~,MLWind] = min(abs(bathy.z(bathyx<200) - MLW));
plot(bathyx(MHWind-5:MHWind+5),repmat(MHW,1,11),'b-');
plot(bathyx(MLWind-5:MLWind+5),repmat(MLW,1,11),'b-');
xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
xlim(xBounds);
ylim([-10,5]);
set(gca,'FontSize',figureFontSize,'YTick',-10:2.5:5);
set(gca,'Xdir','reverse');

% print(fullfile(figureOutDir,'modelBathy'),'-depsc','-r600');

%..........................................................................

%..........................................................................
%% 2.) Figure 1 Map of mean wave energy across the Pacific

% Get some geographic info to plot
waveFile = fullfile(baseDir,'Wavewatch_mean_wave_parameters.mat'); % WaveWatch III mean wave heights, periods, etc
kwajCoord = [167.73,8.73];
land = shaperead('landareas','UseGeoCoords',false); % Load MATLAB-supplied land shapefiles
land(3).X(land(3).X<0) = 360+land(3).X(land(3).X<0);
load(waveFile);
mapc = [0 0 0];

%Initialize figure
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength],'Name','Mean Hs');
fig.Renderer = 'painters';

pcolor(lons,lats,meanHs./count); shading flat; axis equal tight;
hold on;
c = colorbar; ylabel(c,'Mean H_s (m)','FontSize',figureFontSize); caxis([0,6]);
colormap(gca,flipud(colormap('hot')));
mapshow(land,'SymbolSpec',makesymbolspec('Polygon',{'Default','FaceColor',[0,0.35,0],'EdgeColor',mapc,'FaceAlpha',1}));
xlim([min(lons(:)),max(lons(:))]); ylim([min(lats(:)),max(lats(:))]);
scatter(kwajCoord(1),kwajCoord(2),50,'^','MarkerEdgeColor','k','MarkerFaceColor',[0,0.35,0]);
ylabel('Latitude (^{\circ})','FontSize',figureFontSize);
xlabel('Longitude (^{\circ})','Fontsize',figureFontSize);
set(gca,'FontSize',figureFontSize);

% print(fullfile(figureOutDir,'Pacific_wave_height'),'-depsc','-r600');
%..........................................................................

%..........................................................................
%% 3.) Figure 6 Climatic change in mean wave characteristics
%Filenames pointing to data and variables for plotting/processing data
fileName = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File1.mat');
fileName2 = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File2.mat');
bathyfile = fullfile(baseDir,'Kwajalein_transect_model_bathymetry.mat');
bounds = [193,637]; % Indexes of x-grid cells between which to plot
inundThresh = 0.005; % XBeach has standing water of <0.005 m on all 
                     % surfaces, even if not flooded, so we need a
                     % threshold to indicate flooding
names = {'gcmhist','gcm2100_withoutWL','gcmhist_2'}; % Fieldnames in output structs
names2 = {'gcmhist_0p5','gcmhist_1','gcmhist_1p5'};
clusters1 = {histcluster,cluster2100,histcluster}; % We need these values to normalize the output results by proportion of year accounted for by each cluster
clusters2 = {histcluster,histcluster,histcluster};
tide = 't0'; % Consider the MSL runs
tide2 = 't0';

% Load data
load(bathyfile);
bathyx = bathy.x - bathy.x(reefCrestInd); % Normalize transect by distance from reef crest

modelOut = load(fileName);
modelOut = modelOut.modelOut;
modelOut2 = load(fileName2);
modelOut2 = modelOut2.modelOut;

% Initialize output variables
absmeantaubx = struct();
meantaubx = struct();
inundation = struct();
meanHs = struct();
meanEnergy = struct();
meanCorr = struct();
meanHsig = struct();
meanU = struct();
meanStir = struct();
meanSetup = struct();
for x = 1 : length(names)
    fields = fieldnames(modelOut.(names{x}));
    
    % Build struct structure for each output variable
    meantaubx.(names{x}) = struct();
    meantaubx.(names{x}).taubx = zeros(size(bathy.x));
    meantaubx.(names2{x}) = struct();
    meantaubx.(names2{x}).taubx = zeros(size(bathy.x));
    
    meanHs.(names{x}) = struct();
    meanHs.(names{x}).Hs = zeros(size(bathy.x));
    meanHs.(names2{x}) = struct();
    meanHs.(names2{x}).Hs = zeros(size(bathy.x));
    
    meanHsig.(names{x}) = struct();
    meanHsig.(names{x}).Hsig = zeros(size(bathy.x));
    meanHsig.(names2{x}) = struct();
    meanHsig.(names2{x}).Hsig = zeros(size(bathy.x));
    
    meanU.(names{x}).u = zeros(size(bathy.x));
    meanU.(names2{x}).u = zeros(size(bathy.x));
    for y = 1 : length(fields)
        hh1 = modelOut.(names{x}).(sprintf('c%d',y)).(tide).data(HHIND).value;
        hh2 = modelOut2.(names2{x}).(sprintf('c%d',y)).(tide).data(HHIND).value;
        [mod1Hsall,mod1Hrmsall,mod1Hs,mod1Hsig,mod1Hslf,mod1Tp,mod1Tpig,mod1Tplf,mod1Tdom,mod1Hrms,mod1Hrmsig,mod1Hrmslf] = spectral_wave_artificial(squeeze(hh1)',0,1,0,4096,[]);
        [mod2Hsall,mod2Hrmsall,mod2Hs,mod2Hsig,mod2Hslf,mod2Tp,mod2Tpig,mod2Tplf,mod2Tdom,mod2Hrms,mod2Hrmsig,mod2Hrmslf] = spectral_wave_artificial(squeeze(hh2)',0,1,0,4096,[]);
        meanHsig.(names{x}).Hsig = meanHsig.(names{x}).Hsig + sqrt((mod1Hsig.^2+mod1Hslf.^2))' .* clusters1{x}(y,5);
        meanHsig.(names2{x}).Hsig = meanHsig.(names2{x}).Hsig + sqrt((mod2Hsig.^2+mod2Hslf.^2))' .* clusters2{x}(y,5);
        meanHs.(names{x}).Hs = meanHs.(names{x}).Hs + mod1Hs' .* clusters1{x}(y,5);
        meanHs.(names2{x}).Hs = meanHs.(names2{x}).Hs + mod2Hs' .* clusters2{x}(y,5);

        meantaubx.(names{x}).taubx = meantaubx.(names{x}).taubx + squeeze(nanmean((modelOut.(names{x}).(sprintf('c%d',y)).(tide).data(TAUIND).value),1))' * clusters1{x}(y,5);
        meantaubx.(names2{x}).taubx = meantaubx.(names2{x}).taubx + squeeze(nanmean((modelOut2.(names2{x}).(sprintf('c%d',y)).(tide).data(TAUIND).value),1))' * clusters2{x}(y,5);
        meanU.(names{x}).u = meanU.(names{x}).u + squeeze(nanmean((modelOut.(names{x}).(sprintf('c%d',y)).(tide).data(UIND).value),1))' * clusters1{x}(y,5);
        meanU.(names2{x}).u = meanU.(names2{x}).u + squeeze(nanmean((modelOut2.(names2{x}).(sprintf('c%d',y)).(tide).data(UIND).value),1))' * clusters2{x}(y,5);
    end
end

% Initialize figure
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength*0.75],'Name','Mean wave energy');
fig.Renderer = 'painters';

subplot(4,2,1);
plot(bathyx,meanHs.(names{1}).Hs,'Color',pltColors{1});
hold on;
plot(bathyx,meanHs.(names{2}).Hs,'Color',pltColors{2});
ylabel('Mean H_{s,hf} (m)','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([0,1.75]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = []; ax.XDir = 'reverse';
legend({'Historical','Late century','Bathymetry'},'Location','best','FontSize',figureFontSize,'Orientation','horizontal');

subplot(4,2,2);
plot(bathyx,meanHs.(names{1}).Hs,'Color',pltColors{1});
hold on;
plot(bathyx,meanHs.(names2{1}).Hs,'Color',pltColors{2});
plot(bathyx,meanHs.(names2{2}).Hs,'Color',pltColors{3});
plot(bathyx,meanHs.(names2{3}).Hs,'Color',pltColors{4});
plot(bathyx,meanHs.(names{3}).Hs,'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([0,1.75]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.XDir = 'reverse';

subplot(4,2,3);
plot(bathyx,meanHsig.(names{1}).Hsig,'Color',pltColors{1});
hold on;
plot(bathyx,meanHsig.(names{2}).Hsig,'Color',pltColors{2});
ylabel('Mean H_{s,lf} (m)','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([0,0.6]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = []; ax.XDir = 'reverse';

subplot(4,2,4);
plot(bathyx,meanHsig.(names{1}).Hsig,'Color',pltColors{1});
hold on;
plot(bathyx,meanHsig.(names2{1}).Hsig,'Color',pltColors{2});
plot(bathyx,meanHsig.(names2{2}).Hsig,'Color',pltColors{3});
plot(bathyx,meanHsig.(names2{3}).Hsig,'Color',pltColors{4});
plot(bathyx,meanHsig.(names{3}).Hsig,'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([0,0.6]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.XDir = 'reverse';

subplot(4,2,5);
plot(bathyx,meantaubx.(names{1}).taubx,'Color',pltColors{1});
hold on;
plot(bathyx,meantaubx.(names{2}).taubx,'Color',pltColors{2});
plot([min(bathyx),max(bathyx)],[0,0],'k--');
ylabel('Mean \tau_{b} (Nm^{-2})','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-1.5,1]);
ax = gca; ax.FontSize = figureFontSize;
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.YTickLabel = []; ax.XTickLabel = []; ax.XDir = 'reverse';

subplot(4,2,6);
plot(bathyx,meantaubx.(names{1}).taubx,'Color',pltColors{1});
hold on;
plot(bathyx,meantaubx.(names2{1}).taubx,'Color',pltColors{2});
plot(bathyx,meantaubx.(names2{2}).taubx,'Color',pltColors{3});
plot(bathyx,meantaubx.(names2{3}).taubx,'Color',pltColors{4});
plot(bathyx,meantaubx.(names{3}).taubx,'Color',pltColors{5});
plot([min(bathyx),max(bathyx)],[0,0],'k--');
xlim(bathyx(bounds));
ylim([-1.5,1]);
ax = gca; ax.FontSize = figureFontSize; ax.YTickLabel = []; 
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XDir = 'reverse'; ax.XTickLabel = [];

subplot(4,2,7);
plot(bathyx,meanU.(names{1}).u,'Color',pltColors{1});
hold on;
plot(bathyx,meanU.(names{2}).u,'Color',pltColors{2});
plot([min(bathyx),max(bathyx)],[0,0],'k--');
ylabel('Mean u (m s^{-1})','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-0.18,0]);
ax = gca; ax.FontSize = figureFontSize;
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.YTickLabel = []; ax.XDir = 'reverse';

subplot(4,2,8);
plot(bathyx,meanU.(names{1}).u,'Color',pltColors{1});
hold on;
plot(bathyx,meanU.(names2{1}).u,'Color',pltColors{2});
plot(bathyx,meanU.(names2{2}).u,'Color',pltColors{3});
plot(bathyx,meanU.(names2{3}).u,'Color',pltColors{4});
plot(bathyx,meanU.(names{3}).u,'Color',pltColors{5});
plot([min(bathyx),max(bathyx)],[0,0],'k--');
xlim(bathyx(bounds));
ylim([-0.18,0]);
ax = gca; ax.FontSize = figureFontSize; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XDir = 'reverse';

print(fullfile(figureOutDir,'MeanWavePropagation_alfa_HIRES_SLRscenarios_generic'),'-depsc','-r600');

%..........................................................................

%..........................................................................
%% 4.) Figure 7. Compare nonlinearity in wave velocity time series (i.e. skewness 
% (velocity) and asymmetry (acceleration))
% Skewness and Asymmetry calculated as in Ruessink et al. (2012)
% Frequency decomposition performed using EEMD (Wu and Huang, 2009)

fileName = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File1.mat');
fileName2 = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File2.mat');
bathyfile = fullfile(baseDir,'Kwajalein_transect_model_bathymetry.mat');
bounds = [193,637];
load(bathyfile);
bathyx = bathy.x - bathy.x(reefCrestInd);
clusters1 = {histcluster,cluster2100,histcluster};
clusters2 = {histcluster,histcluster,histcluster};
tide = 't0';

% Load in the first output file
modelOut = load(fileName);
modelOut = modelOut.modelOut;

% Initialize the output and calculate average asymmetry for all motions,
% long waves, and short waves
modelOutSu = struct();
modelOutAu = struct();
fields = fieldnames(modelOut);
for ii = 1 : length(fields)
    modelOutSu.(fields{ii}) = struct();
    clusts = fieldnames(modelOut.(fields{ii}));
    modelOutSu.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutSu.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutSu.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    modelOutAu.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutAu.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutAu.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    
    for jj = 1 : length(clusts)
        uSeries = squeeze(modelOut.(fields{ii}).(clusts{jj}).(tide).data(UIND).value)-repmat(nanmean(squeeze(modelOut.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),size(squeeze(modelOut.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),1);
        hiSeries = zeros(size(uSeries));
        loSeries = zeros(size(uSeries));
        resSeries = zeros(size(uSeries));
        uSeries2 = zeros(size(uSeries));
        hiSeries2 = zeros(size(uSeries));
        loSeries2 = zeros(size(uSeries));
        for kk = 1 : size(uSeries,2)
            % Perform EEMD to separate high- and low-frequency components
            [hiSeries(:,kk),loSeries(:,kk),resSeries(:,kk)] = wave_hht(uSeries(:,kk)',modelSampleFreq);
            uSeries2(:,kk) = imag(hilbert(uSeries(:,kk)));
            hiSeries2(:,kk) = imag(hilbert(hiSeries(:,kk)+resSeries(:,kk)));
            loSeries2(:,kk) = imag(hilbert(loSeries(:,kk)));
        end
        modelOutSu.(fields{ii}).u(jj,:) = nanmean(uSeries.^3,1)./std(uSeries,[],1).^3;
        modelOutSu.(fields{ii}).uhi(jj,:) = nanmean((hiSeries+resSeries).^3,1)./std((hiSeries+resSeries),[],1).^3;
        modelOutSu.(fields{ii}).ulo(jj,:) = nanmean((loSeries).^3,1)./std((loSeries),[],1).^3;
        modelOutAu.(fields{ii}).u(jj,:) = nanmean(uSeries2.^3,1)./std(uSeries2,[],1).^3;
        modelOutAu.(fields{ii}).uhi(jj,:) = nanmean((hiSeries2).^3,1)./std((hiSeries2),[],1).^3;
        modelOutAu.(fields{ii}).ulo(jj,:) = nanmean((loSeries2).^3,1)./std((loSeries2),[],1).^3;
    end
end
clear modelOut;

% Load in the first output file
modelOut2 = load(fileName2);
modelOut2 = modelOut2.modelOut;

% Initialize the output and calculate average asymmetry for all motions,
% long waves, and short waves
modelOutSu2 = struct();
modelOutAu2 = struct();
fields = fieldnames(modelOut2);
for ii = 1 : length(fields)
    modelOutSu2.(fields{ii}) = struct();
    clusts = fieldnames(modelOut2.(fields{ii}));
    modelOutSu2.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutSu2.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutSu2.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    modelOutAu2.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutAu2.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutAu2.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    
    for jj = 1 : length(clusts)
        uSeries = squeeze(modelOut2.(fields{ii}).(clusts{jj}).(tide).data(UIND).value)-repmat(nanmean(squeeze(modelOut2.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),size(squeeze(modelOut2.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),1);
        hiSeries = zeros(size(uSeries));
        loSeries = zeros(size(uSeries));
        resSeries = zeros(size(uSeries));
        uSeries2 = zeros(size(uSeries));
        hiSeries2 = zeros(size(uSeries));
        loSeries2 = zeros(size(uSeries));
        for kk = 1 : size(uSeries,2)
            [hiSeries(:,kk),loSeries(:,kk),resSeries(:,kk)] = wave_hht(uSeries(:,kk)',modelSampleFreq);
            uSeries2(:,kk) = imag(hilbert(uSeries(:,kk)));
            hiSeries2(:,kk) = imag(hilbert(hiSeries(:,kk)+resSeries(:,kk)));
            loSeries2(:,kk) = imag(hilbert(loSeries(:,kk)));
        end
        modelOutSu2.(fields{ii}).u(jj,:) = nanmean(uSeries.^3,1)./std(uSeries,[],1).^3;
        modelOutSu2.(fields{ii}).uhi(jj,:) = nanmean((hiSeries+resSeries).^3,1)./std((hiSeries+resSeries),[],1).^3;
        modelOutSu2.(fields{ii}).ulo(jj,:) = nanmean((loSeries).^3,1)./std((loSeries),[],1).^3;
        modelOutAu2.(fields{ii}).u(jj,:) = nanmean(uSeries2.^3,1)./std(uSeries2,[],1).^3;
        modelOutAu2.(fields{ii}).uhi(jj,:) = nanmean((hiSeries2).^3,1)./std((hiSeries2),[],1).^3;
        modelOutAu2.(fields{ii}).ulo(jj,:) = nanmean((loSeries2).^3,1)./std((loSeries2),[],1).^3;
    end
end
clear modelOut2;

% Remove values that are obviously outliers, mainly on the island surface
fields = fieldnames(modelOutSu);
for ii = 1 : length(fields)
    innerFields = fieldnames(modelOutSu.(fields{ii}));
    for jj = 1 : length(innerFields)
        modelOutSu.(fields{ii}).(innerFields{jj})(modelOutSu.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutSu.(fields{ii}).(innerFields{jj})(modelOutSu.(fields{ii}).(innerFields{jj})>100) = nan;
        modelOutSu2.(fields{ii}).(innerFields{jj})(modelOutSu2.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutSu2.(fields{ii}).(innerFields{jj})(modelOutSu2.(fields{ii}).(innerFields{jj})>100) = nan;
        modelOutAu.(fields{ii}).(innerFields{jj})(modelOutAu.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutAu.(fields{ii}).(innerFields{jj})(modelOutAu.(fields{ii}).(innerFields{jj})>100) = nan;
        modelOutAu2.(fields{ii}).(innerFields{jj})(modelOutAu2.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutAu2.(fields{ii}).(innerFields{jj})(modelOutAu2.(fields{ii}).(innerFields{jj})>100) = nan;
    end
end


% Plot the figure
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength*0.75],'Name','Mean wave skewness');
fig.Renderer = 'painters';

subplot(3,2,1);
plot(bathyx,nansum(modelOutSu.gcmhist.u.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutSu.gcm2100_withoutWL.u.*repmat(clusters1{2}(:,5),1,length(bathyx))),'Color',pltColors{2});
ylabel('S_u','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-0.5,1.3]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = []; ax.XDir = 'reverse';
legend({'Historical','Mid-century','Late century','Bathymetry'},'Location','best','FontSize',figureFontSize,'Orientation','horizontal');

subplot(3,2,2);
plot(bathyx,nansum(modelOutSu.gcmhist.u.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutSu2.gcmhist_0p5.u.*repmat(clusters2{1}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot(bathyx,nansum(modelOutSu2.gcmhist_1.u.*repmat(clusters2{2}(:,5),1,length(bathyx))),'Color',pltColors{3});
plot(bathyx,nansum(modelOutSu2.gcmhist_1p5.u.*repmat(clusters2{3}(:,5),1,length(bathyx))),'Color',pltColors{4});
plot(bathyx,nansum(modelOutSu.gcmhist_2.u.*repmat(clusters1{3}(:,5),1,length(bathyx))),'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([-0.5,1.3]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,3);
plot(bathyx,nansum(modelOutSu.gcmhist.uhi.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutSu.gcm2100_withoutWL.uhi.*repmat(clusters1{2}(:,5),1,length(bathyx))),'Color',pltColors{2});
ylabel('S_u','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-0.5,1.3]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,4);
plot(bathyx,nansum(modelOutSu.gcmhist.uhi.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutSu2.gcmhist_0p5.uhi.*repmat(clusters2{1}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot(bathyx,nansum(modelOutSu2.gcmhist_1.uhi.*repmat(clusters2{2}(:,5),1,length(bathyx))),'Color',pltColors{3});
plot(bathyx,nansum(modelOutSu2.gcmhist_1p5.uhi.*repmat(clusters2{3}(:,5),1,length(bathyx))),'Color',pltColors{4});
plot(bathyx,nansum(modelOutSu.gcmhist_2.uhi.*repmat(clusters1{3}(:,5),1,length(bathyx))),'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([-0.5,1.3]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,5);
plot(bathyx,nansum(modelOutSu.gcmhist.ulo.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutSu.gcm2100_withoutWL.ulo.*repmat(clusters1{2}(:,5),1,length(bathyx))),'Color',pltColors{2});
ylabel('S_u','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-0.5,1.3]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
xlabel('Distance from the reef crest (m)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.YTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,6);
plot(bathyx,nansum(modelOutSu.gcmhist.ulo.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutSu2.gcmhist_0p5.ulo.*repmat(clusters2{1}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot(bathyx,nansum(modelOutSu2.gcmhist_1.ulo.*repmat(clusters2{2}(:,5),1,length(bathyx))),'Color',pltColors{3});
plot(bathyx,nansum(modelOutSu2.gcmhist_1p5.ulo.*repmat(clusters2{3}(:,5),1,length(bathyx))),'Color',pltColors{4});
plot(bathyx,nansum(modelOutSu.gcmhist_2.ulo.*repmat(clusters1{3}(:,5),1,length(bathyx))),'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([-0.5,1.3]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
xlabel('Distance from the reef crest (m)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XDir = 'reverse';

% print(fullfile(figureOutDir,'wave_skewness'),'-depsc','-r600');


% Plot the figure for wave asymmetry
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength*0.75],'Name','Mean wave asymmetry');
fig.Renderer = 'painters';

subplot(3,2,1);
plot(bathyx,nansum(modelOutAu.gcmhist.u.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutAu.gcm2100_withoutWL.u.*repmat(clusters1{2}(:,5),1,length(bathyx))),'Color',pltColors{2});
ylabel('A_u','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-1.3,0.5]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = []; ax.XDir = 'reverse';
legend({'Historical','Mid-century','Late century','Bathymetry'},'Location','best','FontSize',figureFontSize,'Orientation','horizontal');

subplot(3,2,2);
plot(bathyx,nansum(modelOutAu.gcmhist.u.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutAu2.gcmhist_0p5.u.*repmat(clusters2{1}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot(bathyx,nansum(modelOutAu2.gcmhist_1.u.*repmat(clusters2{2}(:,5),1,length(bathyx))),'Color',pltColors{3});
plot(bathyx,nansum(modelOutAu2.gcmhist_1p5.u.*repmat(clusters2{3}(:,5),1,length(bathyx))),'Color',pltColors{4});
plot(bathyx,nansum(modelOutAu.gcmhist_2.u.*repmat(clusters1{3}(:,5),1,length(bathyx))),'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([-1.3,0.5]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,3);
plot(bathyx,nansum(modelOutAu.gcmhist.uhi.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutAu.gcm2100_withoutWL.uhi.*repmat(clusters1{2}(:,5),1,length(bathyx))),'Color',pltColors{2});
ylabel('A_u','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-1.3,0.5]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,4);
plot(bathyx,nansum(modelOutAu.gcmhist.uhi.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutAu2.gcmhist_0p5.uhi.*repmat(clusters2{1}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot(bathyx,nansum(modelOutAu2.gcmhist_1.uhi.*repmat(clusters2{2}(:,5),1,length(bathyx))),'Color',pltColors{3});
plot(bathyx,nansum(modelOutAu2.gcmhist_1p5.uhi.*repmat(clusters2{3}(:,5),1,length(bathyx))),'Color',pltColors{4});
plot(bathyx,nansum(modelOutAu.gcmhist_2.uhi.*repmat(clusters1{3}(:,5),1,length(bathyx))),'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([-1.3,0.5]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,5);
plot(bathyx,nansum(modelOutAu.gcmhist.ulo.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutAu.gcm2100_withoutWL.ulo.*repmat(clusters1{2}(:,5),1,length(bathyx))),'Color',pltColors{2});
ylabel('A_u','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-1.3,0.5]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
xlabel('Distance from the reef crest (m)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.YTickLabel = []; ax.XDir = 'reverse';

subplot(3,2,6);
plot(bathyx,nansum(modelOutAu.gcmhist.ulo.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutAu2.gcmhist_0p5.ulo.*repmat(clusters2{1}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot(bathyx,nansum(modelOutAu2.gcmhist_1.ulo.*repmat(clusters2{2}(:,5),1,length(bathyx))),'Color',pltColors{3});
plot(bathyx,nansum(modelOutAu2.gcmhist_1p5.ulo.*repmat(clusters2{3}(:,5),1,length(bathyx))),'Color',pltColors{4});
plot(bathyx,nansum(modelOutAu.gcmhist_2.ulo.*repmat(clusters1{3}(:,5),1,length(bathyx))),'Color',pltColors{5});
xlim(bathyx(bounds));
ylim([-1.3,0.5]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
xlabel('Distance from the reef crest (m)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XDir = 'reverse';

% print(fullfile(figureOutDir,'wave_asymmetry'),'-depsc','-r600');

%..........................................................................

%..........................................................................
%% 5.) Figure 9. Explanatory figure for potential sediment transport
% Bedload transport calculated as in Soulsby and Damgaard (2005)

% Parameters for shear stress distribution
mu = -0.1;
sigma = 2;
skew = 0.5;
kurt = 3;
% Physical coefficients for calculating bedload transport
baseThresh = 0.06; % Critical Shields parameter
A2 = 12; % Coefficient from Soulsby and Damgaard (2005)
criticalShear = baseThresh .* G .* (rhoS - rho) .* 0.002; % Critical shear stress


% Create a skewed distribution of shear stress values
tauDist = pearsrnd(mu,sigma,skew,kurt,10000,1);

% Create a histogram for raw shear stress
[counts,edges] = histcounts(tauDist,50);
edgeInterval = nanmean(diff(edges));
binCenters = (edges(2:end) + edges(1:end-1))/2;
countsPDF = counts ./ sum(counts) ./ edgeInterval;

[~,lind] = min(abs(binCenters+criticalShear));
[~,rind] = min(abs(binCenters-criticalShear));
criticalShearPDFL = countsPDF(lind);
criticalShearPDFR = countsPDF(rind);
xplotL = [binCenters(binCenters <= -criticalShear),-criticalShear,-criticalShear,fliplr(binCenters(binCenters <= -criticalShear))];
yplotL = [countsPDF(binCenters <= -criticalShear),max(countsPDF(binCenters <= -criticalShear)),zeros(1,length(binCenters(binCenters <= -criticalShear))+1)];
xplotR = [criticalShear,binCenters(binCenters >= criticalShear),fliplr(binCenters(binCenters >= criticalShear)),criticalShear];
yplotR = [countsPDF(find(binCenters >= criticalShear,1)),countsPDF(binCenters >= criticalShear),zeros(1,length(binCenters(binCenters >= criticalShear))+1)];

% Create a histogram for bedload transport and excess shear
% Excess shear
excessShearDist = sign(tauDist) .* max(0,abs(tauDist)-criticalShear);
[countse,edgese] = histcounts(excessShearDist(excessShearDist~=0),50);
edgeIntervale = nanmean(diff(edgese));
binCenterse = (edgese(2:end) + edgese(1:end-1))/2;
countsPDFe = countse ./ sum(countse) ./ edgeIntervale;

% Bedload transport
bedloadDist = A2.*sign(tauDist).*sqrt(abs(tauDist)).*max(0,abs(tauDist)-criticalShear)./(G.*rho.^(1.5).*(rhoS./rho-1));
[countsb,edgesb] = histcounts(bedloadDist(bedloadDist~=0),50);
edgeIntervalb = nanmean(diff(edgesb));
binCentersb = (edgesb(2:end) + edgesb(1:end-1))/2;
countsPDFb = countsb ./ sum(countsb) ./ edgeIntervalb;


% Plot the distributions
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,minWidth,minLength*1.5],'Name','ShearToBedload');
fig.Renderer = 'painters';

subplot(3,1,3);
hold on;
fill(xplotL,yplotL,[1,0,0],'EdgeColor','none');
fill(xplotR,yplotR,[1,0,0],'EdgeColor','none');
plot(binCenters(binCenters <= -criticalShear),countsPDF(binCenters <= -criticalShear),'k-');
plot(binCenters(binCenters >= criticalShear),countsPDF(binCenters >= criticalShear),'k-');
plot(binCenters(abs(binCenters)<criticalShear),countsPDF(abs(binCenters)<criticalShear),'k--');
plot([-criticalShear,-criticalShear],[0,criticalShearPDFL],'r-');
plot([criticalShear,criticalShear],[0,criticalShearPDFR],'r-');

xlim([-6,9]);
ax = gca; ax.FontSize = figureFontSize;
xlabel('Bottom shear stress, \tau_b (Nm^{-2})','FontSize',figureFontSize);
ylabel('PDF','FontSize',figureFontSize);

subplot(3,1,2);
hold on;
plot(binCenterse(binCenterse<0),countsPDFe(binCenterse<0),'r-');
plot(binCenterse(binCenterse>0),countsPDFe(binCenterse>0),'r-');
xlim([-6,9]);
ax = gca; ax.XAxisLocation = 'top'; ax.FontSize = figureFontSize;
xlabel('Excess bottom shear stress, MAX(0,|\tau_b| - \tau_c)|\tau_b|/\tau_b (m^3s^{-1})','FontSize',figureFontSize);
ylabel('PDF','FontSize',figureFontSize);

subplot(3,1,1);
hold on;
plot(binCentersb(binCentersb<0),countsPDFb(binCentersb<0),'r-');
plot(binCentersb(binCentersb>0),countsPDFb(binCentersb>0),'r-');
xlim([-3,4.5]*10^-4);
ax = gca; ax.XAxisLocation = 'top'; ax.FontSize = figureFontSize;
xlabel('Bedload transport rate, q_b (m^3s^{-1})','FontSize',figureFontSize);
ylabel('PDF','FontSize',figureFontSize);

% print(fullfile(figureOutDir,'Bedload_explanation'),'-depsc','-r600');
    
%..........................................................................

%..........................................................................
%% 6.) Figure 12. Calculate urep from XBeach using short wave skewness and asymmetry
% Skewness and Asymmetry calculated as in Ruessink et al. (2012)
% Frequency decomposition performed using EEMD (Wu and Huang, 2009)

fileName = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File1.mat');
fileName2 = fullfile(baseDir,'GCM_forecasts_XBEACHX_Results_File2.mat');
bathyfile = fullfile(baseDir,'Kwajalein_transect_model_bathymetry.mat');
bounds = [193,637];
load(bathyfile);
bathyx = bathy.x - bathy.x(reefCrestInd);

clusters1 = {histcluster,cluster2100,histcluster};
clusters2 = {histcluster,histcluster,histcluster};
tide = 't0';
alpha = 0.1; % Coefficient from van Thiel de Vries (2009)
windowSize = 120;   % Moving window (in s) over which to calculate short 
                    % wave characteristics

% Load in the first output file
modelOut = load(fileName);
modelOut = modelOut.modelOut;

% Initialize the output and calculate average asymmetry for all motions,
% long waves, and short waves
modelOutSu = struct();
modelOutAu = struct();
modelOutUrep = struct();
fields = fieldnames(modelOut);
for ii = 1 : length(fields)
    modelOutSu.(fields{ii}) = struct();
    clusts = fieldnames(modelOut.(fields{ii}));
    modelOutUrep.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutSu.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutSu.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutSu.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    modelOutAu.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutAu.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutAu.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    
    for jj = 1 : length(clusts)
        % Just as before, calculate high-frequency skewness and asymmetry
        % within a moving window, then calculate urep as in van Thiel de
        % Vries (2009)
        meanFlow = nanmean(squeeze(modelOut.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1);
        uSeries = squeeze(modelOut.(fields{ii}).(clusts{jj}).(tide).data(UIND).value)-repmat(nanmean(squeeze(modelOut.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),size(squeeze(modelOut.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),1);
        hiSeries = zeros(size(uSeries));
        loSeries = zeros(size(uSeries));
        resSeries = zeros(size(uSeries));
        uSeries2 = zeros(size(uSeries));
        hiSeries2 = zeros(size(uSeries));
        loSeries2 = zeros(size(uSeries));
        urep = zeros(size(uSeries));
        for kk = 1 : size(uSeries,2)
            [hiSeries(:,kk),loSeries(:,kk),resSeries(:,kk)] = wave_hht(uSeries(:,kk)',modelSampleFreq,{{'NumEnsemble',100},{'freqBound',0.04}});
            uSeries2(:,kk) = imag(hilbert(uSeries(:,kk)));
            hiSeries2(:,kk) = imag(hilbert(hiSeries(:,kk)+resSeries(:,kk)));
            loSeries2(:,kk) = imag(hilbert(loSeries(:,kk)));
            Sk = zeros(size(hiSeries,1),1);
            As = zeros(size(hiSeries,1),1);
            urms = zeros(size(hiSeries,1),1);
            for ll = 1 : size(uSeries,1)
                Sk(ll) = nanmean((hiSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)+resSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)).^3)./std((hiSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)+resSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk))).^3;
                As(ll) = nanmean(hiSeries2(max(1,ll-windowSize/2):min(size(hiSeries2,1),ll+windowSize/2),kk).^3)./std(hiSeries2(max(1,ll-windowSize/2):min(size(hiSeries2,1),ll+windowSize/2),kk)).^3;
                urms(ll) = sqrt(nanmean((hiSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)+resSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)).^2));
            end
            modelOutUrep.(fields{ii}).u(jj,kk) = nanmean(meanFlow(kk)+loSeries(:,kk)+(alpha*Sk-alpha*As).*urms);
        end
        modelOutSu.(fields{ii}).u(jj,:) = nanmean(uSeries.^3,1)./std(uSeries,[],1).^3;
        modelOutSu.(fields{ii}).uhi(jj,:) = nanmean((hiSeries+resSeries).^3,1)./std((hiSeries+resSeries),[],1).^3;
        modelOutSu.(fields{ii}).ulo(jj,:) = nanmean((loSeries).^3,1)./std((loSeries),[],1).^3;
        modelOutAu.(fields{ii}).u(jj,:) = nanmean(uSeries2.^3,1)./std(uSeries2,[],1).^3;
        modelOutAu.(fields{ii}).uhi(jj,:) = nanmean((hiSeries2).^3,1)./std((hiSeries2),[],1).^3;
        modelOutAu.(fields{ii}).ulo(jj,:) = nanmean((loSeries2).^3,1)./std((loSeries2),[],1).^3;
        
        
    end
end
clear modelOut;

% Load in the first output file
modelOut2 = load(fileName2);
modelOut2 = modelOut2.modelOut;

% Initialize the output and calculate average asymmetry for all motions,
% long waves, and short waves
modelOutSu2 = struct();
modelOutAu2 = struct();
modelOutUrep2 = struct();
fields = fieldnames(modelOut2);
for ii = 1 : length(fields)
    modelOutSu2.(fields{ii}) = struct();
    clusts = fieldnames(modelOut2.(fields{ii}));
    modelOutUrep2.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutSu2.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutSu2.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutSu2.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    modelOutAu2.(fields{ii}).u = zeros(length(clusts),length(bathyx));
    modelOutAu2.(fields{ii}).uhi = zeros(length(clusts),length(bathyx));
    modelOutAu2.(fields{ii}).ulo = zeros(length(clusts),length(bathyx));
    
    for jj = 1 : length(clusts)
        meanFlow = nanmean(squeeze(modelOut2.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1);
        uSeries = squeeze(modelOut2.(fields{ii}).(clusts{jj}).(tide).data(UIND).value)-repmat(nanmean(squeeze(modelOut2.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),size(squeeze(modelOut2.(fields{ii}).(clusts{jj}).(tide).data(UIND).value),1),1);
        hiSeries = zeros(size(uSeries));
        loSeries = zeros(size(uSeries));
        resSeries = zeros(size(uSeries));
        uSeries2 = zeros(size(uSeries));
        hiSeries2 = zeros(size(uSeries));
        loSeries2 = zeros(size(uSeries));
        for kk = 1 : size(uSeries,2)
            [hiSeries(:,kk),loSeries(:,kk),resSeries(:,kk)] = wave_hht(uSeries(:,kk)',modelSampleFreq,{'NumEnsemble',100});
            uSeries2(:,kk) = imag(hilbert(uSeries(:,kk)));
            hiSeries2(:,kk) = imag(hilbert(hiSeries(:,kk)+resSeries(:,kk)));
            loSeries2(:,kk) = imag(hilbert(loSeries(:,kk)));
            Sk = zeros(size(hiSeries,1),1);
            As = zeros(size(hiSeries,1),1);
            urms = zeros(size(hiSeries,1),1);
            for ll = 1 : size(uSeries,1)
                Sk(ll) = nanmean((hiSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)+resSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)).^3)./std((hiSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)+resSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk))).^3;
                As(ll) = nanmean(hiSeries2(max(1,ll-windowSize/2):min(size(hiSeries2,1),ll+windowSize/2),kk).^3)./std(hiSeries2(max(1,ll-windowSize/2):min(size(hiSeries2,1),ll+windowSize/2),kk)).^3;
                urms(ll) = sqrt(nanmean((hiSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)+resSeries(max(1,ll-windowSize/2):min(size(hiSeries,1),ll+windowSize/2),kk)).^2));
            end
            modelOutUrep2.(fields{ii}).u(jj,kk) = nanmean(meanFlow(kk)+loSeries(:,kk)+(alpha*Sk-alpha*As).*urms);
        end
        modelOutSu2.(fields{ii}).u(jj,:) = nanmean(uSeries.^3,1)./std(uSeries,[],1).^3;
        modelOutSu2.(fields{ii}).uhi(jj,:) = nanmean((hiSeries+resSeries).^3,1)./std((hiSeries+resSeries),[],1).^3;
        modelOutSu2.(fields{ii}).ulo(jj,:) = nanmean((loSeries).^3,1)./std((loSeries),[],1).^3;
        modelOutAu2.(fields{ii}).u(jj,:) = nanmean(uSeries2.^3,1)./std(uSeries2,[],1).^3;
        modelOutAu2.(fields{ii}).uhi(jj,:) = nanmean((hiSeries2).^3,1)./std((hiSeries2),[],1).^3;
        modelOutAu2.(fields{ii}).ulo(jj,:) = nanmean((loSeries2).^3,1)./std((loSeries2),[],1).^3;
    end
end
clear modelOut2;
        
% Remove outlier values of Skewness and asymmetry, generally on island
% surface
fields = fieldnames(modelOutSu);
for ii = 1 : length(fields)
    innerFields = fieldnames(modelOutSu.(fields{ii}));
    for jj = 1 : length(innerFields)
        modelOutSu.(fields{ii}).(innerFields{jj})(modelOutSu.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutSu.(fields{ii}).(innerFields{jj})(modelOutSu.(fields{ii}).(innerFields{jj})>100) = nan;
        modelOutSu2.(fields{ii}).(innerFields{jj})(modelOutSu2.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutSu2.(fields{ii}).(innerFields{jj})(modelOutSu2.(fields{ii}).(innerFields{jj})>100) = nan;
        modelOutAu.(fields{ii}).(innerFields{jj})(modelOutAu.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutAu.(fields{ii}).(innerFields{jj})(modelOutAu.(fields{ii}).(innerFields{jj})>100) = nan;
        modelOutAu2.(fields{ii}).(innerFields{jj})(modelOutAu2.(fields{ii}).(innerFields{jj})<-100) = nan;
        modelOutAu2.(fields{ii}).(innerFields{jj})(modelOutAu2.(fields{ii}).(innerFields{jj})>100) = nan;
    end
end

% Plot the figure
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength*0.25],'Name','Mean wave asymmetry');
fig.Renderer = 'painters';

subplot(1,2,1);
plot(bathyx,nansum(modelOutUrep.gcmhist.u.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutUrep.gcm2100_withoutWL.u.*repmat(clusters1{2}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot([min(bathyx),max(bathyx)],[0,0],'k--');
ylabel('U_{rep}','FontSize',figureFontSize);
xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
xlim(bathyx(bounds));
ylim([-0.08,0.03]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = []; ax.XDir = 'reverse';

subplot(1,2,2);
plot(bathyx,nansum(modelOutUrep.gcmhist.u.*repmat(clusters1{1}(:,5),1,length(bathyx))),'Color',pltColors{1});
hold on;
plot(bathyx,nansum(modelOutUrep2.gcmhist_0p5.u.*repmat(clusters2{1}(:,5),1,length(bathyx))),'Color',pltColors{2});
plot(bathyx,nansum(modelOutUrep2.gcmhist_1.u.*repmat(clusters2{2}(:,5),1,length(bathyx))),'Color',pltColors{3});
plot(bathyx,nansum(modelOutUrep2.gcmhist_1p5.u.*repmat(clusters2{3}(:,5),1,length(bathyx))),'Color',pltColors{4});
plot(bathyx,nansum(modelOutUrep.gcmhist_2.u.*repmat(clusters1{3}(:,5),1,length(bathyx))),'Color',pltColors{5});
plot([min(bathyx),max(bathyx)],[0,0],'k--');
xlim(bathyx(bounds));
xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
ylim([-0.08,0.03]);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.YTickLabel = [];
yyaxis right;
plot(bathyx,bathy.z,'k-');
ylim([-5,5]);
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
ax = gca; ax.FontSize = figureFontSize; ax.XTickLabel = []; ax.XDir = 'reverse';

% print(fullfile(figureOutDir,'RepresentativeTransport'),'-depsc','-r600');
%..........................................................................

%..........................................................................
%% 7.) Figures 10, 11, SI Potential sediment transport rate, figure with 
% all scenarios
% Bedload transport calculated as in Nielsen and Callaghan (2003), assuming
% fully turbulent flow
fileName = fullfile(baseDir,'Modeled_velocity_10hz.mat');
bathyfile = fullfile(baseDir,'Kwajalein_transect_model_bathymetry.mat');
bounds = [49,160]*4-3;
load(bathyfile);
bathyx = bathy.x - bathy.x(reefCrestInd);
posSlope = atan((bathy.z(distInds+1)-bathy.z(distInds-1))./((bathy.x(distInds+1)-bathy.x(distInds-1))));
negSlope = -posSlope;
names = {'gcmhist','gcm2100_withoutWL','gcmhist_2'}; % Fieldnames in output structs
names2 = {'gcmhist_0p5','gcmhist_1','gcmhist_1p5'};
modelSampleFreq = 10;
TS_LENGTH = 40961;
stress = 1;
distcf = mode(distcf) * ones(size(distcf));
phiTaus = [0,10,20,30,40,50,60]*pi/180;

clusters1 = {histcluster,cluster2100,histcluster};
clusters2 = {histcluster,histcluster,histcluster};
baseThresh = 0.05;

% Two different ways to account for bed slope. Either correct the critical
% Shields parameter, or use a constant applied to resultant flux. Turns out
% the results are quite similar either way.
posBaseThresh = baseThresh * sin(phiS+posSlope)/sin(phiS); %Use slope correction of Damgaard et al., 1997
negBaseThresh = baseThresh * sin(phiS+negSlope)/sin(phiS); %Use slope correction of Damgaard et al., 1997
posFluxCoef = 1;%cos(negSlope).*(1+tan(negSlope)./tan(phiS));
negFluxCoef = 1;%cos(posSlope).*(1+tan(posSlope)./tan(phiS));

A2 = 12;
grainSizes = [62,500,2000,16000,64000]/1000000; % D50 in meters
criticalShear = baseThresh .* G .* (rhoS - rho) .* grainSizes;

% Load in the first output file
load(fileName);

% Establish cumulative probability of wave clusters
cumClusters1 = clusters1; cumClusters1{1} = max(ones(size(clusters1{1}(:,5))),round(clusters1{1}(:,5)*1000)); cumClusters1{2} = max(ones(size(clusters1{2}(:,5))),round(clusters1{2}(:,5)*1000)); cumClusters1{3} = max(ones(size(clusters1{3}(:,5))),round(clusters1{3}(:,5)*1000));
cumClusters2 = clusters2; cumClusters2{1} = max(ones(size(clusters2{1}(:,5))),round(clusters2{1}(:,5)*1000)); cumClusters2{2} = max(ones(size(clusters2{2}(:,5))),round(clusters2{2}(:,5)*1000)); cumClusters2{3} = max(ones(size(clusters2{3}(:,5))),round(clusters2{3}(:,5)*1000));

% Now sample the time-series to create distributions of shear stress
fields = names;
sensitivity_results_tau = zeros(length(fields)*2,length(distInds),length(grainSizes),length(phiTaus));
sensitivity_results_qb = zeros(length(fields)*2,length(distInds),length(grainSizes),length(phiTaus));
for ll = 1 : length(phiTaus)
for ii = 1 : length(fields)
    
    modelOutPts.(fields{ii}).dists = struct();
    modelOutPts.(fields{ii}).dists.raw = cell(length(grainSizes),1);
    modelOutPts.(fields{ii}).dists.rawtau = cell(length(grainSizes),1);
    modelOutPts.(fields{ii}).dists.means = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.meansratio = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.taumeans = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.absmeans = zeros(length(distInds),length(grainSizes));
    for kk = 1 : length(grainSizes)
        modelOutPts.(fields{ii}).dists.raw{kk} = [];
    end
    for jj = 1 : length(cumClusters1{ii})
        extractRaw = zeros(TS_LENGTH,length(distInds));
        extractRawDepth = zeros(TS_LENGTH,length(distInds));
        peakFreq = zeros(1,length(distInds));
        uVar = zeros(1,length(distInds));
        alpha = zeros(1,length(distInds));
        for kk = 1 : length(distInds)
            uRaw = modelOutPts.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
            extractRawDepth(:,kk) = modelOutPts.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).rawhh;
            [mod1Hsall,mod1Hrmsall,mod1Hs,mod1Hsig,mod1Hslf,mod1Tp,mod1Tpig,mod1Tplf,mod1Tdom,mod1Hrms,mod1Hrmsig,mod1Hrmslf] = spectral_wave_artificial(squeeze(uRaw)',0,modelSampleFreq,0,4096,[]);
            peakFreq(kk) = 2 * pi / mod1Tp;
            uVar(kk) = var(uRaw);
            extractRaw(:,kk) = modelOutPts.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
            alpha(kk) = nanmean(abs(uRaw)) / (nanmean(abs(uRaw)) + sqrt(2*nanmean((uRaw-nanmean(uRaw)).^2)));
        end
        acc = [];
        for kk = 1 : TS_LENGTH : size(extractRaw,1)
            extractTemp = cat(1,extractRaw(kk,:),extractRaw(kk:min(kk+TS_LENGTH-1,size(extractRaw,1)),:),extractRaw(min(kk+TS_LENGTH-1,size(extractRaw,1)),:));
            acc = cat(1,acc,(extractTemp(3:end,:)-extractTemp(1:end-2,:))./(2/modelSampleFreq*repmat(peakFreq,TS_LENGTH,1)));
        end
        
        for kk = 1 : length(grainSizes)
            Arms = sqrt(2)./peakFreq .* sqrt(uVar);
            f25 = exp(5.5.*(2.5.*grainSizes(kk)./Arms).^0.2-6.3);
            fc = min((0.4./(1+log(2.5.*grainSizes(kk)./30./extractRawDepth))).^2,0.5);
            switch stress
                case 0
                    f = 2*distcf;
                case 1
                    f = f25;
                case 2
                    f = alpha .* fc + (1-alpha) .* f25;
            end
            ustar = sqrt(0.5.*f).*(cos(phiTaus(ll)).*extractRaw + sin(phiTaus(ll)).*acc);
            tau = rho.*sign(ustar).*ustar.^2;
            theta = tau ./ rho ./ ((rhoS/rho-1)*G*grainSizes(kk));
            baseThresh = (ustar<0).*repmat(negBaseThresh,size(extractRaw,1),1) + (ustar>=0).*repmat(posBaseThresh,size(extractRaw,1),1);
            fluxCoef = (ustar<0).*repmat(negFluxCoef,size(extractRaw,1),1) + (ustar>=0).*repmat(posFluxCoef,size(extractRaw,1),1);
            adjCriticalShear = baseThresh .* G .* (rhoS - rho) .* grainSizes(kk);
            modelOutPts.(fields{ii}).dists.rawtau{kk} = cat(1,modelOutPts.(fields{ii}).dists.rawtau{kk},sign(tau).*max(abs(tau)-adjCriticalShear,0).*clusters1{ii}(jj,end)/TS_LENGTH);
            modelOutPts.(fields{ii}).dists.raw{kk} = cat(1,modelOutPts.(fields{ii}).dists.raw{kk},fluxCoef.*A2.*sign(tau).*sqrt(abs(theta)).*max(abs(theta)-baseThresh,0).*sqrt(G.*(rhoS./rho-1).*grainSizes(kk)^3).*clusters1{ii}(jj,end)/TS_LENGTH);
        end
    end
    for jj = 1 : length(distInds)
        for kk = 1 : length(grainSizes)
            modelOutPts.(fields{ii}).dists.means(jj,kk) = nansum(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj));
            if ~isempty(nansum(abs(modelOutPts.(fields{ii}).dists.raw{kk}(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)<0,jj))))
                modelOutPts.(fields{ii}).dists.meansratio(jj,kk) = nansum(modelOutPts.(fields{ii}).dists.raw{kk}(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)>0,jj))/nansum(abs(modelOutPts.(fields{ii}).dists.raw{kk}(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)<0,jj)));
            end
            modelOutPts.(fields{ii}).dists.taumeans(jj,kk) = nansum(modelOutPts.(fields{ii}).dists.rawtau{kk}(:,jj));
            modelOutPts.(fields{ii}).dists.absmeans(jj,kk) = nansum(abs(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)));
        end
    end
    sensitivity_results_tau(ii,:,:,ll) = modelOutPts.(fields{ii}).dists.taumeans;
    sensitivity_results_qb(ii,:,:,ll) = modelOutPts.(fields{ii}).dists.means;
end
end

fields = names2;
for ll = 1 : length(phiTaus)
for ii = 1 : length(fields)
    
    modelOutPts2.(fields{ii}).dists = struct();
    modelOutPts2.(fields{ii}).dists.raw = cell(length(grainSizes),1);
    modelOutPts2.(fields{ii}).dists.rawtau = cell(length(grainSizes),1);
    modelOutPts2.(fields{ii}).dists.means = zeros(length(distInds),length(grainSizes));
    modelOutPts2.(fields{ii}).dists.meansratio = zeros(length(distInds),length(grainSizes));
    modelOutPts2.(fields{ii}).dists.absmeans = zeros(length(distInds),length(grainSizes));
    for kk = 1 : length(grainSizes)
        modelOutPts2.(fields{ii}).dists.raw{kk} = [];
    end
    for jj = 1 : length(cumClusters2{ii})
        extractRaw = zeros(TS_LENGTH,length(distInds));
        extractRawDepth = zeros(TS_LENGTH,length(distInds));
        peakFreq = zeros(1,length(distInds));
        uVar = zeros(1,length(distInds));
        alpha = zeros(1,length(distInds));
        for kk = 1 : length(distInds)
            extractRawDepth(:,kk) = modelOutPts2.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).rawhh;
            uRaw = modelOutPts2.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
            [mod1Hsall,mod1Hrmsall,mod1Hs,mod1Hsig,mod1Hslf,mod1Tp,mod1Tpig,mod1Tplf,mod1Tdom,mod1Hrms,mod1Hrmsig,mod1Hrmslf] = spectral_wave_artificial(squeeze(uRaw)',0,modelSampleFreq,0,4096,[]);
            peakFreq(kk) = 2 * pi / mod1Tp;
            uVar(kk) = var(uRaw);
            extractRaw(:,kk) = modelOutPts2.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
            alpha(kk) = nanmean(abs(uRaw)) / (nanmean(abs(uRaw)) + sqrt(2*nanmean((uRaw-nanmean(uRaw)).^2)));
        end
        acc = [];
        for kk = 1 : TS_LENGTH : size(extractRaw,1)
            extractTemp = cat(1,extractRaw(kk,:),extractRaw(kk:min(kk+TS_LENGTH-1,size(extractRaw,1)),:),extractRaw(min(kk+TS_LENGTH-1,size(extractRaw,1)),:));
            acc = cat(1,acc,(extractTemp(3:end,:)-extractTemp(1:end-2,:))./(2/modelSampleFreq*repmat(peakFreq,TS_LENGTH,1)));
        end
        
        for kk = 1 : length(grainSizes)
            Arms = sqrt(2)./peakFreq .* sqrt(uVar);
            f25 = exp(5.5.*(2.5.*grainSizes(kk)./Arms).^0.2-6.3);
            fc = min((0.4./(1+log(2.5.*grainSizes(kk)./30./extractRawDepth))).^2,0.5);
            switch stress
                case 0
                    f = 2*distcf;
                case 1
                    f = f25;
                case 2
                    f = alpha .* fc + (1-alpha) .* f25;
            end
            ustar = sqrt(0.5.*f).*(cos(phiTaus(ll)).*extractRaw + sin(phiTaus(ll)).*acc);
            tau = rho.*sign(ustar).*ustar.^2;
            theta = tau ./ rho ./ ((rhoS/rho-1)*G*grainSizes(kk));
            baseThresh = (extractRaw<0).*repmat(negBaseThresh,size(extractRaw,1),1) + (extractRaw>=0).*repmat(posBaseThresh,size(extractRaw,1),1);
            fluxCoef = (ustar<0).*repmat(negFluxCoef,size(extractRaw,1),1) + (ustar>=0).*repmat(posFluxCoef,size(extractRaw,1),1);
            adjCriticalShear = baseThresh .* G .* (rhoS - rho) .* grainSizes(kk);
            modelOutPts2.(fields{ii}).dists.rawtau{kk} = cat(1,modelOutPts2.(fields{ii}).dists.rawtau{kk},sign(tau).*max(abs(tau)-adjCriticalShear,0).*clusters2{ii}(jj,end)/TS_LENGTH);
            modelOutPts2.(fields{ii}).dists.raw{kk} = cat(1,modelOutPts2.(fields{ii}).dists.raw{kk},fluxCoef*A2.*sign(tau).*sqrt(abs(theta)).*max(abs(theta)-baseThresh,0).*sqrt(G.*(rhoS./rho-1).*grainSizes(kk)^3).*clusters2{ii}(jj,end)/TS_LENGTH);
        end
    end
    for jj = 1 : length(distInds)
        for kk = 1 : length(grainSizes)
            modelOutPts2.(fields{ii}).dists.means(jj,kk) = nansum(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj));
            if ~isempty(nansum(abs(modelOutPts2.(fields{ii}).dists.raw{kk}(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)<0,jj))))
                modelOutPts2.(fields{ii}).dists.meansratio(jj,kk) = nansum(modelOutPts2.(fields{ii}).dists.raw{kk}(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)>0,jj))/nansum(abs(modelOutPts2.(fields{ii}).dists.raw{kk}(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)<0,jj)));
            end
            modelOutPts2.(fields{ii}).dists.taumeans(jj,kk) = nansum(modelOutPts2.(fields{ii}).dists.rawtau{kk}(:,jj));
            modelOutPts2.(fields{ii}).dists.absmeans(jj,kk) = nansum(abs(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)));
        end
    end
    sensitivity_results_tau(ii+length(fields),:,:,ll) = modelOutPts2.(fields{ii}).dists.taumeans;
    sensitivity_results_qb(ii+length(fields),:,:,ll) = modelOutPts2.(fields{ii}).dists.means;
end
end

ybounds = {[-8,8],[-1.2,1.2],[-1.2,1.2],[-1.2,1.2],[-1.2,1.2], [-100,13000],[-100,410],[-100,410],[-100,410],[-100,410], [-0.6,2.1],[-0.6,2.1],[-0.6,2.1],[-0.6,2.1],[-0.6,2.1], [-100,12000],[-100,425],[-100,425],[-100,425],[-100,425]};
ytick = cell(1,length(ybounds));
for ii = 1 : length(ybounds)
    ytick{ii} = sort(cat(2,ybounds{ii},[0]));
end

% Figures 10, 11. Plot bedload potential transport across grain sizes at
% separate points across the fringing reef, and sensitivity to tau
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength*1.5],'Name','Distributions');
fig.Renderer = 'painters';

% First plot model bathymetry profile
subplot(5,5,2:4);
plot(bathyx,bathy.z,'k-');
hold on;
scatter(bathyx(distInds),bathy.z(distInds),'g^','MarkerFaceColor','g');
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
ax = gca; ax.XDir = 'reverse';

% Now plot one row of excess shear stress for all model scenarios
for ii = 1 : length(distInds)
    subplot(5,5,5+length(distInds)-(ii-1));
    hold on;
    for jj = [1,2,4,5,6,3]
        plot(log(grainSizes*1000),squeeze(sensitivity_results_tau(jj,ii,:,3)),'Color',pltColors{jj},'LineStyle',pltDashes{jj});
    end

    plot([-3,5],[0,0],'k--');
    xlim([-3,5]);
    xticks(log(grainSizes*1000));
    xticklabels({'1','2','3','4','5'});
%     ylim(ybounds{ii}); yticks(ytick{ii});
    ax = gca;
    ax.YScale = 'log';
    ylim([0.001,10]);
    yticks([10^-2,10^-1,10^0,10^1,10^2]);
    if ii == length(distInds)
        ylabel('(\tau - \tau_c) (Nm^{-2})','FontSize',figureFontSize);
    end
end

% Now plot potential bedload transport
for ii = 1 : length(distInds)
    subplot(5,5,10+length(distInds)-(ii-1));
    hold on;
    for jj = [1,2,4,5,6,3]
        plot(log(grainSizes*1000),squeeze(sensitivity_results_qb(jj,ii,:,3))*1e4,'Color',pltColors{jj},'LineStyle',pltDashes{jj});
    end
    plot([-3,5],[0,0],'k--');
    xlim([-3,5]);
    xticks(log(grainSizes*1000));
    xticklabels({'1','2','3','4','5'});
%     ylim(ybounds{ii+length(distInds)}); yticks(ytick{ii+length(distInds)});
    ax = gca;
    ax.YScale = 'log';
    ylim([10^-2,100]);
    yticks([10^-1,10^0,10^1,10^2,10^3]);
    if ii == length(distInds)
        ylabel('Mean q (cm^2s^{-1})','FontSize',figureFontSize);
    end
end

% Now plot changes in tau and qb as functions of SLR and phiTau
for ii = 1 : length(distInds)
    subplot(5,5,15+length(distInds)-(ii-1));
    for jj = 1 : 2 : length(phiTaus)
        hold on;
        plot([0,0.5,1.0,1.5,2.0],sensitivity_results_tau([1,4,5,6,3],ii,3,jj),'-');
%         plot([0,2],[0,0],'k--');
        xlim([0,2]);
        ax = gca;
        ylim([min(min(sensitivity_results_tau([1,4,5,6,3],ii,3,:))),max(max(sensitivity_results_tau([1,4,5,6,3],ii,3,:)))]);
        if ii == length(distInds)
            ylabel('Mean \tau_e (Nm^{-2})','FontSize',figureFontSize);
        end
    end
end

% Now plot changes in tau and qb as functions of SLR and phiTau
for ii = 1 : length(distInds)
    subplot(5,5,20+length(distInds)-(ii-1));
    for jj = 1 : 2 : length(phiTaus)
        hold on;
        plot([0,0.5,1.0,1.5,2.0],sensitivity_results_qb([1,4,5,6,3],ii,3,jj)*1e4,'-');
%         plot([0,2],[0,0],'k--');
        xlim([0,2]);
        ax = gca;
        ylim([min(min(sensitivity_results_qb([1,4,5,6,3],ii,3,:))),max(max(sensitivity_results_qb([1,4,5,6,3],ii,3,:)))]*1e4);
        if ii == length(distInds)
            ylabel('Mean q (cm^2s^{-1})','FontSize',figureFontSize);
        end
    end
end


% print(fullfile(figureOutDir,'Potential_bedload_transport'),'-depsc','-r600');    


% Figure 13. Now plot potential bedload transport across the fringing reef 
% profile
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength*3/4],'Name','Pot_Transport_Distributions');
fig.Renderer = 'painters';
for kk = 2 : length(grainSizes)
    subplot(2,2,kk-1);
    hold on;
    plot(bathyx(distInds),modelOutPts.(names{1}).dists.means(:,kk)*1e4,'Color',pltColors{1});
    plot(bathyx(distInds),modelOutPts.(names{2}).dists.means(:,kk)*1e4,'Color','m');
    plot(bathyx(distInds),modelOutPts2.(names2{1}).dists.means(:,kk)*1e4,'Color',pltColors{2});
    plot(bathyx(distInds),modelOutPts2.(names2{2}).dists.means(:,kk)*1e4,'Color',pltColors{3});
    plot(bathyx(distInds),modelOutPts2.(names2{3}).dists.means(:,kk)*1e4,'Color',pltColors{4});
    plot(bathyx(distInds),modelOutPts.(names{3}).dists.means(:,kk)*1e4,'Color',pltColors{5});
    plot([-50,125],[0,0],'k--');
    xlim([-50,125]); 
    ylim([0.01,500]);
    title(['D_{50}',sprintf(' = %3.2f mm',grainSizes(kk)*1000)],'FontSize',figureFontSize);
    ax = gca; ax.XDir = 'reverse'; ax.YScale = 'log';
    ylim([0.1,100]); yticks([0.1,1,10,100]);
    if kk > 3
        xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
    else
        ax.XTickLabel = '';
    end
    if mod(kk,2)==0
        ylabel('Mean q (cm^2s^{-1})','FontSize',figureFontSize);
    else
        ax.YTickLabel = '';
    end
    yyaxis right;
    plot(bathyx,bathy.z,'Color',[0.2,0.2,0.2]);
    ylim([-5,5]);
    if mod(kk,2) > 0
        ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
    else
        ax.YTickLabel = '';
    end
end

% print(fullfile(figureOutDir,'Sediment_flux_profiles_Nielsen_Callaghan'),'-depsc','-r600');

%..........................................................................

%..........................................................................
%% 8.) Figure 8. Build shear stress distributions
% Shear stress calculated as in Nielsen and Callaghan, 2003
% 
fileName = fullfile(baseDir,'Modeled_velocity_10hz.mat');
bathyfile = fullfile(baseDir,'Kwajalein_transect_model_bathymetry.mat');
bounds = [49,160]*4-3;
load(bathyfile);
bathyx = bathy.x - bathy.x(reefCrestInd);
posSlope = atan((bathy.z(distInds+1)-bathy.z(distInds-1))./((bathy.x(distInds+1)-bathy.x(distInds-1))));
negSlope = -posSlope;
names = {'gcmhist','gcm2100_withoutWL','gcmhist_2'}; % Fieldnames in output structs
names2 = {'gcmhist_0p5','gcmhist_1','gcmhist_1p5'};
modelSampleFreq = 10;
TS_LENGTH = 40961;
phiTaus = [20]*pi/180;%[0,10,20,30,40,50,60]*pi/180;

clusters1 = {histcluster,cluster2100,histcluster};
clusters2 = {histcluster,histcluster,histcluster};
baseThresh = 0.05;
posBaseThresh = baseThresh * sin(phiS+posSlope)/sin(phiS); %Use slope correction of Damgaard et al., 1997
negBaseThresh = baseThresh * sin(phiS+negSlope)/sin(phiS); %Use slope correction of Damgaard et al., 1997
A2 = 12;
grainSizes = [62,500,2000,16000,64000]/1000000; % D50 in meters
criticalShear = baseThresh .* G .* (rhoS - rho) .* grainSizes;

% Load in the first output file
load(fileName);

% Establish cumulative probability of wave clusters
cumClusters1 = clusters1; cumClusters1{1} = max(ones(size(clusters1{1}(:,5))),round(clusters1{1}(:,5)*1000)); cumClusters1{2} = max(ones(size(clusters1{2}(:,5))),round(clusters1{2}(:,5)*1000)); cumClusters1{3} = max(ones(size(clusters1{3}(:,5))),round(clusters1{3}(:,5)*1000));
cumClusters2 = clusters2; cumClusters2{1} = max(ones(size(clusters2{1}(:,5))),round(clusters2{1}(:,5)*1000)); cumClusters2{2} = max(ones(size(clusters2{2}(:,5))),round(clusters2{2}(:,5)*1000)); cumClusters2{3} = max(ones(size(clusters2{3}(:,5))),round(clusters2{3}(:,5)*1000));

% Histogram edges (different per location, but we'll only use 2 mm size sed)
edges = {-20:0.1:20,...
    -20:0.1:20,...
    -20:0.1:20,...
    -20:0.1:20,...
    -20:0.1:20};

% Now sample the time-series to create distributions of shear stress
fields = names;
for ii = 1 : length(fields)
    
    modelOutPts.(fields{ii}).dists = struct();
    modelOutPts.(fields{ii}).dists.raw = cell(length(grainSizes),1);
    modelOutPts.(fields{ii}).dists.rawtau = cell(length(grainSizes),1);
    modelOutPts.(fields{ii}).dists.postau = cell(length(grainSizes),1);
    modelOutPts.(fields{ii}).dists.means = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.meansratio = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.taumeans = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.absmeans = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.posmeans = zeros(length(distInds),length(grainSizes));
    modelOutPts.(fields{ii}).dists.taudist = cell(length(distInds),1);
    for ll = 1 : length(distInds)
        modelOutPts.(fields{ii}).dists.taudist{ll} = cat(1,edges{ll},zeros(size(edges{ll})));
    end
    for kk = 1 : length(grainSizes)
        modelOutPts.(fields{ii}).dists.raw{kk} = [];
    end
    for jj = 1 : length(cumClusters1{ii})
        extractRaw = zeros(TS_LENGTH,length(distInds));
        peakFreq = zeros(1,length(distInds));
        uVar = zeros(1,length(distInds));
        for kk = 1 : length(distInds)
            uRaw = modelOutPts.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
            [mod1Hsall,mod1Hrmsall,mod1Hs,mod1Hsig,mod1Hslf,mod1Tp,mod1Tpig,mod1Tplf,mod1Tdom,mod1Hrms,mod1Hrmsig,mod1Hrmslf] = spectral_wave_artificial(squeeze(uRaw)',0,modelSampleFreq,0,4096,[]);
            peakFreq(kk) = 2 * pi / mod1Tp;
            uVar(kk) = var(uRaw);
            extractRaw(:,kk) = modelOutPts.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
        end
        acc = [];
        for kk = 1 : TS_LENGTH : size(extractRaw,1)
            extractTemp = cat(1,extractRaw(kk,:),extractRaw(kk:min(kk+TS_LENGTH-1,size(extractRaw,1)),:),extractRaw(min(kk+TS_LENGTH-1,size(extractRaw,1)),:));
            acc = cat(1,acc,(extractTemp(3:end,:)-extractTemp(1:end-2,:))./(2/modelSampleFreq*repmat(peakFreq,TS_LENGTH,1)));
        end
        
        for kk = 1 : length(grainSizes)
            Arms = sqrt(2)./peakFreq .* sqrt(uVar);
            f25 = exp(5.5.*(2.5.*grainSizes(kk)./Arms).^0.2-6.3);
            ustar = sqrt(0.5.*f25).*(cos(phiTaus(1)).*extractRaw + sin(phiTaus(1)).*acc);
            tau = rho.*sign(ustar).*ustar.^2;
            theta = tau ./ rho ./ ((rhoS/rho-1)*G*grainSizes(kk));
            baseThresh = (extractRaw<0).*repmat(negBaseThresh,size(extractRaw,1),1) + (extractRaw>=0).*repmat(posBaseThresh,size(extractRaw,1),1);
            if grainSizes(kk)*1000 == 2
                for ll = 1 : length(distInds)
                    counts = histcounts(tau(:,ll),edges{ll});
                    modelOutPts.(fields{ii}).dists.taudist{ll}(2,:) = modelOutPts.(fields{ii}).dists.taudist{ll}(2,:)+cat(2,counts,0)*cumClusters1{ii}(jj,end);
                end
            end
            modelOutPts.(fields{ii}).dists.rawtau{kk} = cat(1,modelOutPts.(fields{ii}).dists.rawtau{kk},sign(tau).*max(abs(tau)-criticalShear(kk),0).*clusters1{ii}(jj,end)/TS_LENGTH);
            modelOutPts.(fields{ii}).dists.postau{kk} = cat(1,modelOutPts.(fields{ii}).dists.postau{kk},(sign(tau).*max(abs(tau)-criticalShear(kk),0)>0).*clusters1{ii}(jj,end)/TS_LENGTH);
            modelOutPts.(fields{ii}).dists.raw{kk} = cat(1,modelOutPts.(fields{ii}).dists.raw{kk},A2.*sign(tau).*sqrt(abs(theta)).*max(abs(theta)-baseThresh,0).*sqrt(G.*(rhoS./rho-1).*grainSizes(kk)^3).*clusters1{ii}(jj,end)/TS_LENGTH);
        end
    end
    for jj = 1 : length(distInds)
        for kk = 1 : length(grainSizes)
            modelOutPts.(fields{ii}).dists.means(jj,kk) = nansum(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj));
            if ~isempty(nansum(abs(modelOutPts.(fields{ii}).dists.raw{kk}(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)<0,jj))))
                modelOutPts.(fields{ii}).dists.meansratio(jj,kk) = nansum(modelOutPts.(fields{ii}).dists.raw{kk}(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)>0,jj))/nansum(abs(modelOutPts.(fields{ii}).dists.raw{kk}(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)<0,jj)));
            end
            modelOutPts.(fields{ii}).dists.taumeans(jj,kk) = nansum(modelOutPts.(fields{ii}).dists.rawtau{kk}(:,jj));
            modelOutPts.(fields{ii}).dists.absmeans(jj,kk) = nansum(abs(modelOutPts.(fields{ii}).dists.raw{kk}(:,jj)));
            modelOutPts.(fields{ii}).dists.posmeans(jj,kk) = nansum(modelOutPts.(fields{ii}).dists.postau{kk}(:,jj));
        end
    end
    
end


fields = names2;
for ii = 1 : length(fields)
    
    modelOutPts2.(fields{ii}).dists = struct();
    modelOutPts2.(fields{ii}).dists.raw = cell(length(grainSizes),1);
    modelOutPts2.(fields{ii}).dists.rawtau = cell(length(grainSizes),1);
    modelOutPts2.(fields{ii}).dists.postau = cell(length(grainSizes),1);
    modelOutPts2.(fields{ii}).dists.means = zeros(length(distInds),length(grainSizes));
    modelOutPts2.(fields{ii}).dists.meansratio = zeros(length(distInds),length(grainSizes));
    modelOutPts2.(fields{ii}).dists.absmeans = zeros(length(distInds),length(grainSizes));
    modelOutPts2.(fields{ii}).dists.posmeans = zeros(length(distInds),length(grainSizes));
    modelOutPts2.(fields{ii}).dists.taudist = cell(length(distInds),1);
    for ll = 1 : length(distInds)
        modelOutPts2.(fields{ii}).dists.taudist{ll} = cat(1,edges{ll},zeros(size(edges{ll})));
    end
    for kk = 1 : length(grainSizes)
        modelOutPts2.(fields{ii}).dists.raw{kk} = [];
    end
    for jj = 1 : length(cumClusters2{ii})
        extractRaw = zeros(TS_LENGTH,length(distInds));
        peakFreq = zeros(1,length(distInds));
        uVar = zeros(1,length(distInds));
        for kk = 1 : length(distInds)
            uRaw = modelOutPts2.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
            [mod1Hsall,mod1Hrmsall,mod1Hs,mod1Hsig,mod1Hslf,mod1Tp,mod1Tpig,mod1Tplf,mod1Tdom,mod1Hrms,mod1Hrmsig,mod1Hrmslf] = spectral_wave_artificial(squeeze(uRaw)',0,modelSampleFreq,0,4096,[]);
            peakFreq(kk) = 2 * pi / mod1Tp;
            uVar(kk) = var(uRaw);
            extractRaw(:,kk) = modelOutPts2.(fields{ii}).(sprintf('c%d',jj)).(sprintf('d%d',distInds(kk))).raw;
        end
        acc = [];
        for kk = 1 : TS_LENGTH : size(extractRaw,1)
            extractTemp = cat(1,extractRaw(kk,:),extractRaw(kk:min(kk+TS_LENGTH-1,size(extractRaw,1)),:),extractRaw(min(kk+TS_LENGTH-1,size(extractRaw,1)),:));
            acc = cat(1,acc,(extractTemp(3:end,:)-extractTemp(1:end-2,:))./(2/modelSampleFreq*repmat(peakFreq,TS_LENGTH,1)));
        end
        
        for kk = 1 : length(grainSizes)
            Arms = sqrt(2)./peakFreq .* sqrt(uVar);
            f25 = exp(5.5.*(2.5.*grainSizes(kk)./Arms).^0.2-6.3);
            ustar = sqrt(0.5.*f25).*(cos(phiTaus(1)).*extractRaw + sin(phiTaus(1)).*acc);
            tau = rho.*sign(ustar).*ustar.^2;
            theta = tau ./ rho ./ ((rhoS/rho-1)*G*grainSizes(kk));
            baseThresh = (extractRaw<0).*repmat(negBaseThresh,size(extractRaw,1),1) + (extractRaw>=0).*repmat(posBaseThresh,size(extractRaw,1),1);
            if grainSizes(kk)*1000 == 2
                for ll = 1 : length(distInds)
                    counts = histcounts(tau(:,ll),edges{ll});
                    modelOutPts2.(fields{ii}).dists.taudist{ll}(2,:) = modelOutPts2.(fields{ii}).dists.taudist{ll}(2,:)+cat(2,counts,0)*cumClusters2{ii}(jj,end);
                end
            end
            modelOutPts2.(fields{ii}).dists.rawtau{kk} = cat(1,modelOutPts2.(fields{ii}).dists.rawtau{kk},sign(tau).*max(abs(tau)-criticalShear(kk),0).*clusters2{ii}(jj,end)/TS_LENGTH);
            modelOutPts2.(fields{ii}).dists.postau{kk} = cat(1,modelOutPts2.(fields{ii}).dists.postau{kk},(sign(tau).*max(abs(tau)-criticalShear(kk),0)>0).*clusters2{ii}(jj,end)/TS_LENGTH);
            modelOutPts2.(fields{ii}).dists.raw{kk} = cat(1,modelOutPts2.(fields{ii}).dists.raw{kk},A2.*sign(tau).*sqrt(abs(theta)).*max(abs(theta)-baseThresh,0).*sqrt(G.*(rhoS./rho-1).*grainSizes(kk)^3).*clusters2{ii}(jj,end)/TS_LENGTH);
        end
    end
    for jj = 1 : length(distInds)
        for kk = 1 : length(grainSizes)
            modelOutPts2.(fields{ii}).dists.means(jj,kk) = nansum(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj));
            if ~isempty(nansum(abs(modelOutPts2.(fields{ii}).dists.raw{kk}(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)<0,jj))))
                modelOutPts2.(fields{ii}).dists.meansratio(jj,kk) = nansum(modelOutPts2.(fields{ii}).dists.raw{kk}(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)>0,jj))/nansum(abs(modelOutPts2.(fields{ii}).dists.raw{kk}(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)<0,jj)));
            end
            modelOutPts2.(fields{ii}).dists.taumeans(jj,kk) = nansum(modelOutPts2.(fields{ii}).dists.rawtau{kk}(:,jj));
            modelOutPts2.(fields{ii}).dists.absmeans(jj,kk) = nansum(abs(modelOutPts2.(fields{ii}).dists.raw{kk}(:,jj)));
            modelOutPts2.(fields{ii}).dists.posmeans(jj,kk) = nansum(modelOutPts2.(fields{ii}).dists.postau{kk}(:,jj));
        end
    end
end


% Now plot the distributions for the sea level rise scenario
cats = {'raw','hi','lo'};

% Set the axis limits for all plots
ybounds = {[0,0.1],[0,0.1],[0,0.1],[0,0.1],[0,0.1], [0,0.11],[0,0.3],[0,0.5],[0,0.6],[0,0.1], [0,0.11],[0,0.3],[0,0.5],[0,0.6],[0,0.1]};
xbounds = {[-500,500],[-20,20],[-10,10],[-7.5,7.5],[-20,20]};
ybounds2 = {[-0.015,0.015],[-0.01,0.06],[-0.01,0.06],[-0.02,0.1],[-0.02,0.1], [-0.01,0.01]/25,[-5e-2,5e-2],[-1e-1,1e-1],[-2e-1,2e-1],[-1e-1,2e-1], [-0.05,0.05]/25,[-2e-2,2e-2],[-3e-2,3e-2],[-7e-2,7e-2],[-1e-1,3e-1]};
xbounds2 = {[-500,500],[-20,20],[-10,10],[-7.5,7.5],[-20,20]};
binWidths = [5,0.1,0.1,0.1,0.1]; % Widths of bins in histogram

%Initialize figure
fig = figure('PaperUnits','centimeters','PaperPosition',[0,0,maxWidth,minLength*1.5],'Name','Distributions');
fig.Renderer = 'painters';

% Plot model bathymetry and locations of distributions
subplot(5,5,2:4);
plot(bathyx,bathy.z,'k-');
hold on;
scatter(bathyx(distInds),bathy.z(distInds),'g^','MarkerFaceColor','g');
xlim(bathyx(bounds));
ylabel('Elevation (m aMSL)','FontSize',figureFontSize);
xlabel('Distance from reef crest (m)','FontSize',figureFontSize);
ax = gca; ax.XDir = 'reverse';

% Plot PDFs of shear stress for historical and +1m SLR scenarios
for ii = 1 : length(distInds)
    subplot(5,5,5+length(distInds)-ii+1);
    hold on;
    histogram('BinEdges',modelOutPts.(names{1}).dists.taudist{ii}(1,:),'BinCounts',modelOutPts.(names{1}).dists.taudist{ii}(2,1:end-1),'Normalization','pdf','DisplayStyle','stairs','EdgeColor',pltColors{1});
    histogram('BinEdges',modelOutPts2.(names2{2}).dists.taudist{ii}(1,:),'BinCounts',modelOutPts2.(names2{2}).dists.taudist{ii}(2,1:end-1),'Normalization','pdf','DisplayStyle','stairs','EdgeColor',pltColors{3});
    ylim(ybounds{ii});
    xlim([min(edges{ii}),max(edges{ll})]);
    
    if (ii == length(distInds))
        ylabel('Probability density function (PDF)','FontSize',figureFontSize);
    end
    ax = gca; ax.FontSize = figureFontSize; ax.XDir = 'reverse';
    
end

% Plot change in shear stress PDFs between historical and +1m SLR scenarios
ii = 1; 
    for jj = 1 : length(distInds)
        subplot(5,5,(ii+1)*5+length(distInds)-jj+1);
        hold on;
        lateBins = modelOutPts2.(names2{2}).dists.taudist{jj}(2,1:end-1);
        histBins = modelOutPts.(names{1}).dists.taudist{jj}(2,1:end-1);
        binWidths = nanmean(diff(edges{ii}));
        binChange = (lateBins ./ (sum(lateBins).*binWidths) - histBins ./ (sum(histBins).*binWidths)); 
        binChange(isnan(binChange)) = 0;
        
        plotEdges = (edges{ii}(1:end-1)+edges{ii}(2:end))/2;
        plot(plotEdges,binChange,'k-');
        plot([min(plotEdges),max(plotEdges)],[0,0],'r:');
        plot([0,0],[min(ybounds2{(ii-1)*5+jj}),max(ybounds2{(ii-1)*5+jj})],'r:');
        ylim(ybounds2{(ii-1)*5+jj});
%         ylim([-0.02,0.02]);
        xlim([min(plotEdges),max(plotEdges)]);
        if (jj == length(distInds)) & (ii == 1)
                ylabel('\Delta PDF','FontSize',figureFontSize);
        end
        if (ii == 1) & (jj == 3)
            xlabel('Bottom shear stress, \tau_{b} (Nm^{-2})','FontSize',figureFontSize);
        end
        ax = gca; ax.FontSize = figureFontSize; ax.XDir = 'reverse';
    end


% print(fullfile(figureOutDir,'ShearStressDistributions_NielsenCallaghan2003'),'-depsc','-r600');

%..........................................................................
