
c$debug
c***********************************************************************
      SUBROUTINE SOLVE
c
c     Subroutine to solve for the S matrix elements from the matching equations
c     at the end points that are outputted from INTEG3.
c     
c     Amplitudes are stored in array D
c***********************************************************************
c
c     KMAX  = Number of integration points
c     NCHAN = Number of sub channels 
c     NCH1  = Number of sub channels in channel 1
c     ispc  = Labeled common sizes
c     nlev  = Number of levels
c
      parameter (kmax = 300, nchan = 64, ispc0 = 16384, ispc1 = 3*16384
     1, ispc4 = 8192, nlkmax = ispc1/2/kmax, nlev = 12)
      parameter (nch1 = 4096/4/nchan, nmatrix = nchan*nchan)
c
      complex*16 amat(nchan,nchan),cmat(nchan,nchan), bmat(nchan,nch1)
     1,          zi(nchan,2), temp
      complex* 8 E(16)
      real*8 dmat(2*nchan,nch1), s(800)
      logical igam(nlev), isym1
c
      COMMON ALPHA(15),HBARC,AMU,CHGSQ,DRT,RMAX
     1,ICON(20),LMAX,KT,MAXD,IOUT,JOUT,NLTR,IDAT(6)
     2,KLAM(nlev),LTRT(nlev),JTRT(nlev),NADD(nlev)
      COMMON EN(nlev),FM(nlev),Z(nlev),FMA(nlev),ZA(nlev),RY(nlev)
     1,FS(nlev),ECM(nlev),FMU(nlev),RC(nlev),FK2(nlev),FK(nlev)
     2,ETA(nlev),AFACT(nlev),PFACT(nlev),DRF(nlev),FN(nlev)
     3,PNLOC(nlev),DUM(nlev),TSIG(nlev)
     4,IS(nlev),NS(nlev),NLI(nlev),IN(nlev),INN(nlev,nlev)
     5,LAM(nlkmax),ISP(nlkmax),JSP(nlkmax),ILK(nlkmax)
     6,LKTAB(2*nlkmax),ICUP(2*nlkmax),JCUP(2*nlkmax)
c
      COMMON /label0/Space0(ispc0)
      COMMON /label1/Space1(ispc1)
      COMMON /label4/Space4(ispc4)
      COMMON /LLL   /NLVAL,NBLKS,NLMIN,NLDEL,MAXTAB(nlev)
     1,INDTAB(nlev+2)
      common/ch3files/input, ioutp
c
      DIMENSION NTAB(nchan),LTAB(nchan),JTAB(nchan),FTAB(nchan)
     1,         D(ispc0),ctab(nchan),stab(nchan),Sig(nchan)
c
      EQUIVALENCE (Space0(1),d(1)) 
     1, (Space1(        1),amat(1,1)), (Space1(4*nmatrix+1),cmat(1,1))
     2, (Space4(        1),bmat(1,1)), (Space4(   1),dmat(1,1))
     3, (Space4(4101+0*nchan),ntab(1)),(Space4(4101+1*nchan),ltab(1))
     4, (Space4(4101+2*nchan),jtab(1)),(Space4(4101+3*nchan),ftab(1))
     5, (Space4(4101+4*nchan),ctab(1)),(Space4(4101+5*nchan),stab(1))
     6, (Space4(4101+6*nchan),zi(1,1))
     7, (Space4(6401)        ,s(1))
c
	data FourPi/12.5663706/
C
      ICOUNT=0
      JMAX=LMAX
      LK=0
      IX=JTRT(1)
      MAX1=JTRT(1)+1
      MXS=MAX1*NS(1)
      JX1=(1.0-PHASEF(JTRT(1)+IS(1)))/2.0
      NPAR=2
      IF(MXS.EQ.1) NPAR=1
