%Transport_calculation_frame.m

% Calculate the longshore transport from bottom frame data, 
% interpolation and extrapolation
% Steps:
% 1)    Interpolate velocity in bottom boundary layer (up to 2.8 mab)
% 2)a    Interpolate SPM to velocity-grid
% 2)b    Extrapolate SPM uniform to seabed
% 3)    transport per (known) z-layer
% 4)a   transport per s-layer
% 4)b   extrapolate transport uniform to surface 

clear
close all
cd('D:\Data\LanderLT\ADCP\Data\mlab\transport')
load('D:\Data\LanderLT\ADCP\Data\mlab\processing\ADCPS')
load('D:\Data\LanderLT\ADCP\Data\mlab\processing\ADCPD')
load('D:\Data\LanderLT\ADCP\Data\mlab\processing\ADCPVf')
load('D:\Data\LanderLT\MEEALL\mlab\dataset\burst2.mat')
load('D:\Data\LanderLT\MEEALL\mlab\dataset\sallspmchla');


load('D:\Data\LanderLT\MEEALL\mlab\dataset\landerDAY') 
[nd,~] = size(ADCPS.sac);  %nz is used for Navicula observations
nt = DAY.recov(2,16);

col1 = erosed;

%% prepare
% constants
pr      = 1:20;     % check-profile part
outp    = 0:0.25:7; % bottom boundary layer size
nl      = 40;       % sigma-layers
    
% pre-allocate
zf          = NaN(11+nd,nt);
velipcf     = NaN(11+nd,nt);
velvvlf     = NaN(11+nd,nt);
% velvvcf     = NaN(11+nd,nt);
spmi        = NaN(11+nd,nt);
sz          = NaN(nl+1,nt);
sh          = NaN(1,nt);

ttfl        = NaN(11+nd,nt);
ttfc        = NaN(11+nd,nt);
trfl        = NaN(nl+1,nt);
trfc        = NaN(nl+1,nt);

%% velocities
% idvc = [];
% idvl = [];
vels = ADCPV.vl;% + ADCPV.vc;

for i = 1:nt %10000%1: 20000% DAY.deploy(2,14):DAY.recov(2,14)%
    % determine if vel+bs-data available, else skip timestep
    if length(find(~isnan(vels(pr,i)))) > 10
%         if ~isnan(nanmean(ADCPS.sac(:,i)))
    
%% 1) Interpolate velocity in bottom boundary layer (up to 2.8 mab)
% calculate interpolation and store per timestep
    % 1) cross-shore: cubic interpolation
    % check NaN in velocity profile
    pr1 = find(~isnan(vels(pr,i)));
    kc(i) = pr1(1)-1;
    profc = ADCPV.vc(pr1,i);
    profc(1) = nanmean(profc(1:4));
    profc(end) = nanmean(profc(end-4:end));
    ws = 4;
