% age2exhume_readInput_v2.m
 
% Calculates the steady exhumation rate associated with different 
% thermochronometers from a an age, crustal properties,
% and assumed thermal properties, including an offset from a regional
% average elevation.

% Input file format: 
% Text file with 1 header row of eight columns, set up as the following: 
% 1-ID number (no letters), 2-latitude(dd), 3-longitude(dd), 4-elevation(km),
% 5-h-value(km), 6-age(Ma), 7-age uncertainty(Ma),
% 8-thermochronologic system (see below)

% System 1=AHe, 2=AFT, 3=ZHe, 4=ZFT, 5=ArAr

% To avoid problems, use a simple count as ID numbers
% and input latitude and longitude in decimal degrees (no spaces).

% ** Note that the first three columns are just used to generate an output
% text file at the end that will keep sample ID location information
% associated with the calculated exhumation rates.

% Results are given in "output" (without headers) or "output2" (with
% headers).
 
% Created by P. van der Beek and T. Schildgen, last modified 21 November
% 2022.
% Based on earlier code published in Schildgen et al. (2018).
 
clear all;
 
% *******************************************************
% DEFINE PARAMETERS
% ******************************************************* 

% Full example dataset: filename = 'Himalaya_mat_input.txt'
filename = 'Himalaya_mat_input.txt';  % Input file name

% Site parameters
T0 = 25;                    % Temperature (deg C) at elevation = 0 km
lapse = 5;                  % Temperature lapse rate (deg C km-1)
 
% Crustal properties
G_init = 25;                % initial geotherm, deg C km-1
kappa = 30;                 % thermal diffusivity, km2 My-1
L = 30;                     % model thickness, km

% *******************************************************
% ******************************************************* 

R = 8.3145;             % Gas law constant, J mol-1 K-1

% Diffusive / thermal properties
A = 27;                 % geometry factor
 
% Thermochronologic systems and associated parameters
% AHe system, from Farley, 2000, JGR
Ea1 = 1.38e5;            % activation energy, J mol-1
omega1 = 5.00e5;         % Do/a2, s-1
Tc_init1 = 70;           % initial closure temp, deg C (not important,
                         % iterations eliminate effect of this value)
 
% AFT system, from Reiners et al., 2006                  
Ea2 = 1.47e5;            % activation energy, J mol-1
omega2 = 3.73e4;         % Do/a2, s-1
Tc_init2 = 120;          % initial closure temp, deg C (not important)                   
                        
% ZHe system, from Reiners et al., 2004                 
Ea3 = 1.69e5;            % activation energy, J mol-1
omega3 = 4.60e3;         % Do/a2, s-1
Tc_init3 = 180;          % initial closure temp, deg C (not important)
                        
% ZFT system, from Reiners et al., 2006                  
Ea4 = 2.08e5;            % activation energy, J mol-1
omega4 = 1.82e6;         % Do/a2, s-1
Tc_init4 = 220;          % initial closure temp, deg C (not important) 

% Muscovite Ar/Ar system, from Hames and Bowring 1994                
Ea5 = 1.80e5;            % activation energy, J mol-1
omega5 = 6.4e-1;         % Do/a2, s-1
Tc_init5 = 350;          % initial closure temp, deg C (not important) 
 
% Read input file and use to define input values
fid = fopen(filename) ;
data = textscan(fid,'%f %f %f %f %f %f %f %f','HeaderLines',1);
data = cell2mat(data);
fclose(fid);

output = data;
zeros = zeros(height(output),1);
output = horzcat(output,zeros,zeros,zeros,zeros,zeros);
elev_all = data(:,4);
elev_avg = mean(elev_all);
T_avg = T0 - (lapse*elev_avg);

