module rw_amplitudes
! calculated sector reduced-width amplitudes
! Time-stamp: "2003-10-06 14:23:21 cjn"
  use io_units, only: fo, fa
  use precisn, only: wp
  use def_sectors, only: nh, nsect
  use mpi_params, only: taskid, numtasks, io_processor
  implicit none

  integer, save     :: ixdro   ! XDR unit number for rwa file

  private
  public start_sa_file, close_sa_file, ampltd

contains


  subroutine start_sa_file (sqno, st, sect)
! open XDR file to hold sector reduced width amplitudes
    use xdr_files, only: open_xdr, xdr_io
    use energy_grid, only: fine_mesh
    use slp, only: qno, val_qno
    use rmx1_in, only: xdr, split_prop
    type(qno), intent(in)   :: sqno        ! scattering q#s
    integer, intent(in)     :: st         ! second target spin
    integer, intent(in)     :: sect ! sector id
    character(len=11)       :: hname
    character(len=2)        :: sp, l
    character(len=1)        :: p
    character(len=4)        :: stem
    character(len=2)        :: isect
    integer                 :: ios, qn(3), status
    integer                 :: nspn ! scattering spin.  

! File naming convention:
!  AMPCSSLLP: coarse mesh, 1st partition (or all channels)
!  AMPFSSLLP: fine mesh, 1st partition (or all channels)
!  AMPCTTLLP: coarse mesh, 2nd partition. TT = target spin multiplicity
!  AMPFTTLLP: fine mesh, 2nd partition, TT = target spin multiplicity
! For no-exchange LS files target spin multiplicity is used
! For J files L is replaced by 2*J and total spin multiplicity by 0 
!                                  or target spin multiplicity by 2*K
! For no-exchange files L is replaced by 2*K 
!                   and total spin multiplicity by 0

    ixdro = fa
    stem = 'AMPC'
    if (fine_mesh) stem = 'AMPF'
    qn = val_qno(sqno)
    if (split_prop) then   ! channel splitting in use
       write (sp,'(i2.2)') st  ! use target spin (splitting parameter)
    else
       if (qn(1) < 0) qn(2) = 0 ! for jpi_coupling
       write (sp,'(i2.2)') abs(qn(2)) ! scattering spin (0 for jpi_coupling)
    end if
    write (l,'(i2.2)') abs(qn(1))
    write (p, '(i1.1)') qn(3)
    write (isect, '(i2.2)') sect
    hname = stem // sp // l // p // isect

    if (xdr) then   ! open XDR output file for Hamiltonian
       ixdro = open_xdr (file=TRIM(hname), action='write')
       call xdr_io (ixdro, nh)
       call xdr_io (ixdro, nsect)
    else
       open (unit=ixdro, file=TRIM(hname), access='sequential',  &
            status='replace', form='unformatted', action='write', &
            iostat=ios)
       if (ios /= 0) then
          write (fo,'(a,i6)') 'start_sa_file: error opening ' // &
               TRIM(hname) // ' iostat = ', ios
          stop
       end if
       write (ixdro) nh, nsect
    end if
    if (xdr) then
       write (fo,'(a)') 'XDR Representation Amplitude file: ' // &
            &TRIM(hname)
    else
       write (fo,'(a)') 'Native Representation Amplitude file: ' // &
            &TRIM(hname)
    end if
  end subroutine start_sa_file

  subroutine start_sa_file_id (sqno, st)
! open XDR file to hold sector reduced width amplitudes
    use xdr_files, only: open_xdr, xdr_io
    use energy_grid, only: fine_mesh
    use slp, only: qno, val_qno
    use rmx1_in, only: xdr, split_prop
    use mpi
    type(qno), intent(in)   :: sqno        ! scattering q#s
    integer, intent(in)     :: st         ! second target spin
    character(len=11)        :: hname
    character(len=2)        :: sp, l, s_mpi
    character(len=1)        :: p
    character(len=4)        :: stem
    integer                 :: ios, qn(3), status
    integer                 :: nspn ! scattering spin.  
    integer ( kind = 4 ) ierr
! File naming convention:
!  AMPCSSLLP: coarse mesh, 1st partition (or all channels)
!  AMPFSSLLP: fine mesh, 1st partition (or all channels)
!  AMPCTTLLP: coarse mesh, 2nd partition. TT = target spin multiplicity
!  AMPFTTLLP: fine mesh, 2nd partition, TT = target spin multiplicity
! For no-exchange LS files target spin multiplicity is used
! For J files L is replaced by 2*J and total spin multiplicity by 0 
!                                  or target spin multiplicity by 2*K
! For no-exchange files L is replaced by 2*K 
!                   and total spin multiplicity by 0

    ixdro = fa
    stem = 'AMPC'
    if (fine_mesh) stem = 'AMPF'
    qn = val_qno(sqno)
    if (split_prop) then   ! channel splitting in use
       write (sp,'(i2.2)') st  ! use target spin (splitting parameter)
    else
       if (qn(1) < 0) qn(2) = 0 ! for jpi_coupling
       write (sp,'(i2.2)') abs(qn(2)) ! scattering spin (0 for jpi_coupling)
    end if
    write (l,'(i2.2)') abs(qn(1))
    write (p, '(i1.1)') qn(3)
    write (s_mpi, '(i2.2)') taskid
    hname = stem // sp // l // p // s_mpi


    if (xdr) then   ! open XDR output file for Hamiltonian
       ixdro = open_xdr (file=TRIM(hname), action='write')
       call xdr_io (ixdro, nh)
       call xdr_io (ixdro, nsect)
    else
       open (unit=ixdro, file=TRIM(hname), access='sequential',  &
            status='replace', form='unformatted', action='write', &
            iostat=ios)
       if (ios /= 0) then
          write (fo,'(a,i6)') 'start_sa_file: error opening ' // &
               TRIM(hname) // ' iostat = ', ios
          stop
       end if
       write (ixdro) nh, nsect
    end if
    if (xdr) then
       write (fo,'(a)') 'XDR Representation Amplitude file: ' // &
            &TRIM(hname)
    else
       write (fo,'(a)') 'Native Representation Amplitude file: ' // &
            &TRIM(hname)
    end if
  end subroutine start_sa_file_id


  subroutine close_sa_file ()
   use xdr_files, only: close_xdr
   call flush(ixdro)
   call close_xdr(ixdro)
   
  end subroutine close_sa_file


  subroutine ampltd (ra1, ra2, nc, nl, nh)
