function [bp, sla_model_adj] = Control_BreakPoints(t,MD,sla_model_adj,bp,sPcum,TGEQ)
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here

% Visual Inspection: 0 = Earthquake; 1 = datum shift (flagged); 2 = datum 
% shift (false detection); 3 = unflagged shifts; 4 = unflagged shift, not
% corrected
bp(:,3) = NaN;
bp(sPcum,3) = 1;
bp(:,4) = TGEQ;
bp(:,5) = NaN;
bp(4,5) = 1;
bp(5,5) = 3;

% HONNGU manual shift Adjustment
i = 5;
isn = find(~isnan(MD(:,i)));
ts1 = MD(isn,i);
t1 = t(isn);
s1 = find(t1<=1995);
s2 = find(t1>1995&t1<2001);
s3 = find(t1>=2001);
[Ts1,yf1] = trend_new(1,t1(s1),ts1(s1),'N');
[Ts2,yf2] = trend_new(1,t1(s3),ts1(s3),'N');
sla_model_adj(isn(s1),i) = yf1(end);  
sla_model_adj(isn(s3),i) = yf2(1); 
sla_model_adj(isn(s2),i) = NaN; 

bp(7,5) = 1; % Earthquakes occcured, but seems to be after the bp detection
bp(8,5) = 2;
bp(9,5) = 0;
bp(16,5) = 0; % River discharge probably a factor; needs to be further evaluated
bp(24,5) = 0;
bp(35,5) = 0;
bp(57,5) = 0;
bp(61,5) = 4;
bp(66,5) = 0;
bp(67,5) = 4;
bp(72,5) = 0;
bp(97,5) = 3;
bp(108,5) = 3;
bp(110,5) = 2;
bp(116,5) = 0;
bp(117,5) = 4;
bp(120,5) = 0;
bp(121,5) = 0;
bp(123,5) = 0;
bp(132,5) = 0;
bp(133,5) = 0;
bp(134,5) = 0;
bp(135,5) = 0;
bp(141,5) = 0;
bp(143,5) = 0;
bp(155,5) = 4;

% Skuru, changepoint analysis doesn't detect anything, but autocorrelation
% parameters and visual inspection indicate 3 changepoints. Taking the
% differences to the neighboring station of Jungfrusund lead to similar
% changepoints, suggesting that the changepoints are indeed datum shifts

bp(159,5) = 3;
i = 159;
ts = MD(:,i);
s1 = find(t<=1911);
s2 = find(t>1911&t<=1930);
s3 = find(t>1930);
mdl = fitlm(t(s1),MD(s1,i),'linear');
yf1 = predict(mdl,t(s1));
mdl = fitlm(t(s2),MD(s2,i),'linear');
yf2 = predict(mdl,t(s2));
mdl = fitlm(t(s3),MD(s3,i),'linear');
yf3 = predict(mdl,t(s3));
sla_model_adj(s1,i) = zeros;  
sla_model_adj(s2,i) = yf2(1)-yf1(end); 
sla_model_adj(s3,i) = yf3(1)-yf2(end)+yf2(1)-yf1(end); 

bp(160,5) = 4;
bp(166,5) = 0;
bp(176,5) = 3;
bp(181,5) = 3;
bp(189,5) = 4;
bp(199,5) = 4;

% Tofino, changepoint analysis doesn't detect anything, because the
% changepoint appears at the time series beginning. It shows up in enhanced
% uncertainty estimates though. PSMSL notes that data before 1920 is
% currently queried.

bp(225,5) = 3;
i = 225;
isn = find(~isnan(MD(:,i)));
ts1 = MD(isn,i);
t1 = t(isn);
s1 = find(t1<=1917);
s2 = find(t1>1917);
[Ts1,yf1] = trend_new(1,t1(s1),ts1(s1),'N');
[Ts2,yf2] = trend_new(1,t1(s2),ts1(s2),'N');
sla_model_adj(isn(s1),i) = yf1(end);  
sla_model_adj(isn(s2),i) = yf2(1); 

bp(239,5) = 0;

% Vlissingen: All Dutch values are 4 values per day, i.e. highs and lows,
% those can be heavily affected by tidal deformations
bp(261,5) = 3;
sla_model_adj(1:20,261) = NaN;

bp(264,5) = 3;
bp(266,5) = 4;
bp(271,5) = 1;
bp(279,5) = 3;

% Kemi, changepoint analysis doesn't detect anything, but inflated
% autocorrelation parameters and visual inspection point toward
% changepoint in 1976/1977, when PSMSL notes that the station was
% relocated. Comparison to nearby location of Raahe/Brahestadt also
% confirms this changepoint.
bp(296,5) = 3;
i = 296;
isn = find(~isnan(MD(:,i)));
ts1 = MD(isn,i);
t1 = t(isn);
s1 = find(t1<=1977);
s2 = find(t1>1977);
[Ts1,yf1] = trend_new(1,t1(s1),ts1(s1),'N');
[Ts2,yf2] = trend_new(1,t1(s2),ts1(s2),'N');
sla_model_adj(isn(s1),i) = yf1(end);  
sla_model_adj(isn(s2),i) = yf2(1); 


bp(318,5) = 4;

i = 319;
bp(319,5) = 3;
load('Astoria_DatumCorrections_Talke2020.mat')
sla_model_adj(:,i) = T20_AstoriaCorr;

bp(322,5) = 0;
bp(331,5) = 4; % Sheerness, undocumented shift before 1940s that needs to be communicated with PSMSL, differences also contain nodal modulations
bp(338,5) = 4; % Narvik, shift before gap, data too short to properly apply shift
bp(347,5) = 3; 