c
c     Setup
c
      INDTAB(1)=0
      NLDEL=0
      IND=1
      DO 1100 I=1,NLTR
      TSIG(I)=0.0
      IND=IND+LMAX
      LKTEMP=(JTRT(I)+1)*NS(I)
      LK=LK+LKTEMP
      MAXTAB(I)=LKTEMP
      INDTAB(I+1)=LK
      NLDEL=MAX0(JTRT(I)+IS(I)+JX1,NLDEL)
      if(fm(i) .eq. 0.0 .and. is(i) .eq. 2) then
        igam(i) = .true.
      else
        igam(i) = .false.
      endif
 1100 CONTINUE
      if(fm(1).eq.fma(1) .and. z(1).eq.za(1)
     1                   .and. is(1).eq.jtrt(1)) then
        isym1 =.true.
      else
        isym1 = .false.
      endif
c
      IS0 = 0
      IS1 = IS(1)
C
C     DETERMINE BLOCK SIZE FOR AMPLITUDES
C
      NLVAL=ispc0/(2*LK*MXS)
      NLMIN=LMAX
      NLVAL=MIN0(NLVAL,LMAX)
      IF(NLVAL.LT.LMAX) NLMIN=NLVAL-NLDEL
      DO 1106 I=1,NLTR
      INDTAB(I+1)=INDTAB(I+1)*(2*MXS*NLVAL)
 1106 CONTINUE
      IF(NLMIN.le.0) then
        IOUT=1
        WRITE(ioutp,'(''0NLMIN.LE.0 ***  D-TABLE TOO SMALL '')')
      endif
 1108 CONTINUE
      NBLKS=(LMAX+NLMIN-1)/NLMIN
      if(icon(19).eq.0) then
      WRITE(ioutp,9011)NBLKS,LMAX,NLMIN,NLVAL
      endif
C
C     CLEAR AMPLITUDE STORAGE
C
      MAX=INDTAB(NLTR+1)
      DO 1110 I=1,MAX
 1110 D(I)=0.0
C
C     J LOOP
C
      if(icon(19).eq.0) then
      if(icon(5).ne.0.or.icon(7).ne.0) then
        WRITE(ioutp,9024)ALPHA,(IDAT(I),I=1,3)
        write(ioutp,9020)
      endif
      endif
c      Correct multiplicity for transverse gamma ray
      if(igam(1)) then
        Sfact = ns(1)-1
      else
        Sfact = ns(1)
      endif
c
      JA1=1
      JA2=LMAX
      IF(NLMIN.LT.LMAX) JA2=NLVAL-NLDEL/2
      JAX=0
      DO 2500 NB=1,NBLKS
      DO 2000 JA=JA1,JA2
        if(icon(19).eq.0) 
     1     write(*,'(''+Solve,  JA = '',i3,''/'',i3)')ja,lmax
      DO 1900 JPAR=1,NPAR
c
      read (3) nll,nli,(ntab(i),ltab(i),jtab(i),i=1,nll)
c
      ncol=NLI(1)
      k=NCOL
      IF(NCOL.EQ.0) GO TO 1900
      MAXA=NLL+NLL
c
c     set up normalization tables
c
      AMFACT=2.0*AMU/HBARC**2
      DO 1120 N1=1,NLL
      I=NTAB(N1)
      L=LTAB(N1)+1
      IND=(I-1)*LMAX
      ARG=S(IND+L)-S(IND+1)
      ctab(n1)=cos(arg)
      stab(n1)=sin(arg)
      if(igam(i)) then
        ftab(n1) = sqrt(((drf(i)/drf(1))**3)*float(2*ltab(n1)+1)
     1   * (fk(1)/(fk(i)*amfact*fmu(1))))/fk(i)
      else
        ftab(n1) = sqrt(((drf(i)/drf(1))**3)*float(2*ltab(n1)+1)
     1   * (fmu(i)*fk(1)/(fk(i)*fmu(1))))/fk(i)
      endif
 1120 CONTINUE
