
#include "fpreproc/f77_dcomplx.h"

c
c  mod DMC: plotting code conditionally compiled, available w/ __DEBUG
c    macro ONLY...
c
#ifdef __DEBUG
      subroutine conturs(plot2d,rmaj2d,zed2d,title,
     >   nchi,npsi,maxchi,maxpsi)

      use rfxqlo_private
      use sglib_mod
c
c dmc 2 March 2005 -- interactivity added...
c
C============
C idecl:  explicitize implicit INTEGER declarations:
      IMPLICIT NONE
      INTEGER npsi,maxchi,maxpsi,nchi,maxchip1
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zth
C============
      integer mcntur
      parameter (mcntur=50)
c
      integer ncntur,npsim1,nchip1,i,j,icneg,iant,inant
      integer k1,k2
      REAL*8 zrav,zzav,rmjmax,rmjmin,pltmin,pltmax,zedmin,zedmax
      REAL*8 dlwall,dgrid,dcont,zrhw,zrhh,zsiz,zpi
      character*(*) title
      character*72 string
      REAL*8 plot2d(maxchi+1,maxpsi),rmaj2d(maxchi+1,maxpsi),
     %     zed2d(maxchi+1,maxpsi),
     %     raymaj(maxchi+maxpsi),rayzed(maxchi+maxpsi),
     %     clevel(mcntur)
      real*8 zrsum,zzsum,zxi,zjac(2,2),zdum
      data ncntur /10/	  ! actual # of contours to plot ( <= mcntur)
c
      integer iup,ileft,iright,idown
c---------------------------------------------------
c  uread stuff added (dmc)
c
      logical lchar,lxout
      integer lnout,iluntty
c---------------------------------------------------
c  interpolation stuff added (dmc)
c
      integer :: ingrid,idollm
      real*8, dimension(:), allocatable :: zrgrid,zzgrid,zwork
      real*8 :: zraxis,zzaxis,zpaxis
      real*8, dimension(:,:), allocatable :: zidata
      real*8 :: zzrmin,zzrmax,zzzmin,zzzmax,zzdmin,zzdmax,zw
      real*8 :: zzrmid,zzzmid,zzrwid,zzzwid,zr1,zr2,zz1,zz2,zxpand
      real*8, parameter :: ZERO = 0.0d0
      real*8, parameter :: ONE = 1.0d0
c---------------------------------------------------
      iant=selected_antenna
c
      npsim1=npsi-1
      nchip1=nchi+1
      zpi=2.d0*abs(acos(0.d0))
c
      zjac=0
c
      iluntty=lnout(0)
      if(iluntty.lt.0) return
c
      zzrmin=minval(rmaj2d(1:nchi,npsi))
      zzrmax=maxval(rmaj2d(1:nchi,npsi))
      zzrmid = (zzrmin+zzrmax)/2
      zzrwid = (zzrmax-zzrmin)/2

      zzzmin=minval(zed2d(1:nchi,npsi))
      zzzmax=maxval(zed2d(1:nchi,npsi))
      zzzmid = (zzzmin+zzzmax)/2
      zzzwid = (zzzmax-zzzmin)/2

      zzdmin=minval(plot2d(1:nchi,1:npsi))
      zzdmax=maxval(plot2d(1:nchi,1:npsi))
c
      zraxis=sum(rmaj2d(1:nchi,1))/nchi
      zzaxis=sum(zed2d(1:nchi,1))/nchi
      zpaxis=sum(plot2d(1:nchi,1))/nchi
c
      idollm=index(title,'$')-1
      if(idollm.le.0) idollm=max(1,len(trim(title)))
c
      do
         if(lxout(0)) then
            write(iluntty,1001) title(1:idollm)
 1001       format(/' subroutine conturs: plot options for: '/6x,a//
     >         ' T -- traditional, with vacuum vessel, antenna.'/
     >         ' I -- interpolated data only, but rescalable'/
     >         ' Q -- quit.')
