!-----------------------------------------------------------------------
!   oreplt -- orbit error plotting:
!      used in debug mode to scan starting values of 2 of the following:
!            xi, vpar/v, and energy of fast ions,
!      follow for a specified no. of midplane crossings, and
!      plot the error in radius, energy, and/or momentum
!
!       john schivell       27 apr 94
!
!
!--   mods
!       jec  25-may-1999  jconboy@jet.uk
!                         f90 std open, ifdef __f90
!       jfs  28 feb 95 -- write out file of error statistics [optional]
!
!       jfs  19 jan 95 -- added plots of xi errors est. from p errors
!
!       jfs  26 aug 94 -- added plot of time-step magnification factor,
!                            turned off making no. of half-orbits even
!
!       jfs  15 oct 94 -- added statistics on xi errors
!
!       jfs   9 dec 94 -- added "data" ztmfst: starting value for tmfcur
!
!=============================================================================
!

subroutine oreplt
#ifdef __defunct
  use nbi_com
  use nbi_dimensions

  !      common block for communication with orbit:
  !xxx      include 'oepcom'

  implicit none
  LOGICAL, EXTERNAL :: LCHAR,LXOUT,LYES,LYESC                !
  !      storage for orbit errors, independent variables
  !============
  integer :: mv1,mv2,kltmfp,klst,lnout,iltmfp,isb,il,isel,ilmxcr
  integer :: jfix,kxi,klam,ken,kcr,iplcx_sv,iplcx,iplfp_sv,iplfp
  integer :: iplcm_sv,iplcm,iorb_loc,ixi,ilam,ien,is1,is2,icr
  integer :: ind_fake,j,lfbadj,i3sd,ilst,kplot
  !============
  real(fp) :: zpi,ztmfst,zthst,zdxicut,zeion,zefix,zxifix,zlamlo
  real(fp) :: zlamhi,zdlam,zelo,zehi,zde,zlamfix,zxilo,zxihi,zdxi
  real(fp) :: zdelt_sv,zdxis0,zdxis1,zdxis2,zs0,zs1,zs2,ztimp1,zen
  real(fp) :: ze0,zam0,zdth1,zthp,zthtol,zdth2,zth0,zxint,zxip
  real(fp) :: zxion0,zexht,zexhm,zwgt,ztimp2,ztimp,txht,txhm,teht
  real(fp) :: tehm,tpht,tphm,zdximn,zdxisd,zmn,zsd,zdenom,position
  !============
  parameter ( mv1 = 51,        mv2 = mv1)
  real(fp) ::  zdxion(mv1,mv2), zdeion(mv1,mv2), zdpion(mv1,mv2)&
       ,  zrat  (mv1,mv2)&
       ,  ztmfpl(mv1,mv2)&
       ,  zedxt (mv1,mv2), zedxm (mv1,mv2)
  real(fp) ::  zv1(mv1), zv2(mv2)

  real(fp) ::  zwork(mv2)

  real(fp) :: zdxis0w, zdxis1w, zdxis2w, zdximnw, zdxisdw
  real(fp) :: zs0w, zs1w, zs2w, zmnw, zsdw

  character*45  ztitle
  character*14  zstfn      ! file name for stats file

  !      arrays for getting ion energy
  logical  ilreact(mireact)
  real(fp) ::     zereact(mireact)

  logical  ifound

  parameter ( zpi = 3.14159265_fp)

  data        ztmfst / 1._fp/

  !-----------------beginning of executable statements--------------------

  zthst = 0._fp
  kltmfp = 0
  klst   = 0

