      program main

      use efvs
      implicit none

c  namelist variables
      integer          :: iunit1,Jmax,ny,ndelta,lrpow,ilevel
      double precision :: ymin,dy,deltamin,ddelta,offset
      logical          :: adiamn,Sdirect,Sindirect,printS,printa,
     1                    printsigpar,printcheck,printS2mat
c  internal loop variables
      integer          :: ibasis,iy,jbasis,idelta,iopen,jopen,nopen,
     1                    ichan,jchan,nopsr,noplr,ntwo,ilev,jlev,
     2                    ifrom,ito,icalc,ift,ifdex,itdex,iclosed,
     3                    kdum
c  internal arrays
      integer          :: Jtotal,nbasis,id,noo,noi,idlev,ifindx
      integer          :: mxcalc,mxstate,mxopen,ncalc,mxlev
      integer          :: levs,nlevs,iblock,iecalc,ielev
      double complex   :: zmat1,zmat2,zmat3
      double complex   :: z2mat1,z2mat2,s2mat
      double complex   :: r_io,r_oi,t_io,t_oi
      double complex   :: Sc
      double complex   ::  xi_R_plus, xi_L_plus, xi_R_minus, xi_L_minus,
     1                    dxi_R_plus,dxi_L_plus,dxi_L_minus,dxi_R_minus
c  internal variables
      integer          :: nlev,jdum(1),idum(1),ldum(1),ifail,ioff,joff,
     1                    ll
      double precision :: y,deltas,edum(1)
      double precision :: abar,pi,phi,awvmax
      double precision :: k,ssquare,simag,sreal,sphase,sqrtk
      double precision :: uj,ujp,un,unp,ratio,ell,darg,xj,xjfac
      double precision :: sigma,dg,dd,wv,sigcur,stemp,slcur,sloss,
     1                    sum_inelastic,sj,sjp,sn,snp,sigpar,slpar
      double precision :: ssum,ssuminel
      double complex   :: iu,phase,pphase,scatlen

c  values read in from iunit1
      integer          :: itype,nstate,nqn,iprogm,jstate,jhalf,ibound,
     1                    nfield,nnrg,nlevel,ndgvl,nconst,nrsq,n,li
      integer          :: jtot,inrg,ib,jfield,iexch,nsig
      double precision :: ured,energy,elevel,rmax,runit
      double precision :: enow,eref,wght,cento,wval
      double precision :: wi
      double precision :: fp0,fpn,gp0,gpn
      double precision :: evecs,evals,degin,degout
      character(80)    :: label
      character(10)    :: runame


      parameter (iu=(0.0,1.0))
      parameter (awvmax=0.01d0)
      parameter (mxcalc=5,mxstate=1003,mxopen=303)
      external zinv

c  arrays used in preparation, or for both methods
      allocatable      :: Jtotal(:),nbasis(:),id(:,:),noo(:),noi(:),
     1                    degin(:),degout(:),wv(:,:),jstate(:,:),
     2                    elevel(:),li(:),cento(:),fp0(:,:),fpn(:,:),
     3                    gp0(:,:),gpn(:,:),wi(:,:),evecs(:,:),
     4                    evals(:),Sc(:,:),energy(:),mxlev(:),
     5                    ielev(:,:),ll(:,:)
c  functions in (O), (O'), (I) and (I') matrices
      allocatable      :: xi_R_plus(:),xi_L_plus(:),xi_R_minus(:),
     1                                                    xi_L_minus(:),
     2                    dxi_R_plus(:),dxi_L_plus(:),dxi_L_minus(:),
     3                                                   dxi_R_minus(:)
c  arrays used for indirect method
      allocatable      :: zmat1(:,:),zmat2(:,:),zmat3(:,:),z2mat1(:,:),
     1                    z2mat2(:,:),s2mat(:,:),r_io(:,:,:),
     2                    r_oi(:,:,:),t_io(:,:,:),t_oi(:,:,:)
c  arrays used for direct method
      double complex, allocatable :: y0(:,:),temp1(:,:),temp2(:,:),
     1                               temp3(:,:),temp4(:,:),temp5(:,:),
     2                               yn(:,:),S(:,:)
c  arrays used in calculation of cross sections
      allocatable :: sigma(:,:,:),levs(:),nlevs(:),sigcur(:,:),
     1               iblock(:),iecalc(:),wght(:),sloss(:,:),
     2               slcur(:),sum_inelastic(:),ifindx(:),
     3               sigpar(:,:,:,:),slpar(:,:,:),
     4               ssuminel(:)

      namelist /input/ adiamn, iunit1, ymin, dy, ny, deltamin,
     1                 ddelta, ndelta, abar, offset, ichan, Sdirect,
     2                 Sindirect,printS,Jmax,lrpow,printa,printsigpar,
     3                 ilevel,printcheck,printS2mat

      pi=2.D0*dacos(0.D0)
c  lazy programming assumes that input is correct
      adiamn=.false.
      iunit1=21
      ymin=0.0
      dy=0.1
      ny=11
      deltamin=0.0
      ddelta=0.01
      ndelta=100
      abar=1.0
