1      subroutine argos_cafe_epme(knode,kfrom,kto,bmod,grid,mgz,energy)
2c
3      implicit none
4c
5#include "argos_cafe_common.fh"
6#include "msgids.fh"
7c
8      integer mgz
9      integer knode(ngz),kfrom(np),kto(np)
10      real*8 bmod(ngmax,3),grid(2,ngx,ngy,mgz),energy
11c
12      real*8 fac,rh1,rh2,rh3,rsq,efac,vfac,ss
13      integer i,j,k,kg,m1,m2,m3
14      integer ngxy,nfx,nfy,nfz
15c
16      fac=(pi/ealpha)**2
17c
18      ngxy=ngx*ngy
19      nfx=ngx/2
20      if(2*nfx.lt.ngx) nfx=nfx+1
21      nfy=ngy/2
22      if(2*nfy.lt.ngy) nfy=nfy+1
23      nfz=ngz/2
24      if(2*nfz.lt.ngz) nfz=nfz+1
25c
26      epme=zero
27c
28      if(kfrom(me+1).gt.0) then
29      do 2 k=kfrom(me+1),kto(me+1)
30      kg=k+1-kfrom(me+1)
31      do 3 j=1,ngy
32      do 4 i=1,ngx
33      if(i+j+k.gt.3) then
34c
35      m1=i-1
36      if(i.gt.nfx) m1=m1-ngx
37      m2=j-1
38      if(j.gt.nfy) m2=m2-ngy
39      m3=k-1
40      if(k.gt.nfz) m3=m3-ngz
41c
42      rh1=recip(1,1)*m1+recip(1,2)*m2+recip(1,3)*m3
43      rh2=recip(2,1)*m1+recip(2,2)*m2+recip(2,3)*m3
44      rh3=recip(3,1)*m1+recip(3,2)*m2+recip(3,3)*m3
45      rsq=rh1*rh1+rh2*rh2+rh3*rh3
46c
47      efac=exp(-fac*rsq)/(pi*volume*bmod(i,1)*bmod(j,2)*bmod(k,3)*rsq)
48      vfac=two*(fac*rsq+one)/rsq
49      ss=grid(1,i,j,kg)**2+grid(2,i,j,kg)**2
50c
51      epme=epme+efac*ss
52      vpme(1)=vpme(1)+efac*ss*(vfac*rh1*rh1-one)
53      vpme(2)=vpme(2)+efac*ss*(vfac*rh1*rh2)
54      vpme(3)=vpme(3)+efac*ss*(vfac*rh1*rh3)
55      vpme(4)=vpme(4)+efac*ss*(vfac*rh2*rh2-one)
56      vpme(5)=vpme(5)+efac*ss*(vfac*rh2*rh3)
57      vpme(6)=vpme(6)+efac*ss*(vfac*rh3*rh3-one)
58c
59      grid(1,i,j,kg)=efac*grid(1,i,j,kg)
60      grid(2,i,j,kg)=efac*grid(2,i,j,kg)
61c
62      endif
63    4 continue
64    3 continue
65    2 continue
66      endif
67c
68      energy=0.5d0*epme
69c
70      return
71      end
72c $Id$
73