anchor.diff:      allocate(presF1(DIMS,grid%nC))
       allocate(condQ(grid%nC))
       allocate(flowRho(grid%nP))
  -    allocate(flowRho1(grid%nP))
       allocate(advRho(grid%nC))
       allocate(advRhou(DIMS,grid%nC))
       allocate(flowRhoH(grid%nP))
  
positive.code:   double precision,allocatable::rho0(:),rhou0(:,:),rhoE0(:),p0(:),u0(:,:),temp0(:),Y0(:,:)
    
    ! auxiliary variables for the PISO iteration
    double precision,allocatable::rho1(:),p1(:),laP1(:),presF1(:,:),temp1(:),rhou1(:,:),flowRho1(:)
    
    ! scales
    double precision,allocatable::rhoScale !< density scale [kg/m^3]
  
negative.code:       end forall
      end if
      call setBC()
      call findMassFlow(grid,rhou,p,presF,dt,flowRho)
      call findAdv(grid,flowRho,advRho)
    end subroutine
    
  