%     if length(pr1) > 9
% %                     lastwarn('')            %clearing the last warning
%     velc(pr1,i) = filtfilt(ones(1,ws)/ws,1,profc);
%     intpcc = interp1([0;0.15;ADCPD.bin(pr1)],[0;burst.uc(i);velc(pr1,i)],outp,'pchip'); % 'spline'
% %                     msgc = lastwarn;
% %                     if ~isempty(msgc)
% %                         idvc = [idvc; i];
% %                     end
%     velipcf(:,i)    = [intpcc(1:11+kc(i))';ADCPV.vc(1+kc(i):end,i)];
%     end
    % longshore:  van veen
%     pr1 = find(~isnan(ADCPV.vl(pr,i)));
%     kl(i) = pr1(1)-1;
%                     lastwarn('')            %clearing the lastwarn
%     al = robustfit(([0.15;ADCPD.bin(1:10)]).^(1/5.2),[burst.ul(i); ADCPV.vl(1:10,i)],'bisquare',4.685,'off');
    al = robustfit(([0.15;ADCPD.bin(pr1(1:8))]).^(1/5.2),[burst.ul(i); ADCPV.vl(pr1(1:8),i)],'bisquare',4.685,'off');
%                     msgc = lastwarn;
%                     if ~isempty(msgc)
%                         idvl = [idvl; i];
%                     end
    intpvl = al*outp'.^(1/5.2); 
% final= 
zf(:,i)         = [outp(1:11+kc(i))'; ADCPD.bin(1+kc(i):end)];
velvvlf(:,i)    = [intpvl(1:11+kc(i));ADCPV.vl(1+kc(i):end,i)];
%         end
    end

end
save('vel_field','zf','velvvlf','velipcf')        

%%       
% idtrfc = [];
% idtrfl = [];
ik = [];
% load('vel_field')
for i = 1:nt %10000%1: 20000% DAY.deploy(2,14):DAY.recov(2,14)%
    % determine if vel+bs-data available, else skip timestep
    if length(find(~isnan(vels(pr,i)))) > 10
%         if ~isnan(nanmean(ADCPS.sac(:,i)))
        ik = [ik; i];
%% 2) intra-/extrapolate SPM
% lineair between known points  (interpolation)
% uniform towards bed           (extrapolation)
% nothing towards surface
rngs = find(~isnan(ADCPS.sac(2:end,i)))+1;
if isempty(rngs)                %no ADCP backscatter data available
    spmi(:,i) = interp1([0;0.3;0.8;1.4;2.0],...
    [sspm(i,1);sspm(i,:)'],...  
    zf(:,i)); % linear    
else
%                     lastwarn('')            %clearing the last warning
    spmi(:,i) = interp1([0;0.3;0.8;1.4;2.0;ADCPD.bin(rngs)],...
    [sspm(i,1);sspm(i,:)';ADCPS.sac(rngs,i)],...  
    zf(:,i)); % linear
end
%                     msgc = lastwarn;
%                     if ~isempty(msgc)
%                         idspm = [idspm; i];
%                     end

%% 3) transport per z-layer
% frame
ttfl(:,i) = spmi(:,i).*velvvlf(:,i);
ttfc(:,i) = spmi(:,i).*velipcf(:,i);

%% 4) interpolate to sigma-grid 
% 40 layers
% depth from ADCPD.dep

sh(i) = ADCPD.dep(i)/nl;
sz(:,i) = 0:sh(i):ADCPD.dep(i);
%                     lastwarn('')            %clearing the last warning
rngl = find(~isnan(ttfl(:,i)));
if ~isempty(rngl)
trfl(:,i) = interp1(zf(rngl,i),ttfl(rngl,i),sz(:,i));
%                     msgc = lastwarn;
%                     if ~isempty(msgc)
%                         idtrfl = [idtrfl; i];
%                     end
%                     lastwarn('')            %clearing the last warning
end
rngc = find(~isnan(ttfc(:,i)));
if~isempty(rngc)
trfc(:,i) = interp1(zf(rngc,i),ttfc(rngc,i),sz(:,i));
%                     msgc = lastwarn;
%                     if ~isempty(msgc)
%                         idtrfc = [idtrfc; i];
%                     end
end


%         end
    end
end

%% 5) extrapolate to surface
for i = ik'
% extrapolate to surface (if data up to xx% of the water column) 
ii = find(~isnan(trfl(:,i)),1,'last');
if ii > 30
trfl(ii+1:36,i) = nanmean(trfl(ii-1:ii,i));
trfc(ii+1:36,i) = nanmean(trfc(ii-1:ii,i));
end
end

%% figure
figure;  subplot(2,1,1); imagesc(DAY.count(1:nt),1:41,trfl); set(gca,'YDir','normal'); clim([-100 100]); plot_deploylines(2,38,'w'); xlim([DAY.deploy(2,2) DAY.recov(2,9)]);  %2:9 // 11:16
colormap(col1);  title('longshore transport'); %datetick('keeplimits');


subplot(2,1,2); plot(DAY.count(1:nt),burst.twf(1:nt),'.'); plot_deploylines(2,30); xlim([DAY.deploy(2,2) DAY.recov(2,9)]); %datetick('keeplimits'); 
title('wave shear stress'); xlabel('Date in 2011')

%subplot(3,1,1);
% wd = 4; fsal = filtfilt(ones(1,wd)/wd,1,diff(CTD.sal(1:nt)))
% subplot(3,1,2); plot(DAY.count(1:nt-1),fsal)
% subplot(3,1,2); plot(DAY.count(1:nt - 1),diff(CTD.sal(1:nt)))
% subplot(3,1,3); plot(burst.tw(1:nt))


