module sector_potl
! generate sector potential
! Time-stamp: "2003-09-19 11:52:40 cjn"
  use precisn, only: wp
  use io_units, only: fo
  use potl_cofs, only: cf, lamd ! asymptotic potential coefficients
  use rmx1_in, only: inc_lrp_prop, n_lambda
  implicit none

  integer, save                 :: ion    ! ionicity
  integer, save                 :: lmx    ! # multipoles
  integer, allocatable, save    :: lc(:)  ! channel orbital a.m.
  integer, save                 :: nc
  real(wp), allocatable, save   :: v(:,:,:)

  private
  public potl, potl_pams, reset_potl
  public v, nc

contains

  subroutine potl_pams (nct, nx, ncp)
! set potential parameters
    use hfile_data, only: lchl
!    use blacs, only: p_error
    use scaling, only: get_ionicity
    use error_out, only: error_check
    integer, intent(in)   :: nct  ! # channels in partition
    integer, intent(in)   :: nx   ! # radial points
    integer, intent(in)   :: ncp  ! channel origin
    integer               :: status

    nc = nct
    allocate (v(nx,nc,nc), lc(nc), stat=status)
    call error_check (status, 'potl_pams: allocation')
    lc(1:nc) = lchl(ncp:ncp+nc-1)
    ion = get_ionicity()
    lmx = n_lambda
  end subroutine potl_pams

  subroutine reset_potl
    integer      :: status
    deallocate (v, lc, stat=status)
    if (status /= 0) then
       write (fo,'(a)') 'reset_potl: deallocation error'
       STOP
    end if
  end subroutine reset_potl

  subroutine potl (nc, nx, r)
! potl: potentials expanded in inverse powers of radial distance;
    use rmx1_in, only: bug7, n_lambda
    integer, intent(in)      :: nc              ! no. channels
    integer, intent(in)      :: nx              ! no. basis fns
    real(wp), intent(in)     :: r(nx)            ! radial points
    real(wp)                 :: dg(nx)
    real(wp)                 :: vt
    integer                  :: i, j, k, locc, lbb, par,  lm

    if (inc_lrp_prop) then  ! asymptotic potential included
       lm = MIN(lmx, n_lambda)
       v = 0.0_wp
       do j = 1, nc
          dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
          locc = j*(j-1)/2
          do i = 1, j
             if (i == j) v(:,j,i) = dg
             do k = 1, lm, 2
                lbb = (k + 1) / 2
                v(:,i,j) = v(:,i,j) + cf(locc+i,lbb) &
                     / r**(lamd(locc+i,lbb)+1)
             end do
          end do
       end do
    else
       v = 0.0_wp
       do j = 1, nc
          dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
          v(:,j,j) = dg
       end do
    end if
    if (bug7 > 0) then
!       do i = 1, nc
       do i = 1, 3
!          do j = 1, nc
          do j = 1, 3
             write (fo,'(a,i4,i4)') ' Potential (v) matrix:, col, &
                  &row ',i,j
             write (fo,'(6f14.7)') v(:,i,j)
          end do
       end do
    end if
  end subroutine potl
end module sector_potl
