
      subroutine ModelBlock(Nout)

      implicit none

      integer, intent(out) :: Nout
      integer N,ctbar(3),Nsphere,ct,ii,jj,kk,DiskRat
      integer fooint
      double precision RadiusFactor,Radius,masX,masXdisk,dx0x,dx0xDisk
      double precision angleFrac,angle,Yr0,dx3(3)
      double precision Zloc,Zloc2,Yloc,Xloc,r,mass,pi,rho_solid
      integer, allocatable, dimension(:) :: FixXY,IsSolid
      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,rho0


c     Start the Specification Parameters

      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.0          ! Forced Density in X direction (Vtest(1))
        write (1,*),100.0        ! Forced Density in Y direction (Vtest(2))
        write (1,*),0.0          ! Forced Density in Z direction (Vtest(3))
        write (1,*),0.3          ! Gamma (blend stress)
        write (1,*),100          ! Number of Recorded Time Steps (Nts)
        write (1,*),1            ! Time Steps between recorded Time Steps (LoopCT)
        write (1,*),1            ! Time Steps between recorded Time Steps (StudyCt)
        write (1,*),5*2.54e-6    ! (StartDefl)
        write (1,*),5*2.54e-6    ! (NewDefl)
        write (1,*),0.05         ! (dTfactor)
        write (1,*),0.005         ! (dTdxfactor)
        write (1,*),15           ! (stop)
        write (1,*),2.15e9       ! (Bulk_water)
      close (1)



      ctbar(1)=101
      ctbar(2)=10
      ctbar(3)=1
      Radius=20.0*(2.54e-4)
      DiskRat=3
      dx0x=0.25*(Radius/20.0)
      pi=ACOS(-1.0)
      rho_solid=7800

      if ((ctbar(1)*dx0x)>(2*Radius)) then
        angleFrac=pi
      else
        angleFrac=2*ASIN((ctbar(1)*dx0x)/(2*Radius))
        angleFrac=angleFrac/2
      endif

      dx0xDisk=dx0x/DiskRat
      masX=(dx0x**3)*rho_solid
      masXdisk=(dx0xDisk**3)*rho_solid

      Nsphere=(angleFrac*Radius/dx0xDisk)
      N=((Nsphere*DiskRat)+(ctbar(1)*ctbar(2)))*ctbar(3)

c _____________________________________________
c                 ALLOCATION
c _____________________________________________



      ALLOCATE (FixXY(N))
      ALLOCATE (IsSolid(N))
      ALLOCATE (X(N,3))
      ALLOCATE (V(N,3))
      ALLOCATE (dVdx(N,3))
      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


      Yr0=(ctbar(2)*dx0x)+Radius
      ct=0
      Zloc=-dx0x
      do kk=1,ctbar(3)
        Zloc=Zloc+dx0x
        Yloc=-dx0x
        do jj=1,ctbar(2)
          Yloc=Yloc+dx0x
          Xloc=-((dx0x*ctbar(1))/2)-(dx0x/2)
          do ii=1,ctbar(1)
            Xloc=Xloc+dx0x
            ct=ct+1
            if (jj==1) then
              FixXY(ct)=1
            else
              FixXY(ct)=0
            endif
            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
            rho0(ct) = rho_solid
            massFct(ct) = masX
            dx0(ct) = (massFct(ct)/rho_solid)**0.333333
            IsSolid(ct) = 1
          enddo
        enddo

        Zloc2=Zloc-(dx0xDisk*(DiskRat+1)/2)
        do ii=1,DiskRat
          Zloc2=Zloc2+dx0xDisk
          angle=(angleFrac/2)+(angleFrac/(2*Nsphere))
          do jj=1,Nsphere
            angle=angle-(angleFrac/Nsphere)
            Xloc=-(SIN(angle))*Radius
            Yloc=Yr0-((COS(angle))*Radius)
            ct=ct+1
            FixXY(ct)=8
            X(ct,1) = Xloc
            X(ct,2) = Yloc
            X(ct,3) = Zloc2
            V(ct,1) = 0
            V(ct,2) = 0
            V(ct,3) = 0
            dVdX(ct,1) = 0
            dVdX(ct,2) = 0
            dVdX(ct,3) = 0
            rho0(ct) = rho_solid
            massFct(ct) = masXdisk
            dx0(ct) = (massFct(ct)/rho_solid)**0.333333
            IsSolid(ct) = 1
          enddo
        enddo
      enddo
      print *,N,ct


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
                  if (FixXY(jj)==FixXY(ii)) 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,*),rho0(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 __________________________________________________________________________
