c     Matthew Marko, PhD PE

      program spam

c     ----------------------------------------------------------------------
      implicit none

      integer ii,jj,kk,jj0,uu,vv,ctbar(3),ctZ,ctZmod
      integer ct,ct2,xx,yy,N,Nsphere,ts,fooint,StudyCt,CenterPt(2)
      integer Nts,LoopCT,LoopCT2,stop,start2,stop2,ctCyl,StopWeight
      integer xxx,yyy,uv(2),break,breakct,abreak,split,StartStop
      integer BSx,cubecyl,ctx(2)
      complex i
      double precision pi,dx3(3),third,Xloc,Yloc,Zloc,StartDefl,NewDefl
      double precision dx00, dxE, dTfactor,dTdxfactor,dTdxfactor2,TTxx
      double precision masX,rho_solid,foo3,Xr0,Yr0,Vydt,dt0x,dx0x
      double precision deflect, deflectWeight, Rred, Ered, bWeight,mass
      double precision h, c02, K_bulk, Vy0,foo,r,dwdx(3),t1,t2
      double precision g(3), foo2, foo3vv, dxJ, Arange, dxEj
      double precision dt,dt0,Dij(3,3),OMGij(3,3),Sij3(3,3)
      double precision dSdt(6),RijAB(3,3),Wij,rho_ij,Radius,RadiusFactor
      double precision Weight0,Weight,Pweight(3),Vtest(3),Tfract,Efrac
      double precision dx,dxd,dE,deltX,deltX2,dxi,a,angleFrac,angle
      double precision lambda,mu,Ey,poisson,TT,Time,Bulk_water
      double precision Csound,rho0min,dV,gamma,dVdx0
      double precision X1(3),X2(3),V1(3),V2(3)
      integer, allocatable, dimension(:) :: Mat, FixXY, IsSolid
      integer, allocatable, dimension(:) :: WeightXY
      integer, allocatable, dimension(:,:) :: Contact,ContactDir
      integer, allocatable, dimension(:,:) :: ContactDirCT
      integer, allocatable, dimension(:,:) :: ContactX,ContactDirX
      integer, allocatable, dimension(:,:) :: ContactDirCTX
      integer, allocatable, dimension(:,:) :: Links
      integer, allocatable, dimension(:) :: BreakContact, RadiusCrap
      double precision, allocatable, dimension(:,:) :: X,X0,X00
      double precision, allocatable, dimension(:,:) :: V,V0,dVdx
      double precision, allocatable, dimension(:,:) :: dVdxFct,dVdyFct
      double precision, allocatable, dimension(:,:) :: dVdzFct
      double precision, allocatable, dimension(:,:) :: V00
      double precision, allocatable, dimension(:,:) :: TcomboOut
      double precision, allocatable, dimension(:,:,:) :: DSDtN,DSDtN0
      double precision, allocatable, dimension(:,:,:) :: Tij,Eij,Eij0
      double precision, allocatable, dimension(:,:,:) :: TTij,Tij0
      double precision, allocatable, dimension(:,:,:) :: TijCombo
      double precision, allocatable, dimension(:) :: rho,rho0,P,P0
      double precision, allocatable, dimension(:) :: printVar,Tvm0
      double precision, allocatable, dimension(:) :: TimeFct
      double precision, allocatable, dimension(:) :: massFct
      double precision, allocatable, dimension(:,:) :: D0,Omg0
      double precision, allocatable, dimension(:,:) :: VonMiss
      double precision, allocatable, dimension(:,:) :: Xout,Yout,Zout
      double precision, allocatable, dimension(:,:) :: VoutX,VoutY,VoutZ
      double precision, allocatable, dimension(:,:) :: dVdy,Tvm
      double precision, allocatable, dimension(:,:) :: Pout,Rout
      double precision, allocatable, dimension(:,:) :: DPDT,DPDT0
      double precision, allocatable, dimension(:,:) :: T22o
      double precision, allocatable, dimension(:,:) :: T11,T12,T13
      double precision, allocatable, dimension(:,:) :: T21,T22,T23
      double precision, allocatable, dimension(:,:) :: T31,T32,T33
      double precision, allocatable, dimension(:,:) :: E22,E11
      double precision, allocatable, dimension(:) :: dx0
      double precision, allocatable, dimension(:) :: zOffset


c     ----------------------------------------------------------------------
c     These are the universal constants
c     ----------------------------------------------------------------------
      pi = 3.141592653589793 ! set the value of pi
      i = (0, 1) ! set the imaginary number i=sqrt(-1)
      g = (/ 0.0, 0.0, 0.0 /)
      third=1.0/3.0