c  crls:  deltas=+/-0.5 (i.e., pi) & y=0 corresponds to Y0=I*infinity
c                                                  (i.e. molscat Y0).
c         offset can be used to shift input values so that deltas corresponds
c         corresponds to input value + offset
      offset=0.5
      ichan=1
      Sdirect=.false.
      Sindirect=.true.
      printS=.false.
      printa=.false.
      printcheck=.false.
      printS2mat=.false.
      Jmax=-1
      lrpow=6
c  bm addition:
      ilevel=1
      printsigpar=.false.

      read(5,input)
      call iaread(iunit1,label,itype,nstate,nqn,ured,iprogm,
     1            jdum,jhalf,ibound,nfield,nnrg,edum,nlevel,edum,
     2            ndgvl,nconst,nrsq,nsig,edum,edum,runit,runame,.true.)
      allocate (jstate(nstate,nqn),elevel(nlevel),energy(nnrg),
     1          degin(nsig),degout(nsig))
      call iaread(iunit1,label,itype,nstate,nqn,ured,iprogm,
     1            jstate,jhalf,ibound,nfield,nnrg,energy,nlevel,elevel,
     2            ndgvl,nconst,nrsq,nsig,degin,degout,runit,runame,
     3            .false.)
      nlevel=nstate

      allocate (r_io(mxstate,mxstate,mxcalc),
     1          r_oi(mxstate,mxstate,mxcalc),
     2          t_io(mxstate,mxstate,mxcalc),
     3          t_oi(mxstate,mxstate,mxcalc))
      allocate (Jtotal(mxcalc),wv(mxstate,mxcalc),id(mxstate,mxcalc),
     1          nbasis(mxcalc),noo(mxcalc),noi(mxcalc),nlevs(mxcalc),
     2          levs(mxstate),iblock(mxcalc),iecalc(mxcalc),
     3          wght(mxcalc),ll(mxstate,mxcalc))

      ifail=0
      icalc=0
      do while (ifail.eq.0)
        icalc=icalc+1
        if (icalc.gt.mxcalc) then
          write(6,*)'increase mxcalc!'
          stop
        endif
c end of header, start of loop
        call ialoop(iunit1,jtot,inrg,ib,jfield,enow,eref,iexch,
     1              wght(icalc),n,idum,ldum,.true.,ifail)
        if (ifail.ne.0) exit
        Jtotal(icalc)=jtot
        iblock(icalc)=ib
        iecalc(icalc)=inrg
        nbasis(icalc)=n
        if (n.gt.mxstate) then
          write(6,*) 'increase mxstate! n =',n
          stop
        endif
        allocate (li(n))
        backspace (iunit1)
        call ialoop(iunit1,jtot,inrg,ib,jfield,enow,eref,iexch,
     1              wght(icalc),n,id(1,icalc),li,.false.,ifail)
        deallocate (li)
        call iawvec(iunit1,wv(1,icalc),noo(icalc),n)

        if (noo(icalc).eq.0) then
          icalc=icalc-1
          cycle
        endif

c  read short-range transformation matrix into evecs and eigenvalues into evals
        allocate (evecs(n,n),evals(n))
        call iasr(iunit1,evecs,evals,n)

c  read centrifugal term into cento
        allocate (cento(n))
        call iacent(iunit1,cento,rmax,n)

c  check number of open levels
        levs(:)=0
        noplr=noo(icalc)
        do iopen=1,noplr
          levs(id(iopen,icalc))=1
        enddo
        nlevs(icalc)=sum(levs)
        if (noplr.gt.mxopen) then
          write(6,*)'increase mxopen! noplr =',noplr
          stop
        endif
        allocate (fp0(n,n),fpn(n,n),gp0(n,n),gpn(n,n))
        call iaypsi(iunit1,fpn,fp0,gpn,gp0,n)

        gp0=matmul(gp0,evecs)
        gpn=matmul(gpn,evecs)
        evecs=transpose(evecs)
        fp0=matmul(evecs,fp0)
        gp0=matmul(evecs,gp0)
        deallocate (evecs)

c  Check that input matrices have the right symmetry (fpn,gp0 are symmetric,
c  and fp0=-[gpn]^T).  Uncomment symmetrisation code if required
        if (printcheck)
     1    write(6,*) '  i1   i2   fpn(i1,i2)-fpn(i2,i1) gpn(i1,i2)+'//
     2               'fp0(i2,i1) gpn(i2,i1)+fp0(i1,i2) gp0(i1,i2)-'//
     3               'gp0(i2,i1)'
        do ibasis=1,n
        do jbasis=ibasis,n
