      subroutine carbo(ncontmax,ncont,nindpndnt,nalpha,nbeta,
     &                 elevelshift,fockalpha,fockbeta,coeffalpha,
     &                 coeffbeta,overlap,wkmat1,wkmat2)

      implicit real*8(a-h,o-z)

      dimension fockalpha(*),fockbeta(*),
     &          coeffalpha(ncontmax,*),coeffbeta(ncontmax,*),
     &          overlap(ncontmax,*)

      dimension wkmat1(ncont,*),wkmat2(ncont,*)

      data zero/0.0d0/

      rewind 45

      do 1001 i=1,ncont
 1001 call fastrd(45,overlap(1,i),ncont)

      do 1002 i=1,ncont
      do 1002 j=1,ncont
      if(j.le.nalpha) wkmat2(j,i)=-coeffalpha(i,j)*elevelshift
      if(j.gt.nalpha) wkmat2(j,i)=+coeffalpha(i,j)*elevelshift
 1002 if(j.gt.nindpndnt) wkmat2(j,i)=zero

      do 1003 j=1,ncont
      do 1004 i=1,ncont
 1004 wkmat1(i,j)=zero
      do 1003 k=1,nindpndnt
      do 1003 i=1,ncont
 1003 wkmat1(i,j)=wkmat1(i,j)+coeffalpha(i,k)*wkmat2(k,j)

      do 1005 j=1,ncont
      do 1006 i=1,ncont
 1006 wkmat2(i,j)=zero
      do 1005 k=1,ncont
      do 1005 i=1,ncont
 1005 wkmat2(i,j)=wkmat2(i,j)+wkmat1(i,k)*overlap(k,j)

      do 1007 j=1,ncont
      istart=j*(j-1)/2
      do 1007 k=1,ncont
      do 1007 i=1,j
 1007 fockalpha(istart+i)=fockalpha(istart+i)+overlap(i,k)*wkmat2(k,j)

      do 1008 i=1,ncont
      do 1008 j=1,ncont
      if(j.le.nbeta) wkmat2(j,i)=-coeffbeta(i,j)*elevelshift
      if(j.gt.nbeta) wkmat2(j,i)=+coeffbeta(i,j)*elevelshift
 1008 if(j.gt.nindpndnt) wkmat2(j,i)=zero

      do 1009 j=1,ncont
      do 1010 i=1,ncont
 1010 wkmat1(i,j)=zero
      do 1009 k=1,nindpndnt
      do 1009 i=1,ncont
 1009 wkmat1(i,j)=wkmat1(i,j)+coeffbeta(i,k)*wkmat2(k,j)

      do 1011 j=1,ncont
      do 1012 i=1,ncont
 1012 wkmat2(i,j)=zero
      do 1011 k=1,ncont
      do 1011 i=1,ncont
 1011 wkmat2(i,j)=wkmat2(i,j)+wkmat1(i,k)*overlap(k,j)

      do 1013 j=1,ncont
      istart=j*(j-1)/2
      do 1013 k=1,ncont
      do 1013 i=1,j
 1013 fockbeta(istart+i)=fockbeta(istart+i)+overlap(i,k)*wkmat2(k,j)

      return
      end
