!-----------------------------------------------------------------------
      real(kind=8) function ylmr(l,ig)
!-----------------------------------------------------------------------
!     real spherical harmonics,  l is combined index for lm  (l=1,..,25)
!     order:  s, p_x, p_z, p_y, d_xy, d_xz, d_z^2, d_yz, d_x^2-y^2,..
!     the states with l>9 should be still checked out

      use gvec
      use cnst
!
      implicit none
      integer l, ig
      real(kind=8) x,y,z,r
!
      if (ig.gt.ng) call error(' ylmr ',' ig.gt.ng ',ig)
      x = gx(ig,1)
      y = gx(ig,2)
      z = gx(ig,3)
!
!     yml is undefined when  g=0 and l>0
!
      r = max(sqrt(x*x+y*y+z*z),1.0d-6)
      x = x/r
      y = y/r
      z = z/r
!
!     only l=1 is ok also when  g=0
!
      if (l.eq.1)  ylmr = sqrt(1.0/fpi)
      if (l.eq.2)  ylmr = sqrt(3.0/fpi)*x
      if (l.eq.3)  ylmr = sqrt(3.0/fpi)*z
      if (l.eq.4)  ylmr = sqrt(3.0/fpi)*y
      if (l.eq.5)  ylmr = sqrt(15.0/fpi)*x*y
      if (l.eq.6)  ylmr = sqrt(15.0/fpi)*x*z
      if (l.eq.7)  ylmr = sqrt(5.0/fpi/4.0)*(2.0*z*z-x*x-y*y)
      if (l.eq.8)  ylmr = sqrt(15.0/fpi)*y*z
      if (l.eq.9)  ylmr = sqrt(15.0/fpi/4.0)*(x*x-y*y)
      if (l.eq.10) ylmr = sqrt(7./fpi)*5./2.*x*(x*x-0.6)
      if (l.eq.11) ylmr = sqrt(7./fpi)*5./2.*y*(y*y-0.6)
      if (l.eq.12) ylmr = sqrt(7.*15./fpi)*x*y*z
      if (l.eq.13) ylmr = sqrt(7./fpi)*5./2.*z*(z*z-0.6)
      if (l.eq.14) ylmr = sqrt(7.*15./fpi)/2.*z*(x*x-y*y)
      if (l.eq.15) ylmr = sqrt(7.*15./fpi)/2.*y*(z*z-x*x)
      if (l.eq.16) ylmr = sqrt(7.*15./fpi)/2.*x*(y*y-z*z)
      if (l.eq.17) ylmr = sqrt(3.*7./fpi)*5./4.*(x**4.+y**4.+z**4.-0.6)
      if (l.eq.18) ylmr = sqrt(9.*35./fpi)/2.*y*z*(y*y-z*z)
      if (l.eq.19) ylmr = sqrt(9.*35./fpi)/2.*x*z*(z*z-x*x)
      if (l.eq.20) ylmr = sqrt(9.*5./fpi)/4.*(x**4.-y**4.-6.*z*z*(x*x-y*y))
      if (l.eq.21) ylmr = sqrt(9.*35./fpi)/2.*x*y*(x*x-y*y)
      if (l.eq.22) ylmr = sqrt(9.*5./fpi)*7./2.*x*y*(z*z-1./7.)
      if (l.eq.23) ylmr = sqrt(9.*5./fpi)*7./2.*x*z*(y*y-1./7.)
      if (l.eq.24) ylmr = sqrt(9.*5./fpi)*7./2.*y*z*(x*x-1./7.)
      if (l.eq.25) ylmr = sqrt(9.*5./fpi/3.)*7./2.*(z**4.-0.5*(x**4.+y**4.)-6./7.*(z*z-0.5*(x*x+y*y))) 
      if (l.ge.26) call error(' ylmr',' higher l not programmed  l=',l)
      return
      end