#### 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).

\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)

#### 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


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
allocate (tarray(tlength), stat=alloc_stat, errmsg=alloc_errmsg)
call errstop(alloc_stat, alloc_errmsg);


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

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
end do


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

do xi = 0, floor(xt / xint)
x = xs + xi*xint
end do
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.

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


Previous ChapterNext Chapter