Subroutine normalize(phidim,nchidim,phi)
use Numeric_Kinds_Module
use Aph_Module
use SurfaceAph_Module
implicit none

INTEGER :: irho,itheta,ichi,aphindex,ir_pt
INTEGER :: phidim,nchidim
REAL(kind=dp)    :: rho_sum, theta_sum, chi_sum
COMPLEX(kind=dp) :: phi(phidim)

!==============================================================================
! Check normalization

  rho_sum=0.d0
  DO irho=1,nrho
   theta_sum=0.d0
    DO itheta=1,ntheta
      chi_sum=0.d0
      DO ichi=1,nchidim
        ir_pt=aphindex(ntheta,nchidim,irho,itheta,ichi,.false.)
        chi_sum=chi_sum+(Conjg(phi(ir_pt))*phi(ir_pt))*deltachi
      ENDDO
      theta_sum=theta_sum+chi_sum!*weights(itheta)
    ENDDO
    rho_sum=rho_sum+theta_sum*deltarho
  ENDDO
  PRINT*,'---------------------------------------------------------------------'
  PRINT*,'Wave packet normalization = ', rho_sum

end subroutine normalize