module hfile_data
! read and store target and channel data from HAM H-file
! Time-stamp: "2003-10-03 17:28:25 cjn"
! reading of jpi_coupling data included april to dec 2006 vmb
  use precisn, only: wp
  use io_units, only: fo
  implicit none

  integer, save                :: iunit = 77

  integer, save                :: nelc, nz, lrang2, lamax, ntarg
  real(wp), save               :: rmatr, bbloch
  integer, allocatable, save   :: jtarg(:), ltarg(:), starg(:), ptarg(:)
  real(wp), allocatable, save  :: etarg(:)
  integer, save                :: nc, nc1, nc2, ns
  integer, allocatable, save   :: lchl(:), tchl(:), kschl(:)
  real(wp), allocatable, save  :: ethr(:)

  private
  public readh1, readh2, cparm, degtar, asc_order
  public readhj1, readhj2
  public get_spins, get_split, split_chls, reorder_chls
  public nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
  public jtarg, ltarg, starg, ptarg, etarg
  public lchl, tchl, kschl, nc, ns, ethr, nc1, nc2

contains

  subroutine cparm
!   cparm : channel parameters eth, lch
    real(wp)                :: e0      ! initial target energy
    integer                  :: i, status

    allocate (ethr(nc), stat=status)
    e0 = MINVAL(etarg(1:ntarg))
    do i = 1, nc
       ethr(i) = etarg(tchl(i)) - e0
    end do
  end subroutine cparm

  subroutine asc_order (en, nc, en_order)
!  define pointer array eorder, giving values of en in ascending order
    integer, intent(in)   :: nc           ! no scattering channels
    real(wp), intent(in)  :: en(nc)       ! array to be ordered
    integer, intent(out)  :: en_order(nc) ! order pointers to en
    logical               :: mask(nc)     ! automatic mask array
    integer :: i, j(1)

    mask(1:nc) = .true.
    do i = 1, nc
       j = MINLOC(en,mask)
       en_order(i) = j(1)
       mask(j(1)) = .false.
    end do
  end subroutine asc_order

  subroutine readh1
! Read part of asymptotic file that is independent of lrgl, spin, parity
    use xdr_files, only: open_xdr, xdr_io
    use rmx1_in, only: filh, xdr, bug2
    use scaling, only: set_charge, scale_radius, scale_etarg
    integer                 :: ihbuf(5)
    real(wp)                :: rhbuf(2), etgr
    integer                 :: i, n, l, status
    real(wp), allocatable   :: cfbut(:)
    real(wp)                :: e0, scale


    write (fo,'(a)') 'Reading H file: ' // TRIM(filh)

    if (xdr) then
       iunit = open_xdr (file=TRIM(filh), action='read')
       call xdr_io (iunit, ihbuf, 5)
       nelc = ihbuf(1)
       nz = ihbuf(2)
       lrang2 = ihbuf(3)
       lamax = ihbuf(4)
       ntarg = ihbuf(5)
       call xdr_io (iunit, rhbuf, 2)
       rmatr = rhbuf(1)
       bbloch = rhbuf(2)

       allocate (etarg(ntarg), ltarg(ntarg), starg(ntarg), &
            ptarg(ntarg), stat=status)

       call xdr_io (iunit, etarg, ntarg)
       call xdr_io (iunit, ltarg, ntarg)
       call xdr_io (iunit, starg, ntarg)
       call xdr_io (iunit, ptarg, ntarg)
       allocate (cfbut(3), stat=status)
       do l = 1, lrang2
          call xdr_io (iunit, cfbut, 3)   ! skip Buttle data for rmx1
       end do
    else          ! normal fortran unformatted output
       read (iunit) nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
       allocate (etarg(ntarg), ltarg(ntarg), starg(ntarg), &
            ptarg(ntarg), cfbut(3*lrang2),  stat=status)
       read (iunit) etarg(1:ntarg)
       read (iunit) ltarg(1:ntarg)
       read (iunit) starg(1:ntarg)
       read (iunit) ptarg(1:ntarg)
       read (iunit) cfbut
    end if

    if (bug2 > 0) then
       write (fo,'(a)') 'Subroutine readh1: sector matrix setup'
       write (fo,'(a)') 'Data on the H file:'
       write (fo,'(5i5,2f12.6)') nelc, nz, lrang2, lamax, ntarg, &
            rmatr, bbloch
       write (fo,'(5f14.6)') etarg
       write (fo,'(5i6)') ltarg
       write (fo,'(5i6)') starg
       write (fo,'(5i6)') ptarg
    end if