%%
% save('transport_vars','zf','velipcf','velvvlf','spmi','ttfl','ttfc','trfl','trfc')
% load('transport_vars')

%% fill slack periods
load('D:\Data\LanderLT\MEEALL\mlab\dataset\T12.mat')

rngpers = [];
for i = [2:9,11:16]
    if i == 9
            rngpers = [rngpers DAY.deploy(2,i):33254];
            
    elseif i == 15
            rngpers = [rngpers DAY.deploy(2,i):67671];
    else
        rngpers = [rngpers DAY.deploy(2,i):DAY.recov(2,i)];
    end
    
end

hloc = NaN(1,nt);
for i = 1:nt
    if ~isnan(nanmean(trfl(:,i)))
        hloc(i) = find(~isnan(trfl(:,i)),1,'last');

    end
end


rngslack1 = intersect([T12.st T12.up],rngpers);

for k = 2:length(rngslack1)
    [~, rngslack2(k)] = nanmin(abs(trfl(4,rngslack1(k)-8:rngslack1(k)+8)));
    rngslack2(k) = rngslack2(k)+rngslack1(k)-9;
end



trfli = trfl;
rngsl = unique(rngslack2); % remove double entries due to non-reversal
for k = 2:length(rngsl)
    rng = rngsl(k)-1:rngsl(k)+1;
    for l = rng
        if hloc(l) < 30
            trfli(hloc(l):30,l)   = trfli(hloc(l)-1,l);
        end
    end
end

        
% figure; imagesc(trfli)
% set(gca,'YDir','normal'); clim([-100 100]); plot_deploylines(2,38,'w');colormap(col1);  xlim([DAY.deploy(2,2) DAY.recov(2,9)]);  %2:9 // 11:16

%% interpolate timewise < 1/4 tide

% trfli1 = trfli;
% trfli2a = trfli;
trfli2b = trfli;
trfli2 = trfli;
% per layer interpolate
for l = 2:36
    intnan = [];
    nanl = find(isnan(trfli(l,:)));
    for k = 1:length(nanl)
        ki = find(~isnan(trfli(l,1:nanl(k))),1,'last');
        kj = find(~isnan(trfli(l,nanl(k):end)),1,'first')+nanl(k)-1;
        if kj-ki < 10
            intnan = [intnan; nanl(k)];  %  
        end
    end
    nonnanl = setxor(1:nt,nanl)';
    filll = intersect(nanl,intnan);
    intl  = sort([nonnanl; filll]);
%     trfli2a(l,intl) = interp1(DAY.count(nonnanl),trfli2(l,nonnanl),DAY.count(intl),'pchip'); %
    trfli2b(l,filll) = trfli2b(l-1,filll);
    for m = 1:numel(filll)
    trfli2(l,filll(m)) = closest(0,[trfli2a(l,filll(m)); trfli2b(l,filll(m))]);  % min(hor-int vert-extr)
    end
    
end

% figure; plot(trfli(5,:),'k'); hold on; plot_deploylines(2,0); plot(trfli1(25,:)); plot(DAY.count(filll),trfli1(25,filll),'*');  plot(DAY.count(filll),trfli2(25,filll),'*')

% extrapolate uniform remaining
trfli3 = trfli2;
hloc = NaN(1,nt);
for i = 1:nt
    if ~isnan(nanmean(trfli2(:,i)))
        hloc(i) = find(~isnan(trfli2(:,i)),1,'last');
        trfli3(hloc(i)+1:end,i) = closest(0,trfli2(hloc(i)-2:hloc(i),i));
    end
end




save('trfl','trfl','trfli','trfli2a','trfli2b','trfli3')

%% 
nrs = [2:9,11:16]; %17?
% trfli = trfl; % back-up

% ij = 2;%nrs;
% rng = (DAY.deploy(2,ij)):DAY.recov(2,ij);
% trfli(trfli == 0) = NaN;

