c$debug
c***********************************************************************
      SUBROUTINE INTEG3
c
c     Suboutine for integration of coupled equations.
c     
c***********************************************************************
C
C     KMAX  = Number of integration points
C     NCHAN = Number of sub channels in INTEG3
c     ispc  = Labeled common sizes
c     nlev  = Number of levels
c
      parameter (kmax = 300, nchan = 64, ispc0 = 16384, ispc1 = 3*16384
     1, ispc4 = 12900, nlkmax = ispc1/2/kmax, nlev = 12)
      parameter (maxij=500, maxsr=1000) 
c
      complex*16 zi(nchan,3),cz(nchan,3),c0(nchan)
      complex*8  W(kmax*nlev), V(kmax*nlev) 
     1          ,XC1(ispc1/2), X(maxij)
      real*8 F(800),FP(800),G(800),GP(800),S(800)
     1, delr, x1, x2, eta_dp
      logical igam(nlev), isym(nlev), iclosed, jpflag, isymp
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 /integ /icode, ii,jj, klii,kljj, fmii,fmjj
     1              ,lx,isx,jx,ix, ly,isy,jy,iy, ll2,is2,jj2, jt
	common/ch3files/input, ioutp
c
      DIMENSION NTAB(nchan),LTAB(nchan),JTAB(nchan),MTAB(nchan)
     1,         NTBL(nchan),FTAB(nchan),VTAB(nchan),mz  (nchan)
     2,         IKTAB(500),JKTAB(500),CTAB(500),LKTLK(500),ICTBL(500)
     3,         DR2(nlev),ITEMP(nlev,nlev)
c
      EQUIVALENCE (Space0(   1),W(  1)),(Space0(nlev*2*kmax+1),V(  1))
     1,           (Space1(   1),XC1(1))
     2,           (Space4(   1),f (1)) ,(Space4(1601),fp(1))
     3,           (Space4(3201),g (1)) ,(Space4(4801),gp(1))
     4,           (Space4(6401),s(1))  
     5,       (Space4( 8001),iktab(1)) ,(Space4( 8501),jktab(1))
     6,       (Space4( 9001),ctab (1)) ,(Space4( 9501),lktlk(1))
     7,       (Space4(10001),ictbl(1))
     8,(Space4(10501+ 0*nchan),ntab(1)),(Space4(10501+ 1*nchan),ltab(1))
     9,(Space4(10501+ 2*nchan),jtab(1)),(Space4(10501+ 3*nchan),mtab(1))
     $,(Space4(10501+ 4*nchan),ntbl(1)),(Space4(10501+ 5*nchan),ftab(1))
     $,(Space4(10501+ 6*nchan),vtab(1)),(Space4(10501+ 7*nchan),mz  (1))
     $,(Space4(10501+ 8*nchan),zi(1,1)),(Space4(10501+20*nchan),cz(1,1))
     $,(Space4(10501+32*nchan),c0(1))
c
      DATA ETA1,SQR4PI/1.0E-20, 3.5449077/
C
      REWIND 3
      NCUP=NLI(2)
      IF(JTRT(1).EQ.0.AND.IS(1).EQ.0) THEN
        NPAR=1
        LDEL=2
      ELSE
        NPAR=2
        LDEL=1
      ENDIF
      IJMAX=0
      NLMAX=0
C
C     CHECK STORAGE SIZE FOR D AMPLITUDE SPACE
C
      LK=0
      NLL=0
      DO 20 I=1,NLTR
      M=(JTRT(I)+1)*NS(I)
      LK=LK+M
      NLL=MAX0(NLL,M)
      DR2(I)=DRF(I)**2
c     set flags for gamma ray channels
      if(fm(i) .eq. 0.0 .and. is(i) .eq. 2) then
        igam(i)=.true.
      else
        igam(i)=.false.
      endif
c     set flags for identical particle
      if(fm(i).eq.fma(i).and.z(i).eq.za(i).and.is(i).eq.jtrt(i)) then
        isym(i) = .true.
      else
        isym(i) = .false.
      endif
   20 CONTINUE
      M=(JTRT(1)+1)*NS(1)
      LK=(LK+LK)*M*LMAX
      NLL=(NLL+NLL)*M
      if(icon(19).eq.0) then
      WRITE(ioutp,'(24H0Number of D ELEMENTS  =,I7)')LK
      endif
C
      IF(NLL.gt.MAXSR) then
        WRITE(ioutp,9003)MAXSR,NLL
 9003   FORMAT('0 CHANNEL MUTIPLICITIES TOO LARGE - REDIMENSION'
     1,  ' SR IN SUBROUTINE XSECT (PRESENT MAXSR = ',I4,'), NEED',I6)
