%* 1-D thermal diffusion model with gradient boundary condition
clear

%** define constant(s)
K=2.8e-7;                   % thermal diffusivity of dry sand, m^2/s 
Ggrad=25/1000;              % geothermal gradient K/m
day=24*60*60;               % seconds in one day
year=365;                   % days in one year 

%** model domain
zsurf=0;                    % surface coordinate, m
zbase=50;                   % base coordinate, m  "deep" 
dz=0.25;                    % vertical interval size, m
z=[zsurf:dz:zbase]';        % vertical coordinates, m
N=length(z)-1;              % intervals in model domain

ti=0;                       % initial time
tf=1;                       % number of years
nD=0.1;             % time step size in days

dt=nD*day;          % time step size in seconds
M=tf*year/nD;          % number of time steps, days


%* boundary conditions
%  start model at coldest part of temperature cycle, -1/2 pi
Tsurfmean=280;                % mean annual temperature, K
Tdiff=30;                     % annual temperature range, K
Tsurf=Tsurfmean + Tdiff*sin(linspace(-1/2*pi, (tf*2-1/2)*pi, M));  

%* variables
T=zeros(N+1,M);           %T history 

%* insert boundary and initial conditions into matrix T 
%  initial condition: use Tsurfmean & monotonic trend to Tbase 
T(1,1:M)=Tsurf;
T(N+1,1)=Tsurfmean+Ggrad*zbase;
T(:,1)=linspace(Tsurfmean,T(N+1,1), N+1)';  

%** numerics

%* set up A matrix (left-hand side of model equations)
A=zeros(N+1, N+1);  

% boundary value rows
A(1,1)= 1;

A(N+1,N)= -2*K*dt/dz^2;
A(N+1, N+1)= 1 + 2*K*dt/dz^2;

% interior rows
for i=2:N
    A(i,i-1)= -K*dt/dz^2;
    A(i,i)= 1 + 2*K*dt/dz^2;
    A(i,i+1)= -K*dt/dz^2;
end

%** let's go!

for n=2:M                  % time steps, n=1 is initial condition
    
    % set up b vector (right-hand side of model equations) using T of time n-1 
    b=zeros(N+1,1);
    b(1)=Tsurf(n);  
    b(N+1)=T(N+1,n-1) + 2*K*dt/dz * Ggrad;

    b(2:N)=T(2:N,n-1);
        
    % solve 
    T(:,n)=A\b;
end


%** plot result

%* a family of curves representing T(z) at selected time steps
Tmin=min(min(T));
Tmax=max(max(T));

figure(3)
clf
axis([Tmin Tmax -zbase zsurf])      % use -z to plot in a perspective that makes physical sense
hold on
plot(T(:,2:30/nD:M-1), -z, 'c-')       % plot every 30th day
plot(T(:,1), -z, 'b-')              % plot initial condition
xlabel('depth (m)')
ylabel('T (K)')
title('temperature cycle in a soil layer, implicit numerical scheme')


%* a plot of the complete T(z,t) space
figure(4)
surf(z,[ti:M-1]*nD,T')
shading interp
xlabel('depth (m)')
ylabel('day of year')
zlabel('T (K)')