paramFhil.m

clc; clear; close all;
set(groot,'defaultAxesXGrid','on');
set(groot,'defaultAxesYGrid','on');
set(groot, 'defaultFigureUnits', 'centimeters', 'defaultFigurePosition', [5 2 40 15]);
set(0,'defaultAxesFontSize',18);

% Add subfolders to path
addpath(genpath(pwd));

% Load datafiles
load('Matfiles/E_hil','Ehil6dm_IC3','Ehil6dm_M1');
load('Data/Mixing_parameterization_fields.mat','decay_scale_bot');
load('Matfiles/N_is.mat','N2_IC3','N2_M1','N_IC3','N_M1','zmid_IC3','zmid_M1','time'); % Contains zmid_M1 and zmid_IC3

% Mooring tags (for image generation)
mooring = ["IC3","M1"];

%% Find the depth-integrated N^2

dz_IC3 = cat(1,10*ones(5,1),25*ones(6,1),50*ones(30,1));
dz_M1 = cat(1,10*ones(5,1),25*ones(6,1),50*ones(31,1));
N2_int_IC3 = nan(length(time),1);
N2_int_M1 = nan(length(time),1);

for i=1:length(time)
    N2_int_IC3(i) = sum(N2_IC3(:,i).*dz_IC3,'omitnan');
    N2_int_M1(i) = sum(N2_M1(:,i).*dz_M1,'omitnan');
end

N2_int_IC3(N2_int_IC3 == 0) = nan;
N2_int_M1(N2_int_M1 == 0) = nan;

N2_int3D_IC3 = repmat(N2_int_IC3,1,1,length(zmid_IC3(:,1)));
N2_int3D_IC3 = squeeze(N2_int3D_IC3(:,1,:))';

N2_int3D_M1 = repmat(N2_int_M1,1,1,length(zmid_M1(:,1)));
N2_int3D_M1 = squeeze(N2_int3D_M1(:,1,:))';

%% IC3: N^2(t) [vertically-integrated N^2(t)]

ax1 = figure;
plot(time,N2_int_IC3);
xlabel('Time');
ylabel('$\int N^2 dz$','Interpreter','latex');
title('IC3: $\int N^2 dz (t)$','Interpreter','latex');

exportgraphics(ax1,'figures/main/HIL/' + mooring(1) + '_N2t.png');

%% M1: N^2(t) [vertically-integrated N^2(t)]

ax2 = figure;
plot(time,N2_int_M1);
xlabel('Time');
ylabel('$\int N^2 dz$','Interpreter','latex');
title('M1: $\int N^2 dz (t)$','Interpreter','latex');

exportgraphics(ax2,'figures/main/HIL/' + mooring(2) + '_N2t.png');

%% IC3: N(z,t) log plot