C
C       ERROR 998
        IOUT=998
        rewind 3
        return
      ENDIF 
C
C     COUNT MAXIMUM NUMBER OF CHANNELS
C
      NLMAX=0
      DO 30 I=1,NLTR
      NLMAX=NLMAX+((JTRT(I)+1)*NS(I)+1)/2
   30 CONTINUE
      IF(NLMAX.GT.NCHAN) THEN
        IOUT=-NLMAX
        GO TO 2000
      ENDIF
C
C     CONSTRUCT EFFECTIVE COUPLING MATRIX
C
      DO 32 I=1,NLTR
      INN(I,I)=INN(I,I)+1
   32 CONTINUE
      DO 36 I=1,NLTR
      DO 34 J=1,NLTR
      DO 34 K=1,NLTR
      MK=0
      DO 33 M=1,NLTR
      MK=MK+INN(J,M)*INN(M,K)
      IF(MK.GT.2) MK=2
   33 CONTINUE
      IF(MK.GT.1) MK=2
      ITEMP(J,K)=MK
   34 CONTINUE
      DO 35 J=1,NLTR
      DO 35 K=1,NLTR
      INN(J,K)=ITEMP(J,K)
   35 CONTINUE
   36 CONTINUE
      DO 37 I=1,NLTR
      INN(I,I)=INN(I,I)-1
   37 CONTINUE
C
C
      JX1=(1.0-PHASEF(JTRT(1)+IS(1)))/2.0
      JMAX=LMAX
      LMAX22=LMAX+LMAX-2
c
c     set up coulomb functions
c
      delr=0.1
      maxl=max0(lmax,3)
      iclosed=.false.
      flclos = 0.
      ind=1
      do 40 i=1,nltr
      x2=drf(i)*float(kt)*fk(i)
      x1=x2-drf(i)*fk(i)
      eta_dp = eta(i)
      if(fk2(i).gt.0.0) then
        call cou_dp(x1,x2,eta_dp,maxl,delr,f(ind),fp(ind),g(ind),gp(ind)
     1, s(ind))
      else
        call han_dp(x1,x2,eta_dp,maxl,delr,f(ind),fp(ind),g(ind),gp(ind)
     1, s(ind))
        iclosed=.true.
        flclos=max(flclos,fk(i))
      endif
      ind=ind+lmax
   40 continue
C
C      JA LOOP
C
        if(icon(19).eq.0) write(*,'(a)')
      DO 1000 JA=1,JMAX
        if(icon(19).eq.0) 
     1      write(*,'(''+Integ,  JA = '',i3,''/'',i3)')ja,lmax
      JT=JA+JA-2+JX1
      LMIN=IABS(IABS(JT-JTRT(1))-IS(1))/2
      DO 1000 JPAR=1,NPAR
C
C     SELECT CHANNELS
C
      NLL=0
      MAXL=0
      DO 80 I=1,NLTR
      NLI(I)=0
      JMXX=     JT+JTRT(I)
      JMNN=IABS(JT-JTRT(I))
      IF(JMXX.LT.JMNN) GO TO 80
      DO 70 J=JMNN,JMXX,2
      LMXX=MIN0(J+IS(I),LMAX22)/2
      LMNN=IABS(J-IS(I))/2
      IF(LMXX.LT.LMNN) GO TO 70
      DO 60 L=LMNN,LMXX
      PHAS=PHASEF(LMIN+JPAR-1-L)*FLOAT(LTRT(1)*LTRT(I))
      IF(PHAS.gt.0.0)then 
        NLL   =NLL   +1
        NLI(I)=NLI(I)+1
        LTAB(NLL)=L
        JTAB(NLL)=J
        NTAB(NLL)=I
        FTAB(NLL)=L*(L+1)
        VTAB(NLL)=0.5*(FLOAT(J*(J+2)-IS(I)*(IS(I)+2))/4.0
     1           -FLOAT(L*(L+1)))
        MTAB(NLL)=(I-1)*KMAX
      ENDIF
   60 CONTINUE
   70 CONTINUE
   80 CONTINUE
c
      ncol=nli(1)
      write(3) nll,nli,(ntab(i),ltab(i),jtab(i),i=1,nll)
      IF(NLI(1).EQ.0) GO TO 999
      IF(NLL.GT.NCHAN) IOUT=NLL
      IF(IOUT.NE.0) GO TO 2000
      MAXA=NLL+NLL
      NLMAX=MAX0(NLMAX,NLL)
