% ProcED II - July, 2010
% This routine is designed to facilitate the handling of a large 
% amount of data derived from ADCP monitoring of estuarine flows in a 
% tidal cycle. The routine makes it possible to eliminate bad ensembles, 
% to extrapolate current velocities  for the surface and bottom blanked areas 
% and to smooth noisy cross profiles, enhancing the visualization of the 
% velocity field. By performing  space and time correlation between cell 
% velocities, ProcED allows for the computation of the residual velocity 
% field, both in the entire cross section and in discrete vertical profiles.
%
% The left margin of the channel (seaward facing) is taken as a reference for 
% cross-section orientation, whereas positive (negative) velocity values 
% indicate seaward (landward) directed flows.
%
% The following routine files must be in the same directory as ProcED_II
% BlWhRe_colormap, rot.m, smartmean.m, isoline.m, filter_5.m
%
% for full information see: 
% Genz,F., Cirano, M., Lessa, G.C. 2010.ProcED: a MATLAB package for 
% Processing ADCP Estuarine Data. Revista Brasileira de Geofísica
% 28(2):183-192.
%
% for information on how to run, data format and initial steps, please
% refer to the manual: "ProcED_II_user_Manual.doc"
%
% - version II - upgraded by Ricardo M. Domingues , ricardormd@gmail.com

% looping indexes:
% i - depth 
% j - cross section
% k - time
%

clc;clear all; close all;
disp(' ')
disp(' ')
disp('          **********************************************************')
disp('          *                       ProcED                           *')
disp('          *                                                        *')
disp('          *        -- tidal cycle ADCP data processing tool --     *')
disp('          **********************************************************')
disp(' ')
disp(' ')
%--------------------------------------------------------------------------
%Processing initial information
dir_path = input('\nEnter the full path of data directory \n ==> ','s');
%Change the current directory to the input path
cd(dir_path);

clear name naminfo;

%---ProcED_processing_options file--------------
%           
%   - refer to the manual
naminfo=dir('ProcED_processing_options.txt');
if isempty(naminfo),
  fprintf('ERROR******* Can''t find file ProcED_processing_options.txt');
  return;
end;
fid_Op=fopen('ProcED_processing_options.txt');
fgetl(fid_Op);files_header=fgetl(fid_Op);
fgetl(fid_Op);label=fgetl(fid_Op);
fgetl(fid_Op);flag=fscanf(fid_Op,'%f');
fgetl(fid_Op);ini_margin=fscanf(fid_Op,'%f');
fgetl(fid_Op);dist_inicio=fscanf(fid_Op,'%f');
fgetl(fid_Op);ang_manual=fscanf(fid_Op,'%f');
fgetl(fid_Op);ctds_inf=fscanf(fid_Op,'%f');
fgetl(fid_Op);dist_ctd=fscanf(fid_Op,'%f');
fclose(fid_Op);

if(ctds_inf>0);
    n_ctd=ctds_inf;ctds_inf=1;
    ctd_pos=1;dist_ctd=dist_ctd';
end

clear fid_Op naminfo

%%%%%% plot type %%%%%%
 fprintf('\nDefine the contours of plot?  (isolines  => 1)  ');
 fprintf('\n                   (pos and neg isolines => 2)  ');
 style=input('\n                                  (color => 3) : ');

 
 opp=input('\nPlot measured velocity field? Yes => 1 :  ');

 
%%%%%%%%%%%%%%%%%%%
% read input file %
%%%%%%%%%%%%%%%%%%%

files=dir([files_header,'*']);
n_perf=length(files);
if(ini_margin==1),inverte=zeros(1,n_perf);else inverte=ones(1,n_perf);end
for fl=1:n_perf,file_names{fl}=files(fl).name;end

clear files fl files_header


% Defining the size of the matrix for each ADCP section

if(length(dist_inicio)==1),dist_inicio(1:n_perf)=dist_inicio;end

