program rmx1
!  Electron-atom R-matrix external region calculation - STAGE 1
! Time-stamp: "2003-09-23 17:06:43 cjn"
! modified for jpi_coupling april to dec 2006 vmb
  use precisn, only: wp
  use io_units, only: fo
  use slp, only: qno, def_qno
  use rmx1_in, only: targ_low_sce, targ_high_sce, rd_nmlst, split_prop,&
       degeny, lrgl1, nspn1, npty1, distribute_indata, ne, bug1,       &
       start_mesh, jpi_coupling
  use hfile_data, only: readh1, degtar, readh2, get_spins, &
       reorder_chls, cparm, distribute_hdata, nc, nc1, nc2, &
       readhj1, readhj2, distribute_hjdata, get_split
  use energy_grid, only: energy_max, distribute_edata
  use potl_cofs, only: rdtmom, rdjtmom
  use pdg_ctl, only: pdiag
  use blacs, only: nprocs, iam, ictxt, io_processor, p, q, ctxt, myrow,&
       mycol, p0, q0, set_cur_ctxt, p_error
  implicit none
  real(wp)                    :: e0
  real(wp), allocatable, save :: ethr(:) ! threshold energies
  integer, pointer            :: en_order(:) ! energy-order chl index
  real(wp)                    :: rafin, rmatr ! prop. limits
  integer                     :: smx = 1
  type(qno)                   :: sqno
  integer                     :: i, j, ion, ncp, ij, status
  integer                     :: st1, st2
  real(wp)                    :: t0, t1
  integer                     :: nmesh, imesh
  integer                     :: imycol, imyrow, ip, iq
  integer                     :: BLACS_PNUM
  integer                     :: c0, c1, cr
  integer, allocatable        :: map(:,:)

  call cpu_time (t0)
  call system_clock (count=c0)

  call BLACS_PINFO (iam, nprocs) ! find process #, total # processors
  call BLACS_GET (-1, 0, ictxt)  ! find default context, ctxt
! define blacs grid to include all processes to broadcast information
  call BLACS_GRIDINIT (ictxt, 'Row-major', 1, nprocs)
  call BLACS_GRIDINFO (ictxt, ip, iq, imyrow, imycol)
  io_processor = (imycol == 0)
  call set_cur_ctxt (1)

  if (io_processor) then ! This is the i/o processor
     write (fo,'(//,15x,a)') '============='
     write (fo,'(15x,a)') 'Program RMX1'
     write (fo,'(15x,a,//)') '============='
     write (fo,'(a)') 'R-Matrix External Region: Sector Diagonalization'
     write (fo,'(a,i6)') 'Number of processors = ', nprocs

     call rd_nmlst    ! read input namelist
     write (fo,'(a,2i6)') 'Blacs process dimensions, p, q = ', p, q

 ! targ_high_sce > 0 signals scattering energies between targets targ_low_sce and
! targ_high_sce. i.e. energies in the resonance region => nmesh = 1
! ne(1) > 0 signals only a single step
     if (targ_high_sce > 0 .or. ne(1) > 0) then
        nmesh = 1
     else
        nmesh = 2
     end if

     if (.not. jpi_coupling) then
        call readh1    !  read first record on H file
     else
        call readhj1 
     end if 

! set near degenerate target states as strictly degenerate: iterate
     degeny = 0.1_wp * degeny
     call degtar (degeny)
     degeny = 5.0_wp * degeny
     call degtar (degeny)
     degeny = 2.0_wp * degeny
     call degtar (degeny)
     call degtar (degeny)

! define required SLp or Jp case:
     sqno = def_qno(lrgl1, nspn1, npty1)
     write (fo,'(/a4,i3,a8,i3,a8,i2)')  'S = ', sqno%nspn,&
          '    L = ', sqno%lrgl, '    P = ', sqno%npty
     if (.not. jpi_coupling) then
        call readh2 (lrgl1, nspn1, npty1)
     else
        call readhj2 (lrgl1, nspn1, npty1)
     end if
     call get_split (st1, st2)

     if (split_prop .and. st2 /= -999) then ! there is channel splitting
        call reorder_chls
        nc2 = nc - nc1
        if (bug1 > 0) then
           write (fo,'(/a/)') 'Channel Splitting Used'
           write (fo,'(a,i6)') 'Number of spin 1 Channels = ', nc1
           write (fo,'(a,i6/)') 'Number of Spin 2 channels = ', nc2
        end if
     else   ! no channel_splitting
        nc1 = nc
        split_prop = .false.
     end if
     call cparm   ! form ethr using reordered channel data
  end if

  call distribute_indata
  if (.not. jpi_coupling ) then
     call rdtmom ! read target transition moments and distribute
     call distribute_hdata (st1, st2, lrgl1, nspn1, npty1, start_mesh, &
       nmesh)
  else
     call rdjtmom 
     call distribute_hjdata (st1, st2, lrgl1, nspn1, npty1, start_mesh, &
       nmesh)
  end if
  sqno = def_qno(lrgl1, nspn1, npty1)
  if (io_processor) io_processor = .false.

  allocate (map(p,q), stat=status)
  call p_error (status, 'rmx1: allocate map')
  ij = 0
  do j = 1, q
     do i = 1, p
        map(i,j) = ij
        ij = ij + 1
     end do
  end do
!  call BLACS_GRIDEXIT (ctxt)     ! kill current mesh
  call BLACS_GET (-1, 0, ctxt)  ! find default context, ctxt
  call BLACS_GRIDMAP (ctxt, map, p, p, q)
! call BLACS_GRIDINIT (ctxt, 'Row-major', p, q)
  call BLACS_GRIDINFO (ctxt, ip, iq, myrow, mycol)
  if (myrow >= 0 .and. myrow < p .and. mycol >= 0 .and. mycol < q) then
     p0 = 0
     q0 = 0
     io_processor = (myrow == p0 .and. mycol == q0)
     call set_cur_ctxt (2)
     mesh_loop: do imesh = start_mesh, nmesh
        if (io_processor) call energy_max (imesh, nmesh)   ! find ebig
        call distribute_edata ! broadcast ebig

 ! initialize full or 1st split partition sector data
        ncp = 1
        call pdiag (nc1, ncp, (/st1,st2/), sqno)
        
        if (nc1 /= nc) then ! initialize 2nd split partition data
           nc2 = nc - nc1
           ncp = nc1 + 1
           call pdiag (nc2, ncp, (/st2, -999/), sqno)
        end if
     end do mesh_loop

     if (io_processor) then
        write (fo,'(/,a,/)') 'end of RMX1'
        call cpu_time (t1)
        write (fo,'(a,f16.4,a)') 'CPU time     = ', t1 - t0, ' secs'
        call system_clock (count=c1, count_rate=cr)
        write (fo,'(a,f16.4,a)') 'Elapsed time = ', REAL(c1-c0,wp) / &
             REAL(cr,wp), ' secs'
     end if
     call BLACS_GRIDEXIT (ctxt)
  end if
  call BLACS_BARRIER (ictxt, 'All')
  call BLACS_EXIT()
  stop
end program rmx1
