subroutine gr(switch)

if (switch.eq.0) then
    ngr=0
    delg =box/(2*nihs)
    do i= 0,nihs
        g(i)=0
  end do
else if (switch.eq.1)then
    ngr=ngr+1
    do i=1,npart-1
        do j=i+1,npart-1
xr=x(i)-x(j)
xr= xr-box*niht(xr/box)
r =sqrt(xr**2)
if(r.lt.box/2)then
ig=int(r/delg)
g(ig)=g(ig)+2
endif
enddo
enddo
else if (switch.eq.2)then
do i=1,nihs
r=delg*(i+0.5)
vb=((i+1)**3-i**3)delg**3
nid = (4/3)*pi*vb*rho
g(i)=g(i)/(ngr*npart*nid)
enddo
endif
return
end