This is the analytical advection-dispersion model from last week. I've modified
it slightly to allow it to read parameters and also to make the code slightly
less terrible (though the main focus this week is on the *numerical* model).

The 1-D advection-dispersion equation is

\frac{\partial c}{\partial t} = D\frac{\partial^{2} c}{\partial x^{2}} - v\frac{\partial c}{\partial x}

where D is the dispersion coefficient, and v is the average linear flow velocity.

For the initial and boundary conditions

c(x,0) = 0 | 0 \le x \le \infty | IC |

c(0,t) = c_o | t \ge 0 | BC |

c(\infty,t) = 0 | t \ge 0 | BC |

the analytical solution is

c = \frac{c_o}{2}\left(\text{erfc}\left(\frac{x-vt}{2\sqrt{Dt}}\right)+e^{vx/D}\text{erfc}\left(\frac{x+vt}{2\sqrt{Dt}}\right)\right)

The variables that are not arguments to c are stored as module variables:

{Variable Declarations 2}

real :: co, v, d

Added to in sections 1:5, 1:6, 1:8, 1:9, 1:10, 1:11, 1:13 and 4

`co`

is a (unitless) concentration. `v`

is in meters per day. `d`

is in square
meters per day.

These are read in as a record from standard input in standard Fortran format:

The program also takes control parameters: start, interval, and total for x (meters) and list of t values (days):

{Variable Declarations 2} +=

real :: xs, xint, xt integer :: tlength real, dimension(:), allocatable :: tarray integer :: alloc_stat character(len=80) :: alloc_errmsg

Added to in sections 1:5, 1:6, 1:8, 1:9, 1:10, 1:11, 1:13 and 4

These are read as three records. The first reads the control parameters for x. The second is the number of t values that will be provided, and the third is the t values:

The implementation is fairly straightforward, basically just transcribing the initial condition and concentration formula.

The program first outputs a record with a zero and then reference x values.
It then outputs records with one t value and then c values. For details on
the output format, see

https://tachibanatech.com/litdoc/hw3.full.d/_book/a_heat_conduct.html#1:5

{Variable Declarations 2} +=

character(len=13), parameter :: real_fmtstr = "(sp,es16.7e3)" real :: x, t integer :: xi, ti

Added to in sections 1:5, 1:6, 1:8, 1:9, 1:10, 1:11, 1:13 and 2

This loop outputs the first record of reference x values:

{Output 4}

write (unit=*, fmt=real_fmtstr, advance="no") 0.0 do xi = 0, floor((xt + epsilon(xt)) / (xint - epsilon(xint))) x = xs + real(xi)*xint write (unit=*, fmt=real_fmtstr, advance="no") x end do write (unit=*, fmt="(a)", advance="yes") ""

Used in section 6

This loop outputs all the subsequent records:

{Output 4} +=

do ti = 1, size(tarray) t = tarray(ti) call at_time(t, xs, xint, xt) end do

Used in section 6

This is the subroutine that outputs a single record for a single point in time:

{Record Output Function 4}

subroutine at_time(t, xs, xint, xt) implicit none real :: t, xs, xint, xt write (unit=*, fmt=real_fmtstr, advance="no") t do xi = 0, floor(xt / xint) x = xs + xi*xint write (unit=*, fmt=real_fmtstr, advance="no") c(x,t) end do write (unit=*, fmt="(a)", advance="yes") "" end subroutine at_time

Used in section 6

This subroutine prints an error message and stops the program if an error has occurred.

This is just the basic outline of the file.

Previous ChapterNext Chapter