C
C     SET UP AND INVERT MATRIX FOR SCATTERING AMPLITUDES
C     EQUATIONS ARE  {A1}*[C] = [F1] + {H1}*[S]
C                    {A2}*[C] = [F2] + {H2}*[S]
C
C                    [C] =  NORM VECTOR
C                    [S] =  S-MATRIX VECTOR
C                    {A1}= N END POINTS AT RMAX-DR
C                    {A2}= N END POINTS AT RMAX
C                    [F1]= REGULAR FCT AT  RMAX-DR
C                    [F2]= REGULAR FCT AT  RMAX
C                    {H1}= HANKEL  FCT AT  RMAX-DR
C                    {H2}= HANKEL  FCT AT  RMAX
C
C     COEFFICIENT MATRIX IS SOLVED BY PARTITIONS
C     (A2*A1(INV)*H1-H2)*SMATRIX=-(A2*A1(INV)*F1-F2)
C
C     FORM A1(TR,INV)*A2(TR) ---> AMAT
C
C     A1(TR) ==> CMAT
C     A2(TR) ==> AMAT
c
      DO 1140 N2=1,NLL
      read (3)(cmat(n2,i),amat(n2,i), i=1,nll)
 1140 CONTINUE
C
C     AMAT = A1(TR,INV)*A2(TR)
C
      CALL INVERT(NLL,NLL,CMAT,AMAT)
C
C     TRANSPOSE AMAT INTO CMAT TO GIVE
C     A2*A1(INV) ---> CMAT
C
      DO 1150 N1=1,NLL
      DO 1150 N2=1,NLL
      CMAT(N2,N1)=AMAT(N1,N2)
 1150 CONTINUE
c	if(ja .eq. 1) then
c	write(20, '(i4,1p8e15.7)') ja,((cmat(n1,n2),n1=1,nll),n2=1,nll)
c	endif
C
C     AMAT =   A2*A1(INV)*H1-H2
C
      do 1170 n2=1,nll
      read (3)(zi(i,1),zi(i,2),i=1,nll)
      do 1165 n1=1,nll
      temp=0.0
      do 1160 i=1,nll
      temp=temp+cmat(n1,i)*zi(i,1)
 1160 continue
      amat(n1,n2)=(temp-zi(n1,2))*ftab(n2)
 1165 continue
 1170 continue
c
C     BMAT = -(A2*A1(INV)*F1-F2)
c
      do 1190 n2=1,ncol
      read (3)(zi(i,1),zi(i,2),i=1,nll)
      do 1185 n1=1,nll
      temp=0.0
      do 1180 i=1,nll
      temp=temp+cmat(n1,i)*zi(i,1)
 1180 continue
      bmat(n1,n2)=-(temp-zi(n1,2))*ftab(n2)
 1185 continue
 1190 continue
C
C
      CALL INVERT(NLL,NCOL,AMAT,BMAT)
C
C     BMAT = SMATRIX
C
      JT=JA+JA-2+JX1
c
      DO 1600 N1=1,NCOL
      NN=2*N1-1
      JX=JTAB(N1)
      LX=2*LTAB(N1)
      FL=LX + 1
      if(isym1 .and. ntab(N1).eq.1 .and. is(1).eq.0) then
        if(phasef(ltab(n1)) .lt. 0.0) then
          bmat(n1,n1) = 0.0
        endif
      endif
      E(N1)=bmat(N1,N1)
c
      DO 1500 N2=1,NLL
      NK=2*N2-1
      JJ=NTAB(N2)
      JY=JTAB(N2)
      LZ=LTAB(N2)
      LY=2*LTAB(N2)
      IY=JTRT(JJ)
      ISJ = IS(JJ)
      MAXI=JTRT(JJ)+1
C
C     Reaction and total inelastic cross section
C
      FACT=float(jt+1)*float(ly+1)/float(lx+1)
      T1 = conjg(bmat(n2,n1))* bmat(n2,n1)
      IF(JJ.EQ.1) THEN
        if(n1.eq.n2) then
        T1 = DIMAG(bmat(n2,n1)) - T1
        else
        T1 =-T1
        endif
c        write(20,'(3(a,i4),1p8e12.4)')' ja =',ja,'  n1 =',n1,'  n2 =',n2
c     1	,bmat(n2,n1),fact,t1,fact*t1,tsig(1)
      if(isym1 .and. ntab(N1).eq.1 .and. is(1).eq.0) then
        bmat(n2,n1) = 2.0*bmat(n2,n1)
        T1 = 4.0*T1
      endif
      ENDIF
      TSIG(JJ)=TSIG(JJ) + T1*FACT
