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:        enddo
  
         do k = ADM_kmin, ADM_kmax
            rhog(suf(ADM_gall_1d,1),k,l) = rhog(suf(ADM_gmax+1,ADM_gmin),k,l)
            rhog(suf(1,ADM_gall_1d),k,l) = rhog(suf(ADM_gmin,ADM_gmax+1),k,l)
         enddo
         rhog(:,ADM_kmin-1,l) = rhog_in(:,ADM_kmin,l)
         rhog(:,ADM_kmax+1,l) = rhog_in(:,ADM_kmax,l)
  
negative.code:           ijm1   = n     - ADM_gall_1d
            ip1j   = n + 1
  
            flux = 0.5D0 * ( (rhovxt(ijm1,TJ)+rhovxt(ij,TI)) * cinterp_HN(ij,l,AI ,1) &
                           + (rhovyt(ijm1,TJ)+rhovyt(ij,TI)) * cinterp_HN(ij,l,AI ,2) &
                           + (rhovzt(ijm1,TJ)+rhovzt(ij,TI)) * cinterp_HN(ij,l,AI ,3) )
  
            flx_h(1,ij  ,k,l) =  flux * cinterp_PRA(ij  ,l) * dt
            flx_h(4,ip1j,k,l) = -flux * cinterp_PRA(ip1j,l) * dt
  