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
  use mpi_params, only: io_processor
  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(:)
  integer, allocatable, save   :: ncf_nonzero(:), ncf2_nonzero(:)
  integer, save                :: ncf_m
  integer, allocatable, save   :: ic_label(:,:), ir_label(:,:)
! these three arrays give the two channel numbers for the condensed cf_nz for each lambda,
! and the number of non-zero elements for each lambda
! nb the matrix is assumed to be symmetric
  real(wp), allocatable, save  :: cf_farm(:,:,:), cf_nz(:,:)
  integer, allocatable, save   :: schl1(:), schl2(:)  ! needed for picking out cf_farm values
  integer, allocatable, save   :: stl1(:), stl2(:)  ! needed for picking out packed cf_nz values
  integer, save                :: spin_v   ! lets cfs know which schl{1,2} to use

  private
  public readh1, readh2, cparm, degtar, asc_order
  public readhj1, readhj2
  public get_spins, get_split, split_chls, reorder_chls
!  public distribute_hdata, distribute_hjdata
  public nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
  public jtarg, ltarg, starg, ptarg, etarg
  public lchl, tchl, kschl, nc, ns, ethr, nc1, nc2
  public schl1, schl2, spin_v, cf_farm
  public ncf_nonzero, ncf2_nonzero, ic_label, ir_label, cf_nz

contains

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

    allocate (ethr(nc), stat=status)
    if (status /= 0) then 
       write(*,*) ' *** Allocation error: Subroutine cparm *** '
       STOP
    end if 
    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_H_in, bug2, farm_format, molecule_format
    use scaling, only: set_charge, scale_radius, scale_etarg
    integer                 :: ihbuf(5)
    real(wp)                :: rhbuf(2), etgr
    integer                 :: i, l, status
    real(wp), allocatable   :: cfbut(:)
    real(wp)                :: e0
    real(wp)                :: r_mol, rmass ! UKRmol dummy parms, needed in prop K/T header output
    real(wp)                :: ezero !UKRmol  partitioned R-matrix parameter
    integer                 :: ibut, iex, nchan_choice ! UKRmol Buttle/partitioned input
    integer                 :: but_dim, iex_d  ! dummy for Buttle array, iex arrays
    real, allocatable       :: sfac(:),ecex(:), rcex(:,:) ! UKRmol partitioned R-matrix arrays

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

    if (xdr_H_in) 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)
       if (status /= 0) then 
         write(*,*) ' *** Allocation error etarg: Subroutine readh1 *** '
         STOP
       end if 
       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)
       if (status /= 0) then 
         write(*,*) ' *** Allocation error cfbut: Subroutine readh1 *** '
         STOP
       end if 
       do l = 1, lrang2
          call xdr_io (iunit, cfbut, 3)   ! skip Buttle data for rmx1
       end do
    else          ! normal fortran unformatted output
       open (unit=iunit, file=TRIM(filh), status='old', action='read', iostat=status, form='unformatted')
       if(status.ne.0) write(*,*) ' Opening ', TRIM(filh), ' Error iostat = ', status
       read (iunit) nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
       ibut = 1
       if (molecule_format) read (iunit) r_mol, rmass, ezero, ibut, iex, nchan_choice
! note that for molecule_format, 'lrang2' is actually the number of channels
       but_dim = lrang2 * 3
       allocate (etarg(ntarg), ltarg(ntarg), starg(ntarg), &
            ptarg(ntarg), cfbut(but_dim),  stat=status)
       if (status /= 0) then 
         write(*,*) ' *** Allocation error etarg: Subroutine readh1 *** '
         STOP
       end if 
       read (iunit) etarg(1:ntarg)
       read (iunit) ltarg(1:ntarg)
       read (iunit) starg(1:ntarg)
       if (.not. farm_format) read (iunit) ptarg(1:ntarg)
! note slight cfbut inconsistency here for farm_format = .false. .
!This corresponds to ham95/fine 95 output as of January 2011.
       if (ibut == 1) then
          read (iunit) cfbut
          deallocate (cfbut, stat=status)
       else if (abs(ibut) == 2) then
