      subroutine coredd(m,n,ncentr,iminj,icontr,inddyy,iloc,coord,
     &                  charge,overlap,core,zeta,zetaa,zetab,axbx,ayby,
     &                  azbz,px,py,pz,pxax,pyay,pzaz,pxbx,pyby,pzbz,eta,
     &                  base,etarr,ss,ss0,ss1,ss2,ss3,ss4,ps,ps0,ps1,
     &                  ps2,ps3,ds0,ds1,ds2,fs0,fs1,gs0,pp0,dpcnd,dp0,
     &                  fpcnd,ddcond,dd0,c,c1,c2,c3,pxcx,pycy,pzcz,u,
     &                  twozeta,eta2)

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(*),icontr(*),inddyy(m,*),iloc(*)

      dimension coord(3,*),charge(*),overlap(*),core(*),zeta(*),
     &          zetaa(*),zetab(*),axbx(*),ayby(*),azbz(*),px(*),py(*),
     &          pz(*),pxax(*),pyay(*),pzaz(*),pxbx(*),pyby(*),pzbz(*),
     &          eta(*),base(*),etarr(*),ss(*),ss0(*),ss1(*),ss2(*),
     &          ss3(*),ss4(*),ps(n,*),ps0(n,*),ps1(n,*),ps2(n,*),
     &          ps3(n,*),ds0(n,*),ds1(n,*),ds2(n,*),fs0(n,*),fs1(n,*),
     &          gs0(n,*),pp0(n,*),dpcnd(n,*),dp0(n,*),fpcnd(n,*),
     &          ddcond(n,*),dd0(n,*),c(*),c1(n,*),c2(n,*),c3(n,*),
     &          pxcx(*),pycy(*),pzcz(*),u(*),twozeta(*),eta2(*)

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

      sqrt3=sqrt(three)

      do 1001 i=1,n
      eta2(i)=two*eta(i)
 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
      ds0(i,1)=pxax(i)*ps0(i,1)+c(i)
      ds0(i,2)=pyay(i)*ps0(i,1)
      ds0(i,3)=pzaz(i)*ps0(i,1)
      ds0(i,4)=pyay(i)*ps0(i,2)+c(i)
      ds0(i,5)=pzaz(i)*ps0(i,2)
 1005 ds0(i,6)=pzaz(i)*ps0(i,3)+c(i)

      do 1006 i=1,n
      pp0(i,1)=pxbx(i)*ps0(i,1)+c(i)
      pp0(i,2)=pyby(i)*ps0(i,1)
      pp0(i,3)=pzbz(i)*ps0(i,1)
      pp0(i,4)=pxbx(i)*ps0(i,2)
      pp0(i,5)=pyby(i)*ps0(i,2)+c(i)
      pp0(i,6)=pzbz(i)*ps0(i,2)
      pp0(i,7)=pxbx(i)*ps0(i,3)
      pp0(i,8)=pyby(i)*ps0(i,3)
 1006 pp0(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
      dp0(i,01)=pxbx(i)*ds0(i,1)+c1(i,1)*two
      dp0(i,02)=pyby(i)*ds0(i,1)
      dp0(i,03)=pzbz(i)*ds0(i,1)
      dp0(i,04)=pxbx(i)*ds0(i,2)+c1(i,2)
      dp0(i,05)=pyby(i)*ds0(i,2)+c1(i,1)
      dp0(i,06)=pzbz(i)*ds0(i,2)
      dp0(i,07)=pxbx(i)*ds0(i,3)+c1(i,3)
      dp0(i,08)=pyby(i)*ds0(i,3)
      dp0(i,09)=pzbz(i)*ds0(i,3)+c1(i,1)
      dp0(i,10)=pxbx(i)*ds0(i,4)
      dp0(i,11)=pyby(i)*ds0(i,4)+c1(i,2)*two
      dp0(i,12)=pzbz(i)*ds0(i,4)
      dp0(i,13)=pxbx(i)*ds0(i,5)
      dp0(i,14)=pyby(i)*ds0(i,5)+c1(i,3)
      dp0(i,15)=pzbz(i)*ds0(i,5)+c1(i,2)
      dp0(i,16)=pxbx(i)*ds0(i,6)
      dp0(i,17)=pyby(i)*ds0(i,6)
 1008 dp0(i,18)=pzbz(i)*ds0(i,6)+c1(i,3)*two

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

      do 1010 j=1,9
      do 1010 i=1,n
 1010 c3(i,j)=pp0(i,j)/twozeta(i)

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

      do 1012 i=1,n
      overlap(inddyy(icontr(i),01))=overlap(inddyy(icontr(i),01))
     &                             +dd0(i,01)
      overlap(inddyy(icontr(i),07))=overlap(inddyy(icontr(i),07))
     &                             +dd0(i,07)
      overlap(inddyy(icontr(i),08))=overlap(inddyy(icontr(i),08))
     &                             +dd0(i,08)
      overlap(inddyy(icontr(i),13))=overlap(inddyy(icontr(i),13))
     &                             +dd0(i,13)
      overlap(inddyy(icontr(i),14))=overlap(inddyy(icontr(i),14))
     &                             +dd0(i,14)
      overlap(inddyy(icontr(i),15))=overlap(inddyy(icontr(i),15))
     &                             +dd0(i,15)
      overlap(inddyy(icontr(i),19))=overlap(inddyy(icontr(i),19))
     &                             +dd0(i,19)
      overlap(inddyy(icontr(i),20))=overlap(inddyy(icontr(i),20))
     &                             +dd0(i,20)
      overlap(inddyy(icontr(i),21))=overlap(inddyy(icontr(i),21))
     &                             +dd0(i,21)
      overlap(inddyy(icontr(i),22))=overlap(inddyy(icontr(i),22))
     &                             +dd0(i,22)
      overlap(inddyy(icontr(i),25))=overlap(inddyy(icontr(i),25))
     &                             +dd0(i,25)
      overlap(inddyy(icontr(i),26))=overlap(inddyy(icontr(i),26))
     &                             +dd0(i,26)
      overlap(inddyy(icontr(i),27))=overlap(inddyy(icontr(i),27))
     &                             +dd0(i,27)
      overlap(inddyy(icontr(i),28))=overlap(inddyy(icontr(i),28))
     &                             +dd0(i,28)
      overlap(inddyy(icontr(i),29))=overlap(inddyy(icontr(i),29))
     &                             +dd0(i,29)
      overlap(inddyy(icontr(i),31))=overlap(inddyy(icontr(i),31))
     &                             +dd0(i,31)
      overlap(inddyy(icontr(i),32))=overlap(inddyy(icontr(i),32))
     &                             +dd0(i,32)
      overlap(inddyy(icontr(i),33))=overlap(inddyy(icontr(i),33))
     &                             +dd0(i,33)
      overlap(inddyy(icontr(i),34))=overlap(inddyy(icontr(i),34))
     &                             +dd0(i,34)
      overlap(inddyy(icontr(i),35))=overlap(inddyy(icontr(i),35))
     &                             +dd0(i,35)
 1012 overlap(inddyy(icontr(i),36))=overlap(inddyy(icontr(i),36))
     &                             +dd0(i,36)

      do 1013 i=1,n
      if(iminj(icontr(i)).ne.0) then

             overlap(inddyy(icontr(i),02))=overlap(inddyy(icontr(i),02))
     &                                    +dd0(i,02)
             overlap(inddyy(icontr(i),03))=overlap(inddyy(icontr(i),03))
     &                                    +dd0(i,03)
             overlap(inddyy(icontr(i),04))=overlap(inddyy(icontr(i),04))
     &                                    +dd0(i,04)
             overlap(inddyy(icontr(i),05))=overlap(inddyy(icontr(i),05))
     &                                    +dd0(i,05)
             overlap(inddyy(icontr(i),06))=overlap(inddyy(icontr(i),06))
     &                                    +dd0(i,06)
             overlap(inddyy(icontr(i),09))=overlap(inddyy(icontr(i),09))
     &                                    +dd0(i,09)
             overlap(inddyy(icontr(i),10))=overlap(inddyy(icontr(i),10))
     &                                    +dd0(i,10)
             overlap(inddyy(icontr(i),11))=overlap(inddyy(icontr(i),11))
     &                                    +dd0(i,11)
             overlap(inddyy(icontr(i),12))=overlap(inddyy(icontr(i),12))
     &                                    +dd0(i,12)
             overlap(inddyy(icontr(i),16))=overlap(inddyy(icontr(i),16))
     &                                    +dd0(i,16)
             overlap(inddyy(icontr(i),17))=overlap(inddyy(icontr(i),17))
     &                                    +dd0(i,17)
             overlap(inddyy(icontr(i),18))=overlap(inddyy(icontr(i),18))
     &                                    +dd0(i,18)
             overlap(inddyy(icontr(i),23))=overlap(inddyy(icontr(i),23))
     &                                    +dd0(i,23)
             overlap(inddyy(icontr(i),24))=overlap(inddyy(icontr(i),24))
     &                                    +dd0(i,24)
             overlap(inddyy(icontr(i),30))=overlap(inddyy(icontr(i),30))
     &                                    +dd0(i,30)

                        endif
 1013 continue

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

      do 1015 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)
 1015 ps(i,3)=pzaz(i)*ss(i)+eta2(i)*ps0(i,3)

      do 1016 i=1,n
 1016 c(i)=ss(i)/twozeta(i)-eta2(i)*ss0(i)/(two*zetaa(i))

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

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

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

      do 1020 j=1,3
      do 1020 i=1,n
 1020 c1(i,j)=ps(i,j)/twozeta(i)

      do 1021 i=1,n
      dp0(i,01)=pxbx(i)*ds1(i,1)+eta2(i)*dp0(i,01)+c1(i,1)*two
      dp0(i,02)=pyby(i)*ds1(i,1)+eta2(i)*dp0(i,02)
      dp0(i,03)=pzbz(i)*ds1(i,1)+eta2(i)*dp0(i,03)
      dp0(i,04)=pxbx(i)*ds1(i,2)+eta2(i)*dp0(i,04)+c1(i,2)
      dp0(i,05)=pyby(i)*ds1(i,2)+eta2(i)*dp0(i,05)+c1(i,1)
      dp0(i,06)=pzbz(i)*ds1(i,2)+eta2(i)*dp0(i,06)
      dp0(i,07)=pxbx(i)*ds1(i,3)+eta2(i)*dp0(i,07)+c1(i,3)
      dp0(i,08)=pyby(i)*ds1(i,3)+eta2(i)*dp0(i,08)
      dp0(i,09)=pzbz(i)*ds1(i,3)+eta2(i)*dp0(i,09)+c1(i,1)
      dp0(i,10)=pxbx(i)*ds1(i,4)+eta2(i)*dp0(i,10)
      dp0(i,11)=pyby(i)*ds1(i,4)+eta2(i)*dp0(i,11)+c1(i,2)*two
      dp0(i,12)=pzbz(i)*ds1(i,4)+eta2(i)*dp0(i,12)
      dp0(i,13)=pxbx(i)*ds1(i,5)+eta2(i)*dp0(i,13)
      dp0(i,14)=pyby(i)*ds1(i,5)+eta2(i)*dp0(i,14)+c1(i,3)
      dp0(i,15)=pzbz(i)*ds1(i,5)+eta2(i)*dp0(i,15)+c1(i,2)
      dp0(i,16)=pxbx(i)*ds1(i,6)+eta2(i)*dp0(i,16)
      dp0(i,17)=pyby(i)*ds1(i,6)+eta2(i)*dp0(i,17)
 1021 dp0(i,18)=pzbz(i)*ds1(i,6)+eta2(i)*dp0(i,18)+c1(i,3)*two

      do 1022 j=1,6
      do 1022 i=1,n
 1022 c2(i,j)=ds1(i,j)/twozeta(i)-eta2(i)*ds0(i,j)/(two*zetab(i))

      do 1023 j=1,9
      do 1023 i=1,n
 1023 c3(i,j)=pp0(i,j)/twozeta(i)

      do 1024 i=1,n
      dd0(i,01)=pxbx(i)*dp0(i,01)+eta2(i)*dd0(i,01)+c2(i,1)+c3(i,1)*two
      dd0(i,02)=pyby(i)*dp0(i,01)+eta2(i)*dd0(i,02)
      dd0(i,03)=pzbz(i)*dp0(i,01)+eta2(i)*dd0(i,03)
      dd0(i,04)=pyby(i)*dp0(i,02)+eta2(i)*dd0(i,04)+c2(i,1)
      dd0(i,05)=pzbz(i)*dp0(i,02)+eta2(i)*dd0(i,05)
      dd0(i,06)=pzbz(i)*dp0(i,03)+eta2(i)*dd0(i,06)+c2(i,1)
      dd0(i,07)=pxbx(i)*dp0(i,04)+eta2(i)*dd0(i,07)+c2(i,2)+c3(i,4)
      dd0(i,08)=pyby(i)*dp0(i,04)+eta2(i)*dd0(i,08)        +c3(i,1)
      dd0(i,09)=pzbz(i)*dp0(i,04)+eta2(i)*dd0(i,09)
      dd0(i,10)=pyby(i)*dp0(i,05)+eta2(i)*dd0(i,10)+c2(i,2)+c3(i,2)
      dd0(i,11)=pzbz(i)*dp0(i,05)+eta2(i)*dd0(i,11)
      dd0(i,12)=pzbz(i)*dp0(i,06)+eta2(i)*dd0(i,12)+c2(i,2)
      dd0(i,13)=pxbx(i)*dp0(i,07)+eta2(i)*dd0(i,13)+c2(i,3)+c3(i,7)
      dd0(i,14)=pyby(i)*dp0(i,07)+eta2(i)*dd0(i,14)
      dd0(i,15)=pzbz(i)*dp0(i,07)+eta2(i)*dd0(i,15)        +c3(i,1)
      dd0(i,16)=pyby(i)*dp0(i,08)+eta2(i)*dd0(i,16)+c2(i,3)
      dd0(i,17)=pzbz(i)*dp0(i,08)+eta2(i)*dd0(i,17)        +c3(i,2)
      dd0(i,18)=pzbz(i)*dp0(i,09)+eta2(i)*dd0(i,18)+c2(i,3)+c3(i,3)
      dd0(i,19)=pxbx(i)*dp0(i,10)+eta2(i)*dd0(i,19)+c2(i,4)
      dd0(i,20)=pyby(i)*dp0(i,10)+eta2(i)*dd0(i,20)        +c3(i,4)*two
      dd0(i,21)=pzbz(i)*dp0(i,10)+eta2(i)*dd0(i,21)
      dd0(i,22)=pyby(i)*dp0(i,11)+eta2(i)*dd0(i,22)+c2(i,4)+c3(i,5)*two
      dd0(i,23)=pzbz(i)*dp0(i,11)+eta2(i)*dd0(i,23)
      dd0(i,24)=pzbz(i)*dp0(i,12)+eta2(i)*dd0(i,24)+c2(i,4)
      dd0(i,25)=pxbx(i)*dp0(i,13)+eta2(i)*dd0(i,25)+c2(i,5)
      dd0(i,26)=pyby(i)*dp0(i,13)+eta2(i)*dd0(i,26)        +c3(i,7)
      dd0(i,27)=pzbz(i)*dp0(i,13)+eta2(i)*dd0(i,27)        +c3(i,4)
      dd0(i,28)=pyby(i)*dp0(i,14)+eta2(i)*dd0(i,28)+c2(i,5)+c3(i,8)
      dd0(i,29)=pzbz(i)*dp0(i,14)+eta2(i)*dd0(i,29)        +c3(i,5)
      dd0(i,30)=pzbz(i)*dp0(i,15)+eta2(i)*dd0(i,30)+c2(i,5)+c3(i,6)
      dd0(i,31)=pxbx(i)*dp0(i,16)+eta2(i)*dd0(i,31)+c2(i,6)
      dd0(i,32)=pyby(i)*dp0(i,16)+eta2(i)*dd0(i,32)
      dd0(i,33)=pzbz(i)*dp0(i,16)+eta2(i)*dd0(i,33)        +c3(i,7)*two
      dd0(i,34)=pyby(i)*dp0(i,17)+eta2(i)*dd0(i,34)+c2(i,6)
      dd0(i,35)=pzbz(i)*dp0(i,17)+eta2(i)*dd0(i,35)        +c3(i,8)*two
 1024 dd0(i,36)=pzbz(i)*dp0(i,18)+eta2(i)*dd0(i,36)+c2(i,6)+c3(i,9)*two

      do 1025 i=1,n
      core(inddyy(icontr(i),01))=core(inddyy(icontr(i),01))+dd0(i,01)
      core(inddyy(icontr(i),07))=core(inddyy(icontr(i),07))+dd0(i,07)
      core(inddyy(icontr(i),08))=core(inddyy(icontr(i),08))+dd0(i,08)
      core(inddyy(icontr(i),13))=core(inddyy(icontr(i),13))+dd0(i,13)
      core(inddyy(icontr(i),14))=core(inddyy(icontr(i),14))+dd0(i,14)
      core(inddyy(icontr(i),15))=core(inddyy(icontr(i),15))+dd0(i,15)
      core(inddyy(icontr(i),19))=core(inddyy(icontr(i),19))+dd0(i,19)
      core(inddyy(icontr(i),20))=core(inddyy(icontr(i),20))+dd0(i,20)
      core(inddyy(icontr(i),21))=core(inddyy(icontr(i),21))+dd0(i,21)
      core(inddyy(icontr(i),22))=core(inddyy(icontr(i),22))+dd0(i,22)
      core(inddyy(icontr(i),25))=core(inddyy(icontr(i),25))+dd0(i,25)
      core(inddyy(icontr(i),26))=core(inddyy(icontr(i),26))+dd0(i,26)
      core(inddyy(icontr(i),27))=core(inddyy(icontr(i),27))+dd0(i,27)
      core(inddyy(icontr(i),28))=core(inddyy(icontr(i),28))+dd0(i,28)
      core(inddyy(icontr(i),29))=core(inddyy(icontr(i),29))+dd0(i,29)
      core(inddyy(icontr(i),31))=core(inddyy(icontr(i),31))+dd0(i,31)
      core(inddyy(icontr(i),32))=core(inddyy(icontr(i),32))+dd0(i,32)
      core(inddyy(icontr(i),33))=core(inddyy(icontr(i),33))+dd0(i,33)
      core(inddyy(icontr(i),34))=core(inddyy(icontr(i),34))+dd0(i,34)
      core(inddyy(icontr(i),35))=core(inddyy(icontr(i),35))+dd0(i,35)
 1025 core(inddyy(icontr(i),36))=core(inddyy(icontr(i),36))+dd0(i,36)

      do 1026 i=1,n
      if(iminj(icontr(i)).ne.0) then

             core(inddyy(icontr(i),02))=core(inddyy(icontr(i),02))
     &                                 +dd0(i,02)
             core(inddyy(icontr(i),03))=core(inddyy(icontr(i),03))
     &                                 +dd0(i,03)
             core(inddyy(icontr(i),04))=core(inddyy(icontr(i),04))
     &                                 +dd0(i,04)
             core(inddyy(icontr(i),05))=core(inddyy(icontr(i),05))
     &                                 +dd0(i,05)
             core(inddyy(icontr(i),06))=core(inddyy(icontr(i),06))
     &                                 +dd0(i,06)
             core(inddyy(icontr(i),09))=core(inddyy(icontr(i),09))
     &                                 +dd0(i,09)
             core(inddyy(icontr(i),10))=core(inddyy(icontr(i),10))
     &                                 +dd0(i,10)
             core(inddyy(icontr(i),11))=core(inddyy(icontr(i),11))
     &                                 +dd0(i,11)
             core(inddyy(icontr(i),12))=core(inddyy(icontr(i),12))
     &                                 +dd0(i,12)
             core(inddyy(icontr(i),16))=core(inddyy(icontr(i),16))
     &                                 +dd0(i,16)
             core(inddyy(icontr(i),17))=core(inddyy(icontr(i),17))
     &                                 +dd0(i,17)
             core(inddyy(icontr(i),18))=core(inddyy(icontr(i),18))
     &                                 +dd0(i,18)
             core(inddyy(icontr(i),23))=core(inddyy(icontr(i),23))
     &                                 +dd0(i,23)
             core(inddyy(icontr(i),24))=core(inddyy(icontr(i),24))
     &                                 +dd0(i,24)
             core(inddyy(icontr(i),30))=core(inddyy(icontr(i),30))
     &                                 +dd0(i,30)

                        endif
 1026 continue

      do 1027 j=1,36
      do 1027 i=1,n
 1027 ddcond(i,j)=zero

      do 2001 nn=1,ncentr

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

      twopichrg=two*pi*charge(nn)

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

      call auggv(5,n,iloc,u,ds0,ss0,ss1,ss2)

      do 2003 i=1,n
      ss0(i)=base(i)*ds0(i,1)*twopichrg/zeta(i)
      ss1(i)=base(i)*ds0(i,2)*twopichrg/zeta(i)
      ss2(i)=base(i)*ds0(i,3)*twopichrg/zeta(i)
      ss3(i)=base(i)*ds0(i,4)*twopichrg/zeta(i)
 2003 ss4(i)=base(i)*ds0(i,5)*twopichrg/zeta(i)

      do 2004 i=1,n

      ps0(i,1)=pxax(i)*ss0(i)-pxcx(i)*ss1(i)
      ps0(i,2)=pyay(i)*ss0(i)-pycy(i)*ss1(i)
      ps0(i,3)=pzaz(i)*ss0(i)-pzcz(i)*ss1(i)

      ps1(i,1)=pxax(i)*ss1(i)-pxcx(i)*ss2(i)
      ps1(i,2)=pyay(i)*ss1(i)-pycy(i)*ss2(i)
      ps1(i,3)=pzaz(i)*ss1(i)-pzcz(i)*ss2(i)

      ps2(i,1)=pxax(i)*ss2(i)-pxcx(i)*ss3(i)
      ps2(i,2)=pyay(i)*ss2(i)-pycy(i)*ss3(i)
      ps2(i,3)=pzaz(i)*ss2(i)-pzcz(i)*ss3(i)

      ps3(i,1)=pxax(i)*ss3(i)-pxcx(i)*ss4(i)
      ps3(i,2)=pyay(i)*ss3(i)-pycy(i)*ss4(i)
      ps3(i,3)=pzaz(i)*ss3(i)-pzcz(i)*ss4(i)

 2004 continue

      do 2005 i=1,n
 2005 c(i)=(ss0(i)-ss1(i))/twozeta(i)

      do 2006 i=1,n
      ds0(i,1)=pxax(i)*ps0(i,1)-pxcx(i)*ps1(i,1)+c(i)
      ds0(i,2)=pyay(i)*ps0(i,1)-pycy(i)*ps1(i,1)
      ds0(i,3)=pzaz(i)*ps0(i,1)-pzcz(i)*ps1(i,1)
      ds0(i,4)=pyay(i)*ps0(i,2)-pycy(i)*ps1(i,2)+c(i)
      ds0(i,5)=pzaz(i)*ps0(i,2)-pzcz(i)*ps1(i,2)
 2006 ds0(i,6)=pzaz(i)*ps0(i,3)-pzcz(i)*ps1(i,3)+c(i)

      do 2007 i=1,n
 2007 c(i)=(ss1(i)-ss2(i))/twozeta(i)

      do 2008 i=1,n
      ds1(i,1)=pxax(i)*ps1(i,1)-pxcx(i)*ps2(i,1)+c(i)
      ds1(i,2)=pyay(i)*ps1(i,1)-pycy(i)*ps2(i,1)
      ds1(i,3)=pzaz(i)*ps1(i,1)-pzcz(i)*ps2(i,1)
      ds1(i,4)=pyay(i)*ps1(i,2)-pycy(i)*ps2(i,2)+c(i)
      ds1(i,5)=pzaz(i)*ps1(i,2)-pzcz(i)*ps2(i,2)
 2008 ds1(i,6)=pzaz(i)*ps1(i,3)-pzcz(i)*ps2(i,3)+c(i)

      do 2009 i=1,n
 2009 c(i)=(ss2(i)-ss3(i))/twozeta(i)

      do 2010 i=1,n
      ds2(i,1)=pxax(i)*ps2(i,1)-pxcx(i)*ps3(i,1)+c(i)
      ds2(i,2)=pyay(i)*ps2(i,1)-pycy(i)*ps3(i,1)
      ds2(i,3)=pzaz(i)*ps2(i,1)-pzcz(i)*ps3(i,1)
      ds2(i,4)=pyay(i)*ps2(i,2)-pycy(i)*ps3(i,2)+c(i)
      ds2(i,5)=pzaz(i)*ps2(i,2)-pzcz(i)*ps3(i,2)
 2010 ds2(i,6)=pzaz(i)*ps2(i,3)-pzcz(i)*ps3(i,3)+c(i)

      do 2011 j=1,3
      do 2011 i=1,n
 2011 c1(i,j)=(ps0(i,j)-ps1(i,j))/twozeta(i)

      do 2012 i=1,n
      fs0(i,01)=pxax(i)*ds0(i,1)-pxcx(i)*ds1(i,1)+c1(i,1)*two
      fs0(i,02)=pyay(i)*ds0(i,1)-pycy(i)*ds1(i,1)
      fs0(i,03)=pzaz(i)*ds0(i,1)-pzcz(i)*ds1(i,1)
      fs0(i,04)=pyay(i)*ds0(i,2)-pycy(i)*ds1(i,2)+c1(i,1)
      fs0(i,05)=pzaz(i)*ds0(i,2)-pzcz(i)*ds1(i,2)
      fs0(i,06)=pzaz(i)*ds0(i,3)-pzcz(i)*ds1(i,3)+c1(i,1)
      fs0(i,07)=pyay(i)*ds0(i,4)-pycy(i)*ds1(i,4)+c1(i,2)*two
      fs0(i,08)=pzaz(i)*ds0(i,4)-pzcz(i)*ds1(i,4)
      fs0(i,09)=pzaz(i)*ds0(i,5)-pzcz(i)*ds1(i,5)+c1(i,2)
 2012 fs0(i,10)=pzaz(i)*ds0(i,6)-pzcz(i)*ds1(i,6)+c1(i,3)*two

      do 2013 j=1,3
      do 2013 i=1,n
 2013 c1(i,j)=(ps1(i,j)-ps2(i,j))/twozeta(i)

      do 2014 i=1,n
      fs1(i,01)=pxax(i)*ds1(i,1)-pxcx(i)*ds2(i,1)+c1(i,1)*two
      fs1(i,02)=pyay(i)*ds1(i,1)-pycy(i)*ds2(i,1)
      fs1(i,03)=pzaz(i)*ds1(i,1)-pzcz(i)*ds2(i,1)
      fs1(i,04)=pyay(i)*ds1(i,2)-pycy(i)*ds2(i,2)+c1(i,1)
      fs1(i,05)=pzaz(i)*ds1(i,2)-pzcz(i)*ds2(i,2)
      fs1(i,06)=pzaz(i)*ds1(i,3)-pzcz(i)*ds2(i,3)+c1(i,1)
      fs1(i,07)=pyay(i)*ds1(i,4)-pycy(i)*ds2(i,4)+c1(i,2)*two
      fs1(i,08)=pzaz(i)*ds1(i,4)-pzcz(i)*ds2(i,4)
      fs1(i,09)=pzaz(i)*ds1(i,5)-pzcz(i)*ds2(i,5)+c1(i,2)
 2014 fs1(i,10)=pzaz(i)*ds1(i,6)-pzcz(i)*ds2(i,6)+c1(i,3)*two

      do 2015 j=1,6
      do 2015 i=1,n
 2015 c2(i,j)=(ds0(i,j)-ds1(i,j))/twozeta(i)

      do 2016 i=1,n
      gs0(i,01)=pxax(i)*fs0(i,01)-pxcx(i)*fs1(i,01)+c2(i,1)*three
      gs0(i,02)=pyay(i)*fs0(i,01)-pycy(i)*fs1(i,01)
      gs0(i,03)=pzaz(i)*fs0(i,01)-pzcz(i)*fs1(i,01)
      gs0(i,04)=pyay(i)*fs0(i,02)-pycy(i)*fs1(i,02)+c2(i,1)
      gs0(i,05)=pzaz(i)*fs0(i,02)-pzcz(i)*fs1(i,02)
      gs0(i,06)=pzaz(i)*fs0(i,03)-pzcz(i)*fs1(i,03)+c2(i,1)
      gs0(i,07)=pyay(i)*fs0(i,04)-pycy(i)*fs1(i,04)+c2(i,2)*two
      gs0(i,08)=pzaz(i)*fs0(i,04)-pzcz(i)*fs1(i,04)
      gs0(i,09)=pzaz(i)*fs0(i,05)-pzcz(i)*fs1(i,05)+c2(i,2)
      gs0(i,10)=pzaz(i)*fs0(i,06)-pzcz(i)*fs1(i,06)+c2(i,3)*two
      gs0(i,11)=pyay(i)*fs0(i,07)-pycy(i)*fs1(i,07)+c2(i,4)*three
      gs0(i,12)=pzaz(i)*fs0(i,07)-pzcz(i)*fs1(i,07)
      gs0(i,13)=pzaz(i)*fs0(i,08)-pzcz(i)*fs1(i,08)+c2(i,4)
      gs0(i,14)=pzaz(i)*fs0(i,09)-pzcz(i)*fs1(i,09)+c2(i,5)*two
 2016 gs0(i,15)=pzaz(i)*fs0(i,10)-pzcz(i)*fs1(i,10)+c2(i,6)*three

      do 2017 i=1,n
      dpcnd(i,01)=axbx(icontr(i))*ds0(i,1)+fs0(i,01)
      dpcnd(i,02)=ayby(icontr(i))*ds0(i,1)+fs0(i,02)
      dpcnd(i,03)=azbz(icontr(i))*ds0(i,1)+fs0(i,03)
      dpcnd(i,04)=axbx(icontr(i))*ds0(i,2)+fs0(i,02)
      dpcnd(i,05)=ayby(icontr(i))*ds0(i,2)+fs0(i,04)
      dpcnd(i,06)=azbz(icontr(i))*ds0(i,2)+fs0(i,05)
      dpcnd(i,07)=axbx(icontr(i))*ds0(i,3)+fs0(i,03)
      dpcnd(i,08)=ayby(icontr(i))*ds0(i,3)+fs0(i,05)
      dpcnd(i,09)=azbz(icontr(i))*ds0(i,3)+fs0(i,06)
      dpcnd(i,10)=axbx(icontr(i))*ds0(i,4)+fs0(i,04)
      dpcnd(i,11)=ayby(icontr(i))*ds0(i,4)+fs0(i,07)
      dpcnd(i,12)=azbz(icontr(i))*ds0(i,4)+fs0(i,08)
      dpcnd(i,13)=axbx(icontr(i))*ds0(i,5)+fs0(i,05)
      dpcnd(i,14)=ayby(icontr(i))*ds0(i,5)+fs0(i,08)
      dpcnd(i,15)=azbz(icontr(i))*ds0(i,5)+fs0(i,09)
      dpcnd(i,16)=axbx(icontr(i))*ds0(i,6)+fs0(i,06)
      dpcnd(i,17)=ayby(icontr(i))*ds0(i,6)+fs0(i,09)
 2017 dpcnd(i,18)=azbz(icontr(i))*ds0(i,6)+fs0(i,10)

      do 2018 i=1,n
      fpcnd(i,01)=axbx(icontr(i))*fs0(i,01)+gs0(i,01)
      fpcnd(i,02)=ayby(icontr(i))*fs0(i,01)+gs0(i,02)
      fpcnd(i,03)=azbz(icontr(i))*fs0(i,01)+gs0(i,03)
      fpcnd(i,04)=axbx(icontr(i))*fs0(i,02)+gs0(i,02)
      fpcnd(i,05)=ayby(icontr(i))*fs0(i,02)+gs0(i,04)
      fpcnd(i,06)=azbz(icontr(i))*fs0(i,02)+gs0(i,05)
      fpcnd(i,07)=axbx(icontr(i))*fs0(i,03)+gs0(i,03)
      fpcnd(i,08)=ayby(icontr(i))*fs0(i,03)+gs0(i,05)
      fpcnd(i,09)=azbz(icontr(i))*fs0(i,03)+gs0(i,06)
      fpcnd(i,10)=axbx(icontr(i))*fs0(i,04)+gs0(i,04)
      fpcnd(i,11)=ayby(icontr(i))*fs0(i,04)+gs0(i,07)
      fpcnd(i,12)=azbz(icontr(i))*fs0(i,04)+gs0(i,08)
      fpcnd(i,13)=axbx(icontr(i))*fs0(i,05)+gs0(i,05)
      fpcnd(i,14)=ayby(icontr(i))*fs0(i,05)+gs0(i,08)
      fpcnd(i,15)=azbz(icontr(i))*fs0(i,05)+gs0(i,09)
      fpcnd(i,16)=axbx(icontr(i))*fs0(i,06)+gs0(i,06)
      fpcnd(i,17)=ayby(icontr(i))*fs0(i,06)+gs0(i,09)
      fpcnd(i,18)=azbz(icontr(i))*fs0(i,06)+gs0(i,10)
      fpcnd(i,19)=axbx(icontr(i))*fs0(i,07)+gs0(i,07)
      fpcnd(i,20)=ayby(icontr(i))*fs0(i,07)+gs0(i,11)
      fpcnd(i,21)=azbz(icontr(i))*fs0(i,07)+gs0(i,12)
      fpcnd(i,22)=axbx(icontr(i))*fs0(i,08)+gs0(i,08)
      fpcnd(i,23)=ayby(icontr(i))*fs0(i,08)+gs0(i,12)
      fpcnd(i,24)=azbz(icontr(i))*fs0(i,08)+gs0(i,13)
      fpcnd(i,25)=axbx(icontr(i))*fs0(i,09)+gs0(i,09)
      fpcnd(i,26)=ayby(icontr(i))*fs0(i,09)+gs0(i,13)
      fpcnd(i,27)=azbz(icontr(i))*fs0(i,09)+gs0(i,14)
      fpcnd(i,28)=axbx(icontr(i))*fs0(i,10)+gs0(i,10)
      fpcnd(i,29)=ayby(icontr(i))*fs0(i,10)+gs0(i,14)
 2018 fpcnd(i,30)=azbz(icontr(i))*fs0(i,10)+gs0(i,15)

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

 2001 continue

      do 2020 i=1,n
      core(inddyy(icontr(i),01))=core(inddyy(icontr(i),01))-ddcond(i,01)
      core(inddyy(icontr(i),07))=core(inddyy(icontr(i),07))-ddcond(i,07)
      core(inddyy(icontr(i),08))=core(inddyy(icontr(i),08))-ddcond(i,08)
      core(inddyy(icontr(i),13))=core(inddyy(icontr(i),13))-ddcond(i,13)
      core(inddyy(icontr(i),14))=core(inddyy(icontr(i),14))-ddcond(i,14)
      core(inddyy(icontr(i),15))=core(inddyy(icontr(i),15))-ddcond(i,15)
      core(inddyy(icontr(i),19))=core(inddyy(icontr(i),19))-ddcond(i,19)
      core(inddyy(icontr(i),20))=core(inddyy(icontr(i),20))-ddcond(i,20)
      core(inddyy(icontr(i),21))=core(inddyy(icontr(i),21))-ddcond(i,21)
      core(inddyy(icontr(i),22))=core(inddyy(icontr(i),22))-ddcond(i,22)
      core(inddyy(icontr(i),25))=core(inddyy(icontr(i),25))-ddcond(i,25)
      core(inddyy(icontr(i),26))=core(inddyy(icontr(i),26))-ddcond(i,26)
      core(inddyy(icontr(i),27))=core(inddyy(icontr(i),27))-ddcond(i,27)
      core(inddyy(icontr(i),28))=core(inddyy(icontr(i),28))-ddcond(i,28)
      core(inddyy(icontr(i),29))=core(inddyy(icontr(i),29))-ddcond(i,29)
      core(inddyy(icontr(i),31))=core(inddyy(icontr(i),31))-ddcond(i,31)
      core(inddyy(icontr(i),32))=core(inddyy(icontr(i),32))-ddcond(i,32)
      core(inddyy(icontr(i),33))=core(inddyy(icontr(i),33))-ddcond(i,33)
      core(inddyy(icontr(i),34))=core(inddyy(icontr(i),34))-ddcond(i,34)
      core(inddyy(icontr(i),35))=core(inddyy(icontr(i),35))-ddcond(i,35)
 2020 core(inddyy(icontr(i),36))=core(inddyy(icontr(i),36))-ddcond(i,36)

      do 2021 i=1,n
      if(iminj(icontr(i)).ne.0) then

              core(inddyy(icontr(i),02))=core(inddyy(icontr(i),02))
     &                                  -ddcond(i,02)
              core(inddyy(icontr(i),03))=core(inddyy(icontr(i),03))
     &                                  -ddcond(i,03)
              core(inddyy(icontr(i),04))=core(inddyy(icontr(i),04))
     &                                  -ddcond(i,04)
              core(inddyy(icontr(i),05))=core(inddyy(icontr(i),05))
     &                                  -ddcond(i,05)
              core(inddyy(icontr(i),06))=core(inddyy(icontr(i),06))
     &                                  -ddcond(i,06)
              core(inddyy(icontr(i),09))=core(inddyy(icontr(i),09))
     &                                  -ddcond(i,09)
              core(inddyy(icontr(i),10))=core(inddyy(icontr(i),10))
     &                                  -ddcond(i,10)
              core(inddyy(icontr(i),11))=core(inddyy(icontr(i),11))
     &                                  -ddcond(i,11)
              core(inddyy(icontr(i),12))=core(inddyy(icontr(i),12))
     &                                  -ddcond(i,12)
              core(inddyy(icontr(i),16))=core(inddyy(icontr(i),16))
     &                                  -ddcond(i,16)
              core(inddyy(icontr(i),17))=core(inddyy(icontr(i),17))
     &                                  -ddcond(i,17)
              core(inddyy(icontr(i),18))=core(inddyy(icontr(i),18))
     &                                  -ddcond(i,18)
              core(inddyy(icontr(i),23))=core(inddyy(icontr(i),23))
     &                                  -ddcond(i,23)
              core(inddyy(icontr(i),24))=core(inddyy(icontr(i),24))
     &                                  -ddcond(i,24)
              core(inddyy(icontr(i),30))=core(inddyy(icontr(i),30))
     &                                  -ddcond(i,30)

                        endif
 2021 continue

      return
      end
