anchor.diff:      double precision,allocatable,intent(inout)::dRho(:) !< net flux of rho
       double precision,allocatable,intent(inout)::dRhou(:,:) !< net flux of rhou
       double precision,allocatable,intent(inout)::dRhoE(:) !< net flux of rhoE
  -    double precision,allocatable::u(:,:),H(:)
  -    double precision,allocatable::f1c(:,:),f2c(:,:,:),f3c(:,:)
  +    double precision,allocatable,save::u(:,:),H(:)
  +    double precision,allocatable,save::f1c(:,:),f2c(:,:,:),f3c(:,:)
       double precision::rhoAvg,uAvg(DIMS),HAvg,cAvg,uNormAvg,rhoJump,pJump,uJump(DIMS),uNormJump
       double precision::dFEntropy(DIMS+2),dFAcoustic1(DIMS+2),dFAcoustic2(DIMS+2),flow(DIMS+2)
       double precision,parameter::GAMM=1.4d0 ! TODO other gamma and real gas
  
positive.code:     dRhoE(:)=0d0
      ! find auxiliary state and flux vectors in cell
      k=minval([size(rho),size(rhou,2),size(rhoE),size(p)])
      allocate(u(DIMS,k))
      allocate(H(k))
      allocate(f1c(DIMS,k))
      allocate(f2c(DIMS,DIMS,k))
      allocate(f3c(DIMS,k))
      !$omp workshare
      forall(i=1:k)
        u(:,i)=rhou(:,i)/rho(i)
  
negative.code:       end if
      end do
      !$end omp parallel do
    end subroutine
    
    !> find advection due to flux f depending on vector s on otGrid
  