i = 347;
sla_model_adj(91:93,i) = NaN;
isn = find(~isnan(MD(:,i)));
ts1 = MD(isn,i);
t1 = t(isn);
s1 = find(t1<=1999);
s2 = find(t1>1999);
[Ts1,yf1] = trend_new(1,t1(s1),ts1(s1),'N');
[Ts2,yf2] = trend_new(1,t1(s2),ts1(s2),'N');
sla_model_adj(isn(s1),i) = yf1(end);  
sla_model_adj(isn(s2),i) = yf2(1); 

bp(350,5) = 4; 
bp(353,5) = 4;
bp(357,5) = 4;
bp(358,5) = 0;
bp(366,5) = 2; % This is located in the Houghly River and likely affected by river changes
bp(367,5) = 2;
bp(369,5) = 3;
bp(370,5) = 4;
bp(392,5) = 3;
bp(409,5) = 4; % Obviously groundwater affected. Therefore the detected bps are unclear and need more attention, they could be earthquake related as well
bp(410,5) = 0;
bp(411,5) = 0;
bp(418,5) = 4;
bp(421,5) = 4; % complex coastal environment

% La Coruna 1 has a changepoint in 1963, which is not detected by any
% changepoint detection algorithm but which shows up in inflated
% autocorrelation parameters. 
% Marcos et al (2005, https://digital.csic.es./bitstream/10261/93204/1/Marcos-JGR-2005-v110-pC03008.pdf)
% indeed also detect a changepoint in 1963. We apply their estimate of -4cm
% as a correction here.

bp(427,5) = 3;
i = 427;
sla_model_adj(s1,i) = zeros;  
sla_model_adj(s2,i) = -0.04; 

bp(434,5) = 3; % Skagway; There are also deformations visible after the 1964 earthquake in Alaska

i = 441; % Puerto Madryn, Argentnia: data has been flagged, but too short to apply CP, so deleted
sla_model_adj(94:97,i) = NaN;
bp(441,5) = 3; 

bp(460,5) = 4; % Davao, Davao Gulf: As noted by PSMSL, this is likley full of datum shifts, communicate with PSMSL; also earthquake prone
bp(464,5) = 3; % Barentsburg, undocumented shift. PSMSL notes suspicion, but no flag
bp(468,5) = 4;
bp(469,5) = 3; % Barentsburg II
bp(472,5) = 3; % Bodo, change in instrumentation, communicate to PSMSL
bp(481,5) = 4;
bp(488,5) = 4;
bp(489,5) = 4;
bp(490,5) = 4;
bp(492,5) = 4;
bp(493,5) = 3; % Dikson, Arctic; needs to be communicated to PSMSL
bp(503,5) = 0; % Onahama
bp(506,5) = 3; % Reykjavik, undocumented shift that needs to be communicated with PSMSL 
bp(511,5) = 3; % Zhelania II, undocumented shift that needs to be communicated with PSMSL
bp(512,5) = 4; % Possible shift in the 1920s, communication with P. Hogarth?
bp(516,5) = 1;
bp(518,5) = 4;
bp(521,5) = 4; % Murmansk, in river location, should be assessed in greater detail
bp(523,5) = 4; % Murmansk, in river location, should be assessed in greater detail

% Hanko, changepoint analysis doesn't detect anything, but autocorrelation
% parameters and visual inspection indicate 2 changepoints. Taking the
% differences to the neighboring station of Turku lead to similar
% changepoints, suggesting that the changepoints are indeed datum shifts

bp(532,5) = 3;
i = 532;
ts = MD(:,i);
s1 = find(t<=1936);
s2 = find(t>1936&t<=1985);
s3 = find(t>1985);
mdl = fitlm(t(s1),MD(s1,i),'linear');
yf1 = predict(mdl,t(s1));
mdl = fitlm(t(s2),MD(s2,i),'linear');
yf2 = predict(mdl,t(s2));
mdl = fitlm(t(s3),MD(s3,i),'linear');
yf3 = predict(mdl,t(s3));
sla_model_adj(s1,i) = zeros;  
sla_model_adj(s2,i) = yf2(1)-yf1(end); 
sla_model_adj(s3,i) = yf3(1)-yf2(end)+yf2(1)-yf1(end); 

bp(540,5) = 2; % Freeport; obviously affected by fluid withdrawals, no shift
bp(547,5) = 0;

i = 549; % Lowestoft, following Hogarth, we delete anomalous values before 1980
sla_model_adj(78:81,i) = NaN;
bp(549,5) = 3;

bp(550,5) = 4;

i = 558;
isn = find(~isnan(MD(:,i)));
ts1 = MD(isn,i);
t1 = t(isn);
s1 = find(t1<=1980);
s2 = find(t1>1980);
[Ts1,yf1] = trend_new(1,t1(s1),ts1(s1),'N');
[Ts2,yf2] = trend_new(1,t1(s2),ts1(s2),'N');
sla_model_adj(isn(s1),i) = yf1(end);  
sla_model_adj(isn(s2),i) = yf2(1); 
bp(558,5) = 3; % Undocumented changepoint; PSMSL noted contriversies, but needs to be communicated

bp(560,5) = 4;
bp(567,5) = 4;
bp(587,5) = 0; % Earthquake site with destructive earthquake and potential post-semic VLM after 1994
bp(601,5) = 0;
bp(611,5) = 4; % Malin Head; possible shifts that need more attention; 1990s values are flagged but not 2000s 

i = 615; % Heysham; data flags and deletion
sla_model_adj(104:105,i) = NaN;
bp(615,5) = 3;

bp(616,5) = 0;
bp(617,5) = 4; % North Shields, needs further investigation

sla_model_adj(:,bp(:,5)==0) = zeros;

end