Annexe 3: Code

Cette annexe présente le code du modèle numérique TRANSCAP-1D ainsi que le cadre qui permet de réaliser les simulations Monte Carlo. Le modèle a été compilé avec le compilateur FORTRAN 90. Le transport de contaminant dans les sédiments est calculé avec la sousroutine "solution". Cette sousroutine est appelé par le cadre du programme "Monte Carlo", qui permet de réaliser des séries simulations stochastiques. Par la suite nous présentons le code de la sousroutine qui choisi les nombres aléatoires.

!

!***********************************************************************

!

! TRANSCAP-1D Finite volume solution of the 1D transport equation in bio-irrigated

! sediments

!

!***********************************************************************

!

! Arrays

!

! a : global left-hand-side matrix

! alfa: vector of irrigation coefficient (input)

! c : vector of nodal concentrations (ug/m**3)

! cfes : vector of concentration of FeS (umol/g)

! d1 : vector with the values of the left diagonal blocks of matrix A

! d2 : vector with the values of the middle diagonal blocks of matrix A

! d3 : vector with the values of the right diagonal blocks of matrix A

! f : global vector of known values (right hand side)

! port: vector of tube porosity

! r : right hand side matrix

! s : right hand side vector (Stu)

! twant : vector of desired time(s) for output

! x : vector of nodal coordinates (m)

!

! Variables

!

! Real

!

! alfa1, alfa2: irrigation coefficients (1/day & 1/m)

! alfa = alfa1*exp(alfa2*(x-xmax))

! area : column area (m**2)

! biox : maximum depth of bioturbation (m) (has to be a multiple of dx!)

! k1, k2 : coefficients of linear dissolution equation

! conc : concentration at the water/sediment interface (ug/m**3)

! cour : courant number

! diff : effective diffusion coefficient (m*2/day)

! ds : diffusion/dispersion coefficient in sediment (m*2/day) d = v*disp+diff

! dm : molecular diffusion

! disp : dispersivity (m)

! dt : time increment (day)

! dx : distance increment (m) (homogeneous grid!)

! epsil : time weighting factor (normally =1)

! fluxlo : flux through lower boundary

! fluxup : flux through upper boundary

! gamma : cinetic coefficient for dissolution reaction

! mins : mass of dissolved contaminant before simulation in sediment

! mint : mass of dissolved contaminant before simulation in tube

! mdiss : mass entering the system by dissolution (tube)

! mends : mass inside the column after simulation in sediment

! mendt : mass inside the column after simulation in tube

! mex: mass exchanged between sediment and tube

! mouts : mass leaving from the boundaries (sediment)

! moutt : mass leaving from the boundaries (tube)

! pec : peclet number

! pors : sediment porosity

! porscap : cap sediment porosity

! port0 : tube porosity at the surface of the sediment

! port=port0*exp(alfa2*(x-xmax))

! qs : advective flow in sediment (m**3/day)

! qq : flux contribution for the inflow boundary equation

! retard : retardation coefficient retard = 1+rho*xkd/pors

! tmax : maximum simulation time (day)

! vs : average water velocity in sediment (m/day) vs = qs/(area*pors)

! vt : average water velocity in tube (m/day)

! xmax : lenght of column (m)

!

!

! Integer

!

! nn : number of nodes

! neq : number of equations per node (tube + sediment)

! nunk : number of unknowns

! ntwant : number of desired time output

!

!

!**********************************************************************

!

! ...define the types of the variables

!

subroutine solution(x, cfes, cs, ct, twant)

!

USE PORTLIB

!

common alfa1, alfa2, area, biox, capdx, &

k1, k2, conc, dm, disp, diff, dp, dt, dx, epsil, &

gamma, pors, porscap, port0, qs, retard, &

tmax, vt, xmax

double precision alfa1, alfa2, area, biox, capdx, &

k1, k2, conc, dm, disp, diff, dp, dt, dx, epsil, &

gamma, pors, porscap, port0, qs, retard, &

tmax, vt, xmax

common sim, nsim, ntwant, nn, AllocateStatus

integer sim, nsim, ntwant, nn, AllocateStatus

common title

character *80 title

!

double precision b, balsed, baltub, cours, courscap, court, &

csurf, ds, dscap, dscheck, dscapcheck, dxmax, &

fluxlos, fluxlot, fluxups, fluxupt, mdiss, mends, mendt, &

mex, mins, mint, mouts, moutt, pecs, pecscap, pect, qq, &

sumex, sumlos, sumlot, sumup, sumups, sumupt, sumdiss, &

tcheck, ttime, vs, vscap, vscheck, vscapcheck

double precision x(nn), cs(nn), ct(nn), cfes(nn),twant(ntwant)

double precision, dimension(:), allocatable :: alfa, c, &

f, port, s, d1, d2, d3

double precision, dimension (:,:), allocatable :: a, r

integer cap, i, ii, j, jj, n, neq, nprint, nunk, nunk2, ucap

parameter (neq=2)

character *1 tab

character(8) whatimeisit

tab=char(9)

!

! ...open the working files

!

open(4, file='mcflux.out')

open(5, file='mcoutput.out')

open(6, file='mccsurf.out')

!

! ...calculate total number of unknown

!

nunk = neq*nn

nunk2 = neq*nunk

!

! ...allocate the memory for the arrays

!

if (sim.eq.1) then

allocate(c(nunk), alfa(nn), port(nn), f(nunk), s(nunk), &

d1(nunk2), d2(nunk2), d3(nunk2), stat = AllocateStatus)

if (AllocateStatus /= 0) stop "***Not enough memory***"

allocate(a(nunk,5), r(nunk,5),stat = AllocateStatus)

if (AllocateStatus /= 0) stop "***Not enough memory***"

endif

!

! ...find the maximum dx to use for the Peclet and Courant

! numbers

!

dxmax = -1.0d30

if(dx.gt.dxmax) dxmax = dx

!

! ...calculate number of nodes corresponding to the capping layer (cap)

! and number of nodes under the capping layer (ucap)

!

cap=int(capdx/dx)

ucap=nn-cap

!

! ...echo back output

!

write(5,405) sim,nn,xmax,area,tmax,dt,epsil,dp,biox,capdx

405 format(//,60('*'),/,5x,'Simulation', t20 ,i4,/,60('*'),//,5x, &

'Number of nodes:',t55,i5,/,5x, &

'Length of the column (m):',t48,d12.5, &

/,5x,'Area of the column (m**2):',t48,d12.5, &

/,5x,'Maximum simulation time (day):',t48,d12.5, &

/,5x,'Time step (day):', t48, d12.5, &

/,5x,'Time weighting factor:',t48,d12.5, &

/,5x,'(DP) Depth of bioirrigation (m):',t48,d12.5, &

/,5x,'(INT) Depth of bioirrigation (m):',t48,d12.5, &

/,5x,'Depth of the capping layer (m):',t48,d12.5)