C
C      COEFFICENT TABLE
C
C     IJ WILL BE CHANNEL COUPLING COUNT
C
      IJ=0
C
C     D(op)*PHI(J) = Sum(I)[C(J,I)*V(J,I)*PHI(I)]
C
      IF(NCUP.ne.0) then
C     I IS *INITIAL* STATE ( ON rhs OF EQUATION )
      DO 340 I=1,NLL
      II=NTAB(I)
C     J IS * FINAL * STATE ( ON lhs OF EQUATION )
      DO 340 J=1,NLL
      JJ=NTAB(J)
      DO 330 K=1,NCUP
      IF(ICUP(K).eq.II.and.JCUP(K).eq.JJ) then
C
      LX =LTAB(I)*2
      JX =JTAB(I)
      IX =JTRT(II)
      ISX=IS  (II)
      KLII=KLAM(II)
      fmii=fma(ii)
C
      LY =LTAB(J)*2
      JY =JTAB(J)
      IY =JTRT(JJ)
      ISY=IS  (JJ)
      KLJJ=KLAM(JJ)
      fmjj=fma(jj)
c
      LK=LKTAB(K)
      LL2=LAM(LK)*2
      JJ2=JSP(LK)
      IS2=ISP(LK)
      icode = ilk(lk)
c
c     gamma channel 
c
      gtemp = 0.0
c     initial channel is gamma
      if(igam(ii)) then
        isi=0
c       calculate effective charges , M/(m+M) factor is in dr scaling
        gtemp = z(jj) + za(jj)*(-fm(jj)/fma(jj))**ltab(i)
        if(icode .ne. 16 .and. ltab(i) .ne. 0)
     1     gtemp = gtemp*sqrt(float(ltab(i)+1)/float(ltab(i)))
     2           *sqr4pi
      else
        isi=isx
        gtemp = 1.0
      endif
c     final   channel is gamma
      if(igam(jj)) then
        isf=0
c       calculate effective charges , M/(m+M) factor is in dr scaling
        gtemp = z(ii) + za(ii)*(-fm(ii)/fma(ii))**ltab(j)
        if(icode .ne. 16 .and. ltab(j) .ne. 0) 
     1     gtemp = gtemp*sqrt(float(ltab(j)+1)/float(ltab(j)))
     2           *sqr4pi
      else
        isf=isy
        gtemp = 1.0
      endif
c
c     set up for identical particles
c
      if(isym(ii)) then
        ispp  =isi
        jmin2 = iabs(jt-ispp)
        jmax2 =      jt+ispp
        jpy   = jy
        jpp   = jx
        lp    = lx
        isymp = .true.
        jpflag= .true.
      elseif(isym(jj).and..not.isym(ii)) then
        ispp  = isf
        jmin2 = iabs(jt-ispp)
        jmax2 =      jt+ispp
        jpx   = jx
        jpp   = jy
        lp    = ly
        isymp = .true.
        jpflag= .false.
      else
        jmin2 = jx
        jmax2 = jx
        jpy   = jy
        isymp = .false.
        jpflag= .true.
      endif
c
      scoeff = 0.0
      do 90 jp=jmin2,jmax2,2
      if(jpflag) then
        jpx = jp
      else
        jpy = jp
      endif
C
      TEMP=gtemp*RACAH(jpx,IX,jpy,IY,JT,JJ2)
     1    *SQRT(FLOAT(jpx+1)*FLOAT(jpy+1)*FLOAT(JJ2+1))
     2    *WINEJ(LX,isi,jpx,LL2,IS2,JJ2,LY,isf,jpy)
     3    *PHASEF((JT-IX-jpy-JJ2)/2)
      if(isymp) then
c       write(20,'(a,5i3,)')' lp, jp, jpp, jt, ispp',jt,lp,jp,jpp,ispp
        tcoeff = sqrt(float((jp+1)*(jpp+1)))*phasef(jt/2)
     1          *Racah(lp,jp,jpp,jt,ispp,ispp)
        if(jp.eq.jpp) tcoeff = tcoeff + 1.0 
