
      subroutine ModelBlock(Nout)

      implicit none

      integer, intent(out) :: Nout
      integer N,ctbar(3),Nsphere,ct,ii,jj,kk,uu
      integer fooint,split
      double precision RadiusFactor,Radius,masX,dx00
      double precision angleFrac,angle,Yr0,dx3(3),third
      double precision Zloc,Yloc,Xloc,r,pi,rho_solid
      integer, allocatable, dimension(:) :: FixXY,IsSolid,BreakContact
      integer, allocatable, dimension(:,:) :: Contact,ContactX
      integer, allocatable, dimension(:,:) :: ContactDir,ContactDirX
      integer, allocatable, dimension(:,:) :: ContactDirCT,ContactDirCTX
      double precision, allocatable, dimension(:,:) :: X,V,dVdx
      double precision, allocatable, dimension(:) :: massFct,dx0
      double precision, allocatable, dimension(:) :: rho0,rho

      open (unit = 1, file = "SimSpec.dat")
        write (1,*),207e9        ! Young's Modulus of Steel
        write (1,*),0.3          ! Poisson's Ratio of Steel
        write (1,*),7800         ! Density of Steel
        write (1,*),0            ! Forced Density in X direction (Vtest(1))
        write (1,*),1e3         ! Forced Density in Y direction (Vtest(2))
        write (1,*),0            ! Forced Density in Z direction (Vtest(3))
        write (1,*),0.3          ! Gamma (blend stress)
        write (1,*),250          ! Number of Recorded Time Steps (Nts)
        write (1,*),6            ! Time Steps between recorded Time Steps (LoopCT)
        write (1,*),1            ! Time Steps between recorded Time Steps (StudyCt)
        write (1,*),0.18         ! (StartDefl)
        write (1,*),0.05         ! (NewDefl)
        write (1,*),0.25         ! (dTfactor)
        write (1,*),0.01         ! (dTdxfactor)
        write (1,*),250          ! (stop)
        write (1,*),2.15e9       ! (Bulk_water)
      close (1)




      ctbar(1)=3     ! 3
      ctbar(2)=11    ! 11
      ctbar(3)=3     ! 3
      masX=0.1
      rho_solid=7800


      third=1.0/3.0
      dx00=(masX/rho_solid)**third
      N=(ctbar(1)*ctbar(2)*ctbar(3))


c _____________________________________________
c                 ALLOCATION
c _____________________________________________


      ALLOCATE (FixXY(N))
      ALLOCATE (IsSolid(N))
      ALLOCATE (BreakContact(N))
      ALLOCATE (X(N,3))
      ALLOCATE (V(N,3))
      ALLOCATE (dVdx(N,3))
      ALLOCATE (rho(N))
      ALLOCATE (rho0(N))
      ALLOCATE (massFct(N))
      ALLOCATE (dx0(N))
      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))

c _____________________________________________




c Set the parameters for the steel bar
      ct = 0
      Zloc=-dx00
      do uu=1,ctbar(3)
        Zloc=Zloc+dx00
        Yloc=-dx00
        do jj=1,ctbar(2)
          Yloc = Yloc + dx00
          Xloc = 0
          do ii=1,ctbar(1)
            ct = ct + 1
            if (jj==ctbar(2)) then
              FixXY(ct)=12
            elseif (jj==1) then
              FixXY(ct)=13
            else
              FixXY(ct)=0
            endif
            Xloc = Xloc + dx00
            X(ct,1) = Xloc
            X(ct,2) = Yloc
            X(ct,3) = Zloc
            V(ct,1) = 0
            V(ct,2) = 0
            V(ct,3) = 0
            dVdX(ct,1) = 0
            dVdX(ct,2) = 0
            dVdX(ct,3) = 0
            rho(ct) = rho_solid
            rho0(ct) = rho_solid
            dx0(ct) = dx00
            IsSolid(ct) = 1
            massFct(ct) = masX
            if (jj==(ctbar(2)/2)) then
              BreakContact(ct)=1
            elseif (jj==(1+(ctbar(2)/2))) then
              BreakContact(ct)=2
            else
              BreakContact(ct)=0
            endif
          enddo
        enddo
      enddo



c Set Contact links
      do jj=1,N
        if (IsSolid(jj)==1) then
          ct=0
          do ii=1,N
            if (IsSolid(ii)==1) then
              if (ii/=jj) then
                r=0
                do kk=1,3
                  dx3(kk)=abs(X(ii,kk)-X(jj,kk))
                  r=r+((X(ii,kk)-X(jj,kk))**2)
                enddo
                r=sqrt(r)
                if (r<(dx0(jj)*1.1)) then
                  split=1
                  if (BreakContact(ii)==1) then
                    if (BreakContact(jj)==2) then
                      split=0
                    endif
                  endif
                  if (BreakContact(ii)==2) then
                    if (BreakContact(jj)==1) then
                      split=0
                    endif
                  endif
                  split=1
                  if (split==1) then
                    ct=ct+1
                    Contact(jj,ct+1)=ii
                    fooint=MAXLOC(dx3,1)
                    ContactDir(jj,ct+1)=fooint
                    ContactDirCT(jj,fooint)=ContactDirCT(jj,fooint)+1
                  endif
                endif
              endif
            endif
          enddo
          Contact(jj,1)=ct
          ContactDir(jj,1)=ct
        endif
      enddo
      
      ContactX=Contact
      ContactDirX=ContactDir
      ContactDirCTX=ContactDirCT



      open (unit = 1, file = "Fixed.dat")
      do jj=1,N
        write (1,*),FixXY(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 = "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 = "ContactDirCTX.dat")
      do jj=1,N
        write (1,*),ContactDirCTX(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,*),dVdx(jj,:)
      enddo
      close (1)

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

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


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


      Nout=N

      end subroutine ModelBlock



c __________________________________________________________________________
