% This is an example of how to create a lookup table between CFC-12 age 
% and TTD parameters gamma and delta. In this example, we use a gamma/delta 
% ratio of 1. This seems to fit many ocean observations.

% Look toward the end of this file for an example of how to use the lookup table 
% with data.

% The reason for using pCFC age is that it makes the lookup table independent 
% of any assumptions you want to make about surface saturation. So when you apply 
% the lookup table to tracer data, you apply those assumptions in the tracer age 
% calculation.

% In the following, I use the function 'f_calc_cfc_age' written by Sabine Mecking 
% (mecking@ocean.washington.edu) to compute pCFC age from pCFC. You can replace 
% this by your function or use the simple hack shown below.

% Observational times at which you (eventually) intend to apply the lookup table
Tpred=[1990:.5:1997]';

% Range of gamma and delta values
Gamma=[.01:.01:.04 .05:.05:1 1.1:.1:2 3:1:400 402:2:800]';
Delta=Gamma;
pCFCage=[.01:.01:.09 .1:.1:2 3:1:65]';

nt=length(Tpred);
ng=length(Gamma);
ntp=length(pCFCage);

% Load CFC-12 data
% add your own code here to read CFC data
code to read data ...
% you should now have 3 variables:
% pcfc_atm_nh: vector of northern hemisphere pcfc history
% pcfc_atm_sh: vector of southern hemisphere pcfc history
% Tcfc:        vector of times (in years) at which above history is given

pcfc_pred=repmat(NaN,[ng nt 2]);
pcfc_age_pred=repmat(NaN,[ng nt 2]);

for it=1:nt
  Tp=Tpred(it);
  for ig=1:ng
    gamma=Gamma(ig)
    delta=Delta(ig);
%   Predict the interior pCFC and pCFC age by convolving the surface pCFC history 
%   with TTD(gamma,delta)
%   Northern hemisphere
    pcfc=ttdConvolveWithG1d(Tp,Tcfc,pcfc_atm_nh,gamma,delta); % NH
%   Compute pCFC age from pCFC. You may replace this with your own function or use the 
%   following hack:
%   pcfc_age=Tpred-interp1(pcfc_atm_nh,Tcfc,pcfc); % this gives the same answer  
    [pcfc_year,pcfc_age]=f_calc_cfc_age(pcfc,Tp,2,'atmos_years',Tcfc,'pcfc_atmos',pcfc_atm_nh);
    pcfc_pred(ig,it,1)=pcfc;
    pcfc_age_pred(ig,it,1)=pcfc_age;
%   Southern hemisphere
    pcfc=ttdConvolveWithG1d(Tp,Tcfc,pcfc_atm_sh,gamma,delta); % SH
%   Compute pCFC age from pCFC. You may replace this with your own function or use the 
%   following hack:
%   pcfc_age=Tpred-interp1(pcfc_atm_sh,Tcfc,pcfc); % this gives the same answer      
    [pcfc_year,pcfc_age]=f_calc_cfc_age(pcfc,Tp,2,'atmos_years',Tcfc,'pcfc_atmos',pcfc_atm_sh);
    pcfc_pred(ig,it,2)=pcfc;
    pcfc_age_pred(ig,it,2)=pcfc_age;
  end
end

% So far we have a lookup table for pCFC age in terms of gamma. We now reverse this 
% to make a lookup table for gamma in terms of pCFC age and Tpred
GAMMA=repmat(NaN,[ntp nt 2]);
for it=1:nt
  GAMMA(:,it,1)=interp1(pcfc_age_pred(:,it,1),Gamma,pCFCage);
  GAMMA(:,it,2)=interp1(pcfc_age_pred(:,it,2),Gamma,pCFCage);
end

save cfc12_age_ttd_lookup_table Tcfc pcfc_atm_nh pcfc_atm_sh ...
     Gamma Delta pcfc_pred pcfc_age_pred Tpred pCFCage GAMMA

% Example of how to use the lookup table given pCFC age data in the vector 'pcfc_age' observed 
% at corresponding times Tobs (in decimal years): 

% agetable=load('cfc12_age_ttd_lookup_table'); % load the age table
% % If the water was mostly ventilated in the northern hemisphere:
% gamma=interpn(agetable.pCFCage,agetable.Tpred,agetable.GAMMA(:,:,1),pcfc_age,Tobs);
% % If the water was mostly ventilated in the southern hemisphere:
% gamma=interpn(agetable.pCFCage,agetable.Tpred,agetable.GAMMA(:,:,2),pcfc_age,Tobs);

% The variable 'gamma' is the predicted value of gamma (assuming a gamma/delta ratio of 1). 
% You can now convolve the 1-d TTD parameterized by gamma and delta(=gamma) with the surface 
% boundary condition for any other tracer.