      subroutine cdftpp(m,n,ncds,ncdspd,ilfunc,icfunc,iminj,icontr,
     &                  iicontr,inddyy,iloc,g,coord,alphaaux,coeffs,
     &                  coeffp,coeffd,tvec,dmat,zeta,axbx,ayby,azbz,px,
     &                  py,pz,base,pxax,pyay,pzaz,u,zetarec,wxcx,wycy,
     &                  wzcz,wxpx,wypy,wzpz,sss0,sss1,sss2,sss3,sss4,
     &                  pss,pss0,pss1,pss2,pss3,dss,dss0,dss1,dss2,ssp1,
     &                  psp,psp0,psp1,dsp,dsp0,dsp1,psd,psd0,dsd,dsd0,
     &                  pps0,ppp0,ppd0,zetac1,zetac2,c,c1,c2,c3,zzeta,
     &                  ppx,ppy,ppz,ppxax,ppyay,ppzaz,aaxbx,aayby,aazbz,
     &                  gather)

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)

      dimension index(9)

      dimension cdfitp(3),cdfitd(6)

      dimension ilfunc(*),icfunc(*),iminj(*),icontr(*),iicontr(*),
     &          inddyy(m,*),iloc(*)

      dimension g(*),coord(3,*),alphaaux(*),coeffs(*),coeffp(*),
     &          coeffd(*),tvec(*),dmat(*),zeta(*),axbx(*),ayby(*),
     &          azbz(*),px(*),py(*),pz(*),base(*),pxax(*),pyay(*),
     &          pzaz(*),u(*),zetarec(*),wxcx(*),wycy(*),wzcz(*),wxpx(*),
     &          wypy(*),wzpz(*),sss0(*),sss1(*),sss2(*),sss3(*),sss4(*),
     &          pss(n,*),pss0(n,*),pss1(n,*),pss2(n,*),pss3(n,*),
     &          dss(n,*),dss0(n,*),dss1(n,*),dss2(n,*),ssp1(n,*),
     &          psp(n,*),psp0(n,*),psp1(n,*),dsp(n,*),dsp0(n,*),
     &          dsp1(n,*),psd(n,*),psd0(n,*),dsd(n,*),dsd0(n,*),
     &          pps0(n,*),ppp0(n,*),ppd0(n,*),zetac1(*),zetac2(*),c(*),
     &          c1(n,*),c2(n,*),c3(n,*),zzeta(*),ppx(*),ppy(*),ppz(*),
     &          ppxax(*),ppyay(*),ppzaz(*),aaxbx(*),aayby(*),aazbz(*),
     &          gather(n,*)

      data one,two,three/1.0d0,2.0d0,3.0d0/
      data five,seven/5.0d0,7.0d0/
      data zero,half,thrhlf,pi/0.0d0,0.5d0,1.5d0,3.1415926535898d0/
      data calcut/1.0d-16/

      twoecnst=two*pi**(five/two)

      sqrt3=sqrt(three)

      dnormcd=seven/sqrt3

      do 1001 naux=1,ncds

      ilaux=ilfunc(naux)

      ccx=coord(1,icfunc(ilaux))
      ccy=coord(2,icfunc(ilaux))
      ccz=coord(3,icfunc(ilaux))

      zetaux=alphaaux(ilaux)

      cdfits=coeffs(ilaux)

      do 1002 i=1,n
 1002 zetarec(i)=one/(zeta(i)+zetaux)

      do 1003 i=1,n
      pxcx=px(i)-ccx
      pycy=py(i)-ccy
      pzcz=pz(i)-ccz
 1003 u(i)=zeta(i)*zetaux*zetarec(i)*(pxcx*pxcx+pycy*pycy+pzcz*pzcz)

      call auggv(3,n,iloc,u,psp0,ppx,ppy,ppz)

      do 1004 i=1,n
      sss0(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,1)/zeta(i)/zetaux
      sss1(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,2)/zeta(i)/zetaux
 1004 sss2(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,3)/zeta(i)/zetaux

      nn=0

      do 1005 i=1,n
      if(abs(sss0(i)).gt.calcut) then
                                       nn=nn+1

                                       iicontr(nn)=icontr(i)

                                       zzeta(nn)=zeta(i)

                                       ppx(nn)=px(i)
                                       ppy(nn)=py(i)
                                       ppz(nn)=pz(i)

                                       ppxax(nn)=pxax(i)
                                       ppyay(nn)=pyay(i)
                                       ppzaz(nn)=pzaz(i)

                                       aaxbx(nn)=axbx(icontr(i))
                                       aayby(nn)=ayby(icontr(i))
                                       aazbz(nn)=azbz(icontr(i))

                                       zetarec(nn)=zetarec(i)

                                       sss0(nn)=sss0(i)
                                       sss1(nn)=sss1(i)
                                       sss2(nn)=sss2(i)
                                 endif
 1005 continue

      do 1006 i=1,nn
 1006 zetac1(i)=zetaux*zetarec(i)

      do 1007 i=1,nn
      wxpx(i)=((zzeta(i)*ppx(i)+zetaux*ccx)*zetarec(i))-ppx(i)
      wypy(i)=((zzeta(i)*ppy(i)+zetaux*ccy)*zetarec(i))-ppy(i)
 1007 wzpz(i)=((zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i))-ppz(i)

      do 1008 i=1,nn

      pss0(i,1)=ppxax(i)*sss0(i)+wxpx(i)*sss1(i)
      pss0(i,2)=ppyay(i)*sss0(i)+wypy(i)*sss1(i)
      pss0(i,3)=ppzaz(i)*sss0(i)+wzpz(i)*sss1(i)

      pss1(i,1)=ppxax(i)*sss1(i)+wxpx(i)*sss2(i)
      pss1(i,2)=ppyay(i)*sss1(i)+wypy(i)*sss2(i)
      pss1(i,3)=ppzaz(i)*sss1(i)+wzpz(i)*sss2(i)

 1008 continue

      do 1009 i=1,nn

      c(i)=(sss0(i)-zetac1(i)*sss1(i))/(two*zzeta(i))

      dss0(i,1)=ppxax(i)*pss0(i,1)+wxpx(i)*pss1(i,1)+c(i)
      dss0(i,2)=ppyay(i)*pss0(i,1)+wypy(i)*pss1(i,1)
      dss0(i,3)=ppzaz(i)*pss0(i,1)+wzpz(i)*pss1(i,1)
      dss0(i,4)=ppyay(i)*pss0(i,2)+wypy(i)*pss1(i,2)+c(i)
      dss0(i,5)=ppzaz(i)*pss0(i,2)+wzpz(i)*pss1(i,2)
      dss0(i,6)=ppzaz(i)*pss0(i,3)+wzpz(i)*pss1(i,3)+c(i)

 1009 continue

      do 1010 i=1,nn
      pps0(i,1)=aaxbx(i)*pss0(i,1)+dss0(i,1)
      pps0(i,2)=aayby(i)*pss0(i,1)+dss0(i,2)
      pps0(i,3)=aazbz(i)*pss0(i,1)+dss0(i,3)
      pps0(i,4)=aaxbx(i)*pss0(i,2)+dss0(i,2)
      pps0(i,5)=aayby(i)*pss0(i,2)+dss0(i,4)
      pps0(i,6)=aazbz(i)*pss0(i,2)+dss0(i,5)
      pps0(i,7)=aaxbx(i)*pss0(i,3)+dss0(i,3)
      pps0(i,8)=aayby(i)*pss0(i,3)+dss0(i,5)
 1010 pps0(i,9)=aazbz(i)*pss0(i,3)+dss0(i,6)

      do 1011 i=1,nn
 1011 gather(i,1)=zero

      do 1012 j=1,9
      do 1012 i=1,nn
 1012 gather(i,1)=gather(i,1)+dmat(inddyy(iicontr(i),j))*pps0(i,j)

      nnn=0

      do 1013 i=1,nn
      if(iminj(iicontr(i)).eq.0) then
                                       nnn=nnn+1
                                       iloc(nnn)=i
                                 endif
 1013 continue

      index(1)=2
      index(2)=3
      index(3)=6

      do 1014 j=1,3
