module def_sectors
! Time-stamp: "2003-10-14 12:51:24 cjn"
  use precisn, only: wp
  use io_units, only: fo
  use error_out, only: error_check
  use mpi_params, only: io_processor
  implicit none

  private
  public reset_sec, sectors
  public nsect, nh, nl, nx, asect, w_sect

  integer, save               :: nsect   ! # sectors
  integer, save               :: nh      ! Hamiltonian dimension (all channels)
  integer, save               :: nl      ! Legendre order
  integer, save               :: nx      ! # Gauss points
  real(wp), allocatable, save :: asect(:)
  real(wp)                    :: w_sect

contains

  subroutine reset_sec
! deallocate memory
    integer            :: status
    deallocate (asect, stat=status)
    call error_check(status, " Error in reset_sec ")
  end subroutine reset_sec

  subroutine sectors (nc)
 ! calculate sector end; rlast lies on the end of last sector
    use energy_grid, only: fine_mesh, ebig
    use rmx1_in, only: nhmax, ncmax, rafin
    use hfile_data, only: rmatr
    integer, intent(in) :: nc    ! # channels
    real(wp)            :: rfirst, rlast ! propagation range
    integer             :: ibuf(4), status, i
    real(wp)            :: dum(1)


       rfirst = rmatr
       rlast = rafin
       if (rafin < rmatr) then
          write (fo,'(a)') 'sectors: propagation range error'
          write (fo,'(a,2f12.6)') 'rmatr, rafin = ', rmatr, rafin
          STOP
       end if

       ! Benchmark version of code allowing easier problem size variation by user
       call bench_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, fine_mesh)

       w_sect = (rlast - rfirst) / REAL(nsect,wp)   ! actual sector width
       allocate (asect(nsect+1), stat=status)
       asect(1) = rfirst
       do i = 1, nsect
          asect(i+1) = asect(1) + i * w_sect
       end do
       nh = nl * nc

! output sector radius
       if(io_processor) then
         write (fo,'(/a,i4)') 'Number of basis functions (Benchmark Case) = ', nl
         write (fo,'(a,i4)') 'Number of abscissae, nx = ', nx
         write (fo,'(a,i8)') 'Hamiltonian dimension = ', nh
         write (fo,'(a,i4)') 'Total Number of Sectors = ', nsect
         write (fo,'(a,f12.6)') 'Sector width   = ', w_sect
         do i = 1, nsect+1
           write (fo,'(i4,5x,a,f12.6)') i, 'asect = ', asect(i)
         end do
         write (fo,'(/)')
         call flush(fo)
       end if
   end subroutine sectors

  subroutine get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
       fine_mesh)
! use PGB formula to determine sector sizes and # of basis functions
! formula:   kmax =  nl * pi   ~=   nl
!                   ---------       --       delta(a) = rlast-rfirst
!                   3*delta(a)    delta(a)
    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
    integer, intent(in)  :: ncmax   ! max # channels for this L
    real(wp), intent(in) :: rfirst ! initial radius
    real(wp), intent(in) :: rlast  ! final radius
    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
    logical, intent(in)  :: fine_mesh ! fine mesh flag
    real(wp)             :: ansect, kmax , anlmax

    kmax = SQRT(ebig)
    if (fine_mesh) then
       anlmax = nhmax / ncmax
       ansect = (rlast - rfirst) * kmax / anlmax
       nsect = MAX(INT(ansect+1.0_wp), 2)
       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 10)
    else ! coarse mesh, therefore set nl to 10
       nl = 10
       ansect = (rlast - rfirst) * kmax / REAL(nl,wp)
       nsect = MAX(INT(ansect+1.0_wp), 2)
    end if
    nx = nl + 5
  end subroutine get_nsect_nl

  subroutine bench_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
       fine_mesh)
    ! **** Benchmark version of routine to allow greater user control of problem size ****
    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
    integer, intent(in)  :: ncmax   ! max # channels for this L
    real(wp), intent(in) :: rfirst ! initial radius
    real(wp), intent(in) :: rlast  ! final radius
    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
    logical, intent(in)  :: fine_mesh ! fine mesh flag
    real(wp)             :: ansect, kmax , anlmax
    CHARACTER(len=255) :: c_nl, c_nsect
    integer :: i_nl, i_nsect, status_val

    kmax = SQRT(ebig)
    if (fine_mesh) then

       nsect = -1
       i_nl = -1
       CALL get_environment_variable("RMX_NL_FINE", c_nl, STATUS=status_val)
       read( c_nl, '(i10)' ) i_nl
       if (i_nl > 0 ) then
          nl = i_nl
       else
          nl = 12
       end if

       CALL get_environment_variable("RMX_NSECT_FINE", c_nsect, STATUS=status_val)
       read( c_nsect, '(i10)' ) i_nsect
       if (i_nsect > 0 ) then
          nsect = i_nsect
       else
          nsect = 5
       end if

       if(io_processor) then
         write (*,'(a,i0)') "BENCHMARK CODE: Number of Sectors in Fine Region = ", nsect
         write (*,'(a,i0)') "BENCHMARK CODE: Number of Basis Functions used in Fine Region = ", nl
       end if

    else ! coarse mesh
       i_nsect=-1
       i_nl=-1
       CALL get_environment_variable("RMX_NL_COARSE", c_nl, STATUS=status_val)
       read( c_nl, '(i10)' ) i_nl
       if (i_nl > 0 ) then
          nl = i_nl
       else
          nl = 6
       end if

       CALL get_environment_variable("RMX_NSECT_COARSE", c_nsect, STATUS=status_val)
       read( c_nsect, '(i10)' ) i_nsect
       if (i_nsect > 0 ) then
          nsect = i_nsect
       else
          nsect = 20
       end if

       if(io_processor) then
         write (*,'(a,i0)') "BENCHMARK CODE: Number of Sectors in Coarse Region = ", nsect
         write (*,'(a,i0)') "BENCHMARK CODE: Number of Basis Functions used in Coarse Region = ", nl
       end if

    end if

    ! number of Gauss points
    nx = nl + 5

  end subroutine bench_get_nsect_nl

end module def_sectors
