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

c----------------------------------------------------------------------
c   Subroutine to determine the information of ghost particles
c   
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), xl, yl,  ! xl,yl-delky stran
     &                 dx, dy
      integer i, j, d, s, mp, np, i2
      
      i2 = 2
     
      
      xl = 0.4
c      xl = 4.e-3
      yl = 0.0125
c      yl = 0.25e-3
      j = 0   ! citac ghost particles
      s = 3   ! s*hsml - polomer support domainu
      do i = 1, ntotal
        if (x(1,i).lt.(-xl+s*hsml(i)).and.x(i2,i).ge.(-yl+s*hsml(i))  !1
     &   .and.x(i2,i).le.(yl-s*hsml(i))) then       
            j = j + 1
            x(1,ntotal+j)  =  -2*xl-x(1,i)
            x(i2,ntotal+j)  =   x(i2,i)
            vx(1,ntotal+j) = -vx(1,i)
            vx(i2,ntotal+j) =  vx(i2,i)
            mass(ntotal+j) =  mass(i)
            rho(ntotal+j)  =   rho(i)
            u(ntotal+j)    =     u(i)
            itype(ntotal+j)=      -1
            hsml(ntotal+j) =  hsml(i)
        endif   
         
        if (x(1,i).gt.(xl-s*hsml(i)).and.x(i2,i).ge.(-yl+s*hsml(i)) !2
     &   .and.x(i2,i).le.(yl-s*hsml(i))) then       
            j = j + 1
            x(1,ntotal+j)  =  2*xl-x(1,i) 
            x(i2,ntotal+j)  =   x(i2,i)
            vx(1,ntotal+j) = -vx(1,i)
            vx(i2,ntotal+j) =  vx(i2,i)
            mass(ntotal+j) =  mass(i)
            rho(ntotal+j)  =   rho(i)
            u(ntotal+j)    =     u(i)
            itype(ntotal+j)=      -1
            hsml(ntotal+j) =  hsml(i)
        endif
         
        
        if (x(i2,i).lt.(-yl+s*hsml(i)).and.x(1,i).ge.(-xl+s*hsml(i)) !3
     &   .and.x(1,i).le.(xl-s*hsml(i))) then       
            j = j + 1
            x(1,ntotal+j)  =   x(1,i) 
            x(i2,ntotal+j)  =  -2*yl-x(i2,i)
            vx(1,ntotal+j) =  vx(1,i)
            vx(i2,ntotal+j) = -vx(i2,i)
            mass(ntotal+j) =  mass(i)
            rho(ntotal+j)  =   rho(i)
            u(ntotal+j)    =     u(i)
            itype(ntotal+j)=      -1
            hsml(ntotal+j) =  hsml(i)
        endif
         
        if (x(i2,i).gt.(yl-s*hsml(i)).and.x(1,i).ge.(-xl+s*hsml(i)) !4
     &   .and.x(1,i).le.(xl-s*hsml(i))) then       
            j = j + 1
            x(1,ntotal+j)  =   x(1,i) 
            x(i2,ntotal+j)  =  2*yl-x(i2,i)
            vx(1,ntotal+j) =  vx(1,i)
            vx(i2,ntotal+j) =  -vx(i2,i)
            mass(ntotal+j) =  mass(i)
            rho(ntotal+j)  =   rho(i)
            u(ntotal+j)    =     u(i)
            itype(ntotal+j)=      -1
            hsml(ntotal+j) =  hsml(i)
        endif
         
      if (x(1,i).lt.(-xl+s*hsml(i)).and.x(i2,i).lt.(-yl+s*hsml(i))) then  !5 
            j = j + 3
            x(1,ntotal+j-2)  =   -2*xl-x(1,i)
            x(1,ntotal+j-1)  =   x(1,i)
            x(1,ntotal+j)    =   -2*xl-x(1,i)
            x(i2,ntotal+j-2)  =  x(i2,i)
            x(i2,ntotal+j-1)  =  -2*yl-x(i2,i)
            x(i2,ntotal+j)    =  -2*yl-x(i2,i)
            vx(1,ntotal+j-2) =  -vx(1,i)
            vx(1,ntotal+j-1) =  vx(1,i)
            vx(1,ntotal+j)   =  -vx(1,i)
            vx(i2,ntotal+j-2) = vx(i2,i)
            vx(i2,ntotal+j-1) = -vx(i2,i)
            vx(i2,ntotal+j)   = -vx(i2,i)
            mass(ntotal+j-2) =  mass(i)
            mass(ntotal+j-1) =  mass(i)
            mass(ntotal+j)   =  mass(i)
            rho(ntotal+j-2)  =   rho(i)
            rho(ntotal+j-1)  =   rho(i)
            rho(ntotal+j)    =   rho(i)
            u(ntotal+j-2)    =     u(i)
            u(ntotal+j-1)    =     u(i)
            u(ntotal+j)      =     u(i)
            itype(ntotal+j-2)=      -1
            itype(ntotal+j-1)=      -1
            itype(ntotal+j)  =      -1
            hsml(ntotal+j-2) =  hsml(i)
            hsml(ntotal+j-1) =  hsml(i)
            hsml(ntotal+j)   =  hsml(i)
        endif
        
       if (x(1,i).gt.(xl-s*hsml(i)).and.x(i2,i).lt.(-yl+s*hsml(i))) then  !6
            j = j + 3
            x(1,ntotal+j-2)  =   2*xl-x(1,i) 
            x(1,ntotal+j-1)  =   x(1,i)
            x(1,ntotal+j)    =   2*xl-x(1,i)
            x(i2,ntotal+j-2)  =  x(i2,i)
            x(i2,ntotal+j-1)  =  -2*yl-x(i2,i)
            x(i2,ntotal+j)    =  -2*yl-x(i2,i)
            vx(1,ntotal+j-2) =  -vx(1,i)
            vx(1,ntotal+j-1) =  vx(1,i)
            vx(1,ntotal+j)   =  -vx(1,i)
            vx(i2,ntotal+j-2) = vx(i2,i)
            vx(i2,ntotal+j-1) = -vx(i2,i)
            vx(i2,ntotal+j)   = -vx(i2,i)
            mass(ntotal+j-2) =  mass(i)
            mass(ntotal+j-1) =  mass(i)
            mass(ntotal+j)   =  mass(i)
            rho(ntotal+j-2)  =   rho(i)
            rho(ntotal+j-1)  =   rho(i)
            rho(ntotal+j)    =   rho(i)
            u(ntotal+j-2)    =     u(i)
            u(ntotal+j-1)    =     u(i)
            u(ntotal+j)      =     u(i)
            itype(ntotal+j-2)=      -1
            itype(ntotal+j-1)=      -1
            itype(ntotal+j)  =      -1
            hsml(ntotal+j-2) =  hsml(i)
            hsml(ntotal+j-1) =  hsml(i)
            hsml(ntotal+j)   =  hsml(i)
        endif
        
       if (x(1,i).lt.(-xl+s*hsml(i)).and.x(i2,i).gt.(yl-s*hsml(i))) then  !7
            j = j + 3
            x(1,ntotal+j-2)  =   -2*xl-x(1,i)
            x(1,ntotal+j-1)  =   x(1,i)
            x(1,ntotal+j)    =   -2*xl-x(1,i)
            x(i2,ntotal+j-2)  =   x(i2,i)
            x(i2,ntotal+j-1)  =  2*yl-x(i2,i)
            x(i2,ntotal+j)    =  2*yl-x(i2,i)
            vx(1,ntotal+j-2) =  -vx(1,i)
            vx(1,ntotal+j-1) =  vx(1,i)
            vx(1,ntotal+j)   =  -vx(1,i)
            vx(i2,ntotal+j-2) =  vx(i2,i)
            vx(i2,ntotal+j-1) = -vx(i2,i)
            vx(i2,ntotal+j)   = -vx(i2,i)
            mass(ntotal+j-2) =  mass(i)
            mass(ntotal+j-1) =  mass(i)
            mass(ntotal+j)   =  mass(i)
            rho(ntotal+j-2)  =   rho(i)
            rho(ntotal+j-1)  =   rho(i)
            rho(ntotal+j)    =   rho(i)
            u(ntotal+j-2)    =     u(i)
            u(ntotal+j-1)    =     u(i)
            u(ntotal+j)      =     u(i)
            itype(ntotal+j-2)=      -1
            itype(ntotal+j-1)=      -1
            itype(ntotal+j)  =      -1
            hsml(ntotal+j-2) =  hsml(i)
            hsml(ntotal+j-1) =  hsml(i)
            hsml(ntotal+j)   =  hsml(i)
        endif
        
        if (x(1,i).gt.(xl-s*hsml(i)).and.x(i2,i).gt.(yl-s*hsml(i))) then  !8
            j = j + 3
            x(1,ntotal+j-2)  =   2*xl-x(1,i)  
            x(1,ntotal+j-1)  =   x(1,i)
            x(1,ntotal+j)    =   2*xl-x(1,i) 
            x(i2,ntotal+j-2)  =   x(i2,i)
            x(i2,ntotal+j-1)  =  2*yl-x(i2,i)
            x(i2,ntotal+j)    =  2*yl-x(i2,i)
            vx(1,ntotal+j-2) =  -vx(1,i)
            vx(1,ntotal+j-1) =  vx(1,i)
            vx(1,ntotal+j)   =  -vx(1,i)
            vx(i2,ntotal+j-2) =  vx(i2,i)
            vx(i2,ntotal+j-1) = -vx(i2,i)
            vx(i2,ntotal+j)   = -vx(i2,i)
            mass(ntotal+j-2) =  mass(i)
            mass(ntotal+j-1) =  mass(i)
            mass(ntotal+j)   =  mass(i)
            rho(ntotal+j-2)  =   rho(i)
            rho(ntotal+j-1)  =   rho(i)
            rho(ntotal+j)    =   rho(i)
            u(ntotal+j-2)    =     u(i)
            u(ntotal+j-1)    =     u(i)
            u(ntotal+j)      =     u(i)
            itype(ntotal+j-2)=      -1
            itype(ntotal+j-1)=      -1
            itype(ntotal+j)  =      -1
            hsml(ntotal+j-2) =  hsml(i)
            hsml(ntotal+j-1) =  hsml(i)
            hsml(ntotal+j)   =  hsml(i)
        endif
         
      enddo 

      nvirt = j    ! pocet ghost particles
     
      if ( wall_particles ) then
      
        mp = 40   
        np = 20
	!xl = 2.0e-3   
	!yl = 1.0e-3        ! pridano
	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, mp+1
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = (i-1)*dx 
          x(i2, ntotal + nvirt) = yl     ! 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, mp+1
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = (i-1)*dx
          x(i2, ntotal + nvirt) = 0.  
          vx(1, ntotal + nvirt) = 0.
	  vx(i2, ntotal + nvirt) = 0.
        enddo

