      subroutine input(x, vx, mass, rho, p, u, itype, hsml, ntotal)
      
c----------------------------------------------------------------------
c     Subroutine for loading or generating initial particle information

c     x-- coordinates of particles                                 [out]
c     vx-- velocities of particles                                 [out]
c     mass-- mass of particles                                     [out]
c     rho-- dnesities of particles                                 [out]
c     p-- pressure  of particles                                   [out]
c     u-- internal energy of particles                             [out]
c     itype-- types of particles                                   [out]
c     hsml-- smoothing lengths of particles                        [out]
c     ntotal-- total particle number                               [out]

      implicit none     
      include 'param.f'

      integer itype(maxn), ntotal       
      double precision x(dim, maxn), vx(dim, maxn), mass(maxn), 
     &                 p(maxn), u(maxn), hsml(maxn), rho(maxn)
      integer i, d, im      

c     load initial particle information from external disk file

      if(config_input) then    
                        
        open(1,file="../data2/f_xv.txt")      !
        open(2,file="../data2/f_state.txt")   !
        open(3,file="../data2/f_other.txt")   !
       
      
        write(*,*)'  **************************************************'
        write(*,*)'      Loading initial particle configuration...   '       
	read (1,*) ntotal 
        write(*,*)'      Total number of particles   ', ntotal    	
        write(*,*)'  **************************************************'	
        do i = 1, ntotal         
          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          
          
      else 
       
        if (shocktube) then  
          open(1,file="../data1/ini_xv.txt")      !
          open(2,file="../data1/ini_state.txt")   !
          open(3,file="../data1/ini_other.txt")   !
        endif  
        
        if (shocktube_2D) then
          open(1,file="../data2/ini_xv.txt")      !
          open(2,file="../data2/ini_state.txt")   !
          open(3,file="../data2/ini_other.txt")   !
        endif  
        
      if (shocktube) call shock_tube(x, vx, mass, rho, p, u, 
     &                    itype, hsml, ntotal)                

      if (shocktube_2D) call shock_tube_2D(x, vx, mass, rho, p, u, 
     &                      itype, hsml, ntotal)
        write(1,*), ntotal,dim,0
        do i = 1, ntotal 
          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, I5, 6(2x, e14.8)) 
1002    format(1x, I5, 7(2x, e14.8)) 
1003    format(1x, I5, 2x, I2, 2x, e14.8) 
        write(*,*)'  **************************************************'
        write(*,*)'      Initial particle configuration generated   '       
        write(*,*)'      Total number of particles   ', ntotal    	
        write(*,*)'  **************************************************' 

      endif

      close(1)
      close(2) 
      close(3) 

      end              
       
       
      subroutine shock_tube(x, vx, mass, rho, p, u, 
     &                        itype, hsml, ntotal)

c----------------------------------------------------------------------     
c     This subroutine is used to generate initial data for the 
c     1 d noh shock tube problem
c     x-- coordinates of particles                                 [out]
c     vx-- velocities of particles                                 [out]
c     mass-- mass of particles                                     [out]
c     rho-- dnesities of particles                                 [out]
c     p-- pressure  of particles                                   [out]
c     u-- internal energy of particles                             [out]
c     itype-- types of particles                                   [out]
c          =1   ideal gas
c     hsml-- smoothing lengths of particles                        [out]
c     ntotal-- total particle number                               [out]

      implicit none     
      include 'param.f'
      
      integer itype(maxn), ntotal
      double precision x(dim, maxn), vx(dim, maxn), mass(maxn),
     &     rho(maxn), p(maxn), u(maxn), hsml(maxn)
      integer i, d
      double precision space_x, mm     

      mm=2; ! nasobek originalniho poctu uzlu
      ntotal=400*mm
      space_x=0.8/(80.*mm) ! zmena z puvodni delky 0.6 na delku 0.8   
      
      do i=1,ntotal
        mass(i)=1./(400.*mm) ! zmena puvodni hmoty 0.75 na 1
        hsml(i)=0.015/mm
        itype(i)=1
        do d = 1, dim
          x(d,i) = 0. 
          vx(d,i) = 0.
        enddo        
      enddo                
                
      do i=1,320*mm
        x(1,i)=-0.8+space_x/4.*(i-1)
      enddo
      
      do i=320*mm+1,ntotal
        x(1,i)=0.+space_x*(i-320*mm)
      enddo   
 
 
