function [Nalpha]=kh_N_alpha_tz(si,sf,ti,tf,pi,pf)
% function [Nalpha]=kh_alpha_tz(si,sf,ti,tf,pi,pf)
%
% Caluclates N'*alpha*ptmp|z (movement of isopycnal*alpha*potential 
% tempe change on pressure surfaces
% All terms need to be entered on neutral density levels.
%
% si=salinity initial (neutral density levels*profiles)
% sf=salinity final etc.
% ti= initial potnetial temperature (needs to be potential temperature)
%
% N=difference in height between two surfaces where + value =upward
% movement of isopycnal.  Note that to be consistent a gradient where the
% water is warm at the surface and cold below is considered a positive
% gradient.
%
% Designed for use in calculating terms from Bindoff and McDougall (1994).
% Equation from B&M (1994) p1143 (equation 14)
%
% Written by Kieran Helm (17-5-2007)

%% Approach 1: Difference approach (EQN 11 - bindoff and mcd)
% [tn]=kh_alpha_tn(si,sf,ti,tf,pi,pf);
% [tz]=kh_alpha_tz(si,sf,ti,tf,pi,pf);
% Nalpha=tn-tz;


%% Approach 2: Full approach (Given in EQN 14)
N=pi-pf; % Therefore + = shoaling of isopycnal (to agree with BM1994 p1143) - added 25-9-07
[sm tm pm]=kh_mid_points(si,sf,ti,tf,pi,pf); % Get midpoints
al=sw_alpha(sm,tm,pm,'ptmp'); % Calculate alpha

    %% Calculate tfz-tiz
    warning off % for divide by zero
    [m n]=size(ti);
    gi=NaN(m,n); gf=NaN(m,n); % gradients of temp/depth slope

    % gradient (dt/dz) intial and final profiles
    % note: p*-1 is to reverse xaxis (as looking from depth to surface)
    for k=1:n
         gi(:,k)=gradient(ti(:,k),pi(:,k)*-1); 
         gf(:,k)=gradient(tf(:,k),pf(:,k)*-1);
    end

%% If gradient too steep then replace with NaN
gi(find(gi>0.5 | gi<-0.5))=NaN;
gf(find(gf>0.5 | gf<-0.5))=NaN;

%% combine all elements of equation 14
Nalpha=0.5*N.*al.*[gf+gi];
 

%% Old script to get t(mp) - not used
% upi=NaN(m,n); upf=NaN(m,n); % overlying pressures
% uti=NaN(m,n); utf=NaN(m,n); % overlying temperature
% for r=2:40 
%     upi(r,:)=pi(r-1,:);        upf(r,:)=pf(r-1,:);
%     uti(r,:)=ti(r-1,:);        utf(r,:)=tf(r-1,:);
% end
% ti_pm=gi.*(pm-upi)+uti; 
% tf_pm=gf.*(pm-upf)+utf; 