anchor.diff:          end if
         end do
       end do
  +    ! $omp end parallel do
       
       ! check whether to stop
       ! Note: ".not.(abs(newe(:))<huge(0d0))" is an alternative to "isnan(new(e))".
  
positive.code:   do i=1,nEle
      call findEleGradScal(i,Temp(1:nEle),gradT(i,:))
      forall(j=1:Ele(i)%SurfNum)
        auxTcorr(i,j)=dot_product(gradT(i,:),Paux(i,j,:))
      end forall
    end do
    !$omp end parallel do
  
negative.code:   use globalvar
    use CNaddvar
    use moduleWrite,only:t
    double precision flux,source,error,buff,Dist,conds,newconds,area
    double precision lookupTab ! functions
    allocate(newTemp(-nFacet:nEle))
    allocate(newe(1:nEle))
  