! partitioned R-matrix, 'lrang2' = number of channels nchan
          deallocate (cfbut, stat=status)
          iex_d = MAX(iex,1)
          allocate (sfac(lrang2), ecex(iex_d), rcex(lrang2,iex_d), stat=status)
          if (status /= 0) then 
            write(*,*) ' *** Allocation error sfac: Subroutine readh1 *** '
            STOP
          end if 
          read (iunit) sfac(1:lrang2)
          if (iex /= 0)then
             read (iunit) ecex(1:iex)
             read (iunit) rcex(1:lrang2,1:iex) 
          end if
          deallocate(rcex, ecex, sfac, stat=status)
       end if
       if (status /= 0) then 
         write(*,*) ' *** Allocation problem: Subroutine readh1 *** '
         STOP
       end if 
    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
       if (farm_format) then
          write (fo,'(a)') 'Farm binary format: etarg, ltarg, starg'
          write (fo,'(5f14.6)') etarg
          write (fo,'(5i6)') ltarg
          write (fo,'(5i6)') starg
       else
          write (fo,'(a)') 'Pfarm preferred xdr/binary format: etarg, ltarg, starg, ptarg'
          write (fo,'(5f14.6)') etarg
          write (fo,'(5i6)') ltarg
          write (fo,'(5i6)') starg
          write (fo,'(5i6)') ptarg
       end if
    end if

! process input data:
    call set_charge (nz, nelc)
    call scale_radius (rmatr)
    call scale_etarg (etarg)     ! etarg now in scaled Ryd
    if(io_processor) then
      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'
    end if  
    e0 = etarg(1)
    do i = 1, ntarg
       etgr = (etarg(i) - e0)           ! convert to Rydbergs
       if(io_processor) write (fo,'(3x,3i10,3x,f12.6)') i, ltarg(i), starg(i), etgr
    end do
  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_H_in, bug2, farm_format, molecule_format
    use scaling, only: set_charge, scale_radius, scale_etarg
    integer                 :: ihbuf(5)
    real(wp)                :: rhbuf(2), etgr
    integer                 :: i, l, status
    real(wp), allocatable   :: cfbut(:)
    real(wp)                :: e0
    real(wp)                :: r_mol, rmass ! UKRmol dummy parms, needed in prop K/T header output
    real(wp)                :: ezero !UKRmol  partitioned R-matrix parameter
    integer                 :: ibut, iex, nchan_choice ! UKRmol Buttle/partitioned input
    integer                 :: but_dim, iex_d  ! dummy for Buttle array, iex arrays
    real, allocatable       :: sfac(:),ecex(:), rcex(:,:) ! UKRmol partitioned R-matrix arrays


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

    if (xdr_H_in) 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)
       if (status /= 0) then 
         write(*,*) ' *** Allocation problem: Subroutine readhj1 *** '
         STOP
       end if 
       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)
       if (status /= 0) then 
         write(*,*) ' *** Allocation problem: Subroutine readhj1 *** '
         STOP
       end if 
       do l = 1, lrang2
          call xdr_io (iunit, cfbut, 3)   ! skip Buttle data for rmx1
       end do
    else          ! normal fortran unformatted output
       open (unit=iunit, file=TRIM(filh), status='old', action='read', iostat=status, form='unformatted')
       if(status.ne.0) write(*,*) ' Opening ', TRIM(filh), ' Error iostat = ', status
       if (status /= 0) then 
         write(*,*) ' *** Allocation problem: Subroutine readhj1 *** '
         STOP
       end if 
       read (iunit) nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
       ibut = 1
       if (molecule_format) read (iunit) r_mol, rmass, ezero, ibut, iex, nchan_choice
! note that for molecule_format, 'lrang2' is actually the number of channels
       but_dim = lrang2 * 3
       allocate (etarg(ntarg), jtarg(ntarg), ltarg(ntarg), starg(ntarg), &
            ptarg(ntarg), cfbut(but_dim),  stat=status)
       if (status /= 0) then 
         write(*,*) ' *** Allocation problem: Subroutine readhj1 *** '
         STOP
       end if 
       read (iunit) etarg(1:ntarg)
       if (farm_format) then
          read (iunit) jtarg(1:ntarg)
          read (iunit) starg(1:ntarg)
       else
          read (iunit) jtarg(1:ntarg)
          read (iunit) ltarg(1:ntarg)
          read (iunit) starg(1:ntarg)
          read (iunit) ptarg(1:ntarg)
       end if
