      subroutine virt_part(itimestep, ntotal,nvirt,hsml,mass,x,vx,
     &           rho,u,p,itype) 

c----------------------------------------------------------------------
c   Subroutine to determine the information of virtual particles
c   Here only the Monaghan type virtual particles for the 2D shear
c   cavity driven problem are generated.
c     itimestep : Current time step                                 [in]
c     ntotal : Number of particles                                  [in]
c     nvirt  : Number of virtual particles                         [out]
c     hsml   : Smoothing Length                                 [in|out]
c     mass   : Particle masses                                  [in|out]
c     x      : Coordinates of all particles                     [in|out]
c     vx     : Velocities of all particles                      [in|out]
c     rho    : Density                                          [in|out]
c     u      : internal energy                                  [in|out]
c     itype   : type of particles                               [in|out]

      implicit none
      include 'param.f'

      integer itimestep, ntotal, nvirt, itype(maxn)
      double precision hsml(maxn),mass(maxn),x(dim,maxn),vx(dim,maxn),
     &                 rho(maxn), u(maxn), p(maxn)
      integer i, j, d, im, mp, np
      double precision xl, yl, dx, dy, v_inf
      integer i2
      i2=2

      if (vp_input) then          
                        
c        open(1,file="../data/xv_vp.dat")
c        open(2,file="../data/state_vp.dat")
c        open(3,file="../data/other_vp.dat")    
        open(1,file="../data2/xv_vp.txt")
        open(2,file="../data2/state_vp.txt")
        open(3,file="../data2/other_vp.txt")        
        read(1,*) nvirt
        do j = 1, nvirt   
          i = ntotal + j      
          read(1,*)im, (x(d, i),d = 1, dim), (vx(d, i),d = 1, dim)                     
          read(2,*)im, mass(i), rho(i), p(i), u(i)        
          read(3,*)im, itype(i), hsml(i)                            
        enddo  
        close(1)
        close(2) 
        close(3) 
      
	else 
       
	nvirt = 0
        mp = 160         ! symericke    
        np = 5           ! symericke
	xl = 0.4           ! symericke
	yl = 0.0125        ! symericke
	dx = xl / mp 
	dy = yl / np       ! pridano
c	v_inf = 1.e-3      ! rychlost pohybujici se horni desky, ted se nehybe

c     Monaghan type virtual particle on the Upper side

        do i = 1, 4*mp+1 + 1 + 1
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = (i-1)*dx/2 - xl - dx/2
          x(i2, ntotal + nvirt) = yl + dx/2   ! xl zmeneno na yl
c          vx(1, ntotal + nvirt) = v_inf
        vx(1, ntotal + nvirt) = 0.
	  vx(i2, ntotal + nvirt) = 0.
        enddo

c     Monaghan type virtual particle on the Lower side

        do i = 1, 4*mp+1 + 1 + 1
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = (i-1)*dx/2 - xl - dx/2
          x(i2, ntotal + nvirt) = -yl - dx/2
          vx(1, ntotal + nvirt) = 0.
	  vx(i2, ntotal + nvirt) = 0.
        enddo

c     Monaghan type virtual particle on the Left side

        do i = 1, 4*np-1 + 2                 ! mp zmeneno na np
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = -xl - dx/2
          x(i2, ntotal + nvirt) = i*dy/2 - yl - dx/2  ! dx zmeneno na dy
          vx(1, ntotal + nvirt) = 0.
	  vx(i2, ntotal + nvirt) = 0.
        enddo

c     Monaghan type virtual particle on the Right side

        do i = 1, 4*np-1 + 2            ! mp zmeneno na np
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = xl + dx/2  ! xl zmeneno na yl 
          x(i2, ntotal + nvirt) = i*dy/2 - yl - dx/2  ! dx zmeneno na dy
          vx(1, ntotal + nvirt) = 0.
	  vx(i2, ntotal + nvirt) = 0.
        enddo

	do i = 1, nvirt
c	  if (x(1, ntotal + i).le.0) then
c	  rho (ntotal + i) = 0.125 ! zmeneno 1000 na 500, 2.2 pro rho = 1
        rho (ntotal + i) = 0.25
c	  endif
c	  rho (ntotal + i) = 0.25
	  mass(ntotal + i) = rho (ntotal + i) * dx * dy  ! druhe dx zmeneno na dy
	  p(ntotal + i) = 0.1795
	  u(ntotal + i) = 1.795
	  itype(ntotal + i) = -1
	  hsml(ntotal + i) = 0  ! bylo tam dx, ale stejne to nema na nic vliv  
        enddo
        
      endif   

c      if (mod(itimestep,save_step).eq.0) then
c        open(1,file="../data/xv_vp.dat")
c        open(2,file="../data/state_vp.dat")
c        open(3,file="../data/other_vp.dat") 

         if (shocktube) then 
          open(1,file="../data1/xv_vp.txt")
          open(2,file="../data1/state_vp.txt")
          open(3,file="../data1/other_vp.txt") 
        endif
        
        if (shocktube_2D) then
          open(1,file="../data2/xv_vp.txt")
          open(2,file="../data2/state_vp.txt")
          open(3,file="../data2/other_vp.txt")
        endif
                        
        write(1,*) nvirt
        do i = ntotal + 1, ntotal + nvirt         
          write(1,1001) i, (x(d, i), d=1,dim), (vx(d, i), d = 1, dim)              
          write(2,1002) i, mass(i), rho(i), p(i), u(i)
          write(3,1003) i, itype(i), hsml(i)                               
        enddo       
1001    format(1x, I6, 6(2x, e14.8))
1002    format(1x, I6, 7(2x, e14.8)) 
1003    format(1x, I6, 2x, I4, 2x, e14.8)
        close(1)
        close(2) 
        close(3) 
c      endif 

      if (mod(itimestep,print_step).eq.0) then
        if (int_stat) then
         print *,' >> Statistics: Virtual boundary particles:'
         print *,'          Number of virtual particles:',NVIRT
        endif     
      endif

      end