c
      T2 = FourPI/(FLOAT(JTRT(1)+1)*Sfact*FK(JJ)**2)
      Sig(n2)=T1*fact*T2
c      write(20,'(12x,1p8e12.4)')bmat(n2,n1),fact,T2,T1,Sig(n2)
c
c     Multiply by coulomb phases
c
      bmat(n2,n1)=bmat(n2,n1)*dcmplx(ctab(n2),stab(n2))
     1                       *dcmplx(ctab(n1),stab(n1))/fk(jj)
    
C
C     STORE S-MATRIX ELEMENTS
C
      JR=NS(1)
      JS=NS(JJ)
      PHASE=PHASEF((LX+IS1+IX-LY-IS(JJ)-IY)/2)
      INDEX=2*(MAXTAB(1)*MAXTAB(JJ)*(LZ-JAX))+INDTAB(JJ)+1
      IND=INDEX+MAX1*MAXI*JR*JS*2-2
C     error trap
      IF(INDEX.le.0.and.icount.le.50) then
      ICOUNT=ICOUNT+1
      if(icon(19).eq.0) then
      WRITE(ioutp,7790)JA,LZ,JA1,INDEX
 7790 FORMAT(' JA,LZ,JA1=',3I5,'  INDEX=',I5)
      endif
C
      endif
      MM=-JTRT(1)
      IF(MM.EQ.0) GO TO 1310
      M1=(MAX1+1)/2
      MI=MAXI
      GO TO 1311
 1310 CONTINUE
      M1=1
      MI=(MAXI+1)/2
 1311 CONTINUE
      DO 1490 MX=1,M1
      MMS=-JTRT(JJ)
      DO 1480 MY=1,MI
      MS=-IS(1)
      DO 1470 M =1,JR
      MSP=-IS(JJ)
      DO 1460 MP=1,JS
      ML=MS+MM-MSP-MMS
      MLT=IABS(ML/2)
      IF(MLT.GT.LZ) GO TO 1458
      PHAS=1.0
      IF(ML.LT.0) PHAS=PHASEF(MLT)
C
         COEFF1 = VCC(LX,IS1,JX,0,MS)*VCC(JX,IX,JT,MS,MM)
     1   *        VCC(LY,IS(JJ),JY,ML,MSP)*VCC(JY,IY,JT,ML+MSP,MMS)
     2   *SQRT(YXFCT(LZ+MLT,LZ-MLT))*PHAS*FLOAT(LY+1)
C
C
      D(INDEX  )=D(INDEX  )+dmat(nk  ,n1)*COEFF1
      D(INDEX+1)=D(INDEX+1)+dmat(nk+1,n1)*COEFF1
C
C     FILL OUT D TABLE
C
      IF(JTRT(1)+JTRT(JJ).EQ.0)   GO TO 1458
      IF(M1.NE.MAX1.AND.MM .EQ.0) GO TO 1458
      IF(MI.NE.MAXI.AND.MMS.EQ.0) GO TO 1458
      PHAS=PHASE*PHAS
      IF(ML.GT.0) PHAS=PHASE*PHASEF(MLT)
      D(IND  )=D(INDEX  )*PHAS
      D(IND+1)=D(INDEX+1)*PHAS
 1458 CONTINUE
      INDEX=INDEX+2
      IND=IND-2
      MSP=MSP+2
 1460 CONTINUE
      MS=MS+2
 1470 CONTINUE
      MMS=MMS+2
 1480 CONTINUE
      MM=MM+2
 1490 CONTINUE
 1500 CONTINUE
c     Write out radial matrix elements
      if(icon(19).eq.0) then
      IF(ICON(7).ne.0) then
          K1=1
          DO 1590 I=1,NLTR
          if(nli(i).ne.0) then
            K2=K1+NLI(I)-1
	    write(ioutp,9022)(ltab(m),jtab(m),m=k1,k2)
            WRITE(ioutp,9021)I,JT,ltab(n1),jtab(n1),(BMAT(M,N1),M=K1,K2)
            if(icon(7) .eq. 2) then
	      write(ioutp,'(12x,1p10e12.4)')(Sig(m),m=k1,k2)
            endif
            K1=K2+1
          endif
 1590   CONTINUE
      endif
      endif
 1600 CONTINUE
