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
  • aeab7a2
  • /
  • src
  • /
  • advection.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:52388ef97ce5c042d140f29501f19c165119af8e
origin badgedirectory badge Iframe embedding
swh:1:dir:fe6372a42dcc6579e797d0a955708be5de150bbf
origin badgerevision badge
swh:1:rev:72984f188f682cc270e18126009c8e56b7aab80d
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: 72984f188f682cc270e18126009c8e56b7aab80d authored by Karline Soetaert on 10 August 2011, 00:00:00 UTC
version 1.3.2
Tip revision: 72984f1
advection.f

      subroutine advection(N, Y, dt, h, hint, v, Bcup, Bcdown,              &
     &       Yup, Ydown, VFint, VF, Aint, A, method, mode, dY, cu, it)

c-----------------------------------------------------------------------------------------
c based on the advection scheme in the GOTM model, code from 2006-11-06
c fluxes are defined on the interfaces, in an upstream-biased way.
c slope-delimeters are applied to obtain monotonic and positive schemes
c also in the presence of large gradients. 
c there are 5 different slope delimeters; first-order upstream, 
c 3rd order upstream-biased polynomial, 3rd order with superbee delimiter, 
c 3rd order with MUSCL limiter
c as described in Pietrzak 98

c Karline: made changes to make it work for negative ww...
c          added volume fraction, surface area; these two properties will generally = 1.
cc-----------------------------------------------------------------------------------------

      IMPLICIT NONE
c  number of vertical layers, time step
      INTEGER                  :: N
      DOUBLE PRECISION         :: dt

c  layer thickness (m), distance from mid to mid of each layer
      DOUBLE PRECISION         :: h(N), hint(0:N)

c  advection speed in the direction of the axis
      DOUBLE PRECISION         :: v(0:N), ww(0:N)  

c  volume fraction and surface at interface and in middle of layers
      DOUBLE PRECISION         :: VFint(0:N), Aint(0:N), VF(N), A(N)  

c  type of upper and lower Boundary Condition (only 1 and 2 used in R)
      INTEGER                  :: Bcdown, Bcup
      integer, parameter       :: Flux     = 1
      integer, parameter       :: Value    = 2
      integer, parameter       :: OneSided = 3
      integer, parameter       :: zeroDivergence = 4

c  value of upper and lower bnd conc
      DOUBLE PRECISION         :: Ydown, Yup

c  type of advection scheme, slope delimeters 
      INTEGER                  :: method
      integer, parameter       :: UPSTREAM =5  
      integer, parameter       :: P2       =4   
      integer, parameter       :: P2_PDM   =3
      integer, parameter       :: Superbee =2
      integer, parameter       :: MUSCL    =1

c  advection mode: 0= non-conservative (e.g.water flow), 1= conservative, e.g. sinking
      INTEGER                  :: mode

c  concentration to be transported
      DOUBLE PRECISION         :: Y(N)

c  rate of change due to advection
      DOUBLE PRECISION         :: dY(N)

      DOUBLE PRECISION         :: one6th=1.0d0/6.0d0
c maximal number of iterations
      INTEGER, parameter       :: itmax=10000

c  LOCAL VARIABLES:
      integer                  :: i,k,it
      DOUBLE PRECISION         :: x,r,Phi,limit=0.d0
      DOUBLE PRECISION         :: Yu,Yc,Yd
      DOUBLE PRECISION         :: c,cmax
      DOUBLE PRECISION         :: cu(0:N)
c-----------------------------------------------------------------------------------------
c  initialize upstream interface fluxes with zero
      cu   = 0.d0

c  initialize maximum Courant number
      cmax = 0.d0

c  copy of current value of state variables 
      do k =1,N
         dy(k) = y(k)
      enddo

c  compute maximum Courant number; estimate nr of iterations 
      do k=0,N
        ww(k) = v(k)*VFint(k)   ! convert to per fraction (in case VFint !=1)
        c = dabs(ww(k))*dt/hint(k)
        if (c.gt.cmax) cmax=c
      enddo

      if (cmax . GT. 1) then 
        it =  min(itmax,int(cmax)+1)   ! WAS: min(itmax,int(cmax)+1)
          it = max(1,it)
      else 
        it = 1
      endif   