! note slight cfbut inconsistency here for farm_format = .false. .
!This corresponds to ham95/fine 95 output as of January 2011.
       if (ibut == 1) then
          read (iunit) cfbut
          deallocate (cfbut, stat=status)
       else if (abs(ibut) == 2) then
! partitioned R-matrix, 'lrang2' = number of channels nchan
          deallocate (cfbut, stat=status)
          iex_d = MAX(iex,1)
          allocate (sfac(lrang2), ecex(iex_d), rcex(lrang2,iex_d), stat=status)
          if (status /= 0) then 
            write(*,*) ' *** Allocation problem sfac: Subroutine readhj1 *** '
            STOP
          end if 
          read (iunit) sfac(1:lrang2)
          if (iex /= 0)then
             read (iunit) ecex(1:iex)
             read (iunit) rcex(1:lrang2,1:iex) 
          end if 
          deallocate(rcex, ecex, sfac, stat=status)
       end if
       if (status /= 0) then 
          write(*,*) ' *** Allocation problem: Subroutine readhj1 *** '
          STOP
       end if
    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
       if (farm_format) then
          write (fo,'(a)') 'Farm binary format: etarg, jtarg, starg'
          write (fo,'(5f14.6)') etarg
          write (fo,'(5i6)') jtarg
          write (fo,'(5i6)') starg
       else
          write (fo,'(a)') 'Pfarm preferred xdr/binary format: etarg, jtarg, ltarg, starg, ptarg'
          write (fo,'(5f14.6)') etarg
          write (fo,'(5i6)') jtarg
          write (fo,'(5i6)') ltarg
          write (fo,'(5i6)') starg
          write (fo,'(5i6)') ptarg
       end if
    end if

! process input data:
    call set_charge (nz, nelc)
    call scale_radius (rmatr)
    call scale_etarg (etarg)     ! etarg now in scaled Ryd
    if(io_processor) write (fo,'(/,a,/)') 'Target states'
    e0 = etarg(1)
    if (farm_format) then
       if(io_processor) write (fo,'(10x,a,5x,a,3x,a,8x,a,/,43x,a)') 'index', &
        '  2j   ', ' (2*s+1)', 'energy', 'scaled ryd'
       do i = 1, ntarg
          etgr = (etarg(i) - e0)           ! convert to Rydbergs
          if(io_processor) write (fo,'(3x,3i10,3x,f12.6)') i, jtarg(i), &
                                       starg(i), etgr
       end do
    else
       if(io_processor) write (fo,'(10x,a,5x,a,3x,a,3x,a,8x,a,/,43x,a)') 'index', &
        '  2j   ', 'total l', '(2*s+1)', 'energy', 'scaled ryd'
       do i = 1, ntarg
          etgr = (etarg(i) - e0)           ! convert to Rydbergs
          if(io_processor) write (fo,'(3x,4i10,3x,f12.6)') i, jtarg(i), ltarg(i), &
                                       starg(i), etgr
       end do
    end if
  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)
  subroutine split_chls_old
! 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, status

    st2 = -999
    st1 = kschl(1)  ! 1st splitting parameter
    do i = 2, nc
       if (kschl(i) == st1) cycle
       st2 = kschl(i)
       exit
    end do
    allocate (schl1(nc), schl2(nc), stat=status)
    allocate (stl1(nc), stl2(nc), stat=status)
    if (status /= 0) then 
       write(*,*) ' *** Allocation problem st1: Subroutine get_split *** '
       STOP
    end if 
    if (st2 == -999) then
       do i = 1, nc
          schl1(i) = i
          stl1(i) = i
       end do
    end if
  end subroutine get_split

!  subroutine split_chls (schl1, schl2)
  subroutine split_chls
     use rmx1_in, only: bug1 