100 call neupag

  !      opportunity to change d_vparallel mode
  write (lnout(0), '(/'' opportunity to change d_vpar mode:''/)' )
  call ureadi( 'd_vpar mode sw, "nchdvp":$$', nchdvp )

  !      opportunity to change time step
  write (lnout(0), '(/'' opportunity to change time step:''/)' )
  call ureadr8( 'time step "dtn":$$', dtn )


  !      opportunity to change starting theta
  write (lnout(0), '(/'' change starting theta?''/)' )
  call ureadr8( 'starting theta:$$', zthst )

  !      opportunity to turn on writing log file of adjustments
  call ureadi( 'log file of adjustments (0=off, 1=on):$$', kltmfp )
  iltmfp = .false.
  if ( kltmfp .eq. 1 ) iltmfp = .true.

  !      get cut value for statistics on non-extreme errors
  zdxicut = 0.01_fp
  call ureadr8( 'xi cut for non-extreme stats:$$', zdxicut )

  !      opportunity to turn on writing file of statistics on errors
  call ureadi(&
       'write file of error statistics (0=off, 1=on):$$', klst )

  !      get sets of ion parameters and select 1 ion
  if (lxout(0)) write (lnout(0),&
       '('' index  fus. prod.? z   a     e (ev)'')' )
  do isb = 1, nsbeam
     lsbeam = isb
     call inispc

     if ( nlfprod(lsbeam) ) then
        if (nlheavy(lsbeam)) &
             write(lnout(0),*) 'oreplt -- heavy ion??'
        call fuslis(lsbeam,ilreact,zereact)
        il = 0
        ifound = .false.
        do while (.not.ifound)
           il = il + 1
           ifound = ilreact(il)
        end do
        zeion = zereact(il)
     else
        zeion = 1.e5_fp  ! arbitrary choice
     end if

     if (lxout(0)) write(lnout(0),&
          '(i5,l10,f8.0,f4.0,1pe11.2)' )&
          isb, nlfprod(lsbeam), xzbeami, abeam, zeion

  end do

  if (lxout(0)) write (lnout(0), '(/)' )
  call ureadi( 'select ion index (0 to quit):$$', isel )

  if ( isel .le. 0 ) then
     ilmxcr = .false.
     return
  end if

  if ( isel .gt. nsbeam ) then
     call uermsg( ' ?oreplt -- bad species index', ' ', 0 )
     go to 100
  end if

  lsbeam = isel
  call inispc
  if (nlfprod(lsbeam)) then
     call fuslis(lsbeam,ilreact,zereact)
     il = 0
     ifound = .false.
     do while (.not.ifound)
        il = il + 1
        ifound = ilreact(il)
     end do
     zefix = zereact(il)
  else
     zefix = 1.e5_fp  ! arbitrary choice
  end if

  !      choose which variable to hold fixed
  if (lxout(0)) write (lnout(0), 9002 )