c
         endif
         call uread('conturs: plot option:$')
         if(lchar('Q')) then
            exit
         else if(lchar('T')) then
c
c CONTUR routine doesn't work unless you have already set up the
c plotting area and the virtual to screen map.  Do this by first
c plotting the vacuum vessel wall (we could have instead called
c the AGRAPH routine to prepare for CONTUR, but AGRAPH doesn't allow
c to set the x and y limits exacatly, not with xneat/yneat, not with
c agixnt/agiynt).
c
            ileft=150
            iright=725
            idown=125
            iup=700

            call initt(960)
            call binitt
            call fpgrlab(1)

c   first fill raymaj and rayzrd with vacuum vessel shape

            zrsum=0
            zzsum=0
            if(rfx_sys%vv_bdy) then
c  DSPRUCE -- use vv bdy
               call rfx_vvxpand(zdum,nchip1,raymaj,rayzed,
     >            nvvmoms,vvrmoms,vvzmoms)
            else
               do i=1,nchip1
                  zth=(i-1)*2*zpi/(1.0d0*nchi)
                  raymaj(i)=0.d0
                  rayzed(i)=0.d0
c  TORIC/SPRUCE -- use bdy + epsilon
                  zxi=xb_qlo(nx_qlo)
                  call rfmomry1(zxi,zth,0,raymaj(i),rayzed(i),zjac)
                  zrsum=zrsum+raymaj(i)
                  zzsum=zzsum+rayzed(i)
               enddo
               zrsum=zrsum/(nchi+1)
               zzsum=zzsum/(nchi+1)
               do i=1,nchi+1
                  raymaj(i)=(zrsum+1.005*(raymaj(i)-zrsum))/100.0
                  rayzed(i)=(zzsum+1.005*(rayzed(i)-zzsum))/100.0
               enddo
            endif
c
            call npts(nchip1)
            rmjmin=1.d30
            rmjmax=-1.d30
            call r8_mnmx(raymaj(1:nchip1),rmjmin,rmjmax)
            zedmin=1.d30
            zedmax=-1.d30
            call r8_mnmx(rayzed(1:nchip1),zedmin,zedmax)

            call xneat(0)
            call yneat(0)
            call slimx(150,725)
            call slimy(125,700)
            zRav=(rmjmax+rmjmin)/2
            zrhw=(rmjmax-rmjmin)/2
            zrhh=(zedmax-zedmin)/2
            zzav=(zedmax+zedmin)/2
            zsiz=1.1d0*max(zrhw,zrhh)
            call r8_dlimx(zRav-zsiz,zRav+zsiz)
            call r8_dlimy(zzav-zsiz,zzav+zsiz)
            call xfrm(2)
            call yfrm(2)
            call r8_check(raymaj(1:nchip1),rayzed(1:nchip1))
            call r8_dsplay(raymaj(1:nchip1),rayzed(1:nchip1))

c  now draw the plasma boundary

            do i=1,nchi
               raymaj(i)=rmaj2d(i,npsi)
               rayzed(i)=zed2d(i,npsi)
            enddo
            raymaj(nchip1)=raymaj(1)
            rayzed(nchip1)=rayzed(1)
            call line(0)
            call r8_cplot(raymaj(1:nchip1),rayzed(1:nchip1))

            dlwall=.005d0*zedmax
            do i=1,nchi
               dgrid=sqrt((rmaj2d(i,npsi)-rmaj2d(i,npsim1))**2+
     %            (zed2d(i,npsi)-zed2d(i,npsim1))**2)
               raymaj(i)=raymaj(i)+dlwall
     %            *(rmaj2d(i,npsi)-rmaj2d(i,npsim1))/dgrid
               rayzed(i)=rayzed(i)+dlwall*
     %            (zed2d(i,npsi)-zed2d(i,npsim1))/dgrid
            enddo
            raymaj(nchip1)=raymaj(1)
            rayzed(nchip1)=rayzed(1)
            call r8_cplot(raymaj(1:nchip1),rayzed(1:nchip1))

