function [tef] = calcTEF(tef,grid,N,Odir)

%inputs:
%1) the number of salinity classes (N)
%2) the Grid output structure from MITgcm mnc package
%3) an output directory to store files (Odir)
%4) and a "tef" structure with the following fields for each TEF transect:
%   Z = coordinates of domain cells (vector)
%   time = model time in seconds (vector)
%   S = salinity of each cell in transect (X by Z by time matrix)
%   HfacC = HFacC field from the Grid MITgcm output (X by Z matrix)
%   eta = surface displacement (X by time matrix)
%   U = u velocity for transect (X by Z by time matrix)
%   V = v velocity for transect (X by Z by time matrix)
%   velDir = indicate which direction is landward (options are: +u,+v,-u,-v)
%   axisDir = which axis transect falls along (options are x or y)

    for i = 1:length(tef)
        %% create s_ik and u_ik matrices
        [nx,nz,nt] = size(tef(i).S);
        s = reshape(tef(i).S,[nx*nz,nt]);
        if tef(i).velDir == '+u' %inward normal velocity (positive landwards)
            u = reshape(tef(i).U,[nx*nz,nt]);
        elseif tef(i).velDir == '+v'
            u = reshape(tef(i).V,[nx*nz,nt]);
        elseif tef(i).velDir == '-u'
            u = reshape(-tef(i).U,[nx*nz,nt]);
         elseif tef(i).velDir == '-v'
            u = reshape(-tef(i).V,[nx*nz,nt]);
        elseif tef(i).velDir == '+r' % use both u and v to project onto a rotated line (defined as Xperp and Yperp)
            for ii = 1:nt
                for iii = 1:nz
                    [perp(:,iii,ii),~] = UPerpPara(tef(i).Xperp,tef(i).Yperp,tef(i).U(:,iii,ii),tef(i).V(:,iii,ii));
                end
            end
            u = reshape(perp,[nx*nz,nt]);
        elseif tef(i).velDir == '-r'
            for ii = 1:nt
                for iii = 1:nz
                    [perp(:,iii,ii),~] = UPerpPara(tef(i).Xperp,tef(i).Yperp,tef(i).U(:,iii,ii),tef(i).V(:,iii,ii));
                end
            end
            u = reshape(-perp,[nx*nz,nt]);
        else
            error('indicate velDir ? u or v direction for velocity transect')
        end

        %% calc cell areas
        dz = abs(diff(grid.Zp1)); %use vertical cell coordinate at cell interface
        DZ = repmat(dz,[1,nx,nt]);
        DZ = permute(DZ,[2,1,3]);
        aa = reshape(tef(i).eta,[nx,1,nt]);
        bb = DZ(:,1,:) + aa;
        DZ(:,1,:) = bb;

        HFacC = repmat(tef(i).HFacC,[1 1 nt]); %partial z cells
        DZ = DZ.*HFacC;

        if tef(i).axisDir == 'y'
            k = dsearchn(grid.Y,tef(i).Y);
            dx = grid.dyG(k);
        elseif tef(i).axisDir == 'x'
            k = dsearchn(grid.X,tef(i).X);
            dx = grid.dxG(k); 
        else
            error('indicate velDir, u or v direction for velocity transect')
        end
        
        DX = repmat(dx,[1, nz,nt]);

        cellArea = DX.*DZ;
        A = reshape(cellArea,[nx*nz,nt]); %cell areas
    
    %% define salinity class vector
    
        Smin = floor(nanmin(s(:)));
        Smax = ceil(nanmax(s(:)));
        S = linspace(Smin,Smax,N);
        deltaS = (Smax - Smin)/N;
        n = floor((s-Smin)/deltaS); %index salt classes
        tef(i).Sclass = S;
    
    %% salt flux per salinity class
        

        for t = 1:nt %one time step at a time
            ss = squeeze(s(:,t));
            uu = squeeze(u(:,t));
            AA = squeeze(A(:,t));
            nn = squeeze(n(:,t));

            for ii = 1:length(S)
                ind = find(nn == ii);   
                qn(ii,t) = nansum(uu(ind).*AA(ind)); %volume flux per class
                sn(ii,t) = nansum(ss(ind)); %salt per class
                qsn(ii,t) = nansum(uu(ind).*AA(ind).*ss(ind)); %salt flux per class
            end
        end
        tef(i).qn = qn;
        tef(i).qsn = qsn;
        tef(i).sn = sn;
        for ii = 1:length(S)
                qnFilt(ii,:) = godinfilt(squeeze(qn(ii,:)),0,1); %tidally-averaged vol flux
                qsnFilt(ii,:) = godinfilt(squeeze(qsn(ii,:)),0,1); %tidally-averaged salt flux
        end
        tef(i).qnFilt = qnFilt;
        tef(i).qsnFilt = qsnFilt;
        
        for t = 1:nt
                tef(i).Q(:,t) = cumsum(qnFilt(:,t),'reverse'); %Integrated flux across class
                tef(i).Qs(:,t) = cumsum(qsnFilt(:,t),'reverse'); %Integrated salt flux across class      
        end

