MODULE basic
!
  USE hashtable
  USE constants
  USE bsplines
  USE mumps_bsplines
  USE futils
  USE mpihelper
  IMPLICIT NONE
!
!   Basic module for time dependent problems
!
  CHARACTER(len=128) :: label1, label2, label3, label4
!
!   BASIC Namelist
!
  LOGICAL          :: nlres = .FALSE.           !< Restart flag
  LOGICAL          :: nlsave = .TRUE.           !< Checkpoint (save) flag
  LOGICAL          :: newres=.FALSE.            !< New result HDF5 file
  LOGICAL          :: nlxg=.FALSE.              !< Show graphical interface Xgrafix
  LOGICAL          :: nlmaxwellsource = .FALSE. !< Activate the maxwell source
  INTEGER          :: nrun=1                    !< Number of time steps to run
  REAL(kind=db) :: job_time=3600.0           !< Time allocated to this job in seconds
  REAL(kind=db) :: tmax=100000.0             !< Maximum simulation time
  REAL(kind=db) :: extra_time=60.0           !< Extra time allocated
  REAL(kind=db) :: dt=1                      !< Time step
  REAL(kind=db) :: time=0                    !< Current simulation time (Init from restart file)
!
!   Other basic global vars and arrays
!
  INTEGER :: jobnum                    !< Job number
  INTEGER :: step                      !< Calculation step of this run
  INTEGER :: cstep=0                   !< Current step number (Init from restart file)
  LOGICAL :: nlend                     !< Signal end of run
  INTEGER :: ierr                      !< Integer used for MPI
  INTEGER :: it0d=1                    !< Number of iterations between 0d values writes to hdf5
  INTEGER :: it2d=100                  !< Number of iterations between 2d values writes to hdf5
  INTEGER :: itparts=1000              !< Number of iterations between particles values writes to hdf5
  INTEGER :: ittext=10                 !< Number of iterations between text outputs in the console
  INTEGER :: itrestart=10000           !< Number of iterations between save of restart.h5 file
  INTEGER :: ittracer=100              !< Number of iterations between save of traced particles position and velocity
  INTEGER :: itcelldiag=100000         !< Number of iterations between save of celldiag diagnostic
  INTEGER :: nbcelldiag=0              !< Number of celldiagnostics
  INTEGER :: itgraph                   !< Number of iterations between graphical interface updates
  INTEGER :: mpirank                   !< MPIrank of the current processus
  INTEGER :: mpisize                   !< Size of the MPI_COMM_WORLD communicator
  INTEGER :: rightproc                 !< Rank of next processor in the z decomposition
  INTEGER :: leftproc                  !< Rank of previous processor in the z decomposition
!
!  List of logical file units
  INTEGER :: lu_in       = 90          !< File duplicated from STDIN
  INTEGER :: lu_stop     = 91          !< stop file, see subroutine TESEND
  INTEGER :: lu_partfile = 120         !< particle loading file, see beam::loadpartfile
!
!  HDF5 file
  CHARACTER(len=256) :: resfile = "results.h5" !< Main result file
  CHARACTER(len=256) :: rstfile = "restart.h5" !< Restart file
  CHARACTER(len=256) :: magnetfile = ""       !< H5 file containing the magnetic field definition
  CHARACTER(len=256) :: partfile(10)=""       !< Particle loading file
  CHARACTER(len=256) :: addedtestspecfile(10)="" !< Particle file list for added particles at restart
  INTEGER           :: fidres                 !< File ID for resfile
  INTEGER           :: fidrst                 !< File ID for restart file
  TYPE(BUFFER_TYPE) :: hbuf0                  !< Hashtable for 0d var
