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, ne, bug1,       &
       start_mesh, jpi_coupling
  use hfile_data, only: readh1, degtar, readh2, get_spins, &
       reorder_chls, cparm, nc, nc1, nc2, &
       readhj1, readhj2, get_split
  use energy_grid, only: energy_max
  use potl_cofs, only: rdtmom, rdjtmom
  use pdg_ctl, only: diag
  use error_out, only: error_check
  use mpi_params, only: taskid, numtasks, io_processor
  use mpi

  !  use blacs, only: nprocs, iam, ictxt, io_processor, p, q, ctxt, myrow,&
  !       mycol, p0, q0, set_cur_ctxt, p_error
  implicit none
  !  include 'mpif.h'
  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(:,:)
  integer ( kind = 4 )        :: errcode, ierr

  call MPI_INIT( ierr )
  call MPI_COMM_RANK( MPI_COMM_WORLD, taskid, ierr )
  call MPI_COMM_SIZE( MPI_COMM_WORLD, numtasks, ierr )
  if(taskid==0) io_processor=.true.

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

  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 MPI Tasks = ',numtasks 
  end if

  call rd_nmlst    ! read input namelist

  ! 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)
  if(io_processor) then
     write (fo,'(/a4,i3,a8,i3,a8,i2)')  'S = ', sqno%nspn,&
          '    L = ', sqno%lrgl, '    P = ', sqno%npty
  end if
  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).and.io_processor) 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)

  ! MPI/LAPACK/GPU version only supports nmesh<=2
  if (nmesh.gt.2) then
    if(io_processor) write(fo,'(a)') '**** Stopping run - MPI/LAPACK/GPU version does not support runs with nmesh > 2 ****'
    call MPI_ABORT(MPI_COMM_WORLD,errcode,ierr)
  end if
  mesh_loop: do imesh = start_mesh, nmesh
     call energy_max (imesh, nmesh)   ! find ebig
     !        call distribute_edata ! broadcast ebig

     ! initialize full or 1st split partition sector data
     ncp = 1
     call diag (nc1, ncp, (/st1,st2/), sqno, imesh)

     if (nc1 /= nc) then ! initialize 2nd split partition data
        nc2 = nc - nc1
        ncp = nc1 + 1
        call diag (nc2, ncp, (/st2, -999/), sqno, imesh)
     end if
  end do mesh_loop

  call cpu_time (t1)
  write (fo,'(a,i0,a,f16.4,a)') 'Task ', taskid, ' ; End of RMX1 - Local CPU time     = ', t1 - t0, ' secs'
  call system_clock (count=c1, count_rate=cr)
  write (fo,'(a,i0,a,f16.4,a)') 'Task ', taskid, ' ; End of RMX1 - Local Elapsed time = ', REAL(c1-c0,wp) / REAL(cr,wp), ' secs'

  call MPI_BARRIER(MPI_COMM_WORLD,ierr)
  if (io_processor) then
     call cpu_time (t1)
     write (fo,'(a,f16.4,a)') 'End of RMX1 - Global CPU time     = ', t1 - t0, ' secs'
     call system_clock (count=c1, count_rate=cr)
     write (fo,'(a,f16.4,a)') 'End of RMX1 - Global Elapsed time = ', REAL(c1-c0,wp) / &
          REAL(cr,wp), ' secs'
  end if
  call MPI_FINALIZE(ierr)
end program rmx1
