program nubeam_test
  
  use rfxqlo_interface
  use nbi_types
  use ezcdf

  !
  !  standalone test driver for NUBEAM
  !  required input data:
  !    1. XPLASMA state file (NetCDF)  ...  written from trfpi/fimain.for:
  !       name of form <runid>_xplasma_nubeam.cdf
  !    2. NUBEAM state file  ...  written from nubeam/nbstart.for:
  !       name of form <runid>_nbi_start.cdf
  !
  !  nb TRANSP users: do `nubeam_test R <runid> A' in the <runid>_fi 
  !                   subdirectory of the directory containing the
  !                   crashed TRANSP run.  This will reproduce the
  !                   MC timestep in this standalone code, with
  !                   extra plotting options.
  !
  ! mg March 2009
  ! Modified  to allow restart and debug on the data set of MPI version of 
  ! NUBEAM as a serial nubeam_test calculation on a sublist of given <myid#>,
  ! which failed, and <numproc#>
  !
  ! bugfix DMC May 2011: when prefix for files is given as <runid>_debug,
  ! extract <runid> rather than using the whole string

  implicit NONE

  type(nbitype_grid) :: zgrid
  type(nbitype_sys) :: zsys
  type(rfxqlo_profile) :: zprof
  character*135 file_xplasma,file_nubeam_fld,file_nubeam_ptcl, &
       file_nubeam_ptcl_local
  integer ifile_xplasma,ifile_nubeam_fld,ifile_nubeam_ptcl,ilen
  character*100 file_path,xplasma_suffix,nubeam_fld_suffix,nubeam_ptcl_suffix
  character*40 runid,runid_pre
  integer lnout
  logical lxout,lchar,lyesc,debug_mpi,nltest_output,key_myidd_0

  integer ierr,ierct,luntrm,lunwrk,id_mcgrid,iorbtot

  integer :: mpierr,idnum
  integer :: istat

  real*8 :: tstart, tend, time0, time1, time2, time3, time4
  real*8 :: time5, time6, time7, time8, time9
  integer :: iorbrzv_option,iorbrzv_rf_option
  !-------------------------------------------------------------------------
  !  RFXQLO data

  !  ***optional***
  !  if filenames are left blank, then, no RF data.

  !  RF data, in case the Monte Carlo RF Operator is active

  character*240 :: rf_idata   ! file containing data that was input to RF
  !  wave field code; $NUBEAM_WORKPATH/trim(rf_idata)

  character*240 :: rf_odata   ! file containing data that was output of RF
  !  wave field code; $NUBEAM_WORKPATH/trim(rf_odata)


#ifndef __MPI
  integer :: icdf, iwarn
                     ! NetCDF channel id
#endif
#ifdef __MPI
  include 'mpif.h'
#endif
  integer :: myid,numprocs,iert,idcod

  !
  !  init filename components
  !

  time0=0.0
  time1=0.0
  time2=0.0
  time3=0.0
  time4=0.0
  time5=0.0
  time6=0.0
  time7=0.0
  time8=0.0
  time9=0.0
  luntrm = lnout(0)

#ifdef __MPI
  call MPI_INIT(mpierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD,myid,mpierr)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,mpierr)
  print *,' '
  print *,'MPI Mode of operation'
  print *,' '
#else
  myid=0
  numprocs=1
  print *,' '
  print *,'NON-MPI Mode of operation'
  print *,' '
#endif

#ifdef __MPI
  tstart=MPI_WTIME()
#else
  tstart=0
#endif


! print *,' '
! print *, 'Process ',myid,' of ',numprocs,' processes alive'
! print *,' '

  runid_pre = '<prefix>'
  file_path = ' '
  xplasma_suffix = '_xplasma_nubeam.cdf' ! see trfpi/fimain.for
  nubeam_fld_suffix = '_nbi_fld_state.cdf' ! see nubeam/nbstart.for
  nubeam_ptcl_suffix = '_nbi_ptcl_state_cpu' ! see nubeam/nbstart.for

  !
  !  init ureadsub
  !

#ifdef __MPI
#else
   call uinitg("NUBEAM_TEST")
#endif


  !  init misc.

  lunwrk = 99   ! fortran i/o unit no.

#ifdef __MPI
#else
  do
#endif

   !
   ! get filenames

     call mkfilenames(file_path, runid_pre, &
          xplasma_suffix, nubeam_fld_suffix, nubeam_ptcl_suffix,   &
          file_xplasma, file_nubeam_fld, file_nubeam_ptcl)

     ifile_xplasma = len_trim(file_xplasma)
     ifile_nubeam_fld = len_trim(file_nubeam_fld)
     ifile_nubeam_ptcl = len_trim(file_nubeam_ptcl)

#ifdef __MPI
#else
          ! menu

     if(lxout(0)) write(luntrm,1001) file_xplasma(1:ifile_xplasma), &
       file_nubeam_fld(1:ifile_nubeam_fld), file_nubeam_ptcl(1:ifile_nubeam_ptcl), &
       trim(runid_pre)
     

1001 format(/ &
          ' XPLASMA state file: "',a,'"'/ &
          ' NUBEAM state file (fld): "',a,'"'// &
          '   to enter "A" successfully, the above must match files that'/ &
          '   exist in the current working directory.'/ &
          ' nubeam_test options:'/ &
          ' P -- reset file path (directory or path)'/ &
          ' R -- reset <runid> or <runid>_debug prefix (currently: "',a,'")'/ &
          ' X -- reset XPLASMA filename suffix <prefix><suffix>'/ &
          ' M -- reset NUBEAM state file (fld) suffix <prefix><suffix>'// &
          ' N -- reset NUBEAM ptcl.state file'/&
          '      <path>/<runid><presuffix><myid#>_<#numproc>.CDF'// &
          ' A -- accept filenames & proceed to read data'/ &
          ' Q -- quit (exit program)'/)

     call uread('nubeam_test: file option:$')


     if(lchar('A')) then

        ! execute

        continue

     else if(lchar('Q')) then

        ! quit

        exit

     else if(lchar('P')) then

        if(lxout(0)) write(luntrm,1002)

        call ureadLw( &
             'directory path for files; "NONE" to use current directory:$', &
             file_path)
        if(file_path.eq.'NONE') file_path=' '

        cycle
     else if(lchar('R')) then


        if(lxout(0)) write(luntrm,1002)

        call ureadLw( &
             'runid:$', runid_pre)

        cycle
     else if(lchar('X')) then

        if(lxout(0)) write(luntrm,1002)

        call ureadLw( &
             'XPLASMA state filename suffix:$', xplasma_suffix)

        cycle
     else if(lchar('M')) then

        if(lxout(0)) write(luntrm,1002)

        call ureadLw( &
             'NUBEAM state (fld) filename suffix:$', nubeam_fld_suffix)

        cycle
     else if(lchar('N')) then


        call ureadLw( &
             'NUBEAM state (ptcl) presuffix only:$', nubeam_ptcl_suffix)

        cycle
     else
        call uermsg(' ?nubeam_test: invalid option code.',' ',0)
        cycle
     endif
     call nubeam_workpath_reset(trim(file_path))

1002 format(/ &
          ' filename syntax:  <path>/<runid><suffix>, where:'/ &
          '   <path> is a directory path; blank means current directory.'/ &
          '   <runid> is a RUNID alphanumeric character string'/ &
          '   <suffix> is the rest of the filename.'// &
          '     example:  path = ".", runid = "12345A06", suffix = ', &
          ' "_xplasma_nubeam.cdf"'/ &
          '     ==> filename = "./12345A06_xplasma_nubeam.cdf"'/&
          '  Please note:'/&
          '    ptcl.filename, as a mandatory imbed CPU# and PROCESSORS#:'/&
          '   <path>/<runid><presuffix><myid#>_<#numproc>.CDF'/&
          '   <presuffix> is everythig between <runid> and <myid#>_<#numproc>.cdf'// &
          '    example:  path = ".", runid = "12345A06"'/&
          '   <presuffix> = "_nubeam_ptcl_state_cpu"'/&
          '    ==> ptcl filename = '/&
          '   "./12345A06_nubeam_ptcl_state_cpu<myid#>_<numproc#>.cdf"'//)
 

#endif
     
     !--------------------------------------------------------
     !  read xplasma data
     !

     ierct=0

     call eq_restore(file_xplasma,ierr)
     if(ierr.ne.0) then
        write(luntrm,*)' ?nubeam_test:  eq_restore failure on:  ',file_xplasma
        call uclear
        ierct = ierct + 1
     endif

#ifdef __MPI
! call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
  time5=MPI_WTIME()
#else
#endif

     !--------------------------------------------------------
     !  read NUBEAM state data
     !



! read field state file from <workpath>/<runid>_<file_nubeam_fld>
! as it's done in nbi_states. Cannot use call nbi_states(..)
! have to retrieve value of debug_mpi and nltest_output
! to ensure that data for debugging had been written
! and have to set myidd and nuprocd to nbi_com first.
!mg 2009
     call nbi_fld_state(file_nubeam_fld,'R',ierr) !only myidd=0 for MPI and serial
                                                  !can read fld stae file
     if(ierr.ne.0) then
        write(luntrm,*)' ?nubeam_test:  nbi_fld_state read failure on:  ',trim(file_nubeam_fld)
        call uclear
        ierct = ierct + 1
     endif

! override state file with local value for: runid
! DMC: this is needed to make source/nubeam_test test procedure work
!      if this is incompatibile with some other nubeam_test use case
!      some new logic may be needed...

! note DMC May 2011: strip off "_debug" suffix if it is detected 
! in runid string

     runid = runid_pre
     ilen=len_trim(runid)
     if(ilen.gt.6) then
        if(runid(ilen-5:ilen).eq.'_debug') then
           runid(ilen-5:ilen)=' '
           ilen = ilen-6
        endif
     endif

     call nubeam_runid_set(runid)


!set RZ -orbit integrator, RF MC operater options

     call nubeam_rf_retrieve(iorbrzv_option,iorbrzv_rf_option)
#ifdef __MPI
#else

     if(iorbrzv_option.eq.1.and.iorbrzv_rf_option.eq.1)then
        !RF data files  need to be provided
        rf_idata=trim(runid)//'_rf_idata.cdf'
        rf_odata=trim(runid)//'_rf_odata.cdf'
        
        if(lxout(0)) write(luntrm,1003) rf_idata,rf_odata
     

1003    format(/ &
          ' file containing wavecode output data: "',a,'"'/ &
          ' file containing wavecode output xplasma state: "',a,'"'/ &
!          ' RF output file: "',a,'"'/ &
!          ' RF output file xplasma state: "',a,'"'/ &
          '   to enter "A" successfully, the above must match files that'/ &
          '   exist in the current working directory.'/ &
          ' nubeam_test options:'/ &
          ' I -- reset wavecode file names'/ &
          ' O -- reset RF output file'/ &
          ' A -- accept filenames & proceed to read data'/ &
          ' Q -- quit (exit program)'/)
        call uread('nubeam_test: file option:$')
             

        if(lchar('A')) then
           
           ! execute
           
           continue
           
        else if(lchar('Q')) then
           
           ! quit
           
           exit
           
        else if(lchar('I')) then
           
           call ureadLw( &
                'file containing wavecode input data:$', &
                rf_idata)
           if(rf_idata.eq.'NONE')then
              rf_idata=' '
              iorbrzv_rf_option=0
              iorbrzv_option=0
           endif
           cycle
        else if(lchar('O')) then
           
           call ureadLw( &
                'file containing wavecode output data:$', &
                rf_odata)
           if(rf_odata.eq.'NONE')then
              rf_odata=' '
              iorbrzv_rf_option=0
              iorbrzv_option=0
           endif
           cycle
        else
           call uermsg(' ?nubeam_test: invalid option code.',' ',0)
           cycle
        endif
        
     endif

     if(iorbrzv_option.eq.1.and.iorbrzv_rf_option.eq.1)then !could be changed, see above...
          ! *** insert xplasma filename in rfxqlo data
        
      
        call rfxqlo_inputfile_read_nbtest(trim(rf_idata),ierr)
        zprof%xplasma_file = file_xplasma
        call rfxqlo_profile_set(zprof,ierr)
        if(ierr.ne.0) then
           write(*,*) &
                ' ?nubeam_test: rfxqlo data recovery error.'
           exit
        endif
        call rfxqlo_outputfile_read(trim(rf_odata),ierr)
         if(ierr.ne.0) then
            write(*,*)' ?nubeam_test:  rfxqlo outputs file read failure:',&
                 &trim(rf_odata)
         endif

     endif
#endif     
#ifdef __MPI
     call nbi_ptcl_read(nubeam_ptcl_suffix)
#else
     if(ierct.gt.0) cycle
#endif
           
     call nubeam_retrieve_debug_mpi(debug_mpi,nltest_output)!retrieve value of debug_mpi
                                                                  !and nltest_output
                                                                  !they are stored in nbi_com
#ifdef __MPI
#else
     if(debug_mpi.and.nltest_output) then
        write(luntrm,2011)
2011          format(/&
                   &' nubeam_test: DEBUG MPI = .TRUE.'/&
                   &'              only sublist will be read'/&
                   &'              example: if ptcl. sublist for myid = 2, from the mpi run'/&
                   &'              with 8 processors will be debugged then enter for '/&
                   &'              number of processors - 8 (last digit of the filename suffix)'/&
                   &'              and for ID number -  2 (from 0 to numprocs-1)'/&
                   &'              enter number processors from the run to be debugged'/&
                   &'              enter sublist ID number which will to be debugged'/)
        call uread('nubeam_test:  enter integer for # of processors:$')
        numprocs=idcod(iert)
22      ierr=0
        call uread('nubeam_test:  enter integer (from 0 to numprocs-1)) for sublist ID# :$')
        myid=idcod(iert)
        if (myid.lt.0.or.myid.gt.numprocs-1) then
           write(luntrm,*)' ?nubeam_test: wrong myid number'
           ierr=1
        endif
        if(ierr.ne.0) goto 22
        call nubeam_debug_set(myid,numprocs)    !set myid and numprocs
                                                !of data which will be 
                                                !debugged
     endif
     call nubeam_workpath_reset(trim(file_path))    ! filepath
                                                     ! of data which will be 
                                                     !debugged
     call nbi_ptcl_read(nubeam_ptcl_suffix)
          
        
     if(ierct.gt.0) cycle
#endif

#ifdef __MPI
! call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
  time6=MPI_WTIME()
#else
#endif

     call nbi_id_restore     ! restore/validate xplasma ids

     !  this cobble tells NUBEAM to regen the MC grid frmt(/ &
     !  in its state file.  It also has the side effect of preventing
     !  state file information from being written by NBSTART.

     call nbi_get_grid(zgrid)
     zgrid%id_mcgrid = -1
     call nbi_update_grid(zgrid)

#ifdef __MPI
  time7=MPI_WTIME()
#else
#endif

     call nbi_interp_profiles(ierr)  ! interpolate profiles to NUBEAM grid
     if(ierr.ne.0) then
        write(luntrm,*) ' ?nubeam_test: NBI_INTERP_PROFILES error.'
#ifdef __MPI
#else
        cycle
#endif
     endif

#ifdef __MPI
  time8=MPI_WTIME()
#else
     do 
              if(lxout(0)) write(luntrm,1009)
1009    format(/' options:'/ &
             ' X: execute beam timestep'/ &
             ' D: debug plotting options'/ &
             ' Q: quit.'/)

        call uread('nubeam_test: execute option:$')

        if(lchar('X')) then
#endif

    ! (re)execute NUBEAM timestep

    !--------------------------------------------------------
    !  ascii printout of NUBEAM control inputs...
 
         call nbi_ascii
 
    !--------------------------------------------------------
    !  perform the NUBEAM calculation

#ifdef __MPI
 time0=MPI_WTIME()
#else
#endif

       print *,' '
       print *, 'Before calling nbstart'
       print *,' '
     call nbstart(ierr)       ! initialize step
       print *,' '
       print *, 'After calling nbstart'
       print *,' '
       if(ierr.ne.0) then
        write(luntrm,*) ' ?nubeam_test:  NBSTART gave ierr=',ierr
       endif

#ifdef __MPI
  time1=MPI_WTIME()
#else
#endif

       write(0,'(a1)')' '
       write(0,'(a6,i6)')'myid= ',myid
       write(0,'(a24,e23.15)') 'Time through nbstart= ',time1-tstart
       write(0,'(a1)')' '

       print *,' '
       print *, 'Before calling depall'
       print *,' '

     call depall(ierr)         ! deposition
       print *,' '
       print *, 'After calling depall'
       print *,' '
       if(ierr.ne.0) then
        write(luntrm,*) ' ?nubeam_test:  DEPALL gave ierr=',ierr
       endif

#ifdef __MPI
 time2=MPI_WTIME()
#else
#endif

       write(0,'(a1)')' '
       write(0,'(a6,i6)')'myid= ',myid
       write(0,'(a17,e23.15)') 'Time for depall= ',time2-time1
       write(0,'(a1)')' '

       key_myidd_0=.FALSE.
       if(myid.eq.0) key_myidd_0=.TRUE.
#ifndef __MPI  
!to allow read serial code when debugging MPI with myidd >0
       key_myidd_0=.TRUE.
#endif
       if(key_myidd_0) then
          print *,' '
          print *, 'Before calling orball'
          print *,' '
          call orball(iorbtot,ierr) ! orbiting
          print *,' '
          print *, 'After calling orball'
          print *,' '
          if(ierr.ne.0) then
             write(luntrm,*) ' ?nubeam_test:  ORBALL gave ierr=',ierr
          endif
       endif

#ifdef __MPI
 time3=MPI_WTIME()
#else
#endif

       write(0,'(a1)')' '
       write(0,'(a6,i6)')'myid= ',myid
       write(0,'(a17,e23.15)') 'Time for orball= ',time3-time2
       write(0,'(a1)')' '

       key_myidd_0=.FALSE.
       if(myid.eq.0) key_myidd_0=.TRUE.
#ifndef __MPI  
!to allow read serial code when debugging MPI with myidd >0
       key_myidd_0=.TRUE.
#endif
       if(key_myidd_0) then
          write(0,'(a1)')' '
          write(0,'(a6,i6)')'myid= ',myid
          write(0,'(a15)') 'Before nbfinish'
          write(0,'(a1)')' '

     call nbfinish(ierr)       ! complete the step
     if(ierr.ne.0) then
      write(luntrm,*) ' ?nubeam_test:  NBFINISH gave ierr=',ierr
     endif

       write(0,'(a1)')' '
       write(0,'(a6,i6)')'myid= ',myid
       write(0,'(a14)') 'After nbfinish'
       write(0,'(a1)')' '

!--------------------------------------------------------
!  ascii printout of scalar outputs...
 
     call nbo_ascii

!--------------------------------------------------------
!  plot output...

!       if (myid.eq.0) then
 
     call eqdbg_plot

   endif

#ifdef __MPI
#else
     exit

  else if(lchar('D')) then

           !--------------------------------------------------------
           !  setup NUBEAM calculation (deposition only)

           call nbstart(ierr)       ! initialize step
           if(ierr.ne.0) then
              if(ierr.ne.111) then
                 write(luntrm,*) ' ?nubeam_test:  NBSTART gave ierr=',ierr
                 exit
              endif
           endif

           write(luntrm,*) ' %nbstart: successful.'

           do
              if(lxout(0)) write(luntrm,1020)
1020          format(/' debug options:'/ &
                   ' C -- check NUBEAM splines'/ &
                   ' N -- NBLTES (beam deposition plot tests)'/ &
                   ' P -- PLOT_FPREACT call (atomic/nuclear rate data)'/ &
                   ' X -- XPLASMA debug plots'/ &
                   ' S -- XSTRALN test plots'/ &
                   ' Q -- quit'/)

              call uread('nubeam_test(debug): option:$')

              if(lchar('C')) then
                 call chkspl
              else if(lchar('N')) then
                 call nbltes
              else if(lchar('P')) then
                 call plot_fpreact
              else if(lchar('X')) then
                 call eqdbg_plot
              else if(lchar('S')) then
                 call xstest(runid)
              else if(lchar('Q')) then
                 exit
              else
                 call uermsg('?nubeam_test(debug): invalid option.',' ',0)
              endif
           enddo

           exit


        else if(lchar('Q')) then
           exit

        else
           call uermsg('nubeam_test: invalid "execute" option.',' ',0)
        endif

     enddo

#ifdef __MPI
#else
  enddo
#endif


  call uexita(0)

#endif

#ifdef __MPI
 tend=MPI_WTIME()

      print*, ' '
      print*, 'myid= ',myid
      print*, 'Time through nbi_ascii= ',time0-tstart
      print*, '('
      print*, 'Time through eq_restore= ',time5-tstart
      print*, 'Time for nbi_states= ',time6-time5
      print*, 'Time for nbi_update_grid= ',time7-time6
      print*, 'Time for nbi_intrp_profiles"= ',time8-time7
      print*, 'Time for nbi_ascii"= ',time0-time8
      print*, ')'
      print*, ' '
      print*, 'Time for nbstart= ',time1-time0
      print*, 'Time for depall= ',time2-time1
      print*, 'Time for orball= ',time3-time2
      print*, 'Time to finish up= ',tend-time3
      print*, 'Total time= ',tend-tstart
      print*, ' '

  call MPI_FINALIZE(mpierr)

#endif

end program nubeam_test

subroutine mkfilenames(path,runid,suffix1,suffix2,suffix3,file1,file2,file3)

  character*(*), intent(in) :: path     ! directory path (can be blank)
  character*(*), intent(in) :: runid    ! common filename prefix

  ! filename suffixes:
  character*(*), intent(in) :: suffix1,suffix2,suffix3

  ! resulting filenames:
  character*(*), intent(out) :: file1,file2,file3

  !----------------------------------

  integer ilpath,ilpre,ilsuf1,ilsuf2,ilsuf3
  integer luntrm,lnout

  !----------------------------------

  ilpath = len_trim(path)
  ilpre = len_trim(runid)
  ilsuf1 = len_trim(suffix1)
  ilsuf2 = len_trim(suffix2)
  ilsuf3 = len_trim(suffix3)

  if(min(ilpre,ilsuf1,ilsuf2).eq.0) then

     file1=' '; file2 = ' '
     luntrm = lnout(0)
     write(luntrm,*) ' ? blank filename component detected:'
     write(luntrm,*) '   runid = "',runid(1:max(1,ilpre)),'"'
     write(luntrm,*) '   suffix1 = "',suffix1(1:max(1,ilsuf1)),'"'
     write(luntrm,*) '   suffix2 = "',suffix2(1:max(1,ilsuf2)),'"'
     write(luntrm,*) '   suffix3 = "',suffix3(1:max(1,ilsuf3)),'"'

  else if(ilpath.eq.0) then

     file1 = runid(1:ilpre)//suffix1(1:ilsuf1)
     file2 = runid(1:ilpre)//suffix2(1:ilsuf2)
     file3 = suffix3(1:ilsuf3)

  else

     if(path(ilpath:ilpath).eq.'/') then

        file1 = path(1:ilpath)//runid(1:ilpre)//suffix1(1:ilsuf1)
        file2 = path(1:ilpath)//runid(1:ilpre)//suffix2(1:ilsuf2)
        file3 = suffix3(1:ilsuf3)

     else

        file1 = path(1:ilpath)//'/'//runid(1:ilpre)//suffix1(1:ilsuf1)
        file2 = path(1:ilpath)//'/'//runid(1:ilpre)//suffix2(1:ilsuf2)
        file3 = suffix3(1:ilsuf3)

     endif

  endif

end subroutine mkfilenames


