       subroutine lookup ( ra, rmin, rmax, ta, tmin, tmax, j1, j2, isen
     1                   , qty, irb, jtb )
c
c    dac:zeus3d.lookup <----------- generates lookup tables for RT -> XY
c                                                            april, 1992
c
c    written by: David Clarke
c    modified 1:
c
c  PURPOSE:  For the purposes of speeding up the RT to XY re-map, this
c  subroutine generates high-resolution look-up vectors ("irb", "jtb")
c  which are used to locate a set of arbitrary coordinates (r,theta) on
c  a discretised grid (ra(i),ta(j)).
c
c  INPUT VARIABLES:
c    ra         face-centred radial coordinate of computational grid
c    ta         face-centred angular coordinate of computational grid
c    j1,j2      starting and ending indices for "ta".
c    rmax,rmin  desired extrema of "ra"
c    tmax,tmin  desired extrema of "ta"
c    isen       >0 => map only the portion of "ta" between "tmin" and
c                     "tmax" to the look-up vector "jtb".
c               <0 => map all of "ta" to the look-up vector "jtb".
c
c  OUTPUT VARIABLES:
c    irb,jtb    look-up vectors to facilitate remapping
c    qty(1,2)   effective range of "ra"  spanned by "irb"
c    qty(3)     resolution of look-up vector "irb"
c    qty(4,5)   effective range of "ta" spanned by "jtb"
c    qty(6)     resolution of look-up vector "jtb"
c
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c------------------------ IMPLICIT STATEMENTS --------------------------
c-----------------------------------------------------------------------
c
       implicit      none
c
c-----------------------------------------------------------------------
c----------------------------- PARAMETERS ------------------------------
c-----------------------------------------------------------------------
c
c      Primary parameters:
c
       integer       in      , jn      , kn      , ijkn    , nxpx
     1             , nypx    , nxrd    , nyrd    , niov    , nios
     2             , ntsl    , mn      , nmd     , isig
       real*8        pi  , tiny  , huge  , one  , mone  , haf  , zro
c
       parameter     ( in = 5005, jn = 1, kn = 10, ijkn = 5005)
       parameter     ( nxpx = 1, nypx = 1, nxrd =   1, nyrd =   1 )
       parameter     ( niov =  20, nios =  20, ntsl =1000 )
       parameter     ( mn   =   1, nmd  =   1 )
       parameter     ( pi   = 3.1415926535898 )
       parameter     ( isig = 9, tiny = 1.0d-99, huge = 1.0d+99 )
       parameter     ( one  = 1.0d+00, mone = -1.0d+00
     1                               , haf  =  0.5d+00, zro = 0.0d+00)
c
c      Secondary parameters:
c
       integer        idim    , jdim
     1              , nt01    , nt02    , nt03    , nt04    , nt05
     2              , nt06    , nt07    , nt08    , nt09    , nt10
     3              , nt11    , nt12    , nt13    , nt14    , nt15
     4              , nt16    , nt17    , nt18    , nt19    , nt20
       parameter     ( idim =   kn, jdim =   in )
       parameter     ( nt01 = 600, nt02 = 600, nt03 =   7, nt04 =   5 )
       parameter     ( nt05 =  18, nt06 =   1, nt07 =  10, nt08 =  12 )
       parameter     ( nt09 =   5, nt10 =  54, nt11 =   8 )
       parameter     ( nt12 =        16, nt17 =         4 )
       parameter     ( nt13 = nt12 +  4, nt18 = nt17 +  1 )
       parameter     ( nt14 = nt13     , nt19 = nt18      )
       parameter     ( nt15 = nt14 + 30, nt20 = nt19 +  3 )
       parameter     ( nt16 = nt14 + 42 )
c
c-----------------------------------------------------------------------
c-------------------- SCRATCH SPACE (WORKER ARRAYS) --------------------
c-----------------------------------------------------------------------
c
c      1-D worker arrays
c
       real*8        wa1d    (ijkn), wb1d    (ijkn), wc1d    (ijkn)
     1             , wd1d    (ijkn), we1d    (ijkn), wf1d    (ijkn)
     2             , wg1d    (ijkn), wh1d    (ijkn), wi1d    (ijkn)
     3             , wj1d    (ijkn), wk1d    (ijkn), wl1d    (ijkn)
     4             , wm1d    (ijkn), wn1d    (ijkn), wo1d    (ijkn)
     5             , wp1d    (ijkn), wq1d    (ijkn), wr1d    (ijkn)
     6             , ws1d    (ijkn), wt1d    (ijkn), wu1d    (ijkn)
     7             , wv1d    (ijkn), ww1d    (ijkn), wx1d    (ijkn)
     8             , wy1d    (ijkn), wz1d    (ijkn)
c
c      2-D worker arrays
c
       real*8        wa2d    (idim,jdim), wa2dt   (jdim,idim)
     1             , wb2d    (idim,jdim), wb2dt   (jdim,idim)
     1             , wc2d    (idim,jdim), wc2dt   (jdim,idim)
     2             , wd2d    (idim,jdim), wd2dt   (jdim,idim)
     1             , we2d    (idim,jdim), we2dt   (jdim,idim)
     2             , wf2d    (idim,jdim), wf2dt   (jdim,idim)
