Crank-Nicolson. Advection - Diffusion
I know nothing about your business logic but you're trying to transpose matlab in Python and it doesn't look like anything: there are some useless parameters, you are building arrays only to take their size. I'll try to provide advices to clean up the code and maybe others will provide runtime improvements.


Remove unused variables like F_bottom or tend which is only an alias to years.
Replaces calls like arr = np.linspace(begin, end, size); var = arr.size with var = size if you don't plan on using the array at all.
Use slice assignment to eliminate visual clutter.
Remove update_bc_dirichlet as you already perform the same kind of operation directly in the loop: stay consistent and see previous point.
Remove linalg_solver as it only indirects a single function call. If you want to avoid lookups in the loop, create an alias as local variable.
Some of the constants should either be extracted out of the function as real constants or defined as parameter with default values.
0 * np.ones(..) is np.zeros(..)
Use docstrings instead of comments to document your functions.


I would also turn adv_diff into a generator, but only because I don't know numpy enough to come up with a more idiomatic code:

import numpy as np 
import scipy as sp 
from scipy.sparse import spdiags 
from scipy import special 
import numba as nb 


@nb.jit 
def adv_diff(D, w, years, L=30, phi=1, dx=0.1, dt=0.001): 
    """What does this function do?"""
    N = int(L / dx) + 1
    C1 = np.zeros((N, 1)) 
    C1[0] = 1 

    linalg_solver = np.linalg.solve
    AL, AR = AL_AR_dirichlet(D, w, phi, dt, dx, N) 

    yield C1
    for _ in range(int(year/dt)): 
        C1 = linalg_solver(AL, AR.dot(C1)) 
        yield C1
        C1[:2, 0] = 1, 1


@nb.jit 
def AL_AR_dirichlet(D, w, phi, dt, dx, N): 
    """AL_AR_dirichlet: creates AL and AR matrices with Dirichlet BC""" 
    s = phi * D * dt / dx / dx # 
    q = phi * w * dt / dx # 
    e1 = np.ones((N, 1)) # 
    AL = spdiags(np.concatenate((e1 * (-s / 2 - q / 4), e1 * (1 + s), e1 * (-s / 2 + q / 4)), axis=1).T, [-1, 0, 1], N, N).toarray() 
    AR = spdiags(np.concatenate((e1 * (s / 2 + q / 4), e1 * (1 - s), e1 * (s / 2 - q / 4)), axis=1).T, [-1, 0, 1], N, N).toarray() 

    AL[0, :2] = 1, 0 
    AL[-1, -2:] = -s, 1 + s 
    AR[0, :2] = 1, 0 
    AR[-1, -2:] = s, 1 - s 

    return AL, AR 


if __name__ == '__main__': 
    D = 0.5 
    w = 3 
    t = 10 
    C = np.asarray(list(adv_diff(D, w, t)))