% ps08_1.m is an m-file that runs and plots the Abyssal Carbon which uses 14C in
%       a 1_D model to estimate upwelling rates and turbulent diffusivities.
%
% Modif'd 2002:11:20 D. Glover for 2002

clear all;
close all;

% load the data

load hydro.mat
load c14.mat

'To see a T-S plot of all the data, press any key.'
pause

figure(1)
plot(s,t,'m*')
title('T-S Plot')
xlabel('Salinity')
ylabel('Temperature')
'To see temperature and salinity plotted separately, press any key.'
pause

figure(2)
subplot(1,3,1)
plot(t,z,'b*')
axis('ij')
title('T vs. Z')
ylabel('Depth')
xlabel('Temperature')
subplot(1,3,2)
plot(s,z,'r*')
axis('ij')
title('S vs. Z')
ylabel('Depth')
xlabel('Salinity')
subplot(1,3,3)
plot(s,t,'m*')
title('T-S Plot')
xlabel('Salinity')
ylabel('Temperature')

'To see the T-S linear sub-range, press any key.'
pause

figure(1)
hold on
all=[z t s o];
dat=all(13:23,1:4);
plot(dat(:,3),dat(:,2),'b*')
pause(1)
text(34.5,5,'Adjusting to linear axis  range...')
pause(1)
figure(3)
plot(dat(:,3),dat(:,2),'b*')
title('T-S Plot. T-S are ~ Linear from 1226.1 m to 3249.7 m')
xlabel('Salinity')
ylabel('Temperature')
'To see the salinity  vs. depth subrange and the fit, press any key.'
pause