c
c      3-D worker arrays
c
       real*8        wa3d    (  in,  jn,  kn), wb3d    (  in,  jn,  kn)
     1             , wc3d    (  in,  jn,  kn), wd3d    (  in,  jn,  kn)
     2             , we3d    (  in,  jn,  kn), wf3d    (  in,  jn,  kn)
     3             , wg3d    (  in,  jn,  kn)
c
       equivalence   ( wa2d    , wa2dt    ), ( wb2d    , wb2dt    )
     1             , ( wc2d    , wc2dt    ), ( wd2d    , wd2dt    )
     2             , ( we2d    , we2dt    ), ( wf2d    , wf2dt    )
       equivalence   ( wa3d    , wa2d     ), ( wb3d    , wb2d     )
     1             , ( wc3d    , wc2d     ), ( wd3d    , wd2d     )
     2             , ( we3d    , we2d     ), ( wf3d    , wf2d     )
c
c---------------------------- Common Blocks ----------------------------
c
       common / scratchr /
     1               wa1d    , wb1d    , wc1d    , wd1d    , we1d
     2             , wf1d    , wg1d    , wh1d    , wi1d    , wj1d
     3             , wk1d    , wl1d    , wm1d    , wn1d    , wo1d
     4             , wp1d    , wq1d    , wr1d    , ws1d    , wt1d
     5             , wu1d    , wv1d    , ww1d    , wx1d    , wy1d
     6             , wz1d
       common / scratchr /
     1               wa3d    , wb3d    , wc3d    , wd3d    , we3d
     2             , wf3d
     3             , wg3d
c
       integer       j1      , j2      , isen
       integer       jmax    , j       , jc      , nmax    , i
     1             , im1     , n       , jm1
       real*8        rmin    , rmax    , tmin    , tmax
       real*8        ramn    , ramx    , dra     , rb      , tamn
     1             , tamx    , dta     , tb      , small
c
       integer       irb     (   *), jtb     (   *)
       real*8        ra      (   *), ta      (   *), qty     (   *)
       real*8        tax   (3*ijkn)
c
c      Careful!  "wa1d" through "wf1d" are used in PIXDMP.
c
       equivalence   ( tax     , wg1d     )
c
c      Data statements
c
       data          small    / 0.001 /
c
c-----------------------------------------------------------------------
c
c      "tax" is the extended angular grid, and stores three branches of
c  the grid specified in "ta" in the same array.  This is only required
c  when the angular grid "ta" spans 2*pi radians, and when the segment
c  to be mapped to a Cartesian grid straddles the branch point.
c  Otherwise, "tax" can be thought of as nothing more than a copy of
c  "ta".
c
       jmax = j2 - j1 + 1
       do 10 j=1,jmax
         jc            = j + j1 - 1
         tax(j       ) = ta(jc) - 2.0 * pi
         tax(j+  jmax) = ta(jc)
         tax(j+2*jmax) = ta(jc) + 2.0 * pi
10     continue
c
c      "irb" is a (long) vector containing the indices of the nearest
c  value of "ra" for a large sample of radial coordinates possible on
c  the Cartesian grid.
c
       nmax = 20 * ijkn
       ramn = 0.0
       ramx = sqrt ( 2.0 ) * rmax
       dra  = ( ramx - ramn ) /  real(nmax)
       rb   = ramn - dra / 2.0
       i    = 1
       im1  = 1
       do 30 n=1,nmax
         rb = rb + dra
         if ( (rb .lt. rmin) .or. (rb .gt. rmax) ) then
           irb(n) = 0
           go to 30
         endif
20       continue
         if (rb .gt. ra(i)) then
           im1 = i
           i   = i + 1
           go to 20
         else
           irb(n) = im1
         endif
30     continue
c
c      "jtb" is a (long) vector containing the indices of the nearest
c  value of "ta" for a large sample of angular coordinates possible on
c  the Cartesian grid.
c
       nmax = 40 * ijkn
       tamn =  real(  int ( 2.0*tmin/pi + 4.0 + small ) ) * pi / 2.0
     1      - 2.0 * pi
       tamx =  real(  int ( 2.0*tmax/pi + 5.0 - small ) ) * pi / 2.0
     1      - 2.0 * pi
       if ( (tamx-tamn .gt. pi+small) .or. (isen .lt. 0) ) then
         tamn = tmax - 2.0 * pi
         tamx = tmax
       endif
       dta  = ( tamx - tamn ) /  real(nmax)
       tb   = tamn - dta / 2.0
       j    = 1
       jm1  = 1
       do 50 n=1,nmax
         tb = tb + dta
         if ( (tb .lt. tmin) .or. (tb .gt. tmax) ) then
           jtb(n) = 0
           go to 50
         endif
40       continue
         if (tb .gt. tax(j)) then
           jm1 = j
           j   = j + 1
           go to 40
         else
           jtb(n) = mod ( jm1, jmax ) + j1 - 1
         endif
50     continue
c
c      Set elements of "qty" to scalars.
c
       qty(1) = ramn
       qty(2) = ramx
       qty(3) = dra
       qty(4) = tamn
       qty(5) = tamx
       qty(6) = dta
c
       return
       end
