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