c symmetrisation code
c         gp0(ibasis,jbasis)=(gp0(ibasis,jbasis)+gp0(jbasis,ibasis))/2d0
c         gp0(jbasis,ibasis)=gp0(ibasis,jbasis)
c         fpn(ibasis,jbasis)=(fpn(ibasis,jbasis)+fpn(jbasis,ibasis))/2d0
c         fpn(jbasis,ibasis)=fpn(ibasis,jbasis)
c         gpn(ibasis,jbasis)=(gpn(ibasis,jbasis)-fp0(jbasis,ibasis))/2d0
c         fp0(jbasis,ibasis)=-gpn(ibasis,jbasis)
c         gpn(jbasis,ibasis)=(gpn(jbasis,ibasis)-fp0(ibasis,jbasis))/2d0
c         fp0(ibasis,jbasis)=-gpn(jbasis,ibasis)
          if (printcheck)
     1    write(6,2000)ibasis,jbasis,fpn(ibasis,jbasis)-
     2                 fpn(jbasis,ibasis),gpn(ibasis,jbasis)+
     3                 fp0(jbasis,ibasis),gpn(jbasis,ibasis)+
     4                 fp0(ibasis,jbasis),gp0(ibasis,jbasis)-
     5                 gp0(jbasis,ibasis)
        enddo
        enddo
 2000   format(2i5,2x,1p,4(1x,g10.3,11x))

        if (Jtotal(icalc).gt.Jmax .and. Jmax.ge.0) then
          deallocate (fp0,fpn,gp0,gpn)
          deallocate (evals)
          deallocate (cento)
          cycle
        endif
        allocate ( xi_R_plus(n), xi_L_plus(n), xi_R_minus(n),
     1                                                    xi_L_minus(n),
     2            dxi_R_plus(n),dxi_L_plus(n),dxi_L_minus(n),
     3                                                   dxi_R_minus(n))
        xi_R_plus  = 0.d0
        xi_L_plus  = 0.d0
        xi_R_minus = 0.d0
        xi_L_minus = 0.d0
        dxi_R_plus = 0.d0
        dxi_L_plus = 0.d0
        dxi_L_minus= 0.d0
        dxi_R_minus= 0.d0

c  long-range (f^o) boundary conditions
        do ibasis = 1,n
          darg=wv(ibasis,icalc)*rmax
          ell=sqrt(cento(ibasis)+0.25)-0.5
          ll(ibasis,icalc)=nint(ell)
          k=sqrt(abs(wv(ibasis,icalc)))
          if (ibasis.le.noplr) then
c  open channels
            call rbesjy(ell, darg, uj, ujp, un, unp)
            sj = uj /k
            sn = un /k
            sjp= ujp*k
            snp= unp*k
            xi_R_plus(ibasis)  = (sj +sn *iu)
            dxi_R_plus(ibasis) = (sjp+snp*iu)
            xi_R_minus(ibasis) = conjg(xi_R_plus(ibasis))
            dxi_R_minus(ibasis)= conjg(dxi_R_plus(ibasis))
          else
c  closed channels
            call rbessk(ell, abs(darg), ratio)
            xi_R_plus(ibasis)  = 1.0D0   /k
            dxi_R_plus(ibasis) = ratio*k
            xi_R_minus(ibasis) = 0.d0  !conjg(xi_R_plus(ibasis))
            dxi_R_minus(ibasis) = 0.d0 !conjg(dxi_R_plus(ibasis))
          endif
          if (printcheck) then
          write(6,2010) 'Wronskian for long range function ',ibasis,
     1                  dxi_R_plus(ibasis)*xi_R_minus(ibasis)-
     2                  dxi_R_minus(ibasis)*xi_R_plus(ibasis)
          write(6,*) 'Flux operator for long range function ',ibasis
          write(6,2020)
     1      'L->R',(dxi_R_plus(ibasis)*conjg(xi_R_plus(ibasis))-
     2              conjg(dxi_R_plus(ibasis))*xi_R_plus(ibasis))
     3             /(2d0*iu),
     4      'L<-R',(dxi_R_minus(ibasis)*conjg(xi_R_minus(ibasis))-
     5              conjg(dxi_R_minus(ibasis))*xi_R_minus(ibasis))
     6             /(2d0*iu)
          endif
        enddo
 2010   format(/1x,a,i3,' = (',g10.3,',',g10.3,')')
 2020   format(1x,2(a,1x,'= (',g10.3,',',g10.3,')',2x))
        write(6,*)
        deallocate (cento)

c  short-range (f^i) boundary conditions
        nopsr=0
        do ibasis=1,n
          wval=evals(ibasis)
          k=sqrt(abs(wval))
          if (wval.le.0) then
          sqrtk=sqrt(k)
c  open channels
            nopsr=nopsr+1
            if (nopsr.ne.ibasis) write(6,*) 'warning ',ibasis,nopsr,
     1                                      icalc
            xi_L_plus(ibasis) = 1.d0/sqrtk
            dxi_L_plus(ibasis)=   iu*sqrtk
          else