%% find maxima and minima (USING LORENZ ET AL. 2O19)
        for t = 1:nt
            %find peaks
            
            [mxPks,mxLocs,mxW,mxP] = findpeaks(tef(i).Q(:,t));
            [mnPks,mnLocs,mnW,mnP] = findpeaks(-tef(i).Q(:,t));
            
            %test with video
            f(t) = figure;
            hold on
            plot(tef(i).Q(:,t),S,'k')
            plot(tef(i).Q(mxLocs,t),S(mxLocs),'^k')
            plot(tef(i).Q(mnLocs,t),S(mnLocs),'vk')
            
            %apply threshold
            Qthresh = 0.25*max(abs(tef(i).Q(:,t))); %recommended value in Lorenz 2019 is 0.1*max(Q)
            ind = find(abs(mxP) < Qthresh);
            mxLocs(ind) = [];
            ind = find(abs(mnP) < Qthresh);
            mnLocs(ind) = [];
            
            Iall = false(size(S)); %all extrema (both Max and Min)
            Iall(1) = 1; %Smin and Smax are dividing salinities
            Iall(end) = 1;
            Iall(mxLocs) = 1;
            Iall(mnLocs) = 1;
            div = find(Iall==1); %dividing indices
            Sdiv = S(div);
            tef(i).Sdiv(t).sdiv = Sdiv;%dividing salinities
            
            for ii = 1:length(Sdiv)
                plot([min(tef(i).Q(:)) max(tef(i).Q(:))],[Sdiv(ii) Sdiv(ii)],'--k','LineWidth',1)
            end
            xlim([min(tef(i).Q(:)) max(tef(i).Q(:))])
            ylim([min(Sdiv) max(Sdiv)])
            set(gca,'ydir','reverse')
            
            movie(t) = getframe(gcf);
            cd(Odir)
            if t == 1
                mkdir('Sdiv')
            end
            cd('Sdiv')
            savefig(f(t),['Timestep ', num2str(t)])
            close(f(t))
         
            for ii = 1:length(div)-1
                tef(i).deltaQj(ii,t) = -(tef(i).Q(div(ii+1),t) - tef(i).Q(div(ii),t)); % this gives wrong sign, but seems to do same thing in Lorenz 2019
                tef(i).deltaQsj(ii,t) = -(tef(i).Qs(div(ii+1),t) - tef(i).Qs(div(ii),t));
            end

%% index and combine

            in = find(tef(i).deltaQj(:,t) > 0);
            out = find(tef(i).deltaQj(:,t) < 0);

            tef(i).Qin(t) = sum(tef(i).deltaQj(in,t));
            tef(i).Qout(t) = sum(tef(i).deltaQj(out,t));
            tef(i).Qs_in(t) = sum(tef(i).deltaQsj(in,t));
            tef(i).Qs_out(t) = sum(tef(i).deltaQsj(out,t));
            tef(i).Sin(t) = tef(i).Qs_in(t)/tef(i).Qin(t);
            tef(i).Sout(t) = tef(i).Qs_out(t)/tef(i).Qout(t);

        end
        
        movi = VideoWriter(['TEFline' num2str(i)]);
        movi.FrameRate = 10;
        open(movi)
        writeVideo(movi,movie)
        close(movi)
   end
end
               