c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine ssor(niter)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   to perform pseudo-time stepping SSOR iterations
c   for five nonlinear pde's.
c---------------------------------------------------------------------

      implicit none
      integer niter

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, k, m, n, lst, lend
      integer istep
      double precision  tmp, tmp2, tv(5*isiz1*isiz2)
      double precision  delunm(5)

      external timer_read
      double precision timer_read

 
c---------------------------------------------------------------------
c   begin pseudo-time stepping iterations
c---------------------------------------------------------------------
      tmp = 1.0d+00 / ( omega * ( 2.0d+00 - omega ) ) 

c---------------------------------------------------------------------
c   initialize a,b,c,d to zero (guarantees that page tables have been
c   formed, if applicable on given architecture, before timestepping).
c---------------------------------------------------------------------
!$omp parallel default(shared) private(m,n,i,j)
!$omp do
      do j=jst,jend
         do i=ist,iend
            do m=1,5
               do n=1,5
                  a(m,n,i,j) = 0.d0
                  b(m,n,i,j) = 0.d0
                  c(m,n,i,j) = 0.d0
                  d(m,n,i,j) = 0.d0
               enddo
            enddo
         enddo
      enddo
!$omp end do nowait
!$omp do
      do j=jend,jst,-1
         do i=iend,ist,-1
            do m=1,5
               do n=1,5
                  au(m,n,i,j) = 0.d0
                  bu(m,n,i,j) = 0.d0
                  cu(m,n,i,j) = 0.d0
                  du(m,n,i,j) = 0.d0
               enddo
            enddo
         enddo
      enddo
!$omp end do nowait
!$omp end parallel
      do i = 1, t_last
         call timer_clear(i)
      end do

c---------------------------------------------------------------------
c   compute the steady-state residuals
c---------------------------------------------------------------------
      call rhs
 
c---------------------------------------------------------------------
c   compute the L2 norms of newton iteration residuals
c---------------------------------------------------------------------
      call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
     >             ist, iend, jst, jend,
     >             rsd, rsdnm )

 
      do i = 1, t_last
         call timer_clear(i)
      end do
      call timer_start(1)
 
c---------------------------------------------------------------------
c   the timestep loop
c---------------------------------------------------------------------
      do istep = 1, niter

         if (mod ( istep, 20) .eq. 0 .or.
     >         istep .eq. itmax .or.
     >         istep .eq. 1) then
            if (niter .gt. 1) write( *, 200) istep
 200        format(' Time step ', i4)
         endif
 
c---------------------------------------------------------------------
c   perform SSOR iteration
c---------------------------------------------------------------------
!$omp parallel default(shared) private(i,j,k,m,tmp2,lst,lend)
!$omp&  shared(ist,iend,jst,jend,nx,ny,nz,nx0,ny0,omega)
!$omp master
         if (timeron) call timer_start(t_rhs)
!$omp end master
         tmp2 = dt
!$omp do
         do k = 2, nz - 1
            do j = jst, jend
               do i = ist, iend
                  do m = 1, 5
                     rsd(m,i,j,k) = tmp2 * rsd(m,i,j,k)
                  end do
               end do
            end do
         end do
!$omp end do nowait
!$omp master
         if (timeron) call timer_stop(t_rhs)
!$omp end master

         lst = ist + jst
         lend = iend + jend
!$omp barrier

         do k = 2, nz -1 
c---------------------------------------------------------------------
c   form the lower triangular part of the jacobian matrix
c---------------------------------------------------------------------
!$omp master
            if (timeron) call timer_start(t_jacld)
!$omp end master
            call jacld(k)
!$omp master
            if (timeron) call timer_stop(t_jacld)
 
c---------------------------------------------------------------------
c   perform the lower triangular solution
c---------------------------------------------------------------------
            if (timeron) call timer_start(t_blts)
!$omp end master
            call blts( isiz1, isiz2, isiz3,
     >                 nx, ny, nz, k,
     >                 omega,
     >                 rsd, 
     >                 a, b, c, d,
     >                 ist, iend, jst, jend, 
     >                 lst, lend )
!$omp master
            if (timeron) call timer_stop(t_blts)
