      subroutine cdftmm(nreal,ninteger,nconts,ncontp,ncontd,ncds,ncdspd,
     &                  nshels,nshelp,nsheld,mtloca,iprimvc,nprimvc,
     &                  icentvc,icdlfunc,icdcfunc,isys,iwkvec,coord,
     &                  alpha,coeff,alphacd,coefscd,coefpcd,coefdcd,
     &                  tvector,density,wkvec)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      include "mpif.h"

      dimension jsend(20)

      dimension istat(mpi_status_size)

      dimension g(7)

      dimension nshels(*),nshelp(*),nsheld(*),mtloca(*),iprimvc(*),
     &          nprimvc(*),icentvc(*),icdlfunc(*),icdcfunc(*),isys(*),
     &          iwkvec(*)

      dimension coord(3,*),alpha(*),coeff(*),alphacd(*),coefscd(*),
     &          coefpcd(*),coefdcd(*),tvector(*),density(*),wkvec(*)

      data zero,half,three,cutoff/0.0d0,0.5d0,3.0d0,1.0d+35/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      nprocesses=numprocs-1

      call igftable

      sqrt3=sqrt(three)

      ndim=nconts+3*ncontp+6*ncontd
      ndim=ndim*(ndim+1)/2

      call mpi_bcast(density,ndim,mpi_real8,0,mpi_comm_world,ierr)

      ndim=ncds+10*ncdspd

      do 1001 i=1,ndim
 1001 tvector(i)=zero

      do 1002 i=1,ncontd
      do 1002 j=1,nconts

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshels(j))
      i4=i2+i3

      density(i4)=density(i4)/sqrt3
      density(i4+i1*3+3)=density(i4+i1*3+3)/sqrt3
      density(i4+i1*5+10)=density(i4+i1*5+10)/sqrt3

 1002 continue

      do 1003 i=1,ncontd
      do 1003 j=1,ncontp

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshelp(j))
      i4=i2+i3

      density(i4)=density(i4)/sqrt3
      density(i4+1)=density(i4+1)/sqrt3
      density(i4+2)=density(i4+2)/sqrt3

      density(i4+i1*3+3)=density(i4+i1*3+3)/sqrt3
      density(i4+i1*3+4)=density(i4+i1*3+4)/sqrt3
      density(i4+i1*3+5)=density(i4+i1*3+5)/sqrt3

      density(i4+i1*5+10)=density(i4+i1*5+10)/sqrt3
      density(i4+i1*5+11)=density(i4+i1*5+11)/sqrt3
      density(i4+i1*5+12)=density(i4+i1*5+12)/sqrt3

 1003 continue

      do 1004 i=1,ncontd
      do 1004 j=1,i

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nsheld(j))
      i4=i2+i3

      density(i4)=density(i4)/three
      if(i.ne.j) density(i4+1)=density(i4+1)/sqrt3
      if(i.ne.j) density(i4+2)=density(i4+2)/sqrt3
      if(i.ne.j) density(i4+3)=density(i4+3)/three
      if(i.ne.j) density(i4+4)=density(i4+4)/sqrt3
      if(i.ne.j) density(i4+5)=density(i4+5)/three

      density(i4+i1)=density(i4+i1)/sqrt3
      if(i.ne.j) density(i4+i1+3)=density(i4+i1+3)/sqrt3
      if(i.ne.j) density(i4+i1+5)=density(i4+i1+5)/sqrt3

      density(i4+i1*2+1)=density(i4+i1*2+1)/sqrt3
      if(i.ne.j) density(i4+i1*2+4)=density(i4+i1*2+4)/sqrt3
      if(i.ne.j) density(i4+i1*2+6)=density(i4+i1*2+6)/sqrt3

      density(i4+i1*3+3)=density(i4+i1*3+3)/three
      density(i4+i1*3+4)=density(i4+i1*3+4)/sqrt3
      density(i4+i1*3+5)=density(i4+i1*3+5)/sqrt3
      density(i4+i1*3+6)=density(i4+i1*3+6)/three
      if(i.ne.j) density(i4+i1*3+7)=density(i4+i1*3+7)/sqrt3
      if(i.ne.j) density(i4+i1*3+8)=density(i4+i1*3+8)/three

      density(i4+i1*4+6)=density(i4+i1*4+6)/sqrt3
      density(i4+i1*4+9)=density(i4+i1*4+9)/sqrt3
      if(i.ne.j) density(i4+i1*4+11)=density(i4+i1*4+11)/sqrt3

      density(i4+i1*5+10)=density(i4+i1*5+10)/three
      density(i4+i1*5+11)=density(i4+i1*5+11)/sqrt3
      density(i4+i1*5+12)=density(i4+i1*5+12)/sqrt3
      density(i4+i1*5+13)=density(i4+i1*5+13)/three
      density(i4+i1*5+14)=density(i4+i1*5+14)/sqrt3
      density(i4+i1*5+15)=density(i4+i1*5+15)/three

 1004 continue

      do 1005 i=1,nconts
      do 1005 j=1,i

      icenta=icentvc(nshels(i))
      icentb=icentvc(nshels(j))

      factor=zero
      if(isys(icenta).eq.1) factor=factor+half
      if(isys(icentb).eq.1) factor=factor+half

      i1=mtloca(nshels(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshels(j))
      i4=i2+i3

      density(i4)=density(i4)*factor

 1005 continue

      do 1006 i=1,ncontp
      do 1006 j=1,nconts

      icenta=icentvc(nshelp(i))
      icentb=icentvc(nshels(j))

      factor=zero
      if(isys(icenta).eq.1) factor=factor+half
      if(isys(icentb).eq.1) factor=factor+half

      i1=mtloca(nshelp(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshels(j))
      i4=i2+i3

      density(i4)=density(i4)*factor
      density(i4+i1)=density(i4+i1)*factor
      density(i4+i1+i1+1)=density(i4+i1+i1+1)*factor

 1006 continue

      do 1007 i=1,ncontp
      do 1007 j=1,i

      icenta=icentvc(nshelp(i))
      icentb=icentvc(nshelp(j))

      factor=zero
      if(isys(icenta).eq.1) factor=factor+half
      if(isys(icentb).eq.1) factor=factor+half

      i1=mtloca(nshelp(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshelp(j))
      i4=i2+i3

      density(i4)=density(i4)*factor

      density(i4+i1)=density(i4+i1)*factor
      density(i4+i1+1)=density(i4+i1+1)*factor

      density(i4+i1+i1+1)=density(i4+i1+i1+1)*factor
      density(i4+i1+i1+2)=density(i4+i1+i1+2)*factor
      density(i4+i1+i1+3)=density(i4+i1+i1+3)*factor

      if(i.ne.j) then
                       density(i4+1)=density(i4+1)*factor
                       density(i4+2)=density(i4+2)*factor

                       density(i4+i1+2)=density(i4+i1+2)*factor
                 endif

 1007 continue

      do 1008 i=1,ncontd
      do 1008 j=1,nconts

      icenta=icentvc(nsheld(i))
      icentb=icentvc(nshels(j))

      factor=zero
      if(isys(icenta).eq.1) factor=factor+half
      if(isys(icentb).eq.1) factor=factor+half

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshels(j))
      i4=i2+i3

      density(i4)=density(i4)*factor
      density(i4+i1)=density(i4+i1)*factor
      density(i4+i1*2+1)=density(i4+i1*2+1)*factor
      density(i4+i1*3+3)=density(i4+i1*3+3)*factor
      density(i4+i1*4+6)=density(i4+i1*4+6)*factor
      density(i4+i1*5+10)=density(i4+i1*5+10)*factor

 1008 continue

      do 1009 i=1,ncontd
      do 1009 j=1,ncontp

      icenta=icentvc(nsheld(i))
      icentb=icentvc(nshelp(j))

      factor=zero
      if(isys(icenta).eq.1) factor=factor+half
      if(isys(icentb).eq.1) factor=factor+half

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshelp(j))
      i4=i2+i3

      density(i4)=density(i4)*factor
      density(i4+1)=density(i4+1)*factor
      density(i4+2)=density(i4+2)*factor

      density(i4+i1)=density(i4+i1)*factor
      density(i4+i1+1)=density(i4+i1+1)*factor
      density(i4+i1+2)=density(i4+i1+2)*factor

      density(i4+i1*2+1)=density(i4+i1*2+1)*factor
      density(i4+i1*2+2)=density(i4+i1*2+2)*factor
      density(i4+i1*2+3)=density(i4+i1*2+3)*factor

      density(i4+i1*3+3)=density(i4+i1*3+3)*factor
      density(i4+i1*3+4)=density(i4+i1*3+4)*factor
      density(i4+i1*3+5)=density(i4+i1*3+5)*factor

      density(i4+i1*4+6)=density(i4+i1*4+6)*factor
      density(i4+i1*4+7)=density(i4+i1*4+7)*factor
      density(i4+i1*4+8)=density(i4+i1*4+8)*factor

      density(i4+i1*5+10)=density(i4+i1*5+10)*factor
      density(i4+i1*5+11)=density(i4+i1*5+11)*factor
      density(i4+i1*5+12)=density(i4+i1*5+12)*factor

 1009 continue

      do 1010 i=1,ncontd
      do 1010 j=1,i

      icenta=icentvc(nsheld(i))
      icentb=icentvc(nsheld(j))

      factor=zero
      if(isys(icenta).eq.1) factor=factor+half
      if(isys(icentb).eq.1) factor=factor+half

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nsheld(j))
      i4=i2+i3

      density(i4)=density(i4)*factor

      density(i4+i1)=density(i4+i1)*factor
      density(i4+i1+1)=density(i4+i1+1)*factor

      density(i4+i1*2+1)=density(i4+i1*2+1)*factor
      density(i4+i1*2+2)=density(i4+i1*2+2)*factor
      density(i4+i1*2+3)=density(i4+i1*2+3)*factor

      density(i4+i1*3+3)=density(i4+i1*3+3)*factor
      density(i4+i1*3+4)=density(i4+i1*3+4)*factor
      density(i4+i1*3+5)=density(i4+i1*3+5)*factor
      density(i4+i1*3+6)=density(i4+i1*3+6)*factor

      density(i4+i1*4+6)=density(i4+i1*4+6)*factor
      density(i4+i1*4+7)=density(i4+i1*4+7)*factor
      density(i4+i1*4+8)=density(i4+i1*4+8)*factor
      density(i4+i1*4+9)=density(i4+i1*4+9)*factor
      density(i4+i1*4+10)=density(i4+i1*4+10)*factor

      density(i4+i1*5+10)=density(i4+i1*5+10)*factor
      density(i4+i1*5+11)=density(i4+i1*5+11)*factor
      density(i4+i1*5+12)=density(i4+i1*5+12)*factor
      density(i4+i1*5+13)=density(i4+i1*5+13)*factor
      density(i4+i1*5+14)=density(i4+i1*5+14)*factor
      density(i4+i1*5+15)=density(i4+i1*5+15)*factor

      if(i.ne.j) then
                       density(i4+1)=density(i4+1)*factor
                       density(i4+2)=density(i4+2)*factor
                       density(i4+3)=density(i4+3)*factor
                       density(i4+4)=density(i4+4)*factor
                       density(i4+5)=density(i4+5)*factor

                       density(i4+i1+2)=density(i4+i1+2)*factor
                       density(i4+i1+3)=density(i4+i1+3)*factor
                       density(i4+i1+4)=density(i4+i1+4)*factor
                       density(i4+i1+5)=density(i4+i1+5)*factor

                       density(i4+i1*2+4)=density(i4+i1*2+4)*factor
                       density(i4+i1*2+5)=density(i4+i1*2+5)*factor
                       density(i4+i1*2+6)=density(i4+i1*2+6)*factor

                       density(i4+i1*3+7)=density(i4+i1*3+7)*factor
                       density(i4+i1*3+8)=density(i4+i1*3+8)*factor

                       density(i4+i1*4+11)=density(i4+i1*4+11)*factor
                 endif

 1010 continue

      ncposs=nconts*(nconts+1)/2
      ncpopp=ncontp*(ncontp+1)/2
      ncpodd=ncontd*(ncontd+1)/2

      ncpops=ncontp*nconts
      ncpods=ncontd*nconts
      ncpodp=ncontd*ncontp

      npposs=0
      nppops=0
      nppopp=0
      nppods=0
      nppodp=0
      nppodd=0

      do 1011 i=1,nconts
      do 1011 j=1,i
 1011 npposs=npposs+nprimvc(nshels(i))*nprimvc(nshels(j))

      do 1012 i=1,ncontp
      do 1012 j=1,nconts
 1012 nppops=nppops+nprimvc(nshelp(i))*nprimvc(nshels(j))

      do 1013 i=1,ncontp
      do 1013 j=1,i
 1013 nppopp=nppopp+nprimvc(nshelp(i))*nprimvc(nshelp(j))

      do 1014 i=1,ncontd
      do 1014 j=1,nconts
 1014 nppods=nppods+nprimvc(nsheld(i))*nprimvc(nshels(j))

      do 1015 i=1,ncontd
      do 1015 j=1,ncontp
 1015 nppodp=nppodp+nprimvc(nsheld(i))*nprimvc(nshelp(j))

      do 1017 i=1,ncontd
      do 1017 j=1,i
 1017 nppodd=nppodd+nprimvc(nsheld(i))*nprimvc(nsheld(j))

      ja=1
      jb=ja+npposs
      jc=jb+ncposs

      if(jc.gt.ninteger) stop 'cdftmm ss integer overflow'

      je=1
      jf=je+npposs
      jg=jf+npposs
      jh=jg+npposs
      ji=jh+npposs
      jj=ji+npposs
      jk=jj+npposs
      jl=jk+npposs
      jm=jl+npposs
      jn=jm+npposs
      jo=jn+npposs
      jp=jo+npposs
      jq=jp+npposs
      jr=jq+npposs*8

      if(jr.gt.nreal) stop 'cdftmm ss real overflow'

      call setupcdss(ncposs,n,nconts,nshels,mtloca,iprimvc,nprimvc,
     &               icentvc,iwkvec(ja),iwkvec(jb),cutoff,coord,alpha,
     &               coeff,wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh),
     &               wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),wkvec(jm),
     &               wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),density)

      jd=jc+n
      j0=jd+n

      if(j0.gt.ninteger) stop 'cdftmm ss integer overflow'

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                    endif

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n

      if(myid.eq.0) then
                          isent=0

                          do 1018 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp

                          call mpi_bsend(jsend,13,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(ja),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop

 1018                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,13,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jq+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+ncposs
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+ncposs*3
      kl=kk+nloop*3
      km=kl+nloop*3
      kn=km+ncposs*6
      ko=kn+nloop*6
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop*10

      if(kv.gt.nreal) stop 'cdftmm ss real'

      call cdftss(ncposs,nloop,ncds,ncdspd,icdlfunc,icdcfunc,iwkvec(ja),
     &            iwkvec(jc),iwkvec(jb),iwkvec(jd),g,coord,alphacd,
     &            coefscd,coefpcd,coefdcd,tvector,density,wkvec(je),
     &            wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jm),wkvec(ka),
     &            wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf),
     &            wkvec(kg),wkvec(kh),wkvec(ki),wkvec(kj),wkvec(kk),
     &            wkvec(kl),wkvec(km),wkvec(kn),wkvec(ko),wkvec(kp),
     &            wkvec(kq),wkvec(kr),wkvec(ks),wkvec(kt),wkvec(ku))

      if(ncontp.eq.0) goto 2001

      ja=1
      jb=ja+nppops
      jc=jb+ncpops*3

      if(jc.gt.ninteger) stop 'cdftmm ps integer overflow'

      je=1
      jf=je+nppops
      jg=jf+nppops
      jh=jg+nppops
      ji=jh+nppops
      jj=ji+nppops
      jk=jj+nppops
      jl=jk+nppops
      jm=jl+nppops
      jn=jm+nppops
      jo=jn+nppops
      jp=jo+nppops
      jq=jp+nppops
      jr=jq+nppops
      js=jr+nppops
      jt=js+nppops
      ju=jt+nppops*11

      if(ju.gt.nreal) stop 'cdftmm ps real overflow'

      call setupcdps(ncpops,n,ncontp,nconts,nshelp,nshels,mtloca,
     &               iprimvc,nprimvc,icentvc,iwkvec(ja),iwkvec(jb),
     &               cutoff,coord,alpha,coeff,wkvec(je),wkvec(jf),
     &               wkvec(jg),wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),
     &               wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),
     &               wkvec(jq),wkvec(jr),wkvec(js),wkvec(jt),density)

      jd=jc+n
      j0=jd+n

      if(j0.gt.ninteger) stop 'cdftmm ps integer overflow'

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                    endif

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n

      if(myid.eq.0) then
                          isent=0

                          do 1019 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js

                          call mpi_bsend(jsend,16,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(ja),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop

 1019                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,16,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jt+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+ncpops*3
      ko=kn+nloop*3
      kp=ko+nloop*3
      kq=kp+nloop*3
      kr=kq+nloop*3
      ks=kr+ncpops*9
      kt=ks+nloop*9
      ku=kt+nloop*9
      kv=ku+ncpops*18
      kw=kv+nloop*18
      kx=kw+nloop
      ky=kx+nloop
      kz=ky+nloop
      la=kz+nloop*3
      lb=la+nloop*3
      lc=lb+nloop
      ld=lc+nloop
      le=ld+nloop
      lf=le+nloop
      lg=lf+nloop
      lh=lg+nloop
      li=lh+nloop
      lj=li+nloop*10

      if(lj.gt.nreal) stop 'cdftmm ps real'

      call cdftps(ncpops,nloop,ncds,ncdspd,icdlfunc,icdcfunc,iwkvec(ja),
     &            iwkvec(jc),iwkvec(jb),iwkvec(jd),g,coord,alphacd,
     &            coefscd,coefpcd,coefdcd,tvector,density,wkvec(je),
     &            wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jm),wkvec(jq),
     &            wkvec(jr),wkvec(js),wkvec(ka),wkvec(kb),wkvec(kc),
     &            wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),
     &            wkvec(ki),wkvec(kj),wkvec(kk),wkvec(kl),wkvec(km),
     &            wkvec(kn),wkvec(ko),wkvec(kp),wkvec(kq),wkvec(kr),
     &            wkvec(ks),wkvec(kt),wkvec(ku),wkvec(kv),wkvec(kw),
     &            wkvec(kx),wkvec(ky),wkvec(kz),wkvec(la),wkvec(lb),
     &            wkvec(lc),wkvec(ld),wkvec(le),wkvec(lf),wkvec(lg),
     &            wkvec(lh),wkvec(li))

      ja=1
      jb=ja+ncpopp
      jc=jb+nppopp
      jd=jc+ncpopp*9

      if(jd.gt.ninteger) stop 'cdftmm pp integer overflow'

      je=1
      jf=je+nppopp
      jg=jf+nppopp
      jh=jg+nppopp
      ji=jh+nppopp
      jj=ji+nppopp
      jk=jj+nppopp
      jl=jk+nppopp
      jm=jl+nppopp
      jn=jm+nppopp
      jo=jn+nppopp
      jp=jo+nppopp
      jq=jp+nppopp
      jr=jq+nppopp
      js=jr+nppopp
      jt=js+nppopp
      ju=jt+nppopp
      jv=ju+nppopp
      jw=jv+nppopp
      jx=jw+ncpopp
      jy=jx+ncpopp
      jz=jy+ncpopp
      j0=jz+ncpopp*11

      if(j0.gt.nreal) stop 'cdftmm pp real overflow'

      call setupcdpp(ncpopp,n,ncontp,nshelp,mtloca,iprimvc,nprimvc,
     &               icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),cutoff,
     &               coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &               wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &               wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &               wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &               wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz),density)

      j0=jd+n
      j1=j0+n

      if(j1.gt.ninteger) stop 'cdftmm pp integer overflow'

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                    endif

      call stringcopy(ncpopp,wkvec(18*n+0*ncpopp+1),wkvec(jw))
      call stringcopy(ncpopp,wkvec(18*n+1*ncpopp+1),wkvec(jx))
      call stringcopy(ncpopp,wkvec(18*n+2*ncpopp+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpopp
      jy=jx+ncpopp
      jz=jy+ncpopp

      if(myid.eq.0) then
                          isent=0

                          do 1020 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          jsend(01)=jb
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv

                          call mpi_bsend(jsend,19,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(jb),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          jb=jb+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop

 1020                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,19,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      jb=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)

      call mpi_recv(iwkvec(jb),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop*3
      kp=ko+nloop*3
      kq=kp+nloop*3
      kr=kq+nloop*3
      ks=kr+nloop*3
      kt=ks+nloop*6
      ku=kt+nloop*6
      kv=ku+nloop*6
      kw=kv+nloop*6
      kx=kw+nloop*3
      ky=kx+nloop*9
      kz=ky+nloop*9
      la=kz+nloop*9
      lb=la+nloop*18
      lc=lb+nloop*18
      ld=lc+nloop*18
      le=ld+nloop*18
      lf=le+nloop*18
      lg=lf+nloop*36
      lh=lg+nloop*36
      li=lh+nloop*9
      lj=li+nloop*27
      lk=lj+nloop*54
      ll=lk+nloop
      lm=ll+nloop
      ln=lm+nloop
      lo=ln+nloop*3
      lp=lo+nloop*6
      lq=lp+nloop*9
      lr=lq+nloop
      ls=lr+nloop
      lt=ls+nloop
      lu=lt+nloop
      lv=lu+nloop
      lw=lv+nloop
      lx=lw+nloop
      ly=lx+nloop
      lz=ly+nloop
      ma=lz+nloop
      mb=ma+nloop*10

      if(mb.gt.nreal) stop 'cdftmm pp real'

      call cdftpp(ncpopp,nloop,ncds,ncdspd,icdlfunc,icdcfunc,iwkvec(ja),
     &            iwkvec(jb),iwkvec(jd),iwkvec(jc),iwkvec(j0),g,coord,
     &            alphacd,coefscd,coefpcd,coefdcd,tvector,density,
     &            wkvec(je),wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jn),
     &            wkvec(jo),wkvec(jp),wkvec(jm),wkvec(jq),wkvec(jr),
     &            wkvec(js),wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),
     &            wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),wkvec(ki),
     &            wkvec(kj),wkvec(kk),wkvec(kl),wkvec(km),wkvec(kn),
     &            wkvec(ko),wkvec(kp),wkvec(kq),wkvec(kr),wkvec(ks),
     &            wkvec(kt),wkvec(ku),wkvec(kv),wkvec(kw),wkvec(kx),
     &            wkvec(ky),wkvec(kz),wkvec(la),wkvec(lb),wkvec(lc),
     &            wkvec(ld),wkvec(le),wkvec(lf),wkvec(lg),wkvec(lh),
     &            wkvec(li),wkvec(lj),wkvec(lk),wkvec(ll),wkvec(lm),
     &            wkvec(ln),wkvec(lo),wkvec(lp),wkvec(lq),wkvec(lr),
     &            wkvec(ls),wkvec(lt),wkvec(lu),wkvec(lv),wkvec(lw),
     &            wkvec(lx),wkvec(ly),wkvec(lz),wkvec(ma))

      if(ncontd.eq.0) goto 2001

      ja=1
      jb=ja+nppods
      jc=jb+ncpods*6

      if(jc.gt.ninteger) stop 'cdftmm ds integer overflow'

      je=1
      jf=je+nppods
      jg=jf+nppods
      jh=jg+nppods
      ji=jh+nppods
      jj=ji+nppods
      jk=jj+nppods
      jl=jk+nppods
      jm=jl+nppods
      jn=jm+nppods
      jo=jn+nppods
      jp=jo+nppods
      jq=jp+nppods
      jr=jq+nppods
      js=jr+nppods
      jt=js+nppods
      ju=jt+nppods*11

      if(ju.gt.nreal) stop 'cdftmm ds real overflow'

      call setupcdds(ncpods,n,ncontd,nconts,nsheld,nshels,mtloca,
     &               iprimvc,nprimvc,icentvc,iwkvec(ja),iwkvec(jb),
     &               cutoff,coord,alpha,coeff,wkvec(je),wkvec(jf),
     &               wkvec(jg),wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),
     &               wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),
     &               wkvec(jq),wkvec(jr),wkvec(js),wkvec(jt),density)

      jd=jc+n
      j0=jd+n

      if(j0.gt.ninteger) stop 'cdftmm ds integer overflow'

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                    endif

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n

      if(myid.eq.0) then
                          isent=0

                          do 1021 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js

                          call mpi_bsend(jsend,16,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(ja),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop

 1021                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,16,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jt+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop*3
      kp=ko+nloop*3
      kq=kp+nloop*3
      kr=kq+nloop*3
      ks=kr+ncpods*6
      kt=ks+nloop*6
      ku=kt+nloop*6
      kv=ku+nloop*6
      kw=kv+nloop*9
      kx=kw+ncpods*18
      ky=kx+nloop*18
      kz=ky+nloop*18
      la=kz+ncpods*36
      lb=la+nloop*36
      lc=lb+nloop
      ld=lc+nloop
      le=ld+nloop
      lf=le+nloop*3
      lg=lf+nloop*6
      lh=lg+nloop*9
      li=lh+nloop
      lj=li+nloop
      lk=lj+nloop
      ll=lk+nloop
      lm=ll+nloop
      ln=lm+nloop
      lo=ln+nloop
      lp=lo+nloop*10

      if(lp.gt.nreal) stop 'cdftmm ds real'

      call cdftds(ncpods,nloop,ncds,ncdspd,icdlfunc,icdcfunc,iwkvec(ja),
     &            iwkvec(jc),iwkvec(jb),iwkvec(jd),g,coord,alphacd,
     &            coefscd,coefpcd,coefdcd,tvector,density,wkvec(je),
     &            wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jm),wkvec(jq),
     &            wkvec(jr),wkvec(js),wkvec(ka),wkvec(kb),wkvec(kc),
     &            wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),
     &            wkvec(ki),wkvec(kj),wkvec(kk),wkvec(kl),wkvec(km),
     &            wkvec(kn),wkvec(ko),wkvec(kp),wkvec(kq),wkvec(kr),
     &            wkvec(ks),wkvec(kt),wkvec(ku),wkvec(kv),wkvec(kw),
     &            wkvec(kx),wkvec(ky),wkvec(kz),wkvec(la),wkvec(lb),
     &            wkvec(lc),wkvec(ld),wkvec(le),wkvec(lf),wkvec(lg),
     &            wkvec(lh),wkvec(li),wkvec(lj),wkvec(lk),wkvec(ll),
     &            wkvec(lm),wkvec(ln),wkvec(lo))

      ja=1
      jb=ja+nppodp
      jc=jb+ncpodp*18

      if(jc.gt.ninteger) stop 'cdftmm dp integer overflow'

      je=1
      jf=je+nppodp
      jg=jf+nppodp
      jh=jg+nppodp
      ji=jh+nppodp
      jj=ji+nppodp
      jk=jj+nppodp
      jl=jk+nppodp
      jm=jl+nppodp
      jn=jm+nppodp
      jo=jn+nppodp
      jp=jo+nppodp
      jq=jp+nppodp
      jr=jq+nppodp
      js=jr+nppodp
      jt=js+nppodp
      ju=jt+nppodp
      jv=ju+nppodp
      jw=jv+nppodp
      jx=jw+ncpodp
      jy=jx+ncpodp
      jz=jy+ncpodp
      j0=jz+ncpodp*11

      if(j0.gt.nreal) stop 'cdftmm dp real overflow'

      call setupcddp(ncpodp,n,ncontd,ncontp,nsheld,nshelp,mtloca,
     &               iprimvc,nprimvc,icentvc,iwkvec(ja),iwkvec(jb),
     &               cutoff,coord,alpha,coeff,wkvec(je),wkvec(jf),
     &               wkvec(jg),wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),
     &               wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),
     &               wkvec(jq),wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),
     &               wkvec(jv),wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz),
     &               density)

      jd=jc+n
      j0=jd+n

      if(j0.gt.ninteger) stop 'cdftmm dp integer overflow'

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                    endif

      call stringcopy(ncpodp,wkvec(18*n+0*ncpodp+1),wkvec(jw))
      call stringcopy(ncpodp,wkvec(18*n+1*ncpodp+1),wkvec(jx))
      call stringcopy(ncpodp,wkvec(18*n+2*ncpodp+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpodp
      jy=jx+ncpodp
      jz=jy+ncpodp

      if(myid.eq.0) then
                          isent=0

                          do 1022 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv

                          call mpi_bsend(jsend,19,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(ja),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop

 1022                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,19,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop*3
      kq=kp+nloop*3
      kr=kq+nloop*3
      ks=kr+nloop*3
      kt=ks+nloop*3
      ku=kt+nloop*6
      kv=ku+nloop*6
      kw=kv+nloop*6
      kx=kw+nloop*6
      ky=kx+nloop*6
      kz=ky+nloop*10
      la=kz+nloop*10
      lb=la+nloop*10
      lc=lb+nloop*10
      ld=lc+nloop*9
      le=ld+nloop*18
      lf=le+nloop*18
      lg=lf+nloop*18
      lh=lg+nloop*30
      li=lh+nloop*30
      lj=li+nloop*30
      lk=lj+nloop*36
      ll=lk+nloop*36
      lm=ll+nloop*60
      ln=lm+nloop*60
      lo=ln+nloop*18
      lp=lo+nloop*54
      lq=lp+nloop*108
      lr=lq+nloop
      ls=lr+nloop
      lt=ls+nloop
      lu=lt+nloop*3
      lv=lu+nloop*6
      lw=lv+nloop*9
      lx=lw+nloop*10
      ly=lx+nloop*18
      lz=ly+nloop
      ma=lz+nloop
      mb=ma+nloop
      mc=mb+nloop
      md=mc+nloop
      me=md+nloop
      mf=me+nloop
      mg=mf+nloop
      mh=mg+nloop
      mi=mh+nloop
      mj=mi+nloop*10

      if(mj.gt.nreal) stop 'cdftmm dp real'

      call cdftdp(ncpodp,nloop,ncds,ncdspd,icdlfunc,icdcfunc,iwkvec(ja),
     &            iwkvec(jc),iwkvec(jb),iwkvec(jd),g,coord,alphacd,
     &            coefscd,coefpcd,coefdcd,tvector,density,wkvec(je),
     &            wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jn),wkvec(jo),
     &            wkvec(jp),wkvec(jm),wkvec(jq),wkvec(jr),wkvec(js),
     &            wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),
     &            wkvec(kf),wkvec(kg),wkvec(kh),wkvec(ki),wkvec(kj),
     &            wkvec(kk),wkvec(kl),wkvec(km),wkvec(kn),wkvec(ko),
     &            wkvec(kp),wkvec(kq),wkvec(kr),wkvec(ks),wkvec(kt),
     &            wkvec(ku),wkvec(kv),wkvec(kw),wkvec(kx),wkvec(ky),
     &            wkvec(kz),wkvec(la),wkvec(lb),wkvec(lc),wkvec(ld),
     &            wkvec(le),wkvec(lf),wkvec(lg),wkvec(lh),wkvec(li),
     &            wkvec(lj),wkvec(lk),wkvec(ll),wkvec(lm),wkvec(ln),
     &            wkvec(lo),wkvec(lp),wkvec(lq),wkvec(lr),wkvec(ls),
     &            wkvec(lt),wkvec(lu),wkvec(lv),wkvec(lw),wkvec(lx),
     &            wkvec(ly),wkvec(lz),wkvec(ma),wkvec(mb),wkvec(mc),
     &            wkvec(md),wkvec(me),wkvec(mf),wkvec(mg),wkvec(mh),
     &            wkvec(mi))

      ja=1
      jb=ja+ncpodd
      jc=jb+nppodd
      jd=jc+ncpodd*36

      if(jd.gt.ninteger) stop 'cdftmm dd integer overflow'

      je=1
      jf=je+nppodd
      jg=jf+nppodd
      jh=jg+nppodd
      ji=jh+nppodd
      jj=ji+nppodd
      jk=jj+nppodd
      jl=jk+nppodd
      jm=jl+nppodd
      jn=jm+nppodd
      jo=jn+nppodd
      jp=jo+nppodd
      jq=jp+nppodd
      jr=jq+nppodd
      js=jr+nppodd
      jt=js+nppodd
      ju=jt+nppodd
      jv=ju+nppodd
      jw=jv+nppodd
      jx=jw+ncpodd
      jy=jx+ncpodd
      jz=jy+ncpodd
      j0=jz+ncpodd*11

      if(j0.gt.nreal) stop 'cdftmm dd real overflow'

      call setupcddd(ncpodd,n,ncontd,nsheld,mtloca,iprimvc,nprimvc,
     &               icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),cutoff,
     &               coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &               wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &               wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &               wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &               wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz),density)

      j0=jd+n
      j1=j0+n

      if(j1.gt.ninteger) stop 'cdftmm dd integer overflow'

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                    endif

      call stringcopy(ncpodd,wkvec(18*n+0*ncpodd+1),wkvec(jw))
      call stringcopy(ncpodd,wkvec(18*n+1*ncpodd+1),wkvec(jx))
      call stringcopy(ncpodd,wkvec(18*n+2*ncpodd+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpodd
      jy=jx+ncpodd
      jz=jy+ncpodd

      if(myid.eq.0) then
                          isent=0

                          do 1023 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          jsend(01)=jb
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv

                          call mpi_bsend(jsend,19,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(jb),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          jb=jb+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop

 1023                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,19,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      jb=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)

      call mpi_recv(iwkvec(jb),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop*3
      kr=kq+nloop*3
      ks=kr+nloop*3
      kt=ks+nloop*3
      ku=kt+nloop*3
      kv=ku+nloop*3
      kw=kv+nloop*6
      kx=kw+nloop*6
      ky=kx+nloop*6
      kz=ky+nloop*6
      la=kz+nloop*6
      lb=la+nloop*6
      lc=lb+nloop*10
      ld=lc+nloop*10
      le=ld+nloop*10
      lf=le+nloop*10
      lg=lf+nloop*10
      lh=lg+nloop*15
      li=lh+nloop*15
      lj=li+nloop*15
      lk=lj+nloop*15
      ll=lk+nloop*9
      lm=ll+nloop*9
      ln=lm+nloop*18
      lo=ln+nloop*18
      lp=lo+nloop*18
      lq=lp+nloop*30
      lr=lq+nloop*30
      ls=lr+nloop*30
      lt=ls+nloop*45
      lu=lt+nloop*45
      lv=lu+nloop*45
      lw=lv+nloop*36
      lx=lw+nloop*36
      ly=lx+nloop*60
      lz=ly+nloop*60
      ma=lz+nloop*90
      mb=ma+nloop*90
      mc=mb+nloop*18
      md=mc+nloop*30
      me=md+nloop*36
      mf=me+nloop*54
      mg=mf+nloop*90
      mh=mg+nloop*108
      mi=mh+nloop*108
      mj=mi+nloop*180
      mk=mj+nloop*216
      ml=mk+nloop
      mm=ml+nloop
      mn=mm+nloop
      mo=mn+nloop*3
      mp=mo+nloop*6
      mq=mp+nloop*10
      mr=mq+nloop*15
      ms=mr+nloop*30
      mt=ms+nloop
      mu=mt+nloop
      mv=mu+nloop
      mw=mv+nloop
      mx=mw+nloop
      my=mx+nloop
      mz=my+nloop
      na=mz+nloop
      nb=na+nloop
      nc=nb+nloop
      nd=nc+nloop*10
                    
      if(mz.gt.nreal) stop 'cdftmm dd real'

      call cdftdd(ncpodd,nloop,ncds,ncdspd,icdlfunc,icdcfunc,iwkvec(ja),
     &            iwkvec(jb),iwkvec(jd),iwkvec(jc),iwkvec(j0),g,coord,
     &            alphacd,coefscd,coefpcd,coefdcd,tvector,density,
     &            wkvec(je),wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jn),
     &            wkvec(jo),wkvec(jp),wkvec(jm),wkvec(jq),wkvec(jr),
     &            wkvec(js),wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),
     &            wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),wkvec(ki),
     &            wkvec(kj),wkvec(kk),wkvec(kl),wkvec(km),wkvec(kn),
     &            wkvec(ko),wkvec(kp),wkvec(kq),wkvec(kr),wkvec(ks),
     &            wkvec(kt),wkvec(ku),wkvec(kv),wkvec(kw),wkvec(kx),
     &            wkvec(ky),wkvec(kz),wkvec(la),wkvec(lb),wkvec(lc),
     &            wkvec(ld),wkvec(le),wkvec(lf),wkvec(lg),wkvec(lh),
     &            wkvec(li),wkvec(lj),wkvec(lk),wkvec(ll),wkvec(lm),
     &            wkvec(ln),wkvec(lo),wkvec(lp),wkvec(lq),wkvec(lr),
     &            wkvec(ls),wkvec(lt),wkvec(lu),wkvec(lv),wkvec(lw),
     &            wkvec(lx),wkvec(ly),wkvec(lz),wkvec(ma),wkvec(mb),
     &            wkvec(mc),wkvec(md),wkvec(me),wkvec(mf),wkvec(mg),
     &            wkvec(mh),wkvec(mi),wkvec(mj),wkvec(mk),wkvec(ml),
     &            wkvec(mm),wkvec(mn),wkvec(mo),wkvec(mp),wkvec(mq),
     &            wkvec(mr),wkvec(ms),wkvec(mt),wkvec(mu),wkvec(mv),
     &            wkvec(mw),wkvec(mx),wkvec(my),wkvec(mz),wkvec(na),
     &            wkvec(nb),wkvec(nc))

 2001 continue

      call mpi_reduce(tvector,wkvec,ndim,mpi_real8,
     &                mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) call stringcopy(ndim,tvector,wkvec)

      return
      end
