module potl_cofs
!  Generate asymptotic coefficients from target multipoles
! Time-stamp: "2003-09-23 09:42:38 cjn"
! modified for jpi_coupling april to dec 2006 vmb
  use precisn, only: wp
  use io_units, only: fo
  implicit none

  private
  public cfs, cfsj, cf, lamd, del_cf, rdtmom, rdjtmom

  logical, save      :: gams = .false.
  integer, parameter :: ngam = 500             ! n entries in gam table
  real(wp), save     :: gam(ngam)              ! gam table

  type tmblock    ! transition moment type definition
     sequence
     integer           :: spn         ! spin value
     integer           :: ntm         ! # transition moments
     integer, pointer  :: lmds(:)     ! multipole index
     integer, pointer  :: lbls(:)     ! locator array
     real(wp), pointer :: tms(:)     ! transition moments
  end type tmblock

  integer, save                    :: ntmb    ! # target spins
  type(tmblock), allocatable, save :: tmb(:)   ! target transition moments

  real(wp), allocatable, save  :: cf(:,:)
  integer, allocatable, save   :: lamd(:,:)

  integer, save                :: nt    ! # targets
  integer, allocatable, save   :: tj(:), tl(:), tspn(:), tpty(:)
  real(wp), allocatable, save  :: tm(:,:) 

contains

  subroutine cfs (nc, tchl, lchl, sqno, spins)
! generate scaled asymptotic potential coefficients
    use slp, only: qno, val_qno
    use hfile_data, only: lamax
!    use blacs, only: ctxt, p_error
    use scaling, only: scale_cf
    use error_out, only: error_check

    integer, intent(in)    :: nc  ! no channels
    integer, intent(in)    :: tchl(nc) ! target sequence of channel
    integer, intent(in)    :: lchl(nc) ! orb. a.m. of channel
    type(qno), intent(in)  :: sqno     ! scattering a. #s
    integer, intent(in)    :: spins(2) ! target spins required
    real(wp)               :: phz, z, drac
    integer                :: i, j, tari, tarj, li, lj, locc, tt, loct, tip
    integer                :: lbx, lbn, lb, lbb, tjl, tjs, tjp, til, tis
    integer                :: nspnt,  qn(3), lrgl, nspn, npty, isp
    integer                :: status, smx

    smx = lamax
    allocate (cf(nc*(nc+1)/2,(smx+1)/2), lamd(nc*(nc+1)/2,(smx+1)/2),&
         stat=status)
    call error_check (status, 'cfs: allocation error')
    cf = 0.0_wp
    lamd = 0

    target_spins: do isp = 1, 2 ! 1 spin value or 2
       nspnt = spins(isp)
       if (nspnt == -999) cycle ! only one spin value
       qn = val_qno(sqno)
       lrgl = qn(1)   ! scattering orb am
       nspn = qn(2)   ! scattering spin multiplicity
       npty = qn(3)   ! scattering parity

       call get_tm (nspnt, smx, nt)

       phz = 1
       if (MOD(lrgl+npty,2) == 1) phz = -1
       do j = 1, nc
          tarj = tchl(j)
          if (tspn(tarj) /= nspnt) cycle
          lj = lchl(j)
          tjl = tl(tarj)
          tjs = tspn(tarj)
          tjp = tpty(tarj)
          locc = j*(j-1)/2
          do i = 1, j
             tari = tchl(i)
             if (tspn(tari) /= nspnt) cycle
             li= lchl(i)
             til = tl(tari)
             tis = tspn(tari)
             tip = tpty(tari)
             if (tari > tarj) then ! ensure input target sequencing
                tt = tarj
                tarj = tari
                tari = tt
             end if
             loct = tarj * (tarj - 1) / 2 + tari
             lbx = MIN(li+lj,til+tjl,smx)    ! min lambda
             lbn = MAX(ABS(li-lj),ABS(til-tjl),1) ! max lambda
             do lb = lbn, lbx, 1
                if (MOD(tjp+tip+lb,2) == 1) cycle
                lbb = (lb + 1) / 2
                if (ABS(tm(loct,lbb)) <= 1.0e-6_wp) cycle
                drac = dracah(2*til, 2*tjl, 2*li, 2*lj, 2*lb, 2*lrgl)
                z = phz * drac * rme(li, lj, lb)
                cf(locc+i,lbb) = tm(loct,lbb) * z
                lamd(locc+i,lbb) = lb
             end do
          end do
       end do
       deallocate (tm, stat=status)
       if (status /= 0) STOP
    end do target_spins
    call scale_cf (cf, lamd)
  end subroutine cfs

  subroutine cfsj (nc, tchl, lchl, kschl, sqno, spins)
