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 diag

contains

  subroutine diag (nc, ncp, spins, sqno, imesh)
   use hfile_data, only: lchl, kschl, tchl
   use rw_amplitudes, only: start_sa_file, start_sa_file_id, close_sa_file, ampltd
   use potl_cofs, only: cfs, cfsj, del_cf
   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: inc_lrp_prop, packed_cf, bug9, jpi_coupling
   use error_out, only: error_check
   use mpi_params, only: io_processor, taskid, numtasks
   use magma
   !  use mpi

   integer, intent(in)      :: nc       ! # channels
   integer, intent(in)      :: ncp
   type(qno), intent(in)    :: sqno     ! scattering q#s
   integer, intent(in)      :: spins(2) ! target spins
   integer                  :: imesh    ! COARSE,FINE mesh - imesh=1,2
   real(wp), allocatable    :: r(:)
   real(wp)                 :: ral, rar
   integer                  :: sect, st, status, nh, pp, qq
   integer, save            :: nsect_fine
   integer                  :: nsect_coarse, mesh_taskid,coarse_mesh_first_taskid
   integer ( kind = 4 ) ierr

   call sectors (nc) ! define sector boundaries
!   if (io_processor) then
!      call start_sa_file (sqno, spins(1),sect)
!   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. packed_cf) 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
    end if

   allocate (r(nx), stat=status)
   call error_check (status, 'pdiag: allocation r')
   
   nh = nc * nl
   if(io_processor)  write (fo,'(a,i4)') 'Total Number of Sectors = ', nsect
 
   ! Allow tasks to cycle through both meshes ensuring better load-balancing
   ! Ideally this would be achieved with a shared counter

   if(imesh==1) then ! fine mesh - original taskid numbering for sectors
      mesh_taskid = taskid
      nsect_fine = nsect
   else ! coarse mesh
      nsect_coarse = nsect
      mesh_taskid = taskid - nsect_fine
      coarse_mesh_first_taskid = MOD(nsect_fine,numtasks)
      mesh_taskid = MOD((taskid - coarse_mesh_first_taskid + numtasks),numtasks)  
   end if

   call magmaf_init()
   do sect = 1, nsect
       if (numtasks.gt.1) then
          if (MOD(sect,numtasks)==MOD(mesh_taskid+1,numtasks)) then
             if(imesh==1) then
                write(fo,'(a,i0,a,i0,a,i0)')'FINE Region MPI Task ', taskid, &
                      ' calculating sector ', sect, ' , Sector Hamiltonian matrix dimension = ',nh
             else
                write(fo,'(a,i0,a,i0,a,i0)')'COARSE Region MPI Task ', taskid, &
                      ' calculating sector ', sect, ' , Sector Hamiltonian matrix dimension = ',nh
             end if
             call flush(fo)
          else
             cycle
          end if
       end if
      call start_sa_file_id (sqno, spins(1), sect)
      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 sh_diag (imesh,sect)      ! diagonalize sector H
      call ampltd (ral, rar, nc, nl, nh)
      call kill_az
      call close_sa_file()
   end do

   call magmaf_finalize()
   call h_reset
   call reset_sec
   call reset_potl
   if (inc_lrp_prop.and.packed_cf) call del_cf
   deallocate (r, stat=status)
   call error_check (status, 'pdiag: deallocation')
   
 end subroutine diag

end module pdg_ctl
