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, bug7, packed_cf
!  use blacs, only:  iam
  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 error_out, only: error_check
    use scaling, only: get_ionicity
    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
!    use blacs, only: ctxt, io_processor
    use error_out, only: error_check
    integer      :: status
    deallocate (v, lc, stat=status)
    if (status /= 0) then
      call error_check (status, 'reset_potl: allocation')
      STOP 
    end if
  end subroutine reset_potl

  subroutine potl (nc, nx, r)
!    use hfile_data, only: cf_nz, ic_label, ir_label, ncf_nonzero, ncf2_nonzero, spin_v, nc1
    use hfile_data, only: cf_nz, ic_label, ir_label, ncf_nonzero, ncf2_nonzero, spin_v
! potl: potentials expanded in inverse powers of radial distance;
    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)
    integer                  :: ncf_temp(n_lambda), icf_st(n_lambda)
    integer                  :: i, j, k, locc, lbb, lm, icf, kp1

    if (inc_lrp_prop) then  ! asymptotic potential included
       lm = MIN(lmx, n_lambda)
       v = 0.0_wp
       if (packed_cf) then
! we just assume a symmetric matrix
          if (spin_v == 1) then
             icf_st(1:lm) = 1
             ncf_temp(1:lm) = ncf_nonzero(1:lm)
          else
             icf_st(1:lm) = ncf_nonzero(1:lm) + 1
             ncf_temp(1:lm) = ncf2_nonzero(1:lm)
          end if 
          do j = 1, nc
             dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
             v(:,j,j) = dg
          end do
          do k = 1, lm
             kp1 = k + 1
             do icf = icf_st(k), ncf_temp(k)
                j = ic_label(icf,k)
                i = ir_label(icf,k)
                v(:,i,j) = v(:,i,j) + cf_nz(icf,k) / r ** kp1
             end do
          end do
! symmetric matrix
          do k = 1, lm
             do icf = icf_st(k), ncf_temp(k)
                j = ic_label(icf,k)
                i = ir_label(icf,k)
                v(:,j,i) = v(:,i,j)
             end do
          end do  
       else
! atomic method with symmetry taken into account
          do j = 1, nc
             dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
             v(:,j,j) = dg
             locc = j*(j-1)/2
             do i = 1, j
                do k = 1, lm, 2
                   lbb = (k + 1) / 2
!                if (i == 1 .and. j == 3) &
!                  write(fo,'(a,2i3,2f14.7)') &
!                '(1,3): lbb, lamd(lbb, rr(1), cf =', lbb, lamd(locc+i,lbb), &
!                  (1.0_wp/ r(1)), cf(locc+i,lbb)           
                   v(:,i,j) = v(:,i,j) + cf(locc+i,lbb) &
                     / r**(lamd(locc+i,lbb)+1)
                   v(:,j,i) = v(:,i,j)
!                 write(fo,*) 'iam', iam, ' i,j,k, lbb, cf:', i,j,k, lbb, cf(locc+i,lbb)
                 end do
             end do
          end do
       end if
    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
!             write (fo,'(a)') ' radial values '
!             write (fo,'(6f14.7)') r(:)
!       do i = 1, nc
       do i = 1, MIN(3,nc)
          do j = 1, nc
!          do j = 1, MIN(3,nc)
!             write (fo,'(a,i4,i4)') ' Potential (v) matrix:, col, &
             write (1025,'(a,i4,i4)') ' Potential (v) matrix:, col, &
                  &row ',i,j
             write (1025,'(6e14.7)') v(:,i,j)
!             write (fo,'(6f14.7)') v(:,i,j)
          end do
       end do
    end if
  end subroutine potl
end module sector_potl
