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:     ! BC for t^n
      call appBCs('new')
      ! solve
      do i=1,nEle
        flux=0d0
        source=0d0
        Pc(:)=Ele(i)%PC(:) ! position of the element center
        
        do j=1,Ele(i)%SurfNum ! find flux
          Ps(:)=Ele(i)%SurfPC(j,:) ! position of the surface center
          Dist=sqrt(dot_product(Pc-Ps,Pc-Ps))
          area=Ele(i)%SurfArea(j)
          if(Ele(i)%Neib(j)<0)then ! for boundary surface
            Dist=2d0*Dist
  
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))".
  