MODULE m_ylm
        IMPLICIT NONE
        PRIVATE
!---> check whether  or not normalizations are needed
        REAL, ALLOCATABLE, SAVE  :: ynorm(:)
        INTEGER,           SAVE  :: lmaxd = -1  ! initial value
        PUBLIC ylm4,ylmnorm_init
#ifdef CPP_GPU 
        PUBLIC ylm4_dev
#endif
      CONTAINS 

#ifdef CPP_GPU
      ATTRIBUTES(device) SUBROUTINE ylm4_dev(lmax,v,ylm)
!************************************************************
!     generate the spherical harmonics for the vector v
!     using a stable upward recursion in l.  (see notes
!     by m. weinert.)
!          m.weinert   january 1982
!     modified by R. Podloucky (added in ynorm); July 1989
!     cleaned up    mw 1995
!
!     modified to make use of f90 constructs. note that
!     the normalization is an internal subroutine and hence
!     can only be called from here. also, no need to dimension
!     arrays for ynorm, done dynamically.          mw 1999
!
!     GPU version added
!             U.Alekseeva      Oktober 2018
!************************************************************
!      USE m_juDFT
      INTEGER, VALUE, INTENT (IN) :: lmax
      REAL,    DEVICE, INTENT (IN) :: v(3)
      COMPLEX, DEVICE, INTENT (OUT):: ylm( (lmax+1)**2 )

      REAL, PARAMETER   :: small = 1.0e-12

      INTEGER  :: l,lm0,m
      REAL     :: fac,x,y,z,xy,r,rxy,cth,sth,cph,sph,cph2
      REAL,allocatable,device     :: c(:),s(:)
      REAL,allocatable,device     :: p(:,:)
      COMPLEX  :: ylms
      REAL,allocatable,device     :: ynorm_dev(:)
      REAL    a,cd,fpi

      allocate(c(0:lmax),s(0:lmax),ynorm_dev((lmax+1)**2))
      allocate(p(0:lmax,0:lmax))

!--->    calculate norm

      fpi = 4.0 * 3.1415926535897932  
      DO l=0,lmax
         lm0 = l*(l+1) + 1
         cd=1.0
         a = sqrt( (2*l+1)/fpi )
         ynorm_dev(lm0) = a
         DO m=1,l
            cd=cd/( (l+1-m)*(l+m) )
            ynorm_dev(lm0+m)  = a*sqrt(cd)
            ynorm_dev(lm0-m) = ( (-1.0)**m )*ynorm_dev(lm0+m)
         ENDDO
      ENDDO

!--->    calculate sin and cos of theta and phi
      x = v(1)
      y = v(2)
      z = v(3)
      xy  = x*x + y*y
      r   = sqrt(xy + z*z)
      rxy = sqrt(xy)

      IF (r.GT.small) THEN
         cth = z/r
         sth = rxy/r
      ELSE
         sth = 0.0
         cth = 1.0
      ENDIF
      IF(rxy.GT.small) THEN
         cph = x/rxy
         sph = y/rxy
      ELSE
         cph = 1.0
         sph = 0.0
      ENDIF

!---> generate associated legendre functions for m.ge.0
      fac = 1.0
!---> loop over m values
      DO m = 0, lmax - 1
         fac = -(m+m-1)*fac
         p(m,m) = fac
         p(m+1,m) = (m+m+1)*cth*fac
!--->    recurse upward in l
         DO l = m+2, lmax
            p(l,m)=((l+l-1)*cth*p(l-1,m)-(l+m-1)*p(l-2,m))/(l-m)
         ENDDO
         fac = fac*sth
      ENDDO
      p(lmax,lmax) = -(lmax+lmax-1)*fac

!--->    determine sin and cos of phi
      s(0) = 0.0
      s(1) = sph
      c(0) = 1.0
      c(1) = cph
      cph2 = cph+cph
      DO m = 2, lmax
         s(m) = cph2*s(m-1)-s(m-2)
         c(m) = cph2*c(m-1)-c(m-2)
      ENDDO

!--->    multiply in the normalization factors
      DO l=0,lmax
         ylm(l*(l+1)+1) = ynorm_dev(l*(l+1)+1)*cmplx(p(l,0),0.0)
      ENDDO
      DO m = 1, lmax
         DO l = m, lmax
            lm0 = l*(l+1)+1
            ylms = p(l,m)*cmplx(c(m),s(m))
            ylm(lm0+m) = ynorm_dev(lm0+m)*ylms
            ylm(lm0-m) = conjg(ylms)*ynorm_dev(lm0-m)
         ENDDO
      ENDDO

      END SUBROUTINE ylm4_dev
