!******************** start file nbplot.for ; group nbplot *************
!  stplot
!
!  control plotting output ("c%" code)
!
subroutine stplot(inew,inu,ict,iplcx,iplfp,isty,ipflr, &
     ztim)
  use iso_c_binding, only: fp => c_double
  !
  !  inew=-1 for deposition plotting only
  !  inew=1 if this is a new ion
  !  inu:  return # of plotting units to send
  !      1 unit= half orbit
  !  ict:  this counter gets cleared to zero
  !  iplfp: =1 to plot fp calculation locations
  !  iplcx: =1 to plot cx calculation locations
  !  isty:  raster location for labels in right margin
  !
  use nbi_com
  use xstraln_calls
  !%      use sglib_mod
  !
  implicit none
  integer :: inew,inu,ict,iplcx,iplfp,isty,ipflr
  real(fp) :: ztim
  !============
  !%      integer :: init,iuntrm,lnout,ierr
  !%      integer :: imomry,iert,ivsp,idcod,ix,i,ip1,imax,inum,izp1,ij,iabs
  !============
  !%      real(fp) :: zfac,zep,zzr,zzz,r8dcod,zphiarg,zxmin,zymin,zlsq
  !%      real(fp) :: zvmax,zxmax,zymax,zth,zx,zdummy,zdtol,zr,zy,zdum
  !%      real, parameter :: zsymsize=0.7    ! symbol size (real not real(fp) ::)
  !============
  !%      logical lxout, lchar, lyesc
  !%      character*80      zbuf
  !%      character*1 ians
  !%      data init/0/
  !     draw limiter
  !
  ! ... flag if momry is used - can modify m.c. ptcl position !
  !%      iuntrm=lnout(0)
  !%      imomry=0
  !
  !%      ict=0
  !%      inu=3
  !%      iplfp=0
  !%      iplcx=0
  !%      if(nbplot.eq.-1) then
  !%         isty=750
  !%         return
  !%      end if
  !%      if(init.gt.0) goto 30
  !%      nbplot=1
  !%      zfac=1._fp
  !%      if(inew.eq.-1) goto 100
  !%      goto 40
  !% 30   continue
  !%      if(inew.eq.-1) goto 100
  !%      if(inew.eq.1) goto 35
  !%      isty=isty-15
  !%      write(zbuf,4003)
  !% 4003 format('*final:$')
  !%      call rmrglb(zbuf,isty)
  !%      write(zbuf,4005)
  !%      call stplab(800,isty,zbuf,ztim)
  !% 35   continue
  !%      if(nbplot.gt.0) call nbplot_bdy(0)
  !% 40   continue
  !------------------------------------------------------------
  !%      zep=1.04414e-12_fp*0.5_fp*vion*vion*abeams(lsbeam)
  !%      call xmomry(xion, th, 0, zzr, zzz, xbjaco)
  !%      if(inew.eq.1) write(nonlin,3004) xion,th,zep,xksid,zzr,zzz
  !% 3004 format(/' new ptcl:  xi,th  = ',2(1pe11.4,1x)/ &
  !%             '      e0(ev) ,ksid = ',2(1pe11.4,1x)/  &
  !%             '              r,z  = ',2(1pe11.4,1x)/ &
  !%             '  enter # of plotting units = 0 to skip ptcl.'/ &
  !%             ' *enter "c" to change orbit parameters ***'/ &
  !%             ' *precede number with "v" for v space plotting',/, &
  !%             ' *enter "e" to toggle to total energy/p-phi display',/, &
  !%             '            (orbs & orbg only, not for orbc)'/ &
  !%             ' *enter "f" to modify the plotting of flr symbols'/ &
  !%           '')
  !% 307  call uread( &
  !%     ' enter # of plotting units to draw (1 unit= one half orbit):$')
  !%      if(lchar('C')) then
  !%         print *, ' *precede number with "r" for r,z coordinate input'
  !%         write(nonlin,3005) 'xi', xion
  !%         write(nonlin,3005) 'r ', zzr
  !% 3005    format(' [ old ', a2, ' value:', 1pe15.5,']')
  !%         call uread('nbplot:  new "xi":$$')
  !%         if (lchar('R')) then
  !%            zzr = r8dcod(iert)
  !%            if (iert .ne. 0) goto 307
  !%            call ureadr8('nbplot:  new "z":$$',zzz)
  !%            zphiarg=0.0_fp
  !%            call xs_fxget(0.01_fp,zphiarg,zzr,zzz, xion, th ,iert)
  !%            if (iert. ne. 0) then
  !%               print *, ' !!error in inverse map xs_fxget'
  !%               goto 307
  !%            end if
  !%         else
  !%            xion = r8dcod(iert)
  !%            if (iert .ne. 0) goto 307
  !%            call ureadr8('nbplot:  new "th":$$',th)
  !%         end if
  !%         call ureadr8('nbplot:  new "kinetic energy":$$',zep)
  !%         call ureadr8('nbplot:  new "ksid":$$',xksid)
  !%         call ureadr8('nbplot:  new "goose":$$',goose)
  !%c   don't adjust for electrostatic potential: zep is kinetic energy
  !%         vion=sqrt(2.0_fp*max(1.0_fp,zep)/ &
  !%                  (1.04414e-12_fp*abeams(lsbeam)))
  !%c   make sure inputs to orbg are defined
  !%         vpll = vion*xksid
  !%         bivpar(1) = vpll
  !%         vperp = sqrt(vion*vion-vpll*vpll)
  !%         bivprp(1) = vperp
  !%c   force magnetic field & jacobian to be correct for orbg & orbs
  !%         bixi(1) = xion
  !%         bithst(1) = th
  !%         call bintrp
  !%         call jintrp
  !%         call gintrp(1)
  !%         goto 40
  !%      else if(lchar('e')) then
  !%         nbplist = -nbplist
  !%         goto 40
  !%      else if(lchar('f')) then
  !%         print *, 'enter 0 for no flr symbol, nonzero integer :: to ' &
  !%           // 'control the frequency'
  !%         print *, 'of flr symbol plotting'
  !%         print *, 'if nlbgflr=.false., use a negative value to also' &
  !%           // ' plot the flr '
  !%         print *, 'returned from a separate call to map_larmor'
  !%         call ureadi('enter flr plot control integer ::$',ipflr)
  !%         goto 40
  !%      else if(lchar('v')) then
  !%         ivsp=1
  !%      else
  !%         ivsp=0
  !%      end if
  !%      if(lchar('@')) call nbpcancl
  !%      inu=idcod(iert)
  !%      if(iert.ne.0) goto 307
  !%      nbplot=1
  !%      if(ivsp.eq.1) nbplot=2
  !%      if(inu.gt.0) goto 308
  !  turn off plotting for this ptcl
  !%      nbplot=0
  !%      goto 9000
  !% 308  call uread('enter scale factor (1.0 ==> std full scale)$')
  !%      zfac=r8dcod(iert)
  !%      zfac=min(2.0_fp,zfac)
  !-----------------------------------/hand-set scale
  !%      if(zfac.le.0.0_fp) then
  !%         if(lxout(0)) write(iuntrm,7059)
  !% 7059    format(/ &
  !%           ' enter "x" not to change a limit on the graph'/)
  !% 71      continue
  !%         if(lxout(0)) write(iuntrm,7101) zxmin
  !% 7101    format('  old xmin (square scale)=',1pe10.3)
  !%         call uread(' sub:grfah2 >> enter new xmin (e10.3)$')
  !%         if(lchar('x')) goto 73
  !%         zxmin=r8dcod(iert)
  !%         if(iert.ne.0) goto 71
  !% 73      if(lxout(0)) write(iuntrm,7301) zymin
  !% 7301    format('  old ymin (square scale)=',1pe10.3)
  !%         call uread(' sub:grfah2 >> enter new ymin$')
  !%         if(lchar('x')) goto 74
  !%         zymin=r8dcod(iert)
  !%         if(iert.ne.0) goto 73
  !% 74      if(lxout(0)) write(iuntrm,7401) zlsq
  !% 7401    format('  old square edge length=',1pe10.3)
  !%         call uread(' sub:grfah2 >> enter new square edge length:$')
  !%         if(lchar('x')) goto 801
  !%         zlsq=r8dcod(iert)
  !%         if(iert.ne.0) goto 74
  !% 801     continue
  !%      end if
  !%      if(iert.ne.0) goto 308
  !%      if(lxout(0)) write(iuntrm,7501)
  !% 7501 format(' ** type "x" to turn off a calculation entirely'/)
  !%      call uread('enter "1" to turn on cx calculation plotting:$')
  !%      iplcx=0
  !%      if(lchar('1')) iplcx=1
  !%      if(lyesc('x')) iplcx=-1
  !%      call uread('enter "1" to turn on fp calculation plotting:$')
  !%      iplfp=0
  !%      if(lchar('1')) iplfp=1
  !%      if(lyesc('x')) iplfp=-1
  !% 100  continue
  !%      call nbplot_bdy(1)
  !%      call ubinit
  !%      call gftics
  !%      call slimx(200,780)
  !%      call slimy(100,700)
  !%      call sizes(zsymsize)
  !%      call newpag
  !%      if(inew.eq.1) call agtmhs(15,720,'new orbit$')
  !%      if(inew.eq.-1) call agtmhs(15,720,'deposition$')
  !%      call agtmhs(15,743,runid(1:lrunid)//':$')
  !%      call agtmhs(15,605,' (cm)$')
  !%      init=init+1
  !%      do 9977 ix=1,5
  !%         if((ix.eq.1).and.(zfac.gt.0.0_fp)) then
  !  scale calculation
  !%            if(nbplot.eq.2) then
  !%               zvmax=1.1_fp*vion
  !%               zxmin=-zvmax
  !%               zxmax=zvmax
  !%               zymin=0.0_fp
  !%               zymax=zvmax
  !%            else
  !%               zxmin=raxis
  !%               zxmax=zxmin
  !%               zymin=yaxis
  !%               zymax=0.0_fp
  !%            end if
  !%         end if
  !%         do 9978 i=1,201
  !%            ip1=i+1
  !  draw plasma boundary
  !%            if(nbplot.eq.2) then
  !%               zth=3.141593_fp*(i-1.0_fp)/200._fp
  !%               zx=1.0_fp-(ix-1)*0.2_fp
  !%               xpos(ip1)=vion*zx*cos(zth)
  !%               ypos(ip1)=vion*zx*sin(zth)
  !%            else
  !%               zx=1.0_fp-(ix-1)*0.2_fp
  !%               zth=3.141593_fp*(i-1._fp)/100._fp
  !%               call xmomry(zx,zth,0,xpos(ip1),ypos(ip1),zdummy)
  !%               imomry=1
  !%               if((ix.eq.1).and.(zfac.gt.0.0_fp)) then
  !  scale calculation
  !%                  zxmin=min(xpos(ip1),zxmin)
  !%                  zxmax=max(xpos(ip1),zxmax)
  !%                  zymin=min(ypos(ip1),zymin)
  !%                  zymax=max(ypos(ip1),zymax)
  !%               end if
  !%            end if
  !% 9978    continue
  !%         xpos(1)=201._fp
  !%         ypos(1)=201._fp
  !%         if((ix.eq.1).and.(zfac.gt.0.0_fp).and.(nbplot.eq.1)) then
  !  scale calculation
  !%            zlsq=max((zymax-zymin),(zxmax-zxmin))*zfac
  !%            zxmin=0.5_fp*(zxmin+zxmax-zlsq)
  !%            zymin=0.5_fp*(zymin+zymax-zlsq)
  !%         end if
  !%         if(ix.eq.1) then
  !%            if(nbplot.eq.1) then
  !%               zxmax=zxmin+zlsq
  !%               zymax=zymin+zlsq
  !%            end if
  !%            call xneat(1)
  !%            call yneat(1)
  !%            call r8_dlimx(zxmin,zxmax)
  !%            call r8_dlimy(zymin,zymax)
  !%            call line(1)
  !%            call r8_check(xpos,ypos)
  !%            call r8_dsplay(xpos,ypos)
  !%            call gfticx
  !%         else
  !%            call line(1)
  !%            call r8_cplot(xpos,ypos)
  !%         end if
  !% 9977 continue
  ! -------------------------
  ! draw xplasma limiters
  !%      imax  = minbplt-1                 ! dimension-1 on xpos,ypos
  !%      zdtol = .001_fp*.01_fp*rmnrmp(lep1,1)
  !%      call eq_limcon(imax, inum, xpos(2), ypos(2), zdtol, ierr)
  !%      if (ierr==0) then
  !%         xpos(2:inum+1)=100._fp*xpos(2:inum+1) ! m to cm
  !%         ypos(2:inum+1)=100._fp*ypos(2:inum+1)
  !%         xpos(1)=inum
  !%         ypos(1)=inum
  !%         call line(0)
  !%         call r8_cplot(xpos,ypos)
  !%      else
  !%         print *, '?nbplot: error calling eq_limcon, ' &
  !%           //'can not draw limiter'
  !%      end if
  ! -------------------------
  !  draw circular limiters
  !c%      if((ncirlm.gt.0).and.(nbplot.eq.1)) then
  !c%        xpos(1)=201.
  !c%         ypos(1)=201.
  !c%        do 703 i=1,ncirlm
  !c%          zx0=crlmr1(i)
  !c%          zy0=crlmy1(i)-plhgt
  !c%          zr=crlmrd(i)
  !c%          call pcircl(zx0,zy0,zr,xpos,ypos,202)
  !c% 703        continue
  !c%      end if
  !  draw line limiters
  !c%      if((nlinlm.gt.0).and.(nbplot.eq.1)) then
  !c%        xpos(1)=3.
  !c%        ypos(1)=3.
  !c%        do 705 i=1,nlinlm
  !c%          zx0=alnlmr(i)
  !c%          zy0=alnlmy(i)-plhgt
  !c%          zcth=alnlmc(i)
  !c%          zsth=alnlms(i)
  !c%          zl=10.*raxis
  !c%          xpos(2)=zx0+zl*zcth
  !c%          xpos(3)=zx0
  !c%          xpos(4)=zx0-zl*zcth
  !c%          ypos(2)=zy0+zl*zsth
  !c%          ypos(3)=zy0
  !c%          ypos(4)=zy0-zl*zsth
  !c%          call line(0)
  !c%          call r8_cplot(xpos,ypos)
  !c% 705        continue
  !c%      end if
  ! ------------------------
  !  draw midplane
  !%      if(nbplot.eq.1) then
  !%         izp1=2*nzones+1
  !%         xpos(1)=izp1
  !%         ypos(1)=izp1
  !%         do i=1,izp1
  !%            ij=lcentr+iabs(nzones+1-i)
  !%            if(i.ge.nzones+1) then
  !%               xpos(i+1)=rmjrmp(ij,1)+rmnrmp(ij,1)
  !%            else
  !%               xpos(i+1)=rmjrmp(ij,1)-rmnrmp(ij,1)
  !%            end if
  !%            ypos(i+1)=ympx(ij)
  !%         end do
  !%         call line(1)
  !%         call r8_cplot(xpos,ypos)
  !%      end if
  !
  !%      if(inew.eq.-1) then
  !%         if(nlfprod(lsbeam)) then
  !%            call agtmhs(200,720,'fusion product fast ion deposition$')
  !%            call agtmhs(100,28,'deposition sce pts (x''s) w/flr gyros$')
  !%         else
  !%            call agtmhs(200,720,'solid lines:  neutral tracks$')
  !%            call agtmhs(200,28,'x''s:  deposition (guiding centers)$')
  !%         end if
  !%         goto 9979
  !%      else
  !%         call agtmhs(200,720, &
  !%           'dashed lines:  neutral tracks$')
  !%         if(iplcx.eq.1) call agtmhs(200,25, &
  !%           'x''s:  charge exchange calculations$')
  !%         if(iplfp.eq.1) call agtmhs(200,5, &
  !%           'stars:  collision calculations$')
  !%      end if
  !%      isty=750
  !%      write(zbuf,4001)
  !% 4001 format('*initial:$')
  !%      call rmrglb(zbuf,isty)
  !%      write(zbuf,4005)
  !% 4005 format('orbit time:$')
  !%      call stplab(800,isty,zbuf,ztim)
  !%      if(nlfprod(lsbeam)) then
  !%         call agtmhs(200,743,'solid line:  fusion ion orbit track$')
  !%      else
  !%         call agtmhs(200,743,'solid line:  beam ion orbit track$')
  !%      end if
  !% 9979 continue
  !% 9000 continue
  ! exit
  !%      if(imomry.eq.1) call xmomry(xion,th,0,zr,zy,zdum)
  !============
  return
end subroutine stplot
!
!%      subroutine nbpcancl
!  call this to cancel the next plot control menu
!%       use iso_c_binding, only: fp => c_double
!%      implicit none
!%      call urdtti('0 ')
!%      return
!%      end
!%c-------
!  rmrglab-- write out values of certain common variables in
!    right margin
!  isty-- raster location to start
!
!%      subroutine rmrglb(zbuf,isty)
!%      use nbi_com
!============
!%      implicit none
!%      integer :: isty,ion,jx
!============
!%      real(fp) :: zf1d,zpsint
!%      real(fp) :: zev,zx0,zh,zkonst,zphi,zpsip,zg,zvpar
!%      real(fp) :: zpphi
!============
!%      parameter (ion=1)
!%      character(len=*)      zbuf
!%
!
!%c==============================================================
!
!%      isty=isty-23
!%      call agtmhs(800,isty,zbuf)
!%
!%      if(nbplist.eq.1) then             !plot kinetic energy
!%
!%         zev=1.04414e-12_fp*0.5_fp*vion*vion*abeam
!%         write(zbuf,4003)
!% 4003    format('k.eion(ev):$')
!%         call stplab(800,isty,zbuf,zev)
!%         write(zbuf,4004)
!% 4004    format('ksid:$')
!%         call stplab(800,isty,zbuf,xksid)
!%         if(nlheavy(lsbeam)) then
!%            write(zbuf,4016)
!% 4016       format('xzbeami:$')
!%            call stplab(800,isty,zbuf,xzbeami)
!%         else
!%            write(zbuf,4006)
!% 4006       format('wght:$')
!%            call stplab(800,isty,zbuf,wght)
!%         end if
!%         write(zbuf,4021)
!% 4021    format('goose:$')
!%         call stplab(800,isty,zbuf,goose)
!%         write(zbuf,4022)
!% 4022    format('fpskip:$')
!%         call stplab(800,isty,zbuf,fpskip)
!%         write(zbuf,4023)
!% 4023    format('cxskip:$')
!%         call stplab(800,isty,zbuf,cxskip)
!%
!%      else                              !plot total energy, p-phi
!%
!%         jx = int(bixi(ion)/dxi(lcentr,2)) + 1
!%         zx0 = real(jx-1,fp)*dxi(lcentr,2)
!%         zh = bixi(ion) - zx0
!%         zkonst = 1.e8_fp*tflux/xpi
!%         zphi  = zf1d(zh,splphi(1,jx),splphi(2,jx), &
!%           splphi(3,jx),splphi(4,jx))
!%         zpsip = zpsint(plflxi(jx),zkonst,zx0,zh,spliot(1,jx), &
!%           spliot(2,jx),spliot(3,jx),spliot(4,jx))
!%         zg = zf1d(zh,splg(1,jx),splg(2,jx), &
!%           splg(3,jx),splg(4,jx))
!%         call bintrp
!     %
!c%  jphe 16-nov-1993 change energy units to ev
!%
!%         zev = (0.5_fp*zmp*abeam*vion*vion + zel*zphi)/1.602e-19_fp
!%         write(zbuf,5003)
!% 5003    format('t.eion:ev:$')
!%         call stplab(800,isty,zbuf,zev)
!%         zvpar = vion*xksid
!%         zpphi = zmp*abeam*zvpar*zg*bi1ovb(ion) - (zel/zc)*zpsip
!%         write(zbuf,5004)
!% 5004    format('pphi:$')
!%         call stplab(800,isty,zbuf,zpphi)
!%         write(zbuf,5005)
!% 5005    format('dtrap:$')
!%         call stplab(800,isty,zbuf,dtrap)
!%         write(zbuf,5021)
!% 5021    format('goose:$')
!%         call stplab(800,isty,zbuf,goose)
!%         write(zbuf,5022)
!% 5022    format('xi:$')
!%         call stplab(800,isty,zbuf,xion)
!%         write(zbuf,5023)
!% 5023    format('theta:$')
!%         call stplab(800,isty,zbuf,th)
!%
!%      end if
!%      return
!%      end
!%c----------/stplab:  generate labels, right margin
!%c
!%      subroutine stplab(ixl,iylf,zlab,zval)
!  iylf decremented by 23=line height (rasters), each call
!============
!%      use iso_c_binding, only: fp => c_double
!%      implicit none
!%      integer :: iylf,ixl
!============
!%      real(fp) :: zval
!============
!%      character(len=*) zlab
!  zlab destroyed by call ***
!%      iylf=iylf-23
!%      call agtmhs(ixl,iylf,zlab)
!%      write(zlab,4001) zval
!% 4001 format(1pe11.4,'$')
!%      iylf=iylf-23
!%      call agtmhs(ixl,iylf,zlab)
!%      return
!%      end
!-----------------------------------------------
!%      subroutine pcircl(x0,y0,r,x,y,n)
!            draw a circle
!%      use iso_c_binding, only: fp => c_double
!%      use sglib_mod
!============
!%      implicit none
!%      integer :: n,i
!============
!%      real(fp) :: y0,r,x0,zth
!============
!%      real(fp) :: x(n),y(n)
!%      x(1)=real(n-1,fp)
!%      y(1)=x(1)
!%      do 10 i=2,n
!%         zth=6.28318_fp*(i-2)/real(n-2,fp)
!%         x(i)=x0+r*cos(zth)
!%         y(i)=y0+r*sin(zth)
!% 10   continue
!%      call line(0)
!%      call r8_cplot(x,y)
!%      return
!%      end
!---------------------------------
!%      subroutine nbplot_bdy(iarg)
!%      use iso_c_binding, only: fp => c_double
!%      implicit none
!%      integer, intent(in) :: iarg  ! 1 for start of new plot, -1 for end
!
!  keep track of ubinit/upaust calls -- set flag to 1 at start of new plot
!  (iarg=1); at end set to zero, but if not zero already, generate
!  upaust(2) call (get user acknl).
!
!%      integer, save :: iflag = 0
!
!%      if(iarg.eq.1) then
!%         iflag = 1
!%      else
!%         if(iflag.eq.1) call upaust(2)
!%         iflag = 0
!%      end if
!
!%      return
!%      end
!********************end file nbplot.for ; group nbplot ***************
! 14sep2004 fgtok