! generate scaled asymptotic potential coefficients
    use slp, only: qno, val_qno
    use hfile_data, only: lamax
!    use blacs, only: ctxt, p_error
    use scaling, only: scale_cf
    use error_out, only: error_check
    integer, intent(in)    :: nc  ! no channels
    integer, intent(in)    :: tchl(nc) ! target sequence of channel
    integer, intent(in)    :: lchl(nc) ! orb. a.m. of channel
    integer, intent(in)    :: kschl(nc) ! k q. no
    type(qno), intent(in)  :: sqno     ! scattering q. #s
    integer, intent(in)    :: spins(2) ! k splits required
    real(wp)               :: phz, z, drac
    integer                :: i, j, tari, tarj, li, lj, locc, tt, loct, tip
    integer                :: lbx, lbn, lb, lbb, tjl, tjs, tjp, til, tis
    integer                :: ki, kj, kbar, tij, tjj
    integer                :: kspnt,  qn(3), lrgl, nspn, npty, isp
    integer                :: status, smx

    smx = lamax
    allocate (cf(nc*(nc+1)/2,(smx+1)/2), lamd(nc*(nc+1)/2,(smx+1)/2),&
         stat=status)
    call error_check (status, 'cfs: allocation error')
    cf = 0.0_wp
    lamd = 0

    channel_k_splits: do isp = 1, 2 ! 1 k value or 2
       kspnt = spins(isp) ! value of ksplit = 2k
       if (kspnt == -999) cycle ! only one k value
       qn = val_qno(sqno)
       lrgl = qn(1)   ! scattering orb am = 2J
       nspn = qn(2)   ! scattering spin multiplicity = 0
       npty = qn(3)   ! scattering parity

       call get_tm (0, smx, nt) ! target moments are not split
       phz = 1
       if (MOD(kspnt, 2) == 0) kbar=kspnt/2
       if (MOD(kspnt, 2) == 1) kbar=(kspnt+1)/2! kbar remains real
       if (MOD(npty+kbar, 2).eq.1) phz = -1
       do j = 1, nc
          kj = kschl(j)
          if (kj /= kspnt) cycle
          tarj = tchl(j)
          lj = lchl(j)
          tjj = tj(tarj)
          tjl = tl(tarj)
          tjs = tspn(tarj)
          tjp = tpty(tarj)
          locc = j*(j-1)/2
          do i = 1, j
             tari = tchl(i)
             li = lchl(i)
             ki = kschl(i)
             if (ki /= kj) cycle
             tij = tj(tari)
             til = tl(tari)
             tis = tspn(tari)
             tip = tpty(tari)
             if (tari > tarj) then ! ensure input target sequencing
                tt = tarj
                tarj = tari
                tari = tt
             end if
             loct = tarj * (tarj - 1) / 2 + tari
             lbx = MIN(li+lj,smx)    ! max lambda
             lbn = MAX(ABS(li-lj),1) ! min lambda
             do lb = lbn, lbx, 1
                if (MOD(tjp+tip+lb,2) == 1) cycle
                lbb = (lb + 1) / 2
                if (ABS(tm(loct,lbb)) <= 1.0e-6_wp) cycle
                drac = dracah(tij, tjj, 2*li, 2*lj, 2*lb, kspnt)
                z = phz * drac * rme(li, lj, lb)
                cf(locc+i,lbb) = tm(loct,lbb) * z
                lamd(locc+i,lbb) = lb
             end do
          end do
       end do
       deallocate (tm, stat=status)
       if (status /= 0) STOP
    end do channel_k_splits
    call scale_cf (cf, lamd)
  end subroutine cfsj

  subroutine del_cf
!    use blacs, only: p_error
    use error_out, only: error_check
    integer        :: status
    deallocate (cf, lamd, stat=status)
    call error_check (status, 'del_cf: deallocation error')
  end subroutine del_cf

  function dracah (i, j, k, l, m, n)
