module pdg_ctl
! control diagonalization of sector Hamiltonians and RWA calculation
! Time-stamp: "2003-10-14 12:49:30 cjn"
! modified for jpi_coupling april to dec 2006 vmb
  use precisn, only: wp
  use io_units, only: fo
  implicit none

  private
  public pdiag

contains

  subroutine pdiag (nc, ncp, spins, sqno)
   use hfile_data, only: lchl, kschl, tchl
   use rw_amplitudes, only: start_sa_file, ampltd
   use potl_cofs, only: cfs, cfsj, del_cf
   use rmx1_in, only: inc_lrp_prop
   use slp, only: qno
   use def_sectors, only: nsect, nl, nx, w_sect, asect, reset_sec, &
        sectors
   use sector_potl, only: potl, potl_pams, reset_potl
   use dist_mat, only: a_fill, sh_diag, kill_az
   use blacs, only: ctxt, p, q, myrow, mycol, io_processor, p0, &
        q0, p_error
   use sec_ham, only: h_reset, gauleg, legndr, def_rb_vals, xi, def_ncq
   use rmx1_in, only: bug9, jpi_coupling
   integer, intent(in)      :: nc       ! # channels
   integer, intent(in)      :: ncp
   type(qno), intent(in)    :: sqno     ! scattering q#s
   integer, intent(in)      :: spins(2) ! target spins
   real(wp), allocatable    :: r(:)
   real(wp)                 :: ral, rar
   integer                  :: sect, st, status, nh, pp, qq

   call sectors (nc) ! define sector boundaries
   if (io_processor) then
      call start_sa_file (sqno, spins(1))
   end if

   call gauleg  ! quadrature formula
   call legndr  ! Legendre functions
   call potl_pams (nc, nx, ncp)
   call def_ncq (ncp) ! pass ethr offset in h_el

! generate asymptotic potential coefficients
   if (inc_lrp_prop) then
   if (.not. jpi_coupling) then
      call cfs (nc, tchl(ncp:), lchl(ncp:), sqno, spins)
   else
      call cfsj (nc, tchl(ncp:), lchl(ncp:), kschl(ncp:), sqno, spins)
   end if
   end if
   allocate (r(nx), stat=status)
   call p_error (status, 'pdiag: allocation r')
   
   nh = nc * nl
   if (io_processor .and. bug9 > 0) then
      write (fo,'(a,i4)') 'Total Number of Sectors = ', nsect
      write (fo,'(a,i4)') 'Partition Hamiltonian dimension = ', nh
   end if

   do sect = 1, nsect
      ral = asect(sect)
      rar = asect(sect+1)
      r = 0.5_wp * ((rar - ral) * xi + rar + ral)
      call def_rb_vals (ral, rar)
      call potl (nc, nx, r)
      call A_fill (nh) ! form distributed sector H
      call BLACS_BARRIER (ctxt, 'A')
      call sh_diag      ! diagonalize sector H
      call ampltd (ral, rar, nc, nl, nh)
      call kill_az
      call BLACS_BARRIER (ctxt, 'A')
   end do
   call h_reset
   call reset_sec
   call reset_potl
   if (inc_lrp_prop) call del_cf
   deallocate (r, stat=status)
   call p_error (status, 'pdiag: deallocation')
 end subroutine pdiag

end module pdg_ctl
