https://github.com/cran/ReacTran
Tip revision: e1145574dc02ebe9e0d96419755e4702aef85172 authored by Karline Soetaert on 16 September 2010, 00:00:00 UTC
version 1.3
version 1.3
Tip revision: e114557
tran.1D.f
C#==============================================================================
C Diffusion in a 1-dimensional finite difference grid
C all inputs are vectors
C#==============================================================================
SUBROUTINE diff1D (N, C, Cup, Cdown, fluxup, fluxdown, aup, adown, &
& BcUp, BcDown, D, VF, A, dx, dxaux, &
& Flux, dC)
IMPLICIT NONE
INTEGER N ! length of C
C input
DOUBLE PRECISION C(N)
C Boundary concentrations (used if Bc..=2,4), fluxes (used if Bc= 1)
C and convection coeff (used if Bc=4)
DOUBLE PRECISION Cup, Cdown, fluxup, fluxdown, aup, adown
C Diffusion, volume fraction, surface Area
DOUBLE PRECISION D(N+1), VF(N+1), A(N+1)
C grid size, distance from mid to mid
DOUBLE PRECISION dx(N), dxaux(N+1)
C boundary concitions (1= flux, 2=conc, 3 = 0-grad, 4=convect)
INTEGER BcUp, BcDown
C output: fluxes and rate of change
DOUBLE PRECISION Flux(N+1), dC(N)
C locals
INTEGER I
DOUBLE PRECISION AVF
C -------------------------------------------------------------------------------
C Flux - first internal cells
DO I = 2,N
Flux(I) = -VF(I)*D(I) * (C(I)-C(I-1)) /dxaux(I)
ENDDO
C Then the outer cells
C upstream boundary
IF (BcUp .EQ. 1) THEN
Flux(1) = fluxup
ELSE IF (BcUp .EQ. 2) THEN
Flux(1) = -VF(1)*D(1) * (C(1)-Cup) /dxaux(1)
ELSE IF (BcUp .EQ. 3) THEN
Flux(1) = 0.D0
ELSE
Flux(1) = aup * (Cup - C(1))
ENDIF
C downstream boundary
IF (BcDown .EQ. 1) THEN
Flux(N+1) = fluxdown
ELSE IF (BcDown .EQ. 2) THEN
Flux(N+1) = -VF(N+1)*D(N+1) * (Cdown-C(N)) /dxaux(N+1)
ELSE IF (BcDown .EQ. 3) THEN
Flux(N+1) =0.D0
ELSE
Flux(N+1) = -adown * (Cdown-C(N))
ENDIF
C Rate of change = negative flux gradient
DO I = 1,N
AVF = 0.25 * (A(I)+A(I+1)) * (VF(I)+VF(I+1))
dC(I) = -(A(I+1)*Flux(I+1) - A(I)*Flux(I)) / AVF / dx(I)
ENDDO
RETURN
END SUBROUTINE diff1D
C#==============================================================================
C Diffusion in a one-dimensional finite difference grid
C all inputs except C are scalars
C#==============================================================================
SUBROUTINE diff1Dscal (N, C, Cup, Cdown, fluxup, fluxdown, &
& aup, adown,BcUp, BcDown, D, VF, dx, &
& Flux, dC)
IMPLICIT NONE
INTEGER N ! length of C
C input
DOUBLE PRECISION C(N), Cup, Cdown, fluxup, fluxdown, aup, adown
DOUBLE PRECISION D, VF, dx
C boundary conditions (1= flux, 2 = conc, 3 = 0 grad, 4= convective)
INTEGER BcUp, BcDown
C output
DOUBLE PRECISION Flux(N+1), dC(N)
C locals
INTEGER I
C -------------------------------------------------------------------------------
C Flux - first internal cells
DO I = 2,N
Flux(I) = -VF*D * (C(I)-C(I-1)) /dx
ENDDO
C Then the outer cells
C upstream
IF (BcUp .EQ. 1) THEN
Flux(1) = fluxup
ELSE IF (BcUp .EQ. 2) THEN
Flux(1) = -VF*D * (C(1)-Cup) /dx
ELSE IF (BcUp .EQ. 3) THEN
Flux(1) = 0.D0
ELSE
Flux(1) = -aup * (C(1)-Cup)
ENDIF
C downstream
IF (BcDown .EQ. 1) THEN
Flux(N+1) = fluxdown
ELSE IF (BcDown .EQ. 2) THEN
Flux(N+1) = -VF*D * (Cdown-C(N)) /dx
ELSE IF (BcDown .EQ. 3) THEN
Flux(N+1) = 0.D0
ELSE
Flux(N+1) = -adown * (Cdown-C(N))
ENDIF
C rate of change = negative flux gradient
DO I = 1,N
dC(I) = -(Flux(I+1) - Flux(I)) /VF /dx
ENDDO
RETURN
END