! 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
    stl1 = 0
    stl2 = 0
    schl1(1) = 1
    stl1(1) = 1
    channels: do i = 2, nc
       if (kschl(i) == st1) then ! first channel split
          nc1 = nc1 + 1
          schl1(nc1) = i
          stl1(i) = nc1
       else
          if (nts == 1) then   ! found a second ks value
             nts = nts + 1
             st2 = kschl(i)
             nc2 = 1
             schl2(nc2) = i
             stl2(i) = nc2
          else if (kschl(i) == st2) then ! second channel split
             nc2 = nc2 + 1
             schl2(nc2) = i
             stl2(i) = nc2
          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
    if (bug1 > 0) then 
      write (fo,*) 'Split Channels Spin 1 = ', (schl1(i),i=1,nc1)
      if (nc2 /= 0) write (fo,*) 'Split Channels Spin 2 = ', (schl2(i),i=1,nc2)
      write (fo,'(/)')
    end if
  end subroutine split_chls

  subroutine reorder_chls
     use rmx1_in, only: packed_cf
! order channel arrays lchl, tchl, kschl according to split parameters
    integer, allocatable      :: tmp_lchl(:), tmp_tchl(:), tmp_kschl(:)
    integer, allocatable      :: tmp_ilabel(:), tmp_jlabel(:)
    integer, allocatable      :: tmp_ic_label(:), tmp_ir_label(:)
    real(wp), allocatable     :: tmp_cf_nz(:), tmp_cf(:)
!    integer, allocatable      :: schl1(:), schl2(:)
    integer                   :: i, ip, jp, nc2, status, ctxt
    integer                   :: st_test, in, jn, ncf, icf, k, icf_1, icf_2, lama1

!    allocate (schl1(nc), schl2(nc), stat=status)
!    call split_chls (schl1, schl2)
    call split_chls
    if (nc1 == nc) then ! only one target spin for this SLp case
       deallocate (schl2, stl2, 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
    if (packed_cf) then
       lama1 = MAX(lamax,1)
       allocate (ncf2_nonzero(lama1), stat=status)
       allocate(tmp_ic_label(ncf_m), tmp_ir_label(ncf_m), tmp_cf_nz(ncf_m), stat=status) 
       allocate(tmp_ilabel(ncf_m), tmp_jlabel(ncf_m), tmp_cf(ncf_m), stat=status) 
! we don't know in advance how many non-zero items belong to nc1
       st_test = kschl(1)
       do k = 1, lamax 
          tmp_ic_label(:) = ic_label(:,k)
          tmp_ir_label(:) = ir_label(:,k)
          tmp_cf_nz(:) = cf_nz(:,k)
          ncf = ncf_nonzero(k)
          icf_1 = 0   ! split ncf into two sets
          icf_2 = 0
          do icf = 1, ncf
             jp = tmp_ic_label(icf)
             ip = tmp_ir_label(icf)  
             if (tmp_kschl(jp) /= tmp_kschl(ip)) then
                write(fo,*) 'problem with split packed coeficients:', &
                'ip, tmp_kschl(ip), jp, tmp_kschl(jp) =', ip, tmp_kschl(ip), jp, tmp_kschl(jp)
                write(fo,*) 'and lambda value = ', k
                status = 1
                STOP
             else if (tmp_kschl(jp) == st_test) then
                icf_1 = icf_1 + 1
                in = stl1(ip) 
                ir_label(icf_1,k) = in
                jn = stl1(jp)
                ic_label(icf_1,k) = jn
                cf_nz(icf_1,k) = tmp_cf_nz(icf)
             else 
                icf_2 = icf_2 + 1
                in = stl2(ip)  
                tmp_ilabel(icf_2) = in
                jn = stl2(jp) 
                tmp_jlabel(icf_2) = jn
                tmp_cf(icf_2) = tmp_cf_nz(icf)
             end if
          end do
          ncf_nonzero(k) = icf_1
          ncf2_nonzero(k) = ncf
          if (icf_1 + icf_2 /= ncf) then
             write(fo,*) 'problem with split packed coeficients:', &
             'icf_1, icf_2 =', icf_1, icf_2
             write(fo,*) 'and lambda value = ', k
             status = 1
             STOP
          end if
          ir_label(icf_1+1:icf_1+icf_2,k) = tmp_ilabel(1:icf_2)
          ic_label(icf_1+1:icf_1+icf_2,k) = tmp_jlabel(1:icf_2)
          cf_nz(icf_1+1:icf_1+icf_2,k) = tmp_cf(1:icf_2)
       end do
       deallocate (tmp_ic_label, tmp_ir_label, tmp_ilabel, tmp_jlabel, tmp_cf_nz, tmp_cf, stat=status)
    end if    
    deallocate (tmp_lchl, tmp_tchl, tmp_kschl, stat=status)
    if (status /= 0) then
       write(fo,*) 'allocation or deallocation failure in reorder_chls'
       STOP
    end if
  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_H_in, bug3, farm_format, nphys, packed_cf, large_wmat
    use scaling, only: scale_cf_farm, scale_cf_farm_pnz

    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
    integer                  :: nphych, nstart, lam
    integer                  :: icf, ncf, lama1

    more2 = 1
    match = .false.
    nt = ntarg
    symmetries: do while (more2 /= 0) ! search for required SLp
       if (xdr_H_in) 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)
          allocate (lchl(nc), nltarg(nt), tchl(nc), kschl(nc), &
                    eig(ns), wmat(nc), 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, end=20) lrgl2, nspn2, npty2, nc, ns, more2
          ncs = nc * ns
          allocate (lchl(nc), nltarg(nt), tchl(nc), kschl(nc), &
                    eig(ns), stat=status)
          read (iunit) nltarg
          read (iunit) lchl
          if (farm_format) then
             if (packed_cf) then
                lama1 = MAX(lamax,1)
                allocate(ncf_nonzero(lama1), stat=status)
                if (status /= 0) call s_error (status, 'ncf_nonzero in readh2: allocation error')
                read (iunit) ncf_nonzero(1:lama1)
                ncf_m = MAX(MAXVAL(ncf_nonzero),1)
                call flush(fo)
                allocate (cf_nz(ncf_m,lama1), ic_label(ncf_m,lama1), ir_label(ncf_m,lama1), stat=status)
                if (status /= 0) call s_error (status, 'cf_nz, label_nz in readh2: allocation error')
                cf_nz = 0.0_wp
                do k = 1, lamax
                   ncf = MAX(ncf_nonzero(k),1)
                   read (iunit) ir_label(1:ncf,k)
                   read (iunit) ic_label(1:ncf,k)