c     ----------------------------------------------------------------------
c     Declare array sizes
c     ----------------------------------------------------------------------

      call ModelBlock(N)

      open (unit=1,file="SimSpec.dat",STATUS='OLD',ACTION='READ')
        read (1,*),Ey
        read (1,*),poisson
        read (1,*),rho_solid
        read (1,*),Vtest(1)
        read (1,*),Vtest(2)
        read (1,*),Vtest(3)
        read (1,*),gamma
        read (1,*),Nts
        read (1,*),LoopCT
        read (1,*),StudyCt
        read (1,*),StartDefl
        read (1,*),NewDefl
        read (1,*),dTfactor
        read (1,*),dTdxfactor
        read (1,*),stop
        read (1,*),Bulk_water
      close (1)  

      Vy0 = Vtest(2)
      K_bulk = Ey/(3*(1-(2*poisson))) ! bulk modulus of iron
      mu=Ey/(2*(1+poisson))
      lambda=K_bulk-(2*mu/3)
      Ered=2/((1-(poisson**2))/Ey)
      
      Nts=Nts*StudyCt

      ALLOCATE (X(N,3))
      ALLOCATE (X0(N,3))
      ALLOCATE (X00(N,3))
      ALLOCATE (V(N,3))
      ALLOCATE (V0(N,3))
      ALLOCATE (V00(N,3))
      ALLOCATE (dx0(N))
      ALLOCATE (dVdx(N,3))
      ALLOCATE (dVdxFct(N,Nts))
      ALLOCATE (dVdyFct(N,Nts))
      ALLOCATE (dVdzFct(N,Nts))
      ALLOCATE (Mat(N))
      ALLOCATE (FixXY(N))
      ALLOCATE (IsSolid(N))
      ALLOCATE (rho(N))
      ALLOCATE (rho0(N))
      ALLOCATE (massFct(N))
      ALLOCATE (P(N))
      ALLOCATE (P0(N))
      ALLOCATE (Tij(N,3,3))
      ALLOCATE (TTij(N,3,3))
      ALLOCATE (Tij0(N,3,3))
      ALLOCATE (TijCombo(N,3,3))
      ALLOCATE (Eij(N,3,3))
      ALLOCATE (Eij0(N,3,3))
      ALLOCATE (DSDtN(N,3,3))
      ALLOCATE (DSDtN0(N,3,3))
      ALLOCATE (DPDT(N,3))
      ALLOCATE (DPDT0(N,3))
      ALLOCATE (D0(N,6))
      ALLOCATE (Omg0(N,6))
      ALLOCATE (Links(N,N+1))
      ALLOCATE (Xout(N,Nts))
      ALLOCATE (Yout(N,Nts))
      ALLOCATE (Zout(N,Nts))
      ALLOCATE (VoutX(N,Nts))
      ALLOCATE (VoutY(N,Nts))
      ALLOCATE (VoutZ(N,Nts))
      ALLOCATE (dVdy(N,Nts))
      ALLOCATE (Pout(N,Nts))
      ALLOCATE (Rout(N,Nts))
      ALLOCATE (TimeFct(Nts))
      ALLOCATE (Tvm0(N))
      ALLOCATE (Tvm(N,Nts))
      ALLOCATE (TcomboOut(N,Nts))
      ALLOCATE (T11(N,Nts))
      ALLOCATE (T12(N,Nts))
      ALLOCATE (T13(N,Nts))
      ALLOCATE (T21(N,Nts))
      ALLOCATE (T22(N,Nts))
      ALLOCATE (T23(N,Nts))
      ALLOCATE (T31(N,Nts))
      ALLOCATE (T32(N,Nts))
      ALLOCATE (T33(N,Nts))
      ALLOCATE (T22o(N,Nts))
      ALLOCATE (E22(N,Nts))
      ALLOCATE (E11(N,Nts))
      ALLOCATE (VonMiss(N,Nts))
      ALLOCATE (Contact(N,(N+1)))
      ALLOCATE (ContactDir(N,(N+1)))
      ALLOCATE (ContactDirCT(N,3))
      ALLOCATE (ContactX(N,(N+1)))
      ALLOCATE (ContactDirX(N,(N+1)))
      ALLOCATE (ContactDirCTX(N,3))
      ALLOCATE (BreakContact(N))

      ALLOCATE (printVar(N))

      ALLOCATE (RadiusCrap(N))