! 
! Plasma parameters
  LOGICAL          :: nlPhis= .TRUE.                 !< Calculate self consistent electric field flag
  LOGICAL          :: nlfreezephi= .FALSE.          !< Freeze the Poisson solver to the field obtained at (re-)start
  LOGICAL          :: nlclassical= .FALSE.           !< If true, solves the equation of motion according to classical 
                                                     !! dynamics
  LOGICAL          :: nlperiod(2)=(/.false.,.false./)!< Set periodic splines on or off
  LOGICAL          :: partperiodic= .TRUE.           !< Sets if the particles boundary conditions are periodic or open
  INTEGER          :: nbaddtestspecies=0             !< On restart number of files to read to add test particles
  INTEGER          :: nplasma                        !< Number of macro-particles on initialisation
  INTEGER          :: nbspecies = 1                  !< Number of particles species also counting tracing particles
  INTEGER          :: npartsalloc = 0                !< Size of particle memory allocated at the begining of the simulation
  INTEGER       :: nblock                         !< Number of slices in Z for stable distribution initialisation
  REAL(kind=db) :: potinn                         !< Electric potential at the inner metallic wall
  REAL(kind=db) :: potout                         !< Electric potential at the outer metallic wall
  REAL(kind=db) :: B0                             !< Max magnitude of magnetic field
  REAL(kind=db), allocatable     :: Bz(:), Br(:)  !< Magnetic field components
  REAL(kind=db), allocatable     :: Athet(:)      !< Theta component of the magnetic vector potential
  TYPE(spline2d), SAVE :: splrz                   !< Spline at r and z for total electric field
  TYPE(spline2d), SAVE :: splrz_ext               !< Spline at r and z for external electric field
  REAL(kind=db), allocatable     :: Ez(:), Er(:)  !< Electric field components ( ext+self )
  REAL(kind=db), allocatable     :: pot(:)        !< Electro static potential ( ext+self )
  REAL(kind=db), allocatable     :: Ezxt(:), Erxt(:)  !< external Electric field components
  REAL(kind=db), allocatable     :: potxt(:)        !< external Electro static potential
  REAL(kind=db)     :: radii(4)                   !< Inner and outer radius of cylinder and radii of fine mesh region
  REAL(kind=db)     :: plasmadim(4)               !< Zmin Zmax Rmin Rmax values for plasma particle loading
  INTEGER           :: distribtype=1              !< Type of distribution function used to load the particles
                                                  !!1: gaussian, 2: Stable as defined in 4.95 of Davidson
  REAL(kind=db)     :: H0=0                       !< Initial value of Hamiltonian for distribution 2
  REAL(kind=db)     :: P0=0                       !< Initial canonical angular momentum for distribution 2
  REAL(kind=db)     :: temprescale = -1.0         !< Factor used for temperature rescaling in case of a restart (<0 -> no rescaling)
  INTEGER           :: samplefactor =-1           !< Factor used for the up-sampling of the particles number 
  REAL(kind=db)     :: lz(2)                      !< Lower and upper cylinder limits in z direction
  REAL(kind=db)     :: n0                         !< Physical plasma density parameter
  REAL(kind=db), DIMENSION(:,:), ALLOCATABLE, SAVE:: moments    !< Moments of the distribution function evaluated every it2d
  REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: rhs    !< right hand side of the poisson equation solver
  REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: volume !< Volume covered by each spline for density calculation
  INTEGER           :: nz                                 !< Number of grid intervals in z
  INTEGER           :: nr                                 !< Total number of grid intervals in r 
  INTEGER           :: nnr(3)                             !< Number of grid intervals in r in each subdomain
  REAL(kind=db)     :: dz                                 !< Cell size in z
  REAL(kind=db)     :: dr(3)                              !< Cell size in r for each region
  REAL(kind=db), ALLOCATABLE :: zgrid(:)                  !< Nodes positions in longitudinal direction
  REAL(kind=db), ALLOCATABLE :: rgrid(:)                  !< Nodes positions in radial direction
  REAL(kind=db)     :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm !< Normalization constants
  REAL(kind=db)     :: qsim                               !< Charge of superparticles
  REAL(kind=db)     :: msim                               !< Mass of superparticles
  REAL(kind=db)     :: partmass=me                        !< Mass of physical particle
  INTEGER           :: femorder(2)                        !< FEM order 
  INTEGER           :: ngauss(2)                          !< Number of gauss points 
  LOGICAL           :: nlppform =.TRUE.                   !< Argument of set_spline
  INTEGER, SAVE     :: nrank(2)                           !< Number of splines in both directions
  REAL(kind=db)     :: omegac                             !< Cyclotronic frequency
  REAL(kind=db)     :: omegap                             !< Plasma frequency
  REAL(kind=db)     :: temp                               !< Initial temperature of plasma
  REAL(kind=db)     :: Rcurv                              !< Magnetic field curvature coefficient
  REAL(kind=db)     :: Width                              !< Distance between two magnetic mirrors
  REAL(kind=db)     :: weights_scale=1.0                  !< Scale factor for the particle weights on restart (only for newres=.true.)
  INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds     !< Index of bounds for local processus in Z direction
  INTEGER           :: bscaling = -1                      !< if >0 rescale the magnetic field read from h5 file before calculating value at grid points, if <0 rescale after interpolation, if = 0 doesn't rescale
  REAL(kind=db):: invdz, invdr(3)
  
CONTAINS
!
!================================================================================
  SUBROUTINE basic_data
!
!   Define basic data
!
    use mpihelper
    IMPLICIT NONE
!
!   Local vars and arrays
    CHARACTER(len=256) :: inputfilename
    INTEGER, EXTERNAL :: OMP_GET_MAX_THREADS
!
    NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg,          &
         &           nplasma, potinn, potout, B0, lz, n0, nz, nnr, femorder, ngauss,             &
         &           nlppform, plasmadim, radii, temp, Rcurv, width, it0d, it2d, itparts, ittext, &
         &           resfile, rstfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0, partperiodic, &
         &           temprescale, samplefactor, nlmaxwellsource, npartsalloc, partfile, partmass, nbspecies, &
         &           ittracer, itcelldiag, nbcelldiag, magnetfile, weights_scale, nlfreezephi, nbaddtestspecies, &
         &           addedtestspecfile, bscaling