! ic: column (1, nc) ; ir: row (ic, nc)
                   read (iunit) cf_nz(1:ncf,k)
                end do
! note: scale_cf_farm_pnz assumes the cf_nz includes the factor of two
                call scale_cf_farm_pnz(cf_nz)
             else
                allocate (cf_farm(nc,nc,lamax), stat=status)
                if (status /= 0) call s_error (status, 'cf in readh2: allocation error')
                read (iunit) cf_farm
! note: scale_cf_farm assumes the cf_farm includes the factor of two
                call scale_cf_farm(cf_farm)
             end if
             read (iunit) eig
             if (large_wmat) then
                allocate (wmat(nc), stat=status)
                do j = 1, ns
                   read(iunit) loc
                   read(iunit) wmat
                end do
             else
                allocate (wmat(ncs), stat=status)
                read (iunit) wmat
             end if
             if (bug3 /= 0) then
                if (packed_cf) then
                   do k = 1, lamax
                      write(fo,*) 'non-zero cf_farm: k =', k
                      do icf = 1, ncf_nonzero(k)
                         j = ic_label(icf,k)
                         i = ir_label(icf,k)
                         write(fo,'((2i3,e14.6))') i, j, cf_nz(icf,k)
                      end do
                   end do  
                else
                   do k = 1, lamax
                      do j = 1, nc
                         write(fo,*) 'cf_farm: j,k =', j, k
                         write(fo,'(6(i3,e14.6))') (i, cf_farm(i,j,k), i = 1, nc)
                      end do
                   end do  
                end if
             end if
          else
             allocate (wmat(nc), stat=status)
             if (status /= 0) call s_error (status, 'wmat in readh2: allocation error')
             read (iunit) eig
             do j = 1, ns
                read (iunit) loc
                read (iunit) wmat(1:nc)
             end do
          end if
       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.
