diff --git a/src/ppm_comp_pp_cell.f b/src/ppm_comp_pp_cell.f
new file mode 100644
index 0000000000000000000000000000000000000000..6c4aa778cb7ac1d86b9ecf2be408c73a973345f3
--- /dev/null
+++ b/src/ppm_comp_pp_cell.f
@@ -0,0 +1,799 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_comp_pp_cell
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_cell_si(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_cell_di(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_cell_sci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_cell_dci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#endif
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_cell_su(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_cell_du(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_cell_scu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_cell_dcu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_cell_st(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_cell_dt(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_cell_sct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_cell_dct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nsublist,Nm,clist,cutoff2,dpd,info)
+#endif
+#endif
+ !!! This routine computes kernel interactions by direct
+ !!! particle-particle interactions using cell lists.
+ !!!
+ !!! [NOTE]
+ !!! Loops over lda have been explicitly unrolled for the cases of
+ !!! `lda.EQ.1` and `lda.EQ.2` to allow vectorization in these cases.
+ !!!
+ !!! [WARNING]
+ !!! clist has to be created by the user program before
+ !!! calling this routine. ppm_neighlist_clist should be
+ !!! used for this. The reason is that this allows the
+ !!! user to have multiple cell lists and call this
+ !!! routine for the appropriate one. Do not forget to
+ !!! call ppm_neighlist_clist_destroy and to DELLOCATE
+ !!! Nm after the cell lists are no longer needed.
+ !!!
+ !!! [NOTE]
+ !!! dpd needs to be allocated to proper size before
+ !!! calling this routine. Also, this routine is not
+ !!! resetting dpd to zero before doing the PP
+ !!! interactions. This allows contributions from
+ !!! different kernels to be accumulated. If needed,
+ !!! set it to zero before calling this routine the
+ !!! first time.
+ !-------------------------------------------------------------------------
+ ! Modules
+ !-------------------------------------------------------------------------
+ USE ppm_module_data
+ USE ppm_module_data_neighlist
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_error
+ USE ppm_module_neighlist_MkNeighIdx
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+ !-------------------------------------------------------------------------
+ ! Includes
+ !-------------------------------------------------------------------------
+#include "ppm_define.h"
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: xp
+ !!! particle co-ordinates
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle strengths used for interaction.
+ !!! Overloaded type: single,double
+ REAL(MK) , DIMENSION(:,:), INTENT( OUT) :: dpd
+ !!! Change of particle data (pdata) due to interaction.
+ !!! This is not initialized by this routine!
+ !!! Overloaded types: single,double.
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle strengths used for interaction.
+ !!! Overloaded types: single complex, double complex.
+ COMPLEX(MK), DIMENSION(:,:), INTENT( OUT) :: dpd
+ !!! Change of particle data (pdata) due to interaction.
+ !!! This is not initialized by this routine!
+ !!! Overloaded types: single complex, double complex.
+#endif
+ INTEGER , INTENT(IN ) :: Np
+ !!! number of particles on this proc.
+ INTEGER , INTENT(IN ) :: lda
+ !!! lading dimension of pdata.
+ INTEGER , INTENT(IN ) :: nsublist
+ !!! number of subdomains on the local processor
+#if __KERNEL == __INTERNAL
+ INTEGER , INTENT(IN ) :: kernel
+ !!! kernel to be used for PP interactions. To use ppm-internal
+ !!! kernels, specify one of:
+ !!!
+ !!! ---------------------------------------
+ !!! ppm_param_kerel_laplace2d_2p
+ !!! (2nd order Laplacian, polynomial in 2D)
+ !!! ppm_param_kerel_laplace3d_2p
+ !!! (2nd order Laplacian, polynomial in 3D)
+ !!! ---------------------------------------
+ !!!
+ !!! To use your own kernel function, pass the function pointer here.
+ !!! Your function should take one argument and return one value.
+ !!! The third possibility is to pass a lookup table with tabulated
+ !!! kernel values. Such a table can be created using
+ !!! ppm_comp_pp_mk_table.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc
+ !!! for description. Type can be single, double
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc
+ !!! for description. Type can be single complex or double complex.
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+ REAL(MK) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. Lookup table version.
+ !!! Pass dxtableinv (the inverse of the table spacing) as a scalar here.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(: ), INTENT(IN ) :: kernel
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(: ), INTENT(IN ) :: kernel
+#endif
+ !!! lookup table with tabulated kernel values. Such a table can
+ !!! be created using <<ppm_comp_pp_mk_table>>.
+#endif
+ LOGICAL , INTENT(IN ) :: lsymm
+ !!! Do symmetric PP interactions?
+ TYPE(ppm_type_ptr_to_clist), DIMENSION(:), INTENT(IN) :: clist
+ !!! Cell list as a list of ptr_to_clist.
+ !!!
+ !!! particle index list of isub: `clist(isub)%lpdx(:)` +
+ !!! pointer to first particle in ibox of isub: ` clist(isub)%lhbx(ibox)`
+ INTEGER , DIMENSION(:,:), INTENT(IN ) :: Nm
+ !!! number of cells in x,y,(z) directions (including the ghosts cells)
+ !!! in each subdomain.
+ !!!
+ !!! 1st index: direction. second index: subid.
+ REAL(MK) , INTENT(IN ) :: cutoff2
+ !!! Squared PP interaction cutoff. Should be .LE. cell size.
+ INTEGER , INTENT( OUT) :: info
+ !!! Returns status, 0 upon success
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ ! counters
+ INTEGER :: i,idom,ibox,jbox,idx
+ INTEGER :: ipart,jpart,ip,jp
+ INTEGER :: cbox,iinter,j,k,ispec
+ ! coordinate differences
+ REAL(MK) :: dx,dy,dz
+ REAL(MK) :: factor,factor2
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) :: summ,summ2
+ REAL(MK) :: eta,dm
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK) :: summ,summ2
+ COMPLEX(MK) :: eta,dm
+#endif
+ ! square of inter particle distance
+ REAL(MK) :: dij,dij2,dij4,dij5
+ ! start and end particle in a box
+ INTEGER :: istart,iend,jstart,jend
+ ! cell neighbor lists
+ INTEGER, DIMENSION(:,:), POINTER :: inp,jnp
+ ! number of interactions for each cell
+ INTEGER, SAVE :: nnp
+ ! cell offsets for box index
+ INTEGER :: n1,n2,nz
+ REAL(MK) :: t0
+ !-------------------------------------------------------------------------
+ ! Externals
+ !-------------------------------------------------------------------------
+#if __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
+ REAL(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
+ COMPLEX(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#endif
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Initialise
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_comp_pp_cell',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Check arguments.
+ !-------------------------------------------------------------------------
+ IF (ppm_debug .GT. 0) THEN
+ IF (.NOT. ppm_initialized) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_ppm_noinit,'ppm_comp_pp_cell', &
+ & 'Please call ppm_init first!',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (Np .LT. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
+ & 'Np must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (lda .LT. 1) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
+ & 'lda must be >0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (nsublist .LT. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
+ & 'nsublist must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (cutoff2 .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
+ & 'cutoff2 must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#if __KERNEL == __LOOKUP_TABLE
+ IF (kpar .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
+ & 'kpar (dxtableinv) must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#endif
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Build interaction index lists
+ !-------------------------------------------------------------------------
+ CALL ppm_neighlist_MkNeighIdx(lsymm,inp,jnp,nnp,info)
+
+ !-------------------------------------------------------------------------
+ ! PARTICLE-PARTICLE INTERACTIONS using symmetry
+ !-------------------------------------------------------------------------
+ IF (lsymm) THEN
+ DO idom=1,nsublist
+ n1 = Nm(1,idom)
+ n2 = Nm(1,idom)*Nm(2,idom)
+ nz = Nm(3,idom)
+ IF (ppm_dim .EQ. 2) THEN
+ n2 = 0
+ nz = 2
+ ENDIF
+ ! loop over all REAL cells (the -2 in the end does this)
+ DO k=0,nz-2
+ DO j=0,Nm(2,idom)-2
+ DO i=0,Nm(1,idom)-2
+ ! index of the center box
+ cbox = i + 1 + n1*j + n2*k
+ ! loop over all box-box interactions
+ DO iinter=1,nnp
+ ! determine box indices for this interaction
+ ibox = cbox+(inp(1,iinter)+n1*inp(2,iinter)+ &
+ & n2*inp(3,iinter))
+ jbox = cbox+(jnp(1,iinter)+n1*jnp(2,iinter)+ &
+ & n2*jnp(3,iinter))
+ !-------------------------------------------------
+ ! Read indices and check if box is empty
+ !-------------------------------------------------
+ istart = clist(idom)%lhbx(ibox)
+ iend = clist(idom)%lhbx(ibox+1)-1
+ IF (iend .LT. istart) CYCLE
+ !-------------------------------------------------
+ ! Within the box itself use symmetry and avoid
+ ! adding the particle itself to its own list
+ !-------------------------------------------------
+ IF (ibox .EQ. jbox) THEN
+ DO ipart=istart,iend
+ ip = clist(idom)%lpdx(ipart)
+ IF (lda .EQ. 1) THEN
+ summ = 0.0_MK
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=(ipart+1),iend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp here... and
+ ! vice versa to use symmetry.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ summ = summ + dm
+ dpd(1,jp) = dpd(1,jp) - dm
+ ENDDO
+ dpd(1,ip) = dpd(1,ip) + summ
+ ELSEIF (lda .EQ. 2) THEN
+ summ = 0.0_MK
+ summ2 = 0.0_MK
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=(ipart+1),iend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp here... and
+ ! vice versa to use symmetry.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ summ = summ + dm
+ dpd(1,jp) = dpd(1,jp) - dm
+ dm = eta*(pdata(2,jp) - &
+ & pdata(2,ip))
+ summ2 = summ2 + dm
+ dpd(2,jp) = dpd(2,jp) - dm
+ ENDDO
+ dpd(1,ip) = dpd(1,ip) + summ
+ dpd(2,ip) = dpd(2,ip) + summ2
+ ELSE
+ DO jpart=(ipart+1),iend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp here... and
+ ! vice versa to use symmetry.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata(ispec,jp)- &
+ & pdata(ispec,ip))
+ dpd(ispec,ip)=dpd(ispec,ip)+dm
+ dpd(ispec,jp)=dpd(ispec,jp)-dm
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ !-------------------------------------------------
+ ! For the other boxes check all particles
+ !-------------------------------------------------
+ ELSE
+ ! get pointers to first and last particle
+ jstart = clist(idom)%lhbx(jbox)
+ jend = clist(idom)%lhbx(jbox+1)-1
+ ! skip this iinter if other box is empty
+ IF (jend .LT. jstart) CYCLE
+ ! loop over all particles inside this cell
+ DO ipart=istart,iend
+ ip = clist(idom)%lpdx(ipart)
+ IF (lda .EQ. 1) THEN
+ summ = 0.0_MK
+ ! check against all particles
+ ! in the other cell
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=jstart,jend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp here... and
+ ! vice versa to use symmetry.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ summ = summ +dm
+ dpd(1,jp) = dpd(1,jp) - dm
+ ENDDO
+ dpd(1,ip) = dpd(1,ip) + summ
+ ELSEIF (lda .EQ. 2) THEN
+ summ = 0.0_MK
+ summ2 = 0.0_MK
+ ! check against all particles
+ ! in the other cell
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=jstart,jend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp here... and
+ ! vice versa to use symmetry.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ summ = summ + dm
+ dpd(1,jp) = dpd(1,jp) - dm
+ dm = eta*(pdata(2,jp) - &
+ & pdata(2,ip))
+ summ2 = summ2 + dm
+ dpd(2,jp) = dpd(2,jp) - dm
+ ENDDO
+ dpd(1,ip) = dpd(1,ip) + summ
+ dpd(2,ip) = dpd(2,ip) + summ2
+ ELSE
+ ! check against all particles
+ ! in the other cell
+ DO jpart=jstart,jend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp here... and
+ ! vice versa to use symmetry.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata(ispec,jp)- &
+ & pdata(ispec,ip))
+ dpd(ispec,ip)=dpd(ispec,ip)+dm
+ dpd(ispec,jp)=dpd(ispec,jp)-dm
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ ENDIF ! ibox .EQ. jbox
+ ENDDO ! iinter
+ ENDDO ! i
+ ENDDO ! j
+ ENDDO ! k
+ ENDDO ! idom
+
+ !-------------------------------------------------------------------------
+ ! PARTICLE-PARTICLE INTERACTIONS not using symmetry
+ !-------------------------------------------------------------------------
+ ELSE
+ DO idom=1,nsublist
+ n1 = Nm(1,idom)
+ n2 = Nm(1,idom)*Nm(2,idom)
+ nz = Nm(3,idom)
+ IF (ppm_dim .EQ. 2) THEN
+ n2 = 0
+ nz = 2
+ ENDIF
+ ! loop over all REAL cells (the -2 in the end does this)
+ DO k=1,nz-2
+ DO j=1,Nm(2,idom)-2
+ DO i=1,Nm(1,idom)-2
+ ! index of the center box
+ cbox = i + 1 + n1*j + n2*k
+ ! loop over all box-box interactions
+ DO iinter=1,nnp
+ ! determine box indices for this interaction
+ ibox = cbox+(inp(1,iinter)+n1*inp(2,iinter)+ &
+ & n2*inp(3,iinter))
+ jbox = cbox+(jnp(1,iinter)+n1*jnp(2,iinter)+ &
+ & n2*jnp(3,iinter))
+ !-------------------------------------------------
+ ! Read indices and check if box is empty
+ !-------------------------------------------------
+ istart = clist(idom)%lhbx(ibox)
+ iend = clist(idom)%lhbx(ibox+1)-1
+ IF (iend .LT. istart) CYCLE
+ !-------------------------------------------------
+ ! Do all interactions within the box itself
+ !-------------------------------------------------
+ IF (ibox .EQ. jbox) THEN
+ DO ipart=istart,iend
+ ip = clist(idom)%lpdx(ipart)
+ IF (lda .EQ. 1) THEN
+ DO jpart=istart,iend
+ jp = clist(idom)%lpdx(jpart)
+ ! No particle interacts with
+ ! itself
+ IF (ip .EQ. jp) CYCLE
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ dpd(1,ip) = dpd(1,ip) + dm
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO jpart=istart,iend
+ jp = clist(idom)%lpdx(jpart)
+ ! No particle interacts with
+ ! itself
+ IF (ip .EQ. jp) CYCLE
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ dpd(1,ip) = dpd(1,ip) + dm
+ dm = eta*(pdata(2,jp) - &
+ & pdata(2,ip))
+ dpd(2,ip) = dpd(2,ip) + dm
+ ENDDO
+ ELSE
+ DO jpart=istart,iend
+ jp = clist(idom)%lpdx(jpart)
+ ! No particle interacts with
+ ! itself
+ IF (ip .EQ. jp) CYCLE
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata(ispec,jp)- &
+ & pdata(ispec,ip))
+ dpd(ispec,ip)=dpd(ispec,ip)+dm
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ !-------------------------------------------------
+ ! Do interactions with all neighboring boxes
+ !-------------------------------------------------
+ ELSE
+ ! get pointers to first and last particle
+ jstart = clist(idom)%lhbx(jbox)
+ jend = clist(idom)%lhbx(jbox+1)-1
+ ! skip if empty
+ IF (jend .LT. jstart) CYCLE
+ ! loop over all particles inside this cell
+ DO ipart=istart,iend
+ ip = clist(idom)%lpdx(ipart)
+ IF (lda .EQ. 1) THEN
+ ! check against all particles
+ ! in the other cell
+ DO jpart=jstart,jend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ dpd(1,ip) = dpd(1,ip) + dm
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ ! check against all particles
+ ! in the other cell
+ DO jpart=jstart,jend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - &
+ & pdata(1,ip))
+ dpd(1,ip) = dpd(1,ip) + dm
+ dm = eta*(pdata(2,jp) - &
+ & pdata(2,ip))
+ dpd(2,ip) = dpd(2,ip) + dm
+ ENDDO
+ ELSE
+ ! check against all particles
+ ! in the other cell
+ DO jpart=jstart,jend
+ jp = clist(idom)%lpdx(jpart)
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx)+(dy*dy)+(dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx)+(dy*dy)
+ ENDIF
+ IF (dij .GT. cutoff2) CYCLE
+ !---------------------------------
+ ! Particle ip interacts with
+ ! particle jp.
+ !---------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata(ispec,jp)- &
+ & pdata(ispec,ip))
+ dpd(ispec,ip)=dpd(ispec,ip)+dm
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ ENDIF ! ibox .EQ. jbox
+ ENDDO ! iinter
+ ENDDO ! i
+ ENDDO ! j
+ ENDDO ! k
+ ENDDO ! idom
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_comp_pp_cell',t0,info)
+ RETURN
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_cell_si
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_cell_di
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_cell_sci
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_cell_dci
+#endif
+
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_cell_su
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_cell_du
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_cell_scu
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_cell_dcu
+#endif
+
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_cell_st
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_cell_dt
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_cell_sct
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_cell_dct
+#endif
+#endif
diff --git a/src/ppm_comp_pp_correct.f b/src/ppm_comp_pp_correct.f
new file mode 100644
index 0000000000000000000000000000000000000000..1c256fa5bbf94bb0b081b6fcd18929e84b5259aa
--- /dev/null
+++ b/src/ppm_comp_pp_correct.f
@@ -0,0 +1,343 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_comp_pp_correct
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_correct_si(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_correct_di(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_correct_sci(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_correct_dci(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#endif
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_correct_su(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_correct_du(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_correct_scu(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_correct_dcu(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_correct_st(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_correct_dt(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_correct_sct(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_correct_dct(dh,cutoff,kernel,eps,kpar,targt, &
+ & kappa,info)
+#endif
+#endif
+ !!! Computes the discretisation correction factor
+ !!! needed to ensure that the discrete second-order
+ !!! moment of the kernel function is equal to some
+ !!! value. The kernel is evaluated on a regular mesh
+ !!! here, so this factor can only be used if the
+ !!! particles are on a regular grid (e.g. after
+ !!! remeshing).
+ !-------------------------------------------------------------------------
+ ! Modules
+ !-------------------------------------------------------------------------
+ USE ppm_module_data
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_error
+ USE ppm_module_write
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+ !-------------------------------------------------------------------------
+ ! Includes
+ !-------------------------------------------------------------------------
+#include "ppm_define.h"
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ REAL(MK), DIMENSION(:), INTENT(IN ) :: dh
+ !!! grid spacings in all directions 1..ppm_dim
+ REAL(MK) , INTENT(IN ) :: cutoff
+ !!! PP interaction cutoff used. If no cutoff is used, pass a negative
+ !!! value here.
+ REAL(MK) , INTENT(IN ) :: eps
+ !!! Kernel size (width)
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , INTENT(IN ) :: targt
+ !!! target value for the second order kernel moment. Overloaded types:
+ !!! single,double
+ REAL(MK) , INTENT( OUT) :: kappa
+ !!! computed discretisation correction factor. Every kernel evaluation has
+ !!! to be multiplied with this in order to correct for discretisation
+ !!! effects. Overloaded types: single,double
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK) , INTENT(IN ) :: targt
+ !!! target value for the second order kernel moment. Overloaded types:
+ !!! single complex, double complex
+ COMPLEX(MK) , INTENT( OUT) :: kappa
+ !!! computed discretisation correction factor. Every kernel evaluation has
+ !!! to be multiplied with this in order to correct for discretisation
+ !!! effects. Overloaded types: single complex, double complex
+#endif
+#if __KERNEL == __INTERNAL
+ INTEGER , INTENT(IN ) :: kernel
+ !!! kernel for which to compute the correction. To use ppm-internal
+ !!! kernels, specify one of:
+ !!!
+ !!! ---------------------------------------
+ !!! ppm_param_kerel_laplace2d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 2D)
+ !!! ppm_param_kerel_laplace3d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 3D)
+ !!! ---------------------------------------
+ !!!
+ !!! To use your own kernel function, pass the function pointer here. Your
+ !!! function should take one argument and return one value.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:), INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single or complex.
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:), INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single complex or double complex.
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+ REAL(MK) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Pass dxtableinv (the inverse of the table spacing) as
+ !!! a scalar here.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(: ), INTENT(IN ) :: kernel
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(: ), INTENT(IN ) :: kernel
+#endif
+ !!! Lookup table with tabulated kernel values.
+ !!! Such a table can be created using <<ppm_comp_pp_mk_table>>.
+#endif
+ INTEGER , INTENT( OUT) :: info
+ !!! Returns 0 upon success
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(MK) :: t0,dij2,dij4,dij5
+ REAL(MK) :: hi,dx,dy,dz,dij,factor,factor2
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) :: eta,summ
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK) :: eta,summ
+#endif
+ INTEGER :: lhx,lhy,lhz,i,j,k,idx
+ CHARACTER(LEN=ppm_char) :: cbuf
+ !-------------------------------------------------------------------------
+ ! Externals
+ !-------------------------------------------------------------------------
+#if __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:), INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
+ REAL(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:), INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
+ COMPLEX(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#endif
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Initialise
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_comp_pp_correct',t0,info)
+ kappa = 1.0_MK
+
+ !-------------------------------------------------------------------------
+ ! Check arguments.
+ !-------------------------------------------------------------------------
+ IF (ppm_debug .GT. 0) THEN
+ IF (.NOT. ppm_initialized) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_ppm_noinit,'ppm_comp_pp_correct', &
+ & 'Please call ppm_init first!',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (eps .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
+ & 'eps must be > 0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (SIZE(dh,1) .LT. ppm_dim) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
+ & 'dh must be at least of length ppm_dim',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ DO i=1,ppm_dim
+ IF (dh(i) .LE. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
+ & 'dh must be >0 in all directions',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ ENDDO
+#if __KERNEL == __LOOKUP_TABLE
+ IF (kpar .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
+ & 'kpar (dxtableinv) must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#endif
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Integration bounds
+ !-------------------------------------------------------------------------
+ IF (cutoff .GE. 0.0_MK) THEN
+ hi = cutoff
+ ELSE
+ hi = 50.0_MK
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Summation bounds
+ !-------------------------------------------------------------------------
+ lhx = INT(hi/dh(1))
+ lhy = INT(hi/dh(2))
+ IF (ppm_dim .GT. 2) lhz = INT(hi/dh(3))
+
+ !-------------------------------------------------------------------------
+ ! Discrete sum
+ !-------------------------------------------------------------------------
+ summ = 0.0_MK
+ IF (ppm_dim .GT. 2) THEN
+ DO i=-lhx,lhx
+ dx = REAL(i,MK)*dh(1)
+ DO j=-lhy,lhy
+ dy = REAL(j,MK)*dh(2)
+ DO k=-lhz,lhz
+ dz = REAL(k,MK)*dh(3)
+ dij = dx*dx + dy*dy + dz*dz
+#include "ppm_comp_pp_kernels.inc"
+ ! second order moment
+ eta = dij*eta
+ summ = summ + eta
+ ENDDO
+ ENDDO
+ ENDDO
+ summ = summ*dh(1)*dh(2)*dh(3)
+ ELSE
+ DO i=-lhx,lhx
+ dx = REAL(i,MK)*dh(1)
+ DO j=-lhy,lhy
+ dy = REAL(j,MK)*dh(2)
+ dij = dx*dx + dy*dy
+#include "ppm_comp_pp_kernels.inc"
+ ! second order moment
+ eta = dij*eta
+ summ = summ + eta
+ ENDDO
+ ENDDO
+ summ = summ*dh(1)*dh(2)
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Correction factor
+ !-------------------------------------------------------------------------
+#if __KERNEL != __LOOKUP_TABLE
+ kappa = (REAL(ppm_dim,MK)*targt*eps)/summ
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Debug output
+ !-------------------------------------------------------------------------
+ IF (ppm_debug .GT. 1) THEN
+ WRITE(cbuf,'(A,F16.10)') 'Discretization correction factor: ',kappa
+ CALL ppm_write(ppm_rank,'ppm_comp_pp_correct',cbuf,info)
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_comp_pp_correct',t0,info)
+ RETURN
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_correct_si
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_correct_di
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_correct_sci
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_correct_dci
+#endif
+
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_correct_su
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_correct_du
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_correct_scu
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_correct_dcu
+#endif
+
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_correct_st
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_correct_dt
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_correct_sct
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_correct_dct
+#endif
+#endif
diff --git a/src/ppm_comp_pp_doring.f b/src/ppm_comp_pp_doring.f
new file mode 100644
index 0000000000000000000000000000000000000000..a5b3bc05257ce7340c9a456e0087890d284ac1d1
--- /dev/null
+++ b/src/ppm_comp_pp_doring.f
@@ -0,0 +1,611 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_comp_pp_doring
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_doring_si(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_doring_di(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_doring_sci(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_doring_dci(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#endif
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_doring_su(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_doring_du(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_doring_scu(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_doring_dcu(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_doring_st(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_doring_dt(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_doring_sct(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_doring_dct(xp,pdata,dpd,Mpart,xp2,pdata2,dpd2, &
+ & Lpart,lda,lsymm,kernel,kpar,mode,info)
+#endif
+#endif
+ !!! Computes direct particle-particle interactions using the
+ !!! ring topology (N**2). This is called by ppm_comp_pp_ring and
+ !!! not by the user directly.
+ !!!
+ !!! [NOTE]
+ !!! The loops over lda are explicitly unrolled for the cases
+ !!! lda.EQ.1 and lda.EQ.2 in order to allow vectorization in these cases.
+ !-------------------------------------------------------------------------
+ ! Modules
+ !-------------------------------------------------------------------------
+ USE ppm_module_data
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_error
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+ !-------------------------------------------------------------------------
+ ! Includes
+ !-------------------------------------------------------------------------
+#include "ppm_define.h"
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: xp
+ !!! particle co-ordinates [Group 1]
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: xp2
+ !!! particle co-ordinates [Group 2]
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle data used for interaction (e.g. vorticity, strength,...)
+ !!! [Group 1]. Overloaded types: single,double
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: pdata2
+ !!! particle data used for interaction (e.g. vorticity, strength,...)
+ !!! [Group 2]. Overloaded types: single,double
+ REAL(MK) , DIMENSION(:,:), INTENT(INOUT) :: dpd
+ !!! Change of particle data due to interaction [Group 1]
+ !!! Overloaded types: single,double
+ REAL(MK) , DIMENSION(:,:), INTENT(INOUT) :: dpd2
+ !!! Change of particle data due to interaction [Group 2]
+ !!! Overloaded types: single,double
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle data used for interaction (e.g. vorticity, strength,...)
+ !!! [Group 1]. Overloaded types: single complex, double complex.
+ COMPLEX(MK), DIMENSION(:,:), INTENT(IN ) :: pdata2
+ !!! particle data used for interaction (e.g. vorticity, strength,...)
+ !!! [Group 2]. Overloaded types: single complex, double complex.
+ COMPLEX(MK), DIMENSION(:,:), INTENT(INOUT) :: dpd
+ !!! Change of particle data due to interaction [Group 1]
+ !!! Overloaded types: single complex,double complex.
+ COMPLEX(MK), DIMENSION(:,:), INTENT(INOUT) :: dpd2
+ !!! Change of particle data due to interaction [Group 2]
+ !!! Overloaded types: single complex,double complex.
+#endif
+ INTEGER , INTENT(IN ) :: Mpart
+ !!! number of particles [Group 1]
+ INTEGER , INTENT(IN ) :: Lpart
+ !!! number of particles [Group 2]
+ INTEGER , INTENT(IN ) :: mode
+ !!! Whether the two groups are the same or not:
+ !!!
+ !!! * 0 not the same group
+ !!! * 1 the same group
+ INTEGER , INTENT(IN ) :: lda
+ !!! leading dimension of pdata,pdata2.
+#if __KERNEL == __INTERNAL
+ INTEGER , INTENT(IN ) :: kernel
+ !!! kernel for which to compute the correction. To use ppm-internal
+ !!! kernels, specify one of:
+ !!!
+ !!! ---------------------------------------
+ !!! ppm_param_kerel_laplace2d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 2D)
+ !!! ppm_param_kerel_laplace3d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 3D)
+ !!! ---------------------------------------
+ !!!
+ !!! To use your own kernel function, pass the function pointer here. Your
+ !!! function should take one argument and return one value.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single or complex.
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single complex or double complex.
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+ REAL(MK) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Pass dxtableinv (the inverse of the table spacing) as
+ !!! a scalar here.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(: ), INTENT(IN ) :: kernel
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(: ), INTENT(IN ) :: kernel
+#endif
+ !!! Lookup table with tabulated kernel values.
+ !!! Such a table can be created using <<ppm_comp_pp_mk_table>>.
+#endif
+ LOGICAL , INTENT(IN ) :: lsymm
+ !!! Whether to use symmetry or not
+ INTEGER , INTENT( OUT) :: info
+ !!! Returns status, 0 upon success
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ ! counters
+ INTEGER :: i,j,ispec,idx
+ ! coordinate differences
+ REAL(MK) :: dx,dy,dz
+ REAL(MK) :: factor,factor2
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) :: summ,summ2
+ REAL(MK) :: eta,dm
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK) :: summ,summ2
+ COMPLEX(MK) :: eta,dm
+#endif
+ ! square of inter particle distance
+ REAL(MK) :: dij,dij2,dij4,dij5
+ REAL(MK) :: t0
+ !-------------------------------------------------------------------------
+ ! Externals
+ !-------------------------------------------------------------------------
+#if __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
+ REAL(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
+ COMPLEX(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#endif
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Initialise
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_comp_pp_doring',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Check arguments
+ !-------------------------------------------------------------------------
+ IF (ppm_debug .GT. 0) THEN
+ IF (Mpart .LT. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_doring', &
+ & 'Np must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (Lpart .LT. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_doring', &
+ & 'Np must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (lda .LT. 1) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_doring', &
+ & 'lda must be >0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF ((mode .NE. 0) .AND. (mode .NE. 1)) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_doring', &
+ & 'mode must be either 0 or 1',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#if __KERNEL == __LOOKUP_TABLE
+ IF (kpar .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_doring', &
+ & 'kpar (dxtableinv) must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#endif
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Check if we are computing the selfinteraction (i.e. Group1 .EQ. Group2)
+ !-------------------------------------------------------------------------
+ IF (mode .EQ. 1) THEN
+ IF (lsymm) THEN
+ !-----------------------------------------------------------------
+ ! SYMMETRY
+ !-----------------------------------------------------------------
+ IF (lda .EQ. 1) THEN
+ DO i = 1,Mpart
+ summ = 0.0_MK
+ DO j = i+1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j and vice versa
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ summ = summ + dm
+ dpd2(1,j) = dpd2(1,j) - dm
+ ENDDO
+ dpd(1,i) = dpd(1,i) + summ
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO i = 1,Mpart
+ summ = 0.0_MK
+ summ2 = 0.0_MK
+ DO j = i+1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j and vice versa
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ summ = summ + dm
+ dpd2(1,j) = dpd2(1,j) - dm
+ dm = eta*(pdata2(2,j) - pdata(2,i))
+ summ2 = summ2 + dm
+ dpd2(2,j) = dpd2(2,j) - dm
+ ENDDO
+ dpd(1,i) = dpd(1,i) + summ
+ dpd(2,i) = dpd(2,i) + summ2
+ ENDDO
+ ELSE
+ DO i = 1,Mpart
+ DO j = i+1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j and vice versa
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata2(ispec,j) - pdata(ispec,i))
+ dpd(ispec,i) = dpd(ispec,i) + dm
+ dpd2(ispec,j) = dpd2(ispec,j) - dm
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ELSE
+ !-----------------------------------------------------------------
+ ! NO SYMMETRY
+ !-----------------------------------------------------------------
+ IF (lda .EQ. 1) THEN
+ DO i = 1,Mpart
+ DO j = 1,Lpart
+ IF (i .EQ. j) CYCLE
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ dpd(1,i) = dpd(1,i) + dm
+ ENDDO
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO i = 1,Mpart
+ DO j = 1,Lpart
+ IF (i .EQ. j) CYCLE
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ dpd(1,i) = dpd(1,i) + dm
+ dm = eta*(pdata2(2,j) - pdata(2,i))
+ dpd(2,i) = dpd(2,i) + dm
+ ENDDO
+ ENDDO
+ ELSE
+ DO i = 1,Mpart
+ DO j = 1,Lpart
+ IF (i .EQ. j) CYCLE
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata2(ispec,j) - pdata(ispec,i))
+ dpd(ispec,i) = dpd(ispec,i) + dm
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! non-self interaction (i.e. Group1 .NE. Group2)
+ !-------------------------------------------------------------------------
+ ELSE
+ IF (lsymm) THEN
+ !-----------------------------------------------------------------
+ ! SYMMETRY
+ !-----------------------------------------------------------------
+ IF (lda .EQ. 1) THEN
+ DO i = 1,Mpart
+ summ = 0.0_MK
+ DO j = 1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j and vice versa
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ summ = summ + dm
+ dpd2(1,j) = dpd2(1,j) - dm
+ ENDDO
+ dpd(1,i) = dpd(1,i) + summ
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO i = 1,Mpart
+ summ = 0.0_MK
+ summ2 = 0.0_MK
+ DO j = 1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j and vice versa
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ summ = summ + dm
+ dpd2(1,j) = dpd2(1,j) - dm
+ dm = eta*(pdata2(2,j) - pdata(2,i))
+ summ2 = summ2 + dm
+ dpd2(2,j) = dpd2(2,j) - dm
+ ENDDO
+ dpd(1,i) = dpd(1,i) + summ
+ dpd(2,i) = dpd(2,i) + summ2
+ ENDDO
+ ELSE
+ DO i = 1,Mpart
+ DO j = 1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j and vice versa
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata2(ispec,j) - pdata(ispec,i))
+ dpd(ispec,i) = dpd(ispec,i) + dm
+ dpd2(ispec,j) = dpd2(ispec,j) - dm
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ELSE
+ !-----------------------------------------------------------------
+ ! NO SYMMETRY
+ !-----------------------------------------------------------------
+ IF (lda .EQ. 1) THEN
+ DO i = 1,Mpart
+ DO j = 1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ dpd(1,i) = dpd(1,i) + dm
+ ENDDO
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO i = 1,Mpart
+ DO j = 1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata2(1,j) - pdata(1,i))
+ dpd(1,i) = dpd(1,i) + dm
+ dm = eta*(pdata2(2,j) - pdata(2,i))
+ dpd(2,i) = dpd(2,i) + dm
+ ENDDO
+ ENDDO
+ ELSE
+ DO i = 1,Mpart
+ DO j = 1,Lpart
+ dx = xp(1,i) - xp2(1,j)
+ dy = xp(2,i) - xp2(2,j)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,i) - xp2(3,j)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !--------------------------------------------------------
+ ! Particle i interacts with j
+ !--------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata2(ispec,j) - pdata(ispec,i))
+ dpd(ispec,i) = dpd(ispec,i) + dm
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDIF
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_comp_pp_doring',t0,info)
+ RETURN
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_doring_si
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_doring_di
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_doring_sci
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_doring_dci
+#endif
+
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_doring_su
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_doring_du
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_doring_scu
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_doring_dcu
+#endif
+
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_doring_st
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_doring_dt
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_doring_sct
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_doring_dct
+#endif
+#endif
diff --git a/src/ppm_comp_pp_mk_table.f b/src/ppm_comp_pp_mk_table.f
new file mode 100644
index 0000000000000000000000000000000000000000..b4963cbbe2897889020880053f2286c6ffa06bb9
--- /dev/null
+++ b/src/ppm_comp_pp_mk_table.f
@@ -0,0 +1,247 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_comp_pp_mk_table
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_mk_table_si(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_mk_table_di(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_mk_table_sci(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_mk_table_dci(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#endif
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_mk_table_su(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_mk_table_du(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_mk_table_scu(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_mk_table_dcu(kernel,kpar,ntable,cutoff2, &
+ & ktab,dxtableinv,info)
+#endif
+#endif
+ !!! Creates a lookup table for a PP interaction kernel.
+ !!! The result is returned in the variables ktab with
+ !!! indices running from 0 to ntable-1. The table
+ !!! contains kernel values eta as a function of the
+ !!! _squared_ distance x**2.
+ !-------------------------------------------------------------------------
+ ! Modules
+ !-------------------------------------------------------------------------
+ USE ppm_module_data
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_error
+ USE ppm_module_alloc
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+ !-------------------------------------------------------------------------
+ ! Includes
+ !-------------------------------------------------------------------------
+#include "ppm_define.h"
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ INTEGER , INTENT(IN ) :: ntable
+ !!! Number of table entries to be created.
+ REAL(MK) , INTENT(IN ) :: cutoff2
+ !!! Square of the cutoff used for PP interactions. The kernel will be
+ !!! tabulated up to this value.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:), POINTER :: ktab
+ !!! Tabulated kernel values 0...ntable-1.
+ !!! Overloaded types: single,double.
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK) , DIMENSION(:), POINTER :: ktab
+ !!! Tabulated kernel values 0...ntable-1.
+ !!! Overloaded types: single complex,double complex.
+#endif
+#if __KERNEL == __INTERNAL
+ INTEGER , INTENT(IN ) :: kernel
+ !!! kernel for which to compute the correction. To use ppm-internal
+ !!! kernels, specify one of:
+ !!!
+ !!! ---------------------------------------
+ !!! ppm_param_kerel_laplace2d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 2D)
+ !!! ppm_param_kerel_laplace3d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 3D)
+ !!! ---------------------------------------
+ !!!
+ !!! To use your own kernel function, pass the function pointer here. Your
+ !!! function should take one argument and return one value.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+#endif
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single, double, single complex or double
+ !!! complex. Omit this argument when using a lookup table.
+#endif
+ REAL(MK) , INTENT( OUT) :: dxtableinv
+ !!! inverse of the dx (kernel evaluation locations) used to create
+ !!! the table. To use the table, look up the value
+ !!! `ktab(INT(dxtableinv*x))`.
+ INTEGER , INTENT( OUT) :: info
+ !!! Returns status, 0 upon success
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ INTEGER :: i,atable,iopt,idx
+ INTEGER, DIMENSION(1) :: ldl,ldu
+ REAL(MK), PARAMETER :: Rmin = 0.0_MK
+ REAL(MK) :: dxtable,Rmax,factor,factor2
+ REAL(MK) :: t0,dij,dij2,dij4,dij5,dx,dy,dz
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) :: eta
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK) :: eta
+#endif
+ !-------------------------------------------------------------------------
+ ! Externals
+ !-------------------------------------------------------------------------
+#if __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
+ REAL(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
+ COMPLEX(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#endif
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Initialise
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_comp_pp_mk_table',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Check arguments.
+ !-------------------------------------------------------------------------
+ IF (ppm_debug .GT. 0) THEN
+ IF (ntable .LT. 1) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_mk_table', &
+ & 'ntable must be >0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (cutoff2 .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_mk_table', &
+ & 'cutoff2 must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Compute Rmax
+ !-------------------------------------------------------------------------
+ Rmax = cutoff2
+
+ !-------------------------------------------------------------------------
+ ! Allocate memory for table lookup
+ !-------------------------------------------------------------------------
+ atable = ntable-1
+
+ iopt = ppm_param_alloc_fit
+ ldl(1) = 0
+ ldu(1) = atable
+ CALL ppm_alloc(ktab,ldl,ldu,iopt,info)
+ IF (info .NE. ppm_param_success) THEN
+ info = ppm_error_fatal
+ CALL ppm_error(ppm_err_alloc,'ppm_comp_pp_mk_table', &
+ & 'kernel lookup table KTAB',__LINE__,info)
+ GOTO 9999
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Compute values
+ !-------------------------------------------------------------------------
+ dxtable = (Rmax - Rmin)/REAL(atable,MK)
+ dxtableinv = 1.0_MK/dxtable
+
+ ! This hack is needed since some kernels use dx,dy,dz explicitly. DO
+ ! NOT MAKE TABLES FOR THOSE KERNELS!
+ dx = 0.0_MK
+ dy = 0.0_MK
+ dz = 0.0_MK
+ DO i=0,atable
+ dij = Rmin + REAL(i,MK)*dxtable
+#include "ppm_comp_pp_kernels.inc"
+ ktab(i) = eta
+ ENDDO
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_comp_pp_mk_table',t0,info)
+ RETURN
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_mk_table_si
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_mk_table_di
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_mk_table_sci
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_mk_table_dci
+#endif
+
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_mk_table_su
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_mk_table_du
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_mk_table_scu
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_mk_table_dcu
+#endif
+#endif
diff --git a/src/ppm_comp_pp_ring.f b/src/ppm_comp_pp_ring.f
new file mode 100644
index 0000000000000000000000000000000000000000..fcf8017e1248bcd7eae2eea7517e3ccbb4f101b9
--- /dev/null
+++ b/src/ppm_comp_pp_ring.f
@@ -0,0 +1,503 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_comp_pp_ring
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_ring_si(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_ring_di(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_ring_sci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_ring_dci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#endif
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_ring_su(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_ring_du(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_ring_scu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_ring_dcu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_ring_st(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_ring_dt(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_ring_sct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_ring_dct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & dpd,info)
+#endif
+#endif
+ !!! Subroutine which computes kernel interactions by
+ !!! direct global (N**2) particle-particle
+ !!! interactions using the ring topology.
+ !!!
+ !!! [NOTE]
+ !!! dpd needs to be allocated to proper size before
+ !!! calling this routine. Also, this routine is not
+ !!! resetting dpd to zero before doing the PP interactions.
+ !!! This allows contributions from different kernels to be
+ !!! accumulated. If needed, set it to zero before calling this
+ !!! routine the first time.
+ !-------------------------------------------------------------------------
+ ! Modules
+ !-------------------------------------------------------------------------
+ USE ppm_module_data
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_error
+ USE ppm_module_alloc
+ USE ppm_module_map_part
+ USE ppm_module_map_part_ring_shift
+ USE ppm_module_comp_pp_doring
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+ !-------------------------------------------------------------------------
+ ! Includes
+ !-------------------------------------------------------------------------
+#include "ppm_define.h"
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: xp
+ !!! particle positions
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle data (e.g. strengths).
+ !!! Overloaded types: single,double
+ REAL(MK) , DIMENSION(:,:), INTENT( OUT) :: dpd
+ !!! Change of particle data (pdata) due to interaction.
+ !!! Overloaded types: single,double.
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle data (e.g. strengths).
+ !!! Overloaded types: single complex,double complex.
+ COMPLEX(MK), DIMENSION(:,:), INTENT( OUT) :: dpd
+ !!! Change of particle data (pdata) due to interaction.
+ !!! Overloaded types: single complex,double complex.
+#endif
+ LOGICAL , INTENT(IN ) :: lsymm
+ !!! Whether to use symmetry or not:
+ !!!
+ !!! * `.FALSE.` PP interaction w/o symmetry
+ !!! * `.TRUE.` PP interaction w/ symmetry
+ INTEGER , INTENT(IN ) :: Np
+ !!! number of particles on local proc.
+ INTEGER , INTENT(IN ) :: lda
+ !!! leading dimension of pdata.
+#if __KERNEL == __INTERNAL
+ INTEGER , INTENT(IN ) :: kernel
+ !!! kernel for which to compute the correction. To use ppm-internal
+ !!! kernels, specify one of:
+ !!!
+ !!! ---------------------------------------
+ !!! ppm_param_kerel_laplace2d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 2D)
+ !!! ppm_param_kerel_laplace3d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 3D)
+ !!! ---------------------------------------
+ !!!
+ !!! To use your own kernel function, pass the function pointer here. Your
+ !!! function should take one argument and return one value.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single or complex.
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single complex or double complex.
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+ REAL(MK) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Pass dxtableinv (the inverse of the table spacing) as
+ !!! a scalar here.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(: ), INTENT(IN ) :: kernel
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(: ), INTENT(IN ) :: kernel
+#endif
+ !!! Lookup table with tabulated kernel values.
+ !!! Such a table can be created using <<ppm_comp_pp_mk_table>>.
+#endif
+ INTEGER , INTENT( OUT) :: info
+ !!! Returns status, 0 upon success
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ INTEGER :: i,j
+ INTEGER :: hops,isource,itarget
+ INTEGER :: ll,lu,rl,ru
+ INTEGER :: Lpart,iopt
+ INTEGER , DIMENSION(2) :: ldu
+ REAL(MK) , DIMENSION(:,:), POINTER :: xp2
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:,:), POINTER :: pdata2,dpd2
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:,:), POINTER :: pdata2,dpd2
+#endif
+ REAL(MK) :: t0
+ !-------------------------------------------------------------------------
+ ! Externals
+ !-------------------------------------------------------------------------
+#if __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
+ REAL(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
+ COMPLEX(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#endif
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Initialise
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_comp_pp_ring',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Check arguments.
+ !-------------------------------------------------------------------------
+ IF (ppm_debug .GT. 0) THEN
+ IF (.NOT. ppm_initialized) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_ppm_noinit,'ppm_comp_pp_ring', &
+ & 'Please call ppm_init first!',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (Np .LT. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_ring', &
+ & 'Np must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (lda .LT. 1) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_ring', &
+ & 'lda must be >0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#if __KERNEL == __LOOKUP_TABLE
+ IF (kpar .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_ring', &
+ & 'kpar (dxtableinv) must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#endif
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Allocate memory for the copy of the local particles
+ !-------------------------------------------------------------------------
+ iopt = ppm_param_alloc_fit
+ ldu(1) = ppm_dim
+ ldu(2) = Np
+ CALL ppm_alloc(xp2,ldu,iopt,info)
+ IF (info .NE. ppm_param_success) THEN
+ info = ppm_error_fatal
+ CALL ppm_error(ppm_err_alloc,'ppm_comp_pp_ring', &
+ & 'local copy of xp XP2',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ ldu(1) = lda
+ ldu(2) = Np
+ CALL ppm_alloc(pdata2,ldu,iopt,info)
+ IF (info .NE. ppm_param_success) THEN
+ info = ppm_error_fatal
+ CALL ppm_error(ppm_err_alloc,'ppm_comp_pp_ring', &
+ & 'local copy of pdata PDATA2',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ CALL ppm_alloc(dpd2,ldu,iopt,info)
+ IF (info .NE. ppm_param_success) THEN
+ info = ppm_error_fatal
+ CALL ppm_error(ppm_err_alloc,'ppm_comp_pp_ring', &
+ & 'local copy of dpd DPD2',__LINE__,info)
+ GOTO 9999
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Make a copy of all local particles including all additional info
+ ! that they carry
+ !-------------------------------------------------------------------------
+ Lpart = Np
+ xp2(1:ppm_dim,1:Np) = xp(1:ppm_dim,1:Np)
+ pdata2(1:lda,1:Np) = pdata(1:lda,1:Np)
+ dpd2(1:lda,1:Np) = dpd(1:lda,1:Np)
+
+ !-------------------------------------------------------------------------
+ ! Calculate the interactions of the local particles with themselves.
+ !-------------------------------------------------------------------------
+ CALL ppm_comp_pp_doring(xp,pdata,dpd,Np,xp2,pdata2,dpd2,Lpart,lda, &
+ & lsymm,kernel,kpar,1,info)
+ IF (info .NE. ppm_param_success) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_sub_failed,'ppm_comp_pp_ring', &
+ & 'Interaction calculation failed.',__LINE__,info)
+ GOTO 9999
+ ENDIF
+
+#ifdef __MPI
+ !-------------------------------------------------------------------------
+ ! Compute whom to send the copy and who to receive the copy from
+ !-------------------------------------------------------------------------
+ itarget = MOD(ppm_nproc+ppm_rank-1,ppm_nproc)
+ isource = MOD(ppm_rank+1,ppm_nproc)
+
+ !-------------------------------------------------------------------------
+ ! Compute how often we have to shift a copy of the particles to the
+ ! neighbour processor
+ !-------------------------------------------------------------------------
+ IF (lsymm) THEN
+ IF (MOD(ppm_nproc,2) .EQ. 0) THEN
+ hops = (ppm_nproc/2)-1
+ ELSE
+ hops = ((ppm_nproc+1)/2)-1
+ ENDIF
+ ELSE
+ hops = ppm_nproc-1
+ ENDIF
+
+ DO i = 1,hops
+ !TODO: check the argument quantity accuracy
+ CALL ppm_map_part_ring_shift(xp2,ppm_dim,Lpart,itarget,isource,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & pdata2,lda,Lpart,Lpart,ppm_param_map_push,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & dpd2,lda,Lpart,Lpart,ppm_param_map_push,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & xp2,ppm_dim,Lpart,Lpart,ppm_param_map_send,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & dpd2,lda,Lpart,Lpart,ppm_param_map_pop,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & pdata2,lda,Lpart,Lpart,ppm_param_map_pop,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & xp2,ppm_dim,Lpart,Lpart,ppm_param_map_pop,info)
+
+ CALL ppm_comp_pp_doring(xp,pdata,dpd,Np,xp2,pdata2,dpd2,Lpart, &
+ & lda,lsymm,kernel,kpar,0,info)
+ IF (info .NE. ppm_param_success) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_sub_failed,'ppm_comp_pp_ring', &
+ & 'Interaction calculation failed.',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ ENDDO
+
+ !-------------------------------------------------------------------------
+ ! If we have symmetry we only have send the particles half the round,
+ ! so we have to check the case of an even number of processors
+ !-------------------------------------------------------------------------
+ IF (lsymm .AND. (MOD(ppm_nproc,2) .EQ. 0)) THEN
+ !TODO: check the argument quantity accuracy
+ CALL ppm_map_part_ring_shift(xp2,ppm_dim,Lpart,itarget,isource,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & pdata2,lda,Lpart,Lpart,ppm_param_map_push,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & dpd2,lda,Lpart,Lpart,ppm_param_map_push,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & xp2,ppm_dim,Lpart,Lpart,ppm_param_map_send,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & dpd2,lda,Lpart,Lpart,ppm_param_map_pop,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & pdata2,lda,Lpart,Lpart,ppm_param_map_pop,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & xp2,ppm_dim,Lpart,Lpart,ppm_param_map_pop,info)
+
+ !---------------------------------------------------------------------
+ ! The processor with the higher ppm_rank computes the upper half of
+ ! the particles and the opposite processor the lower half
+ !---------------------------------------------------------------------
+ IF (ppm_rank .GT. hops) THEN
+ ll = (Np / 2) + 1
+ lu = Np
+ rl = 1
+ ru = Lpart
+ ELSE
+ ll = 1
+ lu = Np
+ rl = 1
+ ru = Lpart / 2
+ ENDIF
+
+ CALL ppm_comp_pp_doring(xp(:,ll:lu),pdata(:,ll:lu),dpd(:,ll:lu),&
+ & lu-ll+1,xp2(:,rl:ru),pdata2(:,rl:ru), &
+ & dpd2(:,rl:ru),ru-rl+1,lda,lsymm,kernel,kpar, &
+ & 0,info)
+ IF (info .NE. ppm_param_success) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_sub_failed,'ppm_comp_pp_ring', &
+ & 'Interaction calculation failed.',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Send the particles where they belong to
+ !-------------------------------------------------------------------------
+ IF (lsymm) THEN
+ itarget = MOD(ppm_rank+(ppm_nproc/2),ppm_nproc)
+ isource = MOD(ppm_nproc+ppm_rank-(ppm_nproc/2),ppm_nproc)
+ ENDIF
+
+ IF (itarget .NE. ppm_rank) THEN
+ !TODO: check the argument quantity accuracy
+ CALL ppm_map_part_ring_shift(xp2,ppm_dim,Lpart,itarget,isource,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & pdata2,lda,Lpart,Lpart,ppm_param_map_push,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & dpd2,lda,Lpart,Lpart,ppm_param_map_push,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & xp2,ppm_dim,Lpart,Lpart,ppm_param_map_send,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & dpd2,lda,Lpart,Lpart,ppm_param_map_pop,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & pdata2,lda,Lpart,Lpart,ppm_param_map_pop,info)
+ CALL ppm_map_part(ppm_param_topo_undefined,ppm_param_topo_undefined, &
+ & xp2,ppm_dim,Lpart,Lpart,ppm_param_map_pop,info)
+
+ IF (Lpart .NE. Np) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_part_lost,'ppm_comp_pp_ring', &
+ & 'Not all particles came back!',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ ENDIF
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Add the particle changes of the local and the long traveled copy
+ !-------------------------------------------------------------------------
+ IF (lda .EQ. 1) THEN
+ DO i = 1,Np
+ dpd(1,i) = dpd(1,i) + dpd2(1,i)
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO i = 1,Np
+ dpd(1,i) = dpd(1,i) + dpd2(1,i)
+ dpd(2,i) = dpd(2,i) + dpd2(2,i)
+ ENDDO
+ ELSE
+ DO i = 1,Np
+ DO j = 1,lda
+ dpd(j,i) = dpd(j,i) + dpd2(j,i)
+ ENDDO
+ ENDDO
+ ENDIF
+
+ 9999 CONTINUE
+
+ !-------------------------------------------------------------------------
+ ! Deallocate memory of copies
+ !-------------------------------------------------------------------------
+ iopt = ppm_param_dealloc
+ CALL ppm_alloc(dpd2,ldu,iopt,info)
+ IF (info .NE. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_dealloc,'ppm_comp_pp_ring', &
+ & 'local copy of dpd DPD2',__LINE__,info)
+ ENDIF
+ CALL ppm_alloc(pdata2,ldu,iopt,info)
+ IF (info .NE. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_dealloc,'ppm_comp_pp_ring', &
+ & 'local copy of pdata PDATA2',__LINE__,info)
+ ENDIF
+ CALL ppm_alloc(xp2,ldu,iopt,info)
+ IF (info .NE. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_dealloc,'ppm_comp_pp_ring', &
+ & 'local copy of xp XP2',__LINE__,info)
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ CALL substop('ppm_comp_pp_ring',t0,info)
+ RETURN
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_ring_si
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_ring_di
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_ring_sci
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_ring_dci
+#endif
+
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_ring_su
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_ring_du
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_ring_scu
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_ring_dcu
+#endif
+
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_ring_st
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_ring_dt
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_ring_sct
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_ring_dct
+#endif
+#endif
diff --git a/src/ppm_comp_pp_verlet.f b/src/ppm_comp_pp_verlet.f
new file mode 100644
index 0000000000000000000000000000000000000000..b611e4498547d772d1a234ecd1683e650a2dcb33
--- /dev/null
+++ b/src/ppm_comp_pp_verlet.f
@@ -0,0 +1,493 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_comp_pp_verlet
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_verlet_si(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_verlet_di(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_verlet_sci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_verlet_dci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#endif
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_verlet_su(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_verlet_du(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_verlet_scu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_verlet_dcu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_comp_pp_verlet_st(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_comp_pp_verlet_dt(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_verlet_sct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ SUBROUTINE ppm_comp_pp_verlet_dct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
+ & nvlist,vlist,dpd,info)
+#endif
+#endif
+ !!! Subroutine which computes kernel interactions by direct
+ !!! particle-particle interactions using Verlet lists.
+ !!!
+ !!! [NOTE]
+ !!! ======================================================================
+ !!! The loops for lda.EQ.1 and lda.EQ.2 are explicitly
+ !!! unfolded to allow vectorization in those cases. For
+ !!! lda.GT.2 this routine will not vectorize.
+ !!!
+ !!! Both nvlist and vlist can be created using
+ !!! ppm_neighlist_vlist. This should be done by the
+ !!! user program since there could be several sets of
+ !!! Verlet lists. This routine can then be called with
+ !!! the appropriate one. Do not forget to DEALLOCATE
+ !!! nvlist and vlist when they are not needed any more.
+ !!!
+ !!! dpd needs to be allocated to proper size before
+ !!! calling this routine. Also, this routine is not
+ !!! resetting dpd to zero before doing the PP
+ !!! interactions. This allows contributions from
+ !!! different kernels to be accumulated. If needed,
+ !!! set it to zero before calling this routine the
+ !!! first time.
+ !!! ======================================================================
+ !!!
+ !-------------------------------------------------------------------------
+ ! Modules
+ !-------------------------------------------------------------------------
+ USE ppm_module_data
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_error
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+ !-------------------------------------------------------------------------
+ ! Includes
+ !-------------------------------------------------------------------------
+#include "ppm_define.h"
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: xp
+ !!! particle coordinates
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle data (strengths).
+ !!! Overloaded types: single,double
+ REAL(MK) , DIMENSION(:,:), INTENT( OUT) :: dpd
+ !!! Change of particle data (strengths).
+ !!! Overloaded types: single,double
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:,:), INTENT(IN ) :: pdata
+ !!! particle data (strengths).
+ !!! Overloaded types: single complex,double complex
+ COMPLEX(MK), DIMENSION(:,:), INTENT( OUT) :: dpd
+ !!! Change of particle data (strengths).
+ !!! Overloaded types: single complex,double complex.
+#endif
+ LOGICAL , INTENT(IN ) :: lsymm
+ !!! using symmetry or not?
+ INTEGER , INTENT(IN ) :: Np
+ !!! number of particles on local proc.
+ INTEGER , INTENT(IN ) :: lda
+ !!! leading dimension of pdata.
+#if __KERNEL == __INTERNAL
+ INTEGER , INTENT(IN ) :: kernel
+ !!! kernel for which to compute the correction. To use ppm-internal
+ !!! kernels, specify one of:
+ !!!
+ !!! ---------------------------------------
+ !!! ppm_param_kerel_laplace2d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 2D)
+ !!! ppm_param_kerel_laplace3d_2p
+ !!! (2nd order Laplacian,
+ !!! polynomial in 3D)
+ !!! ---------------------------------------
+ !!!
+ !!! To use your own kernel function, pass the function pointer here. Your
+ !!! function should take one argument and return one value.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single or complex.
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Type can be single complex or double complex.
+#endif
+#elif __KERNEL == __LOOKUP_TABLE
+ REAL(MK) , INTENT(IN ) :: kpar
+ !!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
+ !!! description. Pass dxtableinv (the inverse of the table spacing) as
+ !!! a scalar here.
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(: ), INTENT(IN ) :: kernel
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(: ), INTENT(IN ) :: kernel
+#endif
+ !!! Lookup table with tabulated kernel values.
+ !!! Such a table can be created using <<ppm_comp_pp_mk_table>>.
+#endif
+ INTEGER , DIMENSION(:,:), INTENT(IN ) :: vlist
+ !!! Verlet lists of all particles. 1st index: particles to interact with
+ !!! (1..nvlist), 2nd: particle (1..Np).
+ INTEGER , DIMENSION( :), INTENT(IN ) :: nvlist
+ !!! length of the Verlet lists of all particles
+ INTEGER , INTENT( OUT) :: info
+ !!! Returns status, 0 upon success
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ INTEGER :: jpart,ip,jp,ispec,idx
+ REAL(MK) :: dx,dy,dz
+ REAL(MK) :: factor,factor2
+ REAL(MK) :: dij,dij2,dij4,dij5
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) :: summ,summ2
+ REAL(MK) :: eta,dm
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK) :: summ,summ2
+ COMPLEX(MK) :: eta,dm
+#endif
+ ! timer
+ REAL(MK) :: t0
+ CHARACTER(LEN=ppm_char) :: mesg
+ !-------------------------------------------------------------------------
+ ! Externals
+ !-------------------------------------------------------------------------
+#if __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
+ REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
+ REAL(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+ COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
+ INTERFACE
+ FUNCTION kernel(x,kpar)
+ IMPLICIT NONE
+#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0E0)
+#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
+ INTEGER, PARAMETER :: MK = KIND(1.0D0)
+#endif
+ REAL(MK), INTENT(IN) :: x
+ COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
+ COMPLEX(MK) :: kernel
+ END FUNCTION kernel
+ END INTERFACE
+#endif
+#endif
+
+ !-------------------------------------------------------------------------
+ ! Initialise
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_comp_pp_verlet',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Check arguments.
+ !-------------------------------------------------------------------------
+ IF (ppm_debug .GT. 0) THEN
+ IF (.NOT. ppm_initialized) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_ppm_noinit,'ppm_comp_pp_verlet', &
+ & 'Please call ppm_init first!',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (Np .LT. 0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_verlet', &
+ & 'Np must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+ IF (lda .LT. 1) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_verlet', &
+ & 'lda must be >0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#if __KERNEL == __LOOKUP_TABLE
+ IF (kpar .LT. 0.0_MK) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_argument,'ppm_comp_pp_verlet', &
+ & 'kpar (dxtableinv) must be >=0',__LINE__,info)
+ GOTO 9999
+ ENDIF
+#endif
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! PARTICLE-PARTICLE INTERACTIONS using symmetry
+ !-------------------------------------------------------------------------
+ IF (lsymm) THEN
+ IF (lda .EQ. 1) THEN
+ DO ip=1,Np
+ summ = 0.0_MK
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=1,nvlist(ip)
+ jp = vlist(jpart,ip)
+ !---------------------------------------------------------
+ ! Calculate the square of the distance between the two
+ ! particles. It will always be .LE. (cutoff+skin)**2 by
+ ! construction of the Verlet list.
+ !---------------------------------------------------------
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !---------------------------------------------------------
+ ! Particle ip interacts with jp and vice versa
+ !---------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - pdata(1,ip))
+ summ = summ + dm
+ dpd(1,jp) = dpd(1,jp) - dm
+ ENDDO
+ dpd(1,ip) = dpd(1,ip) + summ
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO ip=1,Np
+ summ = 0.0_MK
+ summ2 = 0.0_MK
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=1,nvlist(ip)
+ jp = vlist(jpart,ip)
+ !---------------------------------------------------------
+ ! Calculate the square of the distance between the two
+ ! particles. It will always be .LE. (cutoff+skin)**2 by
+ ! construction of the Verlet list.
+ !---------------------------------------------------------
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !---------------------------------------------------------
+ ! Particle ip interacts with jp and vice versa
+ !---------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - pdata(1,ip))
+ summ = summ + dm
+ dpd(1,jp) = dpd(1,jp) - dm
+ dm = eta*(pdata(2,jp) - pdata(2,ip))
+ summ2 = summ2 + dm
+ dpd(2,jp) = dpd(2,jp) - dm
+ ENDDO
+ dpd(1,ip) = dpd(1,ip) + summ
+ dpd(2,ip) = dpd(2,ip) + summ2
+ ENDDO
+ ELSE
+ DO ip=1,Np
+ DO jpart=1,nvlist(ip)
+ jp = vlist(jpart,ip)
+ !---------------------------------------------------------
+ ! Calculate the square of the distance between the two
+ ! particles. It will always be .LE. (cutoff+skin)**2 by
+ ! construction of the Verlet list.
+ !---------------------------------------------------------
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !---------------------------------------------------------
+ ! Particle ip interacts with jp and vice versa
+ !---------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata(ispec,jp) - pdata(ispec,ip))
+ dpd(ispec,ip) = dpd(ispec,ip) + dm
+ dpd(ispec,jp) = dpd(ispec,jp) - dm
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! PARTICLE-PARTICLE INTERACTIONS not using symmetry
+ !-------------------------------------------------------------------------
+ ELSE
+ IF (lda .EQ. 1) THEN
+ DO ip=1,Np
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=1,nvlist(ip)
+ jp = vlist(jpart,ip)
+ !---------------------------------------------------------
+ ! Calculate the square of the distance between the two
+ ! particles. It will always be .LE. (cutoff+skin)**2 by
+ ! construction of the Verlet list.
+ !---------------------------------------------------------
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !---------------------------------------------------------
+ ! Particle ip interacts with jp
+ !---------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - pdata(1,ip))
+ dpd(1,ip) = dpd(1,ip) + dm
+ ENDDO
+ ENDDO
+ ELSEIF (lda .EQ. 2) THEN
+ DO ip=1,Np
+#ifdef __SXF90
+!CDIR NODEP
+#endif
+ DO jpart=1,nvlist(ip)
+ jp = vlist(jpart,ip)
+ !---------------------------------------------------------
+ ! Calculate the square of the distance between the two
+ ! particles. It will always be .LE. (cutoff+skin)**2 by
+ ! construction of the Verlet list.
+ !---------------------------------------------------------
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !---------------------------------------------------------
+ ! Particle ip interacts with jp
+ !---------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ dm = eta*(pdata(1,jp) - pdata(1,ip))
+ dpd(1,ip) = dpd(1,ip) + dm
+ dm = eta*(pdata(2,jp) - pdata(2,ip))
+ dpd(2,ip) = dpd(2,ip) + dm
+ ENDDO
+ ENDDO
+ ELSE
+ DO ip=1,Np
+ DO jpart=1,nvlist(ip)
+ jp = vlist(jpart,ip)
+ !---------------------------------------------------------
+ ! Calculate the square of the distance between the two
+ ! particles. It will always be .LE. (cutoff+skin)**2 by
+ ! construction of the Verlet list.
+ !---------------------------------------------------------
+ dx = xp(1,ip) - xp(1,jp)
+ dy = xp(2,ip) - xp(2,jp)
+ IF (ppm_dim .GT. 2) THEN
+ dz = xp(3,ip) - xp(3,jp)
+ dij = (dx*dx) + (dy*dy) + (dz*dz)
+ ELSE
+ dz = 0.0_MK
+ dij = (dx*dx) + (dy*dy)
+ ENDIF
+ !---------------------------------------------------------
+ ! Particle ip interacts with jp
+ !---------------------------------------------------------
+#include "ppm_comp_pp_kernels.inc"
+ DO ispec=1,lda
+ dm = eta*(pdata(ispec,jp) - pdata(ispec,ip))
+ dpd(ispec,ip) = dpd(ispec,ip) + dm
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_comp_pp_verlet',t0,info)
+ RETURN
+#if __KERNEL == __INTERNAL
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_verlet_si
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_verlet_di
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_verlet_sci
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_verlet_dci
+#endif
+
+#elif __KERNEL == __USER_FUNCTION
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_verlet_su
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_verlet_du
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_verlet_scu
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_verlet_dcu
+#endif
+
+#elif __KERNEL == __LOOKUP_TABLE
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_verlet_st
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_comp_pp_verlet_dt
+#elif __KIND == __SINGLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_verlet_sct
+#elif __KIND == __DOUBLE_PRECISION_COMPLEX
+ END SUBROUTINE ppm_comp_pp_verlet_dct
+#endif
+#endif
diff --git a/src/ppm_module_comp_part.f b/src/ppm_module_comp_part.f
new file mode 100644
index 0000000000000000000000000000000000000000..cc39a2698d4294f027bc78490e74e57e2ac146c7
--- /dev/null
+++ b/src/ppm_module_comp_part.f
@@ -0,0 +1,21 @@
+ !-------------------------------------------------------------------------
+ ! Module : ppm_module_comp_part
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+ MODULE ppm_module_comp_part
+ !!! This module contains all user-callable routines to
+ !!! compute kernel-PP interactions on a set of particles.
+ !----------------------------------------------------------------------
+ ! PPM modules
+ !----------------------------------------------------------------------
+ USE ppm_module_comp_pp_verlet
+ USE ppm_module_comp_pp_cell
+ USE ppm_module_comp_pp_ring
+ USE ppm_module_comp_pp_correct
+ USE ppm_module_comp_pp_mk_table
+
+ END MODULE ppm_module_comp_part
diff --git a/src/ppm_module_comp_pp_cell.f b/src/ppm_module_comp_pp_cell.f
new file mode 100644
index 0000000000000000000000000000000000000000..32a540fa25f6000a27994b7a7f88a5833cd12c7b
--- /dev/null
+++ b/src/ppm_module_comp_pp_cell.f
@@ -0,0 +1,101 @@
+ !-------------------------------------------------------------------------
+ ! Module : ppm_module_comp_pp_cell
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Define types
+ !-------------------------------------------------------------------------
+#define __SINGLE_PRECISION 1
+#define __DOUBLE_PRECISION 2
+#define __SINGLE_PRECISION_COMPLEX 3
+#define __DOUBLE_PRECISION_COMPLEX 4
+
+#define __INTERNAL 5
+#define __USER_FUNCTION 6
+#define __LOOKUP_TABLE 7
+
+ MODULE ppm_module_comp_pp_cell
+ !!! This module provides the routines concerned
+ !!! with computing kernel interactions within a set of particles.
+ !----------------------------------------------------------------------
+ ! Define interfaces to the cell list versions
+ !----------------------------------------------------------------------
+ INTERFACE ppm_comp_pp_cell
+ MODULE PROCEDURE ppm_comp_pp_cell_si
+ MODULE PROCEDURE ppm_comp_pp_cell_di
+ MODULE PROCEDURE ppm_comp_pp_cell_sci
+ MODULE PROCEDURE ppm_comp_pp_cell_dci
+ MODULE PROCEDURE ppm_comp_pp_cell_su
+ MODULE PROCEDURE ppm_comp_pp_cell_du
+ MODULE PROCEDURE ppm_comp_pp_cell_scu
+ MODULE PROCEDURE ppm_comp_pp_cell_dcu
+ MODULE PROCEDURE ppm_comp_pp_cell_st
+ MODULE PROCEDURE ppm_comp_pp_cell_dt
+ MODULE PROCEDURE ppm_comp_pp_cell_sct
+ MODULE PROCEDURE ppm_comp_pp_cell_dct
+ END INTERFACE
+
+ !----------------------------------------------------------------------
+ ! include the source
+ !----------------------------------------------------------------------
+ CONTAINS
+
+#define __KERNEL __INTERNAL
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __USER_FUNCTION
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __LOOKUP_TABLE
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_cell.f"
+#undef __KIND
+#undef __KERNEL
+
+ END MODULE ppm_module_comp_pp_cell
diff --git a/src/ppm_module_comp_pp_correct.f b/src/ppm_module_comp_pp_correct.f
new file mode 100644
index 0000000000000000000000000000000000000000..04b2772293259a834fa3bad4e4eadeb9c29a0d51
--- /dev/null
+++ b/src/ppm_module_comp_pp_correct.f
@@ -0,0 +1,101 @@
+ !-------------------------------------------------------------------------
+ ! Module : ppm_module_comp_pp_correct
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Define types
+ !-------------------------------------------------------------------------
+#define __SINGLE_PRECISION 1
+#define __DOUBLE_PRECISION 2
+#define __SINGLE_PRECISION_COMPLEX 3
+#define __DOUBLE_PRECISION_COMPLEX 4
+
+#define __INTERNAL 5
+#define __USER_FUNCTION 6
+#define __LOOKUP_TABLE 7
+
+ MODULE ppm_module_comp_pp_correct
+ !!! This module provides the routines concerned
+ !!! with computing kernel interactions within a set of particles.
+ !----------------------------------------------------------------------
+ ! Define interfaces to the discretisation correction routine
+ !----------------------------------------------------------------------
+ INTERFACE ppm_comp_pp_correct
+ MODULE PROCEDURE ppm_comp_pp_correct_si
+ MODULE PROCEDURE ppm_comp_pp_correct_di
+ MODULE PROCEDURE ppm_comp_pp_correct_sci
+ MODULE PROCEDURE ppm_comp_pp_correct_dci
+ MODULE PROCEDURE ppm_comp_pp_correct_su
+ MODULE PROCEDURE ppm_comp_pp_correct_du
+ MODULE PROCEDURE ppm_comp_pp_correct_scu
+ MODULE PROCEDURE ppm_comp_pp_correct_dcu
+ MODULE PROCEDURE ppm_comp_pp_correct_st
+ MODULE PROCEDURE ppm_comp_pp_correct_dt
+ MODULE PROCEDURE ppm_comp_pp_correct_sct
+ MODULE PROCEDURE ppm_comp_pp_correct_dct
+ END INTERFACE
+
+ !----------------------------------------------------------------------
+ ! include the source
+ !----------------------------------------------------------------------
+ CONTAINS
+
+#define __KERNEL __INTERNAL
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __USER_FUNCTION
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __LOOKUP_TABLE
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_correct.f"
+#undef __KIND
+#undef __KERNEL
+
+ END MODULE ppm_module_comp_pp_correct
diff --git a/src/ppm_module_comp_pp_doring.f b/src/ppm_module_comp_pp_doring.f
new file mode 100644
index 0000000000000000000000000000000000000000..946f59166dfe1427567ebb3e1691a4c56b3af6e5
--- /dev/null
+++ b/src/ppm_module_comp_pp_doring.f
@@ -0,0 +1,101 @@
+ !-------------------------------------------------------------------------
+ ! Module : ppm_module_comp_pp_doring
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Define types
+ !-------------------------------------------------------------------------
+#define __SINGLE_PRECISION 1
+#define __DOUBLE_PRECISION 2
+#define __SINGLE_PRECISION_COMPLEX 3
+#define __DOUBLE_PRECISION_COMPLEX 4
+
+#define __INTERNAL 5
+#define __USER_FUNCTION 6
+#define __LOOKUP_TABLE 7
+
+ MODULE ppm_module_comp_pp_doring
+ !!! This module provides the routines concerned
+ !!! with computing kernel interactions within a set of particles.
+ !----------------------------------------------------------------------
+ ! Define interfaces to the ring interaction subroutine
+ !----------------------------------------------------------------------
+ INTERFACE ppm_comp_pp_doring
+ MODULE PROCEDURE ppm_comp_pp_doring_si
+ MODULE PROCEDURE ppm_comp_pp_doring_di
+ MODULE PROCEDURE ppm_comp_pp_doring_sci
+ MODULE PROCEDURE ppm_comp_pp_doring_dci
+ MODULE PROCEDURE ppm_comp_pp_doring_su
+ MODULE PROCEDURE ppm_comp_pp_doring_du
+ MODULE PROCEDURE ppm_comp_pp_doring_scu
+ MODULE PROCEDURE ppm_comp_pp_doring_dcu
+ MODULE PROCEDURE ppm_comp_pp_doring_st
+ MODULE PROCEDURE ppm_comp_pp_doring_dt
+ MODULE PROCEDURE ppm_comp_pp_doring_sct
+ MODULE PROCEDURE ppm_comp_pp_doring_dct
+ END INTERFACE
+
+ !----------------------------------------------------------------------
+ ! include the source
+ !----------------------------------------------------------------------
+ CONTAINS
+
+#define __KERNEL __INTERNAL
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __USER_FUNCTION
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __LOOKUP_TABLE
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_doring.f"
+#undef __KIND
+#undef __KERNEL
+
+ END MODULE ppm_module_comp_pp_doring
diff --git a/src/ppm_module_comp_pp_mk_table.f b/src/ppm_module_comp_pp_mk_table.f
new file mode 100644
index 0000000000000000000000000000000000000000..b777d4b0fad861741bd2dfcbbd8ce5d5a638d1f1
--- /dev/null
+++ b/src/ppm_module_comp_pp_mk_table.f
@@ -0,0 +1,79 @@
+ !-------------------------------------------------------------------------
+ ! Module : ppm_module_comp_pp_mk_table
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Define types
+ !-------------------------------------------------------------------------
+#define __SINGLE_PRECISION 1
+#define __DOUBLE_PRECISION 2
+#define __SINGLE_PRECISION_COMPLEX 3
+#define __DOUBLE_PRECISION_COMPLEX 4
+
+#define __INTERNAL 5
+#define __USER_FUNCTION 6
+
+ MODULE ppm_module_comp_pp_mk_table
+ !!! [[ppm_comp_pp_mk_table]]
+ !!! This module provides the routines concerned
+ !!! with computing kernel interactions within a set of particles.
+ !----------------------------------------------------------------------
+ ! Define interfaces to the lookup table creation routine
+ !----------------------------------------------------------------------
+ INTERFACE ppm_comp_pp_mk_table
+ MODULE PROCEDURE ppm_comp_pp_mk_table_si
+ MODULE PROCEDURE ppm_comp_pp_mk_table_di
+ MODULE PROCEDURE ppm_comp_pp_mk_table_sci
+ MODULE PROCEDURE ppm_comp_pp_mk_table_dci
+ MODULE PROCEDURE ppm_comp_pp_mk_table_su
+ MODULE PROCEDURE ppm_comp_pp_mk_table_du
+ MODULE PROCEDURE ppm_comp_pp_mk_table_scu
+ MODULE PROCEDURE ppm_comp_pp_mk_table_dcu
+ END INTERFACE
+
+ !----------------------------------------------------------------------
+ ! include the source
+ !----------------------------------------------------------------------
+ CONTAINS
+
+#define __KERNEL __INTERNAL
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __USER_FUNCTION
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_mk_table.f"
+#undef __KIND
+#undef __KERNEL
+
+ END MODULE ppm_module_comp_pp_mk_table
diff --git a/src/ppm_module_comp_pp_ring.f b/src/ppm_module_comp_pp_ring.f
new file mode 100644
index 0000000000000000000000000000000000000000..61e0e27b2be9871393a7223d5f459305385f8580
--- /dev/null
+++ b/src/ppm_module_comp_pp_ring.f
@@ -0,0 +1,101 @@
+ !-------------------------------------------------------------------------
+ ! Module : ppm_module_comp_pp_ring
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Define types
+ !-------------------------------------------------------------------------
+#define __SINGLE_PRECISION 1
+#define __DOUBLE_PRECISION 2
+#define __SINGLE_PRECISION_COMPLEX 3
+#define __DOUBLE_PRECISION_COMPLEX 4
+
+#define __INTERNAL 5
+#define __USER_FUNCTION 6
+#define __LOOKUP_TABLE 7
+
+ MODULE ppm_module_comp_pp_ring
+ !!! This module provides the routines concerned
+ !!! with computing kernel interactions within a set of particles.
+ !----------------------------------------------------------------------
+ ! Define interfaces to the ring versions
+ !----------------------------------------------------------------------
+ INTERFACE ppm_comp_pp_ring
+ MODULE PROCEDURE ppm_comp_pp_ring_si
+ MODULE PROCEDURE ppm_comp_pp_ring_di
+ MODULE PROCEDURE ppm_comp_pp_ring_sci
+ MODULE PROCEDURE ppm_comp_pp_ring_dci
+ MODULE PROCEDURE ppm_comp_pp_ring_su
+ MODULE PROCEDURE ppm_comp_pp_ring_du
+ MODULE PROCEDURE ppm_comp_pp_ring_scu
+ MODULE PROCEDURE ppm_comp_pp_ring_dcu
+ MODULE PROCEDURE ppm_comp_pp_ring_st
+ MODULE PROCEDURE ppm_comp_pp_ring_dt
+ MODULE PROCEDURE ppm_comp_pp_ring_sct
+ MODULE PROCEDURE ppm_comp_pp_ring_dct
+ END INTERFACE
+
+ !----------------------------------------------------------------------
+ ! include the source
+ !----------------------------------------------------------------------
+ CONTAINS
+
+#define __KERNEL __INTERNAL
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __USER_FUNCTION
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __LOOKUP_TABLE
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_ring.f"
+#undef __KIND
+#undef __KERNEL
+
+ END MODULE ppm_module_comp_pp_ring
diff --git a/src/ppm_module_comp_pp_verlet.f b/src/ppm_module_comp_pp_verlet.f
new file mode 100644
index 0000000000000000000000000000000000000000..88a57e06f2483216ce29baeb5ebbad7dc4298111
--- /dev/null
+++ b/src/ppm_module_comp_pp_verlet.f
@@ -0,0 +1,101 @@
+ !-------------------------------------------------------------------------
+ ! Module : ppm_module_comp_pp_verlet
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! ETH Zurich
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Define types
+ !-------------------------------------------------------------------------
+#define __SINGLE_PRECISION 1
+#define __DOUBLE_PRECISION 2
+#define __SINGLE_PRECISION_COMPLEX 3
+#define __DOUBLE_PRECISION_COMPLEX 4
+
+#define __INTERNAL 5
+#define __USER_FUNCTION 6
+#define __LOOKUP_TABLE 7
+
+ MODULE ppm_module_comp_pp_verlet
+ !!! This module provides the routines concerned
+ !!! with computing kernel interactions within a set of particles.
+ !----------------------------------------------------------------------
+ ! Define interfaces to the Verlet list versions
+ !----------------------------------------------------------------------
+ INTERFACE ppm_comp_pp_verlet
+ MODULE PROCEDURE ppm_comp_pp_verlet_si
+ MODULE PROCEDURE ppm_comp_pp_verlet_di
+ MODULE PROCEDURE ppm_comp_pp_verlet_sci
+ MODULE PROCEDURE ppm_comp_pp_verlet_dci
+ MODULE PROCEDURE ppm_comp_pp_verlet_su
+ MODULE PROCEDURE ppm_comp_pp_verlet_du
+ MODULE PROCEDURE ppm_comp_pp_verlet_scu
+ MODULE PROCEDURE ppm_comp_pp_verlet_dcu
+ MODULE PROCEDURE ppm_comp_pp_verlet_st
+ MODULE PROCEDURE ppm_comp_pp_verlet_dt
+ MODULE PROCEDURE ppm_comp_pp_verlet_sct
+ MODULE PROCEDURE ppm_comp_pp_verlet_dct
+ END INTERFACE
+
+ !----------------------------------------------------------------------
+ ! include the source
+ !----------------------------------------------------------------------
+ CONTAINS
+
+#define __KERNEL __INTERNAL
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __USER_FUNCTION
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+#undef __KERNEL
+
+#define __KERNEL __LOOKUP_TABLE
+#define __KIND __SINGLE_PRECISION
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __SINGLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+
+#define __KIND __DOUBLE_PRECISION_COMPLEX
+#include "ppm_comp_pp_verlet.f"
+#undef __KIND
+#undef __KERNEL
+
+ END MODULE ppm_module_comp_pp_verlet