c  (time) splitting loop
      do i=1,it

c     spatial loop - karline : changed into 1:N
        do k=1,N

c        positive speed 
           if (ww(k) .gt. 0.d0) then

              c=ww(k)/dble(it)*dt/hint(k-1)            ! courant number

              if (k .gt. 1) then
                 Yu=Y(k-1)                             ! upstream value
              else
                 Yu=Y(k)
              end if
              
              Yc=Y(k  )                                ! central value
              
              if (k .lt. N) then
                Yd=Y(k+1)                              ! downstream value
              else
                Yd=Y(k)
              endif
                 
              if (abs(Yd-Yc) .gt. 1e-10) then          ! slope ratio
                r=(Yc-Yu)/(Yd-Yc)
              else
                r=(Yc-Yu)*1.e10
              end if

c        negative speed
           else

             c=-ww(k)/dble(it)*dt/hint(k)              ! courant number

             if (k .lt. N) then
               Yu=Y(k+1)                               ! upstream value
             else
               Yu=Y(k)
             end if
             
             Yc=Y(k)                                   ! central value
             
             if (k .gt. 1) then
               Yd=Y(k-1)                               ! downstream value
             else
               Yd=Y(k)
             endif  

             if (abs(Yc-Yd) .gt. 1e-10) then           ! slope ratio
               r=(Yu-Yc)/(Yc-Yd)
             else
               r=(Yu-Yc)*1.e10                          
             end if

           end if

c        limit the flux according to different suggestions
           select case (method)
             case (UPSTREAM)
               limit=0.d0
             case ((P2),(P2_PDM))
c        the flux-factor phi
               x    =  one6th*(1.-2.0*c)
               Phi  =  (0.5+x)+(0.5-x)*r
               if (method.eq.P2) then
                  limit=Phi
               else
                  limit=max(0.d0,min(Phi,2./(1.d0-c),2.*r/(c+1.e-10)))
               end if
             case (Superbee)
               limit=max(0.d0, min(1.d0, 2.0*r), min(r,2.*1.d0) )
             case (MUSCL)
               limit=max(0.d0,min(2.*1.d0,2.0*r,0.5*(1.d0+r)))
             case default
c               call rerror( 'unkown advection method')  ! should not happen
           end select

c        compute the limited flux KARLINE: changed for negative ww(k) !
         if (ww(k) .gt. 0.) THEN
           cu(k)  = ww(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc))
         else
           cu(k-1)= ww(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc))
         endif           
        end do

c       downstream boundary conditions
        select case (Bcdown)
          case (flux)
            cu(N) = Ydown                ! flux OUT of the domain is positive
          case (value)
            if (ww(N).lt. 0.d0) then
              cu(N) =  ww(N)*Ydown
            else
              cu(N) =  ww(N)*Y(N)
            end if
          case (oneSided)
            if (ww(N).lt. 0.d0) then
              cu(N) =  ww(N)*Y(N)
            else
              cu(N) = 0.d0
            end if
          case (zeroDivergence)
            cu(N) = cu(N-1)
          case default
            call rwarn('unkown downstream boundary condition type')
        end select


c       upstream boundary conditions
        select case (Bcup)
          case (flux)
            cu(0) =   Yup                 ! flux into the domain is positive
          case (value)
            if(ww(0) .gt. 0.d0) then       
              cu(0) =  ww(0)*Yup
            else
              cu(0) =  ww(0)*Y(1)
            end if
          case (oneSided)
            if(ww(0) .ge. 0.d0) then
              cu(0) =  ww(0)*Y(1)
            else
              cu(0) = 0.d0
            end if
          case (zeroDivergence)
           cu(0) = cu(1)
          case default
             call rwarn('unkown upstream boundary condition type')
        end select

c     the advection step

        if (mode.eq.0) then ! non-conservative - KARLINE CHECK...
           do k=1,N
             Y(k)=Y(k)-1.d0/dble(it)*dt*((Aint(k)*cu(k)-                       &
     &                       Aint(k-1)*cu(k-1))/ h(k)/A(k)/VF(k)          &
     &          -Y(k)*(ww(k)-ww(k-1))/h(k)/A(k)/VF(k))
           enddo
        else                ! conservative
           do k=1,N
             Y(k)=Y(k)-1.d0/dble(it)*dt*((Aint(k)*cu(k)-                       &
     &                       Aint(k-1)*cu(k-1))/h(k)/A(k)/VF(k))
           enddo
        end if

      end do ! end of the iteration loop

