clear all
close ALL
% Initialize:
dt=273.15;
g=9.81;
al1=2.605e-5;             % Bouhifd et al. GRL
al2=0.758e-8;
cp1=-402.753e3/77.31;     % Gillet et al. JGR
cp2=74.29e3/77.31;
cp3=87.588e6/77.31;
cp4=-25.913e9/77.31;
cp5=25.374e11/77.31;
dz=100;
nn=3000;
for i=1:nn
  dep(i)=dz*(i-1);
end
tmp=zeros(size(dep));
spi=sqrt(pi);
% Loop twice:
for n=1:2
  disp(sprintf('Geotherm %1i:',n))
  disp('Input variables [Suggested values in brackets]:')
  A0=input('Surface heat production H0 (micro-W/m^2) [1.0]: ');
  A0=A0*1e-6;
  atc=input('Crustal thermal conductivity param A [0.28]: ');
  btc=input('Crustal thermal conductivity param B [0.000316]: ');
  radl=input('Radiogenic length parameter l_rad (km) [10]: ');
  radl=radl*1e3;
  tc=input('Crustal thickness (km) [40]: ');
  tc=tc*1e3;
  nnc=round(tc/dz)+1;
  nnm=nn-nnc;
  lcon=input('Conductive thermal length l_con (km) [70]: ');
  lcon=lcon*1e3;
  Tr=input('Mantle asthenospheric temperature (°C) [1350]: ');
  Tr=Tr+dt;
  Ts=input('Surface temperature (°C) [0]: ');
  Ts=Ts+dt;
  condm=2.3;
  qs=A0*radl*(1-exp((-tc)/radl))*(1-2*radl/(spi*lcon))+2*(Tr-Ts)/(spi*lcon*(atc+Ts*btc));
  for i=1:nnc
    tmp(i)=Ts+(atc*A0*radl*radl*(1-exp(-dep(i)/radl)-erf(dep(i)/lcon))+(Tr-Ts)*erf(dep(i)/lcon))/(1-btc*A0*radl*radl*(1-exp(-dep(i)/radl)-erf(dep(i)/lcon)));
  end
  Tm=tmp(nnc);
  qm=(Tm-tmp(nnc-1))/(dz*(atc+btc*Tm));
  lcm=2*condm*(Tr-Tm)/(spi*qm);
  for i=1:nnm
    tmp(nnc+i)=Tm+(Tr-Tm)*erf(dep(i)/lcm);
  end
% Correct potential temperature to true:
  ttm=0;
  for i=2:nn
    t1=tmp(i-1);
    t2=t1*t1;
    t3=t2*t1;
    ttm=ttm+g*dz*(al1*t1+al2*t2)/(cp1+cp2*log(t1)+cp3/t1+cp4/t2+cp5/t3);
    tmp(i)=tmp(i)+ttm;
  end
% Products:
  disp('    ')
  disp(sprintf('Surface heat flow: %6.1f',qs*1e3))
  disp(sprintf('Mantle heat flow: %6.1f',qm*1e3))
  disp(sprintf('Moho temperature: %7.1f',tmp(nnc)-dt))
  disp('    ')
% Plots:
%AAAAAAAAAAAAA
  if n==1
    tkp=tmp;
    a=plot(tmp-dt*ones(size(tmp)),dep/1000,'-');
    set(a,'LineWidth',5,'Color','Red')
    set(gca,'YDir','reverse')
    hold on
    b=gca;
    set(b,'FontSize',16)
    a=xlabel('Temperature (°C)');
    a=ylabel('Depth (km)');
  else
    a=plot(tmp-dt*ones(size(tmp)),dep/1000,'-');
    set(a,'LineWidth',5,'Color','Blue')
    print -dpng Geotherm.png
    figure(2)
    a=plot(tkp-tmp,dep/1000,'-');
    set(a,'LineWidth',5,'Color','Green')
    set(gca,'YDir','reverse')
    b=gca;
    set(b,'FontSize',16)
    a=xlabel('Temperature Difference (Red Minus Blue, °C)');
    a=ylabel('Depth (km)');
    print -dpng T_Diff.png
  end
end