c     Monaghan type virtual particle on the Left side

        do i = 1, np-1                   ! mp zmeneno na np
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = 0. 
          x(i2, ntotal + nvirt) = i*dy  ! 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, np-1             ! mp zmeneno na np
   	  nvirt = nvirt + 1
	  x(1, ntotal + nvirt) = xl    ! xl zmeneno na yl 
          x(i2, ntotal + nvirt) = i*dy   ! dx zmeneno na dy
          vx(1, ntotal + nvirt) = 0.
	  vx(i2, ntotal + nvirt) = 0.
        enddo

	do i = 1, nvirt - j
	  rho (ntotal + j + i) = 0  ! zmeneno 1000 na 500, 2.2 pro rho = 1
	  mass(ntotal + j + i) = rho (ntotal + i) * dx * dy  ! druhe dx zmeneno na dy
	  p(ntotal + j + i) = 0.
	  u(ntotal + j + i) = 0  ! 357.1
	  itype(ntotal + j + i) = -3
	  hsml(ntotal + j + i) = 0  ! bylo tam dx, ale stejne to nema na nic vliv  
      enddo     
      
      endif

      if ( no_slip )  then
        do i = 1, j
          vx(1, ntotal + i) = -vx(1, ntotal + i)
          vx(i2, ntotal + i) = vx(i2, ntotal + i)
        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