c  Finally draw the antenna slightly thicker

            inant = antenna_array(iant)%nantrz
            call npts(inant)
            call line(0)
            do i=1,inant
               raymaj(i)=antenna_array(iant)%antRgeo(i) !metre
               rayzed(i)=antenna_array(iant)%antZgeo(i) !metre
            enddo
            do j=1,5
               call r8_cplot(raymaj(1:inant),rayzed(1:inant))
               do i=1,inant
                  raymaj(i)=raymaj(i)+0.001d0
               enddo
            enddo

            call frame                  !@

            call agtchs('Rmajor (m)')
            call agtcvs('Z (m)')

c now plot the data:
            maxchip1=maxchi+1
            call npts(nchip1*npsi)
            pltmin=minval(plot2d(1:nchip1,1:npsi))
            pltmax=maxval(plot2d(1:nchip1,1:npsi))
            dcont=(pltmax-pltmin)/AREAL(ncntur)
            if(dcont.eq.0.0d0) dcont=max(abs(pltmin),abs(pltmax))
            if(dcont.eq.0.0d0) dcont=1.0d0
            icneg=-(pltmin/dcont+.5d0)
            clevel=0
            do i=1,ncntur
               clevel(i)=(AREAL(i-icneg)-.5d0)*dcont
            enddo
            k1=ncntur
            k2=icneg+1

            call r8_contur(k1,clevel,k2,
     %         plot2d,maxchip1,rmaj2d,1,nchip1,1,zed2d,1,npsi,1)

            call agtcgt(iup+60,title,0)
            write (string,665) pltmin,pltmax
 665        format('min=',1pe10.3,', max=',1pe9.3,' $')
            call agtcgt(iup+35,string,0)
            call upause(0)
            call finitt(0,0)
         else if(lchar('I')) then
            ingrid=301
            call ureadi(
     >         'conturs: enter interp. gridsize btw 51 and 5001:$',
     >         ingrid)
            ingrid=min(5001,max(51,ingrid))

            allocate(zrgrid(ingrid),zzgrid(ingrid),zwork(ingrid),
     >         zidata(ingrid,ingrid)); zwork=ZERO

            zidata = zzdmin-(zzdmax-zzdmin)/10  ! value for (R,Z) pts outside.

            zxpand=1.0d0 + 1.0d0/ingrid
            zr1=zzrmid-zxpand*zzrwid
            zr2=zzrmid+zxpand*zzrwid
            zz1=zzzmid-zxpand*zzzwid
            zz2=zzzmid+zxpand*zzzwid
            do i=1,ingrid
               zrgrid(i)=((ingrid-i)*zr1 + (i-1)*zr2)/(ingrid-1)
               zzgrid(i)=((ingrid-i)*zz1 + (i-1)*zz2)/(ingrid-1)
            enddo
            do i=1,ingrid
               if((zzrmin.lt.zrgrid(i)).and.(zrgrid(i).lt.zzrmax)) then
                  call interp
               endif
            enddo

            call r8_grf3f1(zrgrid,zidata,zzgrid,zwork,
     >         ingrid,ingrid,ingrid,'R',title(1:idollm),'Z',
     >         'm',' ','m',title(1:idollm),rfx_sys%run_id,0)

            deallocate(zrgrid,zzgrid,zidata,zwork)
         else
            call uermsg(' %conturs: requested option unrecognized.',
     >         ' ',0)
         endif
      enddo
c
      return

      contains
        subroutine interp
c
c  for column of points along line R=zrgrid(i), find interpolation
c  of the data.  Some points maybe outside, but the line R=zrgrid(i)
c  is known to intersect the mapped domain.
c
        integer ichi,ichi_zmin,ichi_zmax
        integer jpsi
        real*8 zr,zrdiff,zzmincep,zzmaxcep,zcep,zxchi,zxpsi,zzcep
        logical :: ifound
