Subroutine spherical_harmonics(ndim,j,m,Yjm)
USE Numeric_Kinds_Module
USE CommonInfo_Module
USE PiFactrs_module
USE Jacobi_Module
USE Complex_Module
Implicit None


!====================================================================
! Routine to calculate spherical harmonics. Requires use of the 
!    routine PLM.f90 to calculate normalized associated legendre
!    polynomials.
!
! Variables:
!   j = angular momentum
!   m = projection of j
!   phi = phi in regular spherical polar coodinates
!   plmn = array of normalized associated legendre polynomials   
!   th = theta in regular spherical polar coordinates      
!   Yjm = array of spherical harmonics
!

Integer :: ndim,j,m,shpt
Real(kind=dp) :: phi(ndim),plmn(0:j),phipart,costh
Complex(kind=dp) :: Yjm(ndim)
!Complex(kind=dp) :: imag

!Print*,'Begin SPHERICAL_HARMONICS'
!Write(Out_Unit,*) 'Begin SPHERICAL_HARMONICS'

! Open(Unit=11,file='spharmpts.txt')

! Error statement for possible m values.

 If(Abs(m).gt.j) Then
       Write(Out_Unit,*) 'Error: m value out of range'
       Write(Out_Unit,*) 'Range = -j ... +j'
       Write(Out_Unit,*) 'j = ',j
       Write(Out_Unit,*) 'm = ',m
       Print*,'ERROR: m is out of range -j ... +j'
       Stop 'spherical_harmonics'
 End If

! imag=Cmplx(0.d0,1.d0,dp)

 phi=0.d0

 Do shpt=1,ndim
         costh=Cos(thetacap(shpt))
         Call plm(j,m,costh,plmn)
         phipart=((1.d0/twopi)**0.5d0)*Exp(imag*m*phi(shpt))
         Yjm(shpt)=plmn(j)*phipart
  !      Write(11,*) thetacap(shpt),Yjm(shpt)
 End Do

! Close(Unit=11)

!Print*,'End SPHERICAL_HARMONICS'
!Write(Out_Unit,*) 'End SPHERICAL_HARMONICS'

Return
End Subroutine spherical_harmonics