anchor.diff:    subroutine clear()
       deallocate(iBC,rho,rhou,rhoE,u,p,temp,c)
       deallocate(y,dRho,dRhou,dRhoE)
  -    deallocate(JacP,JacC)
  +    deallocate(JacP,JacC,JacS)
       deallocate(precPiv,precRhs)
       call grid%clear()
       call bc%clear()
  
positive.code:   !$omp parallel do default(shared)&
    !$omp& private(ier)
    do i=1,grid%nC
      precRhs(:,i)=res([0,1,2,3,4]*grid%nC+i)
      call DGETRS('N',5,1,JacC(:,1:5,i),5,precPiv(:,i),precRhs(:,i),5,ier)
      z([0,1,2,3,4]*grid%nC+i)=precRhs(:,i)
    end do
    !$end omp parallel do
    do l=1,n ! additional Jacobi loops
      !$omp parallel do default(shared)&
      !$omp& private(j,m,ier)
      do i=1,grid%nC
        precRhs(:,i)=res([0,1,2,3,4]*grid%nC+i)
        do j=1,size(grid%neib,1)
          m=grid%neib(j,i)
          if(1<=m.and.m<=grid%nC)then
            precRhs(:,i)=precRhs(:,i)-matmul(JacC(:,j*5+1:j*5+5,i),z([0,1,2,3,4]*grid%nC+m))
          end if
        end do
      end do
  
negative.code:   integer(C_LONG)::nEq !< number of equations
    integer(C_LONG)::iStat(100),iPar(1) !< integer solver outputs and parameters
    double precision::rStat(100),rPar(1) !< real solver outputs and parameters
    real,allocatable::JacS(:,:,:) !< single precision version of JacC
    integer,allocatable::precPiv(:,:) !< pivoting arrays for Lapack
    real,allocatable::precRhs(:,:) !< RHS arrays for Lapack
  
  contains
    
  