c     rate of change due to advection 
      do k =1,N
            dy(k) = (y(k)-dy(k))/dt
          enddo

c Still to do: integrate fluxes in time (now cu = cu of last step )
c            flux = 0. ! at start
c            do k=1,N
c             flux(k)=flux(k)+1.d0/dble(it)*dt*cu(k)
 
      return
      end subroutine advection

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
! ... extensively modified by Karline Soetaert
!-----------------------------------------------------------------------






      subroutine advectvol(N, Y, dt, V, Vint, flow, Bcup, Bcdown,           &
     &                     Yup, Ydown,  method,mode,dY, cu, it)

c-----------------------------------------------------------------------------------------
c Similar as above, but for volumetric transport:
c     use flow   = v*A rather than v
c         volume = h*A rather than h
c         hint  -> c(volume(1),volume) ... for now ... not so important
c-----------------------------------------------------------------------------------------

      IMPLICIT NONE
c  number of vertical layers, time step
      INTEGER                  :: N
      DOUBLE PRECISION         :: dt

c  layer thickness (m), distance from mid to mid
      DOUBLE PRECISION         :: V(N), Vint(0:N)

c  vertical advection speed 
      DOUBLE PRECISION         :: flow(0:N)

c  type of upper and lower Boundary Condition  (1 and 2 used)
      INTEGER                  :: Bcdown, Bcup
      integer, parameter       :: Flux     =1
      integer, parameter       :: Value    =2
      integer, parameter       :: OneSided =3         ! not used
      integer, parameter       :: zeroDivergence =4   ! not used

c  value of upper and lower bnd conc
      DOUBLE PRECISION         :: Ydown, Yup

c  type of advection scheme, slope delimeters 
      INTEGER                  :: method
      integer, parameter       :: UPSTREAM =5  
      integer, parameter       :: P2       =4   
      integer, parameter       :: P2_PDM   =3
      integer, parameter       :: Superbee =2
      integer, parameter       :: MUSCL    =1

c  advection mode 0: non-conservative (e.g. water flow), 1: conservative, eg.g. sinking
      INTEGER                  :: mode

c  concentration to be transported
      DOUBLE PRECISION         :: Y(N)

c  rate of change due to advection
      DOUBLE PRECISION         :: dY(N)

      DOUBLE PRECISION         :: one6th=1.0d0/6.0d0
      INTEGER, parameter       :: itmax=100

c  LOCAL VARIABLES:
      integer                  :: i,k,it
      DOUBLE PRECISION         :: x,r,Phi,limit=0.d0
      DOUBLE PRECISION         :: Yu,Yc,Yd
      DOUBLE PRECISION         :: c,cmax
      DOUBLE PRECISION         :: cu(0:N)
c-----------------------------------------------------------------------------------------

c  initialize interface fluxes with zero
      cu   = 0.d0

c  initialize maximum Courant number
      cmax = 0.d0

c  copy of current value of state variables 
      do k =1,N
            dy(k) = y(k)
          enddo

c  compute maximum Courant number; estimate nr of iterations 
      do k=0,N
        c=dabs(flow(k))*dt/Vint(k)
        if (c.gt.cmax) cmax=c
      enddo

      it=min(itmax,int(cmax)+1)

c  (time) splitting loop
      do i=1,it

c     spatial loop - karline : changed into 1:N
        do k=1,N

c        positive speed 
           if (flow(k) .gt. 0.d0) then

              c=flow(k)/dble(it)*dt/Vint(k-1)          ! courant number

              if (k .gt. 1) then
                 Yu=Y(k-1)                             ! upstream value
              else
                 Yu=Y(k)
              end if
              
              Yc=Y(k  )                                ! central value
              
              if (k .lt. N) then
                Yd=Y(k+1)                              ! downstream value
              else
                Yd=Y(k)
              endif
                 
              if (abs(Yd-Yc) .gt. 1e-10) then          ! slope ratio
                r=(Yc-Yu)/(Yd-Yc)
              else
                r=(Yc-Yu)*1.e10
              end if

