anchor.diff:        if(abs(mFlow(i))<=tiny(1d0))then ! no mass flux
           cycle
         else ! upwinding by mass flux
  -        if(mFlow(i)>=0d0)then
  -          up=m
  -          dn=n
  -        else
  -          up=n
  -          dn=m
  -        end if
  +        sameDirection=mFlow(i)>=0d0
  +        up=merge(m,n,sameDirection)
  +        dn=merge(m,n,.not.sameDirection)
           vf=v(up)
           if(m<=grid%nC.and.n<=grid%nC.and.present(gradV))then ! internal pairs and gradV available
             dV=dot_product(grid%p(:,dn)-grid%p(:,up),gradV(:,up))
  
positive.code:     double precision,intent(in),optional::gradV(:,:,:) !< gradient of v
      double precision::vf(size(v,1)),dV(size(v,1)),r(size(v,1))
      integer::up,dn
      
      call grid%up()
      if(.not.allocated(flow))then
  
negative.code:     double precision,intent(in),optional::gradRho(:,:) !< gradient of rho
      integer::up,dn
      double precision::vMN(DIMS),vMP(DIMS),vPN(DIMS),flux(DIMS),eps,rhoUp,rhoDn,dRho,r,rhof
      logical::sameDirection
      
      call grid%up()
      if(.not.allocated(flow))then
  