function [Aw Af Ah Aw_r2 Af_r2 Ah_r2]=kh_subduction_coeff(s1,s2,t1,t2,p1,p2);
% function [Aw Af Ah Aw_r2 Af_r2 Ah_r2]=kh_subduction_coeff(s1,s2,t1,t2,p1,p2);
% 
% Calculates the dimensionless coefficients Aw, Af, and Ah which represent
% the strengths of warming, freshening, and heave (positive is upwards
% shoaling) using a single process overdetermined solution.
%
% As in Bindoff (1994) Eqn 17 and Bindoff (2000) Eqn A6
%
% s1,s2 are the historical and modern salinity on neutral density
% surfaces*number of profiles
% t1,t2 are insitu density
% p1,p2 are pressure
% w_lat, w_lon=1*number of profiles
% r2 is the proportion of variance explained by each of the processes as in
% equation 18 Bindoff(1994)
%
% Produces solution using Gaussian elimination when xMa=d (x,=1/(Rp-1), m=rp matrix,a=aw
% terms, and d=6'observable quantities').  Uses a=xM\d where M is reduced
% down to a single column as two of Aw,Af,Ah are set to zero (and hence
% becomes an over determined problem).
%
% x bit added where x=1/(Rp-1).  The density bit (p^-1P^-1|z) has been left off
% but is equal to -alphaT'|z+betaS'|z.  Should report Aw, Af and Ah as
% density anomalies as on p1149 bindoff (1994). (24-9-2007) Kieran Helm
clc


[rows cols]=size(t1);
Aw=NaN(rows,cols); Af=NaN(rows,cols); Ah=NaN(rows,cols); 
Aw_r2=NaN(rows,cols); Af_r2=NaN(rows,cols);  Ah_r2=NaN(rows,cols);

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

%% Get terms
tz=kh_alpha_tz(s1,s2,t1,t2,p1,p2);
tn=kh_alpha_tn(s1,s2,t1,t2,p1,p2);
sz=kh_beta_sz(s1,s2,t1,t2,p1,p2);
sn=kh_beta_sn(s1,s2,t1,t2,p1,p2);
% Ntz=tn-tz; % difference method 
% Nsz=sn-sz;
Ntz=kh_N_alpha_tz(s1,s2,t1,t2,p1,p2);% uses full method
Nsz=kh_N_beta_sz(s1,s2,t1,t2,p1,p2);
comment='have terms'

warning off
%% Calc mean Rp
rp1=kh_Rp4(s1,t1,p1);
rp2=kh_Rp4(s2,t2,p2);
rp=(rp1+rp2)/2;
rp(isinf(rp))=NaN;  % remove all infinity point (on Smax and Smin) ( 24-9-2007)
comment='have rp'

%% Calculate subduction coefficients
for c=1:cols % profile at a time
    for r=1:rows % density level at a time
        x=1/(rp(r,c)-1);
        M=[-(rp(r,c)-1) 0 -rp(r,c);...
            1 rp(r,c) 0;...
            rp(r,c) rp(r,c) rp(r,c);...
            0 rp(r,c)-1 -1;...
            1 rp(r,c) 0;...
            1 1 1];
        M=x*M; % added coefficient x out the front (24-9-2007)
        d=[tz(r,c); tn(r,c); Ntz(r,c); sz(r,c); sn(r,c); Nsz(r,c)];
   
        Aw(r,c)=M(:,1)\d;    % (Af=Ah=0)
        Af(r,c)=M(:,2)\d;    % (Aw=Ah=0)
        Ah(r,c)=M(:,3)\d;    % (Aw=Af=0)
        
% Calculate R^2
    % dest=M*[Aw(r,c); Af(r,c); Ah(r,c)]; % for all
    Aw_dest=M*[Aw(r,c); 0; 0];
    Af_dest=M*[0; Af(r,c); 0];
    Ah_dest=M*[0; 0; Ah(r,c)];
    Aw_r2(r,c)=(Aw_dest'*Aw_dest)/(d'*d);
    Af_r2(r,c)=(Af_dest'*Af_dest)/(d'*d);
    Ah_r2(r,c)=(Ah_dest'*Ah_dest)/(d'*d);

    end
     
if  rem(c,1000)==0
    comment=c
end  

end

warning on all