c     write out elastic scattering amplitudes
      if(icon(19).eq.0) then
      IF(ICON(5).ne.0) then
      I=1
      K2=NCOL
      WRITE(ioutp,9023)I,JT,ltab( 1),jtab( 1),(E(M),M=1,K2)
c     write out exp(2i*delta)
      IF(ICON(5).eq.2) then
      DO 1610 M=1,NCOL
      E(M)=1.0 + cmplx(0.0, 2.0)*E(M)
 1610 CONTINUE
      WRITE(ioutp,9026)(E(M),M=1,K2)
      endif
      endif
      endif
 1900 CONTINUE
 2000 CONTINUE
      IF(NBLKS.ne.1) then
      DO 2400 I=1,NLTR
      I1=INDTAB(I)+1
      I2=I1-1+MAXTAB(I)*(NLMIN*MXS*2)
      I3=I1-1+MAXTAB(I)*(NLVAL*MXS*2)
      IF(NB.EQ.NBLKS) I2=I3
      JA3=     JAX+1
      JA4=MIN0(JAX+NLMIN,LMAX)
      WRITE(2)I1,I2,JA3,JA4,(D(IT),IT=I1,I2)
C
C     MOVE AMPLITUDES DOWN IN MEMORY
C
      DO 2370 J=I1,I2
      D(J)=0.0
 2370 CONTINUE
      M=I3-I2
      IF(M.EQ.0) GO TO 2381
      DO 2380 J=1,M
      D(I1-1+J)=D(I2+J)
      D(I2+J)=0.0
 2380 CONTINUE
 2381 CONTINUE
 2400 CONTINUE
      JA1=     JA2+1
      JA2=MIN0(JA2+NLMIN,LMAX)
      JAX=     JAX+NLMIN
      endif
 2500 CONTINUE
      DO 2100 I=1,NLTR
      TSIG(I)=TSIG(I)*FourPi/(FLOAT(JTRT(1)+1)*Sfact*FK(I)**2)
 2100 CONTINUE
      REWIND 2
      REWIND 3
C
      RETURN
 9011 FORMAT('0NBLKS =',I2,3X,'LMAX =',I3,3X,'NLMIN =',I3
     1      ,3X,'NLVAL =',I3)
 9020 format('0ch 2J L1 j1','  (exp(2i*delta) - 1)/(2i)')
 9021 FORMAT(' ',I2,3I3,1P10E12.4,: (/12X,1P10E12.4))
 9022 FORMAT((11x, 5(7h   l2 =,i2,8h  2*j2 =,i2,5x:)))
 9023 FORMAT(' ',I2,3I3,  10F12.6,: (/12X,  10F12.6))
 9024 FORMAT('1SCATTERING AMPLITUDES' 
     1 ,17X,15A4,I4,2(1H/,I2),I4,2(1H.,I2))
 9026 format(12X,10F12.6)
      END


c***********************************************************************
      SUBROUTINE INVERT(MAX,N,A,C)
c
c     Inversion routine used to solve end point matching equations
c     
c***********************************************************************
c
C     KMAX  = Number of integration points
C     NCHAN = Number of sub channels in INTEG2
      parameter (kmax = 300, nchan=64)
      COMPLEX*16 A(nchan,nchan),C(nchan,nchan),B
c
      DO 100 I=1,MAX
      B=1.0/A(I,I)
      DO 10 K=1,N
      C(I,K)=C(I,K)*B
   10 CONTINUE
      DO 20 J=I,MAX
   20 A(I,J)=A(I,J)*B
      DO 50 J=1,MAX
      IF(I.EQ.J) GO TO 50
      B=A(J,I)
      DO 35 K=1,N
      C(J,K)=C(J,K)-C(I,K)*B
   35 CONTINUE
      DO 40 K=I,MAX
   40 A(J,K)=A(J,K)-A(I,K)*B
   50 CONTINUE
  100 CONTINUE
      RETURN
      END