c     ----------------------------------------------------------------------
c     Declare size and velocity of particles
c     ----------------------------------------------------------------------

      open (unit=1,file="Fixed.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),FixXY(ii)
      enddo
      close (1)  

      open (unit=1,file="Contact.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),Contact(ii,:)
      enddo
      close (1)  

      open (unit=1,file="ContactDir.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),ContactDir(ii,:)
      enddo
      close (1)  

      open (unit=1,file="ContactDirCT.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),ContactDirCT(ii,:)
      enddo
      close (1)  

      open (unit=1,file="ContactX.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),ContactX(ii,:)
      enddo
      close (1)  

      open (unit=1,file="ContactDirX.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),ContactDirX(ii,:)
      enddo
      close (1)  

      open (unit=1,file="ContactDirCTX.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),ContactDirCTX(ii,:)
      enddo
      close (1)  

      open (unit=1,file="Mass.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),massFct(ii)
      enddo
      close (1)  

      open (unit=1,file="dX.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),dx0(ii)
      enddo
      close (1)  

      open (unit=1,file="X0.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),X(ii,:)
      enddo
      close (1)  

      open (unit=1,file="V0.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),V(ii,:)
      enddo
      close (1)  

      open (unit=1,file="dVdX.dat",STATUS='OLD',ACTION='READ')
      do ii=1,N
        read (1,*),dVdx(ii,:)
      enddo
      close (1)  

      open (unit=1,file="rho0.dat",STATUS='OLD',ACTION='READ')
      do jj=1,N
        read (1,*),rho0(jj)
      enddo
      close (1)

      open (unit=1,file="rho_init.dat",STATUS='OLD',ACTION='READ')
      do jj=1,N
        read (1,*),rho(jj)
      enddo
      close (1)

      open (unit=1,file="IsSolid.dat",STATUS='OLD',ACTION='READ')
      do jj=1,N
        read (1,*),IsSolid(jj)
      enddo
      close (1)


c --------------------------------------------------------------
c --------------------------------------------------------------

      call SpeedSound(Ey,poisson,rho_solid,Csound)
      dt0=dTfactor*0.5*(0.25*(2*(MINVAL(dx0)))/Csound)

      X00=X
      open (unit = 1, file = "Xoriginal.dat")
        do ii=1,N
          write (1,*),X00(ii,:)
        enddo
      close (1)
      Time=0

      do ii=1,N
        do jj=1,3
          do kk=1,3
            Tij(ii,jj,kk)=0
            TTij(ii,jj,kk)=0
            Eij(ii,jj,kk)=0
          enddo
          if (Vtest(jj)==0) then
            Tij0(:,jj,jj)=0
          else
            Tij0(:,jj,jj)=Vtest(jj)/(abs(Vtest(jj)))
          endif
        enddo
      enddo
      Tij0=Tij0*(1e-9)
      deflect=0

      call cpu_time ( t1 )

      do ts=1,Nts
        call cpu_time ( t2 )
        print *,ts,'/',Nts,', Elapsed CPU time = ', (t2 - t1),' seconds'

        xxx=((ts-1)*StudyCt/Nts)
        if (xxx==0) then
          LoopCT2=3*LoopCT
        else
          LoopCT2=LoopCT
        endif
        deflectWeight=MINVAL(dx0)*(StartDefl+(xxx*NewDefl/StudyCt))

        do yyy=1,LoopCT2

        dt=dTdxfactor*(MINVAL(dx0))/(MAXVAL(abs(V))+(1e-10))
        if (abs(deflect)<deflectWeight) then
          dt0x=abs(dTdxfactor*(MINVAL(dx0))/Vtest(2))
        else
          dt0x=dt0
        endif
        if (dt0x<dt) then
          dt=dt0x
        endif
        if (dt0<dt) then
          dt=dt0
        endif
        Time=Time+dt

        if (ts<=stop) then
          if (abs(deflect)<deflectWeight) then
            deflect=deflect+(Vtest(2)*dt)
          endif
        else
          deflect=deflect
        endif


        call LinkList(N,X,MAXVAL(dx0),Links)
        call ContactStress(N,Tij,Contact,ContactDir,Tij0)

        Contact=ContactX
        ContactDir=ContactDirX
        ContactDirCT=ContactDirCTX

        do jj=1,N
          if (IsSolid(jj)==1) then
            ct=ContactX(jj,1)
            do jj0=1,Links(jj,1)
              ii=Links(jj,jj0+1)

              if (ii/=jj) then
                call kernel(X(ii,:),X(jj,:),2*dx0(jj),Wij,dwdx,r)
                do kk=1,3
                  dx3(kk)=abs(X(ii,kk)-X(jj,kk))
                enddo
                fooint=MAXLOC(dx3,1)
                dxE=dx0(jj)*(1+Eij(jj,fooint,fooint))
                if (dxE<(dx0(jj)*1.05)) then
                  dxE=(dx0(jj)*1.05)
                endif

                if (r<dxE) then
                  kk=1
                  do uu=2,(ContactX(jj,1)+1)
                    vv=ContactX(jj,uu)
                    if (ii==vv) then
                      kk=0
                    endif
                  enddo
                  if (kk==1) then
                    ct=ct+1
                    Contact(jj,ct+1)=ii
                    ContactDir(jj,ct+1)=fooint
                    ContactDirCT(jj,fooint)=ContactDirCT(jj,fooint)+1
                    Tij0(jj,fooint,fooint)=-1
                  endif
                endif
              endif
            enddo
            Contact(jj,1)=ct
            ContactDir(jj,1)=ct
          endif
        enddo

c ____________________________________________


c ____________________________________________

c_______Calculate the new velocity. __________

        X0=X
        V0=V
        do xx=1,3
          do ii=1,N
            dxE=dx0(ii)*(1+Eij(ii,xx,xx))
            foo = g(xx)
            mass=massFct(ii)
            if (IsSolid(ii)==1) then

            if (Contact(ii,1)>0) then
            do uu=2,((Contact(ii,1))+1)
              jj=Contact(ii,uu)
              if (ContactDir(ii,uu)==xx) then
                if (uu>(1+ContactX(ii,1))) then
                  call Compression(xx,ii,jj,N,X,Eij,Ey,dx0,mass,foo2,a)
                else
                  call Tension(xx,ii,jj,N,X,Eij,Ey,dx0,mass,foo,foo2,TT)
                endif
                foo=foo+foo2
              else
                yy=ContactDir(ii,uu)
                call Shear(xx,yy,ii,jj,N,X,Eij,mu,dx0,mass,foo2,TT)
                foo=foo+foo2
              endif
              enddo
            endif
            endif

            if (IsSolid(ii)==1) then
             do jj=1,N
               if (IsSolid(jj)==0) then
                 X1=X(ii,:)
                 X2=X(jj,:)
                 call LJfctSolid(xx,X1,X2,dx0(ii),P(ii),a)
                 foo=foo+(a/mass)
               endif
             enddo
            else
              do jj=1,N
                if (IsSolid(jj)==1) then

                  X1=X(ii,:)
                  X2=X(jj,:)
                  V1=V(ii,:)
                  V2=V(jj,:)
                  call LJfct(xx,X1,X2,V1,V2,dx0(ii),dt,a)

                  foo=foo+a
                endif
              enddo
            endif


            ct=2
            do jj0=1,Links(ii,1)
              jj=Links(ii,jj0+1)
              call kernel(X(ii,:),X(jj,:),2*dx0(ii),Wij,dwdx,r)
              if (IsSolid(jj)==0) then
                foo2=(P(jj)-P(ii))*(dx0(ii)**2)/mass
                foo2=-foo2*((X(ii,xx)-X(jj,xx))/r)
                foo=foo2
              endif
            enddo


            dVdx0=dVdx(ii,xx)
            dVdx(ii,xx) = foo
            V(ii,xx)=V(ii,xx)+(dt*(dVdx0+dVdx(ii,xx))/2)
          enddo
        enddo

        X0=X
        V00=V
        do xx=1,3
          do ii=1,N
            if (IsSolid(ii)==1) then
              if ((Contact(ii,1))>0) then
                foo=0
                do uu=2,((Contact(ii,1))+1)
                  jj=Contact(ii,uu)
                  foo=foo+V00(jj,xx)
                enddo
                foo=foo/(Contact(ii,1))
                foo=(foo+(V00(ii,xx)))/2
                V(ii,xx)=foo
              endif
            endif
          enddo
        enddo


c_______Definition of FixXY___________________

c       FixXY(ii)=0 --> Free Particle
c       FixXY(ii)=1 --> Fixed in space
c       FixXY(ii)=2 --> Fixed in X/1-direction
c       FixXY(ii)=3 --> Fixed in Y/2-direction
c       FixXY(ii)=4 --> Fixed in Z/3-direction
c       FixXY(ii)=5 --> Free in X/1-direction only
c       FixXY(ii)=6 --> Free in Y/2-direction only
c       FixXY(ii)=7 --> Free in Z/3-direction only
c       FixXY(ii)=8 --> Follows V_test
c       FixXY(ii)=9 --> Follows -V_test
c       FixXY(ii)=10 --> Follows V_test in the X/1-direction
c       FixXY(ii)=11 --> Follows -V_test in the X/1-direction
c       FixXY(ii)=12 --> Follows V_test in the Y/2-direction
c       FixXY(ii)=13 --> Follows -V_test in the Y/2-direction
c       FixXY(ii)=14 --> Follows V_test in the Z/3-direction
c       FixXY(ii)=15 --> Follows -V_test in the Z/3-direction

c_______End Definition of FixXY_______________


        foo2=0
c_______Calculate the new location. __________
        do uu=1,3
          do ii=1,N
            if (IsSolid(ii)==1) then
            dV=V(ii,uu)+V0(ii,uu)
            if (dV>(0.1*(dx0(ii))/dt)) then
              dV=0
              V(ii,uu)=0
            elseif (abs(dV)<(1e-6)) then
              dV=0
              V(ii,uu)=0
            endif
            
            if (FixXY(ii)==1) then
              X(ii,uu)=X00(ii,uu)
            elseif (FixXY(ii)==2) then
              if (uu==1) then
                X(ii,uu)=X00(ii,uu)
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==3) then
              if (uu==2) then
                X(ii,uu)=X00(ii,uu)
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==4) then
              if (uu==3) then
                X(ii,uu)=X00(ii,uu)
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==5) then
              if (uu==1) then
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              else
                X(ii,uu)=X00(ii,uu)
              endif
            elseif (FixXY(ii)==6) then
              if (uu==2) then
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              else
                X(ii,uu)=X00(ii,uu)
              endif
            elseif (FixXY(ii)==7) then
              if (uu==3) then
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              else
                X(ii,uu)=X00(ii,uu)
              endif
            elseif (FixXY(ii)==8) then
              if (uu==2) then
                X(ii,uu)=X00(ii,uu)+deflect
              else
                X(ii,uu)=X00(ii,uu)
              endif
            elseif (FixXY(ii)==9) then
              if (uu==2) then
                X(ii,uu)=X00(ii,uu)-deflect
              else
                X(ii,uu)=X00(ii,uu)
              endif
            elseif (FixXY(ii)==10) then
              if (uu==1) then
                X(ii,uu)=X00(ii,uu)+deflect
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==11) then
              if (uu==1) then
                X(ii,uu)=X00(ii,uu)-deflect
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==12) then
              if (uu==2) then
                X(ii,uu)=X00(ii,uu)+deflect
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==13) then
              if (uu==2) then
                X(ii,uu)=X00(ii,uu)-deflect
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==14) then
              if (uu==3) then
                X(ii,uu)=X00(ii,uu)+deflect
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            elseif (FixXY(ii)==15) then
              if (uu==3) then
                X(ii,uu)=X00(ii,uu)-deflect
              else
                X(ii,uu)=X(ii,uu)+(dt*dV/2)
              endif
            else
              X(ii,uu)=X(ii,uu)+(dt*dV/2)
            endif
            else
              X(ii,uu)=X(ii,uu)+(dt*dV/2)
            endif
          enddo
        enddo


c_______Get new Stress. ______________________
        call GetDensity(X,V,N,dx0,dt,massFct,rho,IsSolid,rho)


c Calculate the Tensial Stress
        Tij=Tij-TTij
        do xx=1,3
          do ii=1,N
            mass=massFct(ii)
            if (IsSolid(ii)==1) then

            if (Contact(ii,1)>0) then
              foo=0
              ct2=0
              do uu=2,((Contact(ii,1))+1)
                TTxx=0
                do vv=1,3
                  TTxx=TTxx+Tij0(ii,vv,vv)
                enddo

                jj=Contact(ii,uu)
                if (ContactDir(ii,uu)==xx) then
                 if (uu>(1+ContactX(ii,1))) then
                   call Compression(xx,ii,jj,N,X,Eij,Ey,dx0,mass,a,TT)
                 else
                   call Tension(xx,ii,jj,N,X,Eij,Ey,dx0,mass,TTxx,a,TT)
                 endif
                 foo=foo+TT
                 ct2=ct2+1
                endif
              enddo

              if (FixXY(ii)>0) then
                foo=foo/2
              endif
              if (ct2==0) then
                TTij(ii,xx,xx)=0*TTij(ii,xx,xx)
              else
                TTij(ii,xx,xx)=foo/ct2
              endif
            endif

            else
              TTij(ii,xx,xx)=(Bulk_water/3)*((rho(ii)/rho0(ii))-1)
            endif
          enddo
        enddo



c Calculate the Shear Stress
        do vv=1,3
          call GetUV(vv,uv)
          do ii=1,N
            if (IsSolid(ii)==1) then
              foo=0
              do uu=2,((Contact(ii,1))+1)
                jj=Contact(ii,uu)
                if (ContactDir(ii,uu)==uv(2)) then
                  call Shear(uv(1),uv(2),ii,jj,N,X,Eij,mu,dx0,mass,a,TT)
                  foo=foo+TT
                endif
              enddo

              if (FixXY(ii)==0) then
                if (ContactDirCT(ii,uv(2))==0) then
                  foo=0
                else
                  foo=foo/ContactDirCT(ii,uv(2))
                endif
              else
                foo=foo/2
              endif

              TTij(ii,uv(2),uv(1))=foo
              TTij(ii,uv(1),uv(2))=foo
            endif
          enddo
        enddo

        Tij=Tij+TTij


c _______________Get new strain_______________

        do kk=1,N
          do ii=1,3
            Eij(kk,ii,ii)=(1/Ey)*(Tij(kk,ii,ii))
            call GetUV(ii,uv)
            foo=(poisson/Ey)*(Tij(kk,uv(1),uv(1)))
            foo=foo+(poisson/Ey)*(Tij(kk,uv(2),uv(2)))
            Eij(kk,ii,ii)=Eij(kk,ii,ii)-foo
          enddo
          do ii=1,3
            do jj=1,3
              if (jj/=ii) then
                Eij(kk,ii,jj)=Tij(kk,ii,jj)/(2*mu)
              endif
            enddo
          enddo
        enddo




c _______________Save Data____________________

        do ii=1,N
          foo=(Tij(ii,1,1)-Tij(ii,2,2))**2
          foo=foo+(Tij(ii,2,2)-Tij(ii,3,3))**2
          foo=foo+(Tij(ii,1,1)-Tij(ii,3,3))**2
          foo2=(Tij(ii,1,2)**2)+(Tij(ii,2,3)**2)+(Tij(ii,3,1)**2)
          foo=sqrt((foo+(6*foo2))/2)
          Tvm0(ii)=foo
          P(ii)=(Tij(ii,1,1)+Tij(ii,2,2)+Tij(ii,3,3))/3
        enddo

        TimeFct(ts)=Time
        Xout(:,ts)=X(:,1)
        Yout(:,ts)=X(:,2)
        Zout(:,ts)=X(:,3)
        VoutX(:,ts)=V(:,1)
        VoutY(:,ts)=V(:,2)
        VoutZ(:,ts)=V(:,3)
        Rout(:,ts)=rho
        Pout(:,ts)=P
        T11(:,ts)=Tij(:,1,1)
        T12(:,ts)=Tij(:,1,2)
        T13(:,ts)=Tij(:,1,3)
        T21(:,ts)=Tij(:,2,1)
        T22(:,ts)=Tij(:,2,2)
        T23(:,ts)=Tij(:,2,3)
        T31(:,ts)=Tij(:,3,1)
        T32(:,ts)=Tij(:,3,2)
        T33(:,ts)=Tij(:,3,3)
        Tvm(:,ts)=Tvm0
        E22(:,ts)=Eij(:,2,2)
        E11(:,ts)=Eij(:,1,1)
        dVdxFct(:,ts)=dVdx(:,1)
        dVdyFct(:,ts)=dVdx(:,2)
        dVdzFct(:,ts)=dVdx(:,3)

        enddo
      enddo

      call cpu_time ( t2 )
      print *,'Elapsed CPU time = ', (t2 - t1),' seconds'


c -------------------------------------------------------
c  Print data to dat-files
c -------------------------------------------------------

      open (unit = 1, file = "VarDat.dat")
        write (1,*),Ey
        write (1,*),mu
        write (1,*),StudyCt
      close (1)

      open (unit = 1, file = "X.dat")
      do jj=1,Nts
        printVar=Xout(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "Y.dat")
      do jj=1,Nts
        printVar=Yout(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "Z.dat")
      do jj=1,Nts
        printVar=Zout(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "P.dat")
      do jj=1,Nts
        printVar=Pout(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "rho.dat")
      do jj=1,Nts
        printVar=Rout(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T11.dat")
      do jj=1,Nts
        printVar=T11(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T12.dat")
      do jj=1,Nts
        printVar=T12(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T13.dat")
      do jj=1,Nts
        printVar=T13(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T21.dat")
      do jj=1,Nts
        printVar=T21(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T22.dat")
      do jj=1,Nts
        printVar=T22(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T23.dat")
      do jj=1,Nts
        printVar=T23(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T31.dat")
      do jj=1,Nts
        printVar=T31(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T32.dat")
      do jj=1,Nts
        printVar=T32(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "T33.dat")
      do jj=1,Nts
        printVar=T33(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "E22.dat")
      do jj=1,Nts
        printVar=E22(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "E11.dat")
      do jj=1,Nts
        printVar=E11(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "VonMises.dat")
      do jj=1,Nts
        printVar=Tvm(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "Vx.dat")
      do jj=1,Nts
        printVar=VoutX(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "Vy.dat")
      do jj=1,Nts
        printVar=VoutY(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "Vz.dat")
      do jj=1,Nts
        printVar=VoutZ(:,jj)
        write (1,*),printVar
      enddo
      close(1)

      open (unit = 1, file = "Time.dat")
      do jj=1,Nts
        printVar=TimeFct(jj)
        write (1,*),printVar
      enddo
      close(1)

c ----------------------------

      open (unit = 1, file = "Fixed.dat")
      do jj=1,N
        write (1,*),FixXY(jj)
      enddo
      close (1)  

      open (unit = 1, file = "Contact.dat")
      do jj=1,N
        write (1,*),Contact(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "ContactDir.dat")
      do jj=1,N
        write (1,*),ContactDir(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "ContactDirCT.dat")
      do jj=1,N
        write (1,*),ContactDirCT(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "ContactX.dat")
      do jj=1,N
        write (1,*),ContactX(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "ContactDirX.dat")
      do jj=1,N
        write (1,*),ContactDirX(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "ContactDirCTX.dat")
      do jj=1,N
        write (1,*),ContactDirCTX(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "Links.dat")
      do jj=1,N
        write (1,*),Links(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "Mass.dat")
      do jj=1,N
        write (1,*),massFct(jj)
      enddo
      close (1)  

      open (unit = 1, file = "dX.dat")
      do jj=1,N
        write (1,*),dx0(jj)
      enddo
      close (1)  


      open (unit = 1, file = "X0.dat")
      do jj=1,N
        write (1,*),X(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "V0.dat")
      do jj=1,N
        write (1,*),V(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "dVdX.dat")
      do jj=1,N
        write (1,*),dVdxFct(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "dVdY.dat")
      do jj=1,N
        write (1,*),dVdyFct(jj,:)
      enddo
      close (1)  

      open (unit = 1, file = "dVdZ.dat")
      do jj=1,N
        write (1,*),dVdzFct(jj,:)
      enddo
      close (1)  


      end program 



c __________________________________________________________________________

      subroutine GetUV(ii,uv)
c     ----------------------------------------------------------------------
      implicit none

      integer, intent(in) :: ii
      integer, intent(out) :: uv(2)
      integer jj,ct

      ct=0
      do jj=1,3
        if (jj/=ii) then
          ct=ct+1
          uv(ct)=jj
        endif
      enddo

      END subroutine GetUV


c __________________________________________________________________________

      subroutine FindL(N,h,mu,X,V,P,Tij,DSDtN)
c     ----------------------------------------------------------------------
      implicit none

      integer, intent(in) :: N
      double precision, intent(in) :: Tij(N,3,3)
      double precision, intent(in) :: X(N,3),V(N,3),h
      double precision, intent(in) :: mu
      double precision, intent(out) :: DSDtN(N,3,3),P(N)

      integer ii,jj, kk, uu, vv
      double precision L(3,3),D(3,3),Omg(3,3),r,W,Wtotal
      double precision dV,dx,D0(3,3),Omg0(3,3),DSDt(3,3)
      double precision foo,foo1,foo2,foo3,foo4,Sij(3,3)
      double precision dwdx(3)

      do jj=1,N
        Wtotal = 0
        do vv=1,3
          do uu=1,3
            L(uu,vv)=0
          enddo
        enddo
        do ii=1,N
          call kernel(X(jj,:),X(ii,:),h,W,dwdx,r)
          do vv = 1,3
            do uu = 1,3
              dV=(V(ii,uu)-V(jj,uu))
              dx=(X(ii,vv)-X(jj,vv))
              if (dx==0) then
                L(uu,vv)=L(uu,vv)
              else
                L(uu,vv)=L(uu,vv)+(dwdx(vv)*dV/dx)
                Wtotal=Wtotal+W
              endif
            enddo
          enddo
        enddo
        L=L*0/Wtotal

        do vv=1,3
          do uu=1,3
            D(uu,vv)=(L(uu,vv)+L(vv,uu))/2
            Omg(uu,vv)=(L(uu,vv)-L(vv,uu))/2
          enddo
        enddo


        do vv=1,3
          do uu=1,3
            Sij(uu,vv)=Tij(jj,uu,vv)
            if (uu==vv) then
              Sij(uu,vv)=Sij(uu,vv)-P(jj)
            endif
          enddo
        enddo

        do vv=1,3
          do uu=1,3
            foo1=0
            foo2=0
            do kk=1,3
              foo1=foo1+(Sij(uu,kk)*Omg(vv,kk))
              foo2=foo2+(Omg(uu,kk)*Sij(kk,vv))
            enddo
            foo3=2*mu*D(uu,vv)
            if (uu==vv) then
              foo4=(-2*mu/3)*D(uu,vv)
            else
              foo4=0
            endif
            DSDtN(jj,uu,vv)=foo1+foo2+foo3+foo4
          enddo
        enddo



      enddo

      END subroutine FindL

c __________________________________________________________________________

      subroutine kernel(x1,x2,h,w,dwdx,r)   

c----------------------------------------------------------------------
c   Subroutine to calculate the smoothing kernel wij and its 
c   derivatives dwdxij.

c     r    : Distance between particles i and j                     [in]
c     dx   : x-, y- and z-distance between i and j                  [in]  
c     h : Smoothing length                                       [in]
c     w    : Kernel for all interaction pairs                      [out]
c     dwdx : Derivative of kernel with respect to x, y and z       [out]

      implicit none
      
      double precision, intent(in) :: x1(3),x2(3),h
      double precision, intent(out) :: w,dwdx(3),r
      integer i, j, d, skf      
      double precision q, dw, factor,dx(3),pi,u

      do i=1,3
        dx(i)=abs(X1(i)-X2(i))
      enddo
      r=(dx(1)**2)+(dx(2)**2)+(dx(3)**2)
      r=sqrt(r)

      pi = 3.14159265358979
      u=10/(7*pi)
      q = r/h 
      w = 0.e0
      do d=1,3         
        dwdx(d) = 0.e0
      enddo   

      if (r==0) then
        w=0
        dwdx(1)=0
        dwdx(2)=0
        dwdx(3)=0
      else
        if (q<1) then
          w=1+(-1.5*(q**2))+(0.75*(q**3))
          do i=1,3
            dwdx(i)=((9*r/(4*h))-3)*(dx(i)/(h**2))
          enddo
        else
          if (q<2) then
            w=((2-q)**3)/4
            do i=1,3
              dwdx(i)=(-3/(4*h*r))*((2-q)**2)*dx(i)
            enddo
          else
            w=0
            dwdx(1)=0
            dwdx(2)=0
            dwdx(3)=0
          endif
        endif
      endif
      w=w*u
      do i=1,3
        dwdx(i)=u*dwdx(i)
      enddo
		
      END subroutine kernel

c __________________________________________________________________________

      subroutine GetDPDT(N,X,V,dx,K,massFct,rho,DPDT)
c     ----------------------------------------------------------------------
      implicit none

      integer, intent(in) :: N
      double precision, intent(in) :: X(N,3),V(N,3)
      double precision, intent(in) :: dx(N),massFct(N)
      double precision, intent(in) :: K,rho(N)
      double precision, intent(out) :: DPDT(N,3)

      integer ii,jj, kk
      double precision W0, dwdx(3),foo,dV,r
      double precision h,mass

      do jj=1,N
        h=2*dx(jj)
        mass=massFct(jj)
        do ii=1,N
          call kernel(X(ii,:),X(jj,:),h,W0,dwdx,r)
          foo=0
          do kk=1,3
            dV=V(ii,kk)-V(jj,kk)
            foo=-K*mass*dV*dwdx(kk)/rho(ii)
            DPDT(jj,kk)=DPDT(jj,kk)+foo
          enddo
        enddo
      enddo

      END subroutine GetDPDT

c __________________________________________________________________________

      subroutine GetDensity0(X,N,dx,mass,rho0,IsSolid,rho)
c     ----------------------------------------------------------------------
      implicit none

      integer, intent(in) :: N, IsSolid(N)
      double precision, intent(in) :: X(N,3),rho0(N)
      double precision, intent(in) :: dx(N),mass(N)
      double precision, intent(out) :: rho(N)

      integer ii,jj, kk
      double precision W0, r0, fooN, fooD, h, r, dwdx(3)
      double precision rhoIn(N)

      rhoIn=rho0
      do jj=1,N
        fooN=0
        h = 3*dx(jj)
        do ii=1,N
          if (IsSolid(ii)==IsSolid(jj)) then
            call kernel(X(ii,:),X(jj,:),h,W0,dwdx,r)
            fooN=fooN+(mass(ii)*W0)
          endif
        enddo
        rho(jj)=((3.14159*4/3)*fooN/(h**3))
      enddo

      END subroutine GetDensity0


c __________________________________________________________________________

      subroutine GetDensity(X,V,N,dx,dt,massFct,rho0,IsSolid,rho)
c     ----------------------------------------------------------------------
      implicit none

      integer, intent(in) :: N, IsSolid(N)
      double precision, intent(in) :: X(N,3),V(N,3),rho0(N)
      double precision, intent(in) :: dt,dx(N),massFct(N)
      double precision, intent(out) :: rho(N)

      integer ii,jj, kk
      double precision W,r,dwdx(3),foo1,foo2
      double precision h,mass,rhoIn(N)

      rhoIn=rho0
      do jj=1,N
        h=dx(jj)
        mass=massFct(jj)
        foo1=0
        do ii=1,N
          if (IsSolid(ii)==IsSolid(jj)) then
            call kernel(X(ii,:),X(jj,:),h,W,dwdx,r)
            do kk=1,3
              foo1=foo1-((mass/rhoIn(ii))*V(ii,kk)*dwdx(kk))
            enddo
          else
            foo1=0
          endif
        enddo
        foo2=dt*foo1*rhoIn(jj)
        rho(jj)=rhoIn(jj)+foo2
      enddo


      END subroutine GetDensity

c ______________________________________________________

      subroutine SpeedSound(YM,poisson,rho0min,Csound)

      implicit none

      double precision, intent(in) :: YM,poisson,rho0min
      double precision, intent(out) :: Csound
      double precision CsoundT,CsoundS,G

      CsoundT=(YM*(1-poisson))/rho0min
      CsoundT=CsoundT/((1+poisson)*(1-(2*poisson)))
      CsoundT=sqrt(CsoundT)

      G=(3*(1-poisson)/(1+poisson))-1
      G=G*YM/(4*(1-(2*poisson)))
      CsoundS=sqrt(G/rho0min)

      if (CsoundT>CsoundS) then
        Csound=CsoundT
      else
        Csound=CsoundS
      endif

      end subroutine SpeedSound

c __________________________________________________________________________

      subroutine LJfct(p,X1,X2,V1,V2,r0,dt,a)

      implicit none

      integer, intent(in) :: p
      double precision, intent(in) :: X1(3),X2(3)
      double precision, intent(in) :: V1(3),V2(3)
      double precision, intent(in) :: r0,dt
      double precision, intent(out) :: a

      integer ii
      double precision foo,r,r00,coeff,dir,foo2

      r=0
      do ii=1,3
        r=r+((X1(ii)-X2(ii))**2)
      enddo
      r=sqrt(r)
      if ((abs(X1(p)-X2(p)))==0) then
        dir=0
      else
        dir=(X1(p)-X2(p))/abs(X1(p)-X2(p))
      endif

      if (r>r0) then
        foo=0
      else
        coeff=(((r0/r)**12)-((r0/r)**4))
        foo=(V1(p)-V2(p))*coeff*dir/dt
      endif
      a=foo

      end subroutine LJfct

c __________________________________________________________________________

c __________________________________________________________________________

      subroutine LJfctSolid(p,X1,X2,r0,P2,f)

      implicit none

      integer, intent(in) :: p
      double precision, intent(in) :: X1(3),X2(3)
      double precision, intent(in) :: P2,r0
      double precision, intent(out) :: f

      integer ii
      double precision foo,r,del,trig

      r=0
      do ii=1,3
        r=r+((X1(ii)-X2(ii))**2)
      enddo
      del=(X1(p)-X2(p))
      r=sqrt(r)
      trig=del/r


      if (r>r0) then
        foo=0
      else
        foo=trig*P2*r0*r0
      endif
      f=foo

      end subroutine LJfctSolid

c __________________________________________________________________________



c __________________________________________________________________________

      subroutine Tension(xx,ii,jj,N,X,E,Ey,dx0,mass,TTxx,a,TT)
      
      implicit none

      integer, intent(in) :: xx,ii,jj,N
      double precision, intent(in) :: X(N,3),E(N,3,3)
      double precision, intent(in) :: Ey,dx0(N),mass,TTxx
      double precision, intent(out) :: a,TT
      

      integer ct,uv(2)
      double precision foo,dx,strain,Ai,A0,dir
      double precision dxE,strainE,a2,TT1,TT2,TTxxR

      A0=(dx0(ii))**2
      dx=X(jj,xx)-X(ii,xx)
      dxE=(1+((E(ii,xx,xx)+E(jj,xx,xx))/2))*dx0(ii)

      if (abs(dx)>0) then
        dir=dx/abs(dx)
      else
        dir=0
      endif

      strain=(abs(dx)-dx0(ii))/dx0(ii)
      strainE=(abs(dx)-dxE)/dxE

      a=0
      if (ii==jj) then
        a=0
        TT1=0
      else

        call GetUV(xx,uv)
        Ai=A0*(1+E(ii,uv(1),uv(1)))
        Ai=Ai*(1+E(ii,uv(2),uv(2)))
        a=dir*Ai*Ey*strainE/mass

        if (strain>0) then
          TT1=Ey*strain
        else
          TT1=0
        endif
      endif

      call Compression(xx,ii,jj,N,X,E,Ey,dx0,mass,a2,TT2)

      TTxxR=TTxx
      if (TTxxR>0) then
        TT=TT1
      elseif (TTxxR<0) then
        TT=TT2
      else
        if (strainE<0) then
          TT=TT2
        elseif (strainE>0) then
          TT=TT1
        else
          TT=(TT1+TT2)/2
        endif
      endif
      
      
      end subroutine Tension

c __________________________________________________________________________

      subroutine Compression(xx,ii,jj,N,X,E,Ey,dx0,mass,a,TT)

      implicit none

      integer, intent(in) :: xx,ii,jj,N
      double precision, intent(in) :: X(N,3),E(N,3,3)
      double precision, intent(in) :: Ey,dx0(N),mass
      double precision, intent(out) :: a,TT

      integer ct,uv(2)
      double precision foo,dx,strain,Ai,A0,dir
      double precision dxE,strainE

      A0=(dx0(ii))**2
      dx=X(jj,xx)-X(ii,xx)
      dxE=(1+((E(ii,xx,xx)+E(jj,xx,xx))/2))*dx0(ii)

      if (abs(dx)>0) then
        dir=dx/abs(dx)
      else
        dir=0
      endif

      strain=(abs(dx)-dx0(ii))/dx0(ii)
      strainE=(abs(dx)-dxE)/dxE

      a=0
      if (ii==jj) then
        a=0
        TT=0
      else

        call GetUV(xx,uv)
        Ai=A0*(1+E(ii,uv(1),uv(1)))
        Ai=Ai*(1+E(ii,uv(2),uv(2)))
        a=dir*Ai*Ey*strainE/mass

        if (strain<0) then
          TT=Ey*strain
        else
          TT=0
        endif
      endif
      
      if (TT>0) then
        TT=0
        a=0
      endif

      end subroutine Compression





c __________________________________________________________________________

      subroutine Shear(xx,yy,ii,jj,N,X,E,G,dx0,mass,a,TT)

      implicit none

      integer, intent(in) :: xx,yy,ii,jj,N
      double precision, intent(in) :: X(N,3),E(N,3,3)
      double precision, intent(in) :: G,dx0(N),mass
      double precision, intent(out) :: a,TT

      integer ct,uv(2)
      double precision delta,gamma,L
      double precision foo,Ai,A0,dir

      L=dx0(ii)
      A0=(dx0(ii))**2
      delta=(X(jj,xx)-X(ii,xx))
      if (abs(delta)>0) then
        dir=delta/abs(delta)
      else
        dir=0
      endif
      delta=abs(delta)
      gamma=delta/L


      if (ii==jj) then
        a=0
        TT=0
      else
        if (gamma>0) then
          TT=G*gamma
        else
          TT=0
        endif

        if (gamma>0) then
          call GetUV(yy,uv)
          Ai=A0*(1+E(ii,uv(1),uv(1)))
          Ai=Ai*(1+E(ii,uv(2),uv(2)))
          a=dir*Ai*G*gamma/mass
        else
          a=0
        endif

      endif


      end subroutine Shear


c __________________________________________________________________________

      subroutine ContactStress(N,TijIn,Contact,ContactDir,TijOut)

      implicit none

      integer, intent(in) :: N,Contact(N,N+1),ContactDir(N,N+1)
      double precision, intent(in) :: TijIn(N,3,3)
      double precision, intent(out) :: TijOut(N,3,3)

      integer ii,jj,xx,uu,cycle,ContactDirCTx
      double precision Tij0(N,3,3),Tij(N,3,3)
      double precision foo,Total(3),Avg(3)

      Tij0=TijIn

      do cycle=1,1000
      do ii=1,N
        do xx=1,3
          Avg(xx)=0
          ContactDirCTx=0
          do jj=2,(Contact(ii,1))
            if (ContactDir(ii,jj)==xx) then
              uu=Contact(ii,jj)
              Avg(xx)=Avg(xx)+Tij0(uu,xx,xx)
              ContactDirCTx=ContactDirCTx+1
            endif
          enddo
          if (ContactDirCTx==0) then
            Tij(ii,xx,xx)=Tij0(ii,xx,xx)
          else
            foo=Avg(xx)/ContactDirCTx
            Tij(ii,xx,xx)=(Tij0(ii,xx,xx)+foo)/2
          endif
        enddo
      enddo
      Tij0=Tij
      enddo
      TijOut=Tij


      end subroutine ContactStress


c __________________________________________________________________________

      subroutine LinkList(N,X,h,Links)

      implicit none

      integer, intent(in) :: N
      double precision, intent(in) :: X(N,3),h
      integer, intent(out) :: Links(N,(N+1))

      integer ii,jj,xx,uu,ct,fooInt,LinkCT(3),LLloc(N,3)
      double precision foo,DistCT(3),Wij,dwdx(3),r

      do ii=1,3
        DistCT(ii)=MAXVAL(X(:,ii))-MINVAL(X(:,ii))
        LinkCT(ii)=CEILING(DistCT(ii)/h)
      enddo
      Links(:,:)=0

      do ii=1,N
        do jj=1,3
          foo=(X(ii,jj)-MINVAL(X(:,jj)))/DistCT(jj)
          foo=ceiling(foo*LinkCT(jj))
          if (foo==0) then
            foo=foo+1
          endif
          LLloc(ii,jj)=foo
        enddo
      enddo

      do ii=1,N
        ct=1
        do jj=1,N
          xx=1
          do uu=1,3
            fooInt=abs(LLloc(ii,uu)-LLloc(jj,uu))
            if (fooInt>1) then
              xx=0
            endif
          enddo
          if (xx==1) then
            call kernel(X(ii,:),X(jj,:),h,Wij,dwdx,r)
            if (Wij>0) then
              ct=ct+1
              Links(ii,ct)=jj
            endif
          endif
        enddo
        Links(ii,1)=ct-1
      enddo


      end subroutine LinkList



c __________________________________________________________________________


c ___________________________________________________________________________

