anchor.diff:      call associateVector(c_v,v)
       
       ! find x increment from f increment v and save back to v
  -    !$omp parallel do
  -    do i=1,grid%nC ! TODO: need to figure out why cannot use forall with omp workshare here
  -      v((DIMS+1)*(i-1)+1:(DIMS+1)*(i-1)+DIMS)=v((DIMS+1)*(i-1)+1:(DIMS+1)*(i-1)+DIMS)/rho(i)
  -      v((DIMS+1)*(i-1)+(DIMS+1))=v((DIMS+1)*(i-1)+(DIMS+1))
  -    end do
  -    !$omp end parallel do
  +    !$omp parallel workshare
  +    v(1:(DIMS+1)*grid%nC)=v(1:(DIMS+1)*grid%nC)*pbcDiagPrec(:)
  +    !$omp end parallel workshare
       pbcPsol=0
       foo=dat
     end function
  
positive.code:   double precision,allocatable::rho0(:),rhou0(:,:),rhoE0(:),p0(:),u0(:,:),temp0(:),massFrac0(:,:)
    
    ! state at the beginning of an outer iteration
    double precision,allocatable::rho1(:)
    double precision,allocatable::p1(:)
    
    double precision,allocatable::visc(:) !< viscosity [Pa*s]
    double precision,allocatable::cond(:) !< thermal conductivity [W/m/K]
  
negative.code:   
    !> momentum and pressure preconditioning matrix setup/factor
    function pbcPSet(c_x,c_xScale,c_f,c_fScale,dat)
      use iso_c_binding
      use modSUNDIALS
      use modAdvection
      use modPressure
      type(C_PTR),value::c_x,c_xScale,c_f,c_fScale,dat
      integer(C_INT)::pbcPSet
      double precision,pointer::x(:),xScale(:),f(:),fScale(:)
  