for i = 1:height(data)
    elev = data(i,4);
    delta_h = data(i,5);
    age_1s = data(i,7);
    sys = data(i,8);
    
    if sys == 1
        Ea = Ea1;
        omega = omega1;
        Tc_C = Tc_init1;
    end
    if sys == 2
        Ea = Ea2;
        omega = omega2;
        Tc_C = Tc_init2;
    end
    if sys == 3
        Ea = Ea3;
        omega = omega3;
        Tc_C = Tc_init3;
    end
    if sys == 4
        Ea = Ea4;
        omega = omega4;
        Tc_C = Tc_init4;
    end
    if sys == 5
        Ea = Ea5;
        omega = omega5;
        Tc_C = Tc_init5;
    end

    T_s = T0 - (lapse*(elev-delta_h));  % surface temperature

    % Solve for three cases: age-uncertainty, age, and age+uncertainty
    for j = [-1,0,1]
        age = data(i,6);
        age = age + (j*age_1s);
        if age <= 0
            age = 0.001;
        end
        
        % estimate initial closure depth and exhumation rate
        zc = (Tc_C - T_s)/G_init;   % closure depth; km
        e_rate=(zc + delta_h)/age;  % initial exhumation rate; km/My = mm/y
        dTdt = e_rate*G_init;       % cooling rate, deg C My-1
        TL = T_avg + (L*G_init);    % temperature at base of crust
        Tc_K = Tc_C + 273;          % closure temperature in Kelvin
  
        e_rate_old=0.;     % set initial values
        iter=0;            % iterations
        x = 1.e-3;         % threshold for convergence
        itermax = 1.e4;    % maximum # of iterations

        while abs((e_rate-e_rate_old)/e_rate)>x & iter < itermax % exhumation rate ratio must be within x 
 
            e_rate_old = e_rate;
            zc_old = zc;
            iter = iter+1;
    
            % Perturbation of geotherm due to exhumation
            zc = -(kappa/e_rate)*log(1-(((Tc_C-T_s)/(TL-T_s))*(1-exp(-L*e_rate/kappa))));

            % Solution to the geotherm at the closure temperature:
            G_zc = (((e_rate*(TL-T_s))/(kappa*(1-exp(-L*e_rate/kappa)))))*exp(-(zc*e_rate)/kappa);
            dTdt = abs(G_zc*e_rate);
 
            % Recalculate closure temperature using updated cooling rate
            dTdt_K_s = dTdt/(1000000*365.24*24*3600);
            tau = ((R*Tc_K^2)/(Ea*dTdt_K_s));
            Tc_K = Ea/(R*log(A*tau*omega));
            Tc_C = Tc_K - 273;
    
            % Readjust closure depth following newly caluclated closure temperature
            zc = -(kappa/e_rate)*log(1-(((Tc_C-T_s)/(TL-T_s))*(1-exp(-L*e_rate/kappa))));

            % update exhumation rate
            e_rate=(zc + delta_h)/age;       
            e_rate=(e_rate + e_rate_old)/2;
    
            if e_rate < 0
                delta_age=delta_h/e_rate_old;
                e_rate= zc/(age+delta_age);
            end
        end
        % clean up any lingering imaginary bits
        e_rate = real(e_rate);
        output(i,(10+j)) = e_rate;
        if j == 0
            G_z0 = (((e_rate*(TL-T_s))/(kappa*(1-exp(-L*e_rate/kappa)))))*exp(-(0*e_rate)/kappa);
            output(i,12) = real(G_z0);
            output(i,13) = real(zc);
        end
    end
end


%header = {'ID','lat','long','elev_km','h_val_km','age_Ma','age_unc','system','erate_high','erate_mid','erate_low'};
header = {'ID','lat','long','elev_km','h_val_km','age_Ma','age_unc','system','erate_high','erate_mid','erate_low','G_z0','Zc'};
output2 = [header; num2cell(output)];

% Uncomment to write output (without headers) and/or output2 (with headers) to file
%writecell(output2,'output_cell.txt','Delimiter','space')
%writematrix(output,'output_matrix.txt','Delimiter','space')