c        negative speed
           else

             c=-flow(k)/dble(it)*dt/Vint(k)            ! courant number

             if (k .lt. N) then
               Yu=Y(k+1)                               ! upstream value
             else
               Yu=Y(k)
             end if
             
             Yc=Y(k)                                   ! central value
             
             if (k .gt. 1) then
               Yd=Y(k-1)                               ! downstream value
             else
               Yd=Y(k)
             endif  

             if (dabs(Yc-Yd) .gt. 1e-10) then           ! slope ratio
               r=(Yu-Yc)/(Yc-Yd)
             else
               r=(Yu-Yc)*1.e10                         ! CHECK THIS: not -???
             end if

           end if

c        limit the flux according to different suggestions, phi = flux-factor
           select case (method)
             case (UPSTREAM)
               limit=0.d0
             case ((P2),(P2_PDM))
c         - for quickest
               x    =  one6th*(1.d0-2.d0*c)
               Phi  =  (0.5d0+x)+(0.5d0-x)*r

               if (method.eq.P2) then
                  limit=Phi
               else
                  limit=max(0.d0,min(Phi,2./(1.-c),2.*r/(c+1.e-10)))
               end if
             case (Superbee)
               limit=max(0.d0, min(1.d0, 2.0*r), min(r,2.*1.d0) )
             case (MUSCL)
               limit=max(0.d0,min(2.*1.d0,2.0*r,0.5*(1.0+r)))
             case default
c               call rerror( 'unkown advection method')
           end select

c        compute the limited flux KARLINE: changed for negative flow(k) !
         if (flow(k) .gt. 0.) THEN
           cu(k)  = flow(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc))
         else
           cu(k-1)= flow(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc))
         endif           
        end do

c       downstream boundary conditions
        select case (Bcdown)
          case (flux)
            cu(N) = Ydown                  ! flux OUT of the domain is positive
          case (value)
            if (flow(N).lt. 0.d0) then
              cu(N) =  flow(N)*Ydown
            else
              cu(N) =  flow(N)*Y(N)
            end if
          case (oneSided)
            if (flow(N).lt. 0.d0) then
              cu(N) =  flow(N)*Y(N)
            else
              cu(N) = 0.d0
            end if
          case (zeroDivergence)
            cu(N) = cu(N-1)
          case default
            call rwarn('unkown downstream boundary condition type')
        end select


c       upstream boundary conditions
        select case (Bcup)
          case (flux)
            cu(0) =   Yup                   ! flux into the domain is positive
          case (value)
            if(flow(0) .gt. 0.d0) then      ! Karline: CHECK!
              cu(0) =  flow(0)*Yup
            else
              cu(0) =  flow(0)*Y(1)
            end if
          case (oneSided)
            if(flow(0) .ge. 0.d0) then
              cu(0) =  flow(0)*Y(1)
            else
              cu(0) = 0.d0
            end if
          case (zeroDivergence)
           cu(0) = cu(1)
          case default
             call rwarn('unkown upstream boundary condition type')
        end select

c     the advection step

        if (mode.eq.0) then ! non-conservative - KARLINE CHECK...
           do k=1,N
             Y(k)=Y(k)-1.d0/dble(it)*dt*((cu(k)- cu(k-1))/ V(k)                &
     &          -Y(k)*(flow(k)-flow(k-1))/V(k))
           enddo
        else                ! conservative - this is actually used
           do k=1,N
             Y(k)=Y(k)-1.d0/dble(it)*dt*((cu(k)- cu(k-1))/V(k))
           enddo
        end if

      end do ! end of the iteration loop

c     rate of change due to advection 
      do k =1,N
            dy(k) = (y(k)-dy(k))/dt
          enddo

c Still to do: integrate fluxes in time (now cu = cu of last step )
c            flux = 0. ! at start
c            do k=1,N
c             flux(k)=flux(k)+1.d0/real(it)*dt*cu(k)
 
      return
      end subroutine advectvol

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

back to top