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