9002 format(' choose which variable to hold fixed:'/&
          '      xi (rel flux coord)   1'/&
          '      lambda (vpar/v)       2'/&
          '      energy                3'/ )
  call ureadi( 'index of variable to hold fixed:$$', jfix )

  !      get fixed value and ranges and increments for other variables
  if ( jfix .eq. 1 ) then

     call ureadr8( 'xi:$$', zxifix )
     kxi = 1

     zlamlo = -1._fp
     zlamhi =  1._fp
     zdlam  =  0.04_fp
     call ureadr8( 'lo lambda (vpar/v) :$$', zlamlo )
     call ureadr8( 'hi lambda (vpar/v) :$$', zlamhi )
     call ureadr8( 'increment of lambda:$$', zdlam  )
     klam = (zlamhi-zlamlo)/zdlam + 1.0001_fp
     if (klam.gt.mv1) then
        call uermsg('oreplt: too many steps in lambda',' ',0)
        go to 100
     end if

     zelo   = 1.e4_fp
     zehi   = 5.01e6_fp
     zde    = 1.e5_fp
     call ureadr8( 'lo energy (ev)      :$$', zelo )
     call ureadr8( 'hi energy (ev)      :$$', zehi )
     call ureadr8( 'incr of  energy (ev):$$', zde  )
     ken = (zehi-zelo)/zde + 1.0001_fp
     if (ken.gt.mv1) then
        call uermsg('oreplt: too many steps in energy',' ',0)
        go to 100
     end if

  else if ( jfix .eq. 2 ) then

     call ureadr8( 'lambda:$$', zlamfix )
     klam = 1

     zxilo = -1._fp
     zxihi =  1._fp
     zdxi  =  0.04_fp
     call ureadr8( 'lo xi          :$$', zxilo )
     call ureadr8( 'hi xi          :$$', zxihi )
     call ureadr8( 'increment of xi:$$', zdxi  )
     kxi = (zxihi-zxilo)/zdxi + 1.0001_fp
     if (kxi.gt.mv1) then
        call uermsg('oreplt: too many steps in xi',' ',0)
        go to 100
     end if

     zelo   = 1.e4_fp
     zehi   = 5.01e6_fp
     zde    = 1.e5_fp
     call ureadr8( 'lo energy (ev)      :$$', zelo )
     call ureadr8( 'hi energy (ev)      :$$', zehi )
     call ureadr8( 'incr of  energy (ev):$$', zde  )
     ken = (zehi-zelo)/zde + 1.0001_fp
     if (ken.gt.mv1) then
        call uermsg('oreplt: too many steps in energy',' ',0)
        go to 100
     end if

  else if ( jfix .eq. 3 ) then

     call ureadr8( 'energy (ev):$$', zefix )
     ken = 1

     zxilo = -1._fp
     zxihi =  1._fp
     zdxi  =  0.04_fp
     call ureadr8( 'lo xi          :$$', zxilo )
     call ureadr8( 'hi xi          :$$', zxihi )
     call ureadr8( 'increment of xi:$$', zdxi  )
     kxi = (zxihi-zxilo)/zdxi + 1.0001_fp
     if (kxi.gt.mv1) then
        call uermsg('oreplt: too many steps in xi',' ',0)
        go to 100
     end if

     zlamlo = -1._fp
     zlamhi =  1._fp
     zdlam  =  0.04_fp
     call ureadr8( 'lo lambda (vpar/v) :$$', zlamlo )
     call ureadr8( 'hi lambda (vpar/v) :$$', zlamhi )
     call ureadr8( 'increment of lambda:$$', zdlam  )
     klam = (zlamhi-zlamlo)/zdlam + 1.0001_fp
     if (klam.gt.mv1) then
        call uermsg('oreplt: too many steps in lambda',' ',0)
        go to 100
     end if

  else

     call uermsg( 'oreplt: invalid choice of variable',' ',0)
     go to 100

  end if

  !       get requested no. of midplane crossings
  kcr = 1000
  call ureadi( 'no. times to cross midplane:$$', kcr )

  !      set stop on max no of midplane crossings
  ilmxcr = .true.

  !      save regular time limit and substitute very long one

  zdelt_sv = delt
  delt = 1.e6_fp

  !      turn off perturbations

  iplcx_sv = iplcx
  iplfp_sv = iplfp
  iplcm_sv = iplcm

  iplcx = -1
  iplfp = -1
  iplcm = -1

  !      zero out the statistical sums
  zdxis0 = 0._fp
  zdxis1 = 0._fp
  zdxis2 = 0._fp
  zdxis0w = 0._fp
  zdxis1w = 0._fp
  zdxis2w = 0._fp

  zs0 = 0._fp
  zs1 = 0._fp
  zs2 = 0._fp
  zs0w = 0._fp
  zs1w = 0._fp
  zs2w = 0._fp

  !      zero the count of test orbits
  iorb_loc = 0

  !      do the parameter scans

  call cpu_time( ztimp1 )

  do ixi = 1, kxi

     do ilam = 1, klam

        do ien = 1, ken

           if ( jfix .eq. 1 ) then
              is1 = ilam
              is2 = ien
              zv1(is1) = zlamlo + (ilam-1)*zdlam
              zv2(is2) = zelo   + (ien -1)*zde
              xion  = zxifix
              xksid = zv1(is1)
              zen   = zv2(is2)
           else if ( jfix .eq. 2 ) then
              is1 = ixi
              is2 = ien
              zv1(is1) = zxilo  + (ixi -1)*zdxi
              zv2(is2) = zelo   + (ien -1)*zde
              xion  = zv1(is1)
              xksid = zlamfix
              zen   = zv2(is2)
           else
              is1 = ixi
              is2 = ilam
              zv1(is1) = zxilo  + (ixi -1)*zdxi
              zv2(is2) = zlamlo + (ilam-1)*zdlam
              xion  = zv1(is1)
              xksid = zv2(is2)
              zen   = zefix
           end if

           if ( xion .lt. 0._fp) then
              th = zthst + zpi
           else
              th = zthst
           end if

           !                 make sure th is within normal bounds