write(5,406) qs, vt, pors, porscap, port0, dm, diff, disp

406 format(5x, 'Advective flow in sediment(m**3/day):',t48,d12.5, &

/,5x,'Irrigation velocity (m/day):',t48,d12.5, &

/,5x,'Sediment porosity:',t48,d12.5, &

/,5x,'Capping layer porosity:',t48,d12.5, &

/,5x,'Tube porosity at surface:',t48,d12.5, &

/,5x,'Molecular diffusion (m*2/day):',t48,d12.5, &

/,5x,'Effective diffusion coefficient (m*2/day):',t48,d12.5, &

/,5x,'Dispersivity (m):',t48,d12.5)

write(5,407) k1, k2, gamma, alfa1, alfa2, conc

407 format(5x,'Lin. dissolution eq. coeff. b1 & b2:',t48, &

d12.5,3x,d12.5, &

/,5x,'Cinetic coeff for dissolution (1/day):',t48,d12.5, &

/,5x,'Irrigation coefficient alfa1 (1/day):' ,t48,d12.5, &

/,5x,'Irrigation coefficient alfa2 (1/m):' ,t48,d12.5, &

//,5x,'Conc. at wat./sed. interface (ug/m**3) = ',t48,d12.5)

!

! ...calculate alfa coefficients and tube porosity at each node

!

do 16 i=1,nn

alfa(i)=0.0d0

port(i)=0.0d0

b = xmax - x(i)

if(b.le.biox) then

alfa(i) = alfa1*DEXP(alfa2*(x(i)-xmax))

port(i) = port0*DEXP(alfa2*(x(i)-xmax))

endif

16 continue

!

! ...Calculate mass inside the column (sediment and tubes)

! at the beginning of the simulation

!

mins=0.0d0

mint=0.0d0

do 8 i=1,nn

nc=ucap-i

if (nc.ge.0) then

mins=mins+(cs(i)*pors)*area*dx

elseif (nc.lt.0) then

mins=mins+(cs(i)*porscap)*area*dx

endif

mint=mint+(ct(i)*port(i))*area*dx

8 continue

!

! ...assemble concentration vector c

!

do 12 i=1,nn

j= 2*i-1

n= 2*i

c(j)=cs(i)

c(n)=ct(i)

12 continue

!

! ...compute constants

!

vs = qs/(area*pors)

vscap = qs/(area*porscap)

ds = vs*disp + diff

dscap = vscap*disp + diff

!

! ...echo back calculated parameters

!

write(5,409) vs, vscap, ds, dscap, retard

409 format(/,5x, 'In sediment:',/, &

5x,'Unretarded average linear vel. (in sed.) (m/day):',t48,d12.5,/, &

5x,'Unretarded average linear vel. (in cap) (m/day):',t48,d12.5,/, &

5x,'Unretarded diff./disp. coeff. (in sed.) (m2/day):',t48,d12.5,/, &

5x,'Unretarded diff./disp. coeff. (in cap) (m2/day):',t48,d12.5,/, &

5x,'Retardation coefficient:',t48, d12.5)

!

! ...calculate retarded constants for the sediment

! to check and print Peclet and Courant

!

if(vs.eq.0.0) then

write(5,410)

410 format (/,5x, 'No advection in sediment')

else

vscheck=vs/retard

vscapcheck=vscap/retard

dscheck=ds/retard

dscapcheck=dscap/retard

pecs=vscheck*dxmax/dscheck

pecscap=vscapcheck*dxmax/dscapcheck

cours=vscheck*dt/dxmax

courscap=vscapcheck*dt/dxmax

write(5,411) pecs, pecscap, cours, courscap

411 format(/,5x,'Maximum Peclet number:',t48,d12.5, &

/,5x,'Maximum Peclet number in cap:',t48,d12.5, &

/,5x, 'Maximum Courant number in sediment:', t48,d12.5 &

/,5x, 'Maximum Courant number in cap:', t48,d12.5)

endif

!

! ...calculate retarded constants for the tube

! to check and print Peclet and Courant

!

if(vt.eq.0.0) then

write(5,412)

412 format (/,5x, 'No advection in tube')

else

pect=vt*dxmax/dm

court=vt*dt/dxmax

write(5,413)pect,court

413 format(/,5x,'Maximum Peclet number for tube:',t48,d12.5, &

/,5x, 'Maximum Courant number for tube:', t48,d12.5)

endif

!

! ...assemble the first row of the A matrix (sediment)

!

if (ucap .gt. 1) then

a(1,1)= 0.0d0

a(1,2)= 0.0d0

a(1,3)= pors*dx/dt+epsil*(pors*(vs/retard)+alfa(1)*dx)

a(1,4)= -epsil*alfa(1)*dx

a(1,5)= 0.0d0

elseif (ucap .eq. 1) then

a(1,1)= 0.0d0

a(1,2)= 0.0d0

a(1,3)= pors*dx/dt+epsil*(pors*(vs+vscap)/2/retard+alfa(1)*dx)

a(1,4)= -epsil*alfa(1)*dx

a(1,5)= 0.0d0

elseif (ucap.eq.0) then

a(1,1)= 0.0d0

a(1,2)= 0.0d0

a(1,3)= porscap*dx/dt+epsil*(porscap*(vscap/retard)+alfa(1)*dx)

a(1,4)= -epsil*alfa(1)*dx

a(1,5)= 0.0d0

endif

!

! ...assemble the second row of the A matrix (tube)

!

if (alfa(1) .eq. 0.0d0) then

a(2,1) = 0.0d0

a(2,2) = 0.0d0

a(2,3) = 1.0d0

a(2,4) = 0.0d0

a(2,5) = 0.0d0

else

a(2,1)= 0.0d0

a(2,2)= -epsil*alfa(1)*dx

a(2,3)= port(1)*dx/dt+epsil*(port(1)*vt+(alfa(1)*dx))+gamma*dx

a(2,4)= 0.0d0

a(2,5)= 0.0d0

endif

!

! ...assemble middle nodes of A matrix (sediment)

!

if (pecs.gt.2.0) then

!

! ...if Peclet number is greater than 2 then use upwind weighting

!

do 52 i=3,nunk-3,2

n=(i+1)/2

nc=ucap-n

if (nc.gt.0) then

a(i,1)= -epsil*pors*(ds/retard/dx+vs/retard)

a(i,2)= 0.0d0

a(i,3)= pors*dx/dt+epsil*(pors*2*ds/retard/dx+pors*vs/retard+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*pors*(ds/retard/dx)

elseif (nc.eq.0) then

a(i,1)= -epsil*pors*(ds/retard/dx+vs/retard)

a(i,2)= 0.0d0

a(i,3)= pors*dx/dt+epsil*(pors*(ds+(ds+dscap)/2)/

retard/dx+pors* (vs+vscap)/2/retard+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*porscap*((ds+dscap)/2/retard/dx)

elseif (nc.eq.-1) then

a(i,1)= -epsil*pors*((ds+dscap)/2/retard/dx+(vs+vscap)/2/retard)

a(i,2)= 0.0d0

a(i,3)= porscap*dx/dt+epsil*(porscap*(dscap+(ds+dscap)/2) /

retard/dx+porscap*vscap/retard+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*porscap*(dscap/retard/dx)

elseif (nc.lt.-1) then

a(i,1)= -epsil*porscap*(dscap/retard/dx+vscap/retard)

a(i,2)= 0.0d0

a(i,3)= porscap*dx/dt+epsil*

(porscap*2*dscap/retard/dx+porscap*vscap/retard+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*porscap*(dscap/retard/dx)

endif

52 continue

else

!

! ...else use central weighting

!

do 54 i=3,nunk-3,2

n=(i+1)/2

nc=ucap-n

if (nc.gt.0) then

a(i,1)= -epsil*pors*(ds/retard/dx+vs/retard/2)

a(i,2)= 0.0d0

a(i,3)= pors*dx/dt+epsil*(pors*2*ds/retard/dx+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*pors*(ds/retard/dx-vs/retard/2)

elseif (nc.eq.0) then

a(i,1)= -epsil*pors*(ds/retard/dx+vs/retard/2)

a(i,2)= 0.0d0

a(i,3)= pors*dx/dt+epsil*(pors*(ds+(ds+dscap)/2)/

retard/dx+pors*((vs+vscap)/2-vs)/2/retard+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*porscap*((ds+dscap)/2/retard/dx-(vs+vscap)/2/retard/2)

elseif (nc.eq.-1) then

a(i,1)= -epsil*pors*((ds+dscap)/2/retard/dx+(vs+vscap)/2/retard/2)

a(i,2)= 0.0d0

a(i,3)= porscap*dx/dt+epsil*(porscap*(dscap+(ds+dscap)/2)/

retard/dx+porscap*(vscap-(vs+vscap)/2)/2/retard+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*porscap*(dscap/retard/dx-vscap/retard/2)

elseif (nc.lt.-1) then

a(i,1)= -epsil*porscap*(dscap/retard/dx+vscap/retard/2)

a(i,2)= 0.0d0

a(i,3)= porscap*dx/dt+epsil*(porscap*2*dscap/retard/dx+alfa(n)*dx)

a(i,4)= -epsil*alfa(n)*dx

a(i,5)= -epsil*porscap*(dscap/retard/dx-vscap/retard/2)

endif

54 continue

endif

!

! ...assemble middle rows of A matrix (tube)

!

do 56 i=4,nunk-2,2

n=i/2

if (alfa(n) .eq. 0.0d0) then

a(i,1) = 0.0d0

a(i,2) = 0.0d0

a(i,3) = 1.0d0

a(i,4) = 0.0d0

a(i,5) = 0.0d0

else

a(i,1) = -epsil*port(n-1)*(dm/dx+vt)

a(i,2) = -epsil*(alfa(n)*dx)

a(i,3) = port(n)*dx/dt+epsil*

(dm*2*port(n)/dx+port(n)*vt+alfa(n)*dx)+gamma*dx

a(i,4) = 0.0d0

a(i,5) = -epsil*port(n+1)*(dm/dx)

endif

56 continue

!

! ...set up first-type boundary condition at outlow

! Trick for Dirichlet:

! set the diagonal term of the L.H.S coefficient matrix to a very

! high number and set the flux contribution (qq) in order tu ensure

! that c(nn) goes to the desired value

!

a(nunk-1,1)=0.0d0

a(nunk-1,2)=0.0d0

a(nunk-1,3)=1.0d35

a(nunk-1,4)=0.0d0

a(nunk-1,5)=0.0d0

a(nunk,1)=0.0d0

a(nunk,2)=0.0d0

a(nunk,3)=1.0d35

a(nunk,4)=0.0d0

a(nunk,5)=0.0d0

qq=1.0d35*conc

!

! ...Assemble first row of R matrix (sediment)

!

if (ucap .gt. 1) then

r(1,1)= 0.0d0

r(1,2)= 0.0d0

r(1,3)= pors*dx/dt-(1-epsil)*(pors*vs/retard+alfa(1)*dx)

r(1,4)= (1-epsil)*alfa(1)*dx

r(1,5)= 0.0d0

elseif (ucap.eq.1) then

r(1,1)= 0.0d0

r(1,2)= 0.0d0

r(1,3)= pors*dx/dt-(1-epsil)*(pors*(vs+vscap)/2/retard+alfa(1)*dx)

r(1,4)= (1-epsil)*alfa(1)*dx

r(1,5)= 0.0d0

elseif (ucap.eq.0) then

r(1,1)= 0.0d0

r(1,2)= 0.0d0

r(1,3)= porscap*dx/dt-(1-epsil)*(porscap*vscap/retard+alfa(1)*dx)

r(1,4)= (1-epsil)*alfa(1)*dx

r(1,5)= 0.0d0

endif

!

! ...assemble the second row of the R matrix (tube)

!

r(2,1)= 0.0d0

r(2,2)= (1-epsil)*alfa(1)*dx

r(2,3)= port(1)*dx/dt-(1-epsil)*(port(1)*vt+alfa(1)*dx)

r(2,4)= 0.0d0

r(2,5)= 0.0d0

!

! ...assemble middle rows of R matrix (sediment)

!

if (pecs.gt.2.0) then

!

! ...upwind weighting

!

do 60 i=3,nunk-3,2

n=(i+1)/2

nc=ucap-n

if (nc.gt.0) then

r(i,1)= (1-epsil)*pors*(ds/retard/dx+vs/retard)

r(i,2)= 0.0d0

r(i,3)= pors*dx/dt-(1-epsil)*(pors*(2*ds/retard/dx+vs/retard)+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*pors*(ds/retard/dx)

elseif (nc.eq.0) then

r(i,1)= (1-epsil)*pors*(ds/retard/dx+vs/retard)

r(i,2)= 0.0d0

r(i,3)= pors*dx/dt-(1-epsil)*(pors*((ds+(ds+dscap)/2)/

retard/dx+(vs+vscap)/2/retard)+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*porscap*((ds+dscap)/2/retard/dx)

elseif (nc.eq.-1) then

r(i,1)= (1-epsil)*pors*((ds+dscap)/2/retard/dx+(vs+vscap)/2/retard)

r(i,2)= 0.0d0

r(i,3)= porscap*dx/dt-(1-epsil)*(porscap*((dscap+(ds+dscap)/2)/

retard/dx+vscap/retard)+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*porscap*(dscap/retard/dx)

elseif (nc.lt.-1) then

r(i,1)= (1-epsil)*porscap*(dscap/retard/dx+vscap/retard)

r(i,2)= 0.0d0

r(i,3)= porscap*dx/dt-(1-epsil)*(porscap*(2*dscap/retard/dx+vscap/retard)

+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*porscap*(dscap/retard/dx)

endif

60 continue

else

!

! ...central weighting

!

do 72 i=3,nunk-3,2

n=(i+1)/2

nc=ucap-n

if (nc.gt.0) then

r(i,1)= (1-epsil)*pors*(ds/retard/dx+vs/retard/2)

r(i,2)= 0.0d0

r(i,3)= pors*dx/dt-(1-epsil)*((pors*2*ds/retard/dx)+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*pors*(ds/retard/dx-vs/retard/2)

elseif (nc.eq.0) then

r(i,1)= (1-epsil)*pors*(ds/retard/dx+vs/retard/2)

r(i,2)= 0.0d0

r(i,3)= pors*dx/dt-(1-epsil)*((pors*(ds+(ds+dscap)/2)/retard/dx)+

pors*((vs+vscap)/2-vs)/2/retard+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*porscap*((ds+dscap)/2/retard/dx-(vs+vscap)/2/retard/2)

elseif (nc.eq.-1) then

r(i,1)= (1-epsil)*pors*((ds+dscap)/2/retard/dx+(vs+vscap)/2/retard/2)

r(i,2)= 0.0d0

r(i,3)= porscap*dx/dt-(1-epsil)*((porscap*(dscap+(ds+dscap)/2)/retard/dx)

+porscap*(vscap-(vs+vscap)/2)/2/retard+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*porscap*(dscap/retard/dx-vscap/retard/2)

elseif (nc.lt.-1) then

r(i,1)= (1-epsil)*porscap*(dscap/retard/dx+vscap/retard/2)

r(i,2)= 0.0d0

r(i,3)= porscap*dx/dt-(1-epsil)*((porscap*2*dscap/retard/dx)+alfa(n)*dx)

r(i,4)= (1-epsil)*alfa(n)*dx

r(i,5)= (1-epsil)*porscap*(dscap/retard/dx-vscap/retard/2)

endif

72 continue

endif

!

! ...assemble middle rows of R matrix (tube)

!

do 62 i=4,nunk-2,2

n=i/2

r(i,1)= (1-epsil)*port(n-1)*(dm/dx+vt)

r(i,2)= (1-epsil)*(alfa(n)*dx)

r(i,3)= port(n)*dx/dt-(1-epsil)*(2*dm*port(n)/dx+port(n)*vt+alfa(n)*dx)

r(i,4)= 0.0d0

r(i,5)= (1-epsil)*port(n+1)*(dm/dx)

62 continue

!

! ...assemble last rows of R matrix for Dirichlet condition

!

do 64 j=1,5

r(nunk-1,j)=0.0d0

r(nunk,j)=0.0d0

64 continue

!

! ...Assemble s vector

!

do 66 i=1,nunk-1,2

n=(i+1)/2

b = xmax - x(n)

s(i)=0.0d0

s(i+1)=0.0d0

if(b.le.biox) then

s(i+1)=gamma*1.0d3*(k1*cfes(n)+k2)*dx

endif

66 continue

!

! ...Decompose A matrix in 3 vectors of lenght neq*nn

!

!

do 68 ii=1,nn

i=ii*4-3

j=ii*2-1

d1 (i) = a(j,1)

d1 (i+1) = a(j,2)

d1 (i+2) = 0.0d0

d1 (i+3) = a(j+1,1)

d2 (i) = a(j,3)

d2 (i+1)= a(j,4)

d2 (i+2)= a(j+1,2)

d2 (i+3)= a(j+1,3)

d3 (i) = a(j,5)

d3 (i+1) = 0.0d0

d3 (i+2) = a(j+1,4)

d3 (i+3) = a(j+1,5)

68 continue

!

! ...start loop over time

!

ttime=0.0d0

sumups=0.0d0

sumlos=0.0d0

sumupt=0.0d0

sumlot=0.0d0

mdiss=0.0d0

mex=0.0d0

500 ttime=ttime+dt

tcheck=ttime-tmax

!

! ...if time > maximum time then stop

!

if(tcheck.gt.1.0d-3) goto 999

!

! ...multiply out right-hand side to form the vector of knowns

!

! ...First node

!

f(1)=r(1,3)*c(1)+r(1,4)*c(2)+r(1,5)*c(3)+s(1)

if (alfa(1) .eq. 0.0d0) then

f(2)=0.0d0

else

f(2)=r(2,2)*c(1)+r(2,3)*c(2)+r(2,5)*c(4)+s(2)

endif

!

! Middle nodes

!

do 70 i=3,nunk-3,2

n=(i+1)/2

f(i)=r(i,1)*c(i-2)+r(i,2)*c(i-1)+r(i,3)*c(i)+ &

r(i,4)*c(i+1)+r(i,5)*c(i+2)+s(i)

if (alfa(n) .eq. 0.0d0) then

f(i+1) = 0.0d0

else

f(i+1)=r(i+1,1)*c(i-1)+r(i+1,2)*c(i)+r(i+1,3)*c(i+1)+ &

r(i+1,4)*c(i+2)+r(i+1,5)*c(i+3)+s(i+1)

endif

70 continue

f(nunk-1)=qq

f(nunk)=qq

!

! ...solve the matrix

!

call Thomas(d2,d3,d1,f,c,neq,nn,nunk,nunk2)

!

! ...Calculate contaminant flux through upper boundary of sediment

! depending on the Peclet number

!

if (pecs.gt.2.0) then

if (cap.eq.0) then

fluxups=pors*(vs/retard*(c(nunk-3))-ds/retard*(c(nunk-1)-c(nunk-3))/dx)

elseif (cap.eq.1) then

fluxups=pors*(vs+vscap)/2/retard*(c(nunk-3))-1/retard*

(dscap*porscap*c(nunk-1)-(ds+dscap)/2*pors*c(nunk-3))/dx

elseif (cap.gt.1) then

fluxups=porscap*(vscap/retard*(c(nunk-3))-dscap/retard*

(c(nunk-1)-c(nunk-3))/dx)

endif

else

if (cap.eq.0) then

fluxups=pors*(vs/retard*(c(nunk-1)+c(nunk-3))/2-ds/retard*(c(nunk-1)-

c(nunk-3))/dx)

elseif (cap.eq.1) then

fluxups=1/retard*(vscap*porscap*c(nunk-1)+(vs+vscap)/2*pors*

c(nunk-3))/2-1/retard*(dscap*porscap*c(nunk-1)-(ds+dscap)/2*pors*c(nunk-3))/dx

elseif (cap.gt.1) then

fluxups=porscap*(vscap/retard*(c(nunk-1)+c(nunk-3))/2-

dscap/retard*(c(nunk-1)-c(nunk-3))/dx)

endif

endif

!

! ...Contaminant flux through lower boundary of sediment

!

fluxlos=pors*(-ds/retard*(c(3)-c(1))/dx)

!

! ...Contaminant flux through upper boundary of tube

!

ii= (nunk-2)/2

jj= (nunk)/2

fluxupt=((vt)*port(ii)*c(nunk-2)-dm*(port(jj)*c(nunk)-port(ii)*c(nunk-2))/dx)

!

! ...Contaminant flux through lower boundary (end of bioturbation zone)

!

b=xmax-biox

i=nint(b*2/dx)

j=i+2

ii=i/2

jj=j/2

fluxlot=(-dm*(port(jj)*c(j)-port(ii)*c(i))/dx)

!

! ...Sum of the mass flux to the upper and lower boudaries of sediment

!

sumups=sumups+dabs(fluxups)*area*dt

sumlos=sumlos+dabs(fluxlos)*area*dt

!

! ...Sum of the mass flux to the upper and lower boudaries of tube

!

sumupt=sumupt+dabs(fluxupt)*area*dt

sumlot=sumlot+dabs(fluxlot)*area*dt

!

! ...Contaminant flux exchanged between tube and sediment

!

sumex = 0.0d0

do 76 i=1,nn-1

j=2*i-1

n=2*i

b = xmax - x(i)

if(b.le.biox) then

sumex=sumex+alfa(i)*(c(j)-c(n))*area*dx*dt

endif

76 continue

mex = mex+sumex

!

! ...Contaminant mass "created" by source in the bioturbated zone

!

sumdiss = 0.0d0

do 78 i=1,nn-1

n=2*i

b = xmax - x(i)

if(b.le.biox) then

sumdiss=sumdiss+gamma*(1.0d3*(k1*cfes(i)+k2)-c(n))*dx*area*dt

endif

78 continue

mdiss = mdiss+sumdiss

!

! ...check if ttime is desired for output

!

nprint=0

do 80 i=1,ntwant

tcheck=twant(i)-ttime

tcheck=dabs(tcheck)

if (tcheck.lt.1.0d-4) nprint=1

80 continue

!

! ...print the dephts at which the concentrations are evaluated

!

if(nprint.eq.1) then

write(5,419)

419 format(/,24x,'X(m):',4x,\)

do 82 i=1,nn

write(5,420) x(i)

82 continue

420 format(4x,d12.5,\)

!

! ... print the time and the concentrations in the sediment

!

write(5,434) ttime

434 format(/,'Time(day):',d8.3,6x,'Cs:',6x,\)

do 84 i=1,nunk-1,2

write(5,435) c(i)

84 continue

435 format(4x,d12.5,\)

!

! ...print the concentrations in the tube

!

write (5,436)

436 format (/,24x,'Ct:',6x,\)

do 86 i=2,nunk,2

write(5,437) c(i)

86 continue

437 format(4x,d12.5,\)

end if

!

! ...Go back to beginning of time loop

!

goto 500

!

! ...Print mass fluxes for sediment and tubes

!

999 write(5,440) sumups, sumlos, sumupt, sumlot

440 format(//,5x,'Mass leaving from the upper boundary of the sed.(ug):',t70,d12.5, &

/,5x,'Mass leaving from the lower boundary of the sed.(ug):',t70,d12.5, &

/,5x,'Mass leaving from the upper boundary of the tube (ug):',t70,d12.5, &

/,5x,'Mass leaving from the lower boundary of the tube (ug):',t70,d12.5)

!

! ...Print mass flux that left from the upper bounadry at the end of the simulation

!

sumup=sumups+sumupt

write(4,446) sim, sumup

446 format(//,5x,'Simulation',t20 ,i4, &

/,5x,'Mass leaving from the upper boundary (ug):',t70,d12.5)

!

! ...print surface concentration at the end of the simulation

!

csurf=c(nunk-3)

write(6,439)sim, csurf

439 format(//,5x,'Simulation',t20,i4, &

/,5x,'Surface concentration in the sediment (ug/m3):',t70,d12.5)

!

! Mass balance

!

! ...Mass leaving the sediment by the boundaries

!

mouts=sumlos+sumups

!

! ...Mass leaving the tubes by the boundaries

!

moutt=sumlot+sumupt

!

! ...Mass inside the column after simulation

!

mends=0.0d0

mendt=0.0d0

do 88 i=1, nunk-1,2

n=(i+1)/2

nc=ucap-n

if (nc.ge.0.0d0) then

mends=mends+c(i)*pors*area*dx

else

mends=mends+c(i)*porscap*area*dx

endif

88 continue

do 90 i=2, nunk,2

j=i/2

mendt=mendt+c(i)*port(j)*area*dx

90 continue

balsed=mins-mex-mouts-mends

baltub=mint+mdiss+mex-moutt-mendt

!

! ...Print the mass balance for sediment

!

write(5,442) mins, mends, mouts, mex, balsed

442 format(//,5x,'Diss mass of contaminant in sed. before simulation (ug):',t70,d12.5, &

/,5x,'Diss mass inside the sediment after simulation (ug):',t70,d12.5, &

/,5x,'Diss mass which crossed the sediment boundaries (ug):',t70,d12.5, &

/,5x,'Diss mass exchanged with tube (ug):',t70,d12.5, &

/,5x,'Balance for sediment:',t70,d12.5)

!

! ...Print the mass balance for sediment

!

write(5,444) mint, mendt, moutt, mex, mdiss, baltub

444 format(//,5x,'Diss mass of contaminant in tubes before simulation (ug):',t70,d12.5, &

/,5x,'Diss mass inside the tubes after simulation (ug):',t70,d12.5, &

/,5x,'Diss mass which crossed the tubes boundaries (ug):',t70,d12.5, &

/,5x,'Diss mass exchanged with sediment (ug):',t70,d12.5, &

/,5x,'Mass dissolved in the tube through oxigenation:',t70,d12.5, &

/,5x,'Balance for tube:',t70,d12.5)

!

write(*,*) ' '

whatimeisit = CLOCK ()

write(*,448) sim

448 format(1x,'Simulation ',i5,' endet at ')

write(*,*) whatimeisit

write(*,*) ' '

return

end

!

!**********************************************************************

!

! ...SUBROUTINE Thomas, solves the block tridiagonal matrix.

!

!**********************************************************************

SUBROUTINE Thomas(a,b,c,d,unew,neq,nn,nunk,nunk2)

implicit none

double precision sum

double precision a, b, c, d, unew, w, qq, wtemp, g, z

integer i, ii, ia, ic, icc, icc2, ivar1, ivar2, j, k, kcc, &

kk, neq2, neq, nn, nunk, nunk2

parameter (neq2=4)

dimension a(nunk2), b(nunk2), c(nunk2), d(nunk), unew(nunk), &

w(nunk2), qq(nunk2), wtemp(neq2), g(nunk), z(neq2)

!

! ...First, perform decomposition

!

DO 11 i=1,neq2

w(i)=a(i)

11 CONTINUE

!

! ...Invert w

!

DO 31 ic=2,nn

icc=(ic-2)*neq2

CALL Invert(w,wtemp,neq,icc,nunk2)

!

DO 27 i=1,neq

ii=(i-1)*neq

DO 25 j=1,neq

sum=0.0d0

DO 23 k=1,neq

kk=(k-1)*neq

sum=sum+ wtemp(ii+k)*b(icc+j+kk)

23 CONTINUE

qq(icc+ii+j)=sum

25 CONTINUE

27 CONTINUE

!

DO 37 i=1,neq

ii=(i-1)*neq

DO 35 j=1,neq

sum=0.0d0

DO 33 k=1,neq

kk=(k-1)*neq

sum=sum+ c(icc+neq2+ii+k)*qq(icc+j+kk)

33 CONTINUE

w(icc+neq2+ii+j)=a(icc+neq2+ii+j)-sum

35 CONTINUE

37 CONTINUE

31 CONTINUE

!

! ...Perform forward substitution

!

! ...Invert the first w block.

!

icc=0

CALL Invert(w,wtemp,neq,icc,nunk2)

DO 45 i=1,neq

ii=(i-1)*neq

sum=0.0d0

DO 43 k=1,neq

sum=sum+ wtemp(ii+k)*d(k)

43 CONTINUE

g(i)=sum

45 CONTINUE

DO 131 ic=2,nn

icc=(ic-2)*neq

kcc=(ic-1)*neq2

!

! ...Invert w

!

CALL Invert(w,wtemp,neq,kcc,nunk2)

!

DO 55 i=1,neq

ii=(i-1)*neq

sum=0.0d0

DO 53 k=1,neq

sum=sum+ c(kcc+ii+k)*g(icc+k)

53 CONTINUE

z(i)=d(icc+neq+i)-sum

55 CONTINUE

DO 65 i=1,neq

ii=(i-1)*neq

sum=0.0d0

DO 63 k=1,neq

sum=sum+ wtemp(ii+k)*z(k)

63 CONTINUE

g(icc+neq+i)=sum

65 CONTINUE

131 CONTINUE

!

! ...Perform backward substitution

!

ia=(nn-1)*neq+1

DO 71 i=ia,nn*neq

unew(i)=g(i)

71 CONTINUE

ivar1=(nn-1)*neq

ivar2=(nn-1)*neq2

DO 231 ic=1,nn-1

icc=ivar1-ic*neq

icc2=ivar2-ic*neq2

DO 95 i=1,neq

ii=(i-1)*neq

sum=0.0d0

DO 93 k=1,neq

sum=sum+ qq(icc2+ii+k)*unew(icc+neq+k)

93 CONTINUE

unew(icc+i)=g(icc+i)-sum

95 CONTINUE

231 CONTINUE

RETURN

END

!**********************************************************************

!

! ...SUBROUTINE for matrix inversion

!

!**********************************************************************

SUBROUTINE Invert(w,ainv,neq,icc,nunk2)

implicit none

double precision aa,ainv, w, u, vi, y, l

double precision sum

integer i, ic, ii, in, in1, j, jj, k, kk, neq, neq2, icc, nunk2

parameter (neq2=4)

dimension aa(neq2), ainv(neq2), w(nunk2), u(neq2), vi(neq), &

y(neq), l(neq2)

!

! ...Set aa equal to part of w to invert

!

DO 5 i=1,neq*neq

aa(i)=w(icc+i)

5 CONTINUE

!

! ...Initialize upper part of l and lower part of u to zero.

! and diagonal of u to one.

!

DO 9 i=1,neq

ii=(i-1)*neq

DO 7 j=1,neq

IF(i.EQ.j) u(ii+j)=1.0d0

IF(i.LT.j) l(ii+j)=0.0d0

IF(i.GT.j) u(ii+j)=0.0d0

7 CONTINUE

9 CONTINUE

!

! ...Find matrices l and u.

!

DO 121 j=1,neq

jj=(j-1)*neq

DO 113 i=j,neq

sum=0.0d0

ii=(i-1)*neq

IF(j.EQ.1) THEN

l(ii+j)=aa(ii+j)

ELSE

DO 111 k=1,j-1

kk=k-1

sum=sum + l(ii+k)*u(kk*neq+j)

111 CONTINUE

l(ii+j)=aa(ii+j) - sum

END IF

113 CONTINUE

IF(j.LT.neq) THEN

DO 117 i=j+1,neq

sum=0.0d0

ii=i-1

IF(j.EQ.1) THEN

u(i)=aa(i)/l(1)

ELSE

DO 115 k=1,j-1

kk=k-1

sum=sum + l(jj+k)*u(kk*neq+i)

115 CONTINUE

u(jj+i)=(aa(jj+i)-sum)/l(jj+j)

END IF

117 CONTINUE

END IF

121 CONTINUE

!

! ...Find inverse of aa.

!

! Use vector of identity matrix as r.h.s to find

! consecutive columns of aa'.

! l u a' = i

! l y = vi

! u va' = y

! where va' is a column of aa'.

!

DO 101 ic=1,neq

!

! ...Initialize vector of i.

!

DO 135 j=1,neq

IF(ic.EQ.j) THEN

vi(j)=1.0d0

ELSE

vi(j)=0.0d0

END IF

135 CONTINUE

!

! ...Perform forward substitution

!

y(1)=vi(1)/l(1)

DO 151 i=2,neq

ii=i-1

sum=0.0d0

DO 145 k=1,i-1

sum=sum+ l(ii*neq+k)*y(k)

145 CONTINUE

y(i)=(vi(i)-sum)/l(ii*neq+i)

151 CONTINUE

!

! ...Perform backward substitution

!

ainv(neq*neq-neq+ic)=y(neq)

DO 161 i=2,neq

in=neq-i+1

ii=(neq-i)*neq

sum=0.0d0

DO 155 k=in+1,neq

kk=k-1

sum=sum + u(ii+k)*ainv(kk*neq+ic)

155 CONTINUE

in1=in-1

ainv(in1*neq+ic)=y(in)-sum

161 CONTINUE

101 CONTINUE

!>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<

! itemp=7*neq*neq

! IF(icc.EQ.itemp) THEN

! WRITE(66,1999)

!1999 FORMAT(///,'inversion check, w(i), THEN wtemp(i)',//)

! DO 3000 i=1,neq

! WRITE(66,2000) (w(i,j),j=1,neq)

!3000 CONTINUE

! DO 3500 i=1,neq

! WRITE(66,2000) (wtemp(i,j),j=1,neq)

!3500 CONTINUE

! END IF

!2000 FORMAT(3(3x,d14.7))

!>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<

RETURN

END

!***********************************************************************

!

! Framework for Monte Carlo simulation using TRANSCAP-1D

!

!***********************************************************************

!

!

! Arrays

!

! cfes : vector of concentration of FeS (umol/g)

! twant : vector of desired time(s) for output

! x : vector of nodal coordinates (m)

!

!

! Variables

!

! Real

!

! alfa1, alfa2: irrigation coefficients (1/day & 1/m)

! alfa = alfa1*exp(alfa2*(x-xmax))

! area : column area (m**2)

! biox : maximum depth of bioturbation (m) (has to be a multiple of dx!)

! k1, k2 : coefficients of linear dissolution equation

! conc : concentration at the water/sediment interface (ug/m**3)

! diff : effective diffusion coefficient (m*2/day)

! dm : molecular diffusion

! disp : dispersivity (m)

! dt : time increment (day)

! dx : distance increment (m) (homogeneous grid!)

! epsil : time weighting factor (normally =1)

! gamma : cinetic coefficient for dissolution reaction

! pors : sediment porosity

! porscap : cap sediment porosity

! port0 : tube porosity at the surface of the sediment

! port=port0*exp(alfa2*(x-xmax))

! qs : advective flow in sediment (m**3/day)

! retard : retardation coefficient retard = 1+rho*xkd/pors

! tmax : maximum simulation time (day)

! vt : average water velocity in tube (m/day)

! xmax : lenght of column (m)

! r1 : inner radius of tubes/burrows

! r2 : half the distance between two tubes/burrows

! dp : depth of bioturbation

! d : distance from the burrow axis to the point where the concentration equals

! the horizontally integrated value

!

!

! Distribution parameters:

! ..dist(1) ..dist(2)

! if ..distyp=1=UNIFORM minimum maximum

! if ..distyp=2=NORMAL mean variance

! if ..distyp=3=EXPONENTIAL mean translation

!

!

! Integer

!

! ntub : number of tubes pro square meter

! nn : number of nodes

! neq : number of equations per node (tube + sediment)

! nunk : number of unknowns

! ntwant : number of desired time output

! nsim : number of desired montecarlo simualtions

! mc : if mc=0 then perform only one deterministic simulation

! if mc>0 then perform montecarlo simulation

!

!

! Input parameters

!

! line 1 : tmax, dt, xmax, dx, capdx

! line 2 : area, qs, pors, porscap

! line 3 : disp, diff, dm, epsil,

! line 4 : gamma, k1, k2, r1

! line 5 : ntwant

! line 6 : twant (1) --> twant(ntwant)

! line 7 : conc

! line 8 : cs(1) --> cs(nn)

! line 9 : ct(1) --> ct(nn)

! line 10 : cfes(1) à cfes (nn)

! line 11 : mc

! line 12 : nsim

!

! if montecarlo = 1

!

! line 13 : vt -> disttyp, dist(2), iseed

! line 14 : rt -> disttyp, dist(2), iseed

! line 15 : ntub -> disttyp, dist(2), iseed

! line 16 : dp -> disttyp, dist(2), iseed

!

! if montecarlo = 0

!

! line 13 : vt

! line 14 : retard

! line 15 : alfa1, port0

! line 16 : alfa2, biox

!

!

!***********************************************************************

!

program montecarlo

!

! ...define global variables

!

common alfa1, alfa2, area, biox, capdx, &

k1, k2, conc, dm, disp, diff, dp, dt, dx, epsil, &

gamma, pors, porscap, port0, qs, retard, &

tmax, vt, xmax

double precision alfa1, alfa2, area, biox, capdx, &

k1, k2, conc, dm, disp, diff, dp, dt, dx, epsil, &

gamma, pors, porscap, port0, qs, retard, &

tmax, vt, xmax

common sim, nsim, ntwant, nn, AllocateStatus, mc

integer sim, nsim, ntwant, nn, AllocateStatus, mc

common title

character *80 title

!

! ...define local variables

!

double precision r1, r2, bx, d, pi

double precision, dimension(:), allocatable :: x, cfes, &

cs, ct, twant, vtran, rtran, ntran, dpran, ranvec

double precision dist(2)

integer i, j, ii, disttyp, iseed, ntub

parameter (d = 0.0026)

parameter (pi = 3.14159265358987)

!

! ...open the working files

!

open(3, file='mcinput.dat')

open(4, file='mcflux.out')

open(5, file='mcoutput.out')

open(6, file='mccsurf.out')

!

! ...read the input

!

read(3, 300) title

300 format(a80)

read(3,*) tmax, dt, xmax, dx, capdx

read(3,*) area, qs, pors, porscap

read(3,*) disp, diff, dm, epsil

read(3,*) gamma, k1, k2, r1

!

! ...calculate number of nodes

!

nn = int(xmax/dx)

!

! ...read number for desired output(s)

!

read(3,*) ntwant

!

! ...allocate the memory for the arrays

!

allocate(x(nn), cs(nn), ct(nn), cfes(nn), stat = AllocateStatus)

if (AllocateStatus /= 0) stop "***Not enough memory***"

allocate(twant(ntwant), stat = AllocateStatus)

if (AllocateStatus /= 0) stop "***Not enough memory***"

!

! ...read time(s) for desired output(s)

!

read(3,*) (twant(i), i=1,ntwant)

!

! ...assign the first nodal x-coordinates

!

x(1)= dx/2

!

! ...assign nodal coordinates

!

do 10 i=1,nn-1

x(i+1)=x(i)+dx

10 continue

!

! ...read value of concentration at the sediment-water interface

!

read(3,*) conc

!

! ...read initial concentrations for sediment

!

read(3,*) (cs(i), i=1,nn)

!

! ...read initial concentrations for tube

!

read(3,*) (ct(i), i=1,nn)

!

! ...read concentration of FeS at each node

!

read(3,*) (cfes(i), i=1,nn)

!

! ...choose wheather or not want to do a montecarlo simulation

!

read(3,*) mc

!

! ...read number of desired simulations

!

read(3,*) nsim

!

! ...if not montecarlo simulation then nsim=1

!

if (mc.eq.0) then

nsim=1

endif

!

! ...allocate the memory for the arrays containing random numbers

!

allocate(ranvec(nsim), vtran(nsim), rtran(nsim), &

ntran(nsim), dpran(nsim), stat = AllocateStatus)

if (AllocateStatus /= 0) stop "***Not enough memory***"

!

! ...Write title of the simulation

!

write(4,404) title

write(5,404) title

write(6,404) title

404 format(//,60('*'),/,5x,'1D finite volume solution of the ', &

'transport equation',/,5x, a80,/,60('*'),//)

!

! ... if montecarlo simulation then...

!

if (mc.gt.0) then

!

! ...read distribution parameters

!

do 12 i= 1,4

read (3,*) disttyp

read (3,*) dist(1), dist(2)

read (3,*) iseed

do 14 j=1,nsim

ranvec(j)= 0.0d0

14 continue

!

! ...calculate random numbers for the distribution

! and put them in different vectors

!

call randist(disttyp,dist,iseed,nsim,ranvec)

if (i.eq.1) then

vtran=ranvec

elseif (i.eq.2) then

rtran=ranvec

elseif (i.eq.3) then

ntran=ranvec

elseif (i.eq.4) then

dpran=ranvec

endif

12 continue

!

! ...do montecarlo simulations with random parameters

!

!

do 10000 sim=1,nsim

vt=vtran(sim)

retard=rtran(sim)

ntub=NINT(ntran(sim))

dp=dpran(sim)

!

r2=1/sqrt(REAL(ntub))/2

alfa1=2*diff*r1/((r2**2-r1**2)*(d-r1))

port0=pi*(r1**2)*dx*ntub

bx=dp/dx

biox=NINT(bx)*dx

if (biox.le.0) then

biox=dx

endif

alfa2=-DLOG(alfa1/100)/biox

!

call solution(x, cfes, cs, ct,twant)

10000 continue

!

! ..if the user doesn't choose montecarlo carry out only one simulation

!

!

elseif (mc.eq.0) then

!

! ...read input

!

read (3,*) vt

read (3,*) retard

read (3,*) alfa1, port0

read (3,*) alfa2, biox

!

call solution(x, cfes, cs, ct,twant)

!

else

write(*,*) 'Error: variable montecarlo must be a positive integer or 0'

stop

endif

end

!******************************************************************

!

! This subroutine returns a vector containing a random variable

! with the given distribution type.

!

!******************************************************************

!

subroutine randist(distyp,dist,iseed,nsim,ranvec)

!

implicit real*8(a-h,o-z)

implicit integer*4(i-n)

double precision ranvec(nsim)

!

! ARGUMENTS:

!

! INPUT:

!

! type of distribution

!

integer*4 distyp

!

! distribution parameters:

!

real*8 dist(2)

! dist(1) dist(2)

! if distyp=1=UNIFORM minimum maximum

! if distyp=2=NORMAL mean variance

! if distyp=3=EXPONENTIAL mean translation

!

!******************************************************************

!

! LOCAL DECLARATIONS:

!

!

real*8 ran3

real*8 p

integer*4 i, j

!

!******************************************************************

!

! ...read the input

!

do 10 j=1,nsim

if(distyp .eq. 1) then

ranum=dist(1)+(dist(2)-dist(1))*dble(ran3(iseed))

!was ranum=dist(1)+(dist(2)-dist(1))*dble(rand())

!cc ranum=dist(1)+(dist(2)-dist(1))*ran1(iseed)

elseif(distyp .eq. 2) then

p=0.0d0

do 100 i=1,12

p=p+dble(ran3(iseed))

!cc p=p+ran1(iseed)

!was p=p+dble(rand())

100 continue

ranum=dsqrt(dist(2))*(p-6.0d0)+dist(1)

elseif(distyp .eq. 3) then

p=-dabs(dist(1)-dist(2))*dlog(1.0d0-dble(ran3(iseed)))

!was p=-dabs(dist(1)-dist(2))*dlog(1.0d0-dble(rand()))

!cc p=-dabs(dist(1)-dist(2))*dlog(1.0d0-ran1(iseed))

if(dist(1) .gt. dist(2)) then

ranum=dist(2)+p

else

ranum=dist(2)-p

endif

endif

ranvec(j)=ranum

!

! If you don't want the random number to be negative

!

! ranvec(j)= dabs(ranum)

!

10 continue

return

end

!***********************************************************************

!

! This function returns a random variable

!

!***********************************************************************

real*8 FUNCTION RAN3(IDUM)

implicit real*8(a-h,o-z)

integer*4 MA(55)

real*8 fac

integer*4 mbig,mseed,iff,mz,mj,mk,ii,inext,inextp,idum

save inext,inextp,ma

DATA IFF /0/

MBIG=1000000000

MSEED=161803398

MZ=0

FAC=1.E-9

IF(IDUM.LT.0.OR.IFF.EQ.0)THEN

IFF=1

MJ=MSEED-IABS(IDUM)

MJ=MOD(MJ,MBIG)

MA(55)=MJ

MK=1

DO 11 I=1,54

II=MOD(21*I,55)

MA(II)=MK

MK=MJ-MK

IF(MK.LT.MZ)MK=MK+MBIG

MJ=MA(II)

11 CONTINUE

DO 13 K=1,4

DO 12 I=1,55

MA(I)=MA(I)-MA(1+MOD(I+30,55))

IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG

12 CONTINUE

13 CONTINUE

INEXT=0

INEXTP=31

IDUM=1

ENDIF

INEXT=INEXT+1

IF(INEXT.EQ.56)INEXT=1

INEXTP=INEXTP+1

IF(INEXTP.EQ.56)INEXTP=1

MJ=MA(INEXT)-MA(INEXTP)

IF(MJ.LT.MZ)MJ=MJ+MBIG

MA(INEXT)=MJ

RAN3=MJ*FAC

RETURN

end