! process input data:
    call set_charge (nz, nelc)
    call scale_radius (rmatr)
    call scale_etarg (etarg)     ! etarg now in scaled Ryd
    write (fo,'(/,a,/)') 'Target states'
    write (fo,'(10x,a,5x,a,3x,a,8x,a,/,43x,a)') 'index', 'total l', &
         '(2*s+1)', 'energy', 'scaled ryd'
    e0 = etarg(1)
    do i = 1, ntarg
       etgr = (etarg(i) - e0)           ! convert to Rydbergs
       write (fo,'(3x,3i10,3x,f12.6)') i, ltarg(i), starg(i), etgr
    end do
    deallocate (cfbut, stat=status)
  end subroutine readh1

  subroutine readhj1
! Read part of asymptotic file that is independent of lrgl, spin, parity
! here lrgl corresponds to 2J and spin is zero
    use xdr_files, only: open_xdr, xdr_io
    use rmx1_in, only: filh, xdr, bug2
    use scaling, only: set_charge, scale_radius, scale_etarg
    integer                 :: ihbuf(5)
    real(wp)                :: rhbuf(2), etgr
    integer                 :: i, n, l, status
    real(wp), allocatable   :: cfbut(:)
    real(wp)                :: e0, scale


    write (fo,'(a)') 'Reading H file: ' // TRIM(filh)

    if (xdr) then
       iunit = open_xdr (file=TRIM(filh), action='read')
       call xdr_io (iunit, ihbuf, 5)
       nelc = ihbuf(1)
       nz = ihbuf(2)
       lrang2 = ihbuf(3)
       lamax = ihbuf(4)
       ntarg = ihbuf(5)
       call xdr_io (iunit, rhbuf, 2)
       rmatr = rhbuf(1)
       bbloch = rhbuf(2)

       allocate (etarg(ntarg), jtarg(ntarg), ltarg(ntarg), &
            starg(ntarg), ptarg(ntarg), stat=status)
       call xdr_io (iunit, etarg, ntarg)
       call xdr_io (iunit, jtarg, ntarg)
       call xdr_io (iunit, ltarg, ntarg)
       call xdr_io (iunit, starg, ntarg)
       call xdr_io (iunit, ptarg, ntarg)
       allocate (cfbut(3), stat=status)
       do l = 1, lrang2
          call xdr_io (iunit, cfbut, 3)   ! skip Buttle data for rmx1
       end do
    else          ! normal fortran unformatted output
       read (iunit) nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
       allocate (etarg(ntarg), jtarg(ntarg), ltarg(ntarg), starg(ntarg), &
            ptarg(ntarg), cfbut(3*lrang2),  stat=status)
       read (iunit) etarg(1:ntarg)
       read (iunit) jtarg(1:ntarg)
       read (iunit) ltarg(1:ntarg)
       read (iunit) starg(1:ntarg)
       read (iunit) ptarg(1:ntarg)
       read (iunit) cfbut
    end if

    if (bug2 > 0) then
       write (fo,'(a)') 'Subroutine readhj1: sector matrix setup'
       write (fo,'(a)') 'Data on the H file:'
       write (fo,'(5i5,2f12.6)') nelc, nz, lrang2, lamax, ntarg, &
            rmatr, bbloch
       write (fo,'(5f14.6)') etarg
       write (fo,'(5i6)') jtarg
       write (fo,'(5i6)') ltarg
       write (fo,'(5i6)') starg
       write (fo,'(5i6)') ptarg
    end if

! process input data:
    call set_charge (nz, nelc)
    call scale_radius (rmatr)
    call scale_etarg (etarg)     ! etarg now in scaled Ryd
    write (fo,'(/,a,/)') 'Target states'
    write (fo,'(10x,a,5x,a,3x,a,3x,a,8x,a,/,43x,a)') 'index', &
        '  2j   ', 'total l', '(2*s+1)', 'energy', 'scaled ryd'
    e0 = etarg(1)
    do i = 1, ntarg
       etgr = (etarg(i) - e0)           ! convert to Rydbergs
       write (fo,'(3x,4i10,3x,f12.6)') i, jtarg(i), ltarg(i), &
                                       starg(i), etgr
    end do
    deallocate (cfbut, stat=status)
  end subroutine readhj1

  subroutine get_spins (st1, st2)