150        if ( th .le.  zpi ) goto 155
           th = th - 2._fp*zpi
           goto 150
155        if ( th .gt. -zpi ) goto 160
           th = th + 2._fp*zpi
           goto 155
160        continue

           xion = abs( xion )
           vion  = sqrt( zen / v2toev )
           goose = 1._fp

           ze0    = zen    ! save value before calling orbit

           !                 get starting angular momentum
           call nbmomry( xion, th, 0, rmjion, yion, xbjaco )
           zam0 = rmjion * rmassb * vion

           wght   = 1._fp
           tmfcur = ztmfst
           bperr  = 0._fp
           bpherr = 0._fp

           iorb_loc = iorb_loc + 1

           icr = 0         ! count of midplane crossings
           ind_fake = 1

           call orbit(ind_fake,iorb_loc,ierr)

           if (nbflg.eq.4) then ! completed req. no. xings

              !                   interpolate crossing point
              zdth1 = th - zthp
              if (zdth1 .lt. -zthtol) then
                 zdth1 = zdth1 + 2._fp*zpi
              else if (zdth1 .gt. zthtol) then
                 zdth1 = zdth1 - 2._fp*zpi
              end if

              zdth2 = zth0 - zthp
              if (zdth2 .lt. -zthtol) then
                 zdth2 = zdth2 + 2._fp*zpi
              else if (zdth2 .gt. zthtol) then
                 zdth2 = zdth2 - 2._fp*zpi
              end if

              zxint = zxip + ((xion-zxip)/zdth1) * zdth2

              !                   store errors, time-step magn. factor
              zdxion(is1,is2) = zxint - zxion0
              zdeion(is1,is2) = bperr / ze0
              zdpion(is1,is2) = bpherr / zam0
              ztmfpl(is1,is2) = tmfcur
              zedxt (is1,is2) = zexht
              zedxm (is1,is2) = zexhm

              !                   accumulate sums for xi error statistics
              if ( jfix .gt. 1 ) then

                 zdxis0 = zdxis0 + 1._fp
                 zdxis1 = zdxis1 + zdxion(is1,is2)
                 zdxis2 = zdxis2 + zdxion(is1,is2)**2
                 if (abs(zdxion(is1,is2)).le.zdxicut) then
                    zs0 = zs0 + 1._fp
                    zs1 = zs1 + zdxion(is1,is2)
                    zs2 = zs2 + zdxion(is1,is2)**2
                 end if

                 if ( zxint .le. 1._fp) then
                    j=lcentr+zxion0*nzones
                    zwgt =  rhbs(j,2,lsbeam)
                    zdxis0w = zdxis0w + zwgt
                    zdxis1w = zdxis1w + zwgt * zdxion(is1,is2)
                    zdxis2w = zdxis2w + zwgt * zdxion(is1,is2)**2
                    if (abs(zdxion(is1,is2)).le.zdxicut) then
                       zs0w = zs0w + zwgt
                       zs1w = zs1w + zwgt * zdxion(is1,is2)
                       zs2w = zs2w + zwgt * zdxion(is1,is2)**2
                    end if
                 end if

              end if

           else if (nbflg.eq.1) then ! survived time interval

              write(6, &
                   '(/'' oreplt-- ion did not do req. no. xings''/)')
              call abortr

           else

              zdxion(is1,is2) = 0._fp
              zdeion(is1,is2) = 0._fp
              zdpion(is1,is2) = 0._fp
              ztmfpl(is1,is2) = 0._fp
              zedxt (is1,is2) = 0._fp
              zedxm (is1,is2) = 0._fp

           end if

        end do

     end do

  end do

  call cpu_time( ztimp2 )
  ztimp = ztimp2 - ztimp1

  if (lxout(0)) write ( lnout(0), 9001 )&
       lfbadj, txht, txhm, teht, tehm, tpht, tphm, ztimp