c  closed channels
            xi_L_plus(ibasis) = 1d0/sqrtk
            dxi_L_plus(ibasis)=  iu*sqrtk
          endif
          xi_L_minus(ibasis) =conjg(xi_L_plus(ibasis))
          dxi_L_minus(ibasis)=conjg(dxi_L_plus(ibasis))
          if (printcheck) then
          write(6,2010) 'Wronskian for short range function ',ibasis,
     1                  dxi_L_plus(ibasis)*xi_L_minus(ibasis)-
     2                  dxi_L_minus(ibasis)*xi_L_plus(ibasis)
          write(6,*) 'Flux operator for short range function ',ibasis
          write(6,2020)
     1       'L->R',(dxi_L_plus(ibasis)*conjg(xi_L_plus(ibasis))-
     2              conjg(dxi_L_plus(ibasis))*xi_L_plus(ibasis))
     3             /(2d0*iu),
     4      'L<-R',(dxi_L_minus(ibasis)*conjg(xi_L_minus(ibasis))-
     5              conjg(dxi_L_minus(ibasis))*xi_L_minus(ibasis))
     6             /(2d0*iu)
          endif
        enddo
        noi(icalc)=nopsr
        if (nopsr.gt.mxopen) then
          write(6,*)'increase mxopen! nopsr =',nopsr
          stop
        endif

c  code for indirect method
c
        ntwo=2*n
c  z2mat1 and z2mat2 are the two matrices in equation 44 of my notes -
c  z2mat1 is (eventually) [(FG')(O)-(O')]^{-1}
c  z2mat2 is [(FG')(I)-(I')]
        allocate (z2mat1(ntwo,ntwo),z2mat2(ntwo,ntwo))
        z2mat1=0.d0
        z2mat2=0.d0

        do ibasis=1,n
          do jbasis=1,n
            z2mat1(ibasis,jbasis)=fpn(ibasis,jbasis)*xi_R_plus(jbasis)
            z2mat2(ibasis,jbasis)=fpn(ibasis,jbasis)*xi_R_minus(jbasis)
          enddo
          do jbasis=1,n
            z2mat1(n+jbasis,ibasis)=fp0(jbasis,ibasis)*xi_R_plus(ibasis)
            z2mat2(n+jbasis,ibasis)=fp0(jbasis,ibasis)*
     1                                                xi_R_minus(ibasis)
            z2mat1(ibasis,n+jbasis)=gpn(ibasis,jbasis)*
     1                                                xi_L_minus(jbasis)
            z2mat2(ibasis,n+jbasis)=gpn(ibasis,jbasis)*xi_L_plus(jbasis)
          enddo
          z2mat1(ibasis,ibasis)=z2mat1(ibasis,ibasis)-dxi_R_plus(ibasis)
          z2mat2(ibasis,ibasis)=z2mat2(ibasis,ibasis)-
     1                                               dxi_R_minus(ibasis)
        enddo
        do ibasis=1,n
          ioff=n+ibasis
          do jbasis=1,n
            joff=n+jbasis
            z2mat1(ioff,joff)=gp0(ibasis,jbasis)*xi_L_minus(jbasis)
            z2mat2(ioff,joff)=gp0(ibasis,jbasis)*xi_L_plus(jbasis)
          enddo
          z2mat1(ioff,ioff)=z2mat1(ioff,ioff)-dxi_L_minus(ibasis)
          z2mat2(ioff,ioff)=z2mat2(ioff,ioff)-dxi_L_plus(ibasis)
        enddo
        call zinv(z2mat1,z2mat1,ntwo)
c
c  s2mat is -(rt) using the notation in my notes
        allocate (s2mat(ntwo,ntwo))
        s2mat=matmul(z2mat1,z2mat2)
        if (printS2mat) then
           write(*,*)"Printing complete S-matrix (real, imaginary):"
           do ibasis=1,ntwo
           do jbasis=1,ntwo
              write(*,*)jbasis,ibasis,dble(s2mat(jbasis,ibasis)),
     &        dimag(s2mat(jbasis,ibasis))
           enddo
           enddo
           write(*,*)"End printing complete S-matrix"
           write(*,*)
        endif
        deallocate (z2mat1,z2mat2)

        do jbasis=1,n
          joff=n+jbasis
        do ibasis=1,n
          ioff=n+ibasis
          r_io(ibasis,jbasis,icalc)=-s2mat(ioff,joff)
        enddo
        enddo
        do jbasis=1,noplr
        do ibasis=1,noplr
          r_oi(ibasis,jbasis,icalc)=-s2mat(ibasis,jbasis)
        enddo
        enddo
        do jbasis=1,n
          joff=n+jbasis
        do ibasis=1,noplr
          t_io(ibasis,jbasis,icalc)=-s2mat(ibasis,joff)
          t_oi(jbasis,ibasis,icalc)=-s2mat(joff,ibasis)
        enddo
        enddo
