anchor.diff:      end forall
     end subroutine
     
  +  !> set the boundary conditions
  +  subroutine setBC(x)
  +    double precision,intent(in)::x(*) !< solution vector
  +    double precision::pGst,TGst,ptGst,TtGst,uGst(DIMS),rhoGst
  +    
  +    do i=1,grid%nP
  +    m=grid%iEP(1,i)
  +    n=grid%iEP(2,i)
  +    if(n>grid%nC)then
  +      select case(bc%t(iBC(n)))
  +      case(BC_WALL)
  +        rho(n)=x(m)
  +        rhou(:,n)=x(grid%nC*[1,2,3]+m)&
  +        &         -2d0*dot_product(x(grid%nC*[1,2,3]+m),grid%normP(:,i))*grid%normP(:,i)
  +        rhoE(n)=x(4*grid%nC+m)
  +      case(BC_IN_STATIC)
  +        pGst=bc%p(1,iBC(n))
  +        TGst=bc%p(2,iBC(n))
  +        uGst(:)=bc%p(3:5,iBC(n))
  +        rhoGst=pGst/r/TGst
  +        if(norm2(uGst)/sqrt(gamm*r*TGst)<1d0)then
  +          uGst(:)=x(grid%nC*[1,2,3]+m)/rhoGst
  +        end if
  +        rho(n)=rhoGst
  +        rhou(:,n)=rhoGst*uGst(:)
  +        rhoE(n)=rhoGst*(1d0/(gamm-1d0)*r*TGst+0.5d0*dot_product(uGst,uGst))
  +      case(BC_IN_TOTAL)
  +        ptGst=bc%p(1,iBC(n))
  +        TtGst=bc%p(2,iBC(n))
  +        uGst=bc%p(3:5,iBC(n))
  +        ! TODO: find static values
  +        rhoGst=pGst/r/TGst
  +        if(norm2(uGst)/sqrt(gamm*r*TGst)<1d0)then
  +          uGst(:)=x(grid%nC*[1,2,3]+m)/rhoGst
  +        end if
  +        rho(n)=rhoGst
  +        rhou(:,n)=rhoGst*uGst(:)
  +        rhoE(n)=rhoGst*(1d0/(gamm-1d0)*r*TGst+0.5d0*dot_product(uGst,uGst))
  +      case(BC_OUT)
  +        pGst=bc%p(1,iBC(n))
  +        ! TODO apply pressure
  +        rho(n)=x(m)
  +        rhou(:,n)=x(grid%nC*[1,2,3]+m)&
  +        &         -2d0*dot_product(x(grid%nC*[1,2,3]+m),grid%normP(:,i))*grid%normP(:,i)
  +        rhoE(n)=x(4*grid%nC+m)
  +      case default
  +      end select
  +    end if
  +  end do
  +  end subroutine
  +  
     !> write the state
     subroutine writeState(fName)
       use modFileIO
  
positive.code:   integer,parameter::DIMS=3 !< three dimensions
    double precision,parameter::RTOL=1d-7 !< relative tolerance
    double precision,parameter::ATOL=1d-13 !< absolute tolerance
    
    type(polyFvGrid)::grid !< computational grid
    type(condTab)::bc !< boundary conditions
  
negative.code:   double precision::rPara(*)
    integer::ier
    
    call setBC(x)
    call syncState(x)
    call findAdv(grid,rho,rhou,rhoE,p,dRho,dRhou,dRhoE)
    forall(i=1:grid%nC)
      dxdt(i)=dRho(i)/grid%v(i)
      dxdt(grid%nC+i)=dRhou(1,i)/grid%v(i)
      dxdt(2*grid%nC+i)=dRhou(2,i)/grid%v(i)
      dxdt(3*grid%nC+i)=dRhou(3,i)/grid%v(i)
      dxdt(4*grid%nC+i)=dRhoE(i)/grid%v(i)
    end forall
    ier=0
  end subroutine
  