%metzyme_rates.m
%22 nov 2013
%asantoro@umces.edu
%additional comments and clarifications added for BCO-DMO submission, 25
%mar 2016
%calculates the nitrification rates for three stations in the equatorial
%pacific during a cruise in october 2011. 15NH4 was added to to duplicate bottles,
%plus a no addition control.
%this script must be run separately for each station. the input files are
%hard-coded and must be edited manual as indicated below between lines marked with '*****'.
%this script calls the function nitr_rate_rev.m and requires the MATLAB optimization
%toolbox (for lsqcurvefit.m).
%bottle numbers:
%1-15NH4+ (bottle A)
%2-15NH4+ (bottle B)
%3-no addition (bottle C)
%=====================================
clear all
%load data on initial conditions, columns:
%1-depth (m)
%2-AP15NH4, starting atom fraction 15N (15N/(15N+14N) calculated from mass balance of tracer addition and starting [NH4+]
%3-No_15 (NO3-), the initial concentration (uM) of 15N-NO3-
%4-No_14 (NO3-), the initial concentration (uM) of 14N-NO3-
%********************
%select input files here (initials and isotope data) by uncommenting the
%appropriate files
initials = load('initials_mzst1.csv');
%initials = load('initials_mzst3.csv');
%initials = load('initials_mzst5.csv');
%load isotope data. columns: 1-depth, 2-bottle# (see above), 3-time(h), 4-15R-NO3
data = load('data_mzst1.csv');
% data = load('data_mzst3.csv');
%data = load('data_mzst5.csv');
%***********************
depnum = length(initials); %number of depths
botnum = 3; %number of treatments
%set curve fitting parameters
x0 = [0.00001 0.0001]; %initial guess for k and F
lb = 0; %sets lower bound of zero, i.e. k and F are positive
ub = []; %upper bound set to one, change to empty matrix if there is a fitting error
%set up figure for data and fitting results
figure
%setup matrix for data
metzymerates = zeros(depnum,5);
%break data into matrices by depth and bottle, then fit rate function
for i=1:depnum %for each depth
dep = initials (i,1); %temporary depth
eval(['data_' num2str(dep) '= data (find (data(:,1) == dep),:);']) %break into matrix by depth
tempdata = eval(['data_' num2str(dep)]); %set tempdata to current depth
c = initials (i,2:4); %specify coefficients for rate fitting
for b=1:botnum %for each bottle
bot = tempdata (find (tempdata(:,2) == b),:); %find the data for bottle = botnum
time = bot(:,3); %time data for fitting
ydata = bot(:,4); %15R data for fitting
%use least squares curve fitting to calculate a, the vector of unknown
%coeffs. a(1) = k, a(2) = F. see nitr_rate.m
[a,resnorm,r,exitflag,output,lambda,J] = lsqcurvefit(@(a,time) nitr_rate_rev(a,time,c), x0,time,ydata,lb,ub);
fit = nitr_rate_rev(a,time,c);
subplot(depnum,botnum,b+(i-1)*botnum) %plot
plot(time,ydata,'*') %real data
hold on
plot(time,fit,'r-') %fit
rt = a(2)*1000*24; %convert rate to nM per day
k_removal = a(1) *24; % convert rate constant to per day
% calculate standard error of fit coefficient
% approximate the Hessian using the Jacobian, and use the inverse of that
% as an approximation to the covariance matrix:
%
[Q,R] = qr(J,0); %decomposition of the jacobian
mse = sum(abs(r).^2)/(size(J,1)-size(J,2));
Rinv = inv(R);
Sigma = Rinv*Rinv'*mse;
SE_uM = sqrt(diag(Sigma)); %SE in uM per h
SE_nM = SE_uM(2)*1000*24; %SE in nM per day
%
% where r is the residual vector and J is the Jacobian, both outputs of
% LSQCURVEFIT.
%build matrix
metzymerates (b+(i-1)*botnum,1) = dep; %first column of rate matrix is depth
metzymerates (b+(i-1)*botnum,2) = b; %second column is bottle #
metzymerates (b+(i-1)*botnum,3) = rt; %third column is rate in nM per day
metzymerates (b+(i-1)*botnum,4) = SE_nM; %fourth column is SE of the rate in nM per day
metzymerates (b+(i-1)*botnum,5) = a(1); %fifth column is the removal rate constant in units of per hour
end
end
%output graph
print -painters -dpng metzyme_rates_st1.png
%output rates to csv file
dlmwrite('metzyme_rates_st1.txt', metzymerates, 'delimiter', '\t')