1!
2!---------------------------------------------------------------
3subroutine green(y,lam,e,dvy,chi,vpot,ze2)
4!---------------------------------------------------------------
5   !
6   use kinds, only: DP
7   use ld1inc, only: grid
8   use radial_grids, only: ndmx, series
9   implicit none
10   !
11   ! I/O variables
12   !
13   integer :: lam ! l angular momentum
14   real(DP) :: y(ndmx), chi(ndmx), dvy(ndmx), vpot(ndmx)
15   real(DP) :: e, ze2
16   !
17   ! local variables
18   !
19   integer :: i, l1, ncross, imatch
20   real(DP) :: f(ndmx), g(ndmx)
21   real(DP) :: work(ndmx), int_0_inf_dr
22   real(DP) :: a(0:3), b(0:3), c0, c1, c2, c3, c4, b0e
23   real(DP) :: rr1, rr2, r1, r2
24   real(DP) :: ddx12, sqlhf, xl1, x4l6, x6l12, x8l20
25   real(DP) :: gi, gim1
26   real(DP) :: fac
27! set up constants and initialize
28   ddx12=grid%dx*grid%dx/12.0
29   l1=lam+1
30   sqlhf=(DBLE(lam)+.5)**2
31   xl1=l1
32   x4l6=4*lam+6
33   x6l12=6*lam+12
34   x8l20=8*lam+20
35! series developement of the potential and r.h.s. term near the origin
36   do i=1,4
37      y(i)=vpot(i)-ze2/grid%r(i)
38   end do
39   call series(y,grid%r,grid%r2,b)
40   do i=1,4
41      y(i)= dvy(i)/grid%r(i)**l1
42   end do
43   call series(y,grid%r,grid%r2,a)
44!
45! determine the position of the last change of sign of the f-function
46! f < 0 (approximatively) means classically allowed   region
47! f > 0         "           "        "      forbidden   "
48! set up the f-, g- and y-functions
49!
50   f(1)=ddx12*(grid%r2(1)*(vpot(1)-e)+sqlhf)
51   do i=2,grid%mesh
52      f(i)=ddx12*(grid%r2(i)*(vpot(i)-e)+sqlhf)
53      if( f(i) .ne. sign(f(i),f(i-1)) ) imatch=i
54   end do
55
56   do i=1,grid%mesh
57      f(i)=1-f(i)
58      g(i)=ddx12*dvy(i)*grid%r(i)*grid%sqr(i)
59      y(i)=0.d0
60   end do
61! wave-function in the first two points by series developement
62   b0e=b(0)-e
63   c0 = 1.0d0
64   c1=0.5*ze2*c0/xl1
65   c2=(c1*ze2+b0e*c0+a(0))/x4l6
66   c3=(c2*ze2+c1*b0e+b(1)*c0+a(1))/x6l12
67   c4=(c3*ze2+c2*b0e+c1*b(1)+b(2)*c0+a(2))/x8l20
68   r1=grid%r(1)
69   r2=grid%r(2)
70   rr1=(c0+r1*(c1+r1*(c2+r1*(c3+r1*c4))))*r1**l1
71   rr2=(c0+r2*(c1+r2*(c2+r2*(c3+r2*c4))))*r2**l1
72
73   y(1)=rr1/grid%sqr(1)
74   y(2)=rr2/grid%sqr(2)
75!
76! new set up of the g- function
77   gim1 = g(1)
78   do i=2,grid%mesh-1
79      gi=g(i)
80      g(i)=gim1+10.0d0*gi+g(i+1)
81      gim1=gi
82   end do
83! outward integration: numerov's algorithm
84   call outward(y,f,g,grid%mesh,imatch,ncross)
85
86! inward integration: froese's algorithm
87   call inward(y,f,g,grid%mesh,imatch)
88
89!-orthogonalize to the 0^th order eigenvector
90   do i=1,grid%mesh
91      y(i) = y(i)*grid%sqr(i)
92      work(i)=y(i)*chi(i)
93   end do
94
95   fac = int_0_inf_dr(work,grid,grid%mesh,2*lam+2)
96   do i=1,grid%mesh
97      y(i) = y(i) - fac*chi(i)
98   end do
99!
100   return
101
102end subroutine green
103