to=dat(:,2);
so=dat(:,3);
zsr=(-1*(dat(:,1)-3.2497e3));
figure(4)
plot(so,zsr,'r*')
set(gca,'Ytick',[0:200:2200]);
set(gca,'Xtick',[34.35:.05:34.75]);
axis([34.35 34.75 0 2100])
xlabel('Salinity')
ylabel('Subrange Depth')
title('Salinity vs. Subrange Depth')
pause(2)
m=ones(2,1);
pin=[500];
[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=nlleasqr(zsr,so,pin,'salfunc');
Salp=p;
Salsp = sqrt(diag(covp));
Salcovp=covp;
hold on
plot(f,zsr,'r--')
%text(34.45,1000,sprintf('Z* = %3.2f',p(1)))
text(34.45,1000,sprintf('Z* = %3.2f +- %3.2f',Salp,Salsp))

m(1)=p(1);
'To see the temperature vs depth subrange and the fit, press any key'
pause

figure(5)
plot(to,zsr,'b*')
set(gca,'Ytick',[0:200:2200]);
set(gca,'Xtick',[1.5:.5:4]);
axis([1.5 4 0 2100]);
ylabel('Subrange Depth')
xlabel('Temperature')
title('Temperature vs. Subrange Depth')
pause(2)
figure(5)
pin=[500];
[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2] = nlleasqr(zsr,to,pin,'tempfunc');
Tempp=p;
Tempsp = sqrt(diag(covp));
Tcovp=covp;
hold on
plot(f,zsr,'b--')
text(2.5,1000,sprintf('Z* = %3.2f +- %3.2f',Tempp,Tempsp))
m(2)=p(1);
Zstar=mean(m);
sZstar=sqrt(Salsp^2+Tempsp^2);
text(2.5,800,sprintf('mean Z* = %3.2f +- %3.2f',Zstar,sZstar));
['To plot the oxygen data over the depth subrange and do the fit, press any key.']
pause

figure(6)
oo=dat(:,4);
plot(oo,zsr,'c*')
title('Oxygen vs. Depth Subrange')
xlabel('Oxygen (umol/kg)')
ylabel('Depth Subrange')
pause(2)
pin=[-0.01];
[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2] = nlleasqr(zsr,oo,pin,'oxyfunc');
Oxyp=p;
Oxysp = sqrt(diag(covp));
Ocovp=covp;
hold on
plot(f,zsr,'c--')
text(160,1000,sprintf('J*(O2) = %1.3e +- %1.3e',Oxyp, Oxysp))
JO2=p(1);
JstarO2=Oxyp
JstarCO2=-.62*JstarO2
ErrJstarCO2=sqrt((-.62*Oxysp)^2)
Jstar14CO2=.95*JstarCO2
ErrJstar14CO2=sqrt((.95*ErrJstarCO2)^2)
%ErrOx=Oxysp/Oxyp
%ErrJstar14CO2=Jstar14CO2*ErrOx
['To see a plot of del14C and tCO2 vs. the depth subrange, press any key.']
pause

figure(7)
subplot(1,2,1)
zsrc=2023.6:-224.8:0;
zsrc=zsrc';
plot(del14,zsrc,'r+')
title('Del14C vs Depth Subrange')
xlabel('Del 14C')
ylabel('Depth Subrange')
subplot(1,2,2)
plot(tco2,zsrc,'b+')
title('Total CO2 vs. Depth Subrange')
xlabel('Total CO2')
ylabel('Depth Subrange')
'To see a plot of ERC vs depth subrange, press any key.'
pause

figure(8)
c14c=(1+del14./1000).*tco2;
plot(c14c,zsrc,'g*')
title('Equivalent Radiocarbon Concentration vs. Depth Subrange')
xlabel('Equivalent Radiocarbon Concentration')
ylabel('Depth Subrange')
%['To do a non-linear fit to the data with both Co and Cm relaxed, press any key']
%pause

figure(9)
lam=3.833e-12;
%pin=[70000 1844 1917.4];
%clear p
%[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=nlleasqr(zsrc,c14c,pin,'radiotst1');
%plot(c14c,zsrc,'g*')
%w9p=p;
%w9sp = sqrt(diag(covp));
%w9covp=covp;
%hold on
%plot(f,zsrc,'m--')
%text(1830,2250,sprintf('w* = %1.3e',p(1)))
%w=p(1)*lam;
%J=w*JO2;
%K=w*Zstar;
%text(1850,900,sprintf('w = %1.2e (m/s)',w))
%text(1850,650,sprintf('J = %1.2e (umol/kg/s)',J))
%text(1850,400,sprintf('K = %1.2e (m^2/s)',K))
%title('No constraint')
%pw=p(1);
%clear p
%'To do a non-linear fit with Cm constrained and Co relaxed, press any key.'
%pause

%figure(10)
%[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=nlleasqr(zsrc,c14c,pin,'radiotst2');
%w10p=p;
%w10sp = sqrt(diag(covp));
%w10covp=covp;
%plot(c14c,zsrc,'g*')
%hold on
%plot(f,zsrc,'r--')
%text(1830,2250,sprintf('w* = %1.3e',p(1)))
%w=p(1)*lam;
%J=w*JO2;
%K=w*Zstar;
%text(1850,900,sprintf('w = %1.2e (m/s)',w))
%text(1850,650,sprintf('J = %1.2e (umol/kg/s)',J))
%text(1850,400,sprintf('K = %1.2e (m^2/s)',K))
%title('Cm constrained')
%pm=p(1);
%clear p
%'To do a non-linear fit with Co constrained and Cm relaxed, press any key.'
%pause

%figure(11)
%[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=nlleasqr(zsrc,c14c,pin,'radiotst3');
%w11p=p;
%w11sp = sqrt(diag(covp));
%w11covp=covp;
%plot(c14c,zsrc,'g*')
%hold on
%plot(f,zsrc,'b--')
%text(1830,2250,sprintf('w* = %1.3e',p(1)));
%w=p(1)*lam;
%J=w*JO2;
%K=w*Zstar;
%text(1850,900,sprintf('w = %1.2e (m/s)',w))
%text(1850,650,sprintf('J = %1.2e (umol/kg/s)',J))
%text(1850,400,sprintf('K = %1.2e (m^2/s)',K))
%title('Co Constrained')
%po=p(1);
%clear p
%'To do a non-linear fit with both Co and Cm constrained, press any key.'
%pause

%figure(12)
pin=[80000];
[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=nlleasqr(zsrc,c14c,pin,'radiotst4');
w12p=p;
w12sp = sqrt(diag(covp));
w12covp=covp;
plot(c14c,zsrc,'g*')
hold on
plot(f,zsrc,'c--')
text(1830,2250,sprintf('w* = %1.3e +- %1.3e',w12p, w12sp))
w=p(1)*lam;
J=w*JO2;
K=w*Zstar;
xlabel('Equivalent Radiocarbon Concentration');
ylabel('Depth Subrange');
text(1850,900,sprintf('w = %1.2e (m/s)',w))
text(1850,600,sprintf('J = %1.2e (umol/kg/s)',J))
text(1850,300,sprintf('K = %1.2e (m^2/s)',K))
title('Both Constrained')
pb=p(1);


pin=[80000];
[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=nlleasqr(zsrc,c14c,pin,'radiotst5');
w12p=p;
w12sp = sqrt(diag(covp));
w12covp=covp;
%plot(c14c,zsrc,'g*')
hold on
%plot(f,zsrc,'c--')
text(1830,2100,sprintf('w* = %1.3e +- %1.3e upper limit',w12p, w12sp))
w=p(1)*lam;
J=w*JO2;
K=w*Zstar;
text(1850,800,sprintf('w = %1.2e (m/s) upper limit',w))
text(1850,500,sprintf('J = %1.2e (umol/kg/s) upper limit',J))
text(1850,200,sprintf('K = %1.2e (m^2/s) upper limit',K))
title('Both Constrained')
pb=p(1);


pin=[80000];
[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=nlleasqr(zsrc,c14c,pin,'radiotst6');
w12p=p;
w12sp = sqrt(diag(covp));
w12covp=covp;
%plot(c14c,zsrc,'g*')
hold on
%plot(f,zsrc,'c--')
text(1830,1900,sprintf('w* = %1.3e +- %1.3e lower limit',w12p, w12sp))
w=p(1)*lam;
J=w*JO2;
K=w*Zstar;
text(1850,700,sprintf('w = %1.2e (m/s) lower limit',w))
text(1850,400,sprintf('J = %1.2e (umol/kg/s) lower limit',J))
text(1850,100,sprintf('K = %1.2e (m^2/s) lower limit',K))
title('Both Constrained')
pb=p(1);




% Take ranges of w, J and K divide by 2


% w* +- 1.977e4
% w +- 7.55e-8
% J +- 1.33e-9
% K +- 3.905e-5