for ij = nrs% 15:18;%nrs; %2;%
    rng = (DAY.deploy(2,ij):DAY.recov(2,ij));
    F = figure;
    subplot(3,1,1)
    pcolor(DAY.num(rng),1:41,trfli3(:,rng)), shading flat
    colormap(col1)
    datetick
    clim([-100 100])
    subplot(3,1,2)
    plot(DAY.num(rng),burst.tmax(rng),'.')
    datetick
    subplot(3,1,3) % plot correlation and R2
    plot(DAY.num(rng),trfli3(35,rng)./trfli3(30,rng))
    hold on
    plot(DAY.num(rng),trfli3(30,rng)./trfli3(25,rng)+2)
     plot(DAY.num(rng),trfli3(25,rng)./trfli3(20,rng)+4)
    plot(DAY.num(rng),trfli3(20,rng)./trfli3(15,rng)+6)
    plot(DAY.num(rng),trfli3(15,rng)./trfli3(10,rng)+8)
    plot(DAY.num(rng),trfli3(10,rng)./trfli3(5,rng)+10)
    ylim([-2 14])
    legend('35/30','30/25','25/20','20/15','15/10','10/05')
    hline(0,'k:'),hline(2,'k:'),hline(4,'k:'),hline(6,'k:'),hline(8,'k:'), hline(10,'k:')
    datetick
    set(F,'PaperUnits', 'inches','PaperPosition', [0 0 15 10])
    print(F,'-dtiff','-r300',['transport_2D_per' num2str(ij) '.tif']);
    close

    % first de-deselect entire NaN-time-frames (how come both OBS and ADCP?)
    % How about when interpolation occurs across the non-selected NaN-frames
%     desel = find(nanmean(trfli(:,rng))==0);
%     sel = setxor(rng(desel),rng);


