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(:)=matmul(grid%p(:,dn)-grid%p(:,up),gradV(:,:,up))
  
positive.code:     double precision,intent(in),optional::gradV(:,:) !< gradient of v
      double precision::vf,dV,r
      integer::up,dn
      
      call grid%up()
      if(.not.allocated(flow))then
  
negative.code:     do i=1,grid%nP
        m=grid%iEP(1,i)
        n=grid%iEP(2,i)
        sameDirection=dot_product(u(:,m)+u(:,n),grid%normP(:,i))>=0d0 ! upwind by velocity
        up=merge(m,n,sameDirection)
        dn=merge(m,n,.not.sameDirection)
        rhoUp=rho(up)
        rhoDn=rho(dn)
        if(m<=grid%nC.and.n<=grid%nC)then ! internal pairs
  