9001 format(/   ' adj. on:      ',l7/&
          ' x t tolerance:',1pe11.2/&
          ' x m tolerance:',1pe11.2/&
          ' e t tolerance:',1pe11.2/&
          ' e m tolerance:',1pe11.2/&
          ' p t tolerance:',1pe11.2/&
          ' p m tolerance:',1pe11.2/&
          ' time for scan:',1pe11.2/)
  
  !      calculate and type out statistics on xi errors
  if ( jfix .gt. 1 ) then
     if ( zdxis0 .gt. 0.5_fp) then
        zdximn  = zdxis1  / zdxis0
        zdxisd  = &
             sqrt( max( (zdxis2  / zdxis0  - zdximn**2), 0._fp) )
     else
        zdximn = 0._fp
        zdxisd = 0._fp
     end if

     if ( zdxis0w .gt. 0.5_fp) then
        zdximnw = zdxis1w / zdxis0w
        zdxisdw = &
             sqrt(max((zdxis2w/zdxis0w - zdximnw**2),0._fp))
     else
        zdximnw = 0._fp
        zdxisdw = 0._fp
     end if

     if ( zs0 .gt. 0.5_fp) then
        zmn  = zs1  / zs0
        zsd  = sqrt( max( (zs2  / zs0  - zmn**2), 0._fp) )
     else
        zmn = 0._fp
        zsd = 0._fp
     end if

     if ( zs0w .gt. 0.5_fp) then
        zmnw = zs1w / zs0w
        zsdw = sqrt(max((zs2w/zs0w - zmnw**2),0._fp))
     else
        zmnw = 0._fp
        zsdw = 0._fp
     end if

     if (lxout(0)) write ( lnout(0), 9003 )&
          icr, zdxis0 , zdximn , zdxisd ,&
          zdxis0w, zdximnw, zdxisdw
9003 format(/,&
          1x,'statistics on xi errors (after',i6,' crossings)'/&
          1x,'           number  mean del xi     std dev  '/&
          1x,'unwgtd',  f12.0,  2f12.4/&
          1x,'  wgtd',1pe12.2,0p2f12.4/)
     
     i3sd = 0
     do ixi = 1, kxi
        do ilam = 1, klam
           do ien = 1, ken
              if ( jfix .eq. 2 ) then
                 is1 = ixi
                 is2 = ien
              else
                 is1 = ixi
                 is2 = ilam
              end if
              if (abs(zdxion(is1,is2)-zdximn).gt. &
                   (3._fp*zdxisd)) then
                 i3sd = i3sd + 1
              end if
           end do
        end do
     end do

     if (lxout(0)) write ( lnout(0), 9004 ) i3sd
9004 format(/,&
          1x,'number outside 3 std dev (unwgtd dist)',i5/)

     if (lxout(0)) write ( lnout(0), 9010 )&
          icr, zdxicut, zs0 , zmn , zsd ,&
          zs0w, zmnw, zsdw
9010 format(/,&
          1x,'statistics on xi errors (after',i6,' crossings)'/&
          1x,'  ** extreme (xi err>',f7.4,') points excluded **'/&
          1x,'           number  mean del xi     std dev  '/&
          1x,'unwgtd',  f12.0,  2f12.4/&
          1x,'  wgtd',1pe12.2,0p2f12.4/)

  end if

  ilst = 0                  ! count of writing stats file

  !      display results

