anchor.diff:    use globalvar
     use CNaddvar
     use moduleWrite,only:t
  -  double precision flux,source,error,buff,Pc(3),Ps(3),Dist,conds,newconds,area
  +  double precision flux,source,error,buff,Dist,conds,newconds,area
     double precision lookupTab ! functions
     allocate(newTemp(-nFacet:nEle))
     allocate(newe(1:nEle))
  
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:         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))".
  