cfpp$ nodepchk
      do 1014 i=1,nnn
 1014 gather(iloc(i),1)=gather(iloc(i),1)
     &   -dmat(inddyy(iicontr(iloc(i)),index(j)))*pps0(iloc(i),index(j))

      do 1015 i=1,nn
 1015 gather(i,1)=gather(i,1)*cdfits

      do 1016 i=1,nn
 1016 tvec(naux)=tvec(naux)+gather(i,1)

 1001 continue

      do 2001 naux=1,ncdspd

      ilaux=ilfunc(ncds+naux)

      ccx=coord(1,icfunc(ilaux))
      ccy=coord(2,icfunc(ilaux))
      ccz=coord(3,icfunc(ilaux))

      zetaux=alphaaux(ilaux)

      cdfits=coeffs(ilaux)

      cdfitp(1)=coeffp(ilaux)
      cdfitp(2)=coeffp(ilaux)
      cdfitp(3)=coeffp(ilaux)

      cdfitd(1)=coeffd(ilaux)/dnormcd
      cdfitd(2)=coeffd(ilaux)
      cdfitd(3)=coeffd(ilaux)
      cdfitd(4)=coeffd(ilaux)/dnormcd
      cdfitd(5)=coeffd(ilaux)
      cdfitd(6)=coeffd(ilaux)/dnormcd

      do 2002 i=1,n
 2002 zetarec(i)=one/(zeta(i)+zetaux)

      do 2003 i=1,n
      pxcx=px(i)-ccx
      pycy=py(i)-ccy
      pzcz=pz(i)-ccz
 2003 u(i)=zeta(i)*zetaux*zetarec(i)*(pxcx*pxcx+pycy*pycy+pzcz*pzcz)

      call auggv(5,n,iloc,u,psp0,ppx,ppy,ppz)

      do 2004 i=1,n
      sss0(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,1)/zeta(i)/zetaux
      sss1(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,2)/zeta(i)/zetaux
      sss2(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,3)/zeta(i)/zetaux
      sss3(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,4)/zeta(i)/zetaux
 2004 sss4(i)=base(i)*twoecnst*sqrt(zetarec(i))*psp0(i,5)/zeta(i)/zetaux

      nn=0

      do 2005 i=1,n
      if(abs(sss0(i)).gt.calcut) then
                                       nn=nn+1

                                       iicontr(nn)=icontr(i)

                                       zzeta(nn)=zeta(i)

                                       ppx(nn)=px(i)
                                       ppy(nn)=py(i)
                                       ppz(nn)=pz(i)

                                       ppxax(nn)=pxax(i)
                                       ppyay(nn)=pyay(i)
                                       ppzaz(nn)=pzaz(i)

                                       aaxbx(nn)=axbx(icontr(i))
                                       aayby(nn)=ayby(icontr(i))
                                       aazbz(nn)=azbz(icontr(i))

                                       zetarec(nn)=zetarec(i)

                                       sss0(nn)=sss0(i)
                                       sss1(nn)=sss1(i)
                                       sss2(nn)=sss2(i)
                                       sss3(nn)=sss3(i)
                                       sss4(nn)=sss4(i)
                                 endif
 2005 continue

      do 2006 i=1,nn
      zetac1(i)=zetaux*zetarec(i)
 2006 zetac2(i)=zzeta(i)*zetarec(i)

      do 2007 i=1,nn
      wxcx(i)=((zzeta(i)*ppx(i)+zetaux*ccx)*zetarec(i))-ccx
      wycy(i)=((zzeta(i)*ppy(i)+zetaux*ccy)*zetarec(i))-ccy
 2007 wzcz(i)=((zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i))-ccz

      do 2008 i=1,nn
      wxpx(i)=((zzeta(i)*ppx(i)+zetaux*ccx)*zetarec(i))-ppx(i)
      wypy(i)=((zzeta(i)*ppy(i)+zetaux*ccy)*zetarec(i))-ppy(i)
 2008 wzpz(i)=((zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i))-ppz(i)

      do 2009 i=1,nn

      pss0(i,1)=ppxax(i)*sss0(i)+wxpx(i)*sss1(i)
      pss0(i,2)=ppyay(i)*sss0(i)+wypy(i)*sss1(i)
      pss0(i,3)=ppzaz(i)*sss0(i)+wzpz(i)*sss1(i)

      pss1(i,1)=ppxax(i)*sss1(i)+wxpx(i)*sss2(i)
      pss1(i,2)=ppyay(i)*sss1(i)+wypy(i)*sss2(i)
      pss1(i,3)=ppzaz(i)*sss1(i)+wzpz(i)*sss2(i)

      pss2(i,1)=ppxax(i)*sss2(i)+wxpx(i)*sss3(i)
      pss2(i,2)=ppyay(i)*sss2(i)+wypy(i)*sss3(i)
      pss2(i,3)=ppzaz(i)*sss2(i)+wzpz(i)*sss3(i)

      pss3(i,1)=ppxax(i)*sss3(i)+wxpx(i)*sss4(i)
      pss3(i,2)=ppyay(i)*sss3(i)+wypy(i)*sss4(i)
      pss3(i,3)=ppzaz(i)*sss3(i)+wzpz(i)*sss4(i)

 2009 continue

      do 2010 i=1,nn

      c(i)=(sss0(i)-zetac1(i)*sss1(i))/(two*zzeta(i))

      dss0(i,1)=ppxax(i)*pss0(i,1)+wxpx(i)*pss1(i,1)+c(i)
      dss0(i,2)=ppyay(i)*pss0(i,1)+wypy(i)*pss1(i,1)
      dss0(i,3)=ppzaz(i)*pss0(i,1)+wzpz(i)*pss1(i,1)
      dss0(i,4)=ppyay(i)*pss0(i,2)+wypy(i)*pss1(i,2)+c(i)
      dss0(i,5)=ppzaz(i)*pss0(i,2)+wzpz(i)*pss1(i,2)
      dss0(i,6)=ppzaz(i)*pss0(i,3)+wzpz(i)*pss1(i,3)+c(i)

      c(i)=(sss1(i)-zetac1(i)*sss2(i))/(two*zzeta(i))

      dss1(i,1)=ppxax(i)*pss1(i,1)+wxpx(i)*pss2(i,1)+c(i)
      dss1(i,2)=ppyay(i)*pss1(i,1)+wypy(i)*pss2(i,1)
      dss1(i,3)=ppzaz(i)*pss1(i,1)+wzpz(i)*pss2(i,1)
      dss1(i,4)=ppyay(i)*pss1(i,2)+wypy(i)*pss2(i,2)+c(i)
      dss1(i,5)=ppzaz(i)*pss1(i,2)+wzpz(i)*pss2(i,2)
      dss1(i,6)=ppzaz(i)*pss1(i,3)+wzpz(i)*pss2(i,3)+c(i)

      c(i)=(sss2(i)-zetac1(i)*sss3(i))/(two*zzeta(i))

      dss2(i,1)=ppxax(i)*pss2(i,1)+wxpx(i)*pss3(i,1)+c(i)
      dss2(i,2)=ppyay(i)*pss2(i,1)+wypy(i)*pss3(i,1)
      dss2(i,3)=ppzaz(i)*pss2(i,1)+wzpz(i)*pss3(i,1)
      dss2(i,4)=ppyay(i)*pss2(i,2)+wypy(i)*pss3(i,2)+c(i)
      dss2(i,5)=ppzaz(i)*pss2(i,2)+wzpz(i)*pss3(i,2)
      dss2(i,6)=ppzaz(i)*pss2(i,3)+wzpz(i)*pss3(i,3)+c(i)

 2010 continue

      do 2011 i=1,nn
      ssp1(i,1)=wxcx(i)*sss2(i)
      ssp1(i,2)=wycy(i)*sss2(i)
 2011 ssp1(i,3)=wzcz(i)*sss2(i)

      do 2012 i=1,nn

      c(i)=half*sss1(i)*zetarec(i)

      psp0(i,1)=wxcx(i)*pss1(i,1)+c(i)
      psp0(i,2)=wycy(i)*pss1(i,1)
      psp0(i,3)=wzcz(i)*pss1(i,1)
      psp0(i,4)=wxcx(i)*pss1(i,2)
      psp0(i,5)=wycy(i)*pss1(i,2)+c(i)
      psp0(i,6)=wzcz(i)*pss1(i,2)
      psp0(i,7)=wxcx(i)*pss1(i,3)
      psp0(i,8)=wycy(i)*pss1(i,3)
      psp0(i,9)=wzcz(i)*pss1(i,3)+c(i)

      c(i)=half*sss2(i)*zetarec(i)

      psp1(i,1)=wxcx(i)*pss2(i,1)+c(i)
      psp1(i,2)=wycy(i)*pss2(i,1)
      psp1(i,3)=wzcz(i)*pss2(i,1)
      psp1(i,4)=wxcx(i)*pss2(i,2)
      psp1(i,5)=wycy(i)*pss2(i,2)+c(i)
      psp1(i,6)=wzcz(i)*pss2(i,2)
      psp1(i,7)=wxcx(i)*pss2(i,3)
      psp1(i,8)=wycy(i)*pss2(i,3)
      psp1(i,9)=wzcz(i)*pss2(i,3)+c(i)

 2012 continue

      do 2013 j=1,3
      do 2013 i=1,nn
 2013 c1(i,j)=half*pss1(i,j)*zetarec(i)

      do 2014 i=1,nn
      dsp0(i,01)=wxcx(i)*dss1(i,1)+c1(i,1)*two
      dsp0(i,02)=wycy(i)*dss1(i,1)
      dsp0(i,03)=wzcz(i)*dss1(i,1)
      dsp0(i,04)=wxcx(i)*dss1(i,2)+c1(i,2)
      dsp0(i,05)=wycy(i)*dss1(i,2)+c1(i,1)
      dsp0(i,06)=wzcz(i)*dss1(i,2)
      dsp0(i,07)=wxcx(i)*dss1(i,3)+c1(i,3)
      dsp0(i,08)=wycy(i)*dss1(i,3)
      dsp0(i,09)=wzcz(i)*dss1(i,3)+c1(i,1)
      dsp0(i,10)=wxcx(i)*dss1(i,4)
      dsp0(i,11)=wycy(i)*dss1(i,4)+c1(i,2)*two
      dsp0(i,12)=wzcz(i)*dss1(i,4)
      dsp0(i,13)=wxcx(i)*dss1(i,5)
      dsp0(i,14)=wycy(i)*dss1(i,5)+c1(i,3)
      dsp0(i,15)=wzcz(i)*dss1(i,5)+c1(i,2)
      dsp0(i,16)=wxcx(i)*dss1(i,6)
      dsp0(i,17)=wycy(i)*dss1(i,6)
 2014 dsp0(i,18)=wzcz(i)*dss1(i,6)+c1(i,3)*two

      do 2015 j=1,3
      do 2015 i=1,nn
 2015 c1(i,j)=half*pss2(i,j)*zetarec(i)

      do 2016 i=1,nn
      dsp1(i,01)=wxcx(i)*dss2(i,1)+c1(i,1)*two
      dsp1(i,02)=wycy(i)*dss2(i,1)
      dsp1(i,03)=wzcz(i)*dss2(i,1)
      dsp1(i,04)=wxcx(i)*dss2(i,2)+c1(i,2)
      dsp1(i,05)=wycy(i)*dss2(i,2)+c1(i,1)
      dsp1(i,06)=wzcz(i)*dss2(i,2)
      dsp1(i,07)=wxcx(i)*dss2(i,3)+c1(i,3)
      dsp1(i,08)=wycy(i)*dss2(i,3)
      dsp1(i,09)=wzcz(i)*dss2(i,3)+c1(i,1)
      dsp1(i,10)=wxcx(i)*dss2(i,4)
      dsp1(i,11)=wycy(i)*dss2(i,4)+c1(i,2)*two
      dsp1(i,12)=wzcz(i)*dss2(i,4)
      dsp1(i,13)=wxcx(i)*dss2(i,5)
      dsp1(i,14)=wycy(i)*dss2(i,5)+c1(i,3)
      dsp1(i,15)=wzcz(i)*dss2(i,5)+c1(i,2)
      dsp1(i,16)=wxcx(i)*dss2(i,6)
      dsp1(i,17)=wycy(i)*dss2(i,6)
 2016 dsp1(i,18)=wzcz(i)*dss2(i,6)+c1(i,3)*two

      do 2017 j=1,3
      do 2017 i=1,nn
 2017 c1(i,j)=(pss0(i,j)-zetac2(i)*pss1(i,j))/(two*zetaux)

      do 2018 j=1,3
      do 2018 i=1,nn
 2018 c2(i,j)=half*ssp1(i,j)*zetarec(i)

      do 2019 i=1,nn
      psd0(i,01)=wxcx(i)*psp1(i,1)+c1(i,1)+c2(i,1)
      psd0(i,02)=wycy(i)*psp1(i,1)
      psd0(i,03)=wzcz(i)*psp1(i,1)
      psd0(i,04)=wycy(i)*psp1(i,2)+c1(i,1)
      psd0(i,05)=wzcz(i)*psp1(i,2)
      psd0(i,06)=wzcz(i)*psp1(i,3)+c1(i,1)
      psd0(i,07)=wxcx(i)*psp1(i,4)+c1(i,2)
      psd0(i,08)=wycy(i)*psp1(i,4)        +c2(i,1)
      psd0(i,09)=wzcz(i)*psp1(i,4)
      psd0(i,10)=wycy(i)*psp1(i,5)+c1(i,2)+c2(i,2)
      psd0(i,11)=wzcz(i)*psp1(i,5)
      psd0(i,12)=wzcz(i)*psp1(i,6)+c1(i,2)
      psd0(i,13)=wxcx(i)*psp1(i,7)+c1(i,3)
      psd0(i,14)=wycy(i)*psp1(i,7)
      psd0(i,15)=wzcz(i)*psp1(i,7)        +c2(i,1)
      psd0(i,16)=wycy(i)*psp1(i,8)+c1(i,3)
      psd0(i,17)=wzcz(i)*psp1(i,8)        +c2(i,2)
 2019 psd0(i,18)=wzcz(i)*psp1(i,9)+c1(i,3)+c2(i,3)

      do 2020 j=1,6
      do 2020 i=1,nn
 2020 c2(i,j)=(dss0(i,j)-zetac2(i)*dss1(i,j))/(two*zetaux)

      do 2021 j=1,9
      do 2021 i=1,nn
 2021 c3(i,j)=half*psp1(i,j)*zetarec(i)

      do 2022 i=1,nn
      dsd0(i,01)=wxcx(i)*dsp1(i,01)+c2(i,1)+c3(i,1)*two
      dsd0(i,02)=wycy(i)*dsp1(i,01)
      dsd0(i,03)=wzcz(i)*dsp1(i,01)
      dsd0(i,04)=wycy(i)*dsp1(i,02)+c2(i,1)
      dsd0(i,05)=wzcz(i)*dsp1(i,02)
      dsd0(i,06)=wzcz(i)*dsp1(i,03)+c2(i,1)
      dsd0(i,07)=wxcx(i)*dsp1(i,04)+c2(i,2)+c3(i,4)
      dsd0(i,08)=wycy(i)*dsp1(i,04)        +c3(i,1)
      dsd0(i,09)=wzcz(i)*dsp1(i,04)
      dsd0(i,10)=wycy(i)*dsp1(i,05)+c2(i,2)+c3(i,2)
      dsd0(i,11)=wzcz(i)*dsp1(i,05)
      dsd0(i,12)=wzcz(i)*dsp1(i,06)+c2(i,2)
      dsd0(i,13)=wxcx(i)*dsp1(i,07)+c2(i,3)+c3(i,7)
      dsd0(i,14)=wycy(i)*dsp1(i,07)
      dsd0(i,15)=wzcz(i)*dsp1(i,07)        +c3(i,1)
      dsd0(i,16)=wycy(i)*dsp1(i,08)+c2(i,3)
      dsd0(i,17)=wzcz(i)*dsp1(i,08)        +c3(i,2)
      dsd0(i,18)=wzcz(i)*dsp1(i,09)+c2(i,3)+c3(i,3)
      dsd0(i,19)=wxcx(i)*dsp1(i,10)+c2(i,4)
      dsd0(i,20)=wycy(i)*dsp1(i,10)        +c3(i,4)*two
      dsd0(i,21)=wzcz(i)*dsp1(i,10)
      dsd0(i,22)=wycy(i)*dsp1(i,11)+c2(i,4)+c3(i,5)*two
      dsd0(i,23)=wzcz(i)*dsp1(i,11)
      dsd0(i,24)=wzcz(i)*dsp1(i,12)+c2(i,4)
      dsd0(i,25)=wxcx(i)*dsp1(i,13)+c2(i,5)
      dsd0(i,26)=wycy(i)*dsp1(i,13)        +c3(i,7)
      dsd0(i,27)=wzcz(i)*dsp1(i,13)        +c3(i,4)
      dsd0(i,28)=wycy(i)*dsp1(i,14)+c2(i,5)+c3(i,8)
      dsd0(i,29)=wzcz(i)*dsp1(i,14)        +c3(i,5)
      dsd0(i,30)=wzcz(i)*dsp1(i,15)+c2(i,5)+c3(i,6)
      dsd0(i,31)=wxcx(i)*dsp1(i,16)+c2(i,6)
      dsd0(i,32)=wycy(i)*dsp1(i,16)
      dsd0(i,33)=wzcz(i)*dsp1(i,16)        +c3(i,7)*two
      dsd0(i,34)=wycy(i)*dsp1(i,17)+c2(i,6)
      dsd0(i,35)=wzcz(i)*dsp1(i,17)        +c3(i,8)*two
 2022 dsd0(i,36)=wzcz(i)*dsp1(i,18)+c2(i,6)+c3(i,9)*two

      do 2023 i=1,nn
      pps0(i,1)=aaxbx(i)*pss0(i,1)+dss0(i,1)
      pps0(i,2)=aayby(i)*pss0(i,1)+dss0(i,2)
      pps0(i,3)=aazbz(i)*pss0(i,1)+dss0(i,3)
      pps0(i,4)=aaxbx(i)*pss0(i,2)+dss0(i,2)
      pps0(i,5)=aayby(i)*pss0(i,2)+dss0(i,4)
      pps0(i,6)=aazbz(i)*pss0(i,2)+dss0(i,5)
      pps0(i,7)=aaxbx(i)*pss0(i,3)+dss0(i,3)
      pps0(i,8)=aayby(i)*pss0(i,3)+dss0(i,5)
 2023 pps0(i,9)=aazbz(i)*pss0(i,3)+dss0(i,6)

      do 2024 i=1,nn
      ppp0(i,01)=aaxbx(i)*psp0(i,1)+dsp0(i,01)
      ppp0(i,02)=aaxbx(i)*psp0(i,2)+dsp0(i,02)
      ppp0(i,03)=aaxbx(i)*psp0(i,3)+dsp0(i,03)
      ppp0(i,04)=aayby(i)*psp0(i,1)+dsp0(i,04)
      ppp0(i,05)=aayby(i)*psp0(i,2)+dsp0(i,05)
      ppp0(i,06)=aayby(i)*psp0(i,3)+dsp0(i,06)
      ppp0(i,07)=aazbz(i)*psp0(i,1)+dsp0(i,07)
      ppp0(i,08)=aazbz(i)*psp0(i,2)+dsp0(i,08)
      ppp0(i,09)=aazbz(i)*psp0(i,3)+dsp0(i,09)
      ppp0(i,10)=aaxbx(i)*psp0(i,4)+dsp0(i,04)
      ppp0(i,11)=aaxbx(i)*psp0(i,5)+dsp0(i,05)
      ppp0(i,12)=aaxbx(i)*psp0(i,6)+dsp0(i,06)
      ppp0(i,13)=aayby(i)*psp0(i,4)+dsp0(i,10)
      ppp0(i,14)=aayby(i)*psp0(i,5)+dsp0(i,11)
      ppp0(i,15)=aayby(i)*psp0(i,6)+dsp0(i,12)
      ppp0(i,16)=aazbz(i)*psp0(i,4)+dsp0(i,13)
      ppp0(i,17)=aazbz(i)*psp0(i,5)+dsp0(i,14)
      ppp0(i,18)=aazbz(i)*psp0(i,6)+dsp0(i,15)
      ppp0(i,19)=aaxbx(i)*psp0(i,7)+dsp0(i,07)
      ppp0(i,20)=aaxbx(i)*psp0(i,8)+dsp0(i,08)
      ppp0(i,21)=aaxbx(i)*psp0(i,9)+dsp0(i,09)
      ppp0(i,22)=aayby(i)*psp0(i,7)+dsp0(i,13)
      ppp0(i,23)=aayby(i)*psp0(i,8)+dsp0(i,14)
      ppp0(i,24)=aayby(i)*psp0(i,9)+dsp0(i,15)
      ppp0(i,25)=aazbz(i)*psp0(i,7)+dsp0(i,16)
      ppp0(i,26)=aazbz(i)*psp0(i,8)+dsp0(i,17)
 2024 ppp0(i,27)=aazbz(i)*psp0(i,9)+dsp0(i,18)

      do 2025 i=1,nn
      ppd0(i,01)=aaxbx(i)*psd0(i,01)+dsd0(i,01)
      ppd0(i,02)=aaxbx(i)*psd0(i,02)+dsd0(i,02)
      ppd0(i,03)=aaxbx(i)*psd0(i,03)+dsd0(i,03)
      ppd0(i,04)=aaxbx(i)*psd0(i,04)+dsd0(i,04)
      ppd0(i,05)=aaxbx(i)*psd0(i,05)+dsd0(i,05)
      ppd0(i,06)=aaxbx(i)*psd0(i,06)+dsd0(i,06)
      ppd0(i,07)=aayby(i)*psd0(i,01)+dsd0(i,07)
      ppd0(i,08)=aayby(i)*psd0(i,02)+dsd0(i,08)
      ppd0(i,09)=aayby(i)*psd0(i,03)+dsd0(i,09)
      ppd0(i,10)=aayby(i)*psd0(i,04)+dsd0(i,10)
      ppd0(i,11)=aayby(i)*psd0(i,05)+dsd0(i,11)
      ppd0(i,12)=aayby(i)*psd0(i,06)+dsd0(i,12)
      ppd0(i,13)=aazbz(i)*psd0(i,01)+dsd0(i,13)
      ppd0(i,14)=aazbz(i)*psd0(i,02)+dsd0(i,14)
      ppd0(i,15)=aazbz(i)*psd0(i,03)+dsd0(i,15)
      ppd0(i,16)=aazbz(i)*psd0(i,04)+dsd0(i,16)
      ppd0(i,17)=aazbz(i)*psd0(i,05)+dsd0(i,17)
      ppd0(i,18)=aazbz(i)*psd0(i,06)+dsd0(i,18)
      ppd0(i,19)=aaxbx(i)*psd0(i,07)+dsd0(i,07)
      ppd0(i,20)=aaxbx(i)*psd0(i,08)+dsd0(i,08)
      ppd0(i,21)=aaxbx(i)*psd0(i,09)+dsd0(i,09)
      ppd0(i,22)=aaxbx(i)*psd0(i,10)+dsd0(i,10)
      ppd0(i,23)=aaxbx(i)*psd0(i,11)+dsd0(i,11)
      ppd0(i,24)=aaxbx(i)*psd0(i,12)+dsd0(i,12)
      ppd0(i,25)=aayby(i)*psd0(i,07)+dsd0(i,19)
      ppd0(i,26)=aayby(i)*psd0(i,08)+dsd0(i,20)
      ppd0(i,27)=aayby(i)*psd0(i,09)+dsd0(i,21)
      ppd0(i,28)=aayby(i)*psd0(i,10)+dsd0(i,22)
      ppd0(i,29)=aayby(i)*psd0(i,11)+dsd0(i,23)
      ppd0(i,30)=aayby(i)*psd0(i,12)+dsd0(i,24)
      ppd0(i,31)=aazbz(i)*psd0(i,07)+dsd0(i,25)
      ppd0(i,32)=aazbz(i)*psd0(i,08)+dsd0(i,26)
      ppd0(i,33)=aazbz(i)*psd0(i,09)+dsd0(i,27)
      ppd0(i,34)=aazbz(i)*psd0(i,10)+dsd0(i,28)
      ppd0(i,35)=aazbz(i)*psd0(i,11)+dsd0(i,29)
      ppd0(i,36)=aazbz(i)*psd0(i,12)+dsd0(i,30)
      ppd0(i,37)=aaxbx(i)*psd0(i,13)+dsd0(i,13)
      ppd0(i,38)=aaxbx(i)*psd0(i,14)+dsd0(i,14)
      ppd0(i,39)=aaxbx(i)*psd0(i,15)+dsd0(i,15)
      ppd0(i,40)=aaxbx(i)*psd0(i,16)+dsd0(i,16)
      ppd0(i,41)=aaxbx(i)*psd0(i,17)+dsd0(i,17)
      ppd0(i,42)=aaxbx(i)*psd0(i,18)+dsd0(i,18)
      ppd0(i,43)=aayby(i)*psd0(i,13)+dsd0(i,25)
      ppd0(i,44)=aayby(i)*psd0(i,14)+dsd0(i,26)
      ppd0(i,45)=aayby(i)*psd0(i,15)+dsd0(i,27)
      ppd0(i,46)=aayby(i)*psd0(i,16)+dsd0(i,28)
      ppd0(i,47)=aayby(i)*psd0(i,17)+dsd0(i,29)
      ppd0(i,48)=aayby(i)*psd0(i,18)+dsd0(i,30)
      ppd0(i,49)=aazbz(i)*psd0(i,13)+dsd0(i,31)
      ppd0(i,50)=aazbz(i)*psd0(i,14)+dsd0(i,32)
      ppd0(i,51)=aazbz(i)*psd0(i,15)+dsd0(i,33)
      ppd0(i,52)=aazbz(i)*psd0(i,16)+dsd0(i,34)
      ppd0(i,53)=aazbz(i)*psd0(i,17)+dsd0(i,35)
 2025 ppd0(i,54)=aazbz(i)*psd0(i,18)+dsd0(i,36)

      do 2026 j=1,10
      do 2026 i=1,nn
 2026 gather(i,j)=zero

      do 2027 j=1,9
      do 2027 i=1,nn
 2027 gather(i,1)=gather(i,1)+dmat(inddyy(iicontr(i),j))*pps0(i,j)

      do 2028 k=1,3
      do 2028 j=1,9
      do 2028 i=1,nn
 2028 gather(i,k+1)=gather(i,k+1)
     &             +dmat(inddyy(iicontr(i),j))*ppp0(i,3*(j-1)+k)

      do 2029 k=1,6
      do 2029 j=1,9
      do 2029 i=1,nn
 2029 gather(i,k+4)=gather(i,k+4)
     &             +dmat(inddyy(iicontr(i),j))*ppd0(i,6*(j-1)+k)

      nnn=0

      do 2030 i=1,nn
      if(iminj(iicontr(i)).eq.0) then
                                       nnn=nnn+1
                                       iloc(nnn)=i
                                 endif
 2030 continue

      index(1)=2
      index(2)=3
      index(3)=6

      do 2031 j=1,3