c--------------------------
        jpsi=maxpsi
        zr=zrgrid(i)
c
        zzmincep=zzgrid(ingrid); zzmaxcep=zzgrid(1)
c
        do ichi=1,nchi
           zrdiff=rmaj2d(ichi+1,jpsi)-rmaj2d(ichi,jpsi)
           if(zrdiff.eq.ZERO) cycle
           if((zr-rmaj2d(ichi+1,jpsi))*(zr-rmaj2d(ichi,jpsi)).le.ZERO)
     >        then
              zw=(zr-rmaj2d(ichi,jpsi))/zrdiff
              zcep=zed2d(ichi,jpsi)*(ONE-zw)+zed2d(ichi+1,jpsi)*zw
              if(zcep.lt.zzmincep) then
                 zzmincep=zcep
                 ichi_zmin=ichi
              endif
              if(zcep.gt.zzmaxcep) then
                 zzmaxcep=zcep
                 ichi_zmax=ichi
              endif
           endif
        enddo
c
        ichi=0
        zzcep=zzgrid(1)
        do j=1,ingrid
           if(zzgrid(j).lt.zzmincep) cycle
           if(zzgrid(j).lt.zzcep) cycle
           if(zzgrid(j).gt.zzmaxcep) exit

           if(ichi.eq.0) then
              ichi=ichi_zmin
              jpsi=npsi
           endif

