%% (1) LOAD HYBRID RECONSTRUCTION

Path = 'C:\Users\DangendorfSoenke\Box\Research Works and Working Papers\2023_GMSL Budget KS\Submission Materials\ZENODO\';

load([Path 'KSSLfin.mat'])
load([Path 'KSSL_SEfin.mat'])
RGIAHRSE = RGIAHRwmSE'.*tt; RGIAHRSE = RGIAHRSE-RGIAHRSE(end,:);
GGIAHRSE = GGIAHRwmSE'.*tt; GGIAHRSE = GGIAHRSE-GGIAHRSE(end,:);
s = find(tt>=1900);
t = tt(s); HR = SL_multi(:,s); SDSL = SDSL(:,s); BarySL = BarySL(:,s); 
GIA_Field_KS = GIA_Field_KS(:,s); IBAint = IBAint(:,s); HRrel = HR+GIA_Field_KS;
% HR = SDSL+BarySL+IBAint; HR = SDSL+BarySL+IBAint+GIA_Field_KS
% SDSL = Sterodynamic Sea Level; BarySL = Barystatic Sea Level
% IBAint = Inverse Barometer Effect; GIA_Field_KS = Glacial Isostatic
% Adjustment
SE = SE_SL_multi(:,s); GMSLHRSE = GMSLHRSE(s); GMSSLHRSE = GMSSLHRSE(s);
GBSLHRSE = GBSLHRSE(s); GGIAHRSE = GGIAHRSE(s); RMSLHRSE = RMSLHRSE(s,:); 
RMSSLHRSE = RMSSLHRSE(s,:); RBSLHRSE = RBSLHRSE(s,:); RGIAHRSE = RGIAHRSE(s,:);
SE_SDSL = SE_SDSL(:,s);SE_Bary = SE_Bary(:,s);
wSE_GIA = SE_GIA;
SE_GIA = SE_GIA*t'; SE_GIA = SE_GIA-SE_GIA(:,end);
SE_HR = SE; SE_HRrel = sqrt(SE_HR.^2+SE_GIA.^2); 
Lcg = LALT;

for i = 1:size(t,1);
    GMSLHR(i,1) = sum(HR(:,i).*cosd(Lcg(:,2)))./sum(cosd(Lcg(:,2)));
    GMSLHRrel(i,1) = sum(HRrel(:,i).*cosd(Lcg(:,2)))./sum(cosd(Lcg(:,2)));
    GGIAHR(i,1) = sum(GIA_Field_KS(:,i).*cosd(Lcg(:,2)))./sum(cosd(Lcg(:,2)));
    GMSSLHR(i,1) = sum(SDSL(:,i).*cosd(Lcg(:,2)))./sum(cosd(Lcg(:,2)));
    GBSLHR(i,1) = sum(BarySL(:,i).*cosd(Lcg(:,2)))./sum(cosd(Lcg(:,2)));
    for j = 1:6;
        RMSLHR(i,j) = sum(HR(regionSL==j,i).*cosd(Lcg(regionSL==j,2)))./sum(cosd(Lcg(regionSL==j,2)));
        RMSLHRrel(i,j) = sum(HRrel(regionSL==j,i).*cosd(Lcg(regionSL==j,2)))./sum(cosd(Lcg(regionSL==j,2)));
        RGIAHR(i,j) = sum(GIA_Field_KS(regionSL==j,i).*cosd(Lcg(regionSL==j,2)))./sum(cosd(Lcg(regionSL==j,2)));
        RMSSLHR(i,j) = sum(SDSL(regionSL==j,i).*cosd(Lcg(regionSL==j,2)))./sum(cosd(Lcg(regionSL==j,2)));
        RBSLHR(i,j) = sum(BarySL(regionSL==j,i).*cosd(Lcg(regionSL==j,2)))./sum(cosd(Lcg(regionSL==j,2)));
    end
end

load([Path 'SatelliteAltimetry.mat'])
ALTAc = ALTAc';
ALTAcrel = ALTAc+GIA_Field_KS(:,94:end);
for i = 1:size(ALTAc,2);
    GMSLALT(i,1) = sum(ALTAc(:,i).*cosd(Lcg(:,2)))./sum(cosd(Lcg(:,2)));
    GMSLALTrel(i,1) = sum(ALTAcrel(:,i).*cosd(Lcg(:,2)))./sum(cosd(Lcg(:,2)));