cfpp$ nodepchk
      do 2031 i=1,nnn
 2031 gather(iloc(i),1)=gather(iloc(i),1)
     &                 -dmat(inddyy(iicontr(iloc(i)),index(j)))
     &                 *pps0(iloc(i),index(j))

      do 2032 k=1,3
      do 2032 j=1,3
cfpp$ nodepchk
      do 2032 i=1,nnn
 2032 gather(iloc(i),k+1)=gather(iloc(i),k+1)
     &                   -dmat(inddyy(iicontr(iloc(i)),index(j)))
     &                   *ppp0(iloc(i),3*(index(j)-1)+k)

      do 2033 k=1,6
      do 2033 j=1,3
cfpp$ nodepchk
      do 2033 i=1,nnn
 2033 gather(iloc(i),k+4)=gather(iloc(i),k+4)
     &                   -dmat(inddyy(iicontr(iloc(i)),index(j)))
     &                   *ppd0(iloc(i),6*(index(j)-1)+k)

      do 2034 i=1,nn
      gather(i,01)=gather(i,01)*cdfits
      gather(i,02)=gather(i,02)*cdfitp(1)
      gather(i,03)=gather(i,03)*cdfitp(2)
      gather(i,04)=gather(i,04)*cdfitp(3)
      gather(i,05)=gather(i,05)*cdfitd(1)
      gather(i,06)=gather(i,06)*cdfitd(2)
      gather(i,07)=gather(i,07)*cdfitd(3)
      gather(i,08)=gather(i,08)*cdfitd(4)
      gather(i,09)=gather(i,09)*cdfitd(5)
 2034 gather(i,10)=gather(i,10)*cdfitd(6)

      nc=ncds+(naux-1)*10

      do 2035 i=1,nn
      do 2035 j=1,10
 2035 tvec(nc+j)=tvec(nc+j)+gather(i,j)

 2001 continue

      return
      end