!$omp end master
          end do


          do k = nz - 1, 2, -1
c---------------------------------------------------------------------
c   form the strictly upper triangular part of the jacobian matrix
c---------------------------------------------------------------------
!$omp master
            if (timeron) call timer_start(t_jacu)
!$omp end master
            call jacu(k)
!$omp master
            if (timeron) call timer_stop(t_jacu)

c---------------------------------------------------------------------
c   perform the upper triangular solution
c---------------------------------------------------------------------
            if (timeron) call timer_start(t_buts)
!$omp end master
            call buts( isiz1, isiz2, isiz3,
     >                 nx, ny, nz, k,
     >                 omega,
     >                 rsd, tv,
     >                 du, au, bu, cu,
     >                 ist, iend, jst, jend,
     >                 lst, lend )
!$omp master
            if (timeron) call timer_stop(t_buts)
!$omp end master
          end do

 
c---------------------------------------------------------------------
c   update the variables
c---------------------------------------------------------------------

!$omp master
         if (timeron) call timer_start(t_add)
!$omp end master
         tmp2 = tmp
!$omp do
         do k = 2, nz-1
            do j = jst, jend
               do i = ist, iend
                  do m = 1, 5
                     u( m, i, j, k ) = u( m, i, j, k )
     >                    + tmp2 * rsd( m, i, j, k )
                  end do
               end do
            end do
         end do
!$omp end do nowait
!$omp end parallel
         if (timeron) call timer_stop(t_add)
 
c---------------------------------------------------------------------
c   compute the max-norms of newton iteration corrections
c---------------------------------------------------------------------
         if ( mod ( istep, inorm ) .eq. 0 ) then
            if (timeron) call timer_start(t_l2norm)
            call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
     >                   ist, iend, jst, jend,
     >                   rsd, delunm )
            if (timeron) call timer_stop(t_l2norm)
c            if ( ipr .eq. 1 ) then
c                write (*,1006) ( delunm(m), m = 1, 5 )
c            else if ( ipr .eq. 2 ) then
c                write (*,'(i5,f15.6)') istep,delunm(5)
c            end if
         end if
 
c---------------------------------------------------------------------
c   compute the steady-state residuals
c---------------------------------------------------------------------
         call rhs
 
c---------------------------------------------------------------------
c   compute the max-norms of newton iteration residuals
c---------------------------------------------------------------------
         if ( ( mod ( istep, inorm ) .eq. 0 ) .or.
     >        ( istep .eq. itmax ) ) then
            if (timeron) call timer_start(t_l2norm)
            call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
     >                   ist, iend, jst, jend,
     >                   rsd, rsdnm )
            if (timeron) call timer_stop(t_l2norm)
c            if ( ipr .eq. 1 ) then
c                write (*,1007) ( rsdnm(m), m = 1, 5 )
c            end if
         end if

c---------------------------------------------------------------------
c   check the newton-iteration residuals against the tolerance levels
c---------------------------------------------------------------------
         if ( ( rsdnm(1) .lt. tolrsd(1) ) .and.
     >        ( rsdnm(2) .lt. tolrsd(2) ) .and.
     >        ( rsdnm(3) .lt. tolrsd(3) ) .and.
     >        ( rsdnm(4) .lt. tolrsd(4) ) .and.
     >        ( rsdnm(5) .lt. tolrsd(5) ) ) then
c            if (ipr .eq. 1 ) then
               write (*,1004) istep
c            end if
            go to 900
         end if
 
      end do
  900 continue
 
      call timer_stop(1)
      maxtime= timer_read(1)
 

      return
      
 1001 format (1x/5x,'pseudo-time SSOR iteration no.=',i4/)
 1004 format (1x/1x,'convergence was achieved after ',i4,
     >   ' pseudo-time steps' )
 1006 format (1x/1x,'RMS-norm of SSOR-iteration correction ',
     > 'for first pde  = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for second pde = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for third pde  = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for fourth pde = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for fifth pde  = ',1pe12.5)
 1007 format (1x/1x,'RMS-norm of steady-state residual for ',
     > 'first pde  = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'second pde = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'third pde  = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'fourth pde = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'fifth pde  = ',1pe12.5)
 
      end