for k=1:n_perf

    fprintf('\nfile: %d  \n',k);
    file = [dir_path,file_names{k}];
    T = dlmread(file);T = T';%T = load(file);
    
    
    time=T(1:7,:);
    dx=T(9,:);
    dy=T(8,:);   
    m=(size(T,1)-10)/4;
    
    %Calculating DT for each section
    
    hour(k,1)= time(4,1) + time(5,1)/60;
    hour_i(k,1)=hour(k,1);
    dia = time(3,1);ano=time(1,1);mes=time(2,1);
        
    if(k==1)
        dia_ini = dia;
        data{1}= sprintf('%02g/%02g/%02g',mes,dia,ano);
    elseif(dia>dia_ini)
        data{2}= sprintf('%02g/%02g/%02g',mes,dia,ano);
    end             
    
    if(k>1 && hour(k,1)<hour(k-1,1))
        hour_i(k,1)=hour(k-1,1)+hour(k,1);hour_v = hour(k-1,1);
    elseif(k>1 && hour(k,1)<hour_i(k-1,1)),hour_i(k,1)=hour_i(k-1,1)+ ...
            (hour(k,1)-hour(k-1,1));
    end
    
    DT_sec(k)= (time(3,end)-time(3,1))*24 + [(time(4,end) + time(5,end)/60)...
        - (time(4,1) + time(5,1)/60)];
    hour_end(k)= hour_i(k)+ DT_sec(k);
    
    %Calculating the river angle based on the section angle
    a=mean(diff(dy));b=mean(diff(dx)); ang_sec = (atand(abs(b/a)));

    if(a > 0 && b > 0),ang_sec=ang_sec;end
    if(a < 0 && b > 0),ang_sec=ang_sec + 90;end
    if(a < 0 && b < 0),ang_sec=ang_sec + 180;end
    if(a > 0 && b < 0),ang_sec=ang_sec + 270;end
              
    if(ini_margin==1)
        ang_rio(k)=ang_sec - 90;
    else
        ang_rio(k)=ang_sec + 90;
        if(ang_rio(k)>360),ang_rio(k)=ang_rio(k)-360;end
    end
    %Defining which sections must be inverted
    
    if(k>1 && abs(ang_rio(k)-ang_rio(1))>140)
        if(ang_rio(k)<180)
            ang_rio(k)=ang_rio(k)+180;
        else
            ang_rio(k)=ang_rio(k)-180;
        end
        if(ini_margin==1),inverte(k)=1;else inverte(k)=0;end
    end                     
    
    % sum distance to the bouy in the cross profile beginning plus ensemble
    % avg distance
    dist_tmp=sqrt(dx.^2+dy.^2) + dist_inicio(k);
    dist = dist_tmp + mean(diff(dist_tmp));
    
    u=T(11+2*m:11+3*m-1,:);
    v=T(11+m:11+2*m-1,:);
    [m,n]=size(v);Cont=find(v==flag);
    u(Cont)=nan;v(Cont)=nan;
    
    [v_cross,v_along]=rot(u,v,ang_rio(k));
    
    depth=-1.*T(11:11+m-1,1);
    
    depth_int= mean(diff(depth));
    speed=v_along; % along channel velocity
    ncol1=length(dist);
    
    depth_max=-1*T(10,:)';
    
    clear T m n
        
 % estimating the total distance, max depth, Umax and Umin for each cross profile
   
    max_depth_time(k)=min(depth_max);
    max_dist_time(k)=max(dist);
    min_dist_time(k)=min(dist);
    max_speed_time(k)=max(max(speed));
    min_speed_time(k)=min(min(speed));
    length_time(k)= length(dist);
    length_time_m(k)= nansum(diff(dist));
    ens_size(k)= mean(diff(dist));
    
    clear ind0 ind1 depth_vector speed;
%   
    clear adcp dx dy dist dist_tmp u v v_cross v_along file
end
    clear ncol1 time ang_sec a b Cont ano mes dia dia_ini k

if(ang_manual==9999),ang_rio = mean(ang_rio);
else ang_rio=ang_manual;end


%Calculating the mean DT for all sections
DT_sec = mean(DT_sec);

%Defining which sections mean the same Tidal time step
diff_hour=ones(n_perf-1,1);
for k=2:n_perf
    if((hour_i(k)-hour_end(k-1))<DT_sec)
        diff_hour(k-1)=0;
    end
end

hour=hour_i;

clear hour_i ang_manual
% estimate total distance, max depth, Umax and Umin for the entire tidal
% cycle
    max_dist=max(max_dist_time);
    min_dist=min(min_dist_time);
    max_depth=floor(min(max_depth_time)); 
    max_speed=round(max(max_speed_time)*10)/10;
    min_speed=floor(min(min_speed_time)*10)/10;
    max_length= max(length_time);
    min_length= min(length_time);
    avg_length= mean(length_time);
    avg_length_m= mean(length_time_m);
    avg_ens_size = mean(ens_size);
    CAX = max([abs(min_speed) abs(max_speed)]);
    
    clear max_dist_time max_depth_time ens_size min_dist_time ...
        length_time length_time_m min_speed_time max_speed_time

%%%%%%%%%%%%%%%%
% USER OPTIONS %
%%%%%%%%%%%%%%%%

fprintf('\n ************************************************************************** \n\n\n');
fprintf('Maximum cross profile distance [m] =  %5.1f  \n',max_dist);
fprintf('Minimum number of ensembles =  %4.1f  \n',min_length);
fprintf('Maximum number of ensembles =  %4.1f  \n',max_length);
fprintf('Average ensembles length =  %4.1f  \n',avg_ens_size);
disp(' ')
fprintf('Average length of sections [m] =  %4.1f  \n',avg_length_m);
fprintf('Maximum depth [m] =  %4.2f  \n',max_depth);
fprintf('Depth for the first data cell [m] =  %4.2f  \n',depth(1));
fprintf('ADCP cell heigth [m] =  %4.2f  \n',-depth_int);
disp(' ')
disp(' ---------  STANDARD SECTION SETTINGS  -------------- ')
disp(' ')

n_vert_m = input(' - Set the horizontal resolution (eg. average ensemble length) [m] >>  ');
c_vert = input(' - Set the vertical resolution [m] >>  ');
n_vert = round(avg_length_m/n_vert_m);


clear avg_length avg_length_m avg_ens_size n_vert_m

