********************************************************************************
*                                                                              * 
*  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 ### relcoul1 ###
c
      subroutine relcoul1 (iorb,cw_sor,psi,pot,excp,bpot,d,f3,g,
     &			 lhs,rhs,wk2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      integer*4 cw_sor
      real*8 lhs
      include 'commons8.inc'
      dimension psi(*),pot(*),excp(*),cw_sor(*),
     &     rhs(*),bpot(*),d(*),f3(*),
     &     g(*),lhs(*),wk2(*)

c     dimension of wk2 array need be only (nni+8)*(nmi+4) but
c     (nni+8)*(nmi+8) array is reserved

      if (nel.eq.1) return

c     prepare right-hand side of Poisson's equation

      iborb=i1b(iorb)
      ngrid=i1si(iorb)	
      call prod2 (ngrid,psi(iborb),psi(iborb),wk2)
      call prod2 (ngrid,wk2	,g	,rhs)

c     even symmetry for coulomb potential
      isym=1
      ibpot=i2b(iorb)
      
      do itr1=1,maxsor1	  	
         do ig=1,i1ng(iorb)
            omega=ovfcoul(ig)
            omega1=1.d0-omega
            ioffst=ioffs(ig)
            ioffs1=ioffst+1	  
            ibpot1=ibpot+ioffst	  	  

            do i=1,4
               dmu2t(i)=dmu2(i,ig)
               dmu1t(i)=dmu1(i,ig)
            enddo
	  
c           prepare left-hand side of the poisson equation
c           include the diagonal part of the differentiation 
c           operator in lhs

            do i=1,ngsize(ig)
               lhs(ioffst+i)=f3(ioffst+i)+diag(ig)
            enddo	    
            
            if (i1ng(iorb).eq.1) then 
c              icase=1
               ifill=1
               ngrd1 =ingr1(1,ig)
               ngrd6a=ingr1(2,ig)
               ngrd6b=ingr1(3,ig)
               ngrd7 =ingr1(4,ig)	  	  
              
               call fill (nni,nmu(ig),isym,pot(iborb),wk2)
               do itr2=1,maxsor2/lsor	  		
                  call sor (wk2,lhs(ioffs1),rhs(ioffs1),
     &                 bpot(ioffs1),d(ioffs1),
     &                 cw_sor(iadext(ig)),cw_sor(iadnor(ig)),
     &                 cw_sor(iadex1(ig)),cw_sor(iadex2(ig)),
     &                 cw_sor(iadex3(ig)))
               enddo     
               call refill (nni,nmu(ig),pot(iborb),wk2)     
            else
               if (ig.eq.1) then 
c                 icase=2
                  ifill=2
                  ngrd1 =ingr2(1,ig)
                  ngrd6a=ingr2(2,ig)
                  ngrd6b=ingr2(3,ig)
                  ngrd7 =ingr2(4,ig)	  	  	    
                  
                  call fill2 (nni,nmu(ig),pot(iborb),wk2)

                  do itr2=1,maxsor2/lsor	  		
                     call sor 
     &                    (wk2,lhs(ioffs1),rhs(ioffs1),bpot(ioffs1),
     &                    d(ioffs1),
     &                    cw_sor(iadext(ig+ngrids)),
     &                    cw_sor(iadnor(ig+ngrids)),
     &                    cw_sor(iadex1(ig+ngrids)),
     &                    cw_sor(iadex2(ig+ngrids)),
     &                    cw_sor(iadex3(ig+ngrids)))
                  enddo     
                  call refill (nni,nmu(ig),pot(iborb),wk2)     
               elseif (ig.ne.1.and.ig.ne.i1ng(iorb)) then 
c                 icase=2
                  ifill=3
                  ngrd1 =ingr2(1,ig)
                  ngrd6a=ingr2(2,ig)
                  ngrd6b=ingr2(3,ig)
                  ngrd7 =ingr2(4,ig)	  	  	    
                  
                  muoffs=iemu(ig-1)-1
                  call fill3 (nni,nmu(ig),pot(iborb),wk2)
                  do itr2=1,maxsor2/lsor	  		
                     call sor 
     &                    (wk2,lhs(ioffs1),rhs(ioffs1),bpot(ioffs1),
     &                    d(ioffs1),
     &                    cw_sor(iadext(ig+ngrids)),
     &                    cw_sor(iadnor(ig+ngrids)),
     &                    cw_sor(iadex1(ig+ngrids)),
     &                    cw_sor(iadex2(ig+ngrids)),
     &                    cw_sor(iadex3(ig+ngrids)))
                  enddo     
                  call refill34 (nni,nmu(ig),pot(iborb),wk2)     
                  muoffs=0
               elseif (ig.eq.i1ng(iorb)) then 
                  ifill=4
                  ngrd1 =ingr1(1,ig)
                  ngrd6a=ingr1(2,ig)
                  ngrd6b=ingr1(3,ig)
                  ngrd7 =ingr1(4,ig)	  	  	    
                 
                  muoffs=iemu(ig-1)-1
                  
                  call fill4 (nni,nmu(ig),pot(iborb),wk2)
c                 call dcopy((nni+8)*(nmu(ig)+8),wk2,1,wkref2,1)
                  do itr2=1,maxsor2/lsor	  		
                     call sor 
     &                    (wk2,lhs(ioffs1),rhs(ioffs1),bpot(ioffs1),
     &                    d(ioffs1),
     &                    cw_sor(iadext(ig)),cw_sor(iadnor(ig)),
     &                    cw_sor(iadex1(ig)),cw_sor(iadex2(ig)),
     &                    cw_sor(iadex3(ig)))
                  enddo
                  
                  call refill34 (nni,nmu(ig),pot(iborb),wk2)     
                  muoffs=0	    
               endif
            endif
         enddo
      enddo
01021 format(/,i5,10d12.3,/(5x,10d12.3))
      return
      end