********************************************************************************
*                                                                              *
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### prepscf ####
c
c
      subroutine prepscf (cw_sor,cw_orb,cw_coul,cw_exch,cw_suppl,
     &					        cw_sctch)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      integer*4 cw_sor	
      include 'commons8.inc'
      dimension cw_orb(*),cw_coul(*),cw_exch(*),
     &          cw_suppl(*),cw_sor(*),cw_sctch(*)

c     orbital energies and Lagrange multiplies are initialized
c     or have been retrieved from the disk file

c     asymptotic values of coulomb and exchange potentials are 
c     calculated or have been retrieved fro the disk file

c     see routine inifun for the explanation of the following command

      if (imethod.eq.3) islat=1
c     if (islat.eq.0.and.alphaf.ne.0.d0) islat=1

      if (imethod.eq.2.or.ini.eq.4) nel=1

      if (idbg(495).ne.0) then 
         idump=1
      endif
      if (idump.ne.0) then

c        check orhogonalization and return

         if (idbg(401).ne.0) then
            iprt=2
            do iorb=1,norb
               jorb=norb+1-iorb
               call ortho (jorb,cw_orb,
     &              cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &              cw_sctch(i5b(1)),iprt)
            enddo
         endif

         iprt=0

         if (idbg(400).ne.0) iprt=1

         do iorb=1,norb
	    jorb=norb+1-iorb
       	    call norm (jorb,cw_orb,
     &           cw_suppl(i4b(9)),cw_suppl(i4b(14)),cw_sctch(i5b(1)))
	    call ortho (jorb,cw_orb,cw_suppl(i4b(9)),
     &           cw_suppl(i4b(14)),cw_sctch(i5b(1)),iprt)
         enddo

         do i=1,norb
	    itouch(i)=1
         enddo	  

         if (imethod.ne.2) then 
            call momen(cw_orb,cw_suppl,cw_sctch)

c           determine asymptotic values of coulomb and exchange potentials
c           from multipole expansion

            if (iform.eq.1.or.iform.eq.3) then
               if (iinterp.eq.0) then
                  call asympot (cw_coul,cw_exch)  
c                 write(*,*) '... asympt has been called ...'

               else
                  call asympti (cw_coul,cw_exch)
c                 write(*,*) '... asympti has been called ...'
               endif
            endif
         endif

c  calculate off- and diagonal Lagrange multipliers

         do iorb=norb,1,-1
            if (iform.eq.0.or.iform.eq.2) call rfdexch(iorb,cw_exch)
            call rayl (iorb,cw_orb,cw_coul,cw_exch,
     &           cw_suppl(i4b( 4)),cw_suppl(i4b( 5)),
     &           cw_suppl(i4b(13)),cw_suppl(i4b(14)),
     &           cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),
     &           cw_sctch(i5b( 3)),cw_sctch(i5b( 4)))
            call lagra (iorb,cw_orb,cw_coul,cw_exch,
     &           cw_suppl(i4b(14)),cw_sctch(i5b(1)),cw_sctch(i5b(2)))
            engi(iorb)=eng(iorb)
         enddo
         write(*,*) '... initializing Lagrange multipliers ...'
      endif  
c  ----- if (idump.ne.0) -----
      write(iout6,*)
      call toten (cw_orb,cw_coul,cw_exch,
     &     cw_suppl(i4b( 4)),cw_suppl(i4b( 5)),cw_suppl(i4b(13)),
     &     cw_suppl(i4b(14)),
     &     cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)),
     &     cw_sctch(i5b( 4)))
      write(iout6,6110)  etot
      write(iout6,6100)  evt
      if (islat.eq.0) write(iout6,6120)  virrat

06100 format(1x,'total electronic energy: ',d25.13)
06110 format(1x,'total energy:            ',d25.13)
06120 format(1x,'virial ratio:            ',d25.13)

      return
      end