200 if (lxout(0)) write( lnout(0), 9005 )
9005 format(' plot options:'/&
          '   x   error in xi'/&
          '   e   error in energy'/&
          '   p   error in momentum'/&
          '   r   ratio of xi to momentum errors (only for fixed e)'/&
          '   t   time-step magnification factor'/&
          '   y   orbit total est of xi err from p err (fixed e)'/&
          '   z   orbit max   est of xi err from p err (fixed e)'/&
          '   q   quit plotting'/     )

  call uread( 'enter plot choice:$' )

  if ( lchar('q') ) then
     go to 300
  else if ( lchar('x') ) then
     kplot = 1
  else if ( lchar('e') ) then
     kplot = 2
  else if ( lchar('p') ) then
     kplot = 3
  else if ( lchar('r') .and. jfix.eq.3 ) then
     kplot = 4
  else if ( lchar('t') ) then
     kplot = 5
  else if ( lchar('y') ) then
     kplot = 6
  else if ( lchar('z') ) then
     kplot = 7
  else
     call uermsg( ' ?oreplt -- bad plot code', ' ', 0 )
     go to 200
  end if

  if ( jfix .eq. 1 ) then
     write ( ztitle, 9006 )&
          xzbeami, abeam, kcr, nchdvp, zthst, dtn, zxifix