c      do i=1,ntotal
c        if (x(1,i).le.1.e-8) then
c          u(i)=357.1
c          rho(i)=1000.
c          p(i)=1.
c        endif 
c        if (x(1,i).gt.1.e-8)  then
c          u(i)=300.
c          rho(i)=800
c          p(i)=0.8
c        endif        
c      enddo                           
      do i=1,ntotal
        if (x(1,i).le.1.e-8) then
          u(i)=2.5
          rho(i)=1.
          p(i)=1.
        endif 
        if (x(1,i).gt.1.e-8)  then
          u(i)=1.795
          rho(i)=0.25
          p(i)=0.1795
        endif        
      enddo        
	               
      end
      
      subroutine shock_tube_2D(x, vx, mass, rho, p, u, 
     &                        itype, hsml, ntotal)

c----------------------------------------------------------------------     
c     This subroutine is used to generate initial data for the 
c     2 d shear driven cavity probem with Re = 1
c     x-- coordinates of particles                                 [out]
c     vx-- velocities of particles                                 [out]
c     mass-- mass of particles                                     [out]
c     rho-- dnesities of particles                                 [out]
c     p-- pressure  of particles                                   [out]
c     u-- internal energy of particles                             [out]
c     itype-- types of particles                                   [out]
c          =2   water
c     h-- smoothing lengths of particles                           [out]
c     ntotal-- total particle number                               [out]

      implicit none     
      include 'param.f'
      
      integer dimm
      integer itype(maxn), ntotal, ntotal_1, ntotal_2
      double precision x(dim, maxn), vx(dim, maxn), mass(maxn),
     &     rho(maxn), p(maxn), u(maxn), hsml(maxn)
      integer i, j, d, m, n, mp, np, k, mp_1, mp_2, np_1, np_2, r
      double precision xl, yl, dx, dy ,xl_1, xl_2, dx_1, dx_2, dy_1,dy_2
      integer i2
      i2=2

c     Giving mass and smoothing length as well as other data.

c        %%% Inicializace polohy jednotlivych castic %%%   


        xl = 0.4               ! delka strany x (je to symetricke)
      mp_1 = 160               ! pocet castic v ose x (vlevo)
      np = 10              ! pocet castic v ose y
      ntotal_1 = mp_1 * np  ! celkovy pocet castic (vlevo)
      dx = xl/mp_1            ! rozestup castic v ose x
      dy = dx                 ! rozestup castic v ose y
      yl  = np*dx           ! delka strany y
      

      do i = 1, mp_1
        do j = 1, np
	    k = j + (i-1)*np
	    x(1, k) = -xl + (i-1)*dx + dx/2.
	    x(i2, k) =-yl/2 + (j-1)*dy + dy/2.
        enddo
      enddo
      
      mp_2 = mp_1              ! pocet castic v ose x (vpravo)
      ntotal_2 = mp_2 * np   ! celkovy pocet castic (vpravo)
      
      do i = 1, mp_2           
	  do j = 1, np
	    k = ntotal_1 + j + (i-1)*np
	    x(1, k) = (i-1)*dx + dx/2.
	    x(i2, k) = -yl/2 + (j-1)*dy + dy/2.
        enddo
      enddo
      


      ntotal = ntotal_1 + ntotal_2   ! celkovy pocet castic
      mp = mp_1 + mp_2               ! celkovy pocet castic v ose x


c      %%% Definovani vlastnosti castic %%%

      do i = 1,ntotal
        hsml(i)= dx
        itype(i)=1
        vx(1,i)=0.
        vx(i2,i)=0.
        if (x(1,i).le.0) then
          u(i)=2.5
          rho(i)=1.
          p(i)=1.
        endif
        if (x(1,i).gt.0) then 
          u(i)=1.795
          rho(i)=0.25
          p(i)=0.1795
        endif
        mass(i)=rho(i)/0.80068946e5/2
      enddo
      end  