SUBROUTINE propph1(nt, R_Mat, w, vecnew, tstore, rbegin,          &
                    rmax, nsteps, RE_Mat, zp, usys2, w2,          &
                    zedz, eredev, EDeriv, iprt, zeta,             &
                    omega, new, g0, g1, g2, w0, w1, w2b,          &
                    psi, dmu, energy, gamma,                      &
                    xksq, ncall,smx_time,mcall, sge_time,         &
                    method, diag, cor, asym, ndim, ipiv, ui, work)
USE FileUnits_MODULE
USE dip_MODULE
USE junk_MODULE
USE Boundary_Module
USE BasisRegions_Module
!
!

!             P U R P O S E    O F    S U B R O U T I N E

! Integrates second-order coupled differential equations across a
! region.
! written by G. A. Parker
!
!   Modified by M. Braunstein to treat photodissociation - Spring, 1991
!
!   Stefano's modifications of 9/11/97 added 9/12/97:
!      The variables ipiv, ui, work are now passed in as arguments,
!      which are now dimensioned by the passed-in argument nt.
!    9/16/97  TRP  Added ipiv, ui, work to argument list in CALL
!                     to Manolopous().

!             I N P U T   A R G U M E N T S

!   nt        Number of coupled second-order differential equations.
!   R_Mat      LOG derivative matrix at the begining of the region.
!   w         Temporary storage array.
!   vecnew    Temporary storage array.
!   tstore    Temporary storage array.
!   nsteps    The number of steps in the log-derivative region.
!   rbegin    Beginning distance for the log-derivative region.
!   rmax      Ending distance for the log-derivative region.
!
!   zeta      Gamma vector (inhomogeneous operator) at beginning of
!             region
!

!             O U T P U T     A R G U M E N T S

!   R_Mat    Wigner R-matrix at the end of the region.
!
!   zeta    Gamma vector (inhomogeneous operator) at end of region
!

! <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>

IMPLICIT NONE
CHARACTER(LEN=6) iprt
CHARACTER(LEN=10) method
LOGICAL EDeriv, ifirst
INTEGER nt, isect, nrho, nsteps, isector, i, ncall, mcall, ndim, info
INTEGER ipiv(nt)
REAL(Kind=WP_Kind) rbegin, rmax, rright, rleft, rhov1, xktot, rmax_save
REAL(Kind=WP_Kind) rhov2, usys2, eredev, rhov1g, rhov2g, smx_time
REAL(Kind=WP_Kind) t_smx1, t_smx2, sge_time, t_sge1, t_sge2
REAL(Kind=WP_Kind) work(nt)
!-----------------------------------------------------------------------
! arrays dimensioned as matrices (nt*nt).
!-----------------------------------------------------------------------
REAL(Kind=WP_Kind) w(nt,nt), R_Mat(nt,nt), vecnew(nt), tstore(nt,nt)
REAL(Kind=WP_Kind) RE_Mat(nt,nt), zp(nt,nt), w2(nt,nt), zedz(nt,nt)
REAL(Kind=WP_Kind) g0(nt,nt), g1(nt,nt), g2(nt,nt), dmu(nt,nt), psi(nt)
REAL(Kind=WP_Kind) gamma(nt), rhoval(21), xksq(nt), energy(nt)
REAL(Kind=WP_Kind) ui(nt,nt), w0(nt,nt), w1(nt,nt), w2b(nt,nt)
REAL(Kind=WP_Kind) zeta(nt), omega(nt), new(nt), diag(nt,nt), cor(nt,nt)
REAL(Kind=WP_Kind) asym(nt,nt)
SAVE
COMMON /utotal/ xktot
DATA ifirst/.true./

EXTERNAL popt, mxoutd, trnsfm, sector, overlap
!iprt='little'

t_smx1 = 0.d0
t_smx2 = 0.d0
!-----------------------------------------------------------------------
! Write out intgration parameters.
!-----------------------------------------------------------------------
WRITE(Out_Unit,'(A,I5,A,ES23.15,A,ES23.15)') 'nsteps=',nsteps, ' rbegin=', rbegin, ' rmax=',rmax

