-
Notifications
You must be signed in to change notification settings - Fork 0
/
paramEniw.m
175 lines (137 loc) · 5.86 KB
/
paramEniw.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
clc; clear; close all;
set(groot,'defaultAxesXGrid','on');
set(groot,'defaultAxesYGrid','on');
set(groot, 'defaultFigureUnits', 'centimeters', 'defaultFigurePosition', [5 5 40 15]);
set(0,'defaultAxesFontSize',18);
% Add subfolders to path
addpath(genpath(pwd));
%% update each time with MOORING used - used for saving plots
mooring = "IC3";
%% Load data files
load('merged_hourly_unfiltered_data_20142020.mat', 'lat'); % Load latitude of mooring
load('IC3.mat','t'); % Load time
load('niwEnergy.mat');
%% Normalise the modulation
% Here we modulate the E_wwi signal according to the interpolated 6-month
% energy contained within the combined near-inertial and incoherent
% semidiurnal tidal bands.
% The normalisation factor = niwInterp.
% This must be normalised such that the mean of the time-modulated E_wwi is
% equal to the single-valued E_wwi parametrisation.
tuningFactor = 1./mean((niwInterp - min(niwInterp)) ./ (max(niwInterp) - min(niwInterp)));
normalisation = tuningFactor*(niwInterp - min(niwInterp)) ./ (max(niwInterp) - min(niwInterp));
tuningFactorAtSinglePts = 1./mean((niwEnergy - min(niwEnergy)) ./ (max(niwEnergy) - min(niwEnergy)));
normAtSinglePts = tuningFactorAtSinglePts*(niwEnergy - min(niwEnergy)) ./ (max(niwEnergy) - min(niwEnergy));
clear tuningFactorAtSinglePts;
%% E_niw: Depth-integrated IT dissipation rate due to NIW
% We assume that near-inertial waves (NIW) decay through wave-wave
% interactions (WWI in the parametrisation of De Lavergne et al, 2020).
% Therefore we will use the vertical structure for WWI from that
% parametrisation combined with the estimate for parametrised energy input
% to mixing by near-inertial waves from (Alford, 2020)
alfordIC3 = 3.3e-4; % W/m2
% Eniw_IC3 = normalisation*power_wwi(657,279);
Eniw_IC3 = normalisation*alfordIC3;
% Eniw_IC3_singlePts = normAtSinglePts*power_wwi(657,279);
Eniw_IC3_singlePts = normAtSinglePts*alfordIC3;
Eniw6DM_IC3 = movmean(Eniw_IC3,149);
stdEniw_IC3 = std(Eniw6DM_IC3);
kurEniw_IC3 = kurtosis(Eniw6DM_IC3);
% timePoints = [t(floor(0.5*52873/12)),t(floor(1.5*52873/12)),...
% t(floor(2.5*52873/12)),t(floor(3.5*52873/12)),t(floor(4.5*52873/12)),...
% t(floor(5.5*52873/12)),t(floor(6.5*52873/12)),t(floor(7.5*52873/12)),...
% t(floor(8.5*52873/12)),t(floor(9.5*52873/12)),t(floor(10.5*52873/12)),...
% t(floor(11.5*52873/12))];
% E_wwi(t) incl. mean and std
ax1 = figure;
plot(t,Eniw6DM_IC3,'DisplayName','$E_{niw}(t)$ (6DM)');
hold on
% scatter(timePoints,Eniw_IC3_singlePts,'DisplayName','E_{wwi}(t) (non-interpolated)');
% yline(power_wwi(657,279),'-','DisplayName','E_{wwi}');
yline(alfordIC3,'-','DisplayName','$E_{niw}$');
% yline(mean(Eniw6DM_IC3),'--','DisplayName','Mean E_{hil(t)}','LineWidth',3);
% yline(mean(Eniw6DM_IC3)+stdEniw_IC3,'--','DisplayName','\mu + \sigma','Alpha',0.3);
% yline(mean(Eniw6DM_IC3)-stdEniw_IC3,'--','DisplayName','\mu - \sigma','Alpha',0.3);
hold off
legend('Interpreter','latex');
xlabel('Time');
ylabel('$E_{niw} [W m^{-2}]$','Interpreter','latex');
title('IC3: Depth-integrated dissipation rate due to NIW ($E_{niw}(t)$)','Interpreter','latex');
exportgraphics(ax1,'figures/main/NIW/' + mooring + '_Eniw6DM.png');
%% E_wwi(t): zoomed-in
t1 = 14300;
t2 = 15700;
ax2 = figure;
subplot(1,2,1)
plot(t,Eniw6DM_IC3,'DisplayName','$E_{niw}(t)$ (6DM)','LineWidth',1.5);
legend('Interpreter','latex');
xlabel('Time');
ylabel('$E_{niw} [W m^{-2}]$','Interpreter','latex');
title('$E_{niw}(t)$ (6DM)','Interpreter','latex');
subplot(1,2,2)
plot(t(t1:t2),Eniw6DM_IC3(t1:t2),'DisplayName','$E_{niw}(t)$ (6DM)','LineWidth',1.5);
legend('Interpreter','latex');
xlabel('Time');
ylabel('$E_{niw} [W m^{-2}]$','Interpreter','latex');
title('$E_{niw}(t)$ (6DM): close-up','Interpreter','latex');
exportgraphics(ax2,'figures/main/NIW/' + mooring + '_Eniw6DM_zoom.png');
%% V2. Histogram. Distribution of E.
ax3 = figure;
histogram(Eniw6DM_IC3);
hold on
xline(mean(Eniw6DM_IC3),':','\mu','LineWidth',1.5);
xline(mean(Eniw6DM_IC3)+stdEniw_IC3,':','\mu + \sigma','LineWidth',1.5);
xline(mean(Eniw6DM_IC3)-stdEniw_IC3,':','\mu - \sigma','LineWidth',1.5);
hold off
xlabel('$E_{niw} [W m^{-2}]$','Interpreter','latex');
ylabel('No. of values');
title('Distribution of $E_{niw}(t)$','Interpreter','latex');
exportgraphics(ax3,'figures/main/NIW/' + mooring + '_Eniw_hist.png');
%% Fourier Analysis of E_niw.
Ts = 3600;
fs = 2*pi/Ts;
L = length(t);
fftEniw = fft(Eniw6DM_IC3);
P2 = abs(fftEniw/L);
P1 = P2(1:floor(L/2));
f = (0:L/2-1)*fs/L;
[freqs,names] = frequencies_PF;
Omega = 7.2921e-5;
% f_M2 = freqs(45)/(24*3600);
f_M4 = freqs(76)/(24*3600);
f_S1 = freqs(20)/(24*3600);
% f_S3 = freqs(67)/(24*3600);
fC = 2*Omega*sin(deg2rad(lat(4)));
ax4 = figure;
loglog(f,P1,'DisplayName','$E_{niw}(t)$');
hold on
% xline(f_M2,':',names(45),'DisplayName','M2');
xline(f_M4,':',names(76),'DisplayName','M4');
xline(fC,':','f','DisplayName','f');
xline(f_S1,':',names(20),'DisplayName','S1');
% xline(f_S3,':',names(67),'DisplayName','S3');
hold off
legend('Interpreter','latex');
xlabel('$f [s^{-1}]$','Interpreter','latex');
ylabel('$|E_{niw}(t)|^{2} [W^2 m^{-4}]$','Interpreter','latex');
title('DFT: U component of E');
exportgraphics(ax4,'figures/main/NIW/' + mooring + '_Eniw_FFT.png');
%%
clear f_M2 f_M4 f_S1 f_S3 fC fftEniw freqs names fs Ts Omega P2 L lat;
clear f P1; % comment off/on to see FFT output
%% Evaluate how much of an effect smoothing has on the mean
windowSizes = 6:143:42906;
sad = nan(1,length(windowSizes));
for k = 1 : length(windowSizes)
Eniw_IC3_sm = movmean(Eniw_IC3, windowSizes(k));
sad(k) = sum(abs(Eniw_IC3_sm - Eniw_IC3));
end
ax5 = figure;
plot(windowSizes, sad, 'b*-', 'LineWidth', 2);
xlabel('Window Size');
ylabel('SAD');
title('Summed Absolute Difference (??)');
exportgraphics(ax5,'figures/main/NIW/' + mooring + '_Eniw_SAD.png');
%% Save parameters for the next file
save Matfiles/E_niw.mat Eniw6DM_IC3 Eniw_IC3;
clear ax1 ax2 ax3 ax4 ax5 k;