tran.3D.f
C#==============================================================================
C Diffusion in a 3-dimensional finite difference grid
C all inputs are scalars
C#==============================================================================
SUBROUTINE diff3Dscal (Nx, Ny, Nz, C, &
& Cxup, Cxdown, Cyup, Cydown, Czup, Czdown, &
& Fxup, Fxdown, Fyup, Fydown, Fzup, Fzdown, &
& axup, axdown, ayup, aydown, azup, azdown, &
& BCxUp, BCxDown, BCyUp, BCyDown, BCzUp, BCzDown, &
& D_x, D_y, D_z, VF, dx, dy, dz, &
& FluxXup, FluxXdown, FluxYup, FluxYdown, &
& FluxZup, FluxZdown, dC)
IMPLICIT NONE
INTEGER Nx, Ny, Nz ! dimension of C
C input
DOUBLE PRECISION C(Nx,Ny,Nz)
C when Bc.. = 2,4; boundary concentrations
DOUBLE PRECISION Cxup, Cxdown, Cyup, Cydown, Czup, Czdown
C when Bc.. = 1: boundary fluxes
DOUBLE PRECISION Fxup, Fxdown, Fyup, Fydown, Fzup, Fzdown
C when Bc.. = 4: boundary convective coefficients
DOUBLE PRECISION axup, axdown, ayup, aydown, azup, azdown
C diffusion coefficients, volume fraction, grid sizes
DOUBLE PRECISION D_x, D_y, D_z, VF, dx, dy, dz
C boundaries: 1= flux, 2 = conc, 3 = 0 grad, 4=convection
INTEGER BCxUp, BCxDown, BCyUp, BCyDown, BCzUp, BCzDown
C output: fluxes out and inside the domain..
DOUBLE PRECISION FluxXup(Ny, Nz), FluxXdown(Ny, Nz)
DOUBLE PRECISION FluxYup(Nx, Nz), FluxYdown(Nx, Nz)
DOUBLE PRECISION FluxZup(Ny, Nz), FluxZdown(Ny, Nz)
C rate of change
DOUBLE PRECISION dC(Nx,Ny,Nz)
C locals
INTEGER I, J, K
DOUBLE PRECISION Flux(Nx+1, Ny+1, Nz+1)
C -------------------------------------------------------------------------------
C First diffusion internal cells
DO I = 2, Nx
DO J = 1, Ny
DO K = 1, Nz
Flux(I,J,K) = -VF*D_x*( (C(I,J,K)-C(I-1,J,K)) /dx)
ENDDO
ENDDO
ENDDO
Flux(1,:,:) = 0.D0
Flux(Nx+1,:,:) = 0.D0
C Then the outer cells - depending on boundary
C upstream X-boundary
IF (BCxUp .EQ. 1) THEN
DO J = 1, Ny
DO K = 1, Nz
Flux(1,J,K) = Fxup
ENDDO
ENDDO
ELSE IF (BCxUp .EQ. 2) THEN
DO J = 1, Ny
DO K = 1, Nz
Flux(1,J,K) = -VF*D_x*(C(1,J,K)-Cxup)/dx
ENDDO
ENDDO
ELSE IF (BCxUp .EQ. 3) THEN
Flux(1,:,:) = 0.D0
ELSE IF (BCxUp .EQ. 4) THEN
DO J = 1, Ny
DO K = 1, Nz
Flux(1,J,K) = -axup*(C(1,J,K)-Cxup)
ENDDO
ENDDO
ENDIF
DO J = 1, Ny
DO K = 1, Nz
FluxXup(J,K) = Flux(1,J,K)
ENDDO
ENDDO
C downstream X-boundary
IF (BCxDown .EQ. 1) THEN
DO J = 1, Ny
DO K = 1, Nz
Flux(Nx+1,J,K) = Fxdown
ENDDO
ENDDO
ELSE IF (BCxDown .EQ. 2) THEN
DO J = 1, Ny
DO K = 1, Nz
Flux(Nx+1,J,K) = -VF*D_x*(Cxdown-C(Nx,J,K)) /dx
ENDDO
ENDDO
ELSE IF (BCxDown .EQ. 3) THEN
Flux(Nx+1,:,:) = 0.D0
ELSE IF (BCxDown .EQ. 4) THEN
DO J = 1, Ny
DO K = 1, Nz
Flux(Nx+1,J,K) = -axdown*(Cxdown-C(Nx,J,K))
ENDDO
ENDDO
ENDIF
DO J = 1, Ny
DO K = 1, Nz
FluxXdown(J,K) = Flux(Nx+1,J,K)
ENDDO
ENDDO
C Rate of change = negative flux gradient
DO I = 1,Nx
DO J = 1,Ny
DO K = 1, Nz
dC(I,J,K) = -(Flux(I+1,J,K)-Flux(I,J,K))/VF/dx
ENDDO
ENDDO
ENDDO
C Y-direction
C first internal cells
DO I = 1, Nx
DO J = 2, Ny
DO K = 1, Nz
Flux(I,J,K) = -VF*D_y*( (C(I,J,K)-C(I,J-1,K)) /dy)
ENDDO
ENDDO
ENDDO
Flux(:,1,:) = 0.D0
Flux(:,Ny+1,:) = 0.D0
C upstream Y-boundary
IF (BCyUp .EQ. 1) THEN
DO I = 1, Nx
DO K = 1, Nz
Flux(I,1,K) = Fyup
ENDDO
ENDDO
ELSE IF (BCyUp .EQ. 2) THEN
DO I = 1, Nx
DO K = 1, Nz
Flux(I,1,K) = -VF*D_y*(C(I,1,K)-Cyup)/dy
ENDDO
ENDDO
ELSE IF (BCyUp .EQ. 3) THEN
Flux(:,1,:) = 0.D0
ELSE IF (BCyUp .EQ. 4) THEN
DO I = 1, Nx
DO K = 1, Nz
Flux(I,1,K) = -ayup*(C(I,1,K)-Cyup)
ENDDO
ENDDO
ENDIF
DO I = 1, Nx
DO K = 1, Nz
FluxYup(I,K) = Flux(I,1,K)
ENDDO
ENDDO
C downstream Y-boundary
IF (BCyDown .EQ. 1) THEN
DO I = 1, Nx
DO K = 1, Nz
Flux(I,Ny+1,K) = Fydown
ENDDO
ENDDO
ELSE IF (BCyDown .EQ. 2) THEN
DO I = 1, Nx
DO K = 1, Nz
Flux(I,Ny+1,K) = -VF*D_y*(Cydown-C(I,Ny,K))/dy
ENDDO
ENDDO
ELSE IF (BCyDown .EQ. 3) THEN
Flux(:,Ny+1,:) = 0.D0
ELSE IF (BCyDown .EQ. 4) THEN
DO I = 1, Nx
DO K = 1, Nz
Flux(I,Ny+1,K) = -aydown*(Cydown-C(I,Ny,K))
ENDDO
ENDDO
ENDIF
DO I = 1, Nx
DO K = 1, Nz
FluxYdown(I,K) = Flux(I,Ny+1,K)
ENDDO
ENDDO
C Rate of change = negative flux gradient
DO I = 1,Nx
DO J = 1,Ny
DO K = 1, Nz
dC(I,J,K) = dC(I,J,K) -(Flux(I,J+1,K)-Flux(I,J,K))/VF/dy
ENDDO
ENDDO
ENDDO
C Z-direction
C first internal cells
DO I = 1, Nx
DO J = 1, Ny
DO K = 2, Nz
Flux(I,J,K) = -VF*D_z*( (C(I,J,K)-C(I,J,K-1)) /dz)
ENDDO
ENDDO
ENDDO
Flux(:,:,1) = 0.D0
Flux(:,:,Nz+1) = 0.D0
C upstream Y-boundary
IF (BCzUp .EQ. 1) THEN
DO I = 1, Nx
DO J = 1, Ny
Flux(I,J,1) = Fzup
ENDDO
ENDDO
ELSE IF (BCzUp .EQ. 2) THEN
DO I = 1, Nx
DO J = 1, Ny
Flux(I,J,1) = -VF*D_z*(C(I,J,1)-Czup)/dz
ENDDO
ENDDO
ELSE IF (BCzUp .EQ. 3) THEN
Flux(:,:,1) = 0.D0
ELSE IF (BCyUp .EQ. 4) THEN
DO I = 1, Nx
DO J = 1, Ny
Flux(I,J,1) = -azup*(C(I,J,1)-Czup)
ENDDO
ENDDO
ENDIF
DO I = 1, Nx
DO J = 1, Ny
FluxZup(I,J) = Flux(I,J,1)
ENDDO
ENDDO
C downstream Z-boundary
IF (BCzDown .EQ. 1) THEN
DO I = 1, Nx
DO J = 1, Ny
Flux(I,J,Nz+1) = Fzdown
ENDDO
ENDDO
ELSE IF (BCzDown .EQ. 2) THEN
DO I = 1, Nx
DO J = 1, Ny
Flux(I,J,Nz+1) = -VF*D_z*(Czdown-C(I,J,Nz))/dz
ENDDO
ENDDO
ELSE IF (BCzDown .EQ. 3) THEN
Flux(:,:,Nz+1) = 0.D0
ELSE IF (BCzDown .EQ. 4) THEN
DO I = 1, Nx
DO J = 1, Ny
Flux(I,J,Nz+1) = -azdown*(Czdown-C(I,J,Nz))
ENDDO
ENDDO
ENDIF
DO I = 1, Nx
DO J = 1, Ny
FluxZdown(I,J) = Flux(I,J,Nz+1)
ENDDO
ENDDO
C Rate of change = negative flux gradient
DO I = 1,Nx
DO J = 1,Ny
DO K = 1, Nz
dC(I,J,K) = dC(I,J,K) -(Flux(I,J,K+1)-Flux(I,J,K))/VF/dz
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE Diff3Dscal
C#==============================================================================
C Diffusion in a 3-dimensional spherical coordinate system (r, theta, phi)
C all inputs are vectors/scalars
C CHECK THIS!!!
C#==============================================================================
SUBROUTINE diffspher (Nr, Ntet, Nphi, C, &
& Crup, Crdown, Ctetup, Ctetdown, Cphiup, Cphidown, &
& Frup, Frdown, Ftetup, Ftetdown, Fphiup, Fphidown, &
& Bcrup, Bcrdown, Bctetup, Bctetdown, &
& Bcphiup, Bcphidown, D_r, D_tet, D_phi, &
& r, tet, phi, FluxrUp, FluxrDown, &
& FluxtetUp, FluxtetDown, FluxphiUp, FluxPhiDown, dC)
IMPLICIT NONE
C dimension of
INTEGER Nr, Ntet, Nphi
C concentration
DOUBLE PRECISION C(Nr, Ntet, Nphi)
C used when Bc.. = 2
DOUBLE PRECISION Crup, Crdown, Ctetup, Ctetdown, Cphiup, Cphidown
C used when Bc.. = 1
DOUBLE PRECISION Frup, Frdown, Ftetup, Ftetdown, Fphiup, Fphidown
DOUBLE PRECISION D_r(Nr+1), D_tet(Ntet+1), D_phi(Nphi+1), &
& r(Nr+1), tet(Ntet+1), phi(Nphi+1), &
& rc(Nr), tetc(Ntet), phic(Nphi)
DOUBLE PRECISION dr(Nr),draux(Nr+1),dtet(Ntet),dtetaux(Ntet+1), &
& dphi(Nphi), dphiaux(Nphi+1)
c boundaries: 1= flux, 2 = conc, 3 = zero- grad, 5 = cyclic boundary
INTEGER BcRup, BcRdown, BcTetup, BcTetdown, BcPhiup, BcPhidown
C output
DOUBLE PRECISION FluxrUp(Ntet,Nphi), FluxrDown(Ntet,Nphi)
DOUBLE PRECISION Fluxtetup(Nr,Nphi), FluxtetDown(Nr,Nphi)
DOUBLE PRECISION Fluxphiup(Nr,Ntet), FluxphiDown(Nr,Ntet)
DOUBLE PRECISION dC(Nr,Ntet,Nphi)
C locals
INTEGER I, J, K
DOUBLE PRECISION Flux(Nr+1, Ntet+1, Nphi+1), Cbound
C -------------------------------------------------------------------------------
C The grid
DO I = 1, Nr
rc(I) = 0.5*(r(I)+r(I+1))
dr(I) = r(I+1) - r(I)
ENDDO
DO I = 1, Nr-1
draux(I+1) = (rc(I+1) - rc(I))
ENDDO
draux(1) = rc(1)-r(1)
draux(Nr+1) = r(Nr+1) - rc(Nr)
DO I = 1, Ntet
tetc(I) = 0.5*(tet(I)+tet(I+1))
dtet(I) = tet(I+1) - tet(I)
ENDDO
DO I = 1, Ntet-1
dtetaux(I+1) = (tetc(I+1) - tetc(I))
ENDDO
dtetaux(1) = tetc(1)-tet(1)
dtetaux(Ntet+1) = tet(Ntet+1) - tetc(Ntet)
DO I = 1, Nphi
phic(I) = 0.5*(phi(I)+phi(I+1))
dphi(I) = phi(I+1) - phi(I)
ENDDO
DO I = 1, Nphi-1
dphiaux(I+1) = (phic(I+1) - phic(I))
ENDDO
dphiaux(1) = phic(1)-phi(1)
dphiaux(Nphi+1) = phi(Nphi+1) - phic(Nphi)
C ## transport in r-direction ##
C First diffusion internal cells
DO I = 2, Nr
DO J = 1, Ntet
DO K = 1, Nphi
Flux(I,J,K) = -D_r(I)*(C(I,J,K)-C(I-1,J,K)) /draux(I)
ENDDO
ENDDO
ENDDO
Flux(1,:,:) = 0.D0
Flux(Nr+1,:,:) = 0.D0
C Then the outer cells - depending on boundary type
C upstream r-boundary
IF (Bcrup .EQ. 1) THEN
DO J = 1, Ntet
DO K = 1, Nphi
Flux(1,J,K) = Frup
ENDDO
ENDDO
ELSE IF (Bcrup .EQ. 2) THEN
DO J = 1, Ntet
DO K = 1, Nphi
Flux(1,J,K) = -D_r(1) * (C(1,J,K)-Crup) /draux(1)
ENDDO
ENDDO
ELSE IF (Bcrup .EQ. 3) THEN
Flux(1,:,:) = 0.D0
ELSE IF (Bcrup .EQ. 5) THEN
DO J = 1, Ntet
DO K = 1, Nphi
Cbound = (C(1,J,K)*draux(Nr+1) + C(Nr,J,K)*draux(1)) / &
& (draux(1)+draux(Nr+1))
Flux(1,J,K) = -D_r(1) * (C(1,J,K)-Cbound ) /draux(1)
Flux(Nr+1,J,K) = -D_r(Nr+1) * (Cbound-C(Nr,J,K)) /draux(Nr+1)
ENDDO
ENDDO
ENDIF
DO J = 1, Ntet
DO K = 1, Nphi
Fluxrup(J,K) = Flux(1,J,K)
ENDDO
ENDDO
C downstream r-boundary
IF (Bcrdown .EQ. 1) THEN
DO J = 1, Ntet
DO K = 1, Nphi
Flux(Nr+1,J,K) = Frdown
ENDDO
ENDDO
ELSE IF (Bcrdown .EQ. 2) THEN ! if 5: Crdown already estimated
DO J = 1, Ntet
DO K = 1, Nphi
Flux(Nr+1,J,K) =-D_r(Nr+1) * (Crdown-C(Nr,J,K)) /draux(Nr+1)
ENDDO
ENDDO
ELSE IF (Bcrdown .EQ. 3) THEN
Flux(Nr+1,:,:) = 0.D0
ENDIF
DO J = 1, Ntet
DO K = 1, Nphi
Fluxrdown(J,K) = Flux(Nr+1,J,K)
ENDDO
ENDDO
C Rate of change in r-direction
DO I = 1,Nr
DO J = 1,Ntet
DO K = 1, Nphi
IF(rc(I) .NE. 0.D0) &
& dC(I,J,K) = &
& -(Flux(I+1,J,K)*(r(I+1)**2.)-Flux(I,J,K)*(r(I)**2.)) &
& /dr(I)/(rc(I)**2.)
ENDDO
ENDDO
ENDDO
C ## transport in teta -direction ##
C First diffusion internal cells
DO I = 1, Nr
DO J = 2, Ntet
DO K = 1, Nphi
Flux(I,J,K) = -D_tet(J)*(C(I,J,K)-C(I,J-1,K)) /dtetaux(J) &
& /rc(I)
ENDDO
ENDDO
ENDDO
Flux(:,1,:) = 0.D0
Flux(:,Ntet+1,:) = 0.D0
C upstream teta-boundary
IF (BcTetup .EQ. 1) THEN
DO I = 1, Nr
DO K = 1, Nphi
Flux(I,1,K) = Ftetup
ENDDO
ENDDO
ELSE IF (BcTetup .EQ. 2) THEN
DO I = 1, Nr
DO K = 1, Nphi
Flux(I,1,K) = -D_tet(1)*(C(I,1,K)-Ctetup)/dtetaux(1)/rc(I)
ENDDO
ENDDO
ELSE IF (BcTetup .EQ. 3) THEN
Flux(:,1,:) = 0.D0
ELSE IF (BcTetup .EQ. 5) THEN
DO I = 1, Nr
DO K = 1, Nphi
Cbound = (C(I,1,K)*dtetaux(Ntet+1) + C(I,Ntet,K)*dtetaux(1))/ &
& (dtetaux(1)+dtetaux(Ntet+1))
Flux(I,1,K) = -D_tet(1) *(C(I,1,K)-Cbound) &
& /dtetaux(1)/rc(I)
Flux(I,Ntet+1,K) = -D_tet(Ntet+1)*(Cbound-C(I,Ntet,K)) &
& /dtetaux(Ntet+1)/rc(I)
ENDDO
ENDDO
ENDIF
DO I = 1, Nr
DO K = 1, Nphi
Fluxtetup(I,K) = Flux(I,1,K)
ENDDO
ENDDO
C downstream teta-boundary
IF (BcTetdown .EQ. 1) THEN
DO I = 1, Nr
DO K = 1, Nphi
Flux(I,Ntet+1,K) = Ftetdown
ENDDO
ENDDO
ELSE IF (BcTetdown .EQ. 2) THEN
DO I = 1, Nr
DO K = 1, Nphi
Flux(I,Ntet+1,K) = -D_tet(Ntet+1) * (Ctetdown-C(I,Ntet,K)) &
& /dtetaux(Ntet+1)/rc(I)
ENDDO
ENDDO
ELSE IF (BcTetdown .EQ. 3) THEN
Flux(:,Ntet+1,:) = 0.D0
ENDIF
DO I = 1, Nr
DO K = 1, Nphi
Fluxtetdown(I,K) = Flux(I,Ntet+1,K)
ENDDO
ENDDO
C rate of change, negative flux gradient
DO I = 1,Nr
DO J = 1,Ntet
DO K = 1, Nphi
if (rc(I) .NE. 0.D0 .AND. sin(tetc(J)) .NE. 0.D0) &
& dC(I,J,K) = dC(I,J,K) &
& -(Flux(I,J+1,K)*sin(tet(J+1))-Flux(I,J,K)*sin(tet(J))) &
& /dtet(J)/rc(I)/sin(tetc(J))
ENDDO
ENDDO
ENDDO
C ## transport in phi-direction ##
C First diffusion internal cells
DO I = 1, Nr
DO J = 1, Ntet
DO K = 2, Nphi
Flux(I,J,K) = -D_phi(J)*(C(I,J,K)-C(I,J,K-1)) &
& /dphiaux(K) /rc(I)/sin(tet(j))
ENDDO
ENDDO
ENDDO
Flux(:,:,1) = 0.D0
Flux(:,:,Nphi+1) = 0.D0
C upstream phi-boundary
IF (BcPhiup .EQ. 1) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,1) = Fphiup
ENDDO
ENDDO
ELSE IF (BcPhiup .EQ. 2) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,1) = -D_phi(1)*(C(I,J,1)-Cphiup) &
& /dphiaux(1)/rc(I) /sin(tet(j))
ENDDO
ENDDO
ELSE IF (BcPhiup .EQ. 3) THEN
Flux(:,:,1) = 0.D0
ELSE IF (BcPhiup .EQ. 5) THEN
DO I = 1, Nr
DO J = 1, Ntet
Cbound = (C(I,J,1)*dphiaux(Nphi+1) + C(I,J,Nphi)*dphiaux(1))/ &
& (dphiaux(1)+dphiaux(Nphi+1))
Flux(I,J,1) = -D_phi(1)*(C(I,J,1)-Cbound) &
& /dphiaux(1) /rc(I) /sin(tet(j))
Flux(I,J,Nphi+1) = -D_phi(Nphi+1)*(Cbound-C(I,J,Nphi)) &
& /dphiaux(Nphi+1)/rc(I) /sin(tet(j))
ENDDO
ENDDO
ENDIF
DO I = 1, Nr
DO J = 1, Ntet
Fluxphiup(I,J) = Flux(I,J,1)
ENDDO
ENDDO
C downstream phi-boundary
IF (Bcphidown .EQ. 1) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,Nphi+1) = Fphidown
ENDDO
ENDDO
ELSE IF (Bcphidown .EQ. 2) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,Nphi+1) = -D_phi(Nphi+1) * (Cphidown-C(I,J,Nphi)) &
& /dphiaux(Nphi+1)/rc(I) /sin(tet(j))
ENDDO
ENDDO
ELSE IF (Bcphidown .EQ. 3) THEN
Flux(:,:,Nphi+1) = 0.D0
ENDIF
DO I = 1, Nr
DO J = 1, Ntet
Fluxphidown(I,J) = Flux(I,J,Nphi+1)
ENDDO
ENDDO
C rate of change, negative flux gradient
DO I = 1, Nr
DO J = 1, Ntet
DO K = 1, Nphi
if (rc(I) .NE. 0.D0 .AND. sin(tetc(J)) .NE. 0.D0) &
& dC(I,J,K) = dC(I,J,K) &
& -(Flux(I,J,K+1)-Flux(I,J,K))/dphi(K)/rc(I)/sin(tetc(J))
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE diffspher
C#==============================================================================
C Diffusion in a 3-dimensional cylindrical coordinate system (r, theta, z)
C all inputs are vectors/scalars
C CHECK THIS!!!
C#==============================================================================
SUBROUTINE diffcyl (Nr, Ntet, Nz, C, &
& Crup, Crdown, Ctetup, Ctetdown, Czup, Czdown, &
& Frup, Frdown, Ftetup, Ftetdown, Fzup, Fzdown, &
& Bcrup, Bcrdown, Bctetup, Bctetdown, &
& BCzUp, BCzDown, D_r, D_tet, D_z, &
& r, tet, z, FluxrUp, FluxrDown, &
& FluxtetUp, FluxtetDown, FluxzUp, FluxzDown, dC)
IMPLICIT NONE
C dimension of
INTEGER Nr, Ntet, Nz
C concentration
DOUBLE PRECISION C(Nr, Ntet, Nz)
C used when Bc.. = 2
DOUBLE PRECISION Crup, Crdown, Ctetup, Ctetdown, Czup, Czdown
C used when Bc.. = 1
DOUBLE PRECISION Frup, Frdown, Ftetup, Ftetdown, Fzup, Fzdown
DOUBLE PRECISION D_r(Nr+1), D_tet(Ntet+1), D_z(Nz+1), &
& r(Nr+1), tet(Ntet+1), z(Nz+1), &
& rc(Nr), tetc(Ntet), zc(Nz)
DOUBLE PRECISION dr(Nr), draux(Nr+1), dtet(Ntet),dtetaux(Ntet+1), &
& dz(Nz), dzaux(Nz+1)
c boundaries: 1= flux, 2 = conc, 3 = zero- grad, 5 = cyclic boundary
INTEGER BcRup, BcRdown, BcTetup, BcTetdown, BCzUp, BCzDown
C output
DOUBLE PRECISION FluxrUp(Ntet,Nz), FluxrDown(Ntet,Nz)
DOUBLE PRECISION Fluxtetup(Nr,Nz), FluxtetDown(Nr,Nz)
DOUBLE PRECISION FluxZup(Nr,Ntet), FluxZDown(Nr,Ntet)
DOUBLE PRECISION dC(Nr,Ntet,Nz)
C locals
INTEGER I, J, K
DOUBLE PRECISION Flux(Nr+1, Ntet+1, Nz+1), Cbound
C -------------------------------------------------------------------------------
C The grid
DO I = 1, Nr
rc(I) = 0.5*(r(I)+r(I+1))
dr(I) = r(I+1) - r(I)
ENDDO
DO I = 1, Nr-1
draux(I+1) = (rc(I+1) - rc(I))
ENDDO
draux(1) = rc(1)-r(1)
draux(Nr+1) = r(Nr+1) - rc(Nr)
DO I = 1, Ntet
tetc(I) = 0.5*(tet(I)+tet(I+1))
dtet(I) = tet(I+1) - tet(I)
ENDDO
DO I = 1, Ntet-1
dtetaux(I+1) = (tetc(I+1) - tetc(I))
ENDDO
dtetaux(1) = tetc(1)-tet(1)
dtetaux(Ntet+1) = tet(Ntet+1) - tetc(Ntet)
DO I = 1, Nz
zc(I) = 0.5*(z(I)+z(I+1))
dz(I) = z(I+1) - z(I)
ENDDO
DO I = 1, Nz-1
dzaux(I+1) = (zc(I+1) - zc(I))
ENDDO
dzaux(1) = zc(1)-z(1)
dzaux(Nz+1) = z(Nz+1) - zc(Nz)
C ## transport in r-direction ##
C First diffusion internal cells
DO I = 2, Nr
DO J = 1, Ntet
DO K = 1, Nz
Flux(I,J,K) = -D_r(I)*(C(I,J,K)-C(I-1,J,K)) /draux(I)
ENDDO
ENDDO
ENDDO
Flux(1,:,:) = 0.D0
Flux(Nr+1,:,:) = 0.D0
C Then the outer cells - depending on boundary type
C upstream r-boundary
IF (Bcrup .EQ. 1) THEN
DO J = 1, Ntet
DO K = 1, Nz
Flux(1,J,K) = Frup
ENDDO
ENDDO
ELSE IF (Bcrup .EQ. 2) THEN
DO J = 1, Ntet
DO K = 1, Nz
Flux(1,J,K) = -D_r(1) * (C(1,J,K)-Crup) /draux(1)
ENDDO
ENDDO
ELSE IF (Bcrup .EQ. 3) THEN
Flux(1,:,:) = 0.D0
ELSE IF (Bcrup .EQ. 5) THEN
DO J = 1, Ntet
DO K = 1, Nz
Cbound = (C(1,J,K)*draux(Nr+1) + C(Nr,J,K)*draux(1)) / &
& (draux(1)+draux(Nr+1))
Flux(1,J,K) = -D_r(1) * (C(1,J,K)-Cbound ) /draux(1)
Flux(Nr+1,J,K) = -D_r(Nr+1) * (Cbound-C(Nr,J,K)) /draux(Nr+1)
ENDDO
ENDDO
ENDIF
DO J = 1, Ntet
DO K = 1, Nz
Fluxrup(J,K) = Flux(1,J,K)
ENDDO
ENDDO
C downstream r-boundary
IF (Bcrdown .EQ. 1) THEN
DO J = 1, Ntet
DO K = 1, Nz
Flux(Nr+1,J,K) = Frdown
ENDDO
ENDDO
ELSE IF (Bcrdown .EQ. 2) THEN ! if 5: Crdown already estimated
DO J = 1, Ntet
DO K = 1, Nz
Flux(Nr+1,J,K) =-D_r(Nr+1) * (Crdown-C(Nr,J,K)) /draux(Nr+1)
ENDDO
ENDDO
ELSE IF (Bcrdown .EQ. 3) THEN
Flux(Nr+1,:,:) = 0.D0
ENDIF
DO J = 1, Ntet
DO K = 1, Nz
Fluxrdown(J,K) = Flux(Nr+1,J,K)
ENDDO
ENDDO
C Rate of change in r-direction
DO I = 1,Nr
DO J = 1,Ntet
DO K = 1, Nz
IF (rc(I) .NE. 0.D0) &
& dC(I,J,K) = &
& -(Flux(I+1,J,K)*r(I+1) - Flux(I,J,K)*r(I))/dr(I)/(rc(I))
ENDDO
ENDDO
ENDDO
C ## transport in teta -direction ##
C First diffusion internal cells
DO I = 1, Nr
DO J = 2, Ntet
DO K = 1, Nz
Flux(I,J,K) = -D_tet(J)*(C(I,J,K)-C(I,J-1,K)) /dtetaux(J) &
& /rc(I)
ENDDO
ENDDO
ENDDO
Flux(:,1,:) = 0.D0
Flux(:,Ntet+1,:) = 0.D0
C upstream teta-boundary
IF (BcTetup .EQ. 1) THEN
DO I = 1, Nr
DO K = 1, Nz
Flux(I,1,K) = Ftetup
ENDDO
ENDDO
ELSE IF (BcTetup .EQ. 2) THEN
DO I = 1, Nr
DO K = 1, Nz
Flux(I,1,K) = -D_tet(1)*(C(I,1,K)-Ctetup)/dtetaux(1)/rc(I)
ENDDO
ENDDO
ELSE IF (BcTetup .EQ. 3) THEN
Flux(:,1,:) = 0.D0
ELSE IF (BcTetup .EQ. 5) THEN
DO I = 1, Nr
DO K = 1, Nz
Cbound = (C(I,1,K)*dtetaux(Ntet+1) + C(I,Ntet,K)*dtetaux(1))/ &
& (dtetaux(1)+dtetaux(Ntet+1))
Flux(I,1,K) = -D_tet(1) *(C(I,1,K)-Cbound) &
& /dtetaux(1)/rc(I)
Flux(I,Ntet+1,K) = -D_tet(Ntet+1)*(Cbound-C(I,Ntet,K)) &
& /dtetaux(Ntet+1)/rc(I)
ENDDO
ENDDO
ENDIF
DO I = 1, Nr
DO K = 1, Nz
Fluxtetup(I,K) = Flux(I,1,K)
ENDDO
ENDDO
C downstream teta-boundary
IF (BcTetdown .EQ. 1) THEN
DO I = 1, Nr
DO K = 1, Nz
Flux(I,Ntet+1,K) = Ftetdown
ENDDO
ENDDO
ELSE IF (BcTetdown .EQ. 2) THEN
DO I = 1, Nr
DO K = 1, Nz
Flux(I,Ntet+1,K) = -D_tet(Ntet+1) * (Ctetdown-C(I,Ntet,K)) &
& /dtetaux(Ntet+1)/rc(I)
ENDDO
ENDDO
ELSE IF (BcTetdown .EQ. 3) THEN
Flux(:,Ntet+1,:) = 0.D0
ENDIF
DO I = 1, Nr
DO K = 1, Nz
Fluxtetdown(I,K) = Flux(I,Ntet+1,K)
ENDDO
ENDDO
C rate of change, negative flux gradient
DO I = 1,Nr
DO J = 1,Ntet
DO K = 1, Nz
IF(rc(I) .NE. 0.D0 ) &
& dC(I,J,K) = dC(I,J,K) &
& -(Flux(I,J+1,K)-Flux(I,J,K)) /dtet(J)/rc(I)
ENDDO
ENDDO
ENDDO
C ## transport in z-direction ##
C First diffusion internal cells
DO I = 1, Nr
DO J = 1, Ntet
DO K = 2, Nz
Flux(I,J,K) = -D_z(J)*(C(I,J,K)-C(I,J,K-1)) /dzaux(K)
ENDDO
ENDDO
ENDDO
Flux(:,:,1) = 0.D0
Flux(:,:,Nz+1) = 0.D0
C upstream z-boundary
IF (BCzUp .EQ. 1) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,1) = Fzup
ENDDO
ENDDO
ELSE IF (BCzUp .EQ. 2) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,1) = -D_z(1)*(C(I,J,1)-Czup) /dzaux(1)
ENDDO
ENDDO
ELSE IF (BCzUp .EQ. 3) THEN
Flux(:,:,1) = 0.D0
ELSE IF (BCzUp .EQ. 5) THEN
DO I = 1, Nr
DO J = 1, Ntet
Cbound = (C(I,J,1)*dzaux(Nz+1) + C(I,J,Nz)*dzaux(1)) / &
& (dzaux(1)+dzaux(Nz+1))
Flux(I,J,1) = -D_z(1) *(C(I,J,1)-Cbound) /dzaux(1)
Flux(I,J,Nz+1) = -D_z(Nz+1)*(Cbound-C(I,J,Nz)) /dzaux(Nz+1)
ENDDO
ENDDO
ENDIF
DO I = 1, Nr
DO J = 1, Ntet
FluxZup(I,J) = Flux(I,J,1)
ENDDO
ENDDO
C downstream z-boundary
IF (BCzDown .EQ. 1) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,Nz+1) = Fzdown
ENDDO
ENDDO
ELSE IF (BCzDown .EQ. 2) THEN
DO I = 1, Nr
DO J = 1, Ntet
Flux(I,J,Nz+1) = -D_Z(Nz+1) * (CZdown-C(I,J,Nz)) &
& /dZaux(NZ+1)
ENDDO
ENDDO
ELSE IF (BCzDown .EQ. 3) THEN
Flux(:,:,Nz+1) = 0.D0
ENDIF
DO I = 1, Nr
DO J = 1, Ntet
FluxZdown(I,J) = Flux(I,J,Nz+1)
ENDDO
ENDDO
C rate of change, negative flux gradient
DO I = 1, Nr
DO J = 1, Ntet
DO K = 1, Nz
dC(I,J,K) = dC(I,J,K) &
& -(Flux(I,J,K+1)-Flux(I,J,K))/dz(K)
ENDDO
ENDDO
ENDDO
RETURN
END