end

NRegion = {'South Atlantic','Indo-Pacific','East Pacific','Subpolar North Atlantic','Subtropical North Atlantic','West Pacific'};
for i = 1:size(ALTAc,2);
    for j = 1:6;
        RMSLALT(i,j) = sum(ALTAc(regionSL==j,i).*cosd(Lcg(regionSL==j,2)))./sum(cosd(Lcg(regionSL==j,2)));
        RMSLALTrel(i,j) = sum(ALTAcrel(regionSL==j,i).*cosd(Lcg(regionSL==j,2)))./sum(cosd(Lcg(regionSL==j,2)));
    end
end

GMSLHR = GMSLHR-mean(GMSLHR(94:end));
GMSLHRrel = GMSLHRrel-mean(GMSLHRrel(94:end));
GGIAHR = GGIAHR-mean(GGIAHR(94:end));
GMSSLHR = GMSSLHR-mean(GMSSLHR(94:end));
GBSLHR = GBSLHR-mean(GBSLHR(94:end,:));
RMSLHR = RMSLHR-mean(RMSLHR(94:end,:));
RMSLHRrel = RMSLHRrel-mean(RMSLHRrel(94:end,:));
RGIAHR = RGIAHR-mean(RGIAHR(94:end,:));
RMSSLHR = RMSSLHR-mean(RMSSLHR(94:end,:));
RBSLHR = RBSLHR-mean(RBSLHR(94:end,:));
GMSLALT = detrend(GMSLALT,'constant');
GMSLALTrel = detrend(GMSLALTrel,'constant');
RMSLALT = detrend(RMSLALT,'constant');
RMSLALTrel = detrend(RMSLALTrel,'constant');

%% (2) CALTCULATE GLOBAL AND REGIONAL MSL AND TRENDS

s = find(t>=1993&t<2022);
X = [t(s) ones(size(t(s)))];
m = X\ALTAcrel';
T(:,1) = m(1,:);
m = X\HRrel(:,s)';
T(:,2) = m(1,:);
TD = T(:,2)-T(:,1);

X = [t(s).^2 t(s) ones(size(t(s)))];
m = X\ALTAcrel';
TACC(:,1) = m(1,:)*2;
m = X\HRrel(:,s)';
TACC(:,2) = m(1,:)*2;
TACCD = TACC(:,2)-TACC(:,1);

[figure3] = Plot_TrendMaps(T,TACC,TD,TACCD,Lcg)

K = NaN(size(ALTAc,1),1);
Kpval = NaN(size(ALTAc,1),1);
for i = 1:size(ALTAc,1);
    [K(i,1), Kpval(i)] = corr(detrend(ALTAc(i,:)'),detrend(HR(i,s)'));
end

for i = 1:size(RMSLHR,2);
    KGB(i,1) = corr(detrend(RMSLALT(:,i)),detrend(RMSLHR(s,i)));
end
KGB(7,1) = corr(detrend(GMSLALT),detrend(GMSLHR(s)));