! racah coefficients
    real(wp)           :: dracah
    integer, intent(in) :: i, j, k, l, m, n ! input momenta * 2
    integer :: j1, j2, j3, j4, j5, j6, j7, numin, numax, icount, kk, ki

    if (.not.gams) then
       call gam_tab
       gams = .true.
    end if
    j1 = i + j + m
    j2=  k + l + m
    j3 = i + k + n
    j4 = j + l + n
    if (MOD(j1,2) /= 0 .or. MOD(j2,2) /= 0 .or. MOD(j3,2) /= 0 .or. &
         MOD(j4,2) /= 0) then
       dracah = 0.0_wp
       return
    end if
    if (2*MAX(i,j,m) > j1 .or. 2*MAX(k,l,m) > j2 .or. 2*MAX(i,k,n) > j3 &
         .or. 2*MAX(j,l,n) > j4) then
       dracah = 0.0_wp
       return
    end if
    j1 = j1 / 2
    j2 = j2 / 2
    j3 = j3 / 2
    j4 = j4 / 2
    j5 = (i + j + k + l) / 2
    j6 = (i + l + m + n) / 2
    j7 = (j + k + m + n) / 2
    numin = MAX(j1, j2, j3, j4) + 1
    numax = MIN(j5, j6, j7) + 1
    dracah = 1.0_wp
    icount = 0
    if (numin < numax) then
       do kk = numin+1, numax
          ki = numax - icount
          dracah = 1.0_wp - (dracah * real((ki * (j5-ki+2) * (j6-ki+2) * &
               (j7-ki+2)),wp) / real((ki-1-j1) * (ki-1-j2) * (ki-1-j3) * &
               (ki-1-j4),wp))
          icount = icount + 1
       end do
    end if
    dracah = dracah * (-1.0_wp)**(j5+numin+1) *                     &
         EXP((gam(numin+1) - gam(numin-j1) - gam(numin-j2) -         &
         gam(numin-j3) - gam(numin-j4) - gam(j5+2-numin) -           &
         gam(j6+2-numin) - gam(j7+2-numin)) + ((gam(j1+1-i) +        &
         gam(j1+1-j) + gam(j1+1-m) - gam(j1+2) + gam(j2+1-k) +       &
         gam(j2+1-l) + gam(j2+1-m) - gam(j2+2) + gam(j3+1-i) +       &
         gam(j3+1-k) + gam(j3+1-n) - gam(j3+2) + gam(j4+1-j) +       &
         gam(j4+1-l) + gam(j4+1-n) - gam(j4+2)) / 2.0_wp))
  end function dracah

  function rme (l, lp, k)