! find target spins which contribute to scattering SLp symmetry
    integer, intent(out)   :: st1, st2 ! target spins for SLp case
    integer                :: i

    st2 = -999
    st1 = starg(tchl(1))  ! 1st target spin
    do i = 2, nc
       if (starg(tchl(i)) == st1) cycle
       st2 = starg(tchl(i))
       exit
    end do
  end subroutine get_spins

  subroutine split_chls_old (schl1, schl2)
! split the channels into two possible spins
    integer, intent(out)   :: schl1(nc), schl2(nc) ! channel seqs
    integer                :: ispch(nc) ! new channel order
    integer                :: st1, st2
    integer                :: i, nts

    st1 = starg (tchl(1))    ! first target spin
    nts = 1
    nc1 = 1
    nc2 = 0
    schl1(1) = 1
    channels: do i = 2, nc
       if (starg(tchl(i)) == st1) then ! first target channel
          nc1 = nc1 + 1
          schl1(nc1) = i
       else
          if (nts == 1) then   ! found a second target spin
             nts = nts + 1
             st2 = starg(tchl(i))
             nc2 = 1
             schl2(nc2) = i
          else if (starg(tchl(i)) == st2) then ! second target channel
             nc2 = nc2 + 1
             schl2(nc2) = i
          else
             write (fo,'(a)') 'split_chls: Channel data inconsistent'
             write (fo,'(a,2i6)') 'st1, st2, i, s = ', st1, st2, i, &
                  starg(tchl(i))
             stop
          end if
       end if
    end do channels
    if (nc /= nc1 + nc2) then
       write (fo,'(a,3i4)') 'split_chls: error, nc1, nc2, nc = ', &
            nc1, nc2, nc
       stop
    end if
  end subroutine split_chls_old

  subroutine get_split (st1, st2)
! find splitting parameters which contribute to scattering symmetry
! in SLP symmetries splitting parameters correspond to target spin
! in JP symmetries splitting parameters correspond to channel q no k
    integer, intent(out)   :: st1, st2 ! splitting parameters
    integer                :: i

    st2 = -999
    st1 = kschl(1)  ! 1st splitting parameter
    do i = 2, nc
       if (kschl(i) == st1) cycle
       st2 = kschl(i)
       exit
    end do
  end subroutine get_split

  subroutine split_chls (schl1, schl2)
! this is changed but should be applicable to LS and J
! split the channels into two possible spins
    integer, intent(out)   :: schl1(nc), schl2(nc) ! channel seqs
    integer                :: ispch(nc) ! new channel order
    integer                :: st1, st2
    integer                :: i, nts

    st1 = kschl(1)    ! first ks value, either target spin or k
    nts = 1
    nc1 = 1
    nc2 = 0
    schl1(1) = 1
    channels: do i = 2, nc
       if (kschl(i) == st1) then ! first channel split
          nc1 = nc1 + 1
          schl1(nc1) = i
       else
          if (nts == 1) then   ! found a second ks value
             nts = nts + 1
             st2 = kschl(i)
             nc2 = 1
             schl2(nc2) = i
          else if (kschl(i) == st2) then ! second channel split
             nc2 = nc2 + 1
             schl2(nc2) = i
          else
             write (fo,'(a)') 'split_chls: Channel data inconsistent'
             write (fo,'(a,2i6)') 'st1, st2, i, ks = ', st1, st2, i, &
                    kschl(i)
             stop
          end if
       end if
    end do channels
    if (nc /= nc1 + nc2) then
       write (fo,'(a,3i4)') 'split_chls: error, nc1, nc2, nc = ', &
            nc1, nc2, nc
       stop
    end if
    write (fo,*) ' nc1,nc2 ', nc1,nc2
    write (fo,*) 'schl1 ', (schl1(i),i=1,nc1)
    write (fo,*) 'schl2 ', (schl2(i),i=1,nc2)
  end subroutine split_chls

  subroutine reorder_chls
