      subroutine gradcdpp(m,n,ncds,ilfunc,icfunc,iminj,icontr,inddyy,
     &                    icentr,iicontr,g,coord,density,gradient,
     &                    cdfitc,alphaaux,coeffs,zeta,twozetaa,twozetab,
     &                    px,py,pz,base,pxax,pyay,pzaz,axbx,ayby,azbz,
     &                    zzeta,ttwozetaa,ttwozetab,zetarec,ppx,ppy,ppz,
     &                    ppxax,ppyay,ppzaz,wxpx,wypy,wzpz,zetac1,u,
     &                    sss0,sss1,sss2,sss3,pss0,pss1,pss2,dss0,dss1,
     &                    fss0,c,c1,sscond,pscond,spcond,dscond,ppcond,
     &                    fscond,dpcond,pdcond,grad)

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 ilfunc(*),icfunc(*),iminj(*),icontr(*),inddyy(m,*),
     &          icentr(m,*),iicontr(*)

      dimension g(*),coord(3,*),density(*),gradient(3,*),cdfitc(*),
     &          alphaaux(*),coeffs(*),zeta(*),twozetaa(*),twozetab(*),
     &          px(*),py(*),pz(*),base(*),pxax(*),pyay(*),pzaz(*),
     &          axbx(*),ayby(*),azbz(*),zzeta(*),ttwozetaa(*),
     &          ttwozetab(*),zetarec(*),ppx(*),ppy(*),ppz(*),ppxax(*),
     &          ppyay(*),ppzaz(*),wxpx(*),wypy(*),wzpz(*),zetac1(*),
     &          u(*),sss0(*),sss1(*),sss2(*),sss3(*),pss0(n,*),
     &          pss1(n,*),pss2(n,*),dss0(n,*),dss1(n,*),fss0(n,*),c(*),
     &          c1(n,*),sscond(*),pscond(m,*),spcond(m,*),dscond(m,*),
     &          ppcond(m,*),fscond(m,*),dpcond(m,*),pdcond(m,*),
     &          grad(m,*)

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

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

      do 1001 naux=1,ncds

      ilaux=ilfunc(naux)

      nc=icfunc(ilaux)

      ccx=coord(1,nc)
      ccy=coord(2,nc)
      ccz=coord(3,nc)

      zetaux=alphaaux(ilaux)

      cdfits=cdfitc(naux)*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)

      do 1004 i=1,n
      call augg(4,u(i),g)
      sss0(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(1)/zeta(i)/zetaux
      sss1(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(2)/zeta(i)/zetaux
      sss2(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(3)/zeta(i)/zetaux
 1004 sss3(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(4)/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)

                                       ttwozetaa(nn)=twozetaa(i)
                                       ttwozetab(nn)=twozetab(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)

                                       zetarec(nn)=zetarec(i)

                                       sss0(nn)=sss0(i)
                                       sss1(nn)=sss1(i)
                                       sss2(nn)=sss2(i)
                                       sss3(nn)=sss3(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)
 1008 pss0(i,3)=ppzaz(i)*sss0(i)+wzpz(i)*sss1(i)

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

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

      do 1011 i=1,nn
 1011 c(i)=(sss0(i)-zetac1(i)*sss1(i))/(two*zzeta(i))

      do 1012 i=1,nn
      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)
 1012 dss0(i,6)=ppzaz(i)*pss0(i,3)+wzpz(i)*pss1(i,3)+c(i)

      do 1013 i=1,nn
 1013 c(i)=(sss1(i)-zetac1(i)*sss2(i))/(two*zzeta(i))

      do 1014 i=1,nn
      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)
 1014 dss1(i,6)=ppzaz(i)*pss1(i,3)+wzpz(i)*pss2(i,3)+c(i)

      do 1015 j=1,3
      do 1015 i=1,nn
 1015 c1(i,j)=(pss0(i,j)-zetac1(i)*pss1(i,j))/(two*zzeta(i))

      do 1016 i=1,nn
      fss0(i,01)=ppxax(i)*dss0(i,1)+wxpx(i)*dss1(i,1)+c1(i,1)+c1(i,1)
      fss0(i,02)=ppyay(i)*dss0(i,1)+wypy(i)*dss1(i,1)
      fss0(i,03)=ppzaz(i)*dss0(i,1)+wzpz(i)*dss1(i,1)
      fss0(i,04)=ppyay(i)*dss0(i,2)+wypy(i)*dss1(i,2)+c1(i,1)
      fss0(i,05)=ppzaz(i)*dss0(i,2)+wzpz(i)*dss1(i,2)
      fss0(i,06)=ppzaz(i)*dss0(i,3)+wzpz(i)*dss1(i,3)+c1(i,1)
      fss0(i,07)=ppyay(i)*dss0(i,4)+wypy(i)*dss1(i,4)+c1(i,2)+c1(i,2)
      fss0(i,08)=ppzaz(i)*dss0(i,4)+wzpz(i)*dss1(i,4)
      fss0(i,09)=ppzaz(i)*dss0(i,5)+wzpz(i)*dss1(i,5)+c1(i,2)
 1016 fss0(i,10)=ppzaz(i)*dss0(i,6)+wzpz(i)*dss1(i,6)+c1(i,3)+c1(i,3)

      do 1017 i=1,m
 1017 sscond(i)=zero

      do 1018 j=1,3
      do 1018 i=1,m
 1018 pscond(i,j)=zero

      do 1019 i=1,nn
 1019 sscond(iicontr(i))=sscond(iicontr(i))-sss0(i)

      do 1020 j=1,3
      do 1020 i=1,nn
 1020 pscond(iicontr(i),j)=pscond(iicontr(i),j)-pss0(i,j)

      do 1021 i=1,m
      spcond(i,1)=axbx(i)*sscond(i)+pscond(i,1)
      spcond(i,2)=ayby(i)*sscond(i)+pscond(i,2)
 1021 spcond(i,3)=azbz(i)*sscond(i)+pscond(i,3)

      do 1022 i=1,m
      grad(i,1)=spcond(i,1)*density(inddyy(i,1))*cdfits
      grad(i,2)=spcond(i,2)*density(inddyy(i,2))*cdfits
      grad(i,3)=spcond(i,3)*density(inddyy(i,3))*cdfits
      grad(i,4)=spcond(i,1)*density(inddyy(i,4))*cdfits
      grad(i,5)=spcond(i,2)*density(inddyy(i,5))*cdfits
      grad(i,6)=spcond(i,3)*density(inddyy(i,6))*cdfits
      grad(i,7)=spcond(i,1)*density(inddyy(i,7))*cdfits
      grad(i,8)=spcond(i,2)*density(inddyy(i,8))*cdfits
 1022 grad(i,9)=spcond(i,3)*density(inddyy(i,9))*cdfits

      do 1023 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,2)=zero
                              grad(i,3)=zero
                              grad(i,6)=zero
                        endif
 1023 continue

      do 1024 k=1,3
      do 1024 j=1,3
      do 1024 i=1,m
 1024 gradient(k,nc)=gradient(k,nc)-grad(i,3*(k-1)+j)

      do 1025 k=1,3
      do 1025 j=1,3
      do 1025 i=1,m
 1025 gradient(k,icentr(i,1))=gradient(k,icentr(i,1))+grad(i,3*(k-1)+j)

      do 1026 i=1,m
      grad(i,1)=pscond(i,1)*density(inddyy(i,1))*cdfits
      grad(i,2)=pscond(i,2)*density(inddyy(i,4))*cdfits
      grad(i,3)=pscond(i,3)*density(inddyy(i,7))*cdfits
      grad(i,4)=pscond(i,1)*density(inddyy(i,2))*cdfits
      grad(i,5)=pscond(i,2)*density(inddyy(i,5))*cdfits
      grad(i,6)=pscond(i,3)*density(inddyy(i,8))*cdfits
      grad(i,7)=pscond(i,1)*density(inddyy(i,3))*cdfits
      grad(i,8)=pscond(i,2)*density(inddyy(i,6))*cdfits
 1026 grad(i,9)=pscond(i,3)*density(inddyy(i,9))*cdfits

      do 1027 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,4)=zero
                              grad(i,7)=zero
                              grad(i,8)=zero
                        endif
 1027 continue

      do 1028 k=1,3
      do 1028 j=1,3
      do 1028 i=1,m
 1028 gradient(k,nc)=gradient(k,nc)-grad(i,3*(k-1)+j)

      do 1029 k=1,3
      do 1029 j=1,3
      do 1029 i=1,m
 1029 gradient(k,icentr(i,2))=gradient(k,icentr(i,2))+grad(i,3*(k-1)+j)

      do 1030 j=1,6
      do 1030 i=1,m
 1030 dscond(i,j)=zero

      do 1031 j=1,10
      do 1031 i=1,m
 1031 fscond(i,j)=zero

      do 1032 j=1,6
      do 1032 i=1,nn
 1032 dscond(iicontr(i),j)=dscond(iicontr(i),j)+dss0(i,j)*ttwozetaa(i)

      do 1033 j=1,10
      do 1033 i=1,nn
 1033 fscond(iicontr(i),j)=fscond(iicontr(i),j)+fss0(i,j)*ttwozetaa(i)

      do 1034 i=1,m
      dpcond(i,01)=axbx(i)*dscond(i,1)+fscond(i,01)
      dpcond(i,02)=ayby(i)*dscond(i,1)+fscond(i,02)
      dpcond(i,03)=azbz(i)*dscond(i,1)+fscond(i,03)
      dpcond(i,04)=axbx(i)*dscond(i,2)+fscond(i,02)
      dpcond(i,05)=ayby(i)*dscond(i,2)+fscond(i,04)
      dpcond(i,06)=azbz(i)*dscond(i,2)+fscond(i,05)
      dpcond(i,07)=axbx(i)*dscond(i,3)+fscond(i,03)
      dpcond(i,08)=ayby(i)*dscond(i,3)+fscond(i,05)
      dpcond(i,09)=azbz(i)*dscond(i,3)+fscond(i,06)
      dpcond(i,10)=axbx(i)*dscond(i,4)+fscond(i,04)
      dpcond(i,11)=ayby(i)*dscond(i,4)+fscond(i,07)
      dpcond(i,12)=azbz(i)*dscond(i,4)+fscond(i,08)
      dpcond(i,13)=axbx(i)*dscond(i,5)+fscond(i,05)
      dpcond(i,14)=ayby(i)*dscond(i,5)+fscond(i,08)
      dpcond(i,15)=azbz(i)*dscond(i,5)+fscond(i,09)
      dpcond(i,16)=axbx(i)*dscond(i,6)+fscond(i,06)
      dpcond(i,17)=ayby(i)*dscond(i,6)+fscond(i,09)
 1034 dpcond(i,18)=azbz(i)*dscond(i,6)+fscond(i,10)

      do 1035 i=1,m
      grad(i,01)=dpcond(i,01)*density(inddyy(i,1))*cdfits
      grad(i,02)=dpcond(i,02)*density(inddyy(i,2))*cdfits
      grad(i,03)=dpcond(i,03)*density(inddyy(i,3))*cdfits
      grad(i,04)=dpcond(i,04)*density(inddyy(i,4))*cdfits
      grad(i,05)=dpcond(i,05)*density(inddyy(i,5))*cdfits
      grad(i,06)=dpcond(i,06)*density(inddyy(i,6))*cdfits
      grad(i,07)=dpcond(i,07)*density(inddyy(i,7))*cdfits
      grad(i,08)=dpcond(i,08)*density(inddyy(i,8))*cdfits
      grad(i,09)=dpcond(i,09)*density(inddyy(i,9))*cdfits
      grad(i,10)=dpcond(i,04)*density(inddyy(i,1))*cdfits
      grad(i,11)=dpcond(i,05)*density(inddyy(i,2))*cdfits
      grad(i,12)=dpcond(i,06)*density(inddyy(i,3))*cdfits
      grad(i,13)=dpcond(i,10)*density(inddyy(i,4))*cdfits
      grad(i,14)=dpcond(i,11)*density(inddyy(i,5))*cdfits
      grad(i,15)=dpcond(i,12)*density(inddyy(i,6))*cdfits
      grad(i,16)=dpcond(i,13)*density(inddyy(i,7))*cdfits
      grad(i,17)=dpcond(i,14)*density(inddyy(i,8))*cdfits
      grad(i,18)=dpcond(i,15)*density(inddyy(i,9))*cdfits
      grad(i,19)=dpcond(i,07)*density(inddyy(i,1))*cdfits
      grad(i,20)=dpcond(i,08)*density(inddyy(i,2))*cdfits
      grad(i,21)=dpcond(i,09)*density(inddyy(i,3))*cdfits
      grad(i,22)=dpcond(i,13)*density(inddyy(i,4))*cdfits
      grad(i,23)=dpcond(i,14)*density(inddyy(i,5))*cdfits
      grad(i,24)=dpcond(i,15)*density(inddyy(i,6))*cdfits
      grad(i,25)=dpcond(i,16)*density(inddyy(i,7))*cdfits
      grad(i,26)=dpcond(i,17)*density(inddyy(i,8))*cdfits
 1035 grad(i,27)=dpcond(i,18)*density(inddyy(i,9))*cdfits

      do 1036 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,02)=zero
                              grad(i,03)=zero
                              grad(i,06)=zero
                              grad(i,11)=zero
                              grad(i,12)=zero
                              grad(i,15)=zero
                              grad(i,20)=zero
                              grad(i,21)=zero
                              grad(i,24)=zero
                        endif
 1036 continue

      do 1037 k=1,3
      do 1037 j=1,9
      do 1037 i=1,m
 1037 gradient(k,nc)=gradient(k,nc)-grad(i,9*(k-1)+j)

      do 1038 k=1,3
      do 1038 j=1,9
      do 1038 i=1,m
 1038 gradient(k,icentr(i,1))=gradient(k,icentr(i,1))+grad(i,9*(k-1)+j)

      do 1039 j=1,3
      do 1039 i=1,m
 1039 pscond(i,j)=zero

      do 1040 j=1,6
      do 1040 i=1,m
 1040 dscond(i,j)=zero

      do 1041 j=1,10
      do 1041 i=1,m
 1041 fscond(i,j)=zero

      do 1042 j=1,3
      do 1042 i=1,nn
 1042 pscond(iicontr(i),j)=pscond(iicontr(i),j)+pss0(i,j)*ttwozetab(i)

      do 1043 j=1,6
      do 1043 i=1,nn
 1043 dscond(iicontr(i),j)=dscond(iicontr(i),j)+dss0(i,j)*ttwozetab(i)

      do 1044 j=1,10
      do 1044 i=1,nn
 1044 fscond(iicontr(i),j)=fscond(iicontr(i),j)+fss0(i,j)*ttwozetab(i)

      do 1045 i=1,m
      ppcond(i,1)=axbx(i)*pscond(i,1)+dscond(i,1)
      ppcond(i,2)=ayby(i)*pscond(i,1)+dscond(i,2)
      ppcond(i,3)=azbz(i)*pscond(i,1)+dscond(i,3)
      ppcond(i,4)=axbx(i)*pscond(i,2)+dscond(i,2)
      ppcond(i,5)=ayby(i)*pscond(i,2)+dscond(i,4)
      ppcond(i,6)=azbz(i)*pscond(i,2)+dscond(i,5)
      ppcond(i,7)=axbx(i)*pscond(i,3)+dscond(i,3)
      ppcond(i,8)=ayby(i)*pscond(i,3)+dscond(i,5)
 1045 ppcond(i,9)=azbz(i)*pscond(i,3)+dscond(i,6)

      do 1046 i=1,m
      dpcond(i,01)=axbx(i)*dscond(i,1)+fscond(i,01)
      dpcond(i,02)=ayby(i)*dscond(i,1)+fscond(i,02)
      dpcond(i,03)=azbz(i)*dscond(i,1)+fscond(i,03)
      dpcond(i,04)=axbx(i)*dscond(i,2)+fscond(i,02)
      dpcond(i,05)=ayby(i)*dscond(i,2)+fscond(i,04)
      dpcond(i,06)=azbz(i)*dscond(i,2)+fscond(i,05)
      dpcond(i,07)=axbx(i)*dscond(i,3)+fscond(i,03)
      dpcond(i,08)=ayby(i)*dscond(i,3)+fscond(i,05)
      dpcond(i,09)=azbz(i)*dscond(i,3)+fscond(i,06)
      dpcond(i,10)=axbx(i)*dscond(i,4)+fscond(i,04)
      dpcond(i,11)=ayby(i)*dscond(i,4)+fscond(i,07)
      dpcond(i,12)=azbz(i)*dscond(i,4)+fscond(i,08)
      dpcond(i,13)=axbx(i)*dscond(i,5)+fscond(i,05)
      dpcond(i,14)=ayby(i)*dscond(i,5)+fscond(i,08)
      dpcond(i,15)=azbz(i)*dscond(i,5)+fscond(i,09)
      dpcond(i,16)=axbx(i)*dscond(i,6)+fscond(i,06)
      dpcond(i,17)=ayby(i)*dscond(i,6)+fscond(i,09)
 1046 dpcond(i,18)=azbz(i)*dscond(i,6)+fscond(i,10)

      do 1047 i=1,m
      pdcond(i,01)=axbx(i)*ppcond(i,1)+dpcond(i,01)
      pdcond(i,02)=ayby(i)*ppcond(i,1)+dpcond(i,04)
      pdcond(i,03)=azbz(i)*ppcond(i,1)+dpcond(i,07)
      pdcond(i,04)=ayby(i)*ppcond(i,2)+dpcond(i,05)
      pdcond(i,05)=azbz(i)*ppcond(i,2)+dpcond(i,08)
      pdcond(i,06)=azbz(i)*ppcond(i,3)+dpcond(i,09)
      pdcond(i,07)=axbx(i)*ppcond(i,4)+dpcond(i,04)
      pdcond(i,08)=ayby(i)*ppcond(i,4)+dpcond(i,10)
      pdcond(i,09)=azbz(i)*ppcond(i,4)+dpcond(i,13)
      pdcond(i,10)=ayby(i)*ppcond(i,5)+dpcond(i,11)
      pdcond(i,11)=azbz(i)*ppcond(i,5)+dpcond(i,14)
      pdcond(i,12)=azbz(i)*ppcond(i,6)+dpcond(i,15)
      pdcond(i,13)=axbx(i)*ppcond(i,7)+dpcond(i,07)
      pdcond(i,14)=ayby(i)*ppcond(i,7)+dpcond(i,13)
      pdcond(i,15)=azbz(i)*ppcond(i,7)+dpcond(i,16)
      pdcond(i,16)=ayby(i)*ppcond(i,8)+dpcond(i,14)
      pdcond(i,17)=azbz(i)*ppcond(i,8)+dpcond(i,17)
 1047 pdcond(i,18)=azbz(i)*ppcond(i,9)+dpcond(i,18)

      do 1048 i=1,m
      grad(i,01)=pdcond(i,01)*density(inddyy(i,1))*cdfits
      grad(i,02)=pdcond(i,02)*density(inddyy(i,2))*cdfits
      grad(i,03)=pdcond(i,03)*density(inddyy(i,3))*cdfits
      grad(i,04)=pdcond(i,07)*density(inddyy(i,4))*cdfits
      grad(i,05)=pdcond(i,08)*density(inddyy(i,5))*cdfits
      grad(i,06)=pdcond(i,09)*density(inddyy(i,6))*cdfits
      grad(i,07)=pdcond(i,13)*density(inddyy(i,7))*cdfits
      grad(i,08)=pdcond(i,14)*density(inddyy(i,8))*cdfits
      grad(i,09)=pdcond(i,15)*density(inddyy(i,9))*cdfits
      grad(i,10)=pdcond(i,02)*density(inddyy(i,1))*cdfits
      grad(i,11)=pdcond(i,04)*density(inddyy(i,2))*cdfits
      grad(i,12)=pdcond(i,05)*density(inddyy(i,3))*cdfits
      grad(i,13)=pdcond(i,08)*density(inddyy(i,4))*cdfits
      grad(i,14)=pdcond(i,10)*density(inddyy(i,5))*cdfits
      grad(i,15)=pdcond(i,11)*density(inddyy(i,6))*cdfits
      grad(i,16)=pdcond(i,14)*density(inddyy(i,7))*cdfits
      grad(i,17)=pdcond(i,16)*density(inddyy(i,8))*cdfits
      grad(i,18)=pdcond(i,17)*density(inddyy(i,9))*cdfits
      grad(i,19)=pdcond(i,03)*density(inddyy(i,1))*cdfits
      grad(i,20)=pdcond(i,05)*density(inddyy(i,2))*cdfits
      grad(i,21)=pdcond(i,06)*density(inddyy(i,3))*cdfits
      grad(i,22)=pdcond(i,09)*density(inddyy(i,4))*cdfits
      grad(i,23)=pdcond(i,11)*density(inddyy(i,5))*cdfits
      grad(i,24)=pdcond(i,12)*density(inddyy(i,6))*cdfits
      grad(i,25)=pdcond(i,15)*density(inddyy(i,7))*cdfits
      grad(i,26)=pdcond(i,17)*density(inddyy(i,8))*cdfits
 1048 grad(i,27)=pdcond(i,18)*density(inddyy(i,9))*cdfits

      do 1049 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,02)=zero
                              grad(i,03)=zero
                              grad(i,06)=zero
                              grad(i,11)=zero
                              grad(i,12)=zero
                              grad(i,15)=zero
                              grad(i,20)=zero
                              grad(i,21)=zero
                              grad(i,24)=zero
                        endif
 1049 continue

      do 1050 k=1,3
      do 1050 j=1,9
      do 1050 i=1,m
 1050 gradient(k,nc)=gradient(k,nc)-grad(i,9*(k-1)+j)

      do 1051 k=1,3
      do 1051 j=1,9
      do 1051 i=1,m
 1051 gradient(k,icentr(i,2))=gradient(k,icentr(i,2))+grad(i,9*(k-1)+j)

 1001 continue

      return
      end
