
      subroutine ModelBlock(Nout)

      implicit none

      integer, intent(out) :: Nout
      integer N,ctPV(3),Nsphere,ct,ii,jj,kk,uu
      integer fooint,split
      double precision RadiusFactor,Radius,masX,massL,dx00,Bulk_water
      double precision angleFrac,angle,Yr0,dx3(3),third,RR,pi,P
      double precision Zloc,Yloc,Xloc,rho_solid,rho_liquid,rho_liquidC
      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(:) :: dx0,massFct
      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.0          ! Forced Density in X direction (Vtest(1))
        write (1,*),0.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,*),1000         ! 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,*),0.18         ! (StartDefl)
        write (1,*),0.05         ! (NewDefl)
        write (1,*),0.05         ! (dTfactor)
        write (1,*),0.01         ! (dTdxfactor)
        write (1,*),15           ! (stop)
        write (1,*),2.15e9       ! (Bulk_water)
      close (1)

      ctPV(1)=25     ! 12
      ctPV(2)=20     ! 9
      ctPV(3)=1      ! 1

      masX=0.1
      P=101135*1
      rho_solid=7800
      rho_liquid=1000
      Bulk_water=2.15e9

      rho_liquidC=rho_liquid*((P/Bulk_water)+1)
      print *,rho_liquidC

      third=1.0/3.0
      dx00=((masX/rho_solid)**third)
      massL=rho_liquidC*(dx00**3)

      N=0
      Xloc=(-ctPV(1)*dx00)-(dx00/2)
      do ii=1,(2*ctPV(1))
      Xloc=Xloc+dx00
      Yloc=(-ctPV(1)*dx00)-(dx00/2)
      do jj=1,(2*ctPV(1))
        Yloc=Yloc+dx00
        RR=SQRT((Xloc**2)+(Yloc**2))
        if (RR<(ctPV(1)*dx00)) then
          N=N+1
        endif
      enddo
      enddo
      N=N*ctPV(3)
      print *,N


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 _____________________________________________


      print *,N

c Set the parameters for the steel bar

      ct=0
      Zloc=-dx00
      do uu=1,(ctPV(3))
        Zloc=Zloc+dx00
        Xloc=(-ctPV(1)*dx00)-(dx00/2)
        do ii=1,(2*ctPV(1))
          Xloc=Xloc+dx00
          Yloc=(-ctPV(1)*dx00)-(dx00/2)
          do jj=1,(2*ctPV(1))
            Yloc=Yloc+dx00
            RR=SQRT((Xloc**2)+(Yloc**2))
            if (RR<(ctPV(1)*dx00)) then
              ct = ct + 1
              if (RR<(ctPV(2)*dx00)) then
                IsSolid(ct) = 0
                massFct(ct) = massL
                rho0(ct)=rho_liquid
                rho(ct)=rho_liquidC
              else
                IsSolid(ct) = 1
                massFct(ct) = masX
                rho0(ct)=rho_solid
                rho(ct)=rho_solid
              endif
              FixXY(ct)=0
              X(ct,1) = Xloc
              X(ct,2) = Yloc
              X(ct,3) = Zloc
              
            endif
          enddo
        enddo
      enddo

      dx0 = dx0 + (dx00)


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
                RR=0
                do kk=1,3
                  dx3(kk)=abs(X(ii,kk)-X(jj,kk))
                  RR=RR+((X(ii,kk)-X(jj,kk))**2)
                enddo
                RR=sqrt(RR)
                if (RR<(dx0(jj)*1.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
          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 __________________________________________________________________________