! set 'unphysical' channel pot-cofs to zero
          if (farm_format .and. nphys /= 0 .and. nphys < nc) then
             nphych = 0
             do it  = 1, nphys
                nphych = nphych + nltarg(it)
             enddo
             if (packed_cf) then
                do lam = 1, lamax
                   do icf = 1, ncf_nonzero(lam)
                      i = ic_label(icf,lam)
                      j = ir_label(icf,lam)
                      if ((i > nphych) .or. (j > nphych)) cf_nz(icf,lam) = 0.0_wp
                   end do
                end do
             else
                do i = 1, nc
                   if (i.gt.nphych) then 
                      nstart = 1
                   else
                      nstart = nphych + 1
                   endif
                   do j = nstart, nc 
                      do lam = 1,lamax
                         cf_farm(j,i,lam) = 0.0
                      end do
                   end do
                end do
!                write (99,1234) nphych
!                do lam = 1,lamax
!                it = 0
!                do i = 1,nc
!                   do j = 1,nc
!                     it = it + 1
!                     write (99,1235) lam,i,j,it,cf_farm(j,i,lam)
!                   end do
!                end do
!                end do
 1234           format (' asymptotic coefficients set to zero beyond nphych =', &
     &           I5) 
 1235           format (' j, i1, i2, ii cf(ii,j) = ',4I8,E14.7)
             end if 
          end if
          deallocate (eig, wmat, nltarg, stat=status)
          exit
       end if
       deallocate (lchl, tchl, kschl, eig, wmat, nltarg, stat=status)
       if (farm_format) then
          if (packed_cf) then
             deallocate (cf_nz, ic_label, ir_label, ncf_nonzero, stat=status)
          else
             deallocate (cf_farm, stat=status)
          end if
       end if
! why was this line here?
!       if (more2 == 2) exit  ! too much expected of the H file(?)
!       exit
! a proper farm/H file will have more2 = 0 for the final symmetry
    end do symmetries

 20 continue
    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
       if (farm_format .and. nphys /= 0 .and. nphys < nc) &
           write(*,'(a,i5)') &
           & 'Asymptotic coefficients set to zero beyond nphych =', nphych
    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_H_in, bug3, farm_format, nphys, packed_cf, large_wmat
    use scaling, only: scale_cf_farm, scale_cf_farm_pnz

    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
    integer                  :: nphych, nstart, lam
    integer                  :: icf, ncf, lama1

    more2 = 1
    match = .false.
    nt = ntarg
    symmetries: do while (more2 /= 0) ! search for required SLp
       if (xdr_H_in) 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(nc), stat=status)

!          write(fo,*) ' about to read xdr'
          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
!         write(fo,*) 'xdr read done'
       else
          read (iunit, end=20) lrgl2, nspn2, npty2, nc, ns, more2
          write (fo,*) 'lrgl2, nspn2, npty2, nc, ns, more2', &
                        lrgl2, nspn2, npty2, nc, ns, more2
          ncs = nc * ns
          allocate (lchl(nc), kschl(nc), nltarg(nt), tchl(nc), eig(ns),& 
                   stat=status)

          read (iunit) nltarg
          write(*,*) ' nltarg read'          
          read (iunit) lchl
          write(*,*) 'lchl read'
          if (farm_format) then
             kschl(1:nc) = 0   ! dummy, just need all the same value (?)
             if (packed_cf) then
                lama1 = MAX(lamax,1)
                allocate(ncf_nonzero(lama1), stat=status)
                if (status /= 0) call s_error (status, 'ncf_nonzero in readh2: allocation error')
                read (iunit) ncf_nonzero(1:lama1)
                ncf_m = MAX(MAXVAL(ncf_nonzero),1)
                allocate (cf_nz(ncf_m,lama1), ic_label(ncf_m,lama1), ir_label(ncf_m,lama1), stat=status)
                if (status /= 0) call s_error (status, 'cf_nz, label_nz in readh2: allocation error')
                cf_nz = 0.0_wp
!                write(fo,*) ' about to read, lamax =', lamax
                do k = 1, lamax
                   ncf = MAX(ncf_nonzero(k),1)
                   read (iunit) ir_label(1:ncf,k)
                   read (iunit) ic_label(1:ncf,k)
! ic: column (1, nc) ; ir: row (ic, nc)
!                    write(fo,*) 'labels read, k =', k
                   read (iunit) cf_nz(1:ncf,k)
!                   write(fo,*) 'cf read, k = ', k
                end do
! note: scale_cf_farm_pnz assumes the cf_nz includes the factor of two
                call scale_cf_farm_pnz(cf_nz)