% Put time in datenum meshgrid format for filled contour plot
timeNum_IC3 = meshgrid(datenum(time'),zmid_IC3(:,1));
timeNum_M1 = meshgrid(datenum(time'),zmid_M1(:,1));

ax3 = figure;
contourf(timeNum_IC3,zmid_IC3,log10(N2_IC3),-8.5:0.1:-4.5,'LineColor','none');
ylabel('Depth [m]');
datetick('x','yyyy mmm','keeplimits');
colormap(flipud(cbrewer2('RdPu')));
ylim([-1525 -100]);
% cmocean('tempo');
xlabel('Time'); ylabel('Depth [m]'); title('IC3: $N^2(z,t)$','Interpreter','latex')
c = colorbar;
c.Label.String = 'log_{10}(N^2(z,t)) [s^{-1}]';
exportgraphics(ax3,'figures/main/HIL/' + mooring(1) + '_Nsquared.png');

%% M1: N(z,t) log plot
ax4 = figure;
contourf(timeNum_M1,zmid_M1,log10(N2_M1),-8.5:0.1:-4.5,'LineColor','none');
datetick('x','yyyy mmm','keeplimits');
% cmocean('tempo');
colormap(flipud(cbrewer2('RdPu')));
c = colorbar;
ylim([-1625 -100]);
xlabel('Time'); ylabel('Depth [m]'); title('M1: $N^2(z,t)$','Interpreter','latex')
c.Label.String = 'log_{10}(N^2(z,t)) [s^{-1}]';
exportgraphics(ax4,'figures/main/HIL/' + mooring(2) + '_Nsquared.png');

%% Modulate Hbot

% Hbot for IC3 and M1
Hbot_IC3 = decay_scale_bot(657,279);
Hbot_M1 = decay_scale_bot(660,279);

% Normalisation Factor 'nf'
nf = zeros(2,length(time));
for t = 1:length(time)
    nf(1,t) = 1/sqrt(Ehil6dm_IC3(t));
    nf(2,t) = 1/sqrt(Ehil6dm_M1(t));
end

% Tuning Factor 'tf'
tf = zeros(2,length(time));

normalisation = nan(2,52873);
for i = 1:2
    tf(i,:) = 1./mean((nf(i,:) - min(nf(i,:))) ./ (max(nf(i,:)) - min(nf(i,:))),'omitnan');

    % Normalise such that the mean of nf = E_hil(x,y)
    normalisation(i,:) = tf(i,:) .* (nf(i,:) - min(nf(i,:))) ./ (max(nf(i,:)) - min(nf(i,:)));
end

HbotIC3_norm = normalisation(1,:).*Hbot_IC3;
HbotM1_norm = normalisation(2,:).*Hbot_M1;
meanHBotIC3preNorm = mean(HbotIC3_norm,'omitnan');
meanHBotM1preNorm = mean(HbotM1_norm,'omitnan');

% Apply a floor to the modulated values for Hbot
for i=1:length(time)
    if HbotIC3_norm(i) < 100
        HbotIC3_norm(i) = 100;
    end
    if HbotM1_norm(i) < 100
        HbotM1_norm(i) = 100;
    end
end

ax5 = figure;
plot(time,HbotIC3_norm);
hold on
yline(meanHBotIC3preNorm,'LineWidth',2,'Color','blue');
yline(Hbot_IC3);
hold off
xlabel('Time'); ylabel('$H_{bot}$','Interpreter','latex');
title('IC3: $H_{bot}(t)$','Interpreter','latex');
exportgraphics(ax5,'figures/main/HIL/' + mooring(1) + '_Hbot.png');

ax6 = figure;
plot(time,HbotM1_norm);
hold on
yline(meanHBotM1preNorm,'LineWidth',2,'Color','blue');
yline(Hbot_M1);
hold off
xlabel('Time'); ylabel('$H_{bot}$','Interpreter','latex');
title('M1: $H_{bot}(t)$','Interpreter','latex');
exportgraphics(ax6,'figures/main/HIL/' + mooring(2) + '_Hbot.png');

%% IC3: Compute the vertical structure

% depthsIC3 = zq;
hab_IC3 = nan(1,41);

% From Deployment Overview [De Jong & Fried]
IC3bathy = 1635;
for i = 1:length(zmid_IC3(:,1))
    hab_IC3(i) = IC3bathy + zmid_IC3(i,1);
end
hab_IC3(end) = 0;

if isfile('Matfiles/Fhil_IC3.mat')
    disp('IC3: Time-dependent vertical structure already calculated...');
    load('Matfiles/Fhil_IC3.mat');
else
    disp('Calculating vertical structure...');
    r_bot = 0.86;
    T1_IC3 = nan(41,52873); T2_IC3 = nan(41,52873);
    for t = 1:length(time)
        for k = 1:length(zmid_IC3(:,1))
            T1_IC3(k,t) =  r_bot .* (1./(1 + (hab_IC3(k)./HbotIC3_norm(t))).^2) ...
                            .*(1./IC3bathy + 1./HbotIC3_norm(t));
            T2_IC3(k,t) = (1 - r_bot).*(N2_IC3(k,t)./N2_int3D_IC3(k,t));
        end
    disp(t)
    end
    Fhil_IC3 = T1_IC3 + T2_IC3;
    save Matfiles/Fhil_IC3.mat Fhil_IC3 T1_IC3 T2_IC3;
end


%% M1: Vertical Structure

% M1 coords = [660,279]
hab_M1 = nan(1,42);

% M1 Bathymetry - from Deployment Overview [De Jong & Fried]
M1bathy = 1712;
for i = 1:length(zmid_M1(:,1))
    hab_M1(i) = M1bathy + zmid_M1(i,1);
end
hab_M1(end) = 0;

if isfile('Matfiles/Fhil_M1.mat')
    disp('M1: Time-dependent vertical structure already calculated...');
    load('Matfiles/Fhil_M1.mat');
else
    disp('Calculating vertical structure for M1...');
    r_bot = 0.86;
    T1_M1 = nan(42,52873); T2_M1 = nan(42,52873);
    for t = 1:length(time)
        for k=1:length(zmid_M1(:,1))
            T1_M1(k,t) =  r_bot .* (1./(1 + (hab_M1(k)./HbotM1_norm(t))).^2) ...
                            .*(1./M1bathy + 1./HbotM1_norm(t));
            T2_M1(k,t) = (1 - r_bot).*(N2_M1(k,t)./N2_int3D_M1(k,t));
        end
    disp(t)
    end
    Fhil_M1 = T1_M1 + T2_M1;
    save Matfiles/Fhil_M1.mat Fhil_M1 T1_M1 T2_M1;
end

%% IC3: F(z,t), not binned

ax7 = figure;
contourf(timeNum_IC3,zmid_IC3,log10(Fhil_IC3),-2.5:-0.1:-4.5,'LineColor','none');
datetick('x','yyyy mmm','keeplimits');
colormap(flipud(cbrewer2('RdPu')));
c = colorbar;
c.Label.String = 'log_{10}(F(z,t)) [m^{-1}]';
ylim([-1525 -100]);
xlabel('Time'); ylabel('Depth [m]');
title('IC3: $F_{hil}(t)$','Interpreter','latex');
exportgraphics(ax7,'figures/main/HIL/' + mooring(1) + '_Fhil.png');

%% M1: F(z,t), not binned

ax8 = figure;
contourf(timeNum_M1,zmid_M1,log10(Fhil_M1),-2.5:-0.1:-4.5,'LineColor','none');
datetick('x','yyyy mmm','keeplimits');
colormap(flipud(cbrewer2('RdPu')));
c = colorbar;
c.Label.String = 'log_{10}(F(z,t)) [m^{-1}]';
ylabel('Depth [m]');
ylim([-1625 -100]);
title('M1: $F_{hil}(t)$','Interpreter','latex');
exportgraphics(ax8,'figures/main/HIL/' + mooring(2) + '_Fhil.png');

%% Kurtosis of F(z,t)
% % MAYBE I will come back to this...
% kurtF = kurtosis(Fhil_IC3,0);
% 
% figure
% plot(time,kurtF);
% xlabel('time');
% ylabel('K_F (Kurtosis-value of F(z,t)');
% 
% 
% lessProneToOutliers = 0;
% moreProneToOutliers = 0;
% 
% for i=1:length(time)
%     if kurtF(i) < 3
%         lessProneToOutliers = lessProneToOutliers + 1;
%     end
%     if kurtF(i) > 3
%         moreProneToOutliers = moreProneToOutliers + 1;
%     end
% end
% 
% propProne = 100*moreProneToOutliers./length(time);

%% ADJUST LIMITS
set(groot, 'defaultFigureUnits', 'centimeters', 'defaultFigurePosition', [5 2 25 22]);
set(0,'defaultAxesFontSize',18);

%% IC3: Plot of Vertical Structure components T1 and T2

T1_IC3 = load('Matfiles/Fhil_IC3.mat').T1_IC3;
T2_IC3 = load('Matfiles/Fhil_IC3.mat').T2_IC3;

ax9 = figure;

subplot(2,1,1)
contourf(timeNum_IC3,zmid_IC3,log10(T1_IC3),-2.5:-0.05:-4.5,'LineColor','none');
datetick('x','yyyy mmm','keeplimits');
c = colorbar;
c.Label.String = 'log_{10} (F_{hil}(z,t)) [m^{-1}]';
colormap(flipud(cbrewer2('RdPu')));
ylim([-1525 -100]);
ylabel('Depth [m]');
title('Topography Component');

subplot(2,1,2)
contourf(timeNum_IC3,zmid_IC3,log10(T2_IC3),-2.5:-0.05:-4.5,'LineColor','none');
datetick('x','yyyy mmm','keeplimits');
c = colorbar;
c.Label.String = 'log_{10} (F_{hil}(z,t)) [m^{-1}]';
colormap(flipud(cbrewer2('RdPu')));
ylim([-1525 -100]);
xlabel('Time'); ylabel('Depth [m]');
title('Stratification Component');

sgtitle('IC3: components contributing to $F_{hil}(z,t)$','Interpreter','latex');
exportgraphics(ax9,'figures/main/HIL/' + mooring(1) + '_Fhil_bothTerms.png');

%% M1: Plot of Vertical Structure components T1 and T2

T1_M1 = load('Matfiles/Fhil_M1.mat').T1_M1;
T2_M1 = load('Matfiles/Fhil_M1.mat').T2_M1;

ax10 = figure;

subplot(2,1,1)
contourf(timeNum_M1,zmid_M1,log10(T1_M1),-2.5:-0.1:-4.5,'LineColor','none');
datetick('x','yyyy mmm','keeplimits');
c = colorbar;
c.Label.String = 'log_{10} (F(z,t)) [m^{-1}]';
colormap(flipud(cbrewer2('RdPu')));
ylim([-1625 -100]);
ylabel('Depth [m]');
title('Topography Component');

subplot(2,1,2)
contourf(timeNum_M1,zmid_M1,log10(T2_M1),-2.5:-0.1:-4.5,'LineColor','none');
datetick('x','yyyy mmm','keeplimits');
c = colorbar;
c.Label.String = 'log_{10} (F(z,t)) [m^{-1}]';
colormap(flipud(cbrewer2('RdPu')));
ylim([-1625 -100]);
xlabel('Time'); ylabel('Depth [m]');
title('Stratification Component');

sgtitle('M1: components contributing to $F_{hil}(z,t)$','Interpreter','latex');
exportgraphics(ax10,'figures/main/HIL/' + mooring(2) + '_Fhil_bothTerms.png');

%% Save parameters for next file

save Matfiles/Fhil.mat Fhil_IC3 Fhil_M1;

clear ax1 ax2 ax3 ax4 ax5 ax6 ax7 ax8 ax9 ax10 c decay_scale_bot i t;