% Trends & accelerations since 1993
for i = 1:size(RMSLHR,2);
    TGB93HR(i,:) = trend_new(1,t(s),RMSLHR(s,i),'N')*1000;
    TGB93ALT(i,:) = trend_new(1,t(s),RMSLALT(:,i),'N')*1000;

    for j = 1:numel(t);
        SEsdsl(j,:) = normrnd(0,RMSLHRSE(j,i),1,1000)*1000;
    end

    X = [t(s) ones(size(t(s)))];
    m = X\SEsdsl(s,:);
    TGB93HRSE(i,1) = std(m(1,:)');

    TGB93HRACC(i,:) = trend_new(2,t(s),RMSLHR(s,i),'N')*1000;
    TGB93ALTACC(i,:) = trend_new(2,t(s),RMSLALT(:,i),'N')*1000;

    X = [t(s).^2 t(s) ones(size(t(s)))];
    m = X\SEsdsl(s,:);
    TGB93HRACCSE(i,1) = std(2.*m(1,:)');
end

TGB93HR(7,:) = trend_new(1,t(s),GMSLHR(s),'N')*1000;
TGB93ALT(7,:) = trend_new(1,t(s),GMSLALT,'N')*1000;
for j = 1:numel(t);
    SEsdsl(j,:) = normrnd(0,GMSLHRSE(j,1),1,1000)*1000;
end
X = [t(s) ones(size(t(s)))];
m = X\SEsdsl(s,:);
TGB93HRSE(7,1) = std(m(1,:)');

TGB93HRACC(7,:) = trend_new(2,t(s),GMSLHR(s),'N')*1000;
TGB93ALTACC(7,:) = trend_new(2,t(s),GMSLALT,'N')*1000;
for j = 1:numel(t);
    SEsdsl(j,:) = normrnd(0,GMSLHRSE(j,1),1,1000)*1000;
end
X = [t(s).^2 t(s) ones(size(t(s)))];
m = X\SEsdsl(s,:);
TGB93HRACCSE(7,1) = std(m(1,:)');

TGB93HRSEcomb = sqrt(TGB93HRSE.^2+TGB93HR(:,3).^2)
TGB93HRACCSEcomb = sqrt(TGB93HRACCSE.^2+TGB93HRACC(:,3).^2)

% Plots 
[figure2] = Plot_GlobalandRegional_SeaLevel(t,GMSLHRrel,GMSLALTrel,...
    RMSLHRrel,RMSLALTrel,GMSSLHR,GMSSLHRSE,GBSLHR,GBSLHRSE,GGIAHR,RMSSLHR,...
    RMSSLHRSE,RBSLHR,RBSLHRSE,RGIAHR,GGIAHRSE,RGIAHRSE,Lcg,GMSLHRSE,...
    RMSLHRSE,NRegion,K)
load('GMSL_Palmer_1900-2010.mat')
GMSL_Palmer(:,2:3) = GMSL_Palmer(:,2:3)/1000;
GMSL_PastRec(:,2:end) = GMSL_PastRec(:,2:end)/1000;
[figure7] = Plot_GMSL_PastRec(t,GMSLHR,GMSLHRSE,GMSL_Palmer,GMSL_PastRec)

% Trends since 1900
for i = 1:6;
    TSL(i,:) = trend_new(1,t,RMSLHRrel(:,i),'N')*1000;

    for j = 1:numel(t);
        SEsl(j,:) = normrnd(0,RMSLHRSE(j,i),1,1000)*1000;
    end

    X = [t ones(size(t))];
    m = X\SEsl;
    TSLSE(i,1) = std(m(1,:)');

    TSDSL(i,:) = trend_new(1,t,RMSSLHR(:,i),'N')*1000;

    for j = 1:numel(t);
        SEsdsl(j,:) = normrnd(0,RMSSLHRSE(j,i),1,1000)*1000;
    end

    X = [t ones(size(t))];
    m = X\SEsdsl;
    TSDSLSE(i,1) = std(m(1,:)');

    TBSL(i,:) = trend_new(1,t,RBSLHR(:,i),'N')*1000;

    for j = 1:numel(t);
        SEbsl(j,:) = normrnd(0,RMSSLHRSE(j,i),1,1000)*1000;
    end

    X = [t ones(size(t))];
    m = X\SEbsl;
    TBSLSE(i,1) = std(m(1,:)');

    TGIA(i,:) = trend_new(1,t,RGIAHR(:,i),'N')*1000;
    TGIASE(i,1) = RGIAHRwmSE(i,1)*1000;
end
TSL(7,:) = trend_new(1,t,GMSLHRrel,'N')*1000;
for j = 1:numel(t);
    SEsl(j,:) = normrnd(0,GMSLHRSE(j,1),1,1000)*1000;
end
X = [t ones(size(t))];
m = X\SEsl;
TSLSE(7,1) = std(m(1,:)');
TSDSL(7,:) = trend_new(1,t,GMSSLHR,'N')*1000;
for j = 1:numel(t);
    SEsl(j,:) = normrnd(0,GMSSLHRSE(j,1),1,1000)*1000;
end
X = [t ones(size(t))];
m = X\SEsl;
TSDSLSE(7,1) = std(m(1,:)');
TBSL(7,:) = trend_new(1,t,GBSLHR,'N')*1000;
for j = 1:numel(t);
    SEsl(j,:) = normrnd(0,GBSLHRSE(j,1),1,1000)*1000;
end
X = [t ones(size(t))];
m = X\SEsl;
TBSLSE(7,1) = std(m(1,:)');
TGIA(7,:) = trend_new(1,t,abs(GGIAHR),'N')*1000;
TGIASE(7,1) = GGIAHRwmSE*1000;

TSLSEcomb = sqrt(TSLSE.^2+TSL(:,3).^2);
TSDSLSEcomb = sqrt(TSDSLSE.^2+TSDSL(:,3).^2);
TBSLSEcomb = sqrt(TBSLSE.^2+TBSL(:,3).^2);
TGIASEcomb = sqrt(TGIASE.^2+TGIA(:,3).^2);

%% (3) STERIC VALIDATION

name = {'South Atlantic','Indian Ocean - South Pacific','East Pacific','Subpolar North Atlantic','Subtropical North Atlantic','Northwest Pacific'}
for i = 1:6;
    [RMSSL_F20(:,:,i)] = Read_F20_Basin(name{i});
end
RMSSL_F20 = RMSSL_F20/1000;

for i = 1:6;
    isn = find(~isnan(RMSSL_F20(:,2,i)));
    KST(i,1) = corr(detrend(RMSSLHR(isn,i)),detrend(RMSSL_F20(isn,2,i)));
    KST(i,2) = corr(detrend(RMSSLHR(71:119,i)),detrend(RMSSL_F20(71:119,2,i)));
    TSDKS(i,:) = trend_new(1,t(isn),RMSSLHR(isn,i),'N')*1000;
    TST(i,:) = trend_new(1,t(isn),RMSSL_F20(isn,2,i),'N')*1000;
    
    for j = 1:numel(t);
        SEsdsl(j,:) = normrnd(0,RMSSLHRSE(j,i),1,1000)*1000;
    end

    X = [t(isn) ones(size(t(isn)))];
    m = X\SEsdsl(isn,:);
    TSTKSSE(i,1) = std(m(1,:)');
end
TSDKS(7,:) = trend_new(1,t(isn),GMSSLHR(isn),'N')*1000;
for j = 1:numel(t);
    SEsdsl(j,:) = normrnd(0,GMSSLHRSE(j,1),1,1000)*1000;
end
X = [t(isn) ones(size(t(isn)))];
m = X\SEsdsl(isn,:);
TSTKSSE(7,1) = std(m(1,:)');

TSDKSSEcomb = sqrt(TSTKSSE.^2+TSDKS(:,3).^2)

[figure4, Steric_F20] = Plot_Steric_Comparison(t,GMSSLHR,GMSSLHRSE,RMSSLHR,RMSSLHRSE,RMSSL_F20,NRegion);
TST(7,:) = trend_new(1,t(isn),median(Steric_F20(isn,:),2),'N')*1000;

%% (4) RATES AND REGIONAL BUDGET

% Apply SSA
sm = 15;
yfit = .001*t;yfit = yfit-yfit(end); % SSA sometimes becomes unstable in 
% time series with small trends; adding a trend component (which is removed 
% afterwards again) helps
for i = 1:6;
    [RMSLHRsm(:,i),~,~,~,~,RMSLHRsmEns(:,:,i)] = ssatrend(RMSLHR(:,i),sm,RMSLHRSE(:,i));
    [RMSLALTsm(:,i),~,~,~,~,~] = ssatrend([RMSLHR(1:93,i);RMSLALT(:,i)],sm,RMSLHRSE(1:end,i));
    [RMSSLHRsm(:,i),~,~,~,~,RMSSLHRsmEns(:,:,i)] = ssatrend(RMSSLHR(:,i)+yfit,sm,RMSSLHRSE(:,i));
    [RBSLHRsm(:,i),~,~,~,~,RBSLHRsmEns(:,:,i)] = ssatrend(RBSLHR(:,i),sm,RBSLHRSE(:,i));
end
RMSSLHRsm = RMSSLHRsm-yfit;
RMSSLHRsmEns = RMSSLHRsmEns-yfit;

[GMSLHRsm,~,~,~,~,GMSLHRsmEns] = ssatrend(GMSLHR,sm,GMSLHRSE);
[GMSLALTsm,~,~,~,~,~] = ssatrend([GMSLHR(1:93);GMSLALT],sm,GMSLHRSE(1:end));
[GMSSLHRsm,~,~,~,~,GMSSLHRsmEns] = ssatrend(GMSSLHR,sm,GMSSLHRSE);
[GBSLHRsm,~,~,~,~,GBSLHRsmEns] = ssatrend(GBSLHR,sm,GBSLHRSE);

[figure9] = Plot_GlobalRates(t,GMSLHRsm,GMSSLHRsm,GBSLHRsm,GMSLHRsmEns,GMSSLHRsmEns,GBSLHRsmEns,GMSLALTsm)

% Linear Trends and Accelerations in regional sea levek

s = find(t>=1900);
X = [t(s) ones(size(t(s)))];
m = X\HRrel(:,s)';
TBud(:,1) = m(1,:);
m = X\SDSL(:,s)';
TBud(:,2) = m(1,:);
m = X\BarySL(:,s)';
TBud(:,3) = m(1,:);
m = X\GIA_Field_KS(:,s)';
TBud(:,4) = m(1,:);
m = X\IBAint(:,s)';
TBud(:,5) = m(1,:);

s = find(t>=1970);
X = [t(s).^2 t(s) ones(size(t(s)))];
m = X\HRrel(:,s)';
TACCBud(:,1) = m(1,:)*2;
m = X\SDSL(:,s)';
TACCBud(:,2) = m(1,:)*2;
m = X\BarySL(:,s)';
TACCBud(:,3) = m(1,:)*2;
m = X\GIA_Field_KS(:,s)';
TACCBud(:,4) = m(1,:)*2;
m = X\IBAint(:,s)';
TACCBud(:,5) = m(1,:)*2;

[figure10, TACC_GMSLPastRec] = Plot_Acceleration_1970(t,GMSL_Palmer,GMSLHR,GMSSLHR,GBSLHR,GMSLHRSE,GMSSLHRSE,GBSLHRSE,GMSLALT)

%% (5) TIDE GAUGE COMPARISON

load('tide gauge data.mat')
M = tg.RLRdata(:,21:end)'; L = [tg.lon,tg.lat]; MID = tg.ID; N = tg.N;
regionTG = tg.region;

sF = find(MID==662);
M(:,sF) = []; MID(sF) = []; N(sF) = []; %SF(sF) = [];
L(sF,:) = []; %Flag(:,sF) = [];

s = find(t>=1993&t<2022);
idx = knnsearch(Lcg,L);
MHR = HR(idx,:)';
MSE = SE(idx,:)';
MIBE = IBAint(idx,:)';
MSDSLSE = SE_SDSL(idx,:)';
MBarySE = SE_Bary(idx,:)';
MSDSL = SDSL(idx,:)';
MSDSL = MSDSL-mean(MSDSL(s,:));
MBary = BarySL(idx,:)';
MBary = MBary-mean(MBary(s,:));
% MHRrel = HR(idx,:)'+VLM_KS';
MGIA = GIA_Field_KS(idx,:)';
MGIASE = SE_GIA(idx,:)';
MHRrel = HRrel(idx,:)';
regionTG = regionSL(idx);
for i = 1:size(M,2);
    isn = find(~isnan(M(:,i)));
    MHR(:,i) = MHR(:,i)-mean(MHR(isn,i))+mean(M(isn,i));
    MHRrel(:,i) = MHRrel(:,i)-mean(MHRrel(isn,i))+mean(M(isn,i));
end
MALTrel = ALTAcrel(idx,:)';
MALTrel = detrend(MALTrel,'constant')+mean(MHRrel(s,:));
MALT = ALTAc(idx,:)';
MALT = detrend(MALT,'constant')+mean(MHR(s,:));

[figure1] = Plot_DataAvailability_Map(t,M,L,regionTG,LALT,regionSL)

T1993 = NaN(size(M,2),3,3);
for i = 1:size(M,2);
    isn = find(~isnan(M(s,i)));
    if numel(isn)>=15;
        T1993(i,:,1) = trend_new(1,t(s(isn)),M(s(isn),i),'N');
        T1993(i,:,2) = trend_new(1,t(s(isn)),MALTrel(isn,i),'N');
        T1993(i,:,3) = trend_new(1,t(s(isn)),MHRrel(s(isn),i),'N');
    end
end

[ha] = Plot_Acceleration_LocalRates(LALT,TACCBud,M,t,sF,HRrel,SDSL,BarySL,IBAint,GIA_Field_KS,SE,SE_SDSL,SE_Bary,SE_GIA)

KTG = NaN(size(M,2),2);
for i = 1:size(M,2);
    progressbar(i/size(M,2));
    isn = find(~isnan(M(:,i)));
    KTG(i,1) = corr(detrend(M(isn,i)),detrend(MHR(isn,i)));   
    KTG(i,2) = corr(detrend(MHR(s,i)),detrend(MALT(:,i)));
end

% GPS SONEL
load('SONEL_LINEAR RATES.mat')
GPS_SONEL = table2array(GPS_Table);
idxGPS = knnsearch(GPS_SONEL(:,1:2),L);
GPS_SONEL = GPS_SONEL(idxGPS,:);

load('BSLRad_Field_KS.mat'); % solid earth component of BSL
X = [t(s) ones(size(t(s)))];
m = X\BSLRad_Field_KS';
RAD = m(1,idx)';

%GIA geoid correction for comparison to GPS
F = scatteredInterpolant(LALT(:,1),LALT(:,2),(GIAGEOID_Field_KS(:,end)-GIAGEOID_Field_KS(:,1))/size(GIAGEOID_Field_KS,2),'linear','linear');
GIA5G = F(L(:,1),L(:,2));

wgs84 = wgs84Ellipsoid("km");
for i = 1:size(M,2);
    d(i,1) = distance(L(i,2),L(i,1),GPS_SONEL(i,2),GPS_SONEL(i,1),wgs84);
end
VLM_GPS = (GPS_SONEL(:,5)./1000-RAD)'.*t;
VLM_GPS = VLM_GPS-VLM_GPS(end,:);
GIA5G = (GIA5G-.25/1000)'.*t;
GIA5G = GIA5G-GIA5G(end,:);
MHRrel_SONEL = MHR-VLM_GPS+GIA5G;
for i = 1:size(M,2);
    isn = find(~isnan(M(:,i)));
    MHRrel_SONEL(:,i) = MHRrel_SONEL(:,i)-mean(MHRrel_SONEL(isn,i))+mean(M(isn,i));
end
sd = find(d<=8);

% Trend comparison

Tobs = NaN(size(M,2),3);
THRrel = NaN(size(M,2),3);
THRrel_SONEL = NaN(size(M,2),3);
Tobsfull = NaN(size(M,2),3);
THRrelfull = NaN(size(M,2),3);
THRrelfull_SONEL = NaN(size(M,2),3);
s = find(t>1992);
for i = 1:size(M,2);
    isn = find(~isnan(M(s,i)));
    if numel(isn)>=.75*numel(s);
        [Tobs(i,:)] = trend_new(1,t(s(isn)),M(s(isn),i),'N');
        [THRrel(i,:)] = trend_new(1,t(s(isn)),MHRrel(s(isn),i),'N');
        [THRrel_SONEL(i,:)] = trend_new(1,t(s(isn)),MHRrel_SONEL(s(isn),i),'N');
    end
    isn = find(~isnan(M(:,i)));
    if numel(isn)>=.75*numel(s);
        [Tobsfull(i,:)] = trend_new(1,t(isn),M(isn,i),'N');
        [THRrelfull(i,:)] = trend_new(1,t(isn),MHRrel(isn,i),'N');
        [THRrelfull_SONEL(i,:)] = trend_new(1,t(isn),MHRrel_SONEL(isn,i),'N');
    end
end

TD1 = (Tobs(:,1)-THRrel(:,1))*1000;
TD2 = (Tobs(:,1)-(THRrel_SONEL(:,1)))*1000;

[figure6] = Plot_ReconstructionTrendBias(L, TD1, TD2, sd)

[figure5] = Plot_LocalBudgetsandErrors(L,Tobs,THRrel,M,MHRrel,MSE,MSDSL,MSDSLSE,MBary,MBarySE,MGIA,MGIASE,MIBE,MALTrel,N,t)
