      subroutine gradovdp(ndim,m,n,icontr,inddyy,icentr,density,overlap,
     &                    gradient,zeta,twozetaa,eta,eta2,base,etarr,
     &                    pxax,pyay,pzaz,pxbx,pyby,pzbz,twozeta,ss0,ss,
     &                    ps0,ps,ds,pp,fs,fp,c,c1,c2,ppcond,fpcond,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 icontr(*),inddyy(m,*),icentr(m,*)

      dimension density(*),overlap(ndim,*),gradient(3,*),zeta(*),
     &          twozetaa(*),eta(*),eta2(*),base(*),etarr(*),pxax(*),
     &          pyay(*),pzaz(*),pxbx(*),pyby(*),pzbz(*),twozeta(*),
     &          ss0(*),ss(*),ps0(n,*),ps(n,*),ds(n,*),pp(n,*),fs(n,*),
     &          fp(n,*),c(*),c1(n,*),c2(n,*),ppcond(m,*),fpcond(m,*),
     &          grad(m,*)

      data zero,thrhlf,two,three,pi/0.0d0,1.5d0,
     &                              2.0d0,3.0d0,3.1415926535898d0/

      do 1001 i=1,n
 1001 twozeta(i)=two*zeta(i)

      do 1002 i=1,n
 1002 ss0(i)=((pi/zeta(i))**thrhlf)*base(i)

      do 1003 i=1,n
      ps0(i,1)=pxax(i)*ss0(i)
      ps0(i,2)=pyay(i)*ss0(i)
 1003 ps0(i,3)=pzaz(i)*ss0(i)

      do 1004 i=1,n
 1004 c(i)=ss0(i)/twozeta(i)

      do 1005 i=1,n
      ds(i,1)=pxax(i)*ps0(i,1)+c(i)
      ds(i,2)=pyay(i)*ps0(i,1)
      ds(i,3)=pzaz(i)*ps0(i,1)
      ds(i,4)=pyay(i)*ps0(i,2)+c(i)
      ds(i,5)=pzaz(i)*ps0(i,2)
 1005 ds(i,6)=pzaz(i)*ps0(i,3)+c(i)

      do 1006 i=1,n
      pp(i,1)=pxbx(i)*ps0(i,1)+c(i)
      pp(i,2)=pyby(i)*ps0(i,1)
      pp(i,3)=pzbz(i)*ps0(i,1)
      pp(i,4)=pxbx(i)*ps0(i,2)
      pp(i,5)=pyby(i)*ps0(i,2)+c(i)
      pp(i,6)=pzbz(i)*ps0(i,2)
      pp(i,7)=pxbx(i)*ps0(i,3)
      pp(i,8)=pyby(i)*ps0(i,3)
 1006 pp(i,9)=pzbz(i)*ps0(i,3)+c(i)

      do 1007 j=1,3
      do 1007 i=1,n
 1007 c1(i,j)=ps0(i,j)/twozeta(i)

      do 1008 i=1,n
      fs(i,01)=pxax(i)*ds(i,1)+c1(i,1)+c1(i,1)
      fs(i,02)=pyay(i)*ds(i,1)
      fs(i,03)=pzaz(i)*ds(i,1)
      fs(i,04)=pyay(i)*ds(i,2)+c1(i,1)
      fs(i,05)=pzaz(i)*ds(i,2)
      fs(i,06)=pzaz(i)*ds(i,3)+c1(i,1)
      fs(i,07)=pyay(i)*ds(i,4)+c1(i,2)+c1(i,2)
      fs(i,08)=pzaz(i)*ds(i,4)
      fs(i,09)=pzaz(i)*ds(i,5)+c1(i,2)
 1008 fs(i,10)=pzaz(i)*ds(i,6)+c1(i,3)+c1(i,3)

      do 1009 j=1,6
      do 1009 i=1,n
 1009 c2(i,j)=ds(i,j)/twozeta(i)

      do 1010 i=1,n
      fp(i,01)=pxbx(i)*fs(i,01)+c2(i,1)*three
      fp(i,02)=pyby(i)*fs(i,01)
      fp(i,03)=pzbz(i)*fs(i,01)
      fp(i,04)=pxbx(i)*fs(i,02)+c2(i,2)*two
      fp(i,05)=pyby(i)*fs(i,02)+c2(i,1)
      fp(i,06)=pzbz(i)*fs(i,02)
      fp(i,07)=pxbx(i)*fs(i,03)+c2(i,3)*two
      fp(i,08)=pyby(i)*fs(i,03)
      fp(i,09)=pzbz(i)*fs(i,03)+c2(i,1)
      fp(i,10)=pxbx(i)*fs(i,04)+c2(i,4)
      fp(i,11)=pyby(i)*fs(i,04)+c2(i,2)*two
      fp(i,12)=pzbz(i)*fs(i,04)
      fp(i,13)=pxbx(i)*fs(i,05)+c2(i,5)
      fp(i,14)=pyby(i)*fs(i,05)+c2(i,3)
      fp(i,15)=pzbz(i)*fs(i,05)+c2(i,2)
      fp(i,16)=pxbx(i)*fs(i,06)+c2(i,6)
      fp(i,17)=pyby(i)*fs(i,06)
      fp(i,18)=pzbz(i)*fs(i,06)+c2(i,3)*two
      fp(i,19)=pxbx(i)*fs(i,07)
      fp(i,20)=pyby(i)*fs(i,07)+c2(i,4)*three
      fp(i,21)=pzbz(i)*fs(i,07)
      fp(i,22)=pxbx(i)*fs(i,08)
      fp(i,23)=pyby(i)*fs(i,08)+c2(i,5)*two
      fp(i,24)=pzbz(i)*fs(i,08)+c2(i,4)
      fp(i,25)=pxbx(i)*fs(i,09)
      fp(i,26)=pyby(i)*fs(i,09)+c2(i,6)
      fp(i,27)=pzbz(i)*fs(i,09)+c2(i,5)*two
      fp(i,28)=pxbx(i)*fs(i,10)
      fp(i,29)=pyby(i)*fs(i,10)
 1010 fp(i,30)=pzbz(i)*fs(i,10)+c2(i,6)*three

      do 1011 j=1,9
      do 1011 i=1,m
 1011 ppcond(i,j)=zero

      do 1012 j=1,30
      do 1012 i=1,m
 1012 fpcond(i,j)=zero

      do 1013 j=1,9
      do 1013 i=1,n
 1013 ppcond(icontr(i),j)=ppcond(icontr(i),j)+pp(i,j)

      do 1014 j=1,30
      do 1014 i=1,n
 1014 fpcond(icontr(i),j)=fpcond(icontr(i),j)+fp(i,j)*twozetaa(i)

      do 1015 i=1,m
      overlap(inddyy(i,01),1)=overlap(inddyy(i,01),1)-ppcond(i,1)*two
      overlap(inddyy(i,02),1)=overlap(inddyy(i,02),1)-ppcond(i,2)*two
      overlap(inddyy(i,03),1)=overlap(inddyy(i,03),1)-ppcond(i,3)*two
      overlap(inddyy(i,04),1)=overlap(inddyy(i,04),1)-ppcond(i,4)
      overlap(inddyy(i,05),1)=overlap(inddyy(i,05),1)-ppcond(i,5)
      overlap(inddyy(i,06),1)=overlap(inddyy(i,06),1)-ppcond(i,6)
      overlap(inddyy(i,07),1)=overlap(inddyy(i,07),1)-ppcond(i,7)
      overlap(inddyy(i,08),1)=overlap(inddyy(i,08),1)-ppcond(i,8)
      overlap(inddyy(i,09),1)=overlap(inddyy(i,09),1)-ppcond(i,9)
      overlap(inddyy(i,04),2)=overlap(inddyy(i,04),2)-ppcond(i,1)
      overlap(inddyy(i,05),2)=overlap(inddyy(i,05),2)-ppcond(i,2)
      overlap(inddyy(i,06),2)=overlap(inddyy(i,06),2)-ppcond(i,3)
      overlap(inddyy(i,10),2)=overlap(inddyy(i,10),2)-ppcond(i,4)*two
      overlap(inddyy(i,11),2)=overlap(inddyy(i,11),2)-ppcond(i,5)*two
      overlap(inddyy(i,12),2)=overlap(inddyy(i,12),2)-ppcond(i,6)*two
      overlap(inddyy(i,13),2)=overlap(inddyy(i,13),2)-ppcond(i,7)
      overlap(inddyy(i,14),2)=overlap(inddyy(i,14),2)-ppcond(i,8)
      overlap(inddyy(i,15),2)=overlap(inddyy(i,15),2)-ppcond(i,9)
      overlap(inddyy(i,07),3)=overlap(inddyy(i,07),3)-ppcond(i,1)
      overlap(inddyy(i,08),3)=overlap(inddyy(i,08),3)-ppcond(i,2)
      overlap(inddyy(i,09),3)=overlap(inddyy(i,09),3)-ppcond(i,3)
      overlap(inddyy(i,13),3)=overlap(inddyy(i,13),3)-ppcond(i,4)
      overlap(inddyy(i,14),3)=overlap(inddyy(i,14),3)-ppcond(i,5)
      overlap(inddyy(i,15),3)=overlap(inddyy(i,15),3)-ppcond(i,6)
      overlap(inddyy(i,16),3)=overlap(inddyy(i,16),3)-ppcond(i,7)*two
      overlap(inddyy(i,17),3)=overlap(inddyy(i,17),3)-ppcond(i,8)*two
 1015 overlap(inddyy(i,18),3)=overlap(inddyy(i,18),3)-ppcond(i,9)*two

      do 1016 i=1,m
      overlap(inddyy(i,01),1)=overlap(inddyy(i,01),1)+fpcond(i,01)
      overlap(inddyy(i,02),1)=overlap(inddyy(i,02),1)+fpcond(i,02)
      overlap(inddyy(i,03),1)=overlap(inddyy(i,03),1)+fpcond(i,03)
      overlap(inddyy(i,04),1)=overlap(inddyy(i,04),1)+fpcond(i,04)
      overlap(inddyy(i,05),1)=overlap(inddyy(i,05),1)+fpcond(i,05)
      overlap(inddyy(i,06),1)=overlap(inddyy(i,06),1)+fpcond(i,06)
      overlap(inddyy(i,07),1)=overlap(inddyy(i,07),1)+fpcond(i,07)
      overlap(inddyy(i,08),1)=overlap(inddyy(i,08),1)+fpcond(i,08)
      overlap(inddyy(i,09),1)=overlap(inddyy(i,09),1)+fpcond(i,09)
      overlap(inddyy(i,10),1)=overlap(inddyy(i,10),1)+fpcond(i,10)
      overlap(inddyy(i,11),1)=overlap(inddyy(i,11),1)+fpcond(i,11)
      overlap(inddyy(i,12),1)=overlap(inddyy(i,12),1)+fpcond(i,12)
      overlap(inddyy(i,13),1)=overlap(inddyy(i,13),1)+fpcond(i,13)
      overlap(inddyy(i,14),1)=overlap(inddyy(i,14),1)+fpcond(i,14)
      overlap(inddyy(i,15),1)=overlap(inddyy(i,15),1)+fpcond(i,15)
      overlap(inddyy(i,16),1)=overlap(inddyy(i,16),1)+fpcond(i,16)
      overlap(inddyy(i,17),1)=overlap(inddyy(i,17),1)+fpcond(i,17)
      overlap(inddyy(i,18),1)=overlap(inddyy(i,18),1)+fpcond(i,18)
      overlap(inddyy(i,01),2)=overlap(inddyy(i,01),2)+fpcond(i,04)
      overlap(inddyy(i,02),2)=overlap(inddyy(i,02),2)+fpcond(i,05)
      overlap(inddyy(i,03),2)=overlap(inddyy(i,03),2)+fpcond(i,06)
      overlap(inddyy(i,04),2)=overlap(inddyy(i,04),2)+fpcond(i,10)
      overlap(inddyy(i,05),2)=overlap(inddyy(i,05),2)+fpcond(i,11)
      overlap(inddyy(i,06),2)=overlap(inddyy(i,06),2)+fpcond(i,12)
      overlap(inddyy(i,07),2)=overlap(inddyy(i,07),2)+fpcond(i,13)
      overlap(inddyy(i,08),2)=overlap(inddyy(i,08),2)+fpcond(i,14)
      overlap(inddyy(i,09),2)=overlap(inddyy(i,09),2)+fpcond(i,15)
      overlap(inddyy(i,10),2)=overlap(inddyy(i,10),2)+fpcond(i,19)
      overlap(inddyy(i,11),2)=overlap(inddyy(i,11),2)+fpcond(i,20)
      overlap(inddyy(i,12),2)=overlap(inddyy(i,12),2)+fpcond(i,21)
      overlap(inddyy(i,13),2)=overlap(inddyy(i,13),2)+fpcond(i,22)
      overlap(inddyy(i,14),2)=overlap(inddyy(i,14),2)+fpcond(i,23)
      overlap(inddyy(i,15),2)=overlap(inddyy(i,15),2)+fpcond(i,24)
      overlap(inddyy(i,16),2)=overlap(inddyy(i,16),2)+fpcond(i,25)
      overlap(inddyy(i,17),2)=overlap(inddyy(i,17),2)+fpcond(i,26)
      overlap(inddyy(i,18),2)=overlap(inddyy(i,18),2)+fpcond(i,27)
      overlap(inddyy(i,01),3)=overlap(inddyy(i,01),3)+fpcond(i,07)
      overlap(inddyy(i,02),3)=overlap(inddyy(i,02),3)+fpcond(i,08)
      overlap(inddyy(i,03),3)=overlap(inddyy(i,03),3)+fpcond(i,09)
      overlap(inddyy(i,04),3)=overlap(inddyy(i,04),3)+fpcond(i,13)
      overlap(inddyy(i,05),3)=overlap(inddyy(i,05),3)+fpcond(i,14)
      overlap(inddyy(i,06),3)=overlap(inddyy(i,06),3)+fpcond(i,15)
      overlap(inddyy(i,07),3)=overlap(inddyy(i,07),3)+fpcond(i,16)
      overlap(inddyy(i,08),3)=overlap(inddyy(i,08),3)+fpcond(i,17)
      overlap(inddyy(i,09),3)=overlap(inddyy(i,09),3)+fpcond(i,18)
      overlap(inddyy(i,10),3)=overlap(inddyy(i,10),3)+fpcond(i,22)
      overlap(inddyy(i,11),3)=overlap(inddyy(i,11),3)+fpcond(i,23)
      overlap(inddyy(i,12),3)=overlap(inddyy(i,12),3)+fpcond(i,24)
      overlap(inddyy(i,13),3)=overlap(inddyy(i,13),3)+fpcond(i,25)
      overlap(inddyy(i,14),3)=overlap(inddyy(i,14),3)+fpcond(i,26)
      overlap(inddyy(i,15),3)=overlap(inddyy(i,15),3)+fpcond(i,27)
      overlap(inddyy(i,16),3)=overlap(inddyy(i,16),3)+fpcond(i,28)
      overlap(inddyy(i,17),3)=overlap(inddyy(i,17),3)+fpcond(i,29)
 1016 overlap(inddyy(i,18),3)=overlap(inddyy(i,18),3)+fpcond(i,30)

      do 1017 i=1,n
 1017 ss(i)=eta(i)*(three-two*etarr(i))*ss0(i)

      do 1018 i=1,n
      ps(i,1)=pxax(i)*ss(i)+eta2(i)*ps0(i,1)
      ps(i,2)=pyay(i)*ss(i)+eta2(i)*ps0(i,2)
 1018 ps(i,3)=pzaz(i)*ss(i)+eta2(i)*ps0(i,3)

      do 1019 i=1,n
 1019 c(i)=ss(i)/twozeta(i)-eta2(i)*ss0(i)/twozetaa(i)

      do 1020 i=1,n
      ds(i,1)=pxax(i)*ps(i,1)+eta2(i)*ds(i,1)+c(i)
      ds(i,2)=pyay(i)*ps(i,1)+eta2(i)*ds(i,2)
      ds(i,3)=pzaz(i)*ps(i,1)+eta2(i)*ds(i,3)
      ds(i,4)=pyay(i)*ps(i,2)+eta2(i)*ds(i,4)+c(i)
      ds(i,5)=pzaz(i)*ps(i,2)+eta2(i)*ds(i,5)
 1020 ds(i,6)=pzaz(i)*ps(i,3)+eta2(i)*ds(i,6)+c(i)

      do 1021 i=1,n
 1021 c(i)=ss(i)/twozeta(i)

      do 1022 i=1,n
      pp(i,1)=pxbx(i)*ps(i,1)+eta2(i)*pp(i,1)+c(i)
      pp(i,2)=pyby(i)*ps(i,1)+eta2(i)*pp(i,2)
      pp(i,3)=pzbz(i)*ps(i,1)+eta2(i)*pp(i,3)
      pp(i,4)=pxbx(i)*ps(i,2)+eta2(i)*pp(i,4)
      pp(i,5)=pyby(i)*ps(i,2)+eta2(i)*pp(i,5)+c(i)
      pp(i,6)=pzbz(i)*ps(i,2)+eta2(i)*pp(i,6)
      pp(i,7)=pxbx(i)*ps(i,3)+eta2(i)*pp(i,7)
      pp(i,8)=pyby(i)*ps(i,3)+eta2(i)*pp(i,8)
 1022 pp(i,9)=pzbz(i)*ps(i,3)+eta2(i)*pp(i,9)+c(i)

      do 1023 j=1,3
      do 1023 i=1,n
 1023 c1(i,j)=ps(i,j)/twozeta(i)-eta2(i)*ps0(i,j)/twozetaa(i)

      do 1024 i=1,n
      fs(i,01)=pxax(i)*ds(i,1)+eta2(i)*fs(i,01)+c1(i,1)+c1(i,1)
      fs(i,02)=pyay(i)*ds(i,1)+eta2(i)*fs(i,02)
      fs(i,03)=pzaz(i)*ds(i,1)+eta2(i)*fs(i,03)
      fs(i,04)=pyay(i)*ds(i,2)+eta2(i)*fs(i,04)+c1(i,1)
      fs(i,05)=pzaz(i)*ds(i,2)+eta2(i)*fs(i,05)
      fs(i,06)=pzaz(i)*ds(i,3)+eta2(i)*fs(i,06)+c1(i,1)
      fs(i,07)=pyay(i)*ds(i,4)+eta2(i)*fs(i,07)+c1(i,2)+c1(i,2)
      fs(i,08)=pzaz(i)*ds(i,4)+eta2(i)*fs(i,08)
      fs(i,09)=pzaz(i)*ds(i,5)+eta2(i)*fs(i,09)+c1(i,2)
 1024 fs(i,10)=pzaz(i)*ds(i,6)+eta2(i)*fs(i,10)+c1(i,3)+c1(i,3)

      do 1025 j=1,6
      do 1025 i=1,n
 1025 c2(i,j)=ds(i,j)/twozeta(i)

      do 1026 i=1,n
      fp(i,01)=pxbx(i)*fs(i,01)+eta2(i)*fp(i,01)+c2(i,1)*three
      fp(i,02)=pyby(i)*fs(i,01)+eta2(i)*fp(i,02)
      fp(i,03)=pzbz(i)*fs(i,01)+eta2(i)*fp(i,03)
      fp(i,04)=pxbx(i)*fs(i,02)+eta2(i)*fp(i,04)+c2(i,2)*two
      fp(i,05)=pyby(i)*fs(i,02)+eta2(i)*fp(i,05)+c2(i,1)
      fp(i,06)=pzbz(i)*fs(i,02)+eta2(i)*fp(i,06)
      fp(i,07)=pxbx(i)*fs(i,03)+eta2(i)*fp(i,07)+c2(i,3)*two
      fp(i,08)=pyby(i)*fs(i,03)+eta2(i)*fp(i,08)
      fp(i,09)=pzbz(i)*fs(i,03)+eta2(i)*fp(i,09)+c2(i,1)
      fp(i,10)=pxbx(i)*fs(i,04)+eta2(i)*fp(i,10)+c2(i,4)
      fp(i,11)=pyby(i)*fs(i,04)+eta2(i)*fp(i,11)+c2(i,2)*two
      fp(i,12)=pzbz(i)*fs(i,04)+eta2(i)*fp(i,12)
      fp(i,13)=pxbx(i)*fs(i,05)+eta2(i)*fp(i,13)+c2(i,5)
      fp(i,14)=pyby(i)*fs(i,05)+eta2(i)*fp(i,14)+c2(i,3)
      fp(i,15)=pzbz(i)*fs(i,05)+eta2(i)*fp(i,15)+c2(i,2)
      fp(i,16)=pxbx(i)*fs(i,06)+eta2(i)*fp(i,16)+c2(i,6)
      fp(i,17)=pyby(i)*fs(i,06)+eta2(i)*fp(i,17)
      fp(i,18)=pzbz(i)*fs(i,06)+eta2(i)*fp(i,18)+c2(i,3)*two
      fp(i,19)=pxbx(i)*fs(i,07)+eta2(i)*fp(i,19)
      fp(i,20)=pyby(i)*fs(i,07)+eta2(i)*fp(i,20)+c2(i,4)*three
      fp(i,21)=pzbz(i)*fs(i,07)+eta2(i)*fp(i,21)
      fp(i,22)=pxbx(i)*fs(i,08)+eta2(i)*fp(i,22)
      fp(i,23)=pyby(i)*fs(i,08)+eta2(i)*fp(i,23)+c2(i,5)*two
      fp(i,24)=pzbz(i)*fs(i,08)+eta2(i)*fp(i,24)+c2(i,4)
      fp(i,25)=pxbx(i)*fs(i,09)+eta2(i)*fp(i,25)
      fp(i,26)=pyby(i)*fs(i,09)+eta2(i)*fp(i,26)+c2(i,6)
      fp(i,27)=pzbz(i)*fs(i,09)+eta2(i)*fp(i,27)+c2(i,5)*two
      fp(i,28)=pxbx(i)*fs(i,10)+eta2(i)*fp(i,28)
      fp(i,29)=pyby(i)*fs(i,10)+eta2(i)*fp(i,29)
 1026 fp(i,30)=pzbz(i)*fs(i,10)+eta2(i)*fp(i,30)+c2(i,6)*three

      do 1027 j=1,9
      do 1027 i=1,m
 1027 ppcond(i,j)=zero

      do 1028 j=1,30
      do 1028 i=1,m
 1028 fpcond(i,j)=zero

      do 1029 j=1,9
      do 1029 i=1,n
 1029 ppcond(icontr(i),j)=ppcond(icontr(i),j)+pp(i,j)

      do 1030 j=1,30
      do 1030 i=1,n
 1030 fpcond(icontr(i),j)=fpcond(icontr(i),j)+fp(i,j)*twozetaa(i)

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

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

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

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

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

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

      return
      end