!________________________________________________________________________________
!                   1.   Process Standard Input File
!
  IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN
      WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENT REQUIRED, STOPPING'
      STOP
  ENDIF
  
  CALL GET_COMMAND_ARGUMENT(1,inputfilename)
  
  OPEN(UNIT=lu_in,FILE=trim(inputfilename),ACTION='READ')
  
  IF(mpirank .eq. 0) THEN
!________________________________________________________________________________
!                   1.   Label the run
!
    READ(lu_in,'(a)') label1
    READ(lu_in,'(a)') label2
    READ(lu_in,'(a)') label3
    READ(lu_in,'(a)') label4
!
    WRITE(*,'(12x,a/)') label1(1:len_trim(label1))
    WRITE(*,'(12x,a/)') label2(1:len_trim(label2))
    WRITE(*,'(12x,a/)') label3(1:len_trim(label3))
    WRITE(*,'(12x,a/)') label4(1:len_trim(label4))
!________________________________________________________________________________
!                   2.   Read in basic data specific to run
!
    READ(lu_in,basic)
    WRITE(*,basic)
#if _DEBUG==1
    WRITE(*,*) "Compiled in debug mode"
#endif
  ELSE
    READ(lu_in,basic)
  END IF
  ! Distribute run parameters to all MPI workers
  CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation
  WRITE(*,'(a,i4.2,a,i4.2,a)')"Running on ",mpisize," tasks with", omp_get_max_threads() ," openMP threads"

  IF(samplefactor .gt. 1 .and. .not. newres) THEN
    IF(mpirank.eq.0) WRITE(*,*)"To increase the number of particles, you need to create a new result file (set newres to 1)"
    CALL MPI_abort(MPI_COMM_WORLD,-1,ierr)
  END IF

  IF (npartsalloc .lt. nplasma) THEN
    npartsalloc=nplasma
  END IF
!
  END SUBROUTINE basic_data
  !================================================================================
  SUBROUTINE daytim(str)
    !
    !   Print date and time
    !
        IMPLICIT NONE
    !
        CHARACTER(len=*), INTENT(in) :: str
    !
    !   Local vars and arrays
        CHARACTER(len=16) :: d, t, dat, functime
    !________________________________________________________________________________
    !
        CALL DATE_AND_TIME(d,t)
        dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4)
        functime=t(1:2) // ':' // t(3:4) // ':' // t(5:10)
        WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), functime(1:12)
    !
      END SUBROUTINE daytim
!================================================================================
  SUBROUTINE timera(cntrl, str, eltime)
!
!   Timers (cntrl=0/1 to Init/Update)
!
    IMPLICIT NONE
    INTEGER, INTENT(in)                  :: cntrl
    CHARACTER(len=*), INTENT(in)         :: str
    REAL(kind=db), OPTIONAL, INTENT(out) :: eltime
!
    INTEGER, PARAMETER                    :: ncmax=128
    INTEGER, SAVE                         :: icall=0, nc=0
    REAL(kind=db), SAVE                   :: startt0=0.0
    CHARACTER(len=16), SAVE               :: which(ncmax)
    INTEGER                               :: lstr, found, i
    REAL(kind=db)                         :: seconds
    REAL(kind=db), DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0
!________________________________________________________________________________
    IF( icall .EQ. 0 ) THEN
       icall = icall+1
       startt0 = seconds()
    END IF

    lstr = LEN_TRIM(str)
    IF( lstr .GT. 0 ) found = loc(str)
!________________________________________________________________________________
!
    SELECT CASE (cntrl)
!
    CASE(-1)    !  Current wall time
       IF( PRESENT(eltime) ) THEN
          eltime = seconds() - startt0
       ELSE
          WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', seconds() - startt0
       END IF
!
    CASE(0)    !  Init Timer
       IF( found .EQ. 0 ) THEN  !  Called for the 1st time for 'str'
          nc = nc+1
          which(nc) = str(1:lstr)
          found = nc
       END IF
       startt(found) = seconds()
!
    CASE(1)   !  Update timer
       endt(found) = seconds() - startt(found)
       IF( PRESENT(eltime) ) THEN
          eltime = endt(found)
       ELSE
          WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', endt(found)
       END IF
!
    CASE(2)   !  Update and reset timer
       endt(found) = endt(found) + seconds() - startt(found)
       startt(found) = seconds()
       IF( PRESENT(eltime) ) THEN
          eltime = endt(found)
       END IF
!
    CASE(9)   !  Display all timers
       IF( nc .GT. 0 ) THEN
          WRITE(*,'(a)') "Timer Summary"
          WRITE(*,'(a)') "============="
          DO i=1,nc
             WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endt(i)
          END DO
       END IF
!
    END SELECT
!
  CONTAINS
    INTEGER FUNCTION loc(funcstr)
      CHARACTER(len=*), INTENT(in) :: funcstr
      INTEGER :: j, ind
      loc = 0
      DO j=1,nc
          ind = INDEX(which(j), funcstr(1:lstr))
          IF( ind .GT. 0 .AND. LEN_TRIM(which(j)) .EQ. lstr ) THEN
             loc = j
             EXIT
          END IF
       END DO
    END FUNCTION loc
  END SUBROUTINE timera
!================================================================================
END MODULE basic