%mars_geotherm_for_clathrate_stability

% Simulating the 1-D Diffusion equation (Fourier's equation) by the
...Finite Difference Method(a time march)
% Numerical scheme used is a first order upwind in time and a second order
...central difference in space (both Implicit and Explicit)

%%
%Specifying Parameters
Q_lith = 0.03; %W/m2
k_cond = 2.5; %k_cond
nt=length(annualmeanT);               %Number of time steps 
dt=median(diff(a(:,1)))*3e7;          %Width of each time step
dx=5e2/(nx-1);       %Width of space step
x=0:dx:1e4;            %Range of x (0,2) and specifying the grid points
u=zeros(nx,1);       un=zeros(nx,1);      %Preallocating un
vis=1e-6;            %Thermal diffusivity
beta=vis*dt/(dx*dx); %Stability criterion (0<=beta<=0.5, for explicit)
UL=1;  UR=1;                % Dirichlet B.C
UnL=1; UnR=-(Q_lith/k_cond);                     %Neumann B.C (du/dn=UnL) 

T = zeros(nt,nx);

%%
%Initial Conditions: Geotherm
u(1:end) = annualmeanT(1) + [1:nx].*dx*(Q_lith/k_cond);%dT/dz = Q/k

%%
%B.C vector
bc=zeros(nx-2,1);
bc(1)=vis*dt*UL/dx^2; bc(nx-2)=-UnR*vis*dt/dx;  %Mixed B.Cs
%Calculating the coefficient matrix for the implicit scheme
E=sparse(2:nx-2,1:nx-3,1,nx-2,nx-2);
A=E+E'-2*speye(nx-2);        %Dirichlet B.Cs
A(nx-2,nx-2)=-1;
%A(1,1)=-1; ; %Neumann B.Cs
D=speye(nx-2)-(vis*dt/dx^2)*A;

ULlist = interp1(a(:,1),annualmeanT,startofruninyears+[0:nt].*dt/3e7);

%%
%Calculating for each time step
i=2:nx-1;
T = zeros(nt+1,nx);
for it=0:nt
    if rem(it,1e5)==0;disp(it); end
    UL = ULlist(it+1);%Update UL
    bc(1)=vis*dt*UL/dx^2;
        
    u(1)=[];u(end)=[];
    U=u+bc;
    U=D\U;
    u=[UL;U;U(end)-UnR*dx];                      %Mixed b.c
    T(it+1,:) = u;
end