Skip to main content
  • Home
  • Development
  • Documentation
  • Donate
  • Operational login
  • Browse the archive

swh logo
SoftwareHeritage
Software
Heritage
Archive
Features
  • Search

  • Downloads

  • Save code now

  • Add forge now

  • Help

https://github.com/cran/ReacTran
18 April 2025, 08:44:12 UTC
  • Code
  • Branches (32)
  • Releases (0)
  • Visits
    • Branches
    • Releases
    • HEAD
    • refs/heads/master
    • refs/tags/1.1
    • refs/tags/1.2
    • refs/tags/1.3
    • refs/tags/1.3.1
    • refs/tags/1.3.2
    • refs/tags/1.4.1
    • refs/tags/1.4.2
    • refs/tags/1.4.3.1
    • refs/tags/1.4.3.2
    • refs/tags/R-2.10.0
    • refs/tags/R-2.10.1
    • refs/tags/R-2.11.0
    • refs/tags/R-2.11.1
    • refs/tags/R-2.12.0
    • refs/tags/R-2.12.1
    • refs/tags/R-2.12.2
    • refs/tags/R-2.13.0
    • refs/tags/R-2.13.1
    • refs/tags/R-2.13.2
    • refs/tags/R-2.14.0
    • refs/tags/R-2.14.1
    • refs/tags/R-2.14.2
    • refs/tags/R-2.15.0
    • refs/tags/R-2.15.1
    • refs/tags/R-2.15.2
    • refs/tags/R-2.15.3
    • refs/tags/R-2.9.2
    • refs/tags/R-3.0.0
    • refs/tags/R-3.0.1
    • refs/tags/R-3.0.2
    • refs/tags/R-3.0.3
    No releases to show
  • dc61f12
  • /
  • src
  • /
  • tran3D.f
Raw File Download
Take a new snapshot of a software origin

If the archived software origin currently browsed is not synchronized with its upstream version (for instance when new commits have been issued), you can explicitly request Software Heritage to take a new snapshot of it.

Use the form below to proceed. Once a request has been submitted and accepted, it will be processed as soon as possible. You can then check its processing state by visiting this dedicated page.
swh spinner

Processing "take a new snapshot" request ...

Permalinks

To reference or cite the objects present in the Software Heritage archive, permalinks based on SoftWare Hash IDentifiers (SWHIDs) must be used.
Select below a type of object currently browsed in order to display its associated SWHID and permalink.

  • content
  • directory
  • revision
  • snapshot
origin badgecontent badge Iframe embedding
swh:1:cnt:be38c2a8ec99699d36d2b2361f8d2d241423a4d5
origin badgedirectory badge Iframe embedding
swh:1:dir:efc45e1b032ffb9b87436d1a9121a55fe81e000e
origin badgerevision badge
swh:1:rev:e1d00078210e32e7c6be83abec96f24e1c50b5f9
origin badgesnapshot badge
swh:1:snp:ba7fc766ddf70060bfc9bf13cb3e692cfae73f86
Citations

This interface enables to generate software citations, provided that the root directory of browsed objects contains a citation.cff or codemeta.json file.
Select below a type of object currently browsed in order to generate citations for them.

  • content
  • directory
  • revision
  • snapshot
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Tip revision: e1d00078210e32e7c6be83abec96f24e1c50b5f9 authored by Karline Soetaert on 12 June 2013, 00:00:00 UTC
version 1.4.1
Tip revision: e1d0007
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

back to top

Software Heritage — Copyright (C) 2015–2025, The Software Heritage developers. License: GNU AGPLv3+.
The source code of Software Heritage itself is available on our development forge.
The source code files archived by Software Heritage are available under their own copyright and licenses.
Terms of use: Archive access, API— Contact— JavaScript license information— Web API