function [tz tn sz sn Ntz Nsz rp]=kh_get_all_subduction_terms2(s1,s2,t1,t2,p1,p2)
%
% Gets all terms required for the calculation of pure
% warming/freshening/heave.  Also get mean Rp term.
%
% input =depths*profiles.
% output=depths*profiles.
%       for output definitions see Bindoff (1994) (eqn17)
%
% This version doesn't use alpha or beta
%
% Written by Kieran Helm (23-5-2007)

t1=sw_ptmp(s1,t1,p1,0); % Convert to potential temperature
t2=sw_ptmp(s2,t2,p2,0);

tz=kh_tz(s1,s2,t1,t2,p1,p2);
tn=kh_tn(s1,s2,t1,t2,p1,p2);
sz=kh_sz(s1,s2,t1,t2,p1,p2);
sn=kh_sn(s1,s2,t1,t2,p1,p2);
Ntz=tn-tz; % uses difference method 
Nsz=sn-sz;
% Ntz2=kh_N_alpha_tz(s1,s2,t1,t2,p1,p2);% uses full method (OUTPUT)
% Nsz2=kh_N_beta_sz(s1,s2,t1,t2,p1,p2);

%% Calc mean Rp
warning off
rp=kh_Rp4([(s1+s2)/2],[(t1+t2)/2],[(p1+p2)/2]);


%% Calc difference Ntz from different methods - as a check
% pd=100*(Ntz2-Ntz)./Ntz; % percentage difference from full method
% pd2=reshape(pd(find(pd>-100 & pd<100)),1,[]);
% hist(pd2,200); grid; std(pd2)
% title('Percentage difference between full and subtraction method of N*alpha*tz');
% xlabel('Percentage difference from full method'); ylabel('number of observations')
% plot([230:2:258 260:284]/10,kh_nanmean(pd,2),'r*-'); grid
% title(' Average percentage difference between full and subtraction method of N*alpha*tz by gamma layer')

