      subroutine gradxxdd4(m,n,nn,nc,iminj,inddyy,icentr,iicontr,cdfits,
     &                     cdfitp,cdfitd,gradient,density,ttwozetab,
     &                     dss0,dsp0,dsd0,fss0,fsp0,fsd0,gss0,gsp0,gsd0,
     &                     hss0,hsp0,hsd0,axbx,ayby,azbz,dsscond,
     &                     dspcond,dsdcond,fsscond,fspcond,fsdcond,
     &                     gsscond,gspcond,gsdcond,hsscond,hspcond,
     &                     hsdcond,dpscond,dppcond,dpdcond,fpscond,
     &                     fppcond,fpdcond,gpscond,gppcond,gpdcond,
     &                     ddscond,ddpcond,dddcond,fdscond,fdpcond,
     &                     fddcond,dfscond,dfpcond,dfdcond,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 iminj(*),inddyy(m,*),icentr(m,*),iicontr(*)

      dimension cdfitp(*),cdfitd(*),gradient(3,*),density(*),
     &          ttwozetab(*),dss0(n,*),dsp0(n,*),dsd0(n,*),fss0(n,*),
     &          fsp0(n,*),fsd0(n,*),gss0(n,*),gsp0(n,*),gsd0(n,*),
     &          hss0(n,*),hsp0(n,*),hsd0(n,*),axbx(*),ayby(*),azbz(*),
     &          dsscond(m,*),dspcond(m,*),dsdcond(m,*),fsscond(m,*),
     &          fspcond(m,*),fsdcond(m,*),gsscond(m,*),gspcond(m,*),
     &          gsdcond(m,*),hsscond(m,*),hspcond(m,*),hsdcond(m,*),
     &          dpscond(m,*),dppcond(m,*),dpdcond(m,*),fpscond(m,*),
     &          fppcond(m,*),fpdcond(m,*),gpscond(m,*),gppcond(m,*),
     &          gpdcond(m,*),ddscond(m,*),ddpcond(m,*),dddcond(m,*),
     &          fdscond(m,*),fdpcond(m,*),fddcond(m,*),dfscond(m,*),
     &          dfpcond(m,*),dfdcond(m,*),grad(m,*)

      data zero/0.0d0/

      do 1001 j=1,6
      do 1001 i=1,m
 1001 dsscond(i,j)=zero

      do 1002 j=1,10
      do 1002 i=1,m
 1002 fsscond(i,j)=zero

      do 1003 j=1,15
      do 1003 i=1,m
 1003 gsscond(i,j)=zero

      do 1004 j=1,21
      do 1004 i=1,m
 1004 hsscond(i,j)=zero

      do 1005 j=1,6
      do 1005 i=1,nn
 1005 dsscond(iicontr(i),j)=dsscond(iicontr(i),j)
     &                     +dss0(i,j)*ttwozetab(i)*cdfits

      do 1006 k=1,3
      do 1006 j=1,6
      do 1006 i=1,nn
 1006 dsscond(iicontr(i),j)=dsscond(iicontr(i),j)
     &                     +dsp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1007 k=1,6
      do 1007 j=1,6
      do 1007 i=1,nn
 1007 dsscond(iicontr(i),j)=dsscond(iicontr(i),j)
     &                     +dsd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

      do 1008 j=1,10
      do 1008 i=1,nn
 1008 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fss0(i,j)*ttwozetab(i)*cdfits

      do 1009 k=1,3
      do 1009 j=1,10
      do 1009 i=1,nn
 1009 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fsp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1010 k=1,6
      do 1010 j=1,10
      do 1010 i=1,nn
 1010 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fsd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

      do 1011 j=1,15
      do 1011 i=1,nn
 1011 gsscond(iicontr(i),j)=gsscond(iicontr(i),j)
     &                     +gss0(i,j)*ttwozetab(i)*cdfits

      do 1012 k=1,3
      do 1012 j=1,15
      do 1012 i=1,nn
 1012 gsscond(iicontr(i),j)=gsscond(iicontr(i),j)
     &                     +gsp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1013 k=1,6
      do 1013 j=1,15
      do 1013 i=1,nn
 1013 gsscond(iicontr(i),j)=gsscond(iicontr(i),j)
     &                     +gsd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

      do 1014 j=1,21
      do 1014 i=1,nn
 1014 hsscond(iicontr(i),j)=hsscond(iicontr(i),j)
     &                     +hss0(i,j)*ttwozetab(i)*cdfits

      do 1015 k=1,3
      do 1015 j=1,21
      do 1015 i=1,nn
 1015 hsscond(iicontr(i),j)=hsscond(iicontr(i),j)
     &                     +hsp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1016 k=1,6
      do 1016 j=1,21
      do 1016 i=1,nn
 1016 hsscond(iicontr(i),j)=hsscond(iicontr(i),j)
     &                     +hsd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

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

      do 1018 i=1,m
      fpscond(i,01)=axbx(i)*fsscond(i,01)+gsscond(i,01)
      fpscond(i,02)=ayby(i)*fsscond(i,01)+gsscond(i,02)
      fpscond(i,03)=azbz(i)*fsscond(i,01)+gsscond(i,03)
      fpscond(i,04)=axbx(i)*fsscond(i,02)+gsscond(i,02)
      fpscond(i,05)=ayby(i)*fsscond(i,02)+gsscond(i,04)
      fpscond(i,06)=azbz(i)*fsscond(i,02)+gsscond(i,05)
      fpscond(i,07)=axbx(i)*fsscond(i,03)+gsscond(i,03)
      fpscond(i,08)=ayby(i)*fsscond(i,03)+gsscond(i,05)
      fpscond(i,09)=azbz(i)*fsscond(i,03)+gsscond(i,06)
      fpscond(i,10)=axbx(i)*fsscond(i,04)+gsscond(i,04)
      fpscond(i,11)=ayby(i)*fsscond(i,04)+gsscond(i,07)
      fpscond(i,12)=azbz(i)*fsscond(i,04)+gsscond(i,08)
      fpscond(i,13)=axbx(i)*fsscond(i,05)+gsscond(i,05)
      fpscond(i,14)=ayby(i)*fsscond(i,05)+gsscond(i,08)
      fpscond(i,15)=azbz(i)*fsscond(i,05)+gsscond(i,09)
      fpscond(i,16)=axbx(i)*fsscond(i,06)+gsscond(i,06)
      fpscond(i,17)=ayby(i)*fsscond(i,06)+gsscond(i,09)
      fpscond(i,18)=azbz(i)*fsscond(i,06)+gsscond(i,10)
      fpscond(i,19)=axbx(i)*fsscond(i,07)+gsscond(i,07)
      fpscond(i,20)=ayby(i)*fsscond(i,07)+gsscond(i,11)
      fpscond(i,21)=azbz(i)*fsscond(i,07)+gsscond(i,12)
      fpscond(i,22)=axbx(i)*fsscond(i,08)+gsscond(i,08)
      fpscond(i,23)=ayby(i)*fsscond(i,08)+gsscond(i,12)
      fpscond(i,24)=azbz(i)*fsscond(i,08)+gsscond(i,13)
      fpscond(i,25)=axbx(i)*fsscond(i,09)+gsscond(i,09)
      fpscond(i,26)=ayby(i)*fsscond(i,09)+gsscond(i,13)
      fpscond(i,27)=azbz(i)*fsscond(i,09)+gsscond(i,14)
      fpscond(i,28)=axbx(i)*fsscond(i,10)+gsscond(i,10)
      fpscond(i,29)=ayby(i)*fsscond(i,10)+gsscond(i,14)
 1018 fpscond(i,30)=azbz(i)*fsscond(i,10)+gsscond(i,15)

      do 1019 i=1,m
      gpscond(i,01)=axbx(i)*gsscond(i,01)+hsscond(i,01)
      gpscond(i,02)=ayby(i)*gsscond(i,01)+hsscond(i,02)
      gpscond(i,03)=azbz(i)*gsscond(i,01)+hsscond(i,03)
      gpscond(i,04)=axbx(i)*gsscond(i,02)+hsscond(i,02)
      gpscond(i,05)=ayby(i)*gsscond(i,02)+hsscond(i,04)
      gpscond(i,06)=azbz(i)*gsscond(i,02)+hsscond(i,05)
      gpscond(i,07)=axbx(i)*gsscond(i,03)+hsscond(i,03)
      gpscond(i,08)=ayby(i)*gsscond(i,03)+hsscond(i,05)
      gpscond(i,09)=azbz(i)*gsscond(i,03)+hsscond(i,06)
      gpscond(i,10)=axbx(i)*gsscond(i,04)+hsscond(i,04)
      gpscond(i,11)=ayby(i)*gsscond(i,04)+hsscond(i,07)
      gpscond(i,12)=azbz(i)*gsscond(i,04)+hsscond(i,08)
      gpscond(i,13)=axbx(i)*gsscond(i,05)+hsscond(i,05)
      gpscond(i,14)=ayby(i)*gsscond(i,05)+hsscond(i,08)
      gpscond(i,15)=azbz(i)*gsscond(i,05)+hsscond(i,09)
      gpscond(i,16)=axbx(i)*gsscond(i,06)+hsscond(i,06)
      gpscond(i,17)=ayby(i)*gsscond(i,06)+hsscond(i,09)
      gpscond(i,18)=azbz(i)*gsscond(i,06)+hsscond(i,10)
      gpscond(i,19)=axbx(i)*gsscond(i,07)+hsscond(i,07)
      gpscond(i,20)=ayby(i)*gsscond(i,07)+hsscond(i,11)
      gpscond(i,21)=azbz(i)*gsscond(i,07)+hsscond(i,12)
      gpscond(i,22)=axbx(i)*gsscond(i,08)+hsscond(i,08)
      gpscond(i,23)=ayby(i)*gsscond(i,08)+hsscond(i,12)
      gpscond(i,24)=azbz(i)*gsscond(i,08)+hsscond(i,13)
      gpscond(i,25)=axbx(i)*gsscond(i,09)+hsscond(i,09)
      gpscond(i,26)=ayby(i)*gsscond(i,09)+hsscond(i,13)
      gpscond(i,27)=azbz(i)*gsscond(i,09)+hsscond(i,14)
      gpscond(i,28)=axbx(i)*gsscond(i,10)+hsscond(i,10)
      gpscond(i,29)=ayby(i)*gsscond(i,10)+hsscond(i,14)
      gpscond(i,30)=azbz(i)*gsscond(i,10)+hsscond(i,15)
      gpscond(i,31)=axbx(i)*gsscond(i,11)+hsscond(i,11)
      gpscond(i,32)=ayby(i)*gsscond(i,11)+hsscond(i,16)
      gpscond(i,33)=azbz(i)*gsscond(i,11)+hsscond(i,17)
      gpscond(i,34)=axbx(i)*gsscond(i,12)+hsscond(i,12)
      gpscond(i,35)=ayby(i)*gsscond(i,12)+hsscond(i,17)
      gpscond(i,36)=azbz(i)*gsscond(i,12)+hsscond(i,18)
      gpscond(i,37)=axbx(i)*gsscond(i,13)+hsscond(i,13)
      gpscond(i,38)=ayby(i)*gsscond(i,13)+hsscond(i,18)
      gpscond(i,39)=azbz(i)*gsscond(i,13)+hsscond(i,19)
      gpscond(i,40)=axbx(i)*gsscond(i,14)+hsscond(i,14)
      gpscond(i,41)=ayby(i)*gsscond(i,14)+hsscond(i,19)
      gpscond(i,42)=azbz(i)*gsscond(i,14)+hsscond(i,20)
      gpscond(i,43)=axbx(i)*gsscond(i,15)+hsscond(i,15)
      gpscond(i,44)=ayby(i)*gsscond(i,15)+hsscond(i,20)
 1019 gpscond(i,45)=azbz(i)*gsscond(i,15)+hsscond(i,21)

      do 1020 i=1,m
      ddscond(i,01)=axbx(i)*dpscond(i,01)+fpscond(i,01)
      ddscond(i,02)=ayby(i)*dpscond(i,01)+fpscond(i,04)
      ddscond(i,03)=azbz(i)*dpscond(i,01)+fpscond(i,07)
      ddscond(i,04)=ayby(i)*dpscond(i,02)+fpscond(i,05)
      ddscond(i,05)=azbz(i)*dpscond(i,02)+fpscond(i,08)
      ddscond(i,06)=azbz(i)*dpscond(i,03)+fpscond(i,09)
      ddscond(i,07)=axbx(i)*dpscond(i,04)+fpscond(i,04)
      ddscond(i,08)=ayby(i)*dpscond(i,04)+fpscond(i,10)
      ddscond(i,09)=azbz(i)*dpscond(i,04)+fpscond(i,13)
      ddscond(i,10)=ayby(i)*dpscond(i,05)+fpscond(i,11)
      ddscond(i,11)=azbz(i)*dpscond(i,05)+fpscond(i,14)
      ddscond(i,12)=azbz(i)*dpscond(i,06)+fpscond(i,15)
      ddscond(i,13)=axbx(i)*dpscond(i,07)+fpscond(i,07)
      ddscond(i,14)=ayby(i)*dpscond(i,07)+fpscond(i,13)
      ddscond(i,15)=azbz(i)*dpscond(i,07)+fpscond(i,16)
      ddscond(i,16)=ayby(i)*dpscond(i,08)+fpscond(i,14)
      ddscond(i,17)=azbz(i)*dpscond(i,08)+fpscond(i,17)
      ddscond(i,18)=azbz(i)*dpscond(i,09)+fpscond(i,18)
      ddscond(i,19)=axbx(i)*dpscond(i,10)+fpscond(i,10)
      ddscond(i,20)=ayby(i)*dpscond(i,10)+fpscond(i,19)
      ddscond(i,21)=azbz(i)*dpscond(i,10)+fpscond(i,22)
      ddscond(i,22)=ayby(i)*dpscond(i,11)+fpscond(i,20)
      ddscond(i,23)=azbz(i)*dpscond(i,11)+fpscond(i,23)
      ddscond(i,24)=azbz(i)*dpscond(i,12)+fpscond(i,24)
      ddscond(i,25)=axbx(i)*dpscond(i,13)+fpscond(i,13)
      ddscond(i,26)=ayby(i)*dpscond(i,13)+fpscond(i,22)
      ddscond(i,27)=azbz(i)*dpscond(i,13)+fpscond(i,25)
      ddscond(i,28)=ayby(i)*dpscond(i,14)+fpscond(i,23)
      ddscond(i,29)=azbz(i)*dpscond(i,14)+fpscond(i,26)
      ddscond(i,30)=azbz(i)*dpscond(i,15)+fpscond(i,27)
      ddscond(i,31)=axbx(i)*dpscond(i,16)+fpscond(i,16)
      ddscond(i,32)=ayby(i)*dpscond(i,16)+fpscond(i,25)
      ddscond(i,33)=azbz(i)*dpscond(i,16)+fpscond(i,28)
      ddscond(i,34)=ayby(i)*dpscond(i,17)+fpscond(i,26)
      ddscond(i,35)=azbz(i)*dpscond(i,17)+fpscond(i,29)
 1020 ddscond(i,36)=azbz(i)*dpscond(i,18)+fpscond(i,30)

      do 1021 i=1,m
      fdscond(i,01)=axbx(i)*fpscond(i,01)+gpscond(i,01)
      fdscond(i,02)=ayby(i)*fpscond(i,01)+gpscond(i,04)
      fdscond(i,03)=azbz(i)*fpscond(i,01)+gpscond(i,07)
      fdscond(i,04)=ayby(i)*fpscond(i,02)+gpscond(i,05)
      fdscond(i,05)=azbz(i)*fpscond(i,02)+gpscond(i,08)
      fdscond(i,06)=azbz(i)*fpscond(i,03)+gpscond(i,09)
      fdscond(i,07)=axbx(i)*fpscond(i,04)+gpscond(i,04)
      fdscond(i,08)=ayby(i)*fpscond(i,04)+gpscond(i,10)
      fdscond(i,09)=azbz(i)*fpscond(i,04)+gpscond(i,13)
      fdscond(i,10)=ayby(i)*fpscond(i,05)+gpscond(i,11)
      fdscond(i,11)=azbz(i)*fpscond(i,05)+gpscond(i,14)
      fdscond(i,12)=azbz(i)*fpscond(i,06)+gpscond(i,15)
      fdscond(i,13)=axbx(i)*fpscond(i,07)+gpscond(i,07)
      fdscond(i,14)=ayby(i)*fpscond(i,07)+gpscond(i,13)
      fdscond(i,15)=azbz(i)*fpscond(i,07)+gpscond(i,16)
      fdscond(i,16)=ayby(i)*fpscond(i,08)+gpscond(i,14)
      fdscond(i,17)=azbz(i)*fpscond(i,08)+gpscond(i,17)
      fdscond(i,18)=azbz(i)*fpscond(i,09)+gpscond(i,18)
      fdscond(i,19)=axbx(i)*fpscond(i,10)+gpscond(i,10)
      fdscond(i,20)=ayby(i)*fpscond(i,10)+gpscond(i,19)
      fdscond(i,21)=azbz(i)*fpscond(i,10)+gpscond(i,22)
      fdscond(i,22)=ayby(i)*fpscond(i,11)+gpscond(i,20)
      fdscond(i,23)=azbz(i)*fpscond(i,11)+gpscond(i,23)
      fdscond(i,24)=azbz(i)*fpscond(i,12)+gpscond(i,24)
      fdscond(i,25)=axbx(i)*fpscond(i,13)+gpscond(i,13)
      fdscond(i,26)=ayby(i)*fpscond(i,13)+gpscond(i,22)
      fdscond(i,27)=azbz(i)*fpscond(i,13)+gpscond(i,25)
      fdscond(i,28)=ayby(i)*fpscond(i,14)+gpscond(i,23)
      fdscond(i,29)=azbz(i)*fpscond(i,14)+gpscond(i,26)
      fdscond(i,30)=azbz(i)*fpscond(i,15)+gpscond(i,27)
      fdscond(i,31)=axbx(i)*fpscond(i,16)+gpscond(i,16)
      fdscond(i,32)=ayby(i)*fpscond(i,16)+gpscond(i,25)
      fdscond(i,33)=azbz(i)*fpscond(i,16)+gpscond(i,28)
      fdscond(i,34)=ayby(i)*fpscond(i,17)+gpscond(i,26)
      fdscond(i,35)=azbz(i)*fpscond(i,17)+gpscond(i,29)
      fdscond(i,36)=azbz(i)*fpscond(i,18)+gpscond(i,30)
      fdscond(i,37)=axbx(i)*fpscond(i,19)+gpscond(i,19)
      fdscond(i,38)=ayby(i)*fpscond(i,19)+gpscond(i,31)
      fdscond(i,39)=azbz(i)*fpscond(i,19)+gpscond(i,34)
      fdscond(i,40)=ayby(i)*fpscond(i,20)+gpscond(i,32)
      fdscond(i,41)=azbz(i)*fpscond(i,20)+gpscond(i,35)
      fdscond(i,42)=azbz(i)*fpscond(i,21)+gpscond(i,36)
      fdscond(i,43)=axbx(i)*fpscond(i,22)+gpscond(i,22)
      fdscond(i,44)=ayby(i)*fpscond(i,22)+gpscond(i,34)
      fdscond(i,45)=azbz(i)*fpscond(i,22)+gpscond(i,37)
      fdscond(i,46)=ayby(i)*fpscond(i,23)+gpscond(i,35)
      fdscond(i,47)=azbz(i)*fpscond(i,23)+gpscond(i,38)
      fdscond(i,48)=azbz(i)*fpscond(i,24)+gpscond(i,39)
      fdscond(i,49)=axbx(i)*fpscond(i,25)+gpscond(i,25)
      fdscond(i,50)=ayby(i)*fpscond(i,25)+gpscond(i,37)
      fdscond(i,51)=azbz(i)*fpscond(i,25)+gpscond(i,40)
      fdscond(i,52)=ayby(i)*fpscond(i,26)+gpscond(i,38)
      fdscond(i,53)=azbz(i)*fpscond(i,26)+gpscond(i,41)
      fdscond(i,54)=azbz(i)*fpscond(i,27)+gpscond(i,42)
      fdscond(i,55)=axbx(i)*fpscond(i,28)+gpscond(i,28)
      fdscond(i,56)=ayby(i)*fpscond(i,28)+gpscond(i,40)
      fdscond(i,57)=azbz(i)*fpscond(i,28)+gpscond(i,43)
      fdscond(i,58)=ayby(i)*fpscond(i,29)+gpscond(i,41)
      fdscond(i,59)=azbz(i)*fpscond(i,29)+gpscond(i,44)
 1021 fdscond(i,60)=azbz(i)*fpscond(i,30)+gpscond(i,45)

      do 1022 i=1,m
      dfscond(i,01)=axbx(i)*ddscond(i,01)+fdscond(i,01)
      dfscond(i,02)=ayby(i)*ddscond(i,01)+fdscond(i,07)
      dfscond(i,03)=azbz(i)*ddscond(i,01)+fdscond(i,13)
      dfscond(i,04)=ayby(i)*ddscond(i,02)+fdscond(i,08)
      dfscond(i,05)=azbz(i)*ddscond(i,02)+fdscond(i,14)
      dfscond(i,06)=azbz(i)*ddscond(i,03)+fdscond(i,15)
      dfscond(i,07)=ayby(i)*ddscond(i,04)+fdscond(i,10)
      dfscond(i,08)=azbz(i)*ddscond(i,04)+fdscond(i,16)
      dfscond(i,09)=azbz(i)*ddscond(i,05)+fdscond(i,17)
      dfscond(i,10)=azbz(i)*ddscond(i,06)+fdscond(i,18)
      dfscond(i,11)=axbx(i)*ddscond(i,07)+fdscond(i,07)
      dfscond(i,12)=ayby(i)*ddscond(i,07)+fdscond(i,19)
      dfscond(i,13)=azbz(i)*ddscond(i,07)+fdscond(i,25)
      dfscond(i,14)=ayby(i)*ddscond(i,08)+fdscond(i,20)
      dfscond(i,15)=azbz(i)*ddscond(i,08)+fdscond(i,26)
      dfscond(i,16)=azbz(i)*ddscond(i,09)+fdscond(i,27)
      dfscond(i,17)=ayby(i)*ddscond(i,10)+fdscond(i,22)
      dfscond(i,18)=azbz(i)*ddscond(i,10)+fdscond(i,28)
      dfscond(i,19)=azbz(i)*ddscond(i,11)+fdscond(i,29)
      dfscond(i,20)=azbz(i)*ddscond(i,12)+fdscond(i,30)
      dfscond(i,21)=axbx(i)*ddscond(i,13)+fdscond(i,13)
      dfscond(i,22)=ayby(i)*ddscond(i,13)+fdscond(i,25)
      dfscond(i,23)=azbz(i)*ddscond(i,13)+fdscond(i,31)
      dfscond(i,24)=ayby(i)*ddscond(i,14)+fdscond(i,26)
      dfscond(i,25)=azbz(i)*ddscond(i,14)+fdscond(i,32)
      dfscond(i,26)=azbz(i)*ddscond(i,15)+fdscond(i,33)
      dfscond(i,27)=ayby(i)*ddscond(i,16)+fdscond(i,28)
      dfscond(i,28)=azbz(i)*ddscond(i,16)+fdscond(i,34)
      dfscond(i,29)=azbz(i)*ddscond(i,17)+fdscond(i,35)
      dfscond(i,30)=azbz(i)*ddscond(i,18)+fdscond(i,36)
      dfscond(i,31)=axbx(i)*ddscond(i,19)+fdscond(i,19)
      dfscond(i,32)=ayby(i)*ddscond(i,19)+fdscond(i,37)
      dfscond(i,33)=azbz(i)*ddscond(i,19)+fdscond(i,43)
      dfscond(i,34)=ayby(i)*ddscond(i,20)+fdscond(i,38)
      dfscond(i,35)=azbz(i)*ddscond(i,20)+fdscond(i,44)
      dfscond(i,36)=azbz(i)*ddscond(i,21)+fdscond(i,45)
      dfscond(i,37)=ayby(i)*ddscond(i,22)+fdscond(i,40)
      dfscond(i,38)=azbz(i)*ddscond(i,22)+fdscond(i,46)
      dfscond(i,39)=azbz(i)*ddscond(i,23)+fdscond(i,47)
      dfscond(i,40)=azbz(i)*ddscond(i,24)+fdscond(i,48)
      dfscond(i,41)=axbx(i)*ddscond(i,25)+fdscond(i,25)
      dfscond(i,42)=ayby(i)*ddscond(i,25)+fdscond(i,43)
      dfscond(i,43)=azbz(i)*ddscond(i,25)+fdscond(i,49)
      dfscond(i,44)=ayby(i)*ddscond(i,26)+fdscond(i,44)
      dfscond(i,45)=azbz(i)*ddscond(i,26)+fdscond(i,50)
      dfscond(i,46)=azbz(i)*ddscond(i,27)+fdscond(i,51)
      dfscond(i,47)=ayby(i)*ddscond(i,28)+fdscond(i,46)
      dfscond(i,48)=azbz(i)*ddscond(i,28)+fdscond(i,52)
      dfscond(i,49)=azbz(i)*ddscond(i,29)+fdscond(i,53)
      dfscond(i,50)=azbz(i)*ddscond(i,30)+fdscond(i,54)
      dfscond(i,51)=axbx(i)*ddscond(i,31)+fdscond(i,31)
      dfscond(i,52)=ayby(i)*ddscond(i,31)+fdscond(i,49)
      dfscond(i,53)=azbz(i)*ddscond(i,31)+fdscond(i,55)
      dfscond(i,54)=ayby(i)*ddscond(i,32)+fdscond(i,50)
      dfscond(i,55)=azbz(i)*ddscond(i,32)+fdscond(i,56)
      dfscond(i,56)=azbz(i)*ddscond(i,33)+fdscond(i,57)
      dfscond(i,57)=ayby(i)*ddscond(i,34)+fdscond(i,52)
      dfscond(i,58)=azbz(i)*ddscond(i,34)+fdscond(i,58)
      dfscond(i,59)=azbz(i)*ddscond(i,35)+fdscond(i,59)
 1022 dfscond(i,60)=azbz(i)*ddscond(i,36)+fdscond(i,60)

      do 1023 i=1,m
      grad(i,001)=dfscond(i,01)*density(inddyy(i,01))
      grad(i,002)=dfscond(i,02)*density(inddyy(i,02))
      grad(i,003)=dfscond(i,03)*density(inddyy(i,03))
      grad(i,004)=dfscond(i,04)*density(inddyy(i,04))
      grad(i,005)=dfscond(i,05)*density(inddyy(i,05))
      grad(i,006)=dfscond(i,06)*density(inddyy(i,06))
      grad(i,007)=dfscond(i,11)*density(inddyy(i,07))
      grad(i,008)=dfscond(i,12)*density(inddyy(i,08))
      grad(i,009)=dfscond(i,13)*density(inddyy(i,09))
      grad(i,010)=dfscond(i,14)*density(inddyy(i,10))
      grad(i,011)=dfscond(i,15)*density(inddyy(i,11))
      grad(i,012)=dfscond(i,16)*density(inddyy(i,12))
      grad(i,013)=dfscond(i,21)*density(inddyy(i,13))
      grad(i,014)=dfscond(i,22)*density(inddyy(i,14))
      grad(i,015)=dfscond(i,23)*density(inddyy(i,15))
      grad(i,016)=dfscond(i,24)*density(inddyy(i,16))
      grad(i,017)=dfscond(i,25)*density(inddyy(i,17))
      grad(i,018)=dfscond(i,26)*density(inddyy(i,18))
      grad(i,019)=dfscond(i,31)*density(inddyy(i,19))
      grad(i,020)=dfscond(i,32)*density(inddyy(i,20))
      grad(i,021)=dfscond(i,33)*density(inddyy(i,21))
      grad(i,022)=dfscond(i,34)*density(inddyy(i,22))
      grad(i,023)=dfscond(i,35)*density(inddyy(i,23))
      grad(i,024)=dfscond(i,36)*density(inddyy(i,24))
      grad(i,025)=dfscond(i,41)*density(inddyy(i,25))
      grad(i,026)=dfscond(i,42)*density(inddyy(i,26))
      grad(i,027)=dfscond(i,43)*density(inddyy(i,27))
      grad(i,028)=dfscond(i,44)*density(inddyy(i,28))
      grad(i,029)=dfscond(i,45)*density(inddyy(i,29))
      grad(i,030)=dfscond(i,46)*density(inddyy(i,30))
      grad(i,031)=dfscond(i,51)*density(inddyy(i,31))
      grad(i,032)=dfscond(i,52)*density(inddyy(i,32))
      grad(i,033)=dfscond(i,53)*density(inddyy(i,33))
      grad(i,034)=dfscond(i,54)*density(inddyy(i,34))
      grad(i,035)=dfscond(i,55)*density(inddyy(i,35))
      grad(i,036)=dfscond(i,56)*density(inddyy(i,36))
      grad(i,037)=dfscond(i,02)*density(inddyy(i,01))
      grad(i,038)=dfscond(i,04)*density(inddyy(i,02))
      grad(i,039)=dfscond(i,05)*density(inddyy(i,03))
      grad(i,040)=dfscond(i,07)*density(inddyy(i,04))
      grad(i,041)=dfscond(i,08)*density(inddyy(i,05))
      grad(i,042)=dfscond(i,09)*density(inddyy(i,06))
      grad(i,043)=dfscond(i,12)*density(inddyy(i,07))
      grad(i,044)=dfscond(i,14)*density(inddyy(i,08))
      grad(i,045)=dfscond(i,15)*density(inddyy(i,09))
      grad(i,046)=dfscond(i,17)*density(inddyy(i,10))
      grad(i,047)=dfscond(i,18)*density(inddyy(i,11))
      grad(i,048)=dfscond(i,19)*density(inddyy(i,12))
      grad(i,049)=dfscond(i,22)*density(inddyy(i,13))
      grad(i,050)=dfscond(i,24)*density(inddyy(i,14))
      grad(i,051)=dfscond(i,25)*density(inddyy(i,15))
      grad(i,052)=dfscond(i,27)*density(inddyy(i,16))
      grad(i,053)=dfscond(i,28)*density(inddyy(i,17))
      grad(i,054)=dfscond(i,29)*density(inddyy(i,18))
      grad(i,055)=dfscond(i,32)*density(inddyy(i,19))
      grad(i,056)=dfscond(i,34)*density(inddyy(i,20))
      grad(i,057)=dfscond(i,35)*density(inddyy(i,21))
      grad(i,058)=dfscond(i,37)*density(inddyy(i,22))
      grad(i,059)=dfscond(i,38)*density(inddyy(i,23))
      grad(i,060)=dfscond(i,39)*density(inddyy(i,24))
      grad(i,061)=dfscond(i,42)*density(inddyy(i,25))
      grad(i,062)=dfscond(i,44)*density(inddyy(i,26))
      grad(i,063)=dfscond(i,45)*density(inddyy(i,27))
      grad(i,064)=dfscond(i,47)*density(inddyy(i,28))
      grad(i,065)=dfscond(i,48)*density(inddyy(i,29))
      grad(i,066)=dfscond(i,49)*density(inddyy(i,30))
      grad(i,067)=dfscond(i,52)*density(inddyy(i,31))
      grad(i,068)=dfscond(i,54)*density(inddyy(i,32))
      grad(i,069)=dfscond(i,55)*density(inddyy(i,33))
      grad(i,070)=dfscond(i,57)*density(inddyy(i,34))
      grad(i,071)=dfscond(i,58)*density(inddyy(i,35))
      grad(i,072)=dfscond(i,59)*density(inddyy(i,36))
      grad(i,073)=dfscond(i,03)*density(inddyy(i,01))
      grad(i,074)=dfscond(i,05)*density(inddyy(i,02))
      grad(i,075)=dfscond(i,06)*density(inddyy(i,03))
      grad(i,076)=dfscond(i,08)*density(inddyy(i,04))
      grad(i,077)=dfscond(i,09)*density(inddyy(i,05))
      grad(i,078)=dfscond(i,10)*density(inddyy(i,06))
      grad(i,079)=dfscond(i,13)*density(inddyy(i,07))
      grad(i,080)=dfscond(i,15)*density(inddyy(i,08))
      grad(i,081)=dfscond(i,16)*density(inddyy(i,09))
      grad(i,082)=dfscond(i,18)*density(inddyy(i,10))
      grad(i,083)=dfscond(i,19)*density(inddyy(i,11))
      grad(i,084)=dfscond(i,20)*density(inddyy(i,12))
      grad(i,085)=dfscond(i,23)*density(inddyy(i,13))
      grad(i,086)=dfscond(i,25)*density(inddyy(i,14))
      grad(i,087)=dfscond(i,26)*density(inddyy(i,15))
      grad(i,088)=dfscond(i,28)*density(inddyy(i,16))
      grad(i,089)=dfscond(i,29)*density(inddyy(i,17))
      grad(i,090)=dfscond(i,30)*density(inddyy(i,18))
      grad(i,091)=dfscond(i,33)*density(inddyy(i,19))
      grad(i,092)=dfscond(i,35)*density(inddyy(i,20))
      grad(i,093)=dfscond(i,36)*density(inddyy(i,21))
      grad(i,094)=dfscond(i,38)*density(inddyy(i,22))
      grad(i,095)=dfscond(i,39)*density(inddyy(i,23))
      grad(i,096)=dfscond(i,40)*density(inddyy(i,24))
      grad(i,097)=dfscond(i,43)*density(inddyy(i,25))
      grad(i,098)=dfscond(i,45)*density(inddyy(i,26))
      grad(i,099)=dfscond(i,46)*density(inddyy(i,27))
      grad(i,100)=dfscond(i,48)*density(inddyy(i,28))
      grad(i,101)=dfscond(i,49)*density(inddyy(i,29))
      grad(i,102)=dfscond(i,50)*density(inddyy(i,30))
      grad(i,103)=dfscond(i,53)*density(inddyy(i,31))
      grad(i,104)=dfscond(i,55)*density(inddyy(i,32))
      grad(i,105)=dfscond(i,56)*density(inddyy(i,33))
      grad(i,106)=dfscond(i,58)*density(inddyy(i,34))
      grad(i,107)=dfscond(i,59)*density(inddyy(i,35))
 1023 grad(i,108)=dfscond(i,60)*density(inddyy(i,36))

      do 1024 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,002)=zero
                              grad(i,003)=zero
                              grad(i,004)=zero
                              grad(i,005)=zero
                              grad(i,006)=zero
                              grad(i,009)=zero
                              grad(i,010)=zero
                              grad(i,011)=zero
                              grad(i,012)=zero
                              grad(i,016)=zero
                              grad(i,017)=zero
                              grad(i,018)=zero
                              grad(i,023)=zero
                              grad(i,024)=zero
                              grad(i,030)=zero
                              grad(i,038)=zero
                              grad(i,039)=zero
                              grad(i,040)=zero
                              grad(i,041)=zero
                              grad(i,042)=zero
                              grad(i,045)=zero
                              grad(i,046)=zero
                              grad(i,047)=zero
                              grad(i,048)=zero
                              grad(i,052)=zero
                              grad(i,053)=zero
                              grad(i,054)=zero
                              grad(i,059)=zero
                              grad(i,060)=zero
                              grad(i,066)=zero
                              grad(i,074)=zero
                              grad(i,075)=zero
                              grad(i,076)=zero
                              grad(i,077)=zero
                              grad(i,078)=zero
                              grad(i,081)=zero
                              grad(i,082)=zero
                              grad(i,083)=zero
                              grad(i,084)=zero
                              grad(i,088)=zero
                              grad(i,089)=zero
                              grad(i,090)=zero
                              grad(i,095)=zero
                              grad(i,096)=zero
                              grad(i,102)=zero
                        endif
 1024 continue

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

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

      return
      end
