anchor.diff:  !----------------------------------------------------------------------------- best with 100 columns
   
  -!> global variables for foeuler
  -module gVarEuler
  +!> global variables for fons
  +module gVarNS
     use moduleGrid
     public
     
  
positive.code:   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
  end module
  
  !> foeuler main program
  program foeuler
    use moduleBasicDataStruct
    use moduleFileIO
    use moduleGrid
  
negative.code:   double precision,allocatable::gradRhoE(:,:) !< gradient of rhoE
    double precision pWork !< pressure work done on block surface
    
    ! 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)
        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
      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
  