Analytical Advection-Dispersion

1. Introduction

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) = 00 \le x \le \inftyIC
c(0,t) = c_ot \ge 0BC
c(\infty,t) = 0t \ge 0BC

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)

2. Parameters

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

Used in sections 1:16 and 6

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:

{Input 2}
read *, co, v, d

Added to in section 1:5

Used in sections 1:16 and 6

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

Used in sections 1:16 and 6

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:

{Input 2} +=
read *, xs, xint, xt
read *, tlength
allocate (tarray(tlength), stat=alloc_stat, errmsg=alloc_errmsg)
call errstop(alloc_stat, alloc_errmsg);
read *, tarray

Added to in section 1:5

Used in sections 1:16 and 6

3. Concentration Function

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

{Concentration Function 3}
function c(x,t) result(c_r)
    implicit none
    real, intent(in) :: x, t
    real :: c_r

    if (x == 0) then
        c_r = co
    else
        c_r = co/2 * (erfc((x - v*t)/(2*sqrt(d*t))) &
                + exp(v*x/d)*erfc((x + v*t)/(2*sqrt(d*t))))
    end if
end function c

Used in section 6

4. Output

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

Used in sections 1:16 and 6

\pagebreak

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

5. Error Handling

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

{Stop On Error 5}
subroutine errstop(stat, errmsg)
    integer :: stat
    character(len=*) errmsg

    if (stat > 0) then
        print *, trim(errmsg)
        stop
    end if
end subroutine errstop

Used in section 6

6. Code

This is just the basic outline of the file.

{a_advect_dispers.f90 6}
program a_advect_dispers

    implicit none
    {Variable Declarations, 2}
    {Input, 2}
    {Output, 4}

contains

    {Record Output Function, 4}
    {Concentration Function, 3}
    {Stop On Error, 5}

endprogram

7. Plots


Previous ChapterNext Chapter