c       write(20,'(a,36x,2f10.6)'+',tcoeff,temp
        scoeff = scoeff + temp*tcoeff
      else
        scoeff = temp 
      endif
   90 continue
c     write(20,'(a,f10.6)')' scoeff =',scoeff
      temp = scoeff
c	------------------------------------------
      call icode3(temp)
c	------------------------------------------
      T1= MAX (DRF(II),DRF(JJ))
      T2= MIN (DRF(II),DRF(JJ))
      TEMP=TEMP*SQRT(T1/T2)**3
      IF(TEMP.ne.0.0) then
      IJ=IJ+1
      IF(IJ.gt.MAXIJ) then
C
C       ERROR 888
        IOUT=888
        WRITE(ioutp,9004)MAXIJ
 9004   FORMAT('H0TOO MANY COUPLINGS BETWEEN CHANNELS, MAX =',I4 )
C
        rewind 3
        return
      ENDIF
      IKTAB(IJ)=I
      JKTAB(IJ)=J
      LKTLK(IJ)=(LK-1)*KMAX
      if(igam(jj)) then
        ctab(ij)   =temp        *drf(jj)**3/drf(ii)
      else
        ctab(ij)   =temp*fmu(jj)*drf(jj)**3/drf(ii)
      endif
      endif
      endif
  330 CONTINUE
  340 CONTINUE
      endif
c
      IJMAX=MAX0(IJMAX,IJ)
c
c     start integration at an l-dependent value of r to circumvent the
c     angular momentum barrier problem
c
      DO 430 N2=1,NLL
      II=NTAB(N2)
C     MZ(I) IS INTEGRATION STARTING POINT FOR CHANNEL I
      M1=FLOAT(LTAB(N2)-10)/(2.0*FK(II)*DRF(II))+1.0
      m1=min0(m1,kt-3)
      mz(n2)=min0(max0(m1,1),kmax)
  430 CONTINUE
c     **********************************************************
      if(iclosed) then
        kint = (1.25*fma(1)*0.333333 + 5.0/flclos)/drf(1)
        kint=max(kint,kt)
      else
        kint = kt
      endif
C
C     ITERATION LOOP
C
c     nloop = 1  integrate all channels outward to matching points
c
c     nloop = 2  integrate all channels  inward to matching points
c                  starting from an outgoing hankel function in each
c     nloop = 3  integrate initial channels inward to matching points
c                  starting from a regular function in each
c
      do 900 nloop=1,3
      if    (nloop.eq.1) then
        ma=2
        mb=kint
        mdel=+1
        n1cnt=nll
      elseif(nloop.eq.2) then
        ma=kt-2
        mb=kint-1
        mdel=-1
        n1cnt=nll
      else
        ma=kt-2
        mb=kint-1
        mdel=-1
        n1cnt=ncol
      endif
c
c	if(ja.lt.3) write(20,'(a,4i5)')
c     1	' ma, mb, mdel, nlcnt =',ma,mb,mdel,n1cnt
      DO 800 N1=1,n1cnt
      NJ=N1+N1-1
      JJ=NTAB(N1)
C
C      SET UP INITIAL Z  VALUES
C
      NMM=0
      do 440 n2=1,nll
      ii=ntab(n2)
C     SEARCH DOWN A COLUMN TO TO SEE IF THESE CHANNELS ARE COUPLED
      IF(INN(II,JJ).NE.0.OR.N1.EQ.N2) THEN
        NMM=NMM+1
        NTBL(NMM)=N2
      ENDIF
  440 continue
C
C     CONSTRUCT SECONDARY TABLES
C
C     IK WILL BE CHANNEL COUPLING COUNT FOR THIS ITERATION
C
C
      IK=0
      MZMAX=1
      IF(IJ.NE.0) THEN
      DO 455 I=1,NMM
      I1=NTBL(I)
      DO 450 J=1,NMM
      J1=NTBL(J)
C     SEARCH THRU CHANNEL COUPLINGS FOR A MATCH
      DO 445 K=1,IJ
      IF(JKTAB(K).EQ.NTBL(J).AND.IKTAB(K).EQ.NTBL(I)) THEN
        IK=IK+1
        ICTBL(IK)=K
        if(nloop.eq.1) then
          MZMAX=MAX0(MZMAX,MZ(J1))
C
C         If mj > mi then no coupling until m > mj
C         A negative ictbl value inhibits the coupling
C
          IF(MZ(J1).GT.MZ(I1)) ICTBL(IK)=-K
        endif
      ENDIF
  445 CONTINUE
  450 CONTINUE
  455 CONTINUE
      ENDIF
      MZMAX=MZMAX+1
c
      if(nloop.eq.1) then
C
C     SORT STARTING POINTS IN ASCENDING ORDER FOR CHANNELS
C
      IF(IK.NE.0) THEN
      DO 470 I=1,IK
      M1=KMAX
      K1=I
C     FIND SMALLEST MI
      DO 460 K=I,IK
      K2=IABS(ICTBL(K))
      J= JKTAB(K2)
      M2=MZ(J)
      IF(M2.LT.M1) THEN
        M1=M2
        K1=K
      ENDIF
  460 CONTINUE
C     INTERCHANGE THE ENTRIES
      M2=ICTBL(K1)
      ICTBL(K1)=ICTBL(I )
      ICTBL(I )=M2
  470 CONTINUE
      ENDIF
      ICTBL(IK+1)=0
      endif
c
c     SET UP INITIAL CONDITIONS FOR THIS ITERATION
c
      DO 480 N2 = 1,NLL
      II=NTAB(N2)
      ZI(N2,1)=0.0
      ZI(N2,2)=0.0
      ZI(N2,3)=0.0
      CZ(N2,2)=0.0
      CZ(N2,3)=0.0
      IF(N1.EQ.N2) THEN
c
      if    (nloop.eq.1) then
        mz0=mz(ii)
        mz1=mz0
        zi(n2,3)=eta1
c
c     l = 1 is special case at origin, cz is not zero for r = 0
c
        if(ltab(n2).eq.1) then
          cz(n2  ,2) = 2.0*eta1
        endif
      elseif(nloop.eq.2) then
        mz0=ma+1
        ind=(ntab(n2)-1)*lmax+ltab(n2)+1
        zi(n2,2)=dcmplx(gp(ind),fp(ind))
        zi(n2,3)=dcmplx(g (ind),f (ind))
      else
        mz0=ma+1
        ind=(ntab(n2)-1)*lmax+ltab(n2)+1
        zi(n2,2)=fp(ind)
        zi(n2,3)=f (ind)
      endif
        mk=mtab(n2)+mz0
C
C     COMPUTE DIAGONAL PART OF CZ
C
c	write(*,'(1h+,20x,6i4)')nloop,ii,ma,mz0,mdel,mz0-mdel
        if(nloop.ne.1) then
        CZ(N2,2) =(FTAB(N2)/FLOAT(MZ0-mdel)**2
     1       - DR2(II)*(W(MK-mdel)+V(MK-mdel)*VTAB(N2)))*zi(n2,2)
        endif
        CZ(N2,3) =(FTAB(N2)/FLOAT(MZ0     )**2
     1       - DR2(II)*(W(MK     )+V(MK     )*VTAB(N2)))*zi(n2,3)
      ENDIF
  480 CONTINUE
C
C     COMPUTE NON-DIAGONAL PART OF CZ
C
c	if(ja.lt.3) write(20,'(a,4i6,1p4e12.4)')
c     1	' ja, n1, nloop, mz0 =',ja,n1,nloop,mz0
c
      IF(IK.NE.0) THEN
        DO 490 NY=1,IK
        IC=ICTBL(NY)
        IF(IC.GT.0) THEN
          I1=IKTAB(IC)
          NK=IKTAB(IC)
          NL=JKTAB(IC)
          LK=LKTLK(IC)+mz0
          if(nloop .ne. 1)then
            x(1) = XC1(lk-mdel      )*ctab(ic)
            CZ(NL,2)=CZ(NL,2)-X(1)*ZI(NK,2)
cc
          endif
c
          x(1) = XC1(lk      )*ctab(ic)
          CZ(NL,3)=CZ(NL,3)-X(1)*ZI(NK,3)
cc
        ENDIF
  490   CONTINUE
      ENDIF
c	if(ja.lt.3) write(20,'(a,4i5)')
C
C      INTEGRATION SECTION
C
C      INTEGRATE U(R) BY BAYLIS-PEEL PREDICTOR-CORRECTOR NUMEROV METHOD
C              COMP PHYS COMM 25 (1982) 21.
C
      M1=1
      DO 600 M=ma,mb,mdel
      MK=M+M-1
c
      IF(IK.NE.0) THEN
      DO 510 NY=1,IK
      IC=ICTBL(NY)
      IF(IC.GT.0) THEN
        LK=  LKTLK(IC)+M
        X(NY) = XC1(lk       )*ctab(ic)
      ENDIF
  510 CONTINUE
      ENDIF
C
      IF(M.LT.MZMAX) THEN
  520 CONTINUE
      IF(ICTBL(M1).NE.0) THEN
        K1=IABS(ICTBL(M1))
        J=JKTAB(K1)
        M2=MZ(J)
        IF(M.GE.M2   ) THEN
C         COUPLING IS TURNED ON AFTER M .GE. MZ(J)
C         FOR HIGH PARTIAL WAVES IN SOURCE TERM
          ICTBL(M1)=IABS(ICTBL(M1))
          M1=M1+1
          GO TO 520
        ENDIF
      ENDIF
      ENDIF
      IF(M.LT.MZ1) GO TO 600
C
      DO 530 NY=1,NMM
      N2=NTBL(NY)
C
C     PREDICTION - INTEGRATE U(R) BY N=2 STORMER EQUATION
C
C     STEP INTEGRATION TABLES BACK ONE
C
      CZ(N2,1)=CZ(N2,2)
      CZ(N2,2)=CZ(N2,3)
      ZI(N2,1)=ZI(N2,2)
      ZI(N2,2)=ZI(N2,3)
      ZI(N2,3) = 2.0*ZI(N2,2)-ZI(N2,1)+CZ(N2,2)
  530 CONTINUE
C
      DO 580 ITER=1,2
C
C      COMPUTE NEW CZ VALUE
C
C     COMPUTE DIAGONAL PART OF CZ
C
      DO 540 NY=1,NMM
      N2=NTBL(NY)
      JJ=NTAB(N2)
      NM=MTAB(N2)+M
      CZ(N2,3) =(FTAB(N2)/FLOAT(M)**2
     1         - DR2(JJ)*(W(NM)+V(NM)*VTAB(N2)))*ZI(N2,3)
  540 CONTINUE
C
C     COMPUTE NON-DIAGONAL PART OF CZ
C
      IF(IK.NE.0) THEN
      DO 550 NY=1,IK
      IC=ICTBL(NY)
      IF(IC.GT.0) THEN
        NK=IKTAB(IC)
        NL=JKTAB(IC)
        CZ(NL,3)=CZ(NL,3)-X(NY)*ZI(NK  ,3)
      ENDIF
  550 CONTINUE
      ENDIF
C
      IF(ITER.EQ.1) THEN
C
C     FIRST CORRECTION
C
        DO 560 NY=1,NMM
        NK=NTBL(NY)
        ZI(NK,3)=ZI(NK,3)+(CZ(NK,3)-2.0*CZ(NK,2)+CZ(NK,1))/12.
        C0(NK)  =CZ(NK,3)
  560   CONTINUE
      ELSE
C
C     SECOND CORRECTION
C
        DO 570 NY=1,NMM
        NK=NTBL(NY)
        ZI(NK,3)=ZI(NK,3)+(CZ(NK,3)-C0(NK))/30.0
  570   CONTINUE
      ENDIF
c	if(ja .eq. 1) then
c	write(20, '(3i4,1p8e15.6)')jpar,nloop,m,(zi(1,n),n=1,3)
c	endif
  580 CONTINUE
  600 CONTINUE
C
C     END OF INTEGRATION LOOP
C
C
C     WRITE LAST TWO POINTS ON TAPE
C
      if(nloop.eq.1) then
        write(3)(zi(i,2),zi(i,3),i=1,nll)
c	if(ja .eq. 1) then
c	write(20,'(a)')' **************************'
c	write(20, '(a,2i2,1p8e15.7/(12x,1p8e15.7))')
c     1       ' ja, n2=',ja,n1,(zi(i,2),zi(i,3),i=1,nll)
c	endif
      else
        write(3)(zi(i,3),zi(i,2),i=1,nll)
c	if(ja .eq. 1) then
c	write(20, '(a,2i2,1p8e15.7/(12x,1p8e15.7))')
c     1       ' ja, n2=',ja,n1,(zi(i,3),zi(i,2),i=1,nll)
c	endif
      endif
  800 CONTINUE
C
C     END OF ITERATION LOOP
C
C     CHECK ON CHANNEL COUPLINGS
C
      IF(IK.NE.0) THEN
        DO 890 K=1,IK
        IF(ICTBL(K1).LT.0) THEN
          K1=IABS(ICTBL(K))
          WRITE(6,9020)JA,JPAR,nloop,JKTAB(K1),IKTAB(K1)
        ENDIF
  890   CONTINUE
      ENDIF
  900 continue
c      write(20,'(a)')' **************************'
c      do 920 n2=1,nltr
c      ind=(ntab(n2)-1)*lmax+ltab(n2)+1
c      write(20,'(a,i3,1p10e12.4/(12x,1p10e12.4))')' g, f, gp,fp',ind
c     1    ,g (ind),f (ind),gp(ind),fp(ind)
c  920 continue
c      write(20,'(a)')' **************************'
  999 CONTINUE
 1000 CONTINUE
 2000 CONTINUE
      REWIND 3
      if(icon(19).eq.0) then
      WRITE(ioutp,9005)NLMAX,IJMAX
      endif
c
      RETURN
c
 9005 FORMAT(' MAX. NO. CHANNELS     = ',I6
     1,    '   MAX. NO. CHANNEL COUPLINGS =',I6)
 9020 FORMAT(' WARNING ***, COUPLING FOR JA =',I3,'  JPAR =',I2,
     1  '  Nloop =',i3,'  J =',I3,'  I =',I3,'  NOT REALIZED ***')
      END

      SUBROUTINE COU_dp(R,RP,ETA,L,H,F,FP,G,GP,S)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION F(51),FP(51),G(51),GP(51),S(51)
C
      LL=L
      IF(LL.LT.50) THEN
        ELP=50.
        J=50
      ELSE
        ELP=LL
        J=LL
      ENDIF
      A=ATAN (ETA/ELP)
      B=SQRT (ETA**2+ELP**2)
      Y=A*(ELP-0.5)+ETA*( LOG(B)-1.)-SIN (A)/(12.*B)
     1  +SIN (3.*A)/(360.*B**3)-SIN (5.*A)/(1260.*B**5)
     2  +SIN (7.*A)/(1680.*B**7)-SIN (9.*A)/(1188.*B**9)
C
      K=J-1
      DO 100 I=1,K
      IF(J.LE.LL) S(J)=Y
      ELP=ELP-1.
      J=J-1
      Y=Y-ATAN (ETA/ELP)
  100 CONTINUE
      S(1)=Y
      DEL1=R-2.0*ETA
      RMAX=.41666667*(ETA**2+4.*ETA+3.)
      IF(RMAX.LT.10.0) RMAX=10.0
      DEL=R-RMAX
      IF(ETA.GE.5.0.AND.ABS (DEL1).LT.ABS (DEL)) THEN
      DEL=DEL1
        I1=2
C     Transition line expansion
      X=2.0*ETA
      T1=ETA**2
      T2=T1**2
      T3=ETA** .666666667
      T4=T3**2
      T5=T4**2
      T6=T3*T5
      T7=T4*T6
      T8=T3*T7
      T9=ETA** .166666667
      Y=1.22340402*T9*(1.+.495957017E-1/T4-.888888889E-2/T1
     1  +.245519918E-2/T6-.910895806E-3/T2+.845362E-3/T8)
      Z=-.707881773/T9*(1.-.172826039/T3+.317460317E-3/T1
     1  -.358121485E-2/T5+.311782468E-3/T2-.907396643E-3/T7)
      ELSE
C
      IF(DEL.GE.0.0.OR.ETA.EQ.0.0) THEN
        X=R
        I1=1
      ELSE
        X=RMAX
        I1=2
      ENDIF
C     Asymptotic expansion
      T1=ETA**2
      T2=2.*X
      T4=ETA/T2
      SS=1.
      TS=0.
      SL=0.
      TL=1.-ETA/X
      SSS=1.
      STS=0.
      SSL=0.
      STL=TL
      EN=0.
      DO 700 K=1,25
      IF(ABS (SS/SSS).GT.1.0E-10) THEN
        T5=EN+1.
        T6=T5+EN
        T7=EN*T5
        T8=T6*T4/T5
        T9=(T1-T7)/(T2*T5)
        T5=T8*SS-T9*TS
        TS=T8*TS+T9*SS
        SS=T5
        T5=T8*SL-T9*TL-SS/X
        TL=T8*TL+T9*SL-TS/X
        SL=T5
        SSS=SSS+SS
        STS=STS+TS
        SSL=SSL+SL
        STL=STL+TL
        EN=EN+1.
      ENDIF
  700 CONTINUE
C
      T3=X-ETA* LOG(T2)+S(1)
      T8=SIN (T3)
      T9=COS (T3)
      Y=SSS*T9-STS*T8
      Z=SSL*T9-STL*T8
      ENDIF
C
      DO 810 I=1,I1
      IF(I.EQ.I1) THEN
        G(1) = Y
        W    = Z
        DEL = RP - R
      ENDIF
C     Runge-Kutta integration
      N=ABS (DEL/H)+1.0
      DX=DEL/FLOAT(N)
      T1=DX/2.
      T2=DX/8.
      T3=2.0*ETA
      DO 805 J=1,N
      T4=DX*(T3/X-1.)*Y
      X=X+T1
      T5=DX*(T3/X-1.)*(Y+T1*Z+T2*T4)
      X=X+T1
      T6=DX*(T3/X-1.)*(Y+DX*Z+T1*T5)
      Y=Y+DX*(Z+(T4+2.*T5)/6.)
      Z=Z+(T4+4.*T5+T6)/6.
  805 CONTINUE
  810 CONTINUE
      GP(1)=Y
C
      T1=ETA**2
      T8=SQRT (1.+T1)
      G(2)=((1./R+ETA)*G(1)-W)/T8
      GP(2)=((1./RP+ETA)*Y-Z)/T8
      T2=1.
      T3=2.
C     Recur irregular function upwards
      DO 910 I=3,LL
      T4=T2+T3
      T5=T2*T3
      T6=T3*SQRT (T2**2+T1)
      T7=T2*SQRT (T3**2+T1)
      G (I)=(T4*(ETA+T5/R )*G (I-1)-T6*G (I-2))/T7
      GP(I)=(T4*(ETA+T5/RP)*GP(I-1)-T6*GP(I-2))/T7
      T2=T2+1.
      T3=T3+1.
  910 CONTINUE
      N=MAX0(L+11, INT(2.0*R+11.0) )
      I1=N
      Y =1.0E-20
      YP=Y
      Z =0.
      ZP=Z
C     Recur regular function downwards
      DO 930 I=1,I1
      T2=N
      T3=T2+1.
      T4=T2+T3
      T5=T2*T3
      T6=T2*SQRT (T3**2+T1)
      T7=T3*SQRT (T2**2+T1)
      X =(T4*(ETA+T5/R )*Y -T6*Z )/T7
      XP=(T4*(ETA+T5/RP)*YP-T6*ZP)/T7
      IF(N.LE.LL) THEN
        F (N)=X
        FP(N)=XP
      ELSE
        IF(ABS (X).GT.1.0) THEN
          Y =Y *1.0E-20
          YP=YP*1.0E-20
          X =X *1.0E-20
          XP=XP*1.0E-20
        ENDIF
      ENDIF
      N=N-1
      Z =Y
      ZP=YP
      Y =X
      YP=XP
  930 CONTINUE
      Z =1./(T8*(F (1)*G (2)-F (2)*G (1)))
      ZP=1./(T8*(FP(1)*GP(2)-FP(2)*GP(1)))
      DO 950 I=1,LL
      F (I)=F (I)*Z
  950 FP(I)=FP(I)*ZP
      RETURN
       END

c***********************************************************************
      SUBROUTINE HAN_dp(R ,RP,E,L,H,F,FP,G,GP,S)
c
c     Subroutine to calculate a Hankel function for matching to a
c     closed channel
c     
c***********************************************************************
      implicit real*8(a-h,o-z)
      dimension F(50),FP(50),G(50),GP(50),S(50)
	common/ch3files/input, ioutp
C
      FF(X)=EXP(-X-E* LOG(2.0*X))
      DR=H
      DO 40 J=1,2
      FL=J-1
      T3=1.0
      T4=1.0
      T5=1.0
      T6=1.0
      G (J)=1.0
      GP(J)=1.0
      DO 30 I=1,20
      T3=T3*(FL+1.0-E-T4)*(FL+E+T4)/(2.0*T4)
      T5=T5/R
      T6=T6/RP
      G (J)=G (J)+T5*T3
      GP(J)=GP(J)+T6*T3
      T4=T4+1.0
      TEST=ABS(T5*T3)
      IF(TEST.LT.1.0E-05) GO TO 35
   30 CONTINUE
      WRITE(ioutp,9000)TEST
   35 CONTINUE
      G (J)=G (J)*FF(R )
      GP(J)=GP(J)*FF(RP)
   40 CONTINUE
      F (1)=0.0
      F (2)=0.0
      FP(1)=0.0
      FP(2)=0.0
      S(1)=0.0
      S(2)=0.0
      LL=L
      T1=E*E
      T2=1.0
      T3=2.0
      DO 50 I=3,LL
      T4=T2+T3
      T5=T2*T3
      T6=T3*SQRT(T2**2+T1)
      T7=T2*SQRT(T3**2+T1)
      G (I)=(T4*(E+T5/R )*G (I-1)+T6*G (I-2))/T7
      GP(I)=(T4*(E+T5/RP)*GP(I-1)+T6*GP(I-2))/T7
      T2=T2+1
      T3=T3+1
      F (I)=0.0
      FP(I)=0.0
      S(I)=0.0
   50 CONTINUE
      RETURN
 9000 FORMAT(25H0BOUNDSTATE CONVERGES TO ,1PE12.3,8H IN HAN   )
      END

