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