c  write out and check elements of the matrices,
c  ordered so that pairs of adjacent values should be the same (but problem
c  with sign of fp0/gpn
        if (printcheck) then
        write(6,*)
        write(6,2006) 'roi(i1,i2)','roi(i2,i1)',
     1                'rio(i1,i2)','rio(i2,i1)'
 2006   format(4x,'i1   i2',4(6x,a:,9x))
        do ibasis=1,n
        do jbasis=1,ibasis
        write(6,2005) ibasis,jbasis,s2mat(ibasis,jbasis),
     1                s2mat(jbasis,ibasis),
     2                s2mat(n+ibasis,n+jbasis),
     3                s2mat(n+jbasis,n+ibasis)
        enddo
        enddo
        write(6,*)
        write(6,2006) 'toi(i1,i2)','tio(i2,i1)',
     1                'tio(i1,i2)','toi(i2,i1)'
        do ibasis=1,n
        do jbasis=1,ibasis
        write(6,2005)ibasis,jbasis,s2mat(ibasis,n+jbasis),
     1               s2mat(n+jbasis,ibasis),
     2               s2mat(n+ibasis,jbasis),
     3               s2mat(jbasis,n+ibasis)
        enddo
        enddo
 2005   format(2i5,2x,1p,4(2x,'(',g10.3,',',g10.3,')'))
        write(6,*)
c  this check may not work if noplr/=n, in which case it should be replaced
c  by a more careful calculation, cutting out the long range closed channels
        write(6,*)'n,noplr,nopsr',n,noplr,nopsr
        write(6,*)
        write(6,*)'Testing unitarity of S matrix'
        write(6,2030) matmul(s2mat(1:n+n,1:n+n),
     1                       conjg(s2mat(1:n+n,1:n+n)))
 2030   format(4(1x,'(',g10.3,',',g10.3,')':,2x))
        endif
        deallocate (s2mat)
c  end of code for indirect method

c  start of code for testing via direct method
        if (Sdirect .and. printS) then
          nopsr=noi(icalc)
          noplr=noo(icalc)
          allocate (temp1(n,n),temp2(n,n),
     1              temp3(noplr,n),temp4(noplr,noplr),
     2              temp5(noplr,noplr),y0(n,n),yn(noplr,noplr),
     3              S(noplr,noplr),Sc(n,n))
          do iy=1,ny
            y=ymin+dble(iy-1)*dy
            do idelta=1,ndelta
              deltas=deltamin+dble(idelta-1)*ddelta+offset
              Sc=0.d0
              do iopen=1,nopsr
                Sc(iopen,iopen)=((1.d0-y)/(1.d0+y))*
     1                          exp(2.d0*iu*pi*deltas)
              enddo
              do iclosed=nopsr+1,nbasis(icalc)
                Sc(iclosed,iclosed)=((1.d0-y)/(1.d0+y))
     1                              *exp(2.d0*iu*pi*deltas)
              enddo
c  construct short range Y matrix (equation 57 in my notes)
c  note that f^{i+} and f^{i-} are xi_L_minus and xi_L_plus respectively
              temp1=0.d0
              temp2=0.d0
              do ichan=1,n
                temp1(ichan,ichan)=dxi_L_minus(ichan)
                temp2(ichan,ichan)= xi_L_minus(ichan)
              do jchan=1,n
                temp1(ichan,jchan)=temp1(ichan,jchan)
     1                             +dxi_L_plus(ichan)*Sc(ichan,jchan)
                temp2(ichan,jchan)=temp2(ichan,jchan)
     1                             + xi_L_plus(ichan)*Sc(ichan,jchan)
              enddo
              enddo
              call zinv(temp2,temp2,n)
              y0=matmul(temp1,temp2)
c  construct long range Y matrix (equation 60 in my notes)
              y0=y0-gp0(1:n,1:n)
              call zinv(y0,y0,n)
              temp3=matmul(gpn(1:noplr,1:n),y0)
              yn=fpn+matmul(temp3,fp0(1:n,1:noplr))
c  obtain S matrix from long range boundary conditions and Y matrix (eqn 64)
              temp4=0.d0
              temp5=0.d0
              do iopen=1,noplr
                temp4(iopen,iopen)=-dxi_R_plus(iopen)
                temp5(iopen,iopen)=-dxi_R_minus(iopen)
              do jopen=1,noplr
                temp4(iopen,jopen)=temp4(iopen,jopen)
     1                             +yn(iopen,jopen)*xi_R_plus(jopen)
                temp5(iopen,jopen)=temp5(iopen,jopen)
     1                             +yn(iopen,jopen)*xi_R_minus(jopen)
              enddo
              enddo
              call zinv(temp4,temp4,noplr)
              S=-matmul(temp4,temp5)

c  tabulate the S matrix for the direct method
              write(6,*)
              write(*,*)'direct for y =',y,', delta =',deltas
              write(*,100)'  ROW  COL'
              do iopen=1,noplr
                do jopen=1,noplr
                  sreal=real(S(iopen,jopen))
                  simag=aimag(S(iopen,jopen))
                  ssquare=sreal**2+simag**2
                  if (ssquare.gt.1d0) then
c                   sreal=sreal/sqrt(ssquare)
c                   simag=simag/sqrt(ssquare)
c                   ssquare=1.d0
                  endif
                  sphase=atan2(simag,sreal)/(2d0*pi)
                  if (ssquare.gt.1d-15) 
     1              write(6,110)iopen,jopen,ssquare,sphase,sreal,simag
                enddo
              enddo
            enddo
          enddo
          deallocate (temp1,temp2,temp3,temp4,temp5,y0,yn,S,Sc)
        endif
            
        deallocate (fp0,fpn,gp0,gpn)
        deallocate (xi_R_plus,xi_R_minus,dxi_R_plus,dxi_R_minus)
        deallocate (xi_L_plus,xi_L_minus,dxi_L_plus,dxi_L_minus)
        deallocate (evals)
c  end of code for testing via direct method
      enddo
      deallocate (levs)

c  work out which pair states are included for each energy, and tabulate and
c  order them in ielev
      ncalc=icalc-1
      allocate (ielev(mxopen,nnrg),mxlev(nnrg))
      mxlev=0
c  work through each of the calculations, working out, for each open pair
c  states, which pair level it correlates with
      do icalc=1,ncalc
        inrg=iecalc(icalc)
        iopen_loop: do iopen=1,noo(icalc)
          idlev=id(iopen,icalc)
          do ilev=1,mxlev(inrg)
            if (idlev.eq.ielev(ilev,inrg)) then
              cycle iopen_loop
            endif
          enddo
          mxlev(inrg)=mxlev(inrg)+1
          ielev(mxlev(inrg),inrg)=idlev
        enddo iopen_loop
      enddo
c  then sort them
      do inrg=1,nnrg
        do ilev=1,mxlev(inrg)-1
          do jlev=ilev+1,mxlev(inrg)
            if (ielev(ilev,inrg).gt.ielev(jlev,inrg)) then
              idlev=ielev(ilev,inrg)
              ielev(ilev,inrg)=ielev(jlev,inrg)
              ielev(jlev,inrg)=idlev
            endif
          enddo
        enddo
      enddo

      allocate (sigma(nnrg,nlevel,nlevel))
      allocate (sigpar(nnrg,ncalc,maxval(noo),maxval(noo)))
      allocate (slpar(nnrg,ncalc,maxval(noo)))
      allocate (ssuminel(maxval(noo)))
      allocate (sigcur(nlevel,nlevel))
      allocate (sloss(nnrg,nlevel))
      allocate (slcur(nlevel))
      sigma=0.d0
      sloss=0.d0
      sigpar=0.d0
      slpar=0.d0
      write(6,*)
      do iy=1,ny
        y=ymin+(iy-1)*dy
        do idelta=1,ndelta
          deltas=deltamin +(idelta-1)*ddelta+offset
          write(*,*)'indirect for y =',y,'deltas = ',deltas
          sigma=0.d0
          sloss=0.d0
          sigpar=0.d0
          slpar=0.d0
c  cycle through each calculation, tabulating the S matrix and calculating
c  each contribution to the cross section, also accumulating the elastic,
c  inelastic and short range loss(?)
          do icalc=1,ncalc
            inrg=iecalc(icalc)
            if (Jmax.ge.0 .and. Jtotal(icalc).gt.Jmax) cycle
            xj=pi*wght(icalc)
            if (jhalf.eq.1) then
              xj=xj*dble(2*Jtotal(icalc)+1)
            elseif (jhalf.eq.2) then
              xj=xj*dble(Jtotal(icalc)+1)
            endif
            noplr=noo(icalc)
            nopsr=noi(icalc)
            n=nbasis(icalc)

            allocate (Sc(n,n),zmat1(n,n))
            allocate (zmat2(noplr,n),zmat3(noplr,noplr))
            if (nopsr.eq.0) write (6,*)'no open channels at short range'
            Sc=0.d0
            do iopen=1,nopsr
              Sc(iopen,iopen)=((1.d0-y)/(1.d0+y))*exp(2.d0*iu*pi*deltas)
            enddo
            do iclosed=nopsr+1,n
              Sc(iclosed,iclosed)=((1.d0-y)/(1.d0+y))*
     1                            exp(2.d0*iu*pi*deltas)
            enddo

            zmat1=-matmul(r_io(1:n,1:n,icalc),Sc)
            do ichan=1,n
              zmat1(ichan,ichan)=zmat1(ichan,ichan)+1.d0 !Scinv
            enddo
            call zinv(zmat1,zmat1,n)
            zmat1=matmul(Sc,zmat1)
            zmat2=matmul(t_io(1:noplr,1:n,icalc),zmat1)
            zmat3=matmul(zmat2,t_oi(1:n,1:noplr,icalc))
            zmat3=zmat3+r_oi(1:noplr,1:noplr,icalc)
            sigcur=0.d0
            slcur=0.d0
            if (printa) then
              write(6,605) (trim(runame),kdum=1,3)
  605 FORMAT(/'  K-DEPENDENT SCATTERING LENGTHS/VOLUMES/HYPERVOLUMES',
     1       ' FOR CHANNELS WITH LOW KINETIC ENERGY'/
     2       '  CHAN   L POW',5X,'WVEC*',A10,7X,'RE(A)/',A10,6X,
     3       'IM(A)/',A)
              do iopen=1,noo(icalc)
                !dd=wv(iopen,icalc)**2/runit**2
                dd=wv(iopen,icalc)
                if (dd.gt.awvmax) cycle
!               scatlen=iu*dd**min(lrpow-2,2*ll(iopen,icalc)+1)*
!    1                 (1d0-zmat3(iopen,iopen))/(1d0+zmat3(iopen,iopen))
                scatlen=(1d0-zmat3(iopen,iopen))
     1                 /((1d0+zmat3(iopen,iopen))*
     2                  iu*dd**min(lrpow-2,2*ll(iopen,icalc)+1))
                WRITE(6,604) iopen,LL(iopen,icalc),
     1                     MIN(LRPOW-2,2*LL(iopen,icalc)+1),
     2                     DD,DBLE(SCATLeN),IMAG(SCATLeN)
  604           FORMAT(2I5,I4,4(ES21.13E3,2X))
              enddo
            endif
            if (Sindirect .and. printS)
     1        write(6,*)'S matrix for J=',Jtotal(icalc),'symmetry ',
     2                  'block',iblock(icalc),' at energy = ',
     3                  energy(inrg)
            if (Sindirect .and. printS) write(6,100)'  ROW  COL'
            do iopen=1,noo(icalc)
              ifrom=id(iopen,icalc)
              dd=wv(iopen,icalc)**2/runit**2
              dg=dble(degin(ifrom))
              xjfac=xj/(dg*dd)
              stemp=1.d0
              do jopen=1,noo(icalc)
                ito=id(jopen,icalc)
                simag=aimag(zmat3(jopen,iopen))
                sreal=real(zmat3(jopen,iopen))
                ssquare=sreal**2+simag**2
                if (ssquare.gt.1d0) then
                  sreal=sreal/sqrt(ssquare)
                  simag=simag/sqrt(ssquare)
                  ssquare=1.d0
                endif
                sphase=atan2(simag,sreal)/(2d0*pi)
c  tabulate the S matrix
                if (ssquare.gt.1d-15 .and. Sindirect .and. printS) then
                  write(6,110) iopen,jopen,ssquare,sphase,sreal,simag
                endif
                if (ito.eq.ifrom) stemp=stemp-ssquare
                if (jopen.eq.iopen) ssquare=(1d0-sreal)**2+simag**2
c  calculate (and accumulate) the current contribution to the cross section
                sigcur(ifrom,ito)=sigcur(ifrom,ito)
     1                            +ssquare*xjfac*degout(ito)
                sigma(inrg,ifrom,ito)=sigma(inrg,ifrom,ito)
     1                                +ssquare*xjfac*degout(ito)
c  bm: get partial cross sections                
                sigpar(inrg,icalc,iopen,jopen)=ssquare*xj/dd
              enddo
c  calculate (and accumulate) the loss from each pair state
              slcur(ifrom)=slcur(ifrom)+stemp*xjfac
              sloss(inrg,ifrom)=sloss(inrg,ifrom)+stemp*xjfac
c  bm: get partial loss cross sections                
              slpar(inrg,icalc,iopen)=stemp*xjfac
            enddo
            nlev=nlevs(icalc)
            write(6,*)'current cross section contributions'
            write(6,*)'   Jtot  iblock  energy  from -> to  '//
     1                '              y       delta'
            write(6,120)((Jtotal(icalc),iblock(icalc),energy(inrg),
     1                    ielev(ilev,inrg),ielev(jlev,inrg),
     2                    sigcur(ielev(ilev,inrg),
     3                           ielev(jlev,inrg)),
     4                    y,deltas, ilev=1,nlev), jlev=1,nlev)
            write(6,*)'current loss contributions'
            write(6,*)'     Jtot  iblock  energy    from    '//
     1                '                y       delta'
            write(6,130) (Jtotal(icalc),iblock(icalc),energy(inrg),
     1                    ielev(ilev,inrg),
     2                    slcur(ielev(ilev,inrg)),
     3                    y,deltas, ilev=1,nlev)
            deallocate (Sc,zmat1,zmat2,zmat3)
          enddo
          write(6,*)'-----------------------------------------'

          nopen = maxval(noo)
          do inrg=1,nnrg
            nlev=mxlev(inrg)
c  bm:
c  sum partial cross sections over outgoing L partial waves for given
c  incoming partial wave and print out            
            if(printsigpar) then
              do icalc=1,ncalc
                write(6,*)'accumulated partial cross sections for',
     1                    ' icalc =',icalc         
                write(6,*)'nlevel  iopen  jopen  sigpar  y  dels',
     1                    ' energy'
                do iopen=1,nopen
                  if(.not.(id(iopen,icalc).eq.ilevel)) cycle
                  write(6,180)(ilevel,iopen,jopen,
     1                         sigpar(inrg,icalc,iopen,jopen),
     1                         y,deltas,energy(inrg),' acc xs',
     2                         jopen=1,nopen)
                enddo
                write(*,*)"Diagonal elastic contributions",
     1                    " for a given incoming channel:"
                do iopen=1,nopen
                   if(.not.(id(iopen,icalc).eq.ilevel)) cycle
                   write(*,'(i5,g18.8)')iopen,
     1                   sigpar(inrg,icalc,iopen,iopen)
                enddo
                write(*,*)"Elastic contributions summed over",
     1                    " outgoing channels for a given incoming",
     2                    " channel:"
                do iopen=1,nopen
                   if(.not.(id(iopen,icalc).eq.ilevel)) cycle
                   ssum=0.d0
                   do jopen=1,nopen
                     if(.not.(id(jopen,icalc).eq.ilevel)) cycle
                     ssum=ssum+sigpar(inrg,icalc,iopen,jopen)
                   enddo
                   write(*,'(i5,g18.8)')iopen,ssum
                enddo
                write(*,*)"Inelastic contributions summed over",
     1                    " outgoing channels for a given incoming",
     2                    " channel:"
                do iopen=1,nopen
                   if(.not.(id(iopen,icalc).eq.ilevel)) cycle
                   ssum=0.d0
                   do jopen=1,nopen
                     if(id(jopen,icalc).eq.ilevel) cycle
                     ssum=ssum+sigpar(inrg,icalc,iopen,jopen)
                   enddo
                   ssuminel(iopen)=ssum
                   write(*,'(i5,g18.8)')iopen,ssum
                enddo
                write(*,*)"Total loss contributions summed over",
     1                    " outgoing channels for a given incoming",
     2                    " channel:"
                do iopen=1,nopen
                   if(.not.(id(iopen,icalc).eq.ilevel)) cycle
                   write(*,'(i5,g18.8)')iopen,slpar(inrg,icalc,iopen)
                enddo
                write(*,*)"SR loss contributions summed over",
     1          " outgoing channels for a given incoming",
     2           " channel:"
                do iopen=1,nopen
                   if(.not.(id(iopen,icalc).eq.ilevel)) cycle
                   write(*,'(i5,g18.8)')iopen,
     1             (slpar(inrg,icalc,iopen)-ssuminel(iopen))
                enddo


              enddo ! icalc loop
            endif   ! printsigpar if block
c  bm: partial cross sections printing out ends here

            write(6,*) 'accumulated cross sections/Angstrom^2' !trim(runame)
            write(6,150) ((ielev(ifrom,inrg),ielev(ito,inrg),
     1                     sigma(inrg,ielev(ifrom,inrg),
     2                                ielev(ito,inrg)),
     3                     y,deltas,energy(inrg),' acc xs',
     4                     ifrom=1,nlev),ito=1,nlev)
            allocate (sum_inelastic(mxlev(inrg)),ifindx(mxopen))
            sum_inelastic=0.d0
            do ilev=1,nlev
              ift=ielev(ilev,inrg)
              ifindx(ift)=ilev
            enddo
            do ilev=1,nlev
              ifrom=ielev(ilev,inrg)
              ifdex=ifindx(ifrom)
              do jlev=1,nlev
                ito=ielev(jlev,inrg)
                itdex=ifindx(ito)
                if (ito.eq.ifrom) cycle
                sum_inelastic(ifdex)=sum_inelastic(ifdex)
     1                               +sigma(inrg,ifrom,ito)
              enddo
            enddo
            write(6,*) 'accumulated total inelastic/Angstrom^2' !trim(runame)
            write(6,160) (ielev(ilev,inrg),
     1                    sum_inelastic(ifindx(ielev(ilev,inrg))),
     2                    y,deltas,energy(inrg),'acc inelastic',
     3                                       ilev=1,nlev)
            deallocate (sum_inelastic,ifindx)
            write(6,*)
            write(6,*) 'accumulated losses/Angstrom^2' !trim(runame)
            write(6,170) (ielev(ifrom,inrg),
     1                    sloss(inrg,ielev(ifrom,inrg)),
     2                    y,deltas,energy(inrg),' acc loss',
     3                    ifrom=1,nlev)
          enddo
          write(6,*)'-----------------------------------------'
        enddo
      enddo
  100 format(/a10,7x,'S**2',18x,'PHASE/2PI',14x,'RE (S)',17x,'IM (S)'  )
  110 format(2i5,4(es21.13e3,2X):,f6.2,2x,f10.3:,2i5)
  120 format((1x,'xs',2i5,2x,0pf10.6,2i5,1p,g13.6,0p,2(2x,f6.2)))
  130 format((1x,'loss',2i5,2x,0pf10.6,i5,5x,1p,g13.6,0p,2(2x,f6.2)))
  150 format((0p,2i5,1x,1p,g15.8,2(2x,g16.9),2x,g16.9,1x,a))
  160 format((0p, i5,1x,1p,g15.8,2(2x,g16.9),2x,g16.9,1x,a))
  170 format((0p, i5,1x,1p,g15.8,2(2x,g16.9),2x,g16.9,1x,a))
  180 format((0p,3i5,1x,1p,g15.8,2(2x,g16.9),2x,g16.9,1x,a))    

      deallocate (Jtotal,wv,degin,degout,id,iblock,iecalc)
      deallocate (r_io,r_oi,t_io,t_oi)
      deallocate (sigma,energy,nlevs)
      deallocate (sigcur,sloss)

      end