% SET THE MAIN MATRIXES DIMENSIONS
    
    dist_i=linspace(min_dist,max_dist,n_vert);   % distances and number of columns
    
    depth_i=-c_vert/2:-c_vert:max_depth;         % matrix line number

    z=0:0.1:1;                                   % adimensional Z depths
    masknan=ones(length(depth_i),length(dist_i),n_perf);
    mask=ones(length(depth_i),length(dist_i),n_perf);
    speed_i=ones(length(depth_i),length(dist_i),n_perf);
    V_z=ones(length(z),length(dist_i),n_perf);

    clear min_dist max_dist max_depth min_length max_length;
    
% read ADCP data for the time k

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1st K looping
for k=1:n_perf

    fprintf('\n ******************************************************************** \n\n');    
        file = [dir_path,file_names{k}];
    T = dlmread(file);T = T';
    
    dx=T(9,:);
    dy=T(8,:);
    m=(size(T,1)-10)/4;
    
    % sum distance to the bouy in the cross profile beginning plus ensemble
    % avg distance
    dist_tmp=sqrt(dx.^2+dy.^2) + dist_inicio(k);
    dist = dist_tmp + mean(diff(dist_tmp));
    
    u=T(11+2*m:11+3*m-1,:);
    v=T(11+m:11+2*m-1,:);
    [m,n]=size(v);Cont=find(v==flag);
    u(Cont)=nan;v(Cont)=nan;
    
    [v_cross,v_along]=rot(u,v,ang_rio);

    depth=-1.*T(11:11+m-1,1);
    
    depth_int= mean(diff(depth));
    speed=v_along; % along channel velocity
    
    ncol1=length(dist);      
    nline1=length(depth);
    
    depth_max=-1*T(10,:)';
    
    clear T m n
        
    % check if there is ensemble without bottom track and use interpolation to
    % fill the gap
    ind0=isnan(depth_max);
    ind1=find(ind0==1);
    ind2=find(ind0==0);
    if ind1 > 0
    depth_max(ind1)=interp1(dist(ind2),depth_max(ind2),dist(ind1),'linear','extrap');
    end
    clear file adcp dx dy u v v_cross v_along ind0 ind1 ind2
   
    % check if it is necessary to reverse the cross profile direction

    if inverte(k)==1
        speed_tmp(:,1:ncol1)=speed(:,ncol1:-1:1);
        speed=speed_tmp;
        depth_max_tmp(1:ncol1)=depth_max(ncol1:-1:1);
        depth_max= depth_max_tmp;
        for j=1:ncol1
         dist_tmp_inv(j)=mean(diff(dist_tmp))+dist_tmp(ncol1)...
         -dist_tmp(ncol1-j+1);
        end
        dist=dist_tmp_inv;
        clear dist_tmp_inv speed_tmp dist_tmp depth_max_tmp
    end
    
    % matrix redimension
    for j=1:ncol1
        ind_tmp=isnan(speed(:,j));
        ind0=find(ind_tmp==0);
        if length(ind0)>=1 % check if column has data 
            ind1(j)=1;
        else
            ind1(j)=0;
        end
    end
    ind2= find(ind1==1); % find god data columns
    ncol2= length(ind2); % how many
    speed_z=zeros(length(depth_i),ncol2);
    clear ind2;
    
    % find bad ensemble columns
    n_bad= find(ind1==0);
    if length(n_bad) < 1; n_bad=0; end
    bad_seq= diff(n_bad);
    c=0;
    for kk=1:length(bad_seq)
        if bad_seq(kk) == 1
        c=c+1;
        end
    end
    
    
    fprintf('\nCross profile n. : %d   \n',k);
    fprintf('Bad ensemble id :  %d   \n',n_bad);
    fprintf('N. of sequencial bad ensembles : %d \n\n',c);
    clear ind0 ind1 n_bad bad_seq kk c

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % set matrix speed to the god column data in speed_z   %
    % vertical INTERPOLATION - top and bottom              %
    % Eliminate cells in the Side Lobe Layer Pulse or      %
    % one cell close to bottom                             %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
   Dtop=-depth(1); %depth to the middle of the first cell
   Dsurface= Dtop - (-depth_int)/2; % depth of de Top Discharge Layer

   counter=0; 
   for j=1:ncol1
        ind_tmp=isnan(speed(:,j));
        ind0=find(ind_tmp==0);
        ind1=find(ind_tmp==1);
        speed(ind1,j)=0; % set NaN to ZERO
        if length(ind0)>=1
           counter=counter+1;
           maxdepth_col(counter)=depth_max(j);
           dist_col(counter)=dist(j);
           Dlg_max=-depth(ind0(end)); % WinRiver II ASCII files delete cells after Dlg_max
                                        % jul/2010

           % set bottom layer cells
           ind_pos=find(-depth > Dlg_max & -depth <= -maxdepth_col(counter)); 
        
           if length(ind_pos) >= 1 && ind_pos(1) > 1
               n_pos=length(ind_pos);
               depth_alt(1:ind_pos(1)-1)=depth(1:ind_pos(1)-1); % set cells before Dlg
               speed_alt(1:ind_pos(1)-1,counter)=speed(1:ind_pos(1)-1,j);
               depth_alt(ind_pos(1))=maxdepth_col(counter); % cell after Dlg => bottom depth
               speed_alt(ind_pos(1),counter)=0; % speed = 0 at bottom;
               nn1=length(depth)-ind_pos(end);
                for i=1:nn1 % cell after bottom, modified depth´s
                    depth_alt(ind_pos(1)+i)=depth(ind_pos(end)+i);
                    speed_alt(ind_pos(1)+i,counter)=0; % speed = 0 after bottom;
                end
                
                for i=1:n_pos-1 % set last cells depth and speed (null)
                     depth_alt(nline1-n_pos+1+i)=depth(end)+i*depth_int;
                     speed_alt(nline1-n_pos+1+i,counter)=0;
                end
               
           else
               depth_alt=depth;
               speed_alt(:,counter)=speed(:,j);
           end
           
           % new depth interpolation and to the top extrapolation
           speed_z(:,counter)=interp1(depth_alt,speed_alt(:,counter),depth_i','linear');
           speed_top(:,counter)=interp1(depth_alt,speed_alt(:,counter),depth_i','nearest','extrap');
           dif_depth_top=-depth(1)/c_vert;
           ncell_top=floor(dif_depth_top);
           for dd=1:ncell_top+1
               speed_z(dd,counter)=speed_top(dd,counter); % nearest extrap to top
           end
           
            % mark bottom cells in the new matrix speed_z
           ind_pos_i=find(-depth_i > Dlg_max & -depth_i <= -maxdepth_col(counter)); % set bottom layer cells
           
           if ind_pos_i >= 1 
                Dbottom(counter)=-(maxdepth_col(counter)-depth_i(ind_pos_i(1))) + c_vert/2;
                speed_pos=mean(speed_z(ind_pos_i,counter));
                if ind_pos_i(1) > 1
                    speed_bott(counter)=(speed_z(ind_pos_i(1)-1,counter)+speed_pos+0)/3;  % mean speed - bottom layer
                else
                    speed_bott(counter)=speed_pos/2;  % mean speed - bottom layer
                end
           else
               Dbottom(counter)=0; speed_bott(counter)=0;
           end
           
            % set middle layer cells
           ind_mid=find(-depth_i >= Dsurface & -depth_i <= Dlg_max); % set middle layer cells
           if ind_mid > 0
           speed_mid(counter)=mean(speed_z(ind_mid,counter));       % mean speed - middle layer
           Dmiddle(counter)= Dlg_max - Dsurface; % altered in jul/2010 (old was: length(ind_mid)*c_vert;) 
           else
           speed_mid(counter)=0;
           Dmiddle(counter)= 0;
           end
           
           
            % calc discharge bottom & middle layers
            if counter>1
                Abottom(counter)=Dbottom(counter)*(dist_col(counter)-dist_col(counter-1));
                Qbottom(counter)=Abottom(counter)*speed_bott(counter);
                Amiddle(counter)=Dmiddle(counter)*(dist_col(counter)-dist_col(counter-1));
                Qmiddle(counter)=Amiddle(counter)*speed_mid(counter);
            else
                Abottom(counter)=Dbottom(counter)*dist_col(counter);
                Qbottom(counter)=Abottom(counter)*speed_bott(counter);
                Amiddle(counter)=Dmiddle(counter)*dist_col(counter);
                Qmiddle(counter)=Amiddle(counter)*speed_mid(counter);
            end

%            figure(j) %vertical profile
%            plot(speed(1:ind_pos(end),counter),depth(1:ind_pos(end)),'s',...
%                speed_z(1:ind_pos_i(end),counter),depth_i(1:ind_pos_i(end)),':')
% 
%            pause(4)
%            close

           clear ind_pos depth_alt speed_alt ind_mid ind_pos_i n_pos dd
        end
        clear ind_tmp ind0 ind1 dif_depth_top ncell_top Dlg_max Dsl...
            nn1 speed_pos
   end
   
   ind_tmp=isnan(Qbottom);
   ind_0=find(ind_tmp==0);
   QTOT_bottom=sum(Qbottom(ind_0));
   ATOT_bottom=sum(Abottom(ind_0));
   clear ind_tmp ind_0 Abottom Qbottom;
   
   ind_tmp=isnan(Qmiddle);
   ind_0=find(ind_tmp==0);
   QTOT_middle=sum(Qmiddle(ind_0));
   ATOT_middle=sum(Amiddle(ind_0));
   clear ind_tmp ind_0 Amiddle Qmiddle;
   
   Dtop=-depth(1); %depth to the middle of the first cell
   Dsurface= Dtop - (-depth_int)/2; % depth of de Top Discharge Layer
   A_surface=dist_col(end)*Dsurface;
      
   ind_sf=find(-depth_i < Dtop); % set surface layer cells
   if length(ind_sf) > 1
       speed_sf=mean(mean(speed_z(ind_sf,:))); % mean speed at surface layer
   else
       speed_sf=mean(speed_z(1,:));
   end
   
   QTOT_surface=speed_sf*A_surface;
   
   QTOT=QTOT_middle+QTOT_surface+QTOT_bottom;
   ATOT=ATOT_middle+A_surface+ATOT_bottom;
   Uavg=QTOT/ATOT;
   
    Q_BOT(k)=QTOT_bottom;
    Q_MID(k)=QTOT_middle;
    Q_SUR(k)=QTOT_surface;
    Q_TOTAL(k)=QTOT;

    A_BOT(k)=ATOT_bottom;
    A_MID(k)=ATOT_middle;
    A_SUR(k)=A_surface;
    A_TOTAL(k)=ATOT;
    U_m(k)=Uavg;
    
    %   pause(10)
   
   clear counter maxdepth_tmp depth Dlg_max speed_top Dbottom Dmiddle...
        speed_mid speed_bott ind_sf speed_sf depth_max dist QTOT_bottom...
        QTOT_middle QTOT_surface QTOT ATOT_bottom ATOT_middle A_surface...
         ATOT Uavg;
   

    %%%%%%%%%%%%%%%%%%%%%%%%
    % interpolate the grid %
    %%%%%%%%%%%%%%%%%%%%%%%%

    speed_i(:,:,k) = interp2(dist_col,depth_i',speed_z,dist_i,depth_i','linear',0);
    nline_i= length(depth_i);
    ncol_i=length(speed_i(1,:,k));
    clear speed_z speed

    % max depth interpolation in the new matrix
    maxdepth_i(k,:)=interp1(dist_col,maxdepth_col,dist_i,'linear','extrap');

    izero=find(maxdepth_i(k,:)>0); % check extrapolation
    if izero > 0; maxdepth_i(k,izero)=0; end

    clear maxdepth_col izero dist_col

    % calculate velocity in Z adimensional depths

    for j=1:ncol_i
        ind0=find(depth_i >= maxdepth_i(k,j));
        if sum(ind0 > 1)
        z_tmp=depth_i(ind0)/maxdepth_i(k,j);
        V_z(:,j,k)= interp1(z_tmp',speed_i(ind0,j,k),z','linear','extrap');
        clear ind0 z_tmp;
        else
        V_z(:,j,k)= NaN;
        end
    end

    % define the cross profile mask

    for j=1:ncol_i
    ind0=find(depth_i < maxdepth_i(k,j)); % all cells below max depth are set to NaN
    masknan(ind0,j,k)=NaN;
    clear ind0
    end

    mask(:,:,k)=speed_i(:,:,k).*masknan(:,:,k);


    % plot cross section velocity field

    if opp == 1
        figure(k);
        contourf(dist_i,depth_i,mask(:,:,k),'k');
        caxis([min_speed max_speed]);
        colorbar;
        hold on;
        plot(dist_i,maxdepth_i(k,:),'-k','linewidth',2)
        hold on;

        if ctds_inf == 1
            if ctd_pos == 1
                for i=1:n_ctd
                    plot(dist_ctd(i),0,'k.','MarkerSize',24);
                    if i==1
                    text(dist_ctd(1),0.2,'  1','fontsize',8);
                    else
                    text(dist_ctd(2),0.2,'  2','fontsize',8);
                    end
                end
            end
        end
        title([label,' - velocity field']);
        xlabel('Distance (m)');
        ylabel('Depth (m)');

        pause(2)
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%
    % calculate discharge %
    %%%%%%%%%%%%%%%%%%%%%%%
    
    ind0=isnan(mask(1,:,k)); 
    ind1=find(ind0==0);
    delta_dist=dist_i(2)-dist_i(1);
    area_i=delta_dist*c_vert;
    Q(:,ind1)=area_i*mask(:,ind1,k);
    ind2=isnan(Q);
    ind3=find(ind2==0);
    discharge(k)=sum(Q(ind3));
    area_T(k)=sum(area_i*length(ind3));
    
    clear ind0 ind1 ind2 ind3 area_i Q;  

    % calculate velocities

    mean_speed(k)=discharge(k)/area_T(k); % general average
    mean_col_speed_i(k,1:ncol_i)=smartmean(mask(:,1:ncol_i,k),1); % U column
    vert_speed_profile(1:nline_i,k)=smartmean(mask(1:nline_i,:,k),2); % mean vertical profile
    mean_surface_speed(k)=vert_speed_profile(1,k);% vel close to the surface
    
end %%%%%%%%%%%%%%%%% 1st K Looping %%%%%%%%%%%%%%%%%%%%%

mpp = find(hour>24);
if isempty(mpp) == 0
    hour2=hour;hour2(mpp)=hour2(mpp)-hour_v;
else
    hour2=hour;
end

clear min_speed max_speed ncol1 nline1;

% save 
fout_name=[dir_path,'disch_layers_res.txt'];
ff=fopen(fout_name,'w');
fprintf(ff,'Discharge calculated by ProcED before space and time interpolations\n');
var=[hour2';Q_TOTAL;Q_SUR;Q_MID;Q_BOT;U_m];
fprintf(ff,'%s\n','Time   Qtotal QTop   Qmiddle Qbottom  Uavg');
fprintf(ff,'%2.1f %8.1f %8.1f %8.1f %8.1f %8.3f\n',var);

var1=[hour2';A_TOTAL;A_SUR;A_MID;A_BOT];
fprintf(ff,'%s\n','Time   Atotal ATop   Amiddle Abottom');
fprintf(ff,'%2.1f %8.1f %8.1f %8.1f %8.1f \n',var1);
fclose(ff);
clear fout_name var var1 Cont;

save disch_layers hour2 Q_TOTAL Q_SUR Q_MID Q_BOT U_m A_TOTAL A_SUR A_MID A_BOT

pause(5);
close all;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET AVERAGED CROSS PROFILE VEL FOR THE SAME TIME, DISCHARGE AND VEL     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%diff_hour=diff(hour);
ind_hour=find(diff_hour ~= 0);
last_hour=max(ind_hour)+1;
if ind_hour >=1
    hour_i= hour(ind_hour);
    hour_i(length(ind_hour)+1)=hour(last_hour);
else
    hour_i=hour(1);
end

n_time=length(hour_i); 
n_z=length(z);

% Calculate averaged cross profile vel for the same time (if exists)
    
for k=1:n_time %%%%%%%%%%%%%%%%%%%%% 2nd K Looping %%%%%%%%%%%%%%%%%%%
    
    ind1=find(hour==hour_i(k));

    if length(ind1)>1
    speed_profile(:,:,k)=smartmean(mask(:,:,ind1),3);
    V_z_hour(:,:,k)=smartmean(V_z(:,:,ind1),3);
    maxdepth_hour(k,:)=smartmean(maxdepth_i(ind1,:),1);
    T_discharge(k)=smartmean(discharge(ind1),2);
    T_area(k)=smartmean(area_T(ind1),2);
    else
    speed_profile(:,:,k)=mask(:,:,ind1);
    V_z_hour(:,:,k)=V_z(:,:,ind1);
    maxdepth_hour(k,:)=maxdepth_i(ind1,:);
    T_discharge(k)=discharge(ind1);
    T_area(k)=area_T(ind1);
    end
    T_mean_speed(k)=T_discharge(k)/T_area(k); % average vel
    
end %%%%%%%%%%%%%%%%%%%%% 2nd K Looping %%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE TIDAL CYCLE END AND SAME TIME STEP VEL FIELDS %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if n_time > 1
    
    % define the time when the inicial discharge is reached at end of tidal cycle
    
    time_end=interp1(T_discharge(end-2:end),hour_i(end-2:end),T_discharge(1),'linear','extrap');
    dt=(time_end-hour_i(1))/24; % same time step
    new_hour=[hour_i(1):dt:time_end]; % new time array

    new_n_time=length(new_hour); %%%%%%%
    
    % interpolate U in Zadm and maxdepth to the new time
    % calculate residual cross section field
    clear ind1 ind2;
    for j=1:ncol_i
        V_tmp(1:n_z,1:n_time)= reshape(V_z_hour(:,j,:),[length(V_z_hour(:,1,1)) length(V_z_hour(1,1,:))]);
        new_maxdepth(1:new_n_time,j) = interp1(hour_i,maxdepth_hour(:,j),new_hour','linear','extrap');
        ind1=isnan(V_tmp(1,:));
        ind2=find(ind1==0);
        for i=1:n_z
        new_V_z(i,j,1:new_n_time) = interp1(hour_i(ind2),V_tmp(i,ind2),new_hour','linear','extrap');
        net_speed(i,j)= (0.5*new_V_z(i,j,1)+ sum(new_V_z(i,j,2:end-1))+ 0.5*new_V_z(i,j,end))/24; % residual field
        end
    end

 
    % set bottom cell vel to = 0 in Z adm profiles

    ind0=isnan(new_V_z(11,:,:));
    ind1=find(ind0==0);
    new_V_z(11,ind1)=0;

    % smoothing velocity cross profile fields
    
    fprintf('\n\n # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #\n\n');    
    
    plot_filter = input('Would you like to smooth velocity field for ploting? Yes =>1 :  ');
    
    if plot_filter==1
        nf = input('\n\n Smoothing - How many times?  :  ');    
        fprintf('\n\n ...smoothing');
    end

    % calculates discharge, U, etc
    
    for k2=1:new_n_time  %%%%%% K2 Looping %%%%%%
               
        izero=find(new_maxdepth(k2,:)>0); % check bad bottom depth extrapolation
        if izero > 0; new_maxdepth(k2,izero)=0; end

        ind0=isnan(new_V_z(1,:,k2)); 
        ind1=find(ind0==0);
        delta_dist=dist_i(2)-dist_i(1);
        area_i(:,k2)=delta_dist*-new_maxdepth(k2,:)';
        QT(ind1)=area_i(ind1,k2).*((0.5*new_V_z(1,ind1,k2)+ sum(new_V_z(2:end,ind1,k2),1))/10)';
        ind2=isnan(QT);
        ind3=find(ind2==0);
        new_discharge(k2)=sum(QT(ind3));
        new_area(k2)=sum(area_i(ind3,k2));
        new_mean_speed(k2)=new_discharge(k2)/new_area(k2); %
        new_mean_col_speed_i(k2,1:ncol_i)=smartmean(new_V_z(:,1:ncol_i,k2),1);
        new_vert_speed_profile(1:n_z,k2)=smartmean(new_V_z(1:n_z,:,k2),2); 
        new_mean_surface_speed(k2)= new_vert_speed_profile(1,k2);

        clear ind1 ind2 ind3; 

        %%%%%%%%%%%%% plot velocities fields %%%%%%%%%%%%

        % 5 horizontal cells moving average filter

        fprintf('.');

        if plot_filter==1
            for i=1:nf
                m_size=2;
                new_m=filter_5(new_V_z(:,:,k2),n_z,ncol_i,m_size);
                new_V_z(:,:,k2)=new_m;
                clear new_m m_size
            end
        end

        % reverse depth Z adm to absolute depths

        for j=1:ncol_i
        depth_ret(:,j,k2)=z'.*new_maxdepth(k2,j);
        end
        for i=1:n_z
        dist_ret(i,:)=dist_i;
        end

        % plot        
        isoline(style,new_V_z(:,:,k2),dist_ret,depth_ret(:,:,k2),dist_i,new_maxdepth(k2,:),k2)
        if(style == 3)
            caxis([-1*CAX CAX]);
            colb=colorbar;
            title(colb,'m/s');
        end
       
        if(new_hour(k2)>24)
            hour_plot = num2str(new_hour(k2)-hour_v);
            new_hour(k2)=new_hour(k2)-hour_v;
            hour_plot=[data{2}, ' - ', hour_plot(1:4)];
        else
            hour_plot = num2str(new_hour(k2));hour_plot=hour_plot(1:4);
            hour_plot=[data{1}, ' - ', hour_plot(1:4)];
        end
        
        datetime{k2} = hour_plot;
        
        title([label,' - velocity field - T',num2str(k2),' (',hour_plot,' Hrs)']);
        if ctds_inf == 1
            if ctd_pos == 1
                for i=1:n_ctd
                    plot(dist_ctd(i),0,'k.','MarkerSize',24);
                    if i==1
                    text(dist_ctd(1),0.2,'  1','fontsize',8);
                    else
                    text(dist_ctd(2),0.2,'  2','fontsize',8);
                    end
                end
            end
        end
        hgsave(['vel_field_T',num2str(k2),'.fig']);

        % Calculate data to the CTD position

        if ctds_inf == 1
            depth2ctd(k2,:) = interp1(dist_i,new_maxdepth(k2,:),dist_ctd,'linear','extrap');
            vel2ctd(:,:,k2) = interp2(dist_i,z',new_V_z(:,:,k2),dist_ctd,z','linear');
        end
    
    end %%%%%%%%%%% K2 Looping %%%%%%%%%%%%

    clear speed_temp speed_tmp2 m_tmp_1 m_tmp_2 diff_hour ind_hour last_hour;

    % calculate residual vel profile to the CTD position
    
    if ctds_inf == 1
    net_speed_profile=interp2(dist_i,z',net_speed,dist_ctd,z','linear');
    net_speed_profile(11)=0; % bottom speed == zero
    end
    
    % reverse depth Z adm to absolute depths
    
    for j=1:ncol_i
    depth_ret2(:,j)=z'.*smartmean(maxdepth_hour(:,j),1);
    end
    for i=1:length(z)
    dist_ret2(i,:)=dist_i;
    end
    
    %%%% plot residual vel %%%%

        % 5 horizontal cells moving average filter

    if plot_filter==1
        for i=1:nf
            S=whos('net_speed');
            m_size=length(S.size);
            new_m=filter_5(net_speed,n_z,ncol_i,m_size);
            net_speed=new_m;
            clear new_m m_size S;
        end
    end
    
    % plot
        
    k=40;
    isoline(style,net_speed,dist_ret2,depth_ret2,dist_ret2(11,:),depth_ret2(11,:),k)
    if(style == 3)
            caxis([-1*CAX CAX]);
            colb=colorbar;
            title(colb,'m/s');
    end     

    title([label,' - residual field']);
    if ctds_inf == 1
        if ctd_pos == 1
            for i=1:n_ctd
                plot(dist_ctd(i),0,'k.','MarkerSize',24);
                if i==1
                text(dist_ctd(1),0.2,'  1','fontsize',8);
                else
                text(dist_ctd(2),0.2,'  2','fontsize',8);
                end
            end
        end
    end
    hgsave('residual_field.fig');

    clear ncol_i nline_i;

    % calculate residual discharge

    ind0=isnan(net_speed(1,:));
    ind1=find(ind0==0);
    delta_dist=dist_i(2)-dist_i(1);
    area_i(ind1)= delta_dist*-(smartmean(new_maxdepth(:,ind1),1));
    Q(ind1)= area_i(ind1).*mean(net_speed(:,ind1));

    net_discharge= sum(Q(ind1));      
    net_area= sum(area_i(ind1));
    net_speed_g = net_discharge/net_area; %

    clear ind0 ind1 area_i Q;  

end % %%%%%%%%%%%% IF end %%%%%%%%%%%%%%%

mpp = find(hour_i>24);
if isempty(mpp) == 0
    hour_i(mpp)=hour_i(mpp)-hour_v;
end

%%%%%%%%%%%%%%%%%%%%%
% save output files %
%%%%%%%%%%%%%%%%%%%%%
    
if n_time > 1
    
    % residual field

    fout_name=[dir_path,'residual_field.txt'];
    ff=fopen(fout_name,'w');
    var=[net_discharge;net_speed_g];
    fprintf(ff,'%s\n',' Residual Disch(m3/s) Residual Vel(m/s)');
    fprintf(ff,'%9.2f       %9.3f\n',var);
    fclose(ff);
    clear fout_name var;

    % only 13 profiles to adcp2ctd.m and Estuario.m

    if ctds_inf == 1    
    ind_ctd=1:2:new_n_time; % when new_n_time = 25 ######
    maxdepth_ctd=depth2ctd(ind_ctd,:);
    speed_ctd= vel2ctd(:,:,ind_ctd);
    hour_ctd=new_hour(ind_ctd);
    end
    
    % resume

    fout_name=[dir_path,'run_sets.txt'];
    ff=fopen(fout_name,'w');
    fprintf(ff,'Initialization parameters \n');
    fprintf(ff,'\nProject Title : %s\n',label);
    fprintf(ff,'\nInitial date and time : %s\n',datetime{1});
    fprintf(ff,'\nEnding date and time : %s\n',datetime{end});
    fprintf(ff,'\nPath : %s\n',dir_path);
    fprintf(ff,'\nChannel angle (grads)  : %6.1f\n',ang_rio);
    fprintf(ff,'\nNumber of verticals    : %6.1f\n',n_vert);
    fprintf(ff,'\nVertical resolution [m]: %6.1f\n',c_vert);
    if ctds_inf == 1
        for i=1:n_ctd
        fprintf(ff,'\nCTD distance [m]     : %8.1f\n',dist_ctd(i));
        end
    end
    fprintf(ff,'\nSmooth velocity field for ploting? Yes =>1: %6.0f\n',plot_filter);
    if plot_filter == 1
        fprintf(ff,'\nSmoothing - How many times? %6.0f\n',nf);
    end
    fclose(ff);
    clear fout_name;

    fout_name=[dir_path,'disch_res_interp.txt'];
    ff=fopen(fout_name,'w');
    fprintf(ff,'Discharge calculated by ProcED after space and time interpolation\n');
    var=[new_hour;new_mean_speed;new_area;new_discharge];
    fprintf(ff,'%s\n','Time   U(m/s) Area(m2) Disch(m3/s)');
    fprintf(ff,'%2.1f %6.3f %8.1f %8.1f \n',var);
    fclose(ff);
    clear fout_name var;

    fout_name=[dir_path,'dacv_res.txt'];
    ff=fopen(fout_name,'w');
    fprintf(ff,'Depth averaged velocity calculated by ProcED after space and time interpolation\n');
    fprintf(ff,'1st line is time step\n');
    var1=[new_mean_col_speed_i];
    fprintf(ff,'%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f \n',new_hour);
    fprintf(ff,'%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f \n',var1);
    fclose(ff);
    clear fout_name var1;

    % save matrixs
    save tidal_cycle_matrix depth_i dist_i speed_profile new_V_z depth_ret dist_ret net_speed depth_ret2 dist_ret2 
    save disch_res new_hour new_mean_speed new_area new_discharge data
    if ctds_inf == 1
        ind_ctd=[1:2:new_n_time]; % only 13 profiles to adcp2ctd.m
        save inf2ctd dist_ctd hour_ctd z speed_ctd maxdepth_ctd
    end
else
    save tidal_cycle_matrix depth_i dist_i speed_profile % save matrixs
    save disch_res hour_i T_mean_speed T_area T_discharge data
    if ctds_inf == 1
        save inf2ctd dist_ctd hour_i z speed_ctd maxdepth_ctd
    end

end

% save resume
fout_name=[dir_path,'disch_res2adcp.txt']; %discharge calculated by ProcED before space and time interpolations
ff=fopen(fout_name,'w');
fprintf(ff,'Discharge calculated by ProcED before space and time interpolations\n');
var=[hour2';U_m;A_TOTAL;Q_TOTAL];
fprintf(ff,'%s\n','Time   U(m/s) Area(m2) Disch(m3/s)');
fprintf(ff,'%2.1f %6.3f %8.1f %8.1f \n',var);
fclose(ff);

% save resume
fout_name=[dir_path,'disch_res1.txt']; %discharge calculated by ProcED after cross-section interpolation for each tidal stage
ff=fopen(fout_name,'w');
fprintf(ff,'Discharge calculated by ProcED after cross-section interpolation\n');
fprintf(ff,'(mean value if there were 2 or more measurements to represent one tidal stage)\n');
var=[hour_i';T_mean_speed;T_area;T_discharge];
fprintf(ff,'%s\n','Time   U(m/s) Area(m2) Disch(m3/s)');
fprintf(ff,'%2.1f %6.3f %8.1f %8.1f \n',var);
fclose(ff);

%clear all