! order channel arrays lchl, tchl, kschl according to split parameters
    integer, allocatable      :: tmp_lchl(:), tmp_tchl(:), tmp_kschl(:)
    integer, allocatable      :: schl1(:), schl2(:)
    integer                   :: i, j, ip, jp, nc2, nct, status

    allocate (schl1(nc), schl2(nc), stat=status)
    call split_chls (schl1, schl2)
    if (nc1 == nc) then ! only one target spin for this SLp case
       deallocate (schl1, schl2, stat=status)
       return
    end if
    nc2 = nc - nc1
    allocate (tmp_lchl(nc), tmp_tchl(nc), tmp_kschl(nc), stat=status)
    tmp_lchl = lchl
    tmp_tchl = tchl
    tmp_kschl = kschl
    do i = 1, nc1
       ip = schl1(i)
       lchl(i) = tmp_lchl(ip)
       tchl(i) = tmp_tchl(ip)
       kschl(i) = tmp_kschl(ip)
    end do
    do i = 1, nc2
       jp = nc1 + i
       ip = schl2(i)
       lchl(jp) = tmp_lchl(ip)
       tchl(jp) = tmp_tchl(ip)
       kschl(jp) = tmp_kschl(ip)
    end do
    deallocate (tmp_lchl, tmp_tchl, tmp_kschl, stat=status)
  end subroutine reorder_chls

  subroutine degtar (degeny)
! make target states truly degenerate according to degeny (scaled Rydbergs)
! etarg is in scaled Rydbergs and in energy order
    real(wp), intent(in) :: degeny  ! degeneracy criterion
    real(wp)             :: etdeg(ntarg) ! tmp modified target energies
    integer               :: i, iloop, j, nd
    real(wp)             :: eav, esum, etargi

    i = 1
    do iloop = 1, ntarg
       etargi = etarg(i)
       esum = etargi
       nd = 1
       do j = i+1, ntarg
          if ((etarg(j)- etargi) > degeny) exit
          esum = esum + etarg(j)
          nd = nd + 1
       end do
       eav = esum / nd ! average of near degenerate states
       do j = i, i+nd-1
          etdeg(j) = eav
       end do
       i = i + nd
       if (i > ntarg) exit
    end do
    etarg = etdeg
  end subroutine degtar

  subroutine readh2 (lrgl1, nspn1, npty1)
! read H-file data for a particular scattering symmetry
    use xdr_files, only: xdr_io
    use rmx1_in, only: xdr, bug3
    integer, intent(in)      :: lrgl1   ! scattering orb a.m.
    integer, intent(in)      :: nspn1   ! scattering spin multiplicity
    integer, intent(in)      :: npty1   ! scattering parity
    integer, allocatable     :: nltarg(:)
    real(wp), allocatable    :: wmat(:), eig(:)
    integer                  :: loc, lrgl2, nspn2, npty2
    integer                  :: ibuf(6)
    integer                  :: i, j, k, it, nt, n, more2, nl, status
    integer                  :: ncs
    logical                  :: match

    more2 = 1
    match = .false.
    nt = ntarg
    symmetries: do while (more2 == 1) ! search for required SLp
       if (xdr) then
          call xdr_io (iunit, ibuf, 6)
          lrgl2 = ibuf(1);   nspn2 = ibuf(2);  npty2 = ibuf(3)
          nc = ibuf(4);  ns = ibuf(5);   more2 = ibuf(6)
          ncs = MAX(nc, ns)
          allocate (lchl(nc), nltarg(nt), tchl(nc), kschl(nc), &
                    eig(ns), wmat(ncs), stat=status)
          call xdr_io (iunit, nltarg, nt)
          call xdr_io (iunit, lchl, nc)
          call xdr_io (iunit, eig, ns)  ! skip over eig, wmat
          do j = 1, ns
             call xdr_io (iunit, loc)   ! state number
             call xdr_io (iunit, wmat(1:nc), nc)
          end do
       else
          read (iunit) lrgl2, nspn2, npty2, nc, ns, more2
          ncs = MAX(nc, ns)
          allocate (lchl(nc), nltarg(nt), tchl(nc), kschl(nc), &
                    eig(ns), wmat(ncs), stat=status)
          read (iunit) nltarg
          read (iunit) lchl
          read (iunit) eig
          do j = 1, ns
             read (iunit) loc
             read (iunit) wmat(1:nc)
          end do
       end if

! form tchl array:
       i = 0
       do it  = 1, nt
          nl = nltarg(it)
          do n = 1, nl
             i = i + 1
             tchl(i) = it
          end do
       end do

