module biovmove 

  use constants, only: zero, half, one, hnull

  implicit none
  private
  
  public :: vmove, sedimentation_resuspension, drying_case

contains
  
  subroutine vmove(kb,hi,zooa,wzoo,dt,irc)
    
    !--------------------------------------------------------------------------!
    !
    !   vmove:
    !
    !   calculate vertical movement of plankton
    !
    !   i ke                   number of layers
    !
    !   input:
    !
    !   i kb                   No of active layers at this wet surface point
    !   r hi(ke)               layer thicknesses
    !   r zooa(1:ke)           plankton, old timestep
    !   r wzoo                 plankton vertical eigenvelocity
    !   r dt                   time step
    !
    !   returns:
    !
    !
    !--------------------------------------------------------------------------!
    !
    !  Method: A conservative and stable scheme for determining the biological 
    !          redistribution which, assuming a constant vertical velocity w, 
    !          goes like this:
    !    
    !            zoo_new(k) = zoo_old(k) + (zoo_in - zoo_out)/h(k)
    !
    !          where, according to simple upwinding:
    !
    !            zoo_out = dt*abs(w)*zoo_old(k)    for all w
    !
    !          and 
    !
    !            zoo_in  = dt*w*zoo_old(k+1)       if  w >= 0
    !            zoo_in  = dt*abs(w)*zoo_old(k-1)  if  w <  0
    !
    !          Another formulation, in terms of transports through the upper
    !          and the lower cell faces, goes like this:
    !
    !            zoo_new(k) = zoo_old(k) + (zt_up + zt_lo)/h(k)
    !
    !          where
    !
    !            zt_up = -dt*w*zoo_old(k)   for w >= 0
    !            zt_lo =  dt*w*zoo_old(k+1) for w >= 0
    !
    !            zt_up =  dt*abs(w)*zoo_old(k-1) for w < 0
    !            zt_lo = -dt*abs(w)*zoo_old(k)   for w < 0
    !
    !          or:
    !
    !            zt_up = -zflux(k-1)  
    !            zt_lo =  zflux(k)
    !            zflux(k) = dt*w*zoo_old(k+1) for w >= 0
    !            zflux(k) = dt*w*zoo_old(k)   for w <  0
    !
    !--------------------------------------------------------------------------!

    implicit none

    real(8),    intent(in)    :: dt
    integer(4), intent(in)    :: kb
    integer(4), intent(inout) :: irc
    real(8),    intent(in)    :: hi(:)              ! (1:ke)
    real(8),    intent(inout) :: zooa(1:)           ! (1:ke)
    real(8),    intent(in)    :: wzoo

    real(8)    :: wdt, zoo_old, zoo_new, maxcfl
    integer(4) :: k,ndt

    ! skip trivial case:
    if (kb <= 1) return

    !--- compute the mass flow through the bottom of layer k -------------------
    wdt = wzoo*dt
    maxcfl = zero

    if (wzoo >= 0) then
      zooa(1) = zooa(1) + wdt*zooa(2)/max(hi(1),hnull)
      do k=2,kb-1
        maxcfl  = max(wdt/hi(k),maxcfl)
        zooa(k) = zooa(k) + wdt*(zooa(k+1) - zooa(k))/hi(k) 
      enddo
      maxcfl  = max(wdt/hi(kb),maxcfl)
      zooa(kb) = zooa(kb) - wdt*zooa(kb)/hi(kb)
    else ! sinking
      maxcfl  = max(-wdt/max(hi(1),hnull),maxcfl)
      zoo_old = zooa(1)
      zooa(1) = zooa(1) + wdt*zooa(1)/max(hi(1),hnull)    
      do k=2,kb-1
        maxcfl  = max(-wdt/hi(k),maxcfl)
        zoo_new = zooa(k) + wdt*(zooa(k) - zoo_old)/hi(k)
        zoo_old = zooa(k)
        zooa(k) = zoo_new
      enddo
      zooa(kb) = zooa(kb) - wdt*zoo_old/hi(kb)    
    endif

    ! did we have stability issues?
    ndt = int(maxcfl,4) + 1
    if (ndt > 1) then
      ! Yes, write an error flag'
      irc = 1
    endif
 
  end subroutine vmove

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

  subroutine sedimentation_resuspension(hbot, dt, velkb, eco, ben, kb)
    !- modules
    use bioparam, only : reslim, resuspension, resusprate,         &
                         wdetz              , &
                         idx_det            , &
                         wdetsz              , &
                         idx_dets           , &
                         idx_det_with_ship_N, &
                         idx_det_with_river_N, &
                         bidx_nitr           , &
                         bidx_sili           , &
                         bidx_nitr_with_ship_N, &
                         bidx_nitr_with_river_N, &
                         ntben
    
    implicit none

    !- args
    integer(4), intent(in)    :: kb
    real(8),    intent(in)    :: velkb, hbot, dt
    real(8),    intent(inout) :: eco(1:,0:), ben(0:)

    !- local paramters
    real(8),    parameter :: one     = 1.0_8
    real(8),    parameter :: sqrtcd0 = 0.4_8
    real(8),    parameter :: z0      = 0.0002_8

    !- local variables
    integer(4) :: i, detidx(ntben), idx
    real(8)    :: vstar, fac1, fac2

    ! set up index array:
    detidx(bidx_nitr) = idx_det
    detidx(bidx_sili) = idx_dets
    detidx(bidx_nitr_with_ship_N) = idx_det_with_ship_N
    detidx(bidx_nitr_with_river_N) = idx_det_with_river_N
    ! FIXME: when in the future ntben>2, remember to assign detidx(3:ntben)
    ! NEUMANND: SOLVED

    vstar   = sqrtcd0/abs( one+log(z0/hbot) )*velkb

    ! Use implicit euler
    if (vstar > reslim .and. resuspension) then
      ! resuspension
      fac1 = one/(one + dt*resusprate)
      fac2 = dt*resusprate/hbot
      do i=1,ntben
        idx = detidx(i)
        ben(i)      = ben(i)*fac1
        eco(kb,idx) = eco(kb,idx) + ben(i)*fac2
      enddo
    else
      ! deposition
      ! FIXME NEUMANND: wdetz is always used, remember this when adding other variables
      fac1 = hbot/(hbot - dt*wdetz)
      fac2 = dt*wdetz
      do i=1,ntben
        idx = detidx(i)
        eco(kb,idx) = eco(kb,idx)*fac1
        ben(i)      = ben(i) - eco(kb,idx)*fac2
      enddo
    endif
  end subroutine sedimentation_resuspension

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

  subroutine drying_case(nl, nu, kh, h, hbot, eco, ben, d_kpar, diag)

    !- modules
    use bioparam, only : &
                         idx_det            , &
                         idx_dets           , &
                         idx_det_with_ship_N, &
                         idx_det_with_river_N, &
                         bidx_nitr           , &
                         bidx_sili           , &
                         bidx_nitr_with_ship_N, &
                         bidx_nitr_with_river_N, &
                         ntben

    implicit none

    !- args
    integer(4), intent(in)    :: nl, nu, kh(0:)
    real(8),    intent(in)    :: hbot, h(0:)
    real(8),    intent(inout) :: eco(0:,1:), ben(0:,1:), diag(0:), d_kpar(0:)
  
    !- local variables
    integer(4) :: i, detidx(ntben), n, idx
    real(8)    :: fac
    
    ! set up index array:
    detidx(bidx_nitr) = idx_det
    detidx(bidx_sili) = idx_dets
    detidx(bidx_nitr_with_ship_N) = idx_det_with_ship_N
    detidx(bidx_nitr_with_river_N) = idx_det_with_river_N
    ! FIXME: when in the future ntben>2, remember to assign detidx(3:ntben)
    ! NEUMANND: SOLVED
    
    !- half the detritus becomes benthos (remainder goes to terrestial ecosystem)
    fac = half*hbot
    do i=1,ntben
      idx = detidx(i)
      do n=nl,nu
        if (kh(n) <= 1 .and. h(n) < hnull) then
          ben(n,i)   = ben(n,i) + eco(n,idx)*fac
          eco(n,idx) = zero
          !no light in a drying case
          diag(n)    = -one
          d_kpar(n)  = -one
        endif
      enddo
    enddo

  end subroutine drying_case

end module biovmove