! reduced matrix element (l//c(k)//lp)
! Ref: Fano and Racah, Irreducible Tensorial Sets, Ch 14, P. 81
    real(wp)           :: rme
    integer, intent(in) :: l, lp           ! initial, final orbital a.m.
    integer, intent(in) :: k               ! order of spherical tensor
    integer             :: i2g, ig, i1, i2, i3, imax, ier
    real(wp)           :: qusqrt

    if (.not.gams) then
       call gam_tab
       gams = .true.
    end if
    if (k < ABS(l-lp) .or. k > l+lp) then
       write (fo,'(a,/,a,2i3)') 'rme: a.m. &
            & coupling error', 'rme zeroed for l, lp, k = ', l, lp, k
       rme = 0.0_wp
       return
    end if
    i2g = l + lp + k
    ig = i2g / 2
    if (MOD(i2g,2) /= 0) then
       rme = 0.0_wp
       return
    end if
    if (ig < 0) then
       write (fo,'(a,/,a,2i3)') 'rme: &
            & angle does not match', &
            'rme zeroed for l, lp, k = ', l, lp, k
       rme = 0.0_wp
       return
    elseif (ig == 0) then
       rme = 1.0_wp
       return
    end if
    i1 = ig - l
    i2 = ig - lp
    i3 = ig - k
    imax = MAX(ig+1, i2g+2, 2*i3+1, 2*i2+1, 2*i1+1)
    if (imax > ngam) then
       write (fo,'(a,2i6)') 'rme: table overflow, imax, ngam = ', &
            imax, ngam
       STOP
    end if
    qusqrt = log(real(2*l+1,wp)) + log(real(2*lp+1,wp)) + gam(2*i1+1) +&
         gam(2*i2+1) + gam(2*i3+1) - gam(i2g+2)
    rme = EXP(0.5_wp * qusqrt + gam(ig+1) - gam(i1+1) - gam(i2+1) -     &
         gam(i3+1))
  end function rme

  subroutine gam_tab
! logs of factorials table
    real(wp) :: x
    integer   :: i

    gam(1) = 1.0_wp
    gam(2) = 1.0_wp
    x = 2.0_wp
    do i = 3, 25
       gam(i) = gam(i-1) * x
       x = x + 1.0_wp
    end do
    gam(1:25) = LOG(gam(1:25))
    x = 25.0_wp
    do i = 26, ngam
       gam(i) = gam(i-1) + LOG(x)
       x = x + 1.0_wp
    end do
  end subroutine gam_tab

  subroutine rdtmom
    ! read file of target moments
    !    use blacs, only: ictxt, p_error, io_processor
    use io_units, only: fo, itm
    use xdr_files, only: xdr_io, open_xdr, close_xdr
    use rmx1_in, only: bug4, xdr, fltm
    use error_out, only: error_check
    real(wp), allocatable           :: tmom(:,:)
    integer, allocatable            :: ntmom(:)
    integer                         :: ibuf(2)
    integer                         :: lamind, isp, ijl, ist, jst
    integer                         :: lamlo, lamup, lam, lamn, nijl
    integer                         :: ltmom, status, ios, s, i, j
    integer                         :: nsv, spn, min_lam, max_lam, ntm
    integer                         :: min_spn, max_spn, itmom, lamax, ip

    if (xdr) then    ! output file used XDR format
       itmom = open_xdr (file=TRIM(fltm), action='read')
       call xdr_io (itmom, ibuf, 2)
       nt = ibuf(1)
       lamax = ibuf(2)
       allocate (tl(nt), tspn(nt), tpty(nt), stat=status)
       if (status /= 0) then
          write (fo,'(a,i6)') 'rdtmom: allocation error = ', status
          STOP
       end if
       call xdr_io (itmom, tl, nt)
       call xdr_io (itmom, tspn, nt)
       call xdr_io (itmom, tpty, nt)
       call xdr_io (itmom, ibuf, 2) ! upper, lower limits tgt spin
       min_spn = ibuf(1)
       max_spn = ibuf(2)
    else             ! fortran unformatted output file
       itmom = itm
       open (unit=itmom, file=TRIM(fltm), form='unformatted', &
            status='old', action='read', iostat=ios)
       if (ios /= 0) then
          write (fo,'(a,i6)') 'rdtmom: error opening TARMOM, &
               &iostat = ', ios
          stop
       end if
       read (itmom) nt, lamax
       allocate (tl(nt), tspn(nt), tpty(nt), stat=status)
       call error_check (status, 'rdtmom: allocation error(1)')
       read (itmom) tl(1:nt)
       read (itmom) tspn(1:nt)
       read (itmom) tpty(1:nt)
       read (itmom) min_spn, max_spn
    end if

    lamind = (lamax + 1) / 2
    ltmom = (nt*(nt+1)/2) * ((lamax+1)/2)
    nsv = (max_spn - min_spn) / 2 + 1
    allocate (tmom(ltmom,nsv), ntmom(nsv), stat=status)
    if (status /= 0) then
       write (fo,'(a,i6)') 'rdtmom: allocation error = ', status
       STOP
    end if
    if (bug4 > 0) then ! echo data read from TARMOM
       write (fo,'(a,i4,a,i4)') 'No. of states =', nt, &
            ' max lambda =', lamax
       write (fo,'(a)') 'Target L:'
       write (fo,'(25i3)') tl(1:nt)
       write (fo,'(a)') 'Target S:'
       write (fo,'(25i3)') tspn(1:nt)
       write (fo,'(a)') 'Target P:'
       write (fo,'(25i3)') tpty(1:nt)
       write (fo,'(a)') 'Multiplicity limits:'
       write (fo,'(25i3)') min_spn, max_spn
    end if

    s = 0
    spins: do spn = min_spn, max_spn, 2
       s = s + 1
       if (xdr) then
          call xdr_io (itmom, ibuf, 2)
          isp = ibuf(1)
          ntm = ibuf(2)
          ntmom(s) = ntm
          call xdr_io (itmom, tmom(1:ntm,s), ntm)
       else
          read (itmom) isp, ntm
          ntmom(s) = ntm
          read (itmom) tmom(1:ntm,s)
       end if
       if (bug4 > 0) then
          write (fo,'(a,i4)') 'Target spin = ', isp
          write (fo,'(a,i6)') 'Number of moments = ', ntm
          write (fo,'(5f14.6)') tmom(1:ntm,s)
       end if
    end do spins

    if (xdr) then
       call close_xdr (itmom)
    else
       close (unit=itmom, status='keep')
    end if

    ! unpack tmom:
    ntmb = nsv
    allocate (tmb(nsv), stat=status)
    call error_check (status, 'rdtmom: tmb allocation error')
    s = 0
    do spn = min_spn, max_spn, 2
       s = s + 1
       ntm = ntmom(s)
       allocate (tmb(s)%lbls(ntm), tmb(s)%lmds(ntm), tmb(s)%tms(ntm), &
            stat=status)
       call error_check (status, 'rdtmom: allocation error(5)')
       tmb(s)%spn = spn
       tmb(s)%ntm = ntm
       ijl = 0
       do i = 1, nt
          if (tspn(i) /= spn) cycle
          ip = i * (i - 1) / 2
          do j = 1, i
             if (tspn(j) /= spn) cycle
             min_lam = ABS(tl(i) - tl(j))
             if (MOD(tpty(i)+tpty(j)+min_lam,2) == 1) min_lam = &
                  min_lam + 1
             if (min_lam == 0) min_lam = 2
             max_lam = MIN(tl(i)+tl(j), lamax)
             if ((min_lam+1)/2 <= (max_lam+1)/2) then
                do lam = min_lam, max_lam, 2
                   ijl = ijl + 1
                   tmb(s)%lbls(ijl) = ip + j
                   tmb(s)%lmds(ijl) = lam
                   tmb(s)%tms(ijl) = tmom(ijl,s)
                end do
             end if
          end do
       end do
    end do
    deallocate (tmom, ntmom, stat=status)
    call error_check (status, 'rdtmom: deallocation error')
  end subroutine rdtmom

  subroutine rdjtmom
    ! read file of target moments when jpi_coupling is true
    !    use blacs, only: ictxt, p_error, io_processor
    use io_units, only: fo, itm
    use xdr_files, only: xdr_io, open_xdr, close_xdr
    use rmx1_in, only: bug4, xdr, fltm
    use error_out, only: error_check

    real(wp), allocatable           :: tmom(:,:)
    integer, allocatable            :: ntmom(:)
    integer                         :: ibuf(2)
    integer                         :: lamind, isp, ijl, ist, jst
    integer                         :: lamlo, lamup, lam, lamn, nijl
    integer                         :: ltmom, status, ios, s, i, j
    integer                         :: nsv, spn, min_lam, max_lam, ntm
    integer                         :: min_spn, max_spn, itmom, lamax, ip

    if (xdr) then    ! output file used XDR format
       itmom = open_xdr (file=TRIM(fltm), action='read')
       call xdr_io (itmom, ibuf, 2)
       nt = ibuf(1)
       lamax = ibuf(2)
       allocate (tj(nt), tl(nt), tspn(nt), tpty(nt), stat=status)
       if (status /= 0) then
          write (fo,'(a,i6)') 'rdjtmom: allocation error = ', status
          STOP
       end if
       call xdr_io (itmom, tj, nt)
       call xdr_io (itmom, tl, nt)
       call xdr_io (itmom, tspn, nt)
       call xdr_io (itmom, tpty, nt)
       call xdr_io (itmom, ibuf, 2) ! upper, lower limits tgt spin
       min_spn = ibuf(1)
       max_spn = ibuf(2)
    else             ! fortran unformatted output file
       itmom = itm
       open (unit=itmom, file=TRIM(fltm), form='unformatted', &
            status='old', action='read', iostat=ios)
       if (ios /= 0) then
          write (fo,'(a,i6)') 'rdjtmom: error opening TARMOM, &
               &iostat = ', ios
          stop
       end if
       read (itmom) nt, lamax
       allocate (tj(nt), tl(nt), tspn(nt), tpty(nt), stat=status)
       call error_check (status, 'rdjtmom: allocation error(1)')
       read (itmom) tj(1:nt)
       read (itmom) tl(1:nt)
       read (itmom) tspn(1:nt)
       read (itmom) tpty(1:nt)
       read (itmom) min_spn, max_spn !set zero in fine.f90
    end if

    lamind = (lamax + 1) / 2
    ltmom = (nt*(nt+1)/2) * ((lamax+1)/2)
    nsv = (max_spn - min_spn) / 2 + 1 !nsv=1 here
    allocate (tmom(ltmom,nsv), ntmom(nsv), stat=status)
    if (status /= 0) then
       write (fo,'(a,i6)') 'rdjtmom: allocation error = ', status
       STOP
    end if
    if (bug4 > 0) then ! echo data read from TARMOM
       write (fo,'(a,i4,a,i4)') 'No. of states =', nt, &
            ' max lambda =', lamax
       write (fo,'(a)') 'Target 2J:'
       write (fo,'(25i3)') tj(1:nt)
       write (fo,'(a)') 'Target L:'
       write (fo,'(25i3)') tl(1:nt)
       write (fo,'(a)') 'Target S:'
       write (fo,'(25i3)') tspn(1:nt)
       write (fo,'(a)') 'Target P:'
       write (fo,'(25i3)') tpty(1:nt)
       write (fo,'(a)') 'Multiplicity limits:'
       write (fo,'(25i3)') min_spn, max_spn
    end if

    s = 0
    spins: do spn = min_spn, max_spn, 2! once thro loop
       s = s + 1
       if (xdr) then
          call xdr_io (itmom, ibuf, 2)
          isp = ibuf(1)
          ntm = ibuf(2)
          ntmom(s) = ntm
          call xdr_io (itmom, tmom(1:ntm,s), ntm)
       else
          read (itmom) isp, ntm
          ntmom(s) = ntm
          read (itmom) tmom(1:ntm,s)
       end if
       if (bug4 > 0) then
          write (fo,'(a,i4)') 'Notional Target spin = ', isp
          write (fo,'(a,i6)') 'Number of moments = ', ntm
          write (fo,'(5f14.6)') tmom(1:ntm,s)
       end if
    end do spins

    if (xdr) then
       call close_xdr (itmom)
    else
       close (unit=itmom, status='keep')
    end if

    ! unpack tmom:
    ntmb = nsv
    allocate (tmb(nsv), stat=status)
    call error_check (status, 'rdjtmom: tmb allocation error')
    s = 0
    do spn = min_spn, max_spn, 2
       s = s + 1
       ntm = ntmom(s)
       allocate (tmb(s)%lbls(ntm), tmb(s)%lmds(ntm), tmb(s)%tms(ntm), &
            stat=status)
       call error_check (status, 'rdtmom: allocation error(5)')
       tmb(s)%spn = spn
       tmb(s)%ntm = ntm
       ijl = 0
       do i = 1, nt
          ip = i * (i - 1) / 2
          do j = 1, i
             !             min_lam = ABS(tl(i) - tl(j))
             !             if (MOD(tpty(i)+tpty(j)+min_lam,2) == 1) min_lam = &
             !                  min_lam + 1
             min_lam = MOD(tpty(i)+tpty(j),2)
             if (min_lam == 0) min_lam = 2
             !             max_lam = MIN(tl(i)+tl(j), lamax)
             max_lam = lamax
             if ((min_lam+1)/2 <= (max_lam+1)/2) then
                do lam = min_lam, max_lam, 2
                   ijl = ijl + 1
                   tmb(s)%lbls(ijl) = ip + j
                   tmb(s)%lmds(ijl) = lam
                   tmb(s)%tms(ijl) = tmom(ijl,s)
                end do
             end if
          end do
       end do
    end do
    deallocate (tmom, ntmom, stat=status)
    call error_check (status, 'rdjtmom: deallocation error')
  end subroutine rdjtmom

  subroutine get_tm (spn, smx, nt)
! obtain expanded array of transition moments, tm for spin spn
    use rmx1_in, only: bug5
    integer, intent(in)       :: spn       ! target spin
    integer, intent(in)       :: smx       ! max multipole
    integer, intent(in)       :: nt        ! # targets
    integer                   :: i, j, lp, ntm, l, ier, status, lb

    do i = 1, ntmb
       if (.not.associated(tmb(i)%tms)) exit
       if (tmb(i)%spn /= spn) cycle ! not spin required
       ntm = tmb(i)%ntm      ! # transition moments for spin value
       allocate (tm(nt*(nt+1)/2,(smx+1)/2), stat=status)
       if (status /= 0) STOP
       tm = 0.0_wp
       do j = 1, ntm
          l = tmb(i)%lmds(j)
          lb = tmb(i)%lbls(j)
          tm(lb, (l+1)/2) = tmb(i)%tms(j)
       end do
       if (bug5 > 0) then
          do j = 1, ntm
             l = tmb(i)%lmds(j)
             lb = tmb(i)%lbls(j)
             write (fo,'(3i6,a,f14.6)') j, lb, l, ' tms = ', tmb(i)%tms(j)
          end do
       end if
       return  ! required tm set has been located
    end do
    write (fo,'(a,i6)') 'get_tm: spin not found, spn = ', spn
    STOP
  end subroutine get_tm

end module potl_cofs
