modified: src/electrostatics.f90

	modified:   src/makefile

src/electrostatics.f90

@@ -215,51 +215,74 @@
       use var
       real(kind=pr) :: dis,dis3,disx,disy,disz, &
           Extot,Eytot,Eztot,partmass
-      integer :: n1,n2,ip,jp,lp
+      integer :: n1,n2,ip,jp,lp,m1,m2,n1end
 
 
-      if (elforceScheme.eq.2) then
+      select case (elforceScheme)
+      case (2)
         fx_el=0._pr
         fy_el=0._pr
         fz_el=0._pr
-      endif
 
-      do 1 n1=2,np
-      do 2 n2=1,(n1-1)
-
-      if (n1.eq.n2) cycle
-      if ((elforceScheme.eq.3).and. &
-         ((ip(n1).lt.ip(n2)-1.or.ip(n1).gt.ip(n2)+1).or. &
-          (jp(n1).lt.jp(n2)-1.or.jp(n1).gt.jp(n2)+1).or. &
-          (lp(n1).lt.lp(n2)-1.or.lp(n1).gt.lp(n2)+1))) cycle
-
-! force of n2 acting on n1
-        disx=xp(n1)-xp(n2)
-        disy=yp(n1)-yp(n2)
-        disz=zp(n1)-zp(n2)
-        dis=(disx**2+disy**2+disz**2)**0.5
-        dis3=dis**3
-  
-        partmass=4._pr/3._pr*pi*rhop*partn(n1)*radp(n1)**3
+        do 1 n1=2,np
+        do 2 n2=1,(n1-1)
+          if (n1.eq.n2) cycle
+          call forcesCoulombN1N2(n1,n2)
+2       enddo
+1       enddo
+
+      case (3)
+        do 3 i=1,ii
+        do 3 j=1,jj
+        do 3 l=1,ll
+
+        do 4 m1=2,npic(i,j,l)
+        do 5 m2=1,(m1-1)
+          n1=nic(i,j,l,m1)
+          n2=nic(i,j,l,m2)
+          call forcesCoulombN1N2(n1,n2)
+5       enddo
+4       enddo
+3       enddo
 
-        fx_el(n1)= fx_el(n1) &
-                 + q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
-        fy_el(n1)= fy_el(n1) &
-                 + q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
-        fz_el(n1)= fz_el(n1) &
-                 + q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
+      end select
 
-        partmass=4._pr/3._pr*pi*rhop*partn(n2)*radp(n2)**3
+      return
+      end
 
-        fx_el(n2)= fx_el(n2) &
-                 - q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
-        fy_el(n2)= fy_el(n2) &
-                 - q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
-        fz_el(n2)= fz_el(n2) &
-                 - q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
 
-2     enddo
-1     enddo
+!####################################################################
+!> @author Holger Grosshans
+!> @brief compute Coulomb force of n2 acting on n1
+      subroutine forcesCoulombN1N2(n1,n2)
+      use var
+      real(kind=pr) :: dis,dis3,disx,disy,disz,partmass
+      integer :: n1,n2,ip,jp,lp
+
+
+      disx=xp(n1)-xp(n2)
+      disy=yp(n1)-yp(n2)
+      disz=zp(n1)-zp(n2)
+      dis=(disx**2+disy**2+disz**2)**0.5
+      dis3=dis**3
+  
+      partmass=4._pr/3._pr*pi*rhop*partn(n1)*radp(n1)**3
+
+      fx_el(n1)= fx_el(n1) &
+               + q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
+      fy_el(n1)= fy_el(n1) &
+               + q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
+      fz_el(n1)= fz_el(n1) &
+               + q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
+
+      partmass=4._pr/3._pr*pi*rhop*partn(n2)*radp(n2)**3
+
+      fx_el(n2)= fx_el(n2) &
+               - q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
+      fy_el(n2)= fy_el(n2) &
+               - q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
+      fz_el(n2)= fz_el(n2) &
+               - q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
 
 
       return

src/makefile

@@ -15,7 +15,7 @@ INC =
 # gcc:
 CMP = mpifort
 #O3: optimization
-FLAGS   = -O3 -mcmodel=medium
+FLAGS   = -O3 -mcmodel=medium -pg
 #debugging:
 #FLAGS   = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
 FFLAGS  = -c $(FLAGS)