!

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

!

! 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