anchor.diff:            Qout_max(n,k) = ( q(n,k,l) - CQin_min - qnext_min*(1.D0-Cin-Cout+d(n,k,l)) ) &
                           / ( Cout + zerosw ) * ( 1.D0 - zerosw )                        &
                           + q(n,k,l) * zerosw
  -
  -       enddo
  -       enddo
  -
  -       do k = ADM_kmin, ADM_kmax+1
  +       enddo
  +       enddo
  +
  +       do k = ADM_kmin+1, ADM_kmax
          do n = 1, ADM_gall
  -          mask = 0.5D0 - sign(0.5D0,ck(n,k,l,1)-CNST_EPS_ZERO)
  -
  -          q_h(n,k,l) = (      mask ) * min( max( q_h(n,k,l), Qin_min(n,k  ,1) ), Qin_max(n,k  ,1) ) &
  -                     + ( 1.D0-mask ) * min( max( q_h(n,k,l), Qin_min(n,k-1,2) ), Qin_max(n,k-1,2) )
  -
  -          q_h(n,k,l) = (      mask ) * max( min( q_h(n,k,l), Qout_max(n,k-1) ), Qout_min(n,k-1) ) &
  -                     + ( 1.D0-mask ) * max( min( q_h(n,k,l), Qout_max(n,k  ) ), Qout_min(n,k  ) )
  +          inflagL = 0.5D0 - sign(0.5D0,ck(n,k,l,1)) ! incoming flux: flag=1
  +
  +          q_h(n,k,l) = (      inflagL ) * max( min( q_h(n,k,l), Qout_max(n,k-1) ), Qout_min(n,k-1) ) &
  +                     + ( 1.D0-inflagL ) * max( min( q_h(n,k,l), Qout_max(n,k  ) ), Qout_min(n,k  ) )
          enddo
          enddo
   
  
positive.code:     real(8), intent(in)    :: ck    (ADM_gall,   ADM_kall,ADM_lall   ,2)
      real(8), intent(in)    :: ck_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl,2)
  
      real(8) :: Qin_min    (ADM_gall,   ADM_kall,2)
      real(8) :: Qin_min_pl (ADM_gall_pl,ADM_kall,2)
      real(8) :: Qin_max    (ADM_gall,   ADM_kall,2)
      real(8) :: Qin_max_pl (ADM_gall_pl,ADM_kall,2)
  
      real(8) :: Qout_min   (ADM_gall,   ADM_kall)
      real(8) :: Qout_min_pl(ADM_gall_pl,ADM_kall)
      real(8) :: Qout_max   (ADM_gall,   ADM_kall)
      real(8) :: Qout_max_pl(ADM_gall_pl,ADM_kall)
  
      real(8) :: qnext_min, qnext_min_pl
      real(8) :: qnext_max, qnext_max_pl
      real(8) :: Cin,       Cin_pl
      real(8) :: Cout,      Cout_pl
      real(8) :: CQin_min,  CQin_min_pl
      real(8) :: CQin_max,  CQin_max_pl
  
      real(8) :: mask, mask1, mask2
      real(8) :: zerosw
  
      integer :: n, k, l
  
negative.code:        ADM_kmin,    &
         ADM_kmax
      use mod_cnst, only: &
         BIG => CNST_MAX_REAL, &
         EPS => CNST_EPS_ZERO
      implicit none
  
      real(8), intent(inout) :: q_h   (ADM_gall,   ADM_kall,ADM_lall   )
  