anchor.diff:      integer,parameter::MAX_N_NEIB=25
       integer,parameter::MAX_L_WORK=2000
       integer,parameter::L_IWORK=200
  -    integer::nNeib,iNeib(MAX_N_NEIB),lwork,iwork(L_IWORK),ier
  +    integer::nNeib,iNeib(MAX_N_NEIB),lwork,iwork(L_IWORK),rank,ier
       double precision::dx(MAX_N_NEIB,DIMS),dv(MAX_N_NEIB,size(v,1)),stat(DIMS),work(MAX_L_WORK),rcond
       double precision::tmp
       
  
positive.code:       allocate(gradv(DIMS,m,merge(grid%nC,grid%nE,findCellOnly)))
      end if
      ! find gradient
      do i=1,grid%nE
        if(i>grid%nC.and.findCellOnly)then
          exit
        end if
        ! list of neighbors
        if(i>grid%nC)then ! at non-cell
          if(findCellOnly)then
            exit
          end if
          write(*,*),'findGradPolyVect: support cell element only'
          !TODO:nNeib,lNeib at non-cell
  
negative.code:       call DGELSD(nNeib,DIMS,m,dx,MAX_N_NEIB,dv,MAX_N_NEIB,stat,rcond,rank,work,lwork,iwork,ier)
        gradv(:,:,i)=dv(1:DIMS,:)
      end do
      !$end omp parallel do
      ! limit the gradient
      do i=1,grid%nP
        m=grid%iEP(1,i)
  