%     for il = 1:41       % per layer
% %         selnan = find(isnan(trfli(il,sel)));
% %         selnum = setxor(sel(selnan),sel);
% %         xu = find(diff(selnan)==1);
% %         xv = unique([xu;xu+1]);
% %         sel2 = setxor(1:length(selnan),xv);
% %         trfli(il,sel(selnan(sel2))) = interp1(DAY.num(selnum), trfli(il,selnum),DAY.num(sel(selnan(sel2))));
% 
% % sx = diff(find(~isnan(trfli(il,rng)))); % start and end-profiles to exclude?
%         rx = find(~isnan(trfli(il,rng)));
%         sx = diff(rx); % start and end-profiles to exclude? no (apparently for per2)
%         jj = find(sx>12);
%         kk = [];
%         for kj = 1:length(jj)
%             kk = [kk; (rx(jj(kj))+1:rx(jj(kj)+1)-1)'];
%         end
%         tx = setxor(rng,rng(kk));
%         
%         trfli(il,tx) = interp1(DAY.num(rng(rx)), trfli(il,rng(rx)),DAY.num(tx'));
%     end
    
%     % extrapolate remaining profiles to surface
%     for i = find(isnan(trfli(41,rng)))
%         ii = find(~isnan(trfli(:,rng(i))),1,'last');
%         trfli(ii+1:end,rng(i)) = nanmean(trfli(ii-6:ii,rng(i)));
% %         trfc(ii+1:end,i) = nanmean(trfc(ii-6:ii,i));
%     end

    G = figure;
    subplot(3,1,1)
    plot(rng,trfli3(10,rng)); hold on; plot(rng,trfli3(20,rng)); plot(rng,trfli3(30,rng)); plot(rng,trfli3(40,rng)); legend('10','20','30','40')
    hline(0,'k:')
    ylim([-100 100])
    subplot(3,1,2)
    plot(rng,trfli3(10,rng)-trfli3(20,rng)); hold on;  plot(rng,trfli3(20,rng)-trfli3(30,rng)); plot(rng,trfli3(30,rng)-trfli3(40,rng)); legend('10-20','20-30','30-40')
    hline(0,'k:')
    ylim([-50 50])
    subplot(3,1,3)
    plot(rng,burst.tmax(rng),'.')
    set(G,'PaperUnits', 'inches','PaperPosition', [0 0 7 3])
    print(G,'-dtiff','-r300',['transport_layers_per' num2str(ij) '.tif']);
    close
end


% save('trf','trfl','trfc')
%% per tidal cycle

load('D:\Data\LanderLT\MEEALL\mlab\dataset\T12')

%rngpers at line 196
rngtr = find(~isnan(mean(trfli3)));
st = intersect(T12.st,rngtr); st = st(1:end-1);
en = intersect(T12.en,rngtr); en = en(2:end);
% up = intersect(T12.up,rngpers);
k=0;
for l = 1:length(st)
    if ~isnan(mean(trfli3(4,[st(l)-4:en(l)+4])))
        k = k+1;
    [~, indst(k)] = nanmin(abs(trfli3(4,st(l)-8:st(l)+8)));
     T12tr.st(k) = indst(k)+st(l)-9;
    [~, inden(k)] = nanmin(abs(trfli3(4,en(l)-8:en(l)+8)));
     T12tr.en(k) = inden(k)+en(l)-9;
     [~, indup(k)] = nanmin(abs(trfli3(4,T12tr.st(k)+11:T12tr.en(k)-11)));
     T12tr.up(k) = indup(k)+T12tr.st(k)-1+11;
    end
%     [~, indup(k)] = nanmin(abs(trfl(4,up(k)-8:up(k)+8)));
%      T12tr.up(k-1) = indup(k)+up(k)-9;
end

% save('T12tr','T12tr')


%%

for ij = 1:numel(T12tr.st)%80:229 % periods 2-4
        trng = T12tr.st(ij):T12tr.en(ij);
        ttl(ij) = mean(mean(trfli3(:,trng)));
        ttlu(ij) = mean(mean(trfli3(21:41,trng)));
        ttll(ij) = mean(mean(trfli3(1:21,trng)));
        ttl1(ij) = mean(mean(trfli3(1:11,trng)));
        ttl2(ij) = mean(mean(trfli3(11:21,trng)));
        ttl3(ij) = mean(mean(trfli3(21:31,trng)));
        ttl4(ij) = mean(mean(trfli3(31:41,trng)));
%         indbn = trng(mean(trfli3(:,trng))>0);
%         indbs = setxor(indbn,trng);
        trngn = T12tr.up(ij):T12tr.en(ij);
        trngs = T12tr.st(ij):T12tr.up(ij);
        ttlbn(ij) = mean(mean(trfli3(:,trngn)));
        ttlbs(ij) = mean(mean(trfli3(:,trngs)));
        
    
end
%%
    rngp = 1:numel(T12tr.st);
%     rngpd = T12.st(kk):T12.en(kk);
    figure

subplot(3,1,1)
plot(DAY.num(T12tr.st(rngp)),ttl1(rngp),'k:.')
hold on
% plot(DAY.num(T12.st(rngp)),ttlu(rngp),'k:')
% plot(DAY.num(T12.st(rngp)),ttll(rngp),'k--')
plot(DAY.num(T12tr.st(rngp)),ttl2(rngp),'k--.')
plot(DAY.num(T12tr.st(rngp)),ttl3(rngp),'r:.')
plot(DAY.num(T12tr.st(rngp)),ttl4(rngp),'r--.')
% plot(DAY.num(T12.st(rngp)),ttl(rngp),'k')
xlim([DAY.deploy(1,11) DAY.recov(1,16)])   % 2  9 //  11 16
datetick('keeplimits')
grid on
legend('1','2','3','4')%'UWL','LWL') 'total',
ylabel('net transport')
title('Longshore transport -- 2012') %2011 2012

subplot(3,1,2)
plot(DAY.num(T12tr.st(rngp)),ttlbn(rngp),'r-')
hold on
plot(DAY.num(T12tr.st(rngp)),abs(ttlbs(rngp)),'s-')
xlim([DAY.deploy(1,11) DAY.recov(1,16)])   % 2  9 //  11 16
datetick('keeplimits')
grid on
legend('north','south')
ylabel('gross transport')

subplot(3,1,3)
plot(DAY.num(T12.st),T12.Hsavg)
xlim([DAY.deploy(1,11) DAY.recov(1,16)])   % 2  9 //  11 16
grid on
ylabel('mean Hs')
datetick('keeplimits')
ylim([0 5])
    
% EOF