!                  write(fo,*) 'done scale'
             else
                allocate (cf_farm(nc,nc,lamax), stat=status)
                if (status /= 0) call s_error (status, 'cf in readhj2: allocation error')
                read (iunit) cf_farm
                write(*,*) 'cf_farm read'
! note: scale_cf_farm assumes the cf_farm includes the facotr of two
                call scale_cf_farm(cf_farm)
             end if
             read (iunit) eig
                write(*,*) 'eig read'
             if (large_wmat) then
                allocate (wmat(nc), stat=status)
                do j = 1, ns
                   read(iunit) loc
                   read(iunit) wmat
                end do
             else
                allocate (wmat(ncs), stat=status)
                read (iunit) wmat
                write(*,*) 'wmat read'
             end if
             if (bug3 /= 0) then
                if (packed_cf) then
                   do k = 1, lamax
                      write(fo,*) 'non-zero cf_farm: k =', k
                      do icf = 1, ncf_nonzero(k)
                         j = ic_label(icf,k)
                         i = ir_label(icf,k)
                         write(fo,'((2i3,e14.6))') i, j, cf_nz(icf,k)
                      end do
                   end do  
                else
                   do k = 1, lamax
                      do j = 1, nc
                         write(fo,*) 'cf_farm: j,k =', j, k
                         write(fo,'(6(i3,e14.6))') (i, cf_farm(i,j,k), i = 1, nc)
                      end do
                   end do  
                end if
             end if
 !             write(fo,*) 'read finished' 
          else
             read (iunit) kschl
             allocate (wmat(ncs), stat=status)
             if (status /= 0) call s_error (status, 'wmat in readh2: allocation error')
             read (iunit) eig
             do j = 1, ns
                read (iunit) loc
                read (iunit) wmat(1:nc)
             end do
          end if
       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
!       write(fo,*) 'done tchl'

       if (lrgl1 == lrgl2 .and. npty1 == npty2) then !nspn1 not defined by namelist
          match = .true.
! set 'unphysical' channel pot-cofs to zero
          if (farm_format .and. nphys /= 0 .and. nphys < nc) then
             nphych = 0
             do it  = 1, nphys
                nphych = nphych + nltarg(it)
             enddo
             if (packed_cf) then
                do lam = 1, lamax
                   do icf = 1, ncf_nonzero(lam)
                      i = ic_label(icf,lam)
                      j = ir_label(icf,lam)
                      if ((i > nphych) .or. (j > nphych)) cf_nz(icf,lam) = 0.0_wp
                   end do
                end do
             else
                do  i = 1, nc
                   if (i.gt.nphych) then 
                      nstart = 1
                   else
                      nstart = nphych + 1
                   endif
                   do j = nstart, nc 
                      do lam = 1,lamax
                         cf_farm(j,i,lam) = 0.0
                      end do
                   end do
                end do
             end if 
          end if 
          deallocate (eig, wmat, nltarg, stat=status)
 !      write(fo,*) 'first deallocates done'
          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 (farm_format) then
          if (packed_cf) then
             deallocate (cf_nz, ic_label, ir_label, ncf_nonzero, stat=status)
          else
             deallocate (cf_farm, stat=status)
          end if  
       end if     
! see comments for readh2
!       if (more2 == 2) exit
!       exit
    end do symmetries

 20 continue
    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
       if (farm_format) then 
          write (fo,'(/,a)') 'Channel data (farm binary format):'
          write (fo,'(a)') '  target index       l'
          do i = 1, nc
             write (fo,'(i4,6x,i4,5x,i4,5x,i4)') i, tchl(i), lchl(i)
          end do
       else
          write (fo,'(/,a)') 'Channel data (pfarm preferred xdr/binary format):'
          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
          if (farm_format .and. nphys /= 0 .and. nphys < nc) &
              write(*,'(a,i5)') &
              & 'Asymptotic coefficients set to zero beyond nphych =', nphych
       end if
    end if
  end subroutine readhj2

  subroutine s_error (err, msg)
    use io_units, only: fo
    integer, intent(in)          :: err ! flag
    character(len=*), intent(in) :: msg   ! error message

    if (err /= 0) write (fo, fmt='(a)') msg

  end subroutine s_error

end module hfile_data
