Fortran-Fehler: Unerwarteter-ANWEISUNG, FUNCTION-Anweisung (1)

Ich bin ein fortran-code zu finden, der radial distribution function (RDF) mit harten Kugeln in einem Zell-Modell.
Es ist noch nicht fertig, und jetzt habe ich einen Fehler. Ich bin Implementierung des Histogramms. Das ist mein code.

      implicit double precision (a-h,o-z)
      parameter(npart=3000)
      dimension x(0:npart),y(0:npart),z(0:npart)

c n=Number of particules
c rcel=Radius of the cell
c rpart=Radius of the particules

      pi=3.1415927
      write(*,*)'n,rcel,rpart,dr?'
      read(*,*)n,rcel,rpart,dr
      write(*,*)'nstep,dp'
      read(*,*)nstep,dp
      rpart2=(2*rpart) 

      nfatmax=rcel/dr ! Number of bins
      vtotal=(4/3)*pi*rcel*rcel*rcel
      dentotal=n/vtotal
      write(*,*)'Density of particles, volume and bins = '
      write(*,*)dentotal,vtotal,nfatmax

      x(0)=0
      y(0)=0
      z(0)=0

      write(*,'(a,/)')'Generating start configurations'      
      counter1=0
      counter2=0
      counter3=0

      k=0
      do i=1,n

   21 xx=rcel*(ran()-0.5)*2
      yy=rcel*(ran()-0.5)*2
      zz=rcel*(ran()-0.5)*2

      rr=xx**2+yy**2+zz**2
      dist=sqrt(rr)

      if(dist.gt.(rcel-rpart2))then  !Avoid particles outside the cell
      counter1=counter1+1
      go to 21
      end if

      if(dist.lt.rpart2)then    ! Avoid overlap with central particle
      counter2=counter2+1
      go to 21
      end if   

      if(i.ge.1)then
        do j=1,i-1,1
           sep2=(x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2
           sep=sqrt(sep2)
           if(sep.lt.rpart2)then
           counter3=counter3+1
           go to 21
           end if
        end do
      end if

      k=k+1
      x(k)=xx
      y(k)=yy
      z(k)=zz
      end do

      write(*,*)'Starting config'
      write(*,'(3f8.3)')(x(i),y(i),z(i),i=1,n)  

      counterA=counter1+counter2+counter3
      write(*,*)'Rejection = '
      write(*,*)counterA

c     Monte Carlo loop
      counter4=0
      counter5=0
      counter6=0

      do i = 1,nfatmax
         h(i) = 0       !!!! Error here!!!!!!
      end do
      nobs = 0
      naccept = 0 

      do i=1,nstep
        do j=1,n
         nobs = nobs + 1
         xil=x(j)+dp*(ran()-0.5)
         yil=y(j)+dp*(ran()-0.5)
         zil=z(j)+dp*(ran()-0.5)

         r2=(xil**2)+(yil**2)+(zil**2)
         r=sqrt(r2)

         if(r.gt.(rcel-rpart2))then
           counter4=counter4+1
           go to 444    ! Avoid particles outside the cell
         end if

         if(r.lt.rpart2)then
           counter5=counter5+1
           go to 444    ! Avoid overlap with central particle
         end if 

         do ii=1,j-1
           dist2=(x(ii)-xil)**2+(y(ii)-yil)**2+(z(ii)-zil)**2
           dist=sqrt(dist2)
           if(dist.lt.rpart2)then
             counter6=counter6+1
             go to 444               ! Avoid overlap wit particles
           end if
         end do

c        Accepted configuration 
         x(j)=xil
         y(j)=yil
         z(j)=zil
         naccept = naccept + 1

c        Rejected configuration
 444     continue

         do jj=1,n
           dist2=(x(jj))**2+(y(jj))**2+(z(jj))**2
           dist=sqrt(dist2)
           k=(dist/dr)+1
           h(k) = h(k)+1    !!!!!!!! Error here!!!!!!!!!
         end do
        enddo
      end do

      write(*,*)'Final config'
      write(*,'(3f8.3)')(x(j),y(j),z(j),j=1,n)

      counterB=counter4+counter5+counter6
      write(*,*)'Rejection ='
      write(*,*)counterB
      stop
      end