anchor.diff:      call writeGMSH(13,p,grid,BIND_BLOCK,'p',0,t)
       ! advance in time
       do iStep=1,nStep
  +      call grid%updateDualBlock()
  +      call grid%updateBlockVol()
  +      call grid%updateFacetNorm()
  +      call grid%updateIntfArea()
  +      call grid%updateIntfNorm()
         ! Lagrangian step
  +      tempMass(:)=rho(:)*grid%BlockVol(:)
  +      forall(i=1:grid%nNode)
  +        tempMom(:,i)=rhou(:,i)*grid%NodeVol(i)
  +      end forall
  +      tempEnergy(:)=rhoE(:)*grid%BlockVol(:)
  +      do i=1,grid%nNode
  +        do j=1,size(grid%NodeNeibBlock(i)%dat)
  +          tempMom(:,i)=tempMom(:,i)&
  +          &            -dt*grid%NBAreaVect(i)%dat(:,j)*p(grid%NodeNeibBlock(i)%dat(j))
  +        end do
  +      end do
  +      uIntf=itplNode2Intf(u,grid)
  +      pIntf=itplBlock2Intf(p,grid)
  +      do i=1,grid%nIntf
  +        m=grid%IntfNeibBlock(1,i)
  +        n=grid%IntfNeibBlock(2,i)
  +        pWork=dt*pIntf(i)*grid%IntfArea(i)*dot_product(grid%IntfNorm(:,i),uIntf(:,i))
  +        tempEnergy(m)=tempEnergy(m)-pWork
  +        tempEnergy(n)=tempEnergy(n)+pWork
  +      end do
  +      ! apply BC
  +      do i=1,grid%nFacet
  +        do j=1,grid%Facet(i)%nNode
  +          k=grid%Facet(i)%iNode(j)
  +          tempMom(:,k)=tempMom(:,k)&
  +          &            -grid%FacetNorm(:,i)*dot_product(grid%FacetNorm(:,i),tempMom(:,k))
  +        end do
  +      end do
  +      ! recover intensive state
  +      rhoNode=1d0!TODO:itplBlock2Node(rho,grid)
  +      forall(i=1:grid%nNode)
  +        rhou(:,i)=tempMom(:,i)/grid%NodeVol(i)
  +        u(:,i)=rhou(:,i)/rhoNode(i)
  +      end forall
  +      uBlock=0d0!TODO:itplNode2Block(u,grid)
  +      forall(i=1:grid%nBlock)
  +        rhoE(i)=tempEnergy(i)/grid%BlockVol(i)
  +      end forall
  +      rho(:)=tempMass(:)/grid%BlockVol(:)
  +      call mvGrid(grid,dt*u)
         ! Euler rezoning
  +      call mvGrid(grid,-dt*u)
  +      ! recover state
  +      rhoNode=1d0!TODO:itplBlock2Node(rho,grid)
  +      forall(i=1:grid%nNode)
  +        rhou(:,i)=tempMom(:,i)/grid%NodeVol(i)
  +        u(:,i)=rhou(:,i)/rhoNode(i)
  +      end forall
  +      uBlock=0d0!TODO:itplNode2Block(u,grid)
  +      forall(i=1:grid%nBlock)
  +        rhoE(i)=tempEnergy(i)/grid%BlockVol(i)
  +        E(i)=rhoE(i)/rho(i)
  +        p(i)=(E(i)-dot_product(uBlock(:,i),uBlock(:,i))/2d0)*rho(i)*(gamm-1d0)
  +      end forall
  +      rho(:)=tempMass(:)/grid%BlockVol(:)
         t=t+dt
         ! write results
         call writeGMSH(13,rho,grid,BIND_BLOCK,'rho',iStep,t)
         call writeGMSH(13,u,grid,BIND_NODE,'u',iStep,t)
         call writeGMSH(13,p,grid,BIND_BLOCK,'p',iStep,t)
       end do
  +    deallocate(rhoNode)
  +    deallocate(uBlock)
  +    deallocate(uIntf)
  +    deallocate(pIntf)
  +    deallocate(tempMass)
  +    deallocate(tempMom)
  +    deallocate(tempEnergy)
       close(13)
     else
     end if
  
positive.code:     allocate(E(grid%nBlock))
      allocate(rhou(DIMS,grid%nNode))
      allocate(rhoE(grid%nBlock))
      ! simulation control
      dt=1d-6
      dFinal=1d-4
  
negative.code: 
  !> foeuler main program
  program foeuler
    use moduleBasicDataStruct
    use moduleFileIO
    use moduleGrid
    use moduleGridOperation
    use moduleInterpolation
    use moduleFVMGrad
    use moduleFVMConvect
    use moduleMPIComm
    use gVarEuler
    double precision,allocatable::rhoNode(:) !< density at node
    double precision,allocatable::uBlock(:,:) !< velocity at block
    double precision,allocatable::uIntf(:,:) !< velocity at interface
    double precision,allocatable::pIntf(:) !< pressure at interface
    double precision,allocatable::tempMass(:) !< temporary block mass
    double precision,allocatable::tempMom(:,:) !< temporary extensive node momentum
    double precision,allocatable::tempEnergy(:) !< temporary extensive block energy
    double precision pWork !< pressure work done on block surface
    
    call initMPI()
    if(pidMPI==0)then
      ! read grid
      open(12,file='bin/gridGMSH5.msh',status='old')
      call readGMSH(12,grid)
      call grid%updateIntf()
      close(12)
      ! allocate storage
      allocate(rho(grid%nBlock))
  