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