********************************************************************************
*                                                                              *
*  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 ### momen ###
c
c     Recalculates multipole moment expansion coefficients every time demax(1),
c     i.e. maximum error in orbital energy, is reduced by facmul

c     Warning!
c     Coefficients and then asymptotic values are recalculated only for
c     orbitals which underwent relaxation, i.e. those being touched (itouch=1)

      subroutine momen (cw_orb,cw_suppl,cw_sctch)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension cw_orb(*),cw_suppl(*),cw_sctch(*)

      if (nel.eq.1) return

      call ttime(time1)
      call coulmom (cw_orb,cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &     cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)),
     &     cw_sctch(i5b( 4)),cw_sctch(i5b( 5)),cw_sctch(i5b( 6)),     
     &     cw_sctch(i5b( 7)),
     &     cw_sctch(i5b( 8)),cw_sctch(i5b( 9)),cw_sctch(i5b(10)))

      if (nel.gt.1.and.imethod.eq.1) then
         do  iorb1=1,norb
            do  iorb2=iorb1,norb
               if (idbg(300).eq.0) then
                  call exchmom1 (iorb1,iorb2,cw_orb,
     &                 cw_suppl(i4b(9)),cw_suppl(i4b(14)),	     
     &                 cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),
     &                 cw_sctch(i5b( 3)),cw_sctch(i5b( 4)),
     &                 cw_sctch(i5b( 5)),cw_sctch(i5b( 6)),     
     &                 cw_sctch(i5b( 7)),cw_sctch(i5b( 8)),
     &                 cw_sctch(i5b( 9)),cw_sctch(i5b(10)))
               else
                  call exchmom (iorb1,iorb2,cw_orb,
     &                 cw_suppl(i4b(9)),cw_suppl(i4b(14)),	     
     &                 cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),
     &                 cw_sctch(i5b( 3)),cw_sctch(i5b( 4)),
     &                 cw_sctch(i5b( 5)),cw_sctch(i5b( 6)),     
     &                 cw_sctch(i5b( 7)),cw_sctch(i5b( 8)),
     &                 cw_sctch(i5b( 9)),cw_sctch(i5b(10)))

               endif
            enddo
         enddo
      endif
      call ttime (time2)
      tmomen =tmomen + (time2-time1)
      if (iprtlev.ne.3) then 
         write(iout6,*) 
     &        '... multipole moment expansion coefficients',
     &        ' (re)calculated ...'
      endif
      return
      end