#endif


      SUBROUTINE ylm4(lmax,v,ylm)
!************************************************************
!     generate the spherical harmonics for the vector v
!     using a stable upward recursion in l.  (see notes
!     by m. weinert.)
!          m.weinert   january 1982
!     modified by R. Podloucky (added in ynorm); July 1989
!     cleaned up    mw 1995
!
!     modified to make use of f90 constructs. note that
!     the normalization is an internal subroutine and hence
!     can only be called from here. also, no need to dimension
!     arrays for ynorm, done dynamically.          mw 1999
!************************************************************
      USE m_juDFT
      INTEGER, INTENT (IN) :: lmax
      REAL,    INTENT (IN) :: v(3)
      COMPLEX, INTENT (OUT):: ylm( (lmax+1)**2 )

      REAL,    PARAMETER   :: small = 1.0e-12

      INTEGER l,lm0,m
      REAL    fac,x,y,z,xy,r,rxy,cth,sth,cph,sph,cph2
      REAL    p(0:lmax,0:lmax),c(0:max(1,lmax)),s(0:max(1,lmax))
      COMPLEX ylms


      IF (lmax.GT.lmaxd) THEN
         CALL juDFT_error("lmaxd too small in ylm4")
      ENDIF

!--->    calculate sin and cos of theta and phi
      x = v(1)
      y = v(2)
      z = v(3)
      xy  = x*x + y*y
      r   = sqrt(xy + z*z)
      rxy = sqrt(xy)

      IF (r.GT.small) THEN
         cth = z/r
         sth = rxy/r
      ELSE
         sth = 0.0
         cth = 1.0
      ENDIF
      IF(rxy.GT.small) THEN
         cph = x/rxy
         sph = y/rxy
      ELSE
         cph = 1.0
         sph = 0.0
      ENDIF

!---> generate associated legendre functions for m.ge.0
      fac = 1.0
!---> loop over m values
      DO m = 0, lmax - 1
         fac = -(m+m-1)*fac
         p(m,m) = fac
         p(m+1,m) = (m+m+1)*cth*fac
!--->    recurse upward in l
         DO l = m+2, lmax
            p(l,m)=((l+l-1)*cth*p(l-1,m)-(l+m-1)*p(l-2,m))/(l-m)
         ENDDO
         fac = fac*sth
      ENDDO
      p(lmax,lmax) = -(lmax+lmax-1)*fac

!--->    determine sin and cos of phi
      s(0) = 0.0
      s(1) = sph
      c(0) = 1.0
      c(1) = cph
      cph2 = cph+cph
      DO m = 2, lmax
         s(m) = cph2*s(m-1)-s(m-2)
         c(m) = cph2*c(m-1)-c(m-2)
      ENDDO

!--->    multiply in the normalization factors
      DO l=0,lmax
         ylm(l*(l+1)+1) = ynorm(l*(l+1)+1)*cmplx(p(l,0),0.0)
      ENDDO
      DO m = 1, lmax
         DO l = m, lmax
            lm0 = l*(l+1)+1
            ylms = p(l,m)*cmplx(c(m),s(m))
            ylm(lm0+m) = ynorm(lm0+m)*ylms
            ylm(lm0-m) = conjg(ylms)*ynorm(lm0-m)
         ENDDO
      ENDDO

      RETURN

      END SUBROUTINE ylm4

      SUBROUTINE ylmnorm_init(lmax)
!********************************************************************
!     normalization constants for ylm (internal subroutine has access
!     to lmax and ynorm from above)
!********************************************************************
!$    USE OMP_LIB
      use m_juDFT
      USE m_constants, ONLY : pimach
      INTEGER,INTENT(IN):: lmax
      INTEGER l,lm0,m
      REAL    a,cd,fpi
      IF ( allocated(ynorm) ) DEALLOCATE(ynorm)
      ALLOCATE ( ynorm( (lmax+1)**2 ) )   ! allocate array
      lmaxd = lmax

      fpi = 4.0*pimach()
!$    if (omp_in_parallel()) call juDFT_error(  &
!$          "ylmnorm should not called in parallel",calledby="ylm4.f")
      DO l=0,lmax
         lm0 = l*(l+1) + 1
         cd=1.0
         a = sqrt( (2*l+1)/fpi )
         ynorm(lm0) = a
         DO m=1,l
            cd=cd/( (l+1-m)*(l+m) )
            ynorm(lm0+m)  = a*sqrt(cd)
            ynorm(lm0-m) = ( (-1.0)**m )*ynorm(lm0+m)
         ENDDO
      ENDDO

      RETURN
      END SUBROUTINE ylmnorm_init

      END MODULE m_ylm