! ampltd : computes surface amplitudes and writes to file
    use xdr_files, only: xdr_io
    use sec_ham, only: vbl, vbr
    use dist_mat, only: evals
    use rmx1_in, only: bug8, xdr
!    use blacs, only: p0, q0, ctxt, io_processor, p_error, dlen_, &
!         ctxt_, mb_
    use error_out, only: error_check
    use dist_mat, only: z
    integer, intent(in)      :: nc         ! # channels
    integer, intent(in)      :: nl         ! # Legendre functions
    integer, intent(in)      :: nh         ! sec Hamiltonian dimension
    real(wp), intent(in)     :: ra1, ra2   ! end pt radii of subrange
    real(wp), allocatable    :: ampal(:), ampar(:) ! local amps
    real(wp), allocatable    :: evec(:)
    integer                  :: k, j, ik, i, jj, status
    integer                  :: mm, nn, lwork, mb, nb, rsrc, csrc
    integer                  :: ldd, isize, istart, iend
    integer                  :: jsize, jstart, jend
    real(wp)                 :: alpha, beta, sml, smr
!    integer                  :: descx(dlen_)

!    if (io_processor) then
       allocate (ampal(nc), ampar(nc), stat=status)
       if (status /= 0) then
          write (fo,'(a,i6)') 'ampltd: allocation error = ', status
          STOP
       end if
    
       if (bug8 > 0) then
          write (fo,'(a,e20.12)') 'Sector eigenvalue sum = ', SUM(evals)
       else if (bug8 > 1) then
          write (fo,'(a)') 'R-Matrix Eigenvalues (Ryd)'
          write (fo,'(5f12.5)') evals
       end if

! write eigenvalues to disk
       if (xdr) then
          call xdr_io (ixdro, evals, nh)  ! write R-matrix eigenvalues
       else
          write (ixdro) evals
       end if
!    end if

    allocate (evec(nh), stat=status)
    call error_check (status, 'ampltd: evec allocation')
    lwork = 96
!    lwork = descz(mb_)   ! eigenvectgor block size
!    ctxt = descz(ctxt_)  ! context

! form descriptor for evec
    mm = MAX(1, MIN(nh, lwork))
    nn = MAX(1, INT(lwork / mm))
    mb = mm
    nb = nn
!    rsrc = p0
!    csrc = q0
!    ldd = MAX(1, mm)
!    call descset (descx, mm, nn, mb, nb, rsrc, csrc, ctxt, ldd)

    alpha = 1.0_wp
    beta = 0.0_wp
    sml = 0.0_wp
    smr = 0.0_wp
    evec = z
    jst: do jstart = 1, nh, nn
       jend = MIN(nh, jstart+nn-1)
       jsize = jend - jstart + 1
 !      row_blocks: do istart = 1, nh, mm
 !         iend = MIN(nh, istart+mm-1)
 !         isize = iend - istart + 1
 !         call pdgeadd ('notrans', isize, jsize, alpha, z, istart, &
 !              jstart, descz, beta, evec(istart:iend), 1, 1, descx)
 !      end do row_blocks
 !      call BLACS_BARRIER (ctxt, 'A')
! now have the full eigenvector on the i/o processor:

!       if (io_processor) then
          ampal = 0.0_wp
          ampar = 0.0_wp
          ik = 0
          do i = 1, nc
             do j = 1, nl
                ik = ik + 1
                ampal(i) = ampal(i) + evec(ik) * vbl(j)
                ampar(i) = ampar(i) + evec(ik) * vbr(j)
             end do
          end do
          
          if (bug8 > 0) then
             sml = sml + SUM(ampal)
             smr = smr + SUM(ampar)
          else if (bug8 > 2) then
             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra1,&
                  ' Column: ', jstart
             write (fo,'(5e14.6)') ampal
             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra2,&
                  ' Column: ', jstart
             write (fo,'(5e14.6)') ampar
          end if

! write sector reduced width amplitudes to disk
          if (xdr) then
             call xdr_io (ixdro, ampal, nc) ! write left amplitudes to disk
             call xdr_io (ixdro, ampar, nc) ! write right amplitudes to disk
          else
             write (ixdro) ampal(1:nc)
             write (ixdro) ampar(1:nc)
          end if
 !      end if
    end do jst
!    if (io_processor) then
       if (bug8 > 0) write (fo,'(a,2e20.12)') 'RWAmplitude sums = ', &
            sml, smr
       deallocate (ampal, ampar, evec, stat=status)
!    else
!       deallocate (evec, stat=status)
!    end if
    call error_check (status, 'ampltd: deallocation')
  end subroutine ampltd
end module rw_amplitudes
