% Sets up matfile with subduction coefficients (for faster access).  Two
% versions - one with the top 100m the other with this removed.
%
% This version uses for method for the calculation of Ntz and Nts in
% kh_subduction_coeff script.
%
% Kieran Helm (27-9-2007)
clear; clc
cbar=[0:0.1:1];

%% Load data
d=[230:2:258 260:284];
[t1,t2,w_lat,w_lon]=kh_load_nddata('ns4',d,'t'); %t1=historical, t2=modern
[p1,p2,w_lat,w_lon]=kh_load_nddata('ns4',d,'p');
[s1,s2,w_lat,w_lon]=kh_load_nddata('ns4',d,'s');

%% Get subduction coeff (use in situ temp) - inc top 100 m
[Aw Af Ah Aw_r2 Af_r2 Ah_r2]=kh_subduction_coeff(s1,s2,t1,t2,p1,p2); 
save('subduction_coeff_inc100','Aw','Af','Ah','Aw_r2','Af_r2','Ah_r2');

%% Get subduction coeff 2(use in situ temp) - remove top 100 m
a=find(p1<100); p1(a)=NaN; s1(a)=NaN; t1(a)=NaN; 
a=find(p2<100); p2(a)=NaN; s2(a)=NaN; t2(a)=NaN; 
[Aw Af Ah Aw_r2 Af_r2 Ah_r2]=kh_subduction_coeff(s1,s2,t1,t2,p1,p2); 
save('subduction_coeff','Aw','Af','Ah','Aw_r2','Af_r2','Ah_r2');