module biosolver
  !=============================================================================
  ! Module containing timestepping schemes for the biogeochemical model
  !=============================================================================

  use constants, only : zero, half, one, ten, tinyv

  implicit none

  ! public routines
  public :: solvebio

  ! private routines
  private :: empatankar, bisec, f

contains
  subroutine solvebio(eco, deriv, dt, irc, k, ben, benderiv)
    use bioparam, only : ntbio_notneg, ntben_notneg, ntbio, ntben, ntall, &
                         idx_amm            ,                                 &
                         idx_nit            ,                                 &
                         idx_phos           ,                                 &
                         idx_sil            ,                                 &
                         idx_dia            ,                                 &
                         idx_flag           ,                                 &
                         idx_cyano          ,                                 &
                         idx_mez            ,                                 &
                         idx_miz            ,                                 &
                         idx_det            ,                                 &
                         idx_dets           ,                                 &
                         idx_ldon           ,                                 &
                         idx_oxy            ,                                 &
                         idx_amm_with_ship_N,                                 &
                         idx_nit_with_ship_N,                                 &
                         idx_dia_with_ship_N,                                 &
                         idx_flag_with_ship_N,                                 &
                         idx_cyano_with_ship_N,                                 &
                         idx_mez_with_ship_N,                                 &
                         idx_miz_with_ship_N,                                 &
                         idx_det_with_ship_N,                                 &
                         idx_ldon_with_ship_N,                                 &
                         idx_amm_with_river_N,                                 &
                         idx_nit_with_river_N,                                 &
                         idx_dia_with_river_N,                                 &
                         idx_flag_with_river_N,                                 &
                         idx_cyano_with_river_N,                                 &
                         idx_mez_with_river_N,                                 &
                         idx_miz_with_river_N,                                 &
                         idx_det_with_river_N,                                 &
                         idx_ldon_with_river_N,                                 &
                         bidx_nitr           ,                                &
                         bidx_sili           ,                                &
                         bidx_nitr_with_ship_N,                                &
                         bidx_nitr_with_river_N,                                &
                         impl, impl_cgt, idcs_all_notneg
    use dmi_mpi_global, only : iu06
    implicit none

    ! arguments
    integer(4), intent(inout)           :: irc
    integer(4), intent(in)              :: k
    real(8),    intent(in)              :: deriv(0:), dt
    real(8),    intent(inout)           :: eco(:,0:)
    real(8),    intent(inout), optional :: ben(0:)
    real(8),    intent(in),    optional :: benderiv(0:)
       
    ! locals
    integer(4) :: itr, numimpltr, instablei
    real(8)    :: ecotmp(ntbio_notneg+ntben_notneg), derivtmp(ntbio_notneg+ntben_notneg) 
    real(8)    :: ecotmp_cgt(ntall)
    real(8)    :: instable,ecot,deco,bent,dben
    logical    :: bottom

    ! parameters
    real(8), parameter :: eps  = 1.0e-14_8

    if (present(benderiv)) then
      bottom = .true.
    else
      bottom = .false.
    endif

    if (impl) then
      if (impl_cgt) then
        !----------------------------------------------------------------------
        ! cgt Implicit scheme
        !----------------------------------------------------------------------
        ! Difference to the standard implicite schema:
        !  - we pass all tracers to the (modified-modified) patankar schema
        !  - in the modified-modified patankar schema al dt*deriv are scaled 
        !     by p but p is only calculated on the basis of the have-to-be-
        !     positive tracer
        !  - the first part of the previous bullet point is acutal change
        if ( bottom ) then
          ecotmp_cgt = [eco(k,1:ntbio),ben(1:ntben)]
          call empatankar_cgt_ntall(ecotmp_cgt, [deriv(1:ntbio),benderiv(1:ntben)],ntall, dt,irc) ! update tracers
          eco(k,1:ntbio) = ecotmp_cgt(1:ntbio)
          ben(1:ntben) = ecotmp_cgt(ntbio+1:ntbio+ntben)
        else
          call empatankar_cgt_ntbio(eco(k,1:ntbio), deriv,ntbio, dt,irc) ! update tracers
        endif
        ! the [ ... ] to concatenate arrays is new in Fortran 2003, see
        ! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/292553
        
        if (irc /= 0) return
        
      else
        !-------------------------------------------------------------------------
        ! standard Implicit scheme
        !-------------------------------------------------------------------------
        ! We consider only tracers that have to remain positive. The array
        ! idcs_all_notneg gives us the mapping. In the un-tagged case only 
        ! oxygen is affected.
        DO itr=1,ntbio_notneg
          ecotmp(itr)   = eco(k,idcs_all_notneg(itr))
          derivtmp(itr) = deriv(idcs_all_notneg(itr))
        ENDDO
        if (bottom) then
          ! solution fot bottom point
          DO itr=ntbio_notneg+1,ntbio_notneg+ntben_notneg
            ecotmp(itr)   = ben(idcs_all_notneg(itr))
            derivtmp(itr) = benderiv(idcs_all_notneg(itr))
          enddo
          numimpltr = ntbio_notneg+ntben_notneg ! eco + benthos - maybe-negative
        else 
          ! solution of pelagic
          ecotmp(ntbio_notneg+1:)       = zero 
          derivtmp(ntbio_notneg+1:)     = zero
          numimpltr = ntbio_notneg ! eco - maybe-negative 
        endif

        call empatankar(ecotmp, derivtmp, numimpltr, dt,irc) ! update tracers
        if (irc /= 0) return

        ! copy values back to eco and ben array
        DO itr=1,ntbio_notneg
          eco(k,idcs_all_notneg(itr))   = ecotmp(itr)
        ENDDO

        ! update non-have-to-be-positive tracers (i.e. oxygen) with an euler 
        ! step since they can be negative.
        eco(k,idx_oxy) = eco(k,idx_oxy) + deriv(idx_oxy)*dt
        eco(k,idx_amm_with_ship_N) = eco(k,idx_amm_with_ship_N) + deriv(idx_amm_with_ship_N)*dt
        eco(k,idx_nit_with_ship_N) = eco(k,idx_nit_with_ship_N) + deriv(idx_nit_with_ship_N)*dt
        eco(k,idx_dia_with_ship_N) = eco(k,idx_dia_with_ship_N) + deriv(idx_dia_with_ship_N)*dt
        eco(k,idx_flag_with_ship_N) = eco(k,idx_flag_with_ship_N) + deriv(idx_flag_with_ship_N)*dt
        eco(k,idx_cyano_with_ship_N) = eco(k,idx_cyano_with_ship_N) + deriv(idx_cyano_with_ship_N)*dt
        eco(k,idx_mez_with_ship_N) = eco(k,idx_mez_with_ship_N) + deriv(idx_mez_with_ship_N)*dt
        eco(k,idx_miz_with_ship_N) = eco(k,idx_miz_with_ship_N) + deriv(idx_miz_with_ship_N)*dt
        eco(k,idx_det_with_ship_N) = eco(k,idx_det_with_ship_N) + deriv(idx_det_with_ship_N)*dt
        eco(k,idx_ldon_with_ship_N) = eco(k,idx_ldon_with_ship_N) + deriv(idx_ldon_with_ship_N)*dt
        eco(k,idx_amm_with_river_N) = eco(k,idx_amm_with_river_N) + deriv(idx_amm_with_river_N)*dt
        eco(k,idx_nit_with_river_N) = eco(k,idx_nit_with_river_N) + deriv(idx_nit_with_river_N)*dt
        eco(k,idx_dia_with_river_N) = eco(k,idx_dia_with_river_N) + deriv(idx_dia_with_river_N)*dt
        eco(k,idx_flag_with_river_N) = eco(k,idx_flag_with_river_N) + deriv(idx_flag_with_river_N)*dt
        eco(k,idx_cyano_with_river_N) = eco(k,idx_cyano_with_river_N) + deriv(idx_cyano_with_river_N)*dt
        eco(k,idx_mez_with_river_N) = eco(k,idx_mez_with_river_N) + deriv(idx_mez_with_river_N)*dt
        eco(k,idx_miz_with_river_N) = eco(k,idx_miz_with_river_N) + deriv(idx_miz_with_river_N)*dt
        eco(k,idx_det_with_river_N) = eco(k,idx_det_with_river_N) + deriv(idx_det_with_river_N)*dt
        eco(k,idx_ldon_with_river_N) = eco(k,idx_ldon_with_river_N) + deriv(idx_ldon_with_river_N)*dt
        
        if (bottom) then
          DO itr=ntbio_notneg+1,ntbio_notneg+ntben_notneg
            ben(idcs_all_notneg(itr)) = ecotmp(itr)
          enddo

          ! update non-have-to-be-positive tracers (i.e. oxygen) with an euler 
          ! step since they can be negative.
          ben(bidx_nitr_with_ship_N) = ben(bidx_nitr_with_ship_N) + benderiv(bidx_nitr_with_ship_N)*dt
          ben(bidx_nitr_with_river_N) = ben(bidx_nitr_with_river_N) + benderiv(bidx_nitr_with_river_N)*dt
        endif
      
      endif
      
      ! should we have a check here as well?
      instablei = 0

    else 
      !-------------------------------------------------------------------------
      ! Explicit scheme
      !-------------------------------------------------------------------------
      instable = zero
      do itr=1,ntbio_notneg
        ecot = eco(k,idcs_all_notneg(itr))
        deco = deriv(idcs_all_notneg(itr))*dt
        eco(k,idcs_all_notneg(itr)) = ecot + deco
        instable = max(instable,-deco/(ecot+eps))
      enddo

      ! update non-have-to-be-positive tracers (i.e. oxygen) with an euler 
      ! step since they can be negative.
      eco(k,idx_oxy) = eco(k,idx_oxy) + deriv(idx_oxy)*dt
      eco(k,idx_amm_with_ship_N) = eco(k,idx_amm_with_ship_N) + deriv(idx_amm_with_ship_N)*dt
      eco(k,idx_nit_with_ship_N) = eco(k,idx_nit_with_ship_N) + deriv(idx_nit_with_ship_N)*dt
      eco(k,idx_dia_with_ship_N) = eco(k,idx_dia_with_ship_N) + deriv(idx_dia_with_ship_N)*dt
      eco(k,idx_flag_with_ship_N) = eco(k,idx_flag_with_ship_N) + deriv(idx_flag_with_ship_N)*dt
      eco(k,idx_cyano_with_ship_N) = eco(k,idx_cyano_with_ship_N) + deriv(idx_cyano_with_ship_N)*dt
      eco(k,idx_mez_with_ship_N) = eco(k,idx_mez_with_ship_N) + deriv(idx_mez_with_ship_N)*dt
      eco(k,idx_miz_with_ship_N) = eco(k,idx_miz_with_ship_N) + deriv(idx_miz_with_ship_N)*dt
      eco(k,idx_det_with_ship_N) = eco(k,idx_det_with_ship_N) + deriv(idx_det_with_ship_N)*dt
      eco(k,idx_ldon_with_ship_N) = eco(k,idx_ldon_with_ship_N) + deriv(idx_ldon_with_ship_N)*dt
      eco(k,idx_amm_with_river_N) = eco(k,idx_amm_with_river_N) + deriv(idx_amm_with_river_N)*dt
      eco(k,idx_nit_with_river_N) = eco(k,idx_nit_with_river_N) + deriv(idx_nit_with_river_N)*dt
      eco(k,idx_dia_with_river_N) = eco(k,idx_dia_with_river_N) + deriv(idx_dia_with_river_N)*dt
      eco(k,idx_flag_with_river_N) = eco(k,idx_flag_with_river_N) + deriv(idx_flag_with_river_N)*dt
      eco(k,idx_cyano_with_river_N) = eco(k,idx_cyano_with_river_N) + deriv(idx_cyano_with_river_N)*dt
      eco(k,idx_mez_with_river_N) = eco(k,idx_mez_with_river_N) + deriv(idx_mez_with_river_N)*dt
      eco(k,idx_miz_with_river_N) = eco(k,idx_miz_with_river_N) + deriv(idx_miz_with_river_N)*dt
      eco(k,idx_det_with_river_N) = eco(k,idx_det_with_river_N) + deriv(idx_det_with_river_N)*dt
      eco(k,idx_ldon_with_river_N) = eco(k,idx_ldon_with_river_N) + deriv(idx_ldon_with_river_N)*dt
      
      !update benthos
      if (bottom) then
        do itr=ntbio_notneg+1,ntbio_notneg+ntben_notneg
          bent = ben(idcs_all_notneg(itr))
          dben = benderiv(idcs_all_notneg(itr))*dt
          ben(idcs_all_notneg(itr)) = bent + dben
          instable = max(instable,-dben/(bent+eps))
        enddo

        ! update non-have-to-be-positive sediment tracers with an euler 
        ! step since they can be negative.
        ben(bidx_nitr_with_ship_N) = ben(bidx_nitr_with_ship_N) + benderiv(bidx_nitr_with_ship_N)*dt
        ben(bidx_nitr_with_river_N) = ben(bidx_nitr_with_river_N) + benderiv(bidx_nitr_with_river_N)*dt
      endif


      instablei = int(instable,4) + 1
    endif
 
    if (instablei > 1) then
      ! we will get negative values
      irc = 3
    endif
    
  end subroutine solvebio

  !=============================================================================

  ! 2017-02-21: two major changes by D. Neumann
  !       - removed pmaxfirst and simplified the pmax
  !          setting => same result, less calculations
  !       - removed "if a(itr) close to zero" statement when one is divided by
  !          a(itr). Now, eps is always added to a(itr) when we divide by it.
  !          This has the largest impact, when a(itr) is close to one. Tests 
  !          showed that pdum is reduced in the order of 10^-14 which is 
  !          negligible. However, results might slightly change if the actual
  !          p was between -one/a(itr) and -one/(a(itr)-eps) (which is quite
  !          improbable).
  !       - both changes together reduce the run time by 40 to 50 CPU seconds.
  subroutine empatankar_cgt(tracer, deriv, numtr, dt, irc)
    use bioparam, only :                             &
                         idx_amm            ,             &
                         idx_nit            ,             &
                         idx_phos           ,             &
                         idx_sil            ,             &
                         idx_dia            ,             &
                         idx_flag           ,             &
                         idx_cyano          ,             &
                         idx_mez            ,             &
                         idx_miz            ,             &
                         idx_det            ,             &
                         idx_dets           ,             &
                         idx_ldon           ,             &
                         idx_oxy            ,             &
                         idx_amm_with_ship_N,             &
                         idx_nit_with_ship_N,             &
                         idx_dia_with_ship_N,             &
                         idx_flag_with_ship_N,             &
                         idx_cyano_with_ship_N,             &
                         idx_mez_with_ship_N,             &
                         idx_miz_with_ship_N,             &
                         idx_det_with_ship_N,             &
                         idx_ldon_with_ship_N,             &
                         idx_amm_with_river_N,             &
                         idx_nit_with_river_N,             &
                         idx_dia_with_river_N,             &
                         idx_flag_with_river_N,             &
                         idx_cyano_with_river_N,             &
                         idx_mez_with_river_N,             &
                         idx_miz_with_river_N,             &
                         idx_det_with_river_N,             &
                         idx_ldon_with_river_N,             &
                         notnegall
    implicit none

    ! args
    integer(4), intent(inout) :: irc
    integer(4), intent(in)    :: numtr
    real(8), intent(in)       :: dt,deriv(:)
    real(8), intent(inout)    :: tracer(:)

    ! local
    real(8)    :: pmin,pmax,a(numtr),p,pdt
    integer(4) :: numnegdt,indnegdt(numtr), itr

    integer(4), parameter :: maxit = 50
    real(8), parameter    :: tol   = 0.00001_8
    real(8), parameter    :: eps   = 1.0e-14_8
    
    numnegdt    = 0
    indnegdt(:) = 0
    a(:)        = -eps

    do itr=1,numtr
      if (deriv(itr) < zero .and. notnegall(itr)) then
        numnegdt = numnegdt + 1
        indnegdt(numnegdt) = itr

        ! find pmax
        a(itr) = deriv(itr)*dt/(tracer(itr)+eps)
      endif
    enddo

    pmax = min(-one/minval(a),one)
    pmin = zero

    if (numnegdt > 0) then
      call bisec(a, numnegdt,indnegdt, pmin, pmax,maxit,tol, p,irc )
      if (irc /= 0) return
      pdt = p * dt
    else
      ! p = one
      pdt = dt
    endif

    ! update tracers
    tracer(:) = tracer(:) + deriv(:)*pdt
    !~ do itr=1,numtr
      !~ if( notnegall(itr) ) then
        !~ tracer(itr) = tracer(itr) + deriv(itr)*p*dt
      !~ else
        !~ tracer(itr) = tracer(itr) + deriv(itr)*dt
      !~ endif
    !~ enddo
  
  end subroutine empatankar_cgt

  !=============================================================================

  ! 2017-02-21: two major changes by D. Neumann
  !       - removed pmaxfirst and simplified the pmax
  !          setting => same result, less calculations
  !       - removed "if a(itr) close to zero" statement when one is divided by
  !          a(itr). Now, eps is always added to a(itr) when we divide by it.
  !          This has the largest impact, when a(itr) is close to one. Tests 
  !          showed that pdum is reduced in the order of 10^-14 which is 
  !          negligible. However, results might slightly change if the actual
  !          p was between -one/a(itr) and -one/(a(itr)-eps) (which is quite
  !          improbable).
  !       - both changes together reduce the run time by 40 to 50 CPU seconds.
  subroutine empatankar_cgt_ntbio(tracer, deriv, numtr, dt, irc)
    use exits,    only : exitme
    use bioparam, only : ntbio,                            &
                         idx_amm            ,             &
                         idx_nit            ,             &
                         idx_phos           ,             &
                         idx_sil            ,             &
                         idx_dia            ,             &
                         idx_flag           ,             &
                         idx_cyano          ,             &
                         idx_mez            ,             &
                         idx_miz            ,             &
                         idx_det            ,             &
                         idx_dets           ,             &
                         idx_ldon           ,             &
                         idx_oxy            ,             &
                         idx_amm_with_ship_N,             &
                         idx_nit_with_ship_N,             &
                         idx_dia_with_ship_N,             &
                         idx_flag_with_ship_N,             &
                         idx_cyano_with_ship_N,             &
                         idx_mez_with_ship_N,             &
                         idx_miz_with_ship_N,             &
                         idx_det_with_ship_N,             &
                         idx_ldon_with_ship_N,             &
                         idx_amm_with_river_N,             &
                         idx_nit_with_river_N,             &
                         idx_dia_with_river_N,             &
                         idx_flag_with_river_N,             &
                         idx_cyano_with_river_N,             &
                         idx_mez_with_river_N,             &
                         idx_miz_with_river_N,             &
                         idx_det_with_river_N,             &
                         idx_ldon_with_river_N,             &
                         notnegall
    implicit none

    ! args
    integer(4), intent(inout) :: irc
    integer(4), intent(in)    :: numtr
    real(8), intent(in)       :: dt,deriv(1:ntbio)
    real(8), intent(inout)    :: tracer(1:ntbio)

    ! local
    real(8)    :: pmin,pmax,a(ntbio),p,pdt
    integer(4) :: numnegdt,indnegdt(ntbio), itr

    integer(4), parameter :: maxit = 50
    real(8), parameter    :: tol   = 0.00001_8
    real(8), parameter    :: eps   = 1.0e-14_8
    
    if ( numtr .ne. ntbio ) then
      call exitme(1, 'ERROR: empatankar_cgt_ntbio called with wrong number of tracers NUMTR.')
    endif
    
    numnegdt    = 0
    indnegdt(:) = 0
    a(:)        = -eps

    do itr=1,ntbio
      if (deriv(itr) < zero .and. notnegall(itr)) then
        numnegdt = numnegdt + 1
        indnegdt(numnegdt) = itr

        ! find pmax
        a(itr) = deriv(itr)*dt/(tracer(itr)+eps)
      endif
    enddo

    pmax = min(-one/minval(a),one)
    pmin = zero

    if (numnegdt > 0) then
      call bisec(a, numnegdt,indnegdt, pmin, pmax,maxit,tol, p,irc )
      if (irc /= 0) return
      pdt = p * dt
    else
      ! p = one
      pdt = dt
    endif

    ! update tracers
    tracer(:) = tracer(:) + deriv(:) * pdt
    !~ do itr=1,numtr
      !~ if( notnegall(itr) ) then
        !~ tracer(itr) = tracer(itr) + deriv(itr)*p*dt
      !~ else
        !~ tracer(itr) = tracer(itr) + deriv(itr)*dt
      !~ endif
    !~ enddo
  
  end subroutine empatankar_cgt_ntbio

  !=============================================================================

  ! 2017-02-21: two major changes by D. Neumann
  !       - removed pmaxfirst and simplified the pmax
  !          setting => same result, less calculations
  !       - removed "if a(itr) close to zero" statement when one is divided by
  !          a(itr). Now, eps is always added to a(itr) when we divide by it.
  !          This has the largest impact, when a(itr) is close to one. Tests 
  !          showed that pdum is reduced in the order of 10^-14 which is 
  !          negligible. However, results might slightly change if the actual
  !          p was between -one/a(itr) and -one/(a(itr)-eps) (which is quite
  !          improbable).
  !       - both changes together reduce the run time by 40 to 50 CPU seconds.
  subroutine empatankar_cgt_ntall(tracer, deriv, numtr, dt, irc)
    use exits,    only : exitme
    use bioparam, only : ntall,                            &
                         idx_amm            ,             &
                         idx_nit            ,             &
                         idx_phos           ,             &
                         idx_sil            ,             &
                         idx_dia            ,             &
                         idx_flag           ,             &
                         idx_cyano          ,             &
                         idx_mez            ,             &
                         idx_miz            ,             &
                         idx_det            ,             &
                         idx_dets           ,             &
                         idx_ldon           ,             &
                         idx_oxy            ,             &
                         idx_amm_with_ship_N,             &
                         idx_nit_with_ship_N,             &
                         idx_dia_with_ship_N,             &
                         idx_flag_with_ship_N,             &
                         idx_cyano_with_ship_N,             &
                         idx_mez_with_ship_N,             &
                         idx_miz_with_ship_N,             &
                         idx_det_with_ship_N,             &
                         idx_ldon_with_ship_N,             &
                         idx_amm_with_river_N,             &
                         idx_nit_with_river_N,             &
                         idx_dia_with_river_N,             &
                         idx_flag_with_river_N,             &
                         idx_cyano_with_river_N,             &
                         idx_mez_with_river_N,             &
                         idx_miz_with_river_N,             &
                         idx_det_with_river_N,             &
                         idx_ldon_with_river_N,             &
                         notnegall
    implicit none

    ! args
    integer(4), intent(inout) :: irc
    integer(4), intent(in)    :: numtr
    real(8), intent(in)       :: dt,deriv(1:ntall)
    real(8), intent(inout)    :: tracer(1:ntall)

    ! local
    real(8)    :: pmin,pmax,a(ntall),p,pdt
    integer(4) :: numnegdt,indnegdt(ntall), itr

    integer(4), parameter :: maxit = 50
    real(8), parameter    :: tol   = 0.00001_8
    real(8), parameter    :: eps   = 1.0e-14_8
    
    if ( numtr .ne. ntall ) then
      call exitme(1, 'ERROR: empatankar_cgt_ntall called with wrong number of tracers NUMTR.')
    endif
    
    numnegdt    = 0
    indnegdt(:) = 0
    a(:)        = -eps

    do itr=1,ntall
      if (deriv(itr) < zero .and. notnegall(itr)) then
        numnegdt = numnegdt + 1
        indnegdt(numnegdt) = itr

        ! find pmax
        a(itr) = deriv(itr)*dt/(tracer(itr)+eps)
      endif
    enddo

    pmax = min(-one/minval(a),one)
    pmin = zero

    if (numnegdt > 0) then
      call bisec(a, numnegdt,indnegdt, pmin, pmax,maxit,tol, p,irc )
      if (irc /= 0) return
      pdt = p * dt
    else
      ! p = one
      pdt = dt
    endif

    ! update tracers
    tracer(:) = tracer(:) + deriv(:)*pdt
    !~ do itr=1,numtr
      !~ if( notnegall(itr) ) then
        !~ tracer(itr) = tracer(itr) + deriv(itr)*p*dt
      !~ else
        !~ tracer(itr) = tracer(itr) + deriv(itr)*dt
      !~ endif
    !~ enddo
  
  end subroutine empatankar_cgt_ntall

  !=============================================================================

  subroutine empatankar(tracer, deriv, numtr, dt, irc)
    implicit none

    ! args
    integer(4), intent(inout) :: irc
    integer(4), intent(in)    :: numtr
    real(8), intent(in)       :: dt,deriv(numtr)
    real(8), intent(inout)    :: tracer(numtr)

    ! local
    real(8)    :: pmin,pmax,a(numtr),p,pdt
    integer(4) :: numnegdt,indnegdt(numtr), itr

    integer(4), parameter :: maxit = 50
    real(8), parameter    :: tol   = 0.00001_8
    real(8), parameter    :: eps   = 1.0e-14_8

    numnegdt    = 0
    indnegdt(:) = 0
    a(:)        = -eps

    do itr=1,numtr
      if (deriv(itr) < zero) then
        numnegdt = numnegdt + 1
        indnegdt(numnegdt) = itr

        ! find pmax
        a(itr) = deriv(itr)*dt/(tracer(itr)+eps)
      endif
    enddo

    pmax = min(-one/minval(a),one)
    pmin = zero

    if (numnegdt > 0) then
      call bisec(a, numnegdt,indnegdt, pmin, pmax,maxit,tol, p,irc )
      if (irc /= 0) return
      pdt = p * dt
    else
      pdt = dt
    endif

    ! update tracers
    do itr=1,numtr
      tracer(itr) = tracer(itr) + deriv(itr)*pdt
    enddo
  
  end subroutine empatankar

  !============================================================================
  subroutine bisec(a,numnegdt,indnegdt,x1,x2,maxit,eps_step,c,irc )
    !--------
    ! simple bisection routine for finding the root of function f in the 
    ! low-high p-space
    ! remember dg/dp is always negative, thus we can simplify the bisection
    !--------
    
    implicit none

    !args
    real(8), intent(in)       :: a(:), eps_step
    real(8), intent(inout)    :: c, x1, x2
    integer(4), intent(in)    :: numnegdt, indnegdt(:), maxit
    integer(4), intent(inout) :: irc

    ! local
    real(8)    :: rtb,fc, ff,dx
    integer(4) :: k

    ! Usually the root is located close to the upper bound, make a quick look
    ! around this point
    fc = f(x2,a,numnegdt,indnegdt)
    ff = f(x2-ten*eps_step,a,numnegdt,indnegdt)
    if ( ff > zero) then
      x1 = x2-ten*eps_step
    else
      ff = f(x2-ten*ten*eps_step,a,numnegdt,indnegdt)
      if (ff > zero) then
        x1 = x2-ten*ten*eps_step
      else
        ff = f(x1,a,numnegdt,indnegdt)
      endif
    endif
    
    ! Does the bracketing contains the root?
    if (ff*fc >= zero) then
      ! Root is not bracketed'
      irc = 4
      return
    endif

    ! we know the orientation already since the dg/dp is negative
    dx  = x1-x2
    rtb = x2
    do k=1,maxit
      dx = dx*half
      c  = rtb+dx
      fc = f(c,a,numnegdt,indnegdt)

      if (fc < zero) rtb = c
!FIXME: Do we need to check for "fc == zero"?
!       This will probably never be the case.
      if (abs(dx) < eps_step .or. fc == zero) return
    enddo

    !Implicit bioscheme did not converge'
    irc = 5
  end subroutine bisec

  !============================================================================

  real(8) function f(p,a,numnegdt,indnegdt)
    implicit none

    ! args
    real(8), intent(in)    :: p, a(:)
    integer(4), intent(in) :: numnegdt, indnegdt(:)

    ! local
    integer(4)        :: ii, i
    real(8)           :: y

    ii = indnegdt(1)
    y = (one+a(ii)*p)
    do i=2,numnegdt
      ii = indnegdt(i)
      y = y * (one+a(ii)*p)
    enddo
    y = y - p

    f = y
  end function f

  !============================================================================

end module biosolver