c  search, from box [jpsi-1:jpsi]x[ichi:ichi+1]

           call find(zr,zzgrid(j),ichi,jpsi,zxchi,zxpsi,ifound,zzcep)
           if(.not.ifound) cycle

           if(jpsi.gt.1) then
              zidata(i,j) = (ONE-zxchi)*
     >           ((ONE-zxpsi)*plot2d(ichi,jpsi-1) +
     >           zxpsi*plot2d(ichi,jpsi)) +
     >           zxchi*
     >           ((ONE-zxpsi)*plot2d(ichi+1,jpsi-1) +
     >           zxpsi*plot2d(ichi+1,jpsi))
           else
              zidata(i,j) = (ONE-zxchi)*
     >           ((ONE-zxpsi)*zpaxis +
     >           zxpsi*plot2d(ichi,jpsi)) +
     >           zxchi*
     >           ((ONE-zxpsi)*zpaxis +
     >           zxpsi*plot2d(ichi+1,jpsi))
           endif
        enddo
        return
        end subroutine interp

        subroutine find(zr,zz,ichi,jpsi,zxchi,zxpsi,ifound,zzcep)

        ! find the containing cell and the displacements within the cell
        ! the caller has assured that the zz value of the target along the
        ! line zr=const, lies between the min and max intercept of the
        ! mapped region with that line.

        real*8, intent(in) :: zr,zz  ! target point
        integer, intent(inout) :: ichi,jpsi  ! outer segment of zone
                                     ! guess on input, actual on output
        real*8, intent(out) :: zxchi,zxpsi   ! displacements within zone
        logical, intent(out) :: ifound ! .TRUE. if cell is found.
        real*8, intent(inout) :: zzcep ! re-entry intercept, if any

        !------------------------------
        real*8 :: zrbox(5),zzbox(5),ztest,ztest2,zzside(2),zxside(2)
        real*8 :: zxchix(2),zxpsix(2),zztol
        integer :: iside(2),jside,isides,idir
        integer :: icorner(2)
        logical :: icorner1
        !------------------------------

        zztol=1.0d-12*zzzwid

        ifound=.TRUE.
        do
           zrbox(1)=rmaj2d(ichi,jpsi)
           zzbox(1)=zed2d(ichi,jpsi)
           zrbox(2)=rmaj2d(ichi+1,jpsi)
           zzbox(2)=zed2d(ichi+1,jpsi)
           if(jpsi.gt.1) then
              zrbox(3)=rmaj2d(ichi+1,jpsi-1)
              zzbox(3)=zed2d(ichi+1,jpsi-1)
              zrbox(4)=rmaj2d(ichi,jpsi-1)
              zzbox(4)=zed2d(ichi,jpsi-1)
           else
              zrbox(3)=zraxis
              zzbox(3)=zzaxis
              zrbox(4)=zraxis
              zzbox(4)=zzaxis
           endif
           zrbox(5)=zrbox(1); zzbox(5)=zzbox(1)

           isides=0
           do jside=1,4
              if(zrbox(jside).eq.zrbox(jside+1)) cycle

              ztest=(zr-zrbox(jside))*(zr-zrbox(jside+1))
              if(ztest.ge.0) cycle

              isides=isides+1
              iside(isides)=0
              icorner(isides)=0
              zxside(isides)=ZERO

              if(ztest.eq.ZERO) then
                 ! hit a corner
                 if((zr-zrbox(jside)).eq.ZERO) then
                    icorner(isides)=jside
                    zzside(isides)=zzbox(jside)
                 else
                    if(jside.eq.4) then
                       icorner(isides)=1
                    else
                       icorner(isides)=jside+1
                    endif
                    zzside(isides)=zzbox(icorner(isides))
                 endif
                 if(isides.gt.1) then
                    if(icorner(isides).eq.icorner(isides-1)) then
                       ! skip duplicate corner count
                       isides=isides-1
                       cycle
                    endif
                 endif
              else
                 ! hit an edge
                 iside(isides)=jside
                 zxside(isides)=(zr-zrbox(jside))/
     >              (zrbox(jside+1)-zrbox(jside))
                 zzside(isides)=zzbox(jside)+zxside(isides)*
     >              (zzbox(jside+1)-zzbox(jside))
              endif

              if(isides.eq.2) exit
           enddo

           icorner1=.FALSE.
           if(isides.eq.1) then
              if(icorner(1).eq.0) then
                 ! if only hit one point, it *must* be a corner...
                 call errmsg_exit('?conturs search error 1.')
              else
                 icorner1=.TRUE.
                 isides=2
                 iside(2)=0
                 icorner(2)=icorner(1)
                 zzside(2)=zzside(1)
              endif
           endif

           ztest=(zz-zzside(1))*(zz-zzside(2))
           ztest2=min(abs(zz-zzside(1)),abs(zz-zzside(2)))
           if((ztest.le.ZERO).or.(ztest2.lt.zztol)) then
              do jside=1,isides
                 if(icorner(jside).eq.0) then
                    if(iside(jside).eq.1) then
                       zxpsix(jside)=ONE
                       zxchix(jside)=zxside(jside)
                    else if(iside(jside).eq.2) then
                       zxpsix(jside)=ONE-zxside(jside)
                       zxchix(jside)=ONE
                    else if(iside(jside).eq.3) then
                       zxpsix(jside)=ZERO
                       zxchix(jside)=ONE-zxside(jside)
                    else if(iside(jside).eq.4) then
                       zxpsix(jside)=zxside(jside)
                       zxchix(jside)=ZERO
                    else
                       call errmsg_exit('?conturs search error 2.')
                    endif
                 else if(icorner(jside).eq.1) then
                    zxpsix(jside)=ONE
                    zxchix(jside)=ZERO
                 else if(icorner(jside).eq.2) then
                    zxpsix(jside)=ONE
                    zxchix(jside)=ONE
                 else if(icorner(jside).eq.3) then
                    zxpsix(jside)=ZERO
                    zxchix(jside)=ONE
                 else if(icorner(jside).eq.4) then
                    zxpsix(jside)=ZERO
                    zxchix(jside)=ZERO
                 endif
              enddo
              if(zzside(1).eq.zzside(2)) then
                 zw=0.5d0
              else
                 zw=(zz-zzside(1))/(zzside(2)-zzside(1))
                 zw=max(ZERO,min(ONE,zw))
              endif
              zxpsi=zxpsix(1)+zw*(zxpsix(2)-zxpsix(1))
              zxchi=zxchix(1)+zw*(zxchix(2)-zxchix(1))
              exit                      ! *** answer found ***
           endif

           ! not found; have to move to next cell...
           ! convert corners to sides...

           if(icorner1) then
              if(icorner(1).eq.1) then
                 iside(1)=4
                 if(zzbox(4).gt.zzbox(2)) then
                    zzside(1)=zzbox(1)
                    zzside(2)=zzbox(1)-zztol
                 else
                    zzside(1)=zzbox(1)-zztol
                    zzside(2)=zzbox(1)
                 endif
              else
                 iside(1)=icorner(1)-1
                 if(zzbox(icorner(1)-1).gt.zzbox(icorner(1)+1)) then
                    zzside(1)=zzbox(1)
                    zzside(2)=zzbox(1)-zztol
                 else
                    zzside(1)=zzbox(1)-zztol
                    zzside(2)=zzbox(1)
                 endif
              endif
              zxside(1)=ONE
              iside(2)=icorner(1)
              zxside(2)=ZERO
           else
              do jside=1,2
                 if(icorner(jside).gt.0) then
                    iside(jside)=icorner(jside)
                    zxside(jside)=ZERO
                 endif
              enddo
           endif

           ! exit via the side that has the larger z value
           ! exception: if jpsi=maxpsi, we have to find re-entry point...

           if(zzside(1).gt.zzside(2)) then
              jside=1
           else
              jside=2
           endif

           if(iside(jside).eq.2) then
              ichi=ichi+1
              if(ichi.eq.maxchi+1) ichi=1
           else if(iside(jside).eq.4) then
              ichi=ichi-1
              if(ichi.eq.0) ichi=maxchi
           else if(iside(jside).eq.3) then
              if(jpsi.eq.1) then
                 ! cannot exit through origin, must choose ichi neighbor
                 if(zzbox(1).gt.zzbox(2)) then
                    ichi=ichi-1
                    if(ichi.eq.0) ichi=maxchi
                 else
                    ichi=ichi+1
                    if(ichi.eq.maxchi+1) ichi=1
                 endif
              else
                 jpsi=jpsi-1
              endif
           else
              if(jpsi.lt.maxpsi) then
                 jpsi=jpsi+1
              else
                 ! crossed bdy; must look for re-entry point; may miss
                 ! target point as a result
                 if(zzbox(1).gt.zzbox(2)) then
                    idir=-1
                 else
                    idir=1
                 endif
                 do
                    ichi=ichi-idir
                    if(ichi.eq.0) ichi=maxchi
                    if(ichi.eq.maxchi+1) ichi=1
                    if(rmaj2d(ichi,jpsi).eq.rmaj2d(ichi+1,jpsi)) cycle
                    ztest=(zr-rmaj2d(ichi,jpsi))*
     >                 (zr-rmaj2d(ichi+1,jpsi))
                    if(ztest.gt.ZERO) cycle

                    ! zr intercept found

                    zw=(zr-rmaj2d(ichi,jpsi))/
     >                 (rmaj2d(ichi+1,jpsi)-rmaj2d(ichi,jpsi))
                    ztest=zed2d(ichi,jpsi)+zw*
     >                 (zed2d(ichi+1,jpsi)-zed2d(ichi,jpsi))
                 enddo
                 if(ztest.gt.zz) then
                    ! target was in gap
                    ifound=.FALSE.
                    zxchi=ZERO
                    zxpsi=ZERO
                    zzcep=ztest  ! other targets may also be in gap...
                    exit
                 endif
              endif  ! (jpsi.eq.maxpsi)
           endif  ! (iside test)

        enddo  ! keep searching...
        return
        end subroutine find
      end

#else
      subroutine conturs

      write(6,*) ' %NOTE: dummy CONTURS call: compiled w/o __DEBUG '

      return
      end
#endif
