tran3D.f

C#==============================================================================
C Diffusion in a 3-dimensional spherical coordinate system (r, theta, phi)
C all inputs are vectors/scalars
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)*D_r(1)*draux(Nr+1) +                                &
&              C(Nr,J,K)*D_r(Nr+1)*draux(1))  /                             &
&              (draux(1)*D_r(Nr+1)+draux(Nr+1)*D_r(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) THEN
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.)
ELSE
dC(I,J,K) = 0.D0
ENDIF
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
IF(rc(I) .NE. 0.D0) THEN
Flux(I,J,K) = -D_tet(J)*(C(I,J,K)-C(I,J-1,K)) /dtetaux(J)            &
&                    /rc(I)
ELSE
Flux(I,J,K) = 0.D0
ENDIF
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
IF(rc(I) .NE. 0.D0) THEN
Flux(I,1,K) = -D_tet(1)*(C(I,1,K)-Ctetup)/dtetaux(1)/rc(I)
ENDIF
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

IF(rc(I) .NE. 0.D0) THEN

Cbound = (C(I,1,K)   *D_tet(1)     *dtetaux(Ntet+1) +              &
&               C(I,Ntet,K)*D_tet(Ntet+1)*dtetaux(1)        )/              &
&              (D_tet(Ntet+1)*dtetaux(1)+D_tet(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)
ENDIF

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
IF(rc(I) .NE. 0.D0) THEN
Flux(I,Ntet+1,K) = -D_tet(Ntet+1) * (Ctetdown-C(I,Ntet,K))           &
&                                           /dtetaux(Ntet+1)/rc(I)
ENDIF
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) * 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
IF(rc(I) * sin(tet(J)) .NE. 0.D0) THEN
Flux(I,J,K) = -D_phi(K)*(C(I,J,K)-C(I,J,K-1))                      &
&                           /dphiaux(K)  /rc(I)/sin(tet(j))
ELSE
Flux(I,J,K) = 0.D0
ENDIF
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
IF (rc(I) * sin(tet(J)) .NE. 0.D0)                                   &
&        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
IF(rc(I) * sin(tet(j)) .NE. 0.D0) THEN

Cbound = (C(I,J,1)   *D_phi(1)     *dphiaux(Nphi+1)                    &
&            + C(I,J,Nphi)*D_phi(Nphi+1)*dphiaux(1))/                       &
&              (D_phi(Nphi+1)*dphiaux(1)+D_phi(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))
ENDIF
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
if (rc(I) * sin(tet(J)) .NE. 0.D0)                                   &
&      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#==============================================================================

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)*D_r(1)*draux(Nr+1) +                                &
&              C(Nr,J,K)*D_r(Nr+1)*draux(1))  /                             &
&              (draux(1)*D_r(Nr+1)+draux(Nr+1)*D_r(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)   *D_tet(1)     *dtetaux(Ntet+1) +                  &
&              C(I,Ntet,K)*D_tet(Ntet+1)*dtetaux(1))/                       &
&              (D_tet(Ntet+1)*dtetaux(1)+D_tet(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 - Karline:changed from D_Z(j) to D_Z(K) 16/04/2013

DO I = 1, Nr
DO J = 1, Ntet
DO K = 2, Nz
Flux(I,J,K) = -D_z(K)*(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 )*D_z(1)   *dzaux(Nz+1) +                            &
&              C(I,J,Nz)*D_z(Nz+1)*dzaux(1)) /                              &
&              (D_z(Nz+1)*dzaux(1)+D_z(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