9006 format(2f2.0,i5,i2,f6.2,1x,1pe8.2,' xi=',0pf6.3)

     if ( kplot .eq. 1 ) then
        call r8_grf3f1( zv1, zdxion, zv2, zwork, klam, ken, mv1,&
             'lambda (vpar/v)', 'xi error', 'energy',&
             'vpar/v', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 2 ) then
        call r8_grf3f1( zv1, zdeion, zv2, zwork, klam, ken, mv1,&
             'lambda (vpar/v)', 'frac. energy error', 'energy',&
             'vpar/v', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 3 ) then
        call r8_grf3f1( zv1, zdpion, zv2, zwork, klam, ken, mv1,&
             'lambda (vpar/v)', 'frac. momentum error', 'energy',&
             'vpar/v', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 5 ) then
        call r8_grf3f1( zv1, ztmfpl, zv2, zwork, klam, ken, mv1,&
             'lambda (vpar/v)', 'timestep magn. factor', 'energy',&
             'vpar/v', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     end if

  else if ( jfix .eq. 2 ) then
     write ( ztitle, 9007 )&
          xzbeami, abeam, kcr, nchdvp, zthst, dtn, zlamfix
9007 format(2f2.0,i5,i2,f6.2,1x,1pe8.2,' vp/v=',0pf5.2)

     if ( kplot .eq. 1 ) then
        call r8_grf3f1( zv1, zdxion, zv2, zwork, kxi, ken, mv1,&
             'xi', 'xi error', 'energy',&
             'xi', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 2 ) then
        call r8_grf3f1( zv1, zdeion, zv2, zwork, kxi, ken, mv1,&
             'xi', 'frac. energy error', 'energy',&
             'xi', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 3) then
        call r8_grf3f1( zv1, zdpion, zv2, zwork, kxi, ken, mv1,&
             'xi', 'frac. momentum error', 'energy',&
             'xi', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 5) then
        call r8_grf3f1( zv1, ztmfpl, zv2, zwork, kxi, ken, mv1,&
             'xi', 'timestep magn. factor', 'energy',&
             'xi', ' ', 'ev',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     end if

  else
     write ( ztitle, 9008 )&
          xzbeami, abeam, kcr, nchdvp, zthst, dtn, zefix
9008 format(2f2.0,i5,i2,f6.2,1x,1pe8.2,' e=',1pe8.2)

     if ( kplot .eq. 1 ) then
        call r8_grf3f1( zv1, zdxion, zv2, zwork, kxi, klam, mv1,&
             'xi', 'xi error', 'lambda (vpar/v)',&
             'xi', ' ', 'vpar/v',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 2 ) then
        call r8_grf3f1( zv1, zdeion, zv2, zwork, kxi, klam, mv1,&
             'xi', 'frac. energy error', 'lambda (vpar/v)',&
             'xi', ' ', 'vpar/v',&
             ztitle, 'run '//runid(1:lrunid), 0 )
        
     else if ( kplot .eq. 3 ) then
        call r8_grf3f1( zv1, zdpion, zv2, zwork, kxi, klam, mv1,&
             'xi', 'frac. momentum error', 'lambda (vpar/v)',&
             'xi', ' ', 'vpar/v',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 4 ) then
        do ixi = 1, kxi
           do ilam = 1, klam
              if ( zdpion(ixi,ilam) .ge. 0._fp) then
                 zdenom = max( 1.e-3_fp, zdpion(ixi,ilam) )
              else
                 zdenom = min(-1.e-3_fp, zdpion(ixi,ilam) )
              end if
              zrat(ixi,ilam) = zdxion(ixi,ilam) / zdenom
           end do
        end do
        call r8_grf3f1( zv1, zrat, zv2, zwork, kxi, klam, mv1,&
             'xi', 'ratio of xi to p errors', 'lambda (vpar/v)',&
             'xi', ' ', 'vpar/v',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 5 ) then
        call r8_grf3f1( zv1, ztmfpl, zv2, zwork, kxi, klam, mv1,&
             'xi', 'timestep magn. factor', 'lambda (vpar/ev)',&
             'xi', ' ', 'vpar/v',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 6 ) then
        call r8_grf3f1( zv1, zedxt , zv2, zwork, kxi, klam, mv1,&
             'xi', 't xi err est from p', 'lambda (vpar/v)',&
             'xi', ' ', 'vpar/v',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     else if ( kplot .eq. 7 ) then
        call r8_grf3f1( zv1, zedxm , zv2, zwork, kxi, klam, mv1,&
             'xi', 'm xi err est from p', 'lambda (vpar/v)',&
             'xi', ' ', 'vpar/v',&
             ztitle, 'run '//runid(1:lrunid), 0 )

     end if

  end if

  !           write out file of stats
  if ( jfix.gt.1 .and. klst.eq.1 .and. ilst.eq.0 ) then
     zstfn = runid(1:lrunid)//'st.dat'
#if __f90 || __rs6000
     open ( unit    = 21&
          ,file    = zstfn&
          ,status  ='unknown'&
          ,position='append' )
#else
     open ( unit    = 21&
          ,file    = zstfn&
          ,status  ='unknown'&
          ,access  ='append'&           ! - __f90
          ,carriagecontrol='fortran' )   ! - __f90
#endif
     write ( 21, '(1x,''run '',a10/)' )  runid(1:lrunid)
     write ( 21, '(1x,a45/)' )  ztitle
     write ( 21, 9001 )&
          lfbadj, txht, txhm, teht, tehm, tpht, tphm, ztimp
     write ( 21, 9003 )&
          icr, zdxis0 , zdximn , zdxisd ,&
          zdxis0w, zdximnw, zdxisdw
     write ( 21, 9004 ) i3sd
     write ( 21, 9010 )&
          icr, zdxicut, zs0 , zmn , zsd ,&
          zs0w, zmnw, zsdw
     close ( 21 )
     ilst = ilst + 1
     if (lxout(0)) write ( lnout(0), 9009 )
9009 format (//1x,'*** error statistics file written ***'//)
  end if

  go to 200

  !      loop back to do again

300 delt = zdelt_sv

  iplcx = iplcx_sv
  iplfp = iplfp_sv
  iplcm = iplcm_sv

  go to 100
#else
  write(6,*) ' %nbdebug: oreplt has been disabled.'
#endif

end subroutine oreplt
