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