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:     ! BC for t^n
      call appBCs('new')
      ! solve
      ! $omp parallel do &
      ! $omp default(shared) &
      ! $omp private(flux,source,Dist,area,conds,newconds,buff,i,j,k) &
      ! $omp reduction(max:error)
      do i=1,nEle
        flux=0d0
        source=0d0
        do j=1,Ele(i)%SurfNum ! find flux
          Dist=sqrt(dot_product(Ele(i)%PC(:)-Ele(i)%SurfPC(j,:),Ele(i)%PC(:)-Ele(i)%SurfPC(j,:)))
          area=Ele(i)%SurfArea(j)
          if(Ele(i)%Neib(j)<0)then ! for boundary surface
            Dist=2d0*Dist
  