********************************************************************************
*                                                                              * 
*  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 ### sor ###
c
c     Performs one iteration of the sor process applied to
c     the equation (diagonal part of the differentiation operator, i.e.
c     the coefficient multiplying fun(ni,mi) is included in lhs array)
c
c    { d2(mi) + b(mi)*d1(mi) + d2(ni) + d(ni)*d1(ni)
c		      + lhs(ni,mi) } fun(ni,mi) = rhs(ni,mi)
c
c     warning! this routine uses mesh defined by mesh
c
      subroutine sor (fun,lhs,rhs,b,d,
     &     indx1,indx2,indx6a,indx6b,indx7)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      real*8	lhs
      include 'commons8.inc'
      
      dimension fun(*),lhs(*),rhs(*),b(*),d(*),
     &     indx1(*),indx2(*),indx6a(*),indx6b(*),indx7(*)	

c     indx1 (indx2) array elements point to the particular (consequtive)
c     element of fun (lhs, rhs, b and d) to be relaxed

c     first relax the inner points- region i (formula a1)
      
c     ipoiss=1 - approximately 20 mflops (Cray XMP)

      do i=1,ngrd1
         ij=indx1(i)
         ik=indx2(i)

         ddmi2=dmu2t(1) * ( fun(ij-nni4) + fun(ij+nni4) )+
     &        dmu2t(2) * ( fun(ij-nni3) + fun(ij+nni3) )+
     &        dmu2t(3) * ( fun(ij-nni2) + fun(ij+nni2) )+
     &        dmu2t(4) * ( fun(ij-nni1) + fun(ij+nni1) )
        
         ddmi1=dmu1t(1) * ( fun(ij-nni4) - fun(ij+nni4) )+
     &        dmu1t(2) * ( fun(ij-nni3) - fun(ij+nni3) )+
     &        dmu1t(3) * ( fun(ij-nni2) - fun(ij+nni2) )+
     &        dmu1t(4) * ( fun(ij-nni1) - fun(ij+nni1) )

         ddni2=dni2(1) * ( fun(ij-   4) + fun(ij+   4) )+
     &        dni2(2) * ( fun(ij-   3) + fun(ij+   3) )+
     &        dni2(3) * ( fun(ij-   2) + fun(ij+   2) )+
     &        dni2(4) * ( fun(ij-   1) + fun(ij+   1) )

         ddni1=dni1(1) * ( fun(ij-   4) - fun(ij+   4) )+
     &        dni1(2) * ( fun(ij-   3) - fun(ij+   3) )+
     &        dni1(3) * ( fun(ij-   2) - fun(ij+   2) )+
     &        dni1(4) * ( fun(ij-   1) - fun(ij+   1) )

         cc=ddmi2 + b(ik)*ddmi1 + ddni2 + d(ik)*ddni1
         fun(ij) = omega * (rhs(ik)-cc)/lhs(ik)+ omega1 * fun(ij)
      enddo

c     determine values at the bondary points from extrapolation
      
      if (isym.eq.1) then
         if (ifill.eq.1.or.ifill.eq.2) then
            do i=1,ngrd7
               ij=indx7(i)
               fun(ij)=exeven(1)*fun(ij+nni1)+exeven(2)*fun(ij+nni2)+
     &              exeven(3)*fun(ij+nni3)+exeven(4)*fun(ij+nni4)+
     &              exeven(5)*fun(ij+nni5)
            enddo
         endif     

         do i=1,ngrd6a
            ij=indx6a(i)
            fun(ij)=exeven(1)*fun(ij+1)+exeven(2)*fun(ij+2)+
     &           exeven(3)*fun(ij+3)+exeven(4)*fun(ij+4)+
     &           exeven(5)*fun(ij+5)
         enddo

         do i=1,ngrd6b
            ij=indx6b(i)
            fun(ij)=exeven(1)*fun(ij-1)+exeven(2)*fun(ij-2)+
     &           exeven(3)*fun(ij-3)+exeven(4)*fun(ij-4)+
     &           exeven(5)*fun(ij-5)
         enddo
      else
         if (ifill.eq.1.or.ifill.eq.2) then
            do i=1,ngrd7
               ij=indx7(i)
               fun(ij)=0.d0
            enddo
         endif
        
         do i=1,ngrd6a
            ij=indx6a(i)
            fun(ij)=0.d0
         enddo

         do i=1,ngrd6b
            ij=indx6b(i)
            fun(ij)=0.d0
         enddo
      endif

      return
      end