1*
2* $Id$
3*
4
5      subroutine kerker_G_init()
6      implicit none
7
8#include "bafdecls.fh"
9#include "errquit.fh"
10
11*     **** local variables ****
12      integer npack0,nfft3d,G(3)
13      integer i,ccode
14      real*8  gg,gg0,g0
15      integer tmp1(2)
16
17*     **** external functions ****
18      real*8   control_kerker_g0
19      external control_kerker_g0
20      integer  G_indx,c_G_indx,control_code
21      external G_indx,c_G_indx,control_code
22
23*     **** common block used for kerker_G.F ****
24c     real*8 tg(nfft3d)
25      logical isband,dokerker
26      integer tg_indx,tg_hndl
27      common / kerker_g_block / tg_indx,tg_hndl,isband,dokerker
28
29
30      g0 = control_kerker_g0()
31      dokerker =  (g0.gt.0.0d0)
32      if (dokerker) then
33
34      ccode = control_code()
35      isband =((ccode.eq.5).or.(ccode.eq.13).or.(ccode.eq.14))
36
37      if (isband) then
38         call C3dB_nfft3d(1,nfft3d)
39         call Cram_npack(0,npack0)
40         G(1)= c_G_indx(1)
41         G(2)= c_G_indx(2)
42         G(3)= c_G_indx(3)
43      else
44         call D3dB_nfft3d(1,nfft3d)
45         call Pack_npack(0,npack0)
46         G(1)= G_indx(1)
47         G(2)= G_indx(2)
48         G(3)= G_indx(3)
49      end if
50
51      if (.not.BA_alloc_get(mt_dbl,npack0,
52     >                     'tg',tg_hndl,tg_indx))
53     >   call errquit('kerker_G_init:out of heap memory',0,MA_ERR)
54
55      if (.not.BA_push_get(mt_dbl,nfft3d,'tmp1',tmp1(2),tmp1(1)))
56     >   call errquit('kerker_G_init:out of stack memory',0,MA_ERR)
57
58
59      gg0 = g0*g0
60      do i = 1,nfft3d
61         gg  = ( dbl_mb(G(1)+i-1)*dbl_mb(G(1)+i-1)
62     >         + dbl_mb(G(2)+i-1)*dbl_mb(G(2)+i-1)
63     >         + dbl_mb(G(3)+i-1)*dbl_mb(G(3)+i-1))
64         dbl_mb(tmp1(1)+i-1) = gg/(gg+gg0)
65      end do
66
67      if (isband) then
68         call Cram_r_pack(0,dbl_mb(tmp1(1)))
69         call Cram_r_Copy(0,dbl_mb(tmp1(1)),dbl_mb(tg_indx))
70      else
71         call Pack_t_pack(0,dbl_mb(tmp1(1)))
72         call Pack_t_Copy(0,dbl_mb(tmp1(1)),dbl_mb(tg_indx))
73      end if
74      if (.not.BA_pop_stack(tmp1(2)))
75     >   call errquit('kerker_G_init:popping stack memory',0,MA_ERR)
76
77      end if
78      return
79      end
80
81      subroutine kerker_G_end()
82      implicit none
83
84#include "bafdecls.fh"
85#include "errquit.fh"
86
87*     **** common block used for kerker_G.F ****
88      logical isband,dokerker
89      integer tg_indx,tg_hndl
90      common / kerker_G_block / tg_indx,tg_hndl,isband,dokerker
91
92      if (dokerker) then
93         if (.not.BA_free_heap(tg_hndl))
94     >   call errquit('error freeing heap',0, MA_ERR)
95      end if
96      return
97      end
98
99
100      subroutine kerker_G(V)
101      implicit none
102      real*8     V(*)
103
104#include "bafdecls.fh"
105#include "errquit.fh"
106
107*     **** common block used for kerker_G.F ****
108      logical isband,dokerker
109      integer tg_indx,tg_hndl
110      common / kerker_g_block / tg_indx,tg_hndl,isband,dokerker
111
112*     **** local variables ****
113      integer nfft3d,n2ft3d,tmp1(2)
114      integer nx,ny,nz
115      real*8  scal1
116
117      if (dokerker) then
118
119      if (isband) then
120         call C3dB_nx(1,nx)
121         call C3dB_ny(1,ny)
122         call C3dB_nz(1,nz)
123         call C3dB_nfft3d(1,nfft3d)
124         n2ft3d = nfft3d
125      else
126         call D3dB_nx(1,nx)
127         call D3dB_ny(1,ny)
128         call D3dB_nz(1,nz)
129         call D3dB_nfft3d(1,nfft3d)
130         n2ft3d = 2*nfft3d
131      end if
132      scal1 = 1.0d0/dble(nx*ny*nz)
133
134
135      if (.not.BA_push_get(mt_dcpl,nfft3d,'tmp1',tmp1(2),tmp1(1)))
136     >   call errquit('kerker_G: out of stack memory',0, MA_ERR)
137      if (isband) then
138         call C3dB_rc_SMul(1,scal1,V,dcpl_mb(tmp1(1)))
139         call C3dB_rc_fft3f(1,dcpl_mb(tmp1(1)))
140         call Cram_c_pack(0,dcpl_mb(tmp1(1)))
141         call Cram_rc_Mul2(0,dbl_mb(tg_indx),dcpl_mb(tmp1(1)))
142         call Cram_c_unpack(0,dcpl_mb(tmp1(1)))
143         call C3dB_cr_fft3b(1,dcpl_mb(tmp1(1)))
144         call C3dB_cr_real(1,dcpl_mb(tmp1(1)),V)
145      else
146         call D3dB_r_SMul(1,scal1,V,dcpl_mb(tmp1(1)))
147         call D3dB_r_Zero_Ends(1,dcpl_mb(tmp1(1)))
148         call D3dB_rc_fft3f(1,dcpl_mb(tmp1(1)))
149         call Pack_c_pack(0,dcpl_mb(tmp1(1)))
150         call Pack_tc_Mul2(0,dbl_mb(tg_indx),dcpl_mb(tmp1(1)))
151         call Pack_c_unpack(0,dcpl_mb(tmp1(1)))
152         call D3dB_cr_fft3b(1,dcpl_mb(tmp1(1)))
153         call D3dB_r_Zero_Ends(1,dcpl_mb(tmp1(1)))
154         call dcopy(n2ft3d,dcpl_mb(tmp1(1)),1,V,1)
155      end if
156      if (.not.BA_pop_stack(tmp1(2)))
157     >   call errquit('kerker_G: popping stack memory',0,MA_ERR)
158
159
160      endif
161      return
162      end
163
164