!-----------------------------------------------------------------------
! calculate lower and upper bounds for energy-labeled
! channels which are to be propagated
!-----------------------------------------------------------------------
WRITE(Out_Unit,*) 'nt=',nt, '   xktot=', xktot
rmax_save=rmax
rleft = rbegin
rright = rmax
rhov1 = rbegin
isect = 1
READ(Ovr_Unit)ISect,Rhov1,Rhov2
nrho = NSectors
Rhov1=Basis_Dist(Isect)
Rhov2=Basis_Dist(Isect+1)
Rleft=Sector_Boundaries(Isect)
Rright=Sector_Boundaries(Isect+1)
Rmax=Sector_Boundaries(Nsectors+1)
Nrho=Nsectors
IF(lphoto)THEN
   READ(fph_unit)rhov1g,rhov2g
ENDIF

DO isector = 1, NSectors
   IF(isector/=1)newsect=.true.
   IF(isector/=1)THEN
      READ(Ovr_Unit)ISect,Rhov1,rhov2
      Rhov1=Basis_Dist(Isector)
      Rhov2=Basis_Dist(Isector+1)
      Rleft=Sector_Boundaries(Isector)
      Rright=Sector_Boundaries(Isector+1)
      Rmax=Sector_Boundaries(Nsectors+1)
      IF(isector==NSectors)THEN
         RRight=rmax_save
         Rmax=rmax_save
      ENDIF
      Nrho=Nsectors
      IF(lphoto)THEN
         READ(fph_unit)rhov1g,rhov2g
      ENDIF
      nrho = NSectors
   ENDIF
   !-----------------------------------------------------------------------
   ! integrate from rleft to rright using the log-derivative integrator.
   !-----------------------------------------------------------------------
   IF(ifirst.AND.(iprt/='little'))THEN
       WRITE(Out_Unit,'(A,es23.15,A,es23.15,A,I5)')'rleft=',rleft,' rright=',rright,' nsteps=',nsteps
   ENDIF
   IF(method=='new_logder')THEN
      !----------------------------------------------------------------------
      ! This method does not use matrix inversions but uses matrix solutions.
      !----------------------------------------------------------------------
      IF(isector==1)THEN
         WRITE(Out_Unit,'(A)')'Using new_logder'
      ENDIF
      CALL logderph1(tstore, R_Mat, w, nt, rleft, rright, nsteps, RE_Mat,  &
                     zp, usys2, w2, zedz, EDeriv, zeta, omega, new,        &
                     g0, g1, g2, dmu, psi, gamma,                          &
                     rhoval, xksq, energy, w0, w1, w2b,                    &
                     diag, cor, asym, ndim, ipiv, ui, work)
   ELSEIF(method=='man_method')THEN
      !----------------------------------------------------------------------
      ! This is the Manolopous Procedure
      !----------------------------------------------------------------------
      IF(isector==1)THEN
         WRITE(Out_Unit,*)'Using Manolopous Method'
         WRITE(Msg_Unit,*)'Using Manolopous Method'
      ENDIF
      CALL Manolopous(tstore, R_Mat, w, nt, rleft, rright, nsteps, RE_Mat,  &
                     zp, usys2, w2, zedz, EDeriv, zeta, omega, new,         &
                     g0, g1, g2, dmu, psi, gamma,                           &
                     rhoval, xksq, energy, w0, w1, w2b,                     &
                     diag, cor, asym, ndim, ipiv, ui, work)
   ELSE
      WRITE(Msg_Unit,*)'Error variable method not set correctly'
      WRITE(Msg_Unit,*)method
   ENDIF
   IF(iprt=='medium')THEN
      WRITE(Out_Unit,*) 'logder_matrix at r=', rright
      CALL MxOut(R_Mat, nt, nt)
      IF(EDeriv)THEN
         WRITE(Out_Unit,*) 'Derivative of logder_matrix', ' with respect to E at r=', rright
         CALL MxOut(RE_Mat, nt, nt)
      ENDIF
   ELSEIF(iprt=='full  ')THEN
      WRITE(Out_Unit,*) 'logder_matrix at r=', rright
      CALL MxOutD(R_Mat, nt, nt, 0)
      IF(EDeriv)THEN
         WRITE(Out_Unit,*) 'Derivative of logder_matrix', ' with respect to E at r=', rright
         CALL MxOutD(RE_Mat, nt, nt, 0)
      ENDIF
   ENDIF
   !-----------------------------------------------------------------------
   ! If sector adiabatic basis functions are used we need to obtain
   ! the overlap matrix and transform to the new basis.
   !-----------------------------------------------------------------------
   IF(iprt=='medium')THEN
      WRITE(Out_Unit,'(1x,A,ES15.7,A,ES15.7,A,ES15.7)')'RLeft=',RLeft,'  RRight=',RRight, ' Ovlp(1,1)=',TStore(1,1)
   ENDIF
   IF(isector /= nrho)THEN
      tstore=0.d0
      CALL overlap2(rhov1, rhov2, tstore, vecnew, nt, ifirst)
      IF(iprt=='medium')THEN
         WRITE(Out_Unit,*)'Overlap Matrix'
         CALL MxOut(tstore, nt, nt)
      ENDIF
      w=MATMUL(R_Mat,TStore)
      R_Mat=MATMUL(TRANSPOSE(TStore),w)
      IF(EDeriv)THEN
         w=MATMUL(RE_Mat,TStore)
         RE_Mat=MATMUL(TRANSPOSE(TStore),w)
      ENDIF
      !------------------------------------------------------------------
      ! transfrom gamma - half transfrom of vector
      !------------------------------------------------------------------
      IF(lphoto)THEN
         IF(iprt=='medium')THEN
            WRITE(Out_Unit,*)'gamma vector'
            WRITE(Out_Unit,*)(zeta(i),i=1,nt)          ! WRITE intermediate gamma to file
         ENDIF  
         CALL dgemv('T',nt,nt,1.d0,tstore,nt,zeta,1,0.d0,w,1)
         DO i=1,nt
            zeta(i)=w(i,1)
         ENDDO
      ENDIF
      IF(iprt=='medium')THEN
         WRITE(Out_Unit,*) 'logder_matrix at r=', rright
         WRITE(Out_Unit,*)'transformed with an overlap',  ' matrix calculated',   &
                          ' between functions at rho = ',  rhov1,' and ',rhov2
         CALL MxOut(R_Mat, nt, nt)
         IF(EDeriv)THEN
            WRITE(Out_Unit,*) 'Derivative of logder_matrix',  ' with respect to E at r=', rright
            WRITE(Out_Unit,*)'transformed with an overlap',  ' matrix calculated',  &
                             'between functions at rho = ',  rhov1,' and ',rhov2
            CALL MxOut(RE_Mat, nt, nt)
         ENDIF
      ELSEIF(iprt=='full  ')THEN
         WRITE(Out_Unit,*) 'logder_matrix at r=', rright
         WRITE(Out_Unit,*)'transformed with an overlap',  ' matrix calculated',  &
                          'between functions at rho = ',  rhov1,' and ',rhov2
         CALL MxOutD(R_Mat, nt, nt, 0)
         IF(EDeriv)THEN
            WRITE(Out_Unit,*) 'Derivative of logder_matrix',  ' with respect to E at r=', rright
            WRITE(Out_Unit,*)'transformed with an overlap',  ' matrix calculated', &
                            'between functions at rho = ',  rhov1,' and ',rhov2
            CALL MxOutD(RE_Mat, nt, nt, 0)
         ENDIF
      ENDIF
   ENDIF
ENDDO

!-----------------------------------------------------------------------
! calculate R = Y**(-1)
!-----------------------------------------------------------------------
CALL cputime(t_smx1)
CALL unit (nt, ui)
CALL dsysv ('l', nt, nt, R_Mat, nt, ipiv, ui, nt, work, nt, info)
R_Mat=ui

CALL cputime(t_smx2)
ncall = ncall + 1
smx_time = smx_time + t_smx2-t_smx1
!-----------------------------------------------------------------------
! calculate RE = - R * YE * R
!-----------------------------------------------------------------------
IF(EDeriv)THEN
   RE_Mat=MATMUL(R_Mat,RE_Mat)
   RE_Mat=-MATMUL(RE_Mat,R_Mat)
ENDIF

rmax = rright
ifirst = .false.
RETURN
END