Use a first order upwind (for the convection component) and a second order central difference (for the diffusion component). So the end result would be equivalent to discretising the equation,
$$
\frac{\partial u}{\partial t} = \frac{\partial \boldsymbol{v}}{\partial x} + D\frac{\partial^2 u}{\partial x^2}
$$
So using the $\theta$-method you will end up with,
$$
\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} =
\boldsymbol{v} \left[ \frac{1-\theta}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \right) + \frac{\theta}{2\Delta x} \left( u_{j+1}^{n+1} - u_{j-1}^{n+1} \right) \right] + \\ D \left[ \frac{1 - \theta}{(\Delta x)^2} \left( u_{j-1}^{n} - 2u_{j}^{n} + u_{j+1}^{n} \right) + \frac{\theta}{(\Delta x)^2} \left( u_{j-1}^{n+1} - 2u_{j}^{n+1} + u_{j+1}^{n+1} \right) \right]
$$
where $u^n$ and $u^{n+1}$ terms are the present and future time step, respectively. So a fully implicit scheme would be recovered by setting $\theta=1$ and the Crank-Nicolson by setting $\theta=1/2$.
To solve this equation you need to write as a linear system so group in terms of the solution variable $u$ and move the unknowns and known to different sides and you should be able to solve.
For example, the end result would have the form,
$$
\underbrace{\boldsymbol{A}}_{A}\cdot\underbrace{\boldsymbol{u^{n+1}}}_{x} = \underbrace{\boldsymbol{M}\cdot\boldsymbol{u^{n}}}_{d}
$$
where,
$$
\boldsymbol{A}=\begin{pmatrix}
1+2s\theta & -\theta(s + r) & & 0 \\
\theta(r-s) & 1+2s\theta & -\theta (s + r) & \\
& \ddots & \ddots & \ddots \\
& \theta(r-s) & 1+2s\theta & -\theta (s + r) \\
0 & & \theta(r-s) & 1+2s\theta \\
\end{pmatrix}
$$
and
$$
\boldsymbol{M}=\begin{pmatrix}
1-2s(1-\theta) & (1 - \theta)(s + r) & & 0 \\
(1 - \theta)(s - r) & 1-2s(1-\theta) & (1 - \theta)(s + r) & \\
& \ddots & \ddots & \ddots \\
& (1 - \theta)(s - r) & 1-2s(1-\theta) & (1 - \theta)(s + r) \\
0 & & (1 - \theta)(s - r) & 1-2s(1-\theta) \\
\end{pmatrix}
$$
Finally $\boldsymbol{u^n}$ and $\boldsymbol{u^{n+1}}$ are just column vectors of the solution variable discretised over the space. And very finally, $s=D\frac{\Delta t}{(\Delta x)^2}$ and $r=\boldsymbol{v}\frac{\Delta t}{2 \Delta x}$.
A side note.
For 1D conservation problems you could also consider the finite volume method, it is just as easy to implement. It has some advantages in the fact that the conservative property is preserved to the discretisation level. You have to be a little be careful because discretisation will be come unstable for advection dominated problems so you need to introduce adaptive-upwinding or exponential fitting. For a step-by-step guide see my notes.
