c---------------------------------------------
c  Relaxation routine.
c
c  Wed Feb 14 12:55:03 PST 2001
c  
c  For the Black String Coalition.
c---------------------------------------------

      subroutine relax(alpha, a, c, br, bz, krr,
     &                 phi, pi, phi_r, phi_z,
     &                 Nr, Nz, r, z,
     &                 order, tol, maxiter, 
     &                 mg_rhs1,mg_rhs2, mg_rhs3,
     &                 mg_rhs4)
      
      implicit      none 

      logical       ltrace
      parameter   ( ltrace = .true. )

      logical       relax_trace
      parameter   ( relax_trace = .true. )


      character*5   cdnm
      parameter   ( cdnm = 'relax' )

      integer       Nr,   Nz    
      integer       i,    j    

      real*8        dr,     dz 

      real*8        alpha(Nr, Nz)
      real*8        a(Nr, Nz),     c(Nr,Nz)
      real*8        br(Nr, Nz),    bz(Nr, Nz)

      real*8        krr(Nr, Nz)
      real*8        phi(Nr, Nz),   pi(Nr, Nz)
      real*8        phi_r(Nr, Nz), phi_z(Nr, Nz)
      
      real*8        r(Nr),         z(Nz)         

      real*8        alpha1d(Nr)
      real*8        a1d(Nr)
      real*8        br1d(Nr)


c---------------------------------------------
c   tol      ---> tolerance. 
c   maxiter  ---> maximum number of iterations. 
c---------------------------------------------

      real*8        tol
      integer       maxiter,     iter    

c---------------------------------------------
c    order ---> 0 for Lexicographic.
c               1 for Red-Black.
c---------------------------------------------

      integer       order

      real*8        norm_res           

c---------------------------------------------
c   Functions.
c---------------------------------------------
      
      real*8        res_relax

      integer       shape(2)
      real*8        bbox(4)
      integer       rank

      integer       out 
      real*8        time

c---------------------------------------------
c  Multigrid rhs   
c   mg_rhs1 for the slicing condition.  
c   mg_rhs1 for the HC.  
c   mg_rhs1 for the r-MC.  
c   mg_rhs1 for the z-MC.  
c---------------------------------------------
      
      real*8         mg_rhs1(Nr,Nz)
      real*8         mg_rhs2(Nr,Nz)
      real*8         mg_rhs3(Nr,Nz)
      real*8         mg_rhs4(Nr,Nz)

c=============================================
c=============================================
      dr = r(2) - r(1)
      dz = z(2) - z(1)
      if ( ltrace ) then
         write(0,*) cdnm,':                  Nr =', Nr 
         write(0,*) cdnm,':                  Nz =', Nz 
         write(0,*) cdnm,':                  dr =', dr    
         write(0,*) cdnm,':                  dz =', dz    
         write(0,*) cdnm,':  order(0=lexi/1=R-B)=', order
         write(0,*) cdnm,':            maxiter  =', maxiter   
         write(0,*) cdnm,':            tol      =', tol   
      endif

      if ( order .gt. 1 )  goto 9999 
      if ( order .lt. 0 )  goto 9999 
      if ( Nr .le. 0 )     goto 9999 
      if ( Nz .le. 0 )     goto 9999 
      if ( dr .le. 0.0d0 ) goto 9999 
      if ( dz .le. 0.0d0 ) goto 9999 
      shape(1) = Nr   
      shape(2) = Nz   
      rank = 2  
      bbox(1) = r(1)
      bbox(2) = r(Nr)
      bbox(3) = z(1)
      bbox(4) = z(Nz)
      out = 1
c-----------------------------------------------------------
c     R E L A X A T I O N    L O O P
c-----------------------------------------------------------
      time = 1.0d0 
      do iter = 1, maxiter 
         norm_res = 0.0d0
         norm_res = res_relax(alpha, a, c, br, bz, krr,
     &                 phi, pi, phi_r, phi_z,
     &                 Nr, Nz, r, z, order, mg_rhs1,
     &                 mg_rhs2, mg_rhs3, mg_rhs4)
         if ( tol .ne. 0.0d0  ) then
              write(0,*) cdnm,':  Nr*Nz=',Nr*Nz,', iter=',iter,
     &                        ', norm_res=',norm_res
         endif
         out = out + 1 
         if ( out .eq. 10 ) out = 0 
         do i = 1 , Nr
             a1d(i) = a(i,1)
             alpha1d(i) = alpha(i,1)
             br1d(i) = br(i,1)
         enddo 
         if ( relax_trace .and. out .eq. 1) then
            call gft_out_bbox('alpha', time, shape, rank, bbox, alpha)
            call gft_out_bbox('a', time, shape, rank, bbox, a)
            call gft_out_bbox('br', time, shape, rank, bbox, br)
            call gft_out_bbox('bz', time, shape, rank, bbox, bz)
            call vsxynt('a1d',time,r,a1d,nr) 
            call vsxynt('alpha1d',time,r,alpha1d,nr) 
            call vsxynt('br1d',time,r,br1d,nr) 
         endif
         time = time + 1.0d0
         if ( norm_res .le. tol ) goto 1000
      enddo

      if ( tol .ne. 0.0d0 ) then
         write(0,*) cdnm,':  iter=maxiter' 
      endif 

1000  continue
      return
9999  continue 
      write(0,*) cdnm,':********ATTENTION*****************'
      write(0,*) cdnm, ' failed !!!'
      write(0,*)
      write(0,*) cdnm,':                  Nr =', Nr 
      write(0,*) cdnm,':                  Nz =', Nz 
      write(0,*) cdnm,':                  dr =', dr    
      write(0,*) cdnm,':                  dz =', dz    
      write(0,*) cdnm,':  order(0=lexi/1=R-B)=', order
      write(0,*)
      write(0,*) cdnm,':********ATTENTION*****************'
      return
      end