anchor.diff:    double precision,allocatable::gradRhoE(:,:) !< gradient of rhoE
     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))
  -    allocate(u(DIMS,grid%nNode))
  -    allocate(p(grid%nBlock))
  -    allocate(E(grid%nBlock))
  -    allocate(rhou(DIMS,grid%nNode))
  -    allocate(rhoE(grid%nBlock))
  -    allocate(rhoNode(grid%nNode))
  -    allocate(uBlock(DIMS,grid%nBlock))
  -    allocate(uIntf(DIMS,grid%nIntf))
  -    allocate(pIntf(grid%nIntf))
  -    allocate(tempMass(grid%nBlock))
  -    allocate(tempMom(DIMS,grid%nNode))
  -    allocate(tempEnergy(grid%nBlock))
  -    allocate(gradRho(DIMS,grid%nBlock))
  -    allocate(gradRhou(DIMS,DIMS,grid%nNode))
  -    allocate(gradRhoE(DIMS,grid%nBlock))
  -    ! simulation control
  -    dt=1d-5
  -    tFinal=0.0005d0
  -    nStep=ceiling(tFinal/dt)
  -    ! initial value of variables
  -    call grid%updateBlockPos()
  -    gamm=1.4d0
  -    u(:,:)=0d0
  -    rhou(:,:)=0d0
  +  ! 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))
  +  allocate(u(DIMS,grid%nNode))
  +  allocate(p(grid%nBlock))
  +  allocate(E(grid%nBlock))
  +  allocate(rhou(DIMS,grid%nNode))
  +  allocate(rhoE(grid%nBlock))
  +  allocate(rhoNode(grid%nNode))
  +  allocate(uBlock(DIMS,grid%nBlock))
  +  allocate(uIntf(DIMS,grid%nIntf))
  +  allocate(pIntf(grid%nIntf))
  +  allocate(tempMass(grid%nBlock))
  +  allocate(tempMom(DIMS,grid%nNode))
  +  allocate(tempEnergy(grid%nBlock))
  +  allocate(gradRho(DIMS,grid%nBlock))
  +  allocate(gradRhou(DIMS,DIMS,grid%nNode))
  +  allocate(gradRhoE(DIMS,grid%nBlock))
  +  ! simulation control
  +  dt=1d-5
  +  tFinal=0.0005d0 !TODO: read simulation times
  +  ! initial value of variables
  +  call grid%updateBlockPos()
  +  gamm=1.4d0
  +  u(:,:)=0d0
  +  rhou(:,:)=0d0
  +  forall(i=1:grid%nBlock)
  +    rho(i)=merge(1d0,0.125d0,grid%BlockPos(1,i)<0.5d0)
  +    p(i)=merge(1d5,1d4,grid%BlockPos(1,i)<0.5d0)
  +    E(i)=p(i)/rho(i)/(gamm-1d0)
  +    rhoE(i)=rho(i)*E(i)
  +  end forall
  +  t=0d0
  +  iWrite=0
  +  ! write initial states
  +  open(13,file='rstNS.msh',status='replace')
  +  call writeGMSH(13,grid)
  +  call writeGMSH(13,rho,grid,BIND_BLOCK,'rho',iWrite,t)
  +  call writeGMSH(13,u,grid,BIND_NODE,'u',iWrite,t)
  +  call writeGMSH(13,p,grid,BIND_BLOCK,'p',iWrite,t)
  +  ! advance in time
  +  do while(t<tFinal)
  +    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=itplBlock2Node(rho,grid)
  +    forall(i=1:grid%nNode)
  +      rhou(:,i)=tempMom(:,i)/grid%NodeVol(i)
  +      u(:,i)=rhou(:,i)/rhoNode(i)
  +    end forall
  +    uBlock=itplNode2Block(u,grid)
       forall(i=1:grid%nBlock)
  -      rho(i)=merge(1d0,0.125d0,grid%BlockPos(1,i)<0.5d0)
  -      p(i)=merge(1d5,1d4,grid%BlockPos(1,i)<0.5d0)
  -      E(i)=p(i)/rho(i)/(gamm-1d0)
  -      rhoE(i)=rho(i)*E(i)
  +      rhoE(i)=tempEnergy(i)/grid%BlockVol(i)
       end forall
  -    t=0d0
  -    ! write initial states
  -    open(13,file='rstEuler.msh',status='replace')
  -    call writeGMSH(13,grid)
  -    call writeGMSH(13,rho,grid,BIND_BLOCK,'rho',0,t)
  -    call writeGMSH(13,u,grid,BIND_NODE,'u',0,t)
  -    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=itplBlock2Node(rho,grid)
  -      forall(i=1:grid%nNode)
  -        rhou(:,i)=tempMom(:,i)/grid%NodeVol(i)
  -        u(:,i)=rhou(:,i)/rhoNode(i)
  -      end forall
  -      uBlock=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
  -      gradRho=findGrad(rho,grid,BIND_BLOCK)
  -      gradRhou=findGrad(rhou,grid,BIND_NODE)
  -      gradRhoE=findGrad(rhoE,grid,BIND_BLOCK)
  -      tempMass=tempMass+findDispConvect(rho,BIND_BLOCK,-dt*u,grid,gradRho,limiter=vanLeer)
  -      tempMom=tempMom+findDispConvect(rhou,BIND_NODE,-dt*u,grid,gradRhou,limiter=vanLeer)
  -      tempEnergy=tempEnergy+findDispConvect(rhoE,BIND_BLOCK,-dt*u,grid,gradRhoE,limiter=vanLeer)
  -      call mvGrid(grid,-dt*u)
  -      ! recover state
  -      call grid%updateDualBlock()
  -      call grid%updateBlockVol()
  -      rhoNode=itplBlock2Node(rho,grid)
  -      forall(i=1:grid%nNode)
  -        rhou(:,i)=tempMom(:,i)/grid%NodeVol(i)
  -        u(:,i)=rhou(:,i)/rhoNode(i)
  -      end forall
  -      uBlock=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)
  -      call showProg(t/tFinal)
  -    end do
  -    write(*,*),''
  -    deallocate(rhoNode)
  -    deallocate(uBlock)
  -    deallocate(uIntf)
  -    deallocate(pIntf)
  -    deallocate(tempMass)
  -    deallocate(tempMom)
  -    deallocate(tempEnergy)
  -    close(13)
  -  else
  -  end if
  -  call finalMPI()
  +    rho(:)=tempMass(:)/grid%BlockVol(:)
  +    call mvGrid(grid,dt*u)
  +    ! Euler rezoning
  +    gradRho=findGrad(rho,grid,BIND_BLOCK)
  +    gradRhou=findGrad(rhou,grid,BIND_NODE)
  +    gradRhoE=findGrad(rhoE,grid,BIND_BLOCK)
  +    tempMass=tempMass+findDispConvect(rho,BIND_BLOCK,-dt*u,grid,gradRho,limiter=vanLeer)
  +    tempMom=tempMom+findDispConvect(rhou,BIND_NODE,-dt*u,grid,gradRhou,limiter=vanLeer)
  +    tempEnergy=tempEnergy+findDispConvect(rhoE,BIND_BLOCK,-dt*u,grid,gradRhoE,limiter=vanLeer)
  +    call mvGrid(grid,-dt*u)
  +    ! recover state
  +    call grid%updateDualBlock()
  +    call grid%updateBlockVol()
  +    rhoNode=itplBlock2Node(rho,grid)
  +    forall(i=1:grid%nNode)
  +      rhou(:,i)=tempMom(:,i)/grid%NodeVol(i)
  +      u(:,i)=rhou(:,i)/rhoNode(i)
  +    end forall
  +    uBlock=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
  +    iWrite=iWrite+1
  +    call writeGMSH(13,rho,grid,BIND_BLOCK,'rho',iWrite,t)
  +    call writeGMSH(13,u,grid,BIND_NODE,'u',iWrite,t)
  +    call writeGMSH(13,p,grid,BIND_BLOCK,'p',iWrite,t)
  +    call showProg(t/tFinal)
  +  end do
  +  write(*,*),''
  +  deallocate(rhoNode)
  +  deallocate(uBlock)
  +  deallocate(uIntf)
  +  deallocate(pIntf)
  +  deallocate(tempMass)
  +  deallocate(tempMom)
  +  deallocate(tempEnergy)
  +  close(13)
   end program
  
positive.code: !----------------------------------------------------------------------------- best with 100 columns
  
  !> global variables for foeuler
  module gVarEuler
    use moduleGrid
    public
    
  
negative.code:   use moduleInterpolation
    use moduleFVMGrad
    use moduleFVMConvect
    use moduleCLIO
    use gVarNS
    double precision,allocatable::rhoNode(:) !< density at node
    double precision,allocatable::uBlock(:,:) !< velocity at block
    double precision,allocatable::uIntf(:,:) !< velocity at interface
  