C  Electrostatic lens calculation, based on chapter 7.
C  Uses block SOR with red/black orderings.
C  Auxilliary routines in lensaux.f
C  SOR function in lensor.f
C
C  File:  lens.f          Included file:  lens.h
C
C  Domain is represented by two groups of arrays
C  corresponding to:
C
C           U linear
C      Nx+------------+
C        |            |
C        |            |       U = potential
C        |            |
C     U=0|     B      |U=1000
C   ^    |            |
C   |    |            |      U=1000
C  k|  N | - - - - -  +-----------------+
C   |    |            |                 |
C        |                              |
C     U=0|              A               |U=1000
C        |                              |
C        |            |                 |
C      0 +------------------------------+
C        0    Axis    M     Axis        Mx
C                   i ---->
C
C  Parameters and global entities are defined in
C  the file  lens.h  using common blocks
C
C  /LENSFE/ - scalars       /LENSPE/ - arrays
C
C  Ax, Bx, Dx, Rx - Parameters for array dim's
C
C   M, Mx,  N, NX - Problem geometry
C  IA, KA, IB, KB - Problem coordinate arrays
C          UA, UB - The potential
C     Ared, Bred  - Mark red interior entries
C     Ablk, Bblk  - Mark black interior entries
C     Xred, Xblk  - Mark red/black on axis
C            W    - SOR parameter

C  Manage tests and examination of result.

      program testlens
      INCLUDE 'lens.h'

        logical checklens
        external checklens

      INTERFACE
        real function SORlens(it)
        integer it
      END INTERFACE

      character*50 buff     ! For terminal input
      integer iter, titer
      real Up, eps, err
      integer eltm    ! For elapsed time on mpp
      logical more

      more = .TRUE.
      do while ( more )
        print  '(" M, Mx, N, Nx > " $)'
        read *, M, Mx, N, Nx
        more = .not. checklens()                
      end do
      call setuplens()
      print '(" potential, iterations, epsilon, W > "$)'
      read *, Up, iter, eps, W
      call initlens(Up)
      err = 1.0 + 2.0*eps   ! To assure iterations
      titer = 0
       do while ( err > eps )
        call mpTimerStart()
        err = SORlens(iter)
        eltm = mpTimerElapsed()
        titer = titer + iter
        print 21, titer, err, eltm/1000.0
 21     FORMAT (1X, I4, " iterations,  Err = ", E13.6,
     *          "  Time = ", F8.3, " sec." )
      end do

C  Browse results until user decides to quit:
      more = .TRUE.
      print *, "Specify ranges in problem coordinates."
      do while ( more )
        print '(" Browse? ( A | B | N | n ) > " $)'
        read *, buff
        if ((buff(1:1) == 'N') .or. (buff(1:1) == 'n'))
     *      more = .FALSE.
        if ( (buff(1:1) == 'A')
     *      .or. (buff(1:1) == 'B')) then
          print '( " i1, i2, k1, k2 > " $)'
          read *, i1, i2, k1, k2
        end if
        if ( buff(1:1) == 'A' )
     *    call prA(i1, i2, k1, k2)
        if ( buff(1:1) == 'B' )
     *     call prb(i1, i2, k1, k2)
      end do
      end program testlens
