      subroutine gradcddd(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,sss4,sss5,pss0,pss1,pss2,
     &                    pss3,pss4,dss0,dss1,dss2,dss3,fss0,fss1,fss2,
     &                    gss0,gss1,hss0,c,c1,pscond,dscond,ppcond,
     &                    fscond,dpcond,pdcond,gscond,fpcond,ddcond,
     &                    hscond,gpcond,fdcond,dfcond,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(*),sss4(*),sss5(*),
     &          pss0(n,*),pss1(n,*),pss2(n,*),pss3(n,*),pss4(n,*),
     &          dss0(n,*),dss1(n,*),dss2(n,*),dss3(n,*),fss0(n,*),
     &          fss1(n,*),fss2(n,*),gss0(n,*),gss1(n,*),hss0(n,*),c(*),
     &          c1(n,*),pscond(m,*),dscond(m,*),ppcond(m,*),fscond(m,*),
     &          dpcond(m,*),pdcond(m,*),gscond(m,*),fpcond(m,*),
     &          ddcond(m,*),hscond(m,*),gpcond(m,*),fdcond(m,*),
     &          dfcond(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(6,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
      sss3(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(4)/zeta(i)/zetaux
      sss4(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(5)/zeta(i)/zetaux
 1004 sss5(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(6)/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)
                                       sss4(nn)=sss4(i)
                                       sss5(nn)=sss5(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
      pss3(i,1)=ppxax(i)*sss3(i)+wxpx(i)*sss4(i)
      pss3(i,2)=ppyay(i)*sss3(i)+wypy(i)*sss4(i)
 1011 pss3(i,3)=ppzaz(i)*sss3(i)+wzpz(i)*sss4(i)

      do 1012 i=1,nn
      pss4(i,1)=ppxax(i)*sss4(i)+wxpx(i)*sss5(i)
      pss4(i,2)=ppyay(i)*sss4(i)+wypy(i)*sss5(i)
 1012 pss4(i,3)=ppzaz(i)*sss4(i)+wzpz(i)*sss5(i)

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

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

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

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

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

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

      do 1019 i=1,nn
 1019 c(i)=(sss3(i)-zetac1(i)*sss4(i))/(two*zzeta(i))

      do 1020 i=1,nn
      dss3(i,1)=ppxax(i)*pss3(i,1)+wxpx(i)*pss4(i,1)+c(i)
      dss3(i,2)=ppyay(i)*pss3(i,1)+wypy(i)*pss4(i,1)
      dss3(i,3)=ppzaz(i)*pss3(i,1)+wzpz(i)*pss4(i,1)
      dss3(i,4)=ppyay(i)*pss3(i,2)+wypy(i)*pss4(i,2)+c(i)
      dss3(i,5)=ppzaz(i)*pss3(i,2)+wzpz(i)*pss4(i,2)
 1020 dss3(i,6)=ppzaz(i)*pss3(i,3)+wzpz(i)*pss4(i,3)+c(i)

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

      do 1022 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)
 1022 fss0(i,10)=ppzaz(i)*dss0(i,6)+wzpz(i)*dss1(i,6)+c1(i,3)+c1(i,3)

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

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

      do 1025 j=1,3
      do 1025 i=1,nn
 1025 c1(i,j)=(pss2(i,j)-zetac1(i)*pss3(i,j))/(two*zzeta(i))

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

      do 1027 j=1,6
      do 1027 i=1,nn
 1027 c1(i,j)=(dss0(i,j)-zetac1(i)*dss1(i,j))/(two*zzeta(i))

      do 1028 i=1,nn
      gss0(i,01)=ppxax(i)*fss0(i,01)+wxpx(i)*fss1(i,01)+c1(i,1)*three
      gss0(i,02)=ppyay(i)*fss0(i,01)+wypy(i)*fss1(i,01)
      gss0(i,03)=ppzaz(i)*fss0(i,01)+wzpz(i)*fss1(i,01)
      gss0(i,04)=ppyay(i)*fss0(i,02)+wypy(i)*fss1(i,02)+c1(i,1)
      gss0(i,05)=ppzaz(i)*fss0(i,02)+wzpz(i)*fss1(i,02)
      gss0(i,06)=ppzaz(i)*fss0(i,03)+wzpz(i)*fss1(i,03)+c1(i,1)
      gss0(i,07)=ppyay(i)*fss0(i,04)+wypy(i)*fss1(i,04)+c1(i,2)*two
      gss0(i,08)=ppzaz(i)*fss0(i,04)+wzpz(i)*fss1(i,04)
      gss0(i,09)=ppzaz(i)*fss0(i,05)+wzpz(i)*fss1(i,05)+c1(i,2)
      gss0(i,10)=ppzaz(i)*fss0(i,06)+wzpz(i)*fss1(i,06)+c1(i,3)*two
      gss0(i,11)=ppyay(i)*fss0(i,07)+wypy(i)*fss1(i,07)+c1(i,4)*three
      gss0(i,12)=ppzaz(i)*fss0(i,07)+wzpz(i)*fss1(i,07)
      gss0(i,13)=ppzaz(i)*fss0(i,08)+wzpz(i)*fss1(i,08)+c1(i,4)
      gss0(i,14)=ppzaz(i)*fss0(i,09)+wzpz(i)*fss1(i,09)+c1(i,5)*two
 1028 gss0(i,15)=ppzaz(i)*fss0(i,10)+wzpz(i)*fss1(i,10)+c1(i,6)*three

      do 1029 j=1,6
      do 1029 i=1,nn
 1029 c1(i,j)=(dss1(i,j)-zetac1(i)*dss2(i,j))/(two*zzeta(i))

      do 1030 i=1,nn
      gss1(i,01)=ppxax(i)*fss1(i,01)+wxpx(i)*fss2(i,01)+c1(i,1)*three
      gss1(i,02)=ppyay(i)*fss1(i,01)+wypy(i)*fss2(i,01)
      gss1(i,03)=ppzaz(i)*fss1(i,01)+wzpz(i)*fss2(i,01)
      gss1(i,04)=ppyay(i)*fss1(i,02)+wypy(i)*fss2(i,02)+c1(i,1)
      gss1(i,05)=ppzaz(i)*fss1(i,02)+wzpz(i)*fss2(i,02)
      gss1(i,06)=ppzaz(i)*fss1(i,03)+wzpz(i)*fss2(i,03)+c1(i,1)
      gss1(i,07)=ppyay(i)*fss1(i,04)+wypy(i)*fss2(i,04)+c1(i,2)*two
      gss1(i,08)=ppzaz(i)*fss1(i,04)+wzpz(i)*fss2(i,04)
      gss1(i,09)=ppzaz(i)*fss1(i,05)+wzpz(i)*fss2(i,05)+c1(i,2)
      gss1(i,10)=ppzaz(i)*fss1(i,06)+wzpz(i)*fss2(i,06)+c1(i,3)*two
      gss1(i,11)=ppyay(i)*fss1(i,07)+wypy(i)*fss2(i,07)+c1(i,4)*three
      gss1(i,12)=ppzaz(i)*fss1(i,07)+wzpz(i)*fss2(i,07)
      gss1(i,13)=ppzaz(i)*fss1(i,08)+wzpz(i)*fss2(i,08)+c1(i,4)
      gss1(i,14)=ppzaz(i)*fss1(i,09)+wzpz(i)*fss2(i,09)+c1(i,5)*two
 1030 gss1(i,15)=ppzaz(i)*fss1(i,10)+wzpz(i)*fss2(i,10)+c1(i,6)*three

      do 1031 j=1,10
      do 1031 i=1,nn
 1031 c1(i,j)=(fss0(i,j)-zetac1(i)*fss1(i,j))/(two*zzeta(i))

      do 1032 i=1,nn
      hss0(i,01)=ppxax(i)*gss0(i,01)+wxpx(i)*gss1(i,01)+c1(i,01)*four
      hss0(i,02)=ppyay(i)*gss0(i,01)+wypy(i)*gss1(i,01)
      hss0(i,03)=ppzaz(i)*gss0(i,01)+wzpz(i)*gss1(i,01)
      hss0(i,04)=ppyay(i)*gss0(i,02)+wypy(i)*gss1(i,02)+c1(i,01)
      hss0(i,05)=ppzaz(i)*gss0(i,02)+wzpz(i)*gss1(i,02)
      hss0(i,06)=ppzaz(i)*gss0(i,03)+wzpz(i)*gss1(i,03)+c1(i,01)
      hss0(i,07)=ppyay(i)*gss0(i,04)+wypy(i)*gss1(i,04)+c1(i,02)*two
      hss0(i,08)=ppzaz(i)*gss0(i,04)+wzpz(i)*gss1(i,04)
      hss0(i,09)=ppzaz(i)*gss0(i,05)+wzpz(i)*gss1(i,05)+c1(i,02)
      hss0(i,10)=ppzaz(i)*gss0(i,06)+wzpz(i)*gss1(i,06)+c1(i,03)*two
      hss0(i,11)=ppyay(i)*gss0(i,07)+wypy(i)*gss1(i,07)+c1(i,04)*three
      hss0(i,12)=ppzaz(i)*gss0(i,07)+wzpz(i)*gss1(i,07)
      hss0(i,13)=ppzaz(i)*gss0(i,08)+wzpz(i)*gss1(i,08)+c1(i,04)
      hss0(i,14)=ppzaz(i)*gss0(i,09)+wzpz(i)*gss1(i,09)+c1(i,05)*two
      hss0(i,15)=ppzaz(i)*gss0(i,10)+wzpz(i)*gss1(i,10)+c1(i,06)*three
      hss0(i,16)=ppyay(i)*gss0(i,11)+wypy(i)*gss1(i,11)+c1(i,07)*four
      hss0(i,17)=ppzaz(i)*gss0(i,11)+wzpz(i)*gss1(i,11)
      hss0(i,18)=ppzaz(i)*gss0(i,12)+wzpz(i)*gss1(i,12)+c1(i,07)
      hss0(i,19)=ppzaz(i)*gss0(i,13)+wzpz(i)*gss1(i,13)+c1(i,08)*two
      hss0(i,20)=ppzaz(i)*gss0(i,14)+wzpz(i)*gss1(i,14)+c1(i,09)*three
 1032 hss0(i,21)=ppzaz(i)*gss0(i,15)+wzpz(i)*gss1(i,15)+c1(i,10)*four

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

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

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

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

      do 1037 j=1,6
      do 1037 i=1,nn
 1037 dscond(iicontr(i),j)=dscond(iicontr(i),j)-dss0(i,j)

      do 1038 j=1,10
      do 1038 i=1,nn
 1038 fscond(iicontr(i),j)=fscond(iicontr(i),j)-fss0(i,j)

      do 1039 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)
 1039 ppcond(i,9)=azbz(i)*pscond(i,3)+dscond(i,6)

      do 1040 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)
 1040 dpcond(i,18)=azbz(i)*dscond(i,6)+fscond(i,10)

      do 1041 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)
 1041 pdcond(i,18)=azbz(i)*ppcond(i,9)+dpcond(i,18)

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

      do 1043 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,02)=zero
                              grad(i,03)=zero
                              grad(i,04)=zero
                              grad(i,05)=zero
                              grad(i,06)=zero
                              grad(i,09)=zero
                              grad(i,10)=zero
                              grad(i,11)=zero
                              grad(i,12)=zero
                              grad(i,16)=zero
                              grad(i,17)=zero
                              grad(i,18)=zero
                              grad(i,21)=zero
                              grad(i,22)=zero
                              grad(i,23)=zero
                              grad(i,24)=zero
                              grad(i,29)=zero
                              grad(i,30)=zero
                              grad(i,36)=zero
                              grad(i,40)=zero
                              grad(i,41)=zero
                              grad(i,42)=zero
                              grad(i,48)=zero
                        endif
 1043 continue

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

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

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

      do 1047 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,02)=zero
                              grad(i,03)=zero
                              grad(i,06)=zero
                              grad(i,19)=zero
                              grad(i,20)=zero
                              grad(i,21)=zero
                              grad(i,23)=zero
                              grad(i,24)=zero
                              grad(i,26)=zero
                              grad(i,27)=zero
                              grad(i,30)=zero
                              grad(i,37)=zero
                              grad(i,38)=zero
                              grad(i,39)=zero
                              grad(i,40)=zero
                              grad(i,41)=zero
                              grad(i,42)=zero
                              grad(i,44)=zero
                              grad(i,45)=zero
                              grad(i,47)=zero
                              grad(i,48)=zero
                              grad(i,51)=zero
                        endif
 1047 continue

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

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

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

      do 1051 j=1,15
      do 1051 i=1,m
 1051 gscond(i,j)=zero

      do 1052 j=1,21
      do 1052 i=1,m
 1052 hscond(i,j)=zero

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

      do 1054 j=1,15
      do 1054 i=1,nn
 1054 gscond(iicontr(i),j)=gscond(iicontr(i),j)+gss0(i,j)*ttwozetaa(i)

      do 1055 j=1,21
      do 1055 i=1,nn
 1055 hscond(iicontr(i),j)=hscond(iicontr(i),j)+hss0(i,j)*ttwozetaa(i)

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

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

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

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

      do 1060 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
 1060 continue

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

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

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

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

      do 1065 j=1,15
      do 1065 i=1,m
 1065 gscond(i,j)=zero

      do 1066 j=1,21
      do 1066 i=1,m
 1066 hscond(i,j)=zero

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

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

      do 1069 j=1,15
      do 1069 i=1,nn
 1069 gscond(iicontr(i),j)=gscond(iicontr(i),j)+gss0(i,j)*ttwozetab(i)

      do 1070 j=1,21
      do 1070 i=1,nn
 1070 hscond(iicontr(i),j)=hscond(iicontr(i),j)+hss0(i,j)*ttwozetab(i)

      do 1071 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)
 1071 dpcond(i,18)=azbz(i)*dscond(i,6)+fscond(i,10)

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

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

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

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

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

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

      do 1078 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
 1078 continue

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

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

 1001 continue

      return
      end