! form kschl array:  ! this array will be used for channel splitting
       do i = 1, nc
          kschl(i) = starg(tchl(i))
       end do

       if (lrgl1 == lrgl2 .and. nspn1 == nspn2 .and. npty1 == npty2) &
            then
          match = .true.
          deallocate (eig, wmat, nltarg, stat=status)
          exit
       end if
       deallocate (lchl, tchl, kschl, eig, wmat, nltarg, stat=status)
       if (more2 == 2) exit
       exit
    end do symmetries

    if (.NOT.match) then
       write (fo,'(a,3i5)') 'readh2: unmatched symmetry, SLp = ', &
            nspn1, lrgl1, npty1
       stop
    end if
    if (bug3 /= 0) then
       write (fo,'(a,i6,a,i6)') 'nc = ', nc, '    ns = ', ns
       write (fo,'(/,a)') 'Channel data:'
       write (fo,'(a)') '         target sequence   orbital a.m.'
       do i = 1, nc
          write (fo,'(i4,15x,i4,5x,i4)') i, tchl(i), lchl(i)
       end do
    end if
  end subroutine readh2

  subroutine readhj2 (lrgl1, nspn1, npty1)
! read H-file data for a particular scattering symmetry
! here lrgl corresponds to 2J and spin multiplicity is zero
    use xdr_files, only: xdr_io
    use rmx1_in, only: xdr, bug3
    integer, intent(in)      :: lrgl1   ! scattering orb a.m.=2J or -2K
    integer, intent(inout)   :: nspn1   ! 0 or weighting factor info for K-file
    integer, intent(in)      :: npty1   ! scattering parity
    integer, allocatable     :: nltarg(:)
    real(wp), allocatable    :: wmat(:), eig(:)
    integer                  :: loc, lrgl2, nspn2, npty2
    integer                  :: ibuf(6)
    integer                  :: i, j, k, it, nt, n, more2, nl, status
    integer                  :: ncs
    logical                  :: match

    more2 = 1
    match = .false.
    nt = ntarg
    symmetries: do while (more2 == 1) ! search for required SLp
       if (xdr) then
          call xdr_io (iunit, ibuf, 6)
          lrgl2 = ibuf(1);   nspn2 = ibuf(2);  npty2 = ibuf(3)
          nc = ibuf(4);  ns = ibuf(5);   more2 = ibuf(6)
          ncs = MAX(nc, ns)
          allocate (lchl(nc), kschl(nc), nltarg(nt), tchl(nc), eig(ns),& 
                   wmat(ncs), stat=status)

          call xdr_io (iunit, nltarg, nt)
          call xdr_io (iunit, lchl, nc)
          call xdr_io (iunit, kschl, nc)
          call xdr_io (iunit, eig, ns)  ! skip over eig, wmat
          do j = 1, ns
             call xdr_io (iunit, loc)   ! state number
             call xdr_io (iunit, wmat(1:nc), nc)
          end do
       else
          read (iunit) lrgl2, nspn2, npty2, nc, ns, more2
          ncs = MAX(nc, ns)
          allocate (lchl(nc), kschl(nc), nltarg(nt), tchl(nc), eig(ns),& 
                   wmat(ncs), stat=status)

          read (iunit) nltarg
          read (iunit) lchl
          read (iunit) kschl
          read (iunit) eig
          do j = 1, ns
             read (iunit) loc
             read (iunit) wmat(1:nc)
          end do
       end if

! form tchl array:
       i = 0
       do it  = 1, nt
          nl = nltarg(it)
          do n = 1, nl
             i = i + 1
             tchl(i) = it
          end do
       end do

       if (lrgl1 == lrgl2 .and. npty1 == npty2) then !nspn1 not defined by namelist
          match = .true.
          deallocate (eig, wmat, nltarg, stat=status)
          exit
       end if
       nspn1 = nspn2 !use the value on the h-file (weighting factor for K-files)
       deallocate (lchl, kschl, tchl, eig, wmat, nltarg, stat=status)
       if (more2 == 2) exit
       exit
    end do symmetries

    if (.NOT.match) then
       write (fo,'(a,3i5)') 'readhj2: unmatched symmetry, SLp = ', &
            nspn1, lrgl1, npty1
       stop
    end if
    if (bug3 /= 0) then
       write (fo,'(a,i6,a,i6)') 'nc = ', nc, '    ns = ', ns
       write (fo,'(/,a)') 'Channel data:'
       write (fo,'(a)') '  target index       l          2k'
       do i = 1, nc
          write (fo,'(i4,6x,i4,5x,i4,5x,i4)') i, tchl(i), lchl(i),&
                 kschl(i)
       end do
    end if
  end subroutine readhj2

end module hfile_data
