       subroutine rttoxy ( iar, ar, rvecb, drvecbi, tvecb, dtvecbi
     1                   , j1, j2, nxz, nyz, xb, yb, qvar, irb, jtb
     2                   , th0, iqty, qty )
c
c    dac:zeus3d.rttoxy <------------------- Map 2-D RT array to XY array
c                                                         february, 1992
c
c    written by: David Clarke
c    modified 1:
c
c  PURPOSE:  Maps a (not necessarily uniformly gridded) spherical polar
c  2-D array "ar" onto a (not necessarily uniformly gridded) Cartesian
c  2-D array "qty".  The spherical polar image is rotated by an angle
c  "th0" before it is remapped.
c
c  INPUT VARIABLES:
c    iar         leading dimension of "ar"
c    ar          input array (a 2-D polar slice of a 3-D datacube)
c    rvecb       zone-centred radial coordinate of computational grid
c    drvecbi     inverse of radial b-grid differences
c    tvecb       zone-centred angular coordinate of computational grid
c    dtvecbi     inverse of angular b-grid differences
c    j1,j2       starting and ending indices for "tvecb"
c    nxz,nyz     number of Cartesian zones in x- and y-directions.
c    xb,yb       zone-centred Cartesian grid.
c    qvar(1,2)   effective range of radial  coordinates spanned by "irb"
c    qvar(3)     resolution of look-up vector "irb"
c    qvar(4,5)   effective range of angular coordinates spanned by "jtb"
c    qvar(6)     resolution of look-up vector "jtb"
c    irb,jtb     look-up vectors to facilitate remapping
c    th0         CC angle about centre to rotate image before remapping
c    iqty        leading dimension of "qty"
c
c  The quantities "qvar", "irb", and "jtb" are computed by LOOKUP.
c
c  OUTPUT VARIABLES:
c    qty         remapped array
c
c  EXTERNALS: [NONE]
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
       integer       iar     , j1      , j2      , nxz     , nyz
     1             , iqty    , j       , nx      , ny      , ic
     2             , icmp1   , jc      , jcmp1
       real*8        th0     , sinth0  , costh0  , rbmn    , drai
     1             , rb0     , drb     , tamn    , tamx    , tbmn
     2             , dtai    , tb0     , dtb     , sintb0  , xb0
     3             , yb0     , sgnxb   , sgnyb   , sgnmn   , sgnmx
     4             , q1      , q2      , q3      , darr    , darl
     5             , d2ar    , dqx1    , dqx2    , qmax    , qmin
c
       integer       irb     (   *), jtb     (   *)
       real*8        rvecb   (   *), drvecbi (   *), tvecb   (   *)
     1             , dtvecbi (   *), xb      (   *), yb      (   *)
     2             , qvar    (   *), sintb   (ijkn), costb   (ijkn)
c
       real*8        ar      ( iar,   *), qty     (iqty,   *)
c
c  WARNING:  Do not equivalence the temporary vectors to the worker
c            arrays.  Conflicts result when this routine is called in
c            parallel loops in routines BININT and ZRPTOXYZ.
c
c-----------------------------------------------------------------------
c
       do 10 j=j1,j2
         sintb(j) = sin(tvecb(j))
         costb(j) = cos(tvecb(j))
10     continue
       sinth0 = sin(th0)
       costh0 = cos(th0)
c
       drai = 1.0 / qvar(3)
       rbmn = qvar(1) + qvar(3) / 2.0
       tamn = qvar(4)
       tamx = qvar(5)
       dtai = 1.0 / qvar(6)
       tbmn = qvar(4) + qvar(6) / 2.0
c
c      Re-map the data ("ar" -> "qty").  For each point on the pixel
c  grid, the coordinates of the nearest computational zone are
c  determined using the look-up tables ("irb" and "jtb").  The value
c  assigned to the pixel is then interpolated by use of a van-Leer
c  differencing scheme which has been determined (by trial and error) to
c  yield a much smoother (less jagged) re-map than, say, a weighted
c  average.
c
       do 30 ny=1,nyz
         q1 = yb(ny) * sinth0
         q2 = yb(ny) * costh0
         q3 = yb(ny)**2
         do 20 nx=1,nxz
           xb0    = q1 + xb(nx) * costh0
           sgnxb  = sign ( haf, xb0 )
           yb0    = q2 - xb(nx) * sinth0
           sgnyb  = sign ( one, yb0 )
           rb0    = sqrt ( xb(nx)**2 + q3 )
           ic     = irb (  int ( (rb0-rbmn) * drai ) + 1 )
           sintb0 =   max ( mone,   min ( one, xb0/(rb0+tiny) ) )
           tb0    = sgnyb * asin(sintb0) + ( 1.0 - sgnyb ) * sgnxb * pi
           sgnmn  = sign ( haf, tb0-tamn )
           sgnmx  = sign ( haf, tb0-tamx )
           tb0    = tb0 - ( sgnmn + sgnmx ) * 2.0 * pi
           jc     = jtb (  int ( (tb0-tbmn) * dtai ) + 1 )
           if (ic*jc .eq. 0) then
             qty(nx,ny) = 0.0
           else
             darr       = ( ar(ic+1,jc) - ar(ic  ,jc) ) * drvecbi(ic+1)
             darl       = ( ar(ic  ,jc) - ar(ic-1,jc) ) * drvecbi(ic  )
             d2ar       = darr * darl
             drb        = rb0 - rvecb(ic)
             dqx1       = drb * ( 1.0 + sign ( one, d2ar ) )
     1                  * d2ar / ( ( darr + darl ) + pi*tiny )
             darr       = ( ar(ic,jc+1) - ar(ic,jc  ) ) * dtvecbi(jc+1)
             darl       = ( ar(ic,jc  ) - ar(ic,jc-1) ) * dtvecbi(jc  )
             d2ar       = darr * darl
             dtb        = ( sin(tb0)*costb(jc) - sintb(jc)*cos(tb0) )
             dqx2       = dtb * ( 1.0 + sign ( one, d2ar ) )
     1                  * d2ar / ( ( darr + darl ) + pi*tiny )
             qty(nx,ny) = ar(ic,jc) + dqx1 + dqx2
             icmp1      = ic +  int ( sign ( one, drb ) )
             jcmp1      = jc +  int ( sign ( one, dtb ) )
             qmax       =   max ( ar(ic   ,jc), ar(ic   ,jcmp1)
     1                          , ar(icmp1,jc), ar(icmp1,jcmp1) )
             qmin       =   min ( ar(ic   ,jc), ar(ic   ,jcmp1)
     1                          , ar(icmp1,jc), ar(icmp1,jcmp1) )
             qty(nx,ny) =   max ( qmin,   min ( qmax, qty(nx,ny) ) )
           endif
20       continue
30     continue
c
       return
       end
