anchor.diff:    double precision t !< current time
     double precision dt !< time step size
     double precision tFinal !< final time
  -  integer iStep !< current time step number
  -  integer nStep !< total number of time steps
  +  integer iWrite !< index of result output snapshot
   end module
   
  -!> foeuler main program
  -program foeuler
  +!> fons main program
  +program fons
     use moduleBasicDataStruct
     use moduleFileIO
     use moduleGrid
  
positive.code:   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
      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
      ! 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()
  end program
  
negative.code: !----------------------------------------------------------------------------- best with 100 columns
  
  !> global variables for fons
  module gVarNS
    use moduleGrid
    public
    
  