anchor.diff:      do i=1,grid%nP
         m=grid%iEP(1,i)
         n=grid%iEP(2,i)
  -      if(m<=grid%nC.and.n<=grid%nC)then ! internal pairs
  +      if(m<=size(s,2).and.m<=size(f,3).and.n<=size(s,2).and.n<=size(f,3))then
           do j=1,size(s,1)
             if(abs(s(j,m)-s(j,n))<=tiny(1d0))then
               up=m
  
positive.code:     double precision,intent(in)::f(:,:,:) !< fluxes
      double precision,allocatable,intent(inout)::adv(:,:) !< advection output
      integer::up,dn
      double precision::flow(size(s,1)),fUp,fDn,df,r
      double precision,allocatable::gradF(:,:,:)
      
      call grid%up()
  
negative.code:           end if
            fUp=dot_product(f(:,j,up),grid%normP(:,i))
            fDn=dot_product(f(:,j,dn),grid%normP(:,i))
            if(m<=grid%nC.and.n<=grid%nC)then ! internal pairs
              df=dot_product(matmul(grid%p(:,dn)-grid%p(:,up),gradF(:,j*DIMS-(DIMS-1):j*DIMS,up)),&
              &              grid%normP(:,i))
              r=merge(2d0*df/(fDn-fUp)-1d0,0d0,abs(fDn-fUp)>tiny(1d0))
              flow=-grid%aP(i)*(fUp+0.5d0*vanAlbada(r)*(fDn-fUp))
              adv(j,m)=adv(j,m)+flow
              adv(j,n)=adv(j,n)-flow
            else ! boundary pairs
              flow=-grid%aP(i)*fUp
              adv(j,m)=adv(j,m)+flow
            end if
          end do
        end if
      end do
      deallocate(gradF)
  