          module bigstuff
          implicit real*8(a-h,o-z)
          real*8,allocatable :: CF(:,:,:)
          end module bigstuff
          module angouter
          real*8,allocatable :: AC(:,:),BVC(:,:),BLC(:,:)
          end module angouter
          module workham
          real*8,allocatable :: wrkcb(:,:,:)
          real*8,allocatable :: wrkbb(:)
          end module workham
          module workcont
          real*8,allocatable :: wrkcc(:,:,:)
          end module workcont
          module workdipole
          real*8,allocatable :: wrkd(:,:,:)
          end module workdipole
          module datamod
          implicit real*8(a-h,o-z)
          integer :: MJS=1,jfile=32,kfile=33,iprint=1,ICODE=25,IZESP=0
          integer :: JRELOP=1,maxnc=0,isat=0,NSPN=0
          logical :: outer=.true.
          real*8 :: ABUTL=0.0d0,ABUTV=0.0d0,BBUTL=0.0d0,BBUTV=0.0d0
          real*8 :: eta=0.0d0
          end module datamod
C CPB                                                        14/01/06
C NRB                          PROGRAM DTO3                  14/01/06
C
c based-on kab's d2h of Dec 2002.
c
c  Converts DARC to RMBP file format: DSTG2.DAT to RECUPH.DAT and RECUPD.DAT.

c  ifile = unit number for DSTG2.DAT input
c  INAST = number of N+1 electron symmetries on file: needs to be defined
C          but can be set large, e.g. 999, to obtain all (stg3 moans a bit) NRB.
c  nproc = number of RECUPH.DAT files to be written (=1 for serial codes).

      logical ok
      data ifile/30/, INAST/0/

      print*,'Enter number of N+1 electron symmetries on DSTG2.DAT and'
      print*,'number of RECUPH.DAT files to write (=1 for serial codes)'
      read(5,*) INAST,nproc
C      INAST=1
C      nproc=1
      if( INAST.le.0 .or. nproc.le.0 ) stop

c Open dstg2.f file for reading:

      INQUIRE (FILE='DSTG2.DAT',EXIST=ok)
      if (.not.ok) then
        print*,'STOP: file DSTG2.DAT does not exist'
        stop
      endif
      OPEN(UNIT=ifile,STATUS='OLD',FORM='UNFORMATTED',FILE='DSTG2.DAT')

      CALL DTO3(INAST,nproc,ifile)
      stop
      end

c-----------------------------------------------------------------------
      SUBROUTINE DTO3(INAST,nproc,IFILE)
      use bigstuff
      use angouter
      use workham
      use workcont
      use workdipole
      use datamod
c  Converts DARC to RMBP file format: DSTG2.DAT to RECUPH.DAT and RECUPD.DAT.

c  The following need to be predefined:
c  INAST = number of N+1 electron symmetries on file
c  ifile = unit number for DSTG2.DAT input
c  nproc = number of RECUPH.DAT files to be written (=1 for serial codes).

      IMPLICIT REAL*8(A-H,O-Z)

c  DSTG2.DAT is the file from the DARC program dstg2.f, 
c  and can contain both Hamiltonian and dipole matrix elements.

c  The DTO3 routine reads this file and creates file RECUPH.DAT with
c  the Hamiltonians, and RECUPD.DAT with the dipoles if present,
c  in the format required by the BPRM program stg3r.f
c  (the format was defined by Berrington et al 1995 in CPC 92,290-420).

c  It also converts the Buttle correction to 'MJS' style.
c  It also detects any external-region contributions to the dipole matrix
c  (these can now be calculated in dstg2.f with some mods by kab).

c   DARC style Buttle correction written to DBUT.DAT (case MJS=0) for
c   automatic read by stgf.

c Include the dimension parameters MZLR2,MZNR2,MZNC2,MZLMX,MZCHF,MZSLP,MZTAR
c using the Badnell (stg3r.f etc) convention:

      INCLUDE 'PARAM'

c This could be done using the DARC dstg2.f dimension plants as follows:
c      PARAMETER (MZLR2={NK})
c      PARAMETER (MZNR2={NB})
c      PARAMETER (NXNC2={NC})
c      PARAMETER (MZCHF={CH})
c      PARAMETER (MZTAR={NL})
c      PARAMETER (MZLMX=4)
c      PARAMETER (MZSLP=100)

c Derive further parameters to set dimensions.
c The biggest array here is work, whose size=MZNR2*MZNC2*MZCHF*MZCHF/2.

      PARAMETER (INML=MZNR2/MZNC2, JNML=MZNC2/MZNR2, KNML=INML+JNML)
      PARAMETER (ND24=MZNR2*INML/KNML+MZNC2*JNML/KNML)
c      PARAMETER (kccc=(MZCHF*(MZCHF+1))/2)
      PARAMETER (NBUT=2*(MZNR2-1))
      parameter (zero=0.0)
c
      REAL*8,allocatable :: ARRAY(:),workcc(:,:)
c
      CHARACTER IHED*80,RECORD*20
      CHARACTER*1 NUM(0:9)

      logical bnew

      DIMENSION MAXNQN(MZLR2), MINNQN(MZLR2),
     *          MAXFUL(MZLR2), LAG(MZLR2), KAG(MZLR2), JAG(MZLR2),
     *          NCONAT(MZTAR), K2P(MZCHF), L2P(MZCHF),
     *          MBUT(MZLR2), ksym(MZCHF,MZCHF),
     *          ILSPI(3,MZSLP), temp(MZLR2)
      DIMENSION EBUT(NBUT,MZLR2),CBUT(NBUT,MZLR2),
     *          e(NBUT), f(NBUT), eigmax(MZLR2)
C     *          wrkcb(MZNR2,MZNC2,MZCHF),wrkbb(MZNC2)

      COMMON/COM1/ENDS(MZNR2,MZLR2),ENAT(MZTAR),
     *            EIGENS(MZNR2,MZLR2),COEFF(3,MZLR2)
      COMMON/COM2/JAT(MZTAR),LPTY(MZTAR), MAXNLG(MZLR2)
      COMMON/COM3/NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,LAM,
     *            NAST,NIX,RA,BSTO,HINT

c  MJS   = 1 (default) for producing Buttle fit;
c        = 0 for dump on DBUT.DAT (this also occurs if fit fails).
c  outer = true (default) for outer region contributions, if any;
c        = false for inserting zeros.
c  jfile = 32 (default) unit number for RECUPH.DAT output.
c  kfile = 33 (default) unit number for RECUPD.DAT output.
c  iprint= 1 (default) for printing useful stuff on standard output;
c        = 0 for emergency printing only.
c  The other parameters set here should be acceptable for BPRM.

      DATA NUM/'0','1','2','3','4','5','6','7','8','9'/

c  Check dimensions for number of symmetries input

      if( INAST.gt.MZSLP ) then
        print*,'STOP! increase dimension MZLSP to:',MZSLP
        stop
      endif

c
c nrec=number of symmetries per file
c
      nrec=inast/nproc
      if(nrec.le.0)nrec=1

c  Open Hamiltonian output file and info file
      if(nproc.eq.1)then
        OPEN( UNIT=jfile, STATUS='UNKNOWN', FORM='UNFORMATTED',
     *      FILE='RECUPH.DAT')
      endif
      OPEN(99,file='sizeBP.dat',status='unknown',form='formatted')

c-----------------------------------------------------------------------
c Read header from DSTG2.DAT file

c Some of the variables are DARC-specific and do not seem to correspond
c to anything in BPRM, eg speed of light CL, so we will try to ignore them.

c On the other hand, there are some variables that BPRM would like,
c which the DARC file does not have, eg. LAMAX, so we will try to derive them.

c Along the way, we will also check dimensions.

      rewind ifile
      READ (ifile) IHED,RECORD
      if(iprint.ge.1)print*,  IHED,RECORD
      READ (ifile) KBMAX,KCMIN,KCMAX,KX,NRANG2,NELC,NZ
      READ (ifile) LAMBB,LAMBC,LAMCC

      LRANG1 = KX
      LRANG2 = KCMAX                ! but starts at KCMIN
      LAMAX  = LAMBB
      LAM    = 1
      if( LAMBC.gt.0 ) LAM = 2
      if( LAMCC.gt.0 ) LAM = 3

      if( LRANG1.gt.MZLR2 .or. LRANG2.gt.MZLR2  ) then
        print*,'STOP! Increase dimension MZLR2 to:',max(LRANG1,LRANG2)
        stop
      endif
      if( LAMAX.gt.MZLMX ) then
        print*,'STOP! Increase dimension MZLMX to:',LAMAX
        stop
      endif
      if( NRANG2.gt.MZNR2 ) then
        print*,'STOP! Increase dimension MZNR2 to:',NRANG2
        stop
      endif

      IF( KBMAX.GT.0 ) THEN
        READ (ifile) (MINNQN(L),L=1,KX), (MAXNQN(L),L=1,KX),
     *               (MAXNLG(L),L=1,KX), (MAXFUL(L),L=1,KBMAX),
     *               (LAG(L),L=1,KX), (KAG(L),L=1,KX), (JAG(L),L=1,KX)
      ELSE
        READ (ifile) (MINNQN(L),L=1,KX), (MAXNQN(L),L=1,KX),
     *               (MAXNLG(L),L=1,KX), (MAXFUL(L),L=1,KX),
     *               (LAG(L),L=1,KX),(KAG(L),L=1,KX), (JAG(L),L=1,KX)
      ENDIF

      READ (ifile) RA,BSTO,HINT,CL,DR1
      READ (ifile) NIX
      READ (ifile) (IRX,I=1,NIX)
      READ (ifile) NUMORB

c Write header to standard RECUPH.DAT hamilitonian output file

      if(nproc.eq.1)then
        WRITE(jfile) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM,IZESP,
     *               JRELOP,JRELOP,JRELOP
        WRITE(jfile) (MAXNLG(L),L=1,LRANG1), (MAXNLG(L),L=1,LRANG1),
     *               (maxnc,L=1,LRANG1)
      endif

      IF (KCMIN.GT.1) THEN
        DO L = 1,KCMIN-1
          DO N=1,NRANG2
            ENDS(N,L)=ZERO
            EIGENS(N,L)=ZERO
          ENDDO
          if(nproc.eq.1)then
            WRITE(jfile) (EIGENS(N,L),N=1,NRANG2)
            WRITE(jfile) (ENDS(N,L),N=1,NRANG2)
          endif
        ENDDO
      ENDIF
c Finish reading header

      DO L = KCMIN,KCMAX
        READ (ifile) (EIGENS(N,L),N=1,NRANG2)
        if(nproc.eq.1)WRITE(jfile) (2*EIGENS(N,L),N=1,NRANG2) !a.u. to Ryd
        eigmax(L) = EIGENS(NRANG2,L)
        READ (ifile) (ENDS(N,L),N=1,NRANG2)
        if(nproc.eq.1)WRITE(jfile) (ENDS(N,L),N=1,NRANG2)
        READ (ifile) MBUT(L)
        READ (ifile) (EBUT(I,L),I=1,MBUT(L))
        READ (ifile) (CBUT(I,L),I=1,MBUT(L))
      ENDDO

      READ (ifile) NAST
      if(NAST.gt.MZTAR)then
        print *,'STOP! Increase dimension MZTAR to:',NAST
        stop
      endif
      READ (ifile) (JAT(I),I=1,NAST)
      READ (ifile) (ENAT(I),I=1,NAST)

c Print target states

      if( iprint.ge.1 ) then
        print*,'NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,RA,NAST:',
     *          NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,RA,NAST
        print*,'Target states:'
        do i = 1,nast
          LPTY(i) = 0
          if( JAT(I).lt.0 ) LPTY(i) = 1
          JAT(i) = abs(JAT(i))-1
          aj = 0.5*JAT(i)
          if( LPTY(i).eq.1 ) aj = -aj
          ee = 2.0*(enat(i)-enat(1))
          if( i.eq.1 ) write(6,1001) aj,enat(1)
          if( i.gt.1 ) write(6,1002) aj,ee,ee*13.6057
        enddo
 1001   format(f5.1,'    0.0',1pe17.9,' au')
 1002   format(f5.1,2f12.6,' eV')
      endif

C  FIT BUTTLE CORRECTION AND STORE FIT PARAMETERS IN  COEFF(I,L),I=1,3

      IF (KCMIN.GT.1) THEN
        DO L = 1,KCMIN-1
          COEFF(1,L) = zero
          COEFF(2,L) = zero
          COEFF(3,L) = zero
        ENDDO
      ENDIF

      IF (BSTO.NE.ZERO .OR. MJS.EQ.0) GOTO 90

C  IF MJS=1 and BSTO=0 USE IMPROVED FITTING PROCEDURE 
c  (M J SEATON, J.PHYS.B20(1987)L69-72) - remember to convert to Ryd!

      COEFF(3,1)=-30000.      !FLAG NEWBUT CASE KCMIN.GT.1
      DO L = KCMIN,KCMAX
        imax = MBUT(L)
        do i=1,imax
          e(i) = 2*EBUT(i,L)  !BECAUSE BUTFIT REDEFINES E
          f(i) = CBUT(i,L)    !UNNECESSARY AS F NOT REDEFINED
        enddo
        emax = 2*eigmax(L)
        CALL BUTFIT(imax,e,f,RA,emax,ALPHA,BETA,nNBUT,DELTA)
        if(iprint.gt.1)print*,'Buttle',L,nNBUT,ALPHA,BETA
        IF (nNBUT.le.0) THEN
          print*,'NEW BUTFIT DIVERGING, USE SUBROUTINE INTBUT IN STGF/B'
          goto 90
        ENDIF
        COEFF(1,L) = ALPHA
        COEFF(2,L) = BETA
        COEFF(3,L) = -10000*nNBUT
      ENDDO
      goto 120

C  IF MJS=0 OR IF BUTFIT DIVERGING, dump Buttle corrections into DBUT.DAT 
C  (STGF CANNOT HANDLE A MIXTURE OF OLD AND NEW BUTTLE FITS): keep as a.u.

   90 OPEN(31,FILE='DBUT.DAT',FORM='UNFORMATTED')
C
      IF (KCMIN.GT.1) THEN
        DO L = 1,KCMIN-1
          MBUT(L) = 1
          EBUT(1,L) = zero
          CBUT(1,L) = zero
          eigmax(L) = zero
        ENDDO
      ENDIF
C
      DO L = 1, LRANG2
        COEFF(1,L) = zero
        COEFF(2,L) = zero
        COEFF(3,L) = zero
        WRITE (31) MBUT(L)
        WRITE (31) (EBUT(N,L),N=1,MBUT(L))
        WRITE (31) (CBUT(N,L),N=1,MBUT(L))
      ENDDO
c 1031 format(1p100e20.12)
      if(iprint.ge.0)print*,
     *  'Note, Buttle correction for INTBUT routine put in DBUT.DAT'
  120 continue

c Finish writing header
      NIX=0                                       !TO SKIP READ OF WERNERS'
      if(nproc.eq.1)then
        WRITE(jfile) RA,BSTO,HINT,eta,eta,NIX
        WRITE(jfile) ((COEFF(I,L),I=1,3),L=1,LRANG2)
        WRITE(jfile) NAST
        WRITE(jfile) (ENAT(I),I=1,NAST), (JAT(I),I=1,NAST),
     *               (isat,I=1,NAST), (LPTY(I),I=1,NAST)
      endif
       

      IPOLPH = 1
      if(iprint.ge.1)print*,'N+1 electron symmetries:'
c-----------------------------------------------------------------------
c  Loop over each total symmetry:
      iam=-1
      do loop = 1, INAST

        bnew=nproc.gt.1.and.mod(loop-1,nrec).eq.0
        if(bnew)then
          if(iam.ge.0)close(jfile)
          iam=iam+1
          i1=iam/100
          i2=(iam-100*(iam/100))/10
          i3=iam-(100*(iam/100))-i2*10
          OPEN(UNIT=jfile, STATUS='UNKNOWN', FORM='UNFORMATTED',
     *    FILE='RECUPH'//NUM(i1)//NUM(i2)//NUM(i3))
          CALL HEAD(jfile)                       !write header to each file
        endif


        READ(ifile,end=99) J2P1,NPTY,NCFGP

c  The following IF structure allows automatic detection of where the
c  dipole matrices start: it is only used if you have set INAST too big.

        if( IPOLPH.eq.2 .and. J2P1.ge.0 ) then
          if(iprint.ge.1)print*, J2P1,NPTY,NCFGP
          backspace ifile        
          go to 10
        endif

        if( NCFGP.gt.MZNC2 ) then
          print*,'STOP! Increase dimension MZNC2 to:', NCFGP
          stop
        endif

        READ(ifile) MNP1,NCHAN
        if(iprint.ge.1)print*,loop,') J2P,NPTY,NCFGP,MNP1,NCHAN:',
     *                             J2P1-1,NPTY,NCFGP,MNP1,NCHAN
        if( NCHAN.gt.MZCHF ) then
          print*,'STOP! Increase dimension MZCHF to:', NCHAN
          stop
        endif

C   CPB  create  sizeBP.dat for pstg3r allocations
        write(99,*)NCHAN,MNP1-NCFGP,MNP1,
     A            ' 2J= ',abs(J2P1)-1,' PI ',NPTY

c Initialize ksym pointer to symmetric hamilitonian matrix elements

        kcc = 0
        do i = 1, NCHAN
          do j = 1,i
            kcc = kcc + 1
            ksym(i,j) = kcc
            ksym(j,i) = kcc
          enddo
        enddo

        READ(ifile) (NCONAT(I),I=1,NAST)
        READ(ifile) (K2P(I),I=1,NCHAN)          !kappa

        MORE   = 1
        if( loop.eq.INAST ) MORE = 0            !case INAST/nproc.ne.integer
        if(mod(loop,nrec).eq.0)MORE=0

        if (J2P1.lt.0) then
          IPOLPH = 2
          J2P1 = -J2P1
        endif
        do i=1,NCHAN                            !NRB
          L2P(I) = 2*K2P(I)                     !2*kappa
          if(K2P(I).lt.0) L2P(I) = -L2P(I)-1    !K, analogue of L
          L2P(I)=L2P(I)-1                       !K-1
        enddo
        if(iprint.gt.1)print*,'L2P:',(L2P(i),i=1,NCHAN)

        ILSPI(1,loop) = J2P1*NPTY
        ILSPI(2,loop) = MNP1
        ILSPI(3,loop) = NCHAN
        if( NPTY.gt.0 ) NPTY = 0
        if( NPTY.lt.0 ) NPTY = 1

        WRITE(jfile) J2P1-1,NSPN,NPTY,NCFGP,IPOLPH
        WRITE(jfile) MNP1,NRANG2*NCHAN,NCHAN
        write(6,*)MNP1,NRANG2*NCHAN,NCHAN
        WRITE(jfile) (NCONAT(I),I=1,NAST)
        WRITE(jfile) (L2P(I),I=1,NCHAN),(K2P(I),I=1,NCHAN) !kappa not used
        WRITE(jfile) MORE

c Asymptotic coeffs

        READ (ifile) IKMAX
        allocate(CF(NCHAN,NCHAN,LAMAX))
c
c
c
        IF (IKMAX.GT.0) THEN
          DO LAM = 1,IKMAX
            READ (ifile) ((CF(J,I,LAM),J=1,NCHAN),I=1,NCHAN)
          ENDDO
        ENDIF
        IF(IKMAX.LT.LAMAX)THEN
          DO LAM=IKMAX+1,LAMAX
            DO I=1,NCHAN
              DO J=I,NCHAN
                CF(J,I,LAM)=ZERO
              ENDDO
            ENDDO
          ENDDO
        ENDIF
c
c
        if( LAMAX.GT.0 ) then
      OPEN(UNIT=2233,FILE='CFSCRAT',FORM='UNFORMATTED',STATUS='UNKNOWN')
c
          do I = 1, NCHAN
       WRITE(2233)((2*CF(J,I,K),K=1,LAMAX),J=I,NCHAN)  !2* NRB
          enddo
        deallocate(CF)
        CLOSE(2233)
        endif
c
c  Hamiltonian matrix

C  Continuum-continuum blocks
C
C       The following direct access route does slow things 
C       down but scales the size of the executable by orders 
C       of magnitude   ! CPB 2006
C
C
        NRANG2SQ=NRANG2*NRANG2
        kccc=(NCHAN*(NCHAN+1))/2
        allocate(array(NRANG2SQ),stat=ierr)
        CALL DA2(2,0,199,NRANG2SQ,ARRAY)
C
        allocate(workcc(NRANG2,NRANG2))
C
C
        DO ICH = 1,NCHAN
          DO JCH = ICH,NCHAN
            READ (ifile) I,J
            kcc = ksym(I,J)
            IF (I.EQ.J) THEN
              READ (ifile) ((workcc(M,N),N=M,NRANG2),M=1,NRANG2)
c  Note, the diagonal block needs to be symmetrised
              do M=1,NRANG2
                do N=M,NRANG2
                  workcc(N,M)=workcc(M,N)
                enddo
              enddo
C
             KK=0
             DO II=1,NRANG2
              DO JJ=1,NRANG2
               KK=KK+1
               ARRAY(KK)=workcc(JJ,II)
              ENDDO
             ENDDO
C
             CALL DA2(2,kcc,199,NRANG2SQ,ARRAY)
C              print*,NRANG2SQ,ARRAY
C              print*,kcc,i,j,ich,jch,(wrkcc(M,M,kcc),M=1,NRANG2)
            ELSE
              READ (ifile) ((workcc(M,N),N=1,NRANG2),M=1,NRANG2)
C
             KK=0
             DO II=1,NRANG2
              DO JJ=1,NRANG2
               KK=KK+1
               ARRAY(KK)=workcc(JJ,II)
              ENDDO
             ENDDO
C
             CALL DA2(2,kcc,199,NRANG2SQ,ARRAY)
            ENDIF
          ENDDO
        ENDDO
c
c  Note, transpose is written out, because input is on one side of the
c  triangle (JCH = ICH,NCHAN) and output is on other side (JCH = 1,ICH)
c
c
c
        DO ICH = 1,NCHAN
          DO JCH = 1,ICH
            kcc = ksym(JCH,ICH)               ! unsure why this is needed  CPB
            CALL DA2(1,kcc,199,NRANG2SQ,ARRAY)
C
              KKK=0
              do III=1,NRANG2
               do JJJ=1,NRANG2
                KKK=KKK+1
                workcc(JJJ,III)=ARRAY(KKK)
               enddo
              enddo
C
C
        WRITE(jfile) ((workcc(III,JJJ),III=1,NRANG2),JJJ=1,NRANG2)
          ENDDO
        ENDDO
c
C
        CALL DA2(3,0,199,NRANG2SQ,ARRAY) !!! nuke WCC.DAT direct access file CPB
        deallocate(array)
        deallocate(workcc)
c
c

        IF (NCFGP.GT.0) THEN

C  Continuum-bound blocks
          allocate(wrkcb(NRANG2,NCFGP,NCHAN))
          DO ICH = 1,NCHAN
            READ (ifile) I,((wrkcb(M,N,I),N=1,NCFGP),M=1,NRANG2)
          ENDDO
          DO ICH = 1,NCHAN
            WRITE(jfile) ((wrkcb(M,N,ICH),N=1,NCFGP),M=1,NRANG2)
          ENDDO
          deallocate(wrkcb)

C  Bound-bound elements
          allocate(wrkbb(NCFGP))
          DO M = 1,NCFGP
            READ (ifile) (wrkbb(N),N=M,NCFGP)
            WRITE(jfile) (wrkbb(N),N=M,NCFGP)
          ENDDO
          deallocate(wrkbb)
        ENDIF

c  Asymptotic coefficients

        if( LAMAX.GT.0 ) then
         allocate(CF(NCHAN,NCHAN,LAMAX))
c
      OPEN(UNIT=2233,STATUS='UNKNOWN',FORM='UNFORMATTED',FILE='CFSCRAT')
c
          do I = 1, NCHAN
            READ(2233) ((CF(J,I,K),K=1,LAMAX),J=I,NCHAN) 
          enddo
c
          do I = 1, NCHAN
            WRITE(jfile) ((CF(J,I,K),K=1,LAMAX),J=I,NCHAN) 
          enddo
        deallocate(CF)
        close(2233)
        endif

        INAST2 = loop
      enddo

c-----------------------------------------------------------------------
c  Dipole matrix

      if( IPOLPH.le.1 ) goto 99
   10 continue

c    To be implemented...
      if(nproc.gt.1)then
        print *,'Using a single RECUPD.DAT dipole file....'
        print *,'Re-run with 1 RECUPH.DAT file if necessary.'
      endif

c  Open dipole matrix output file
 
      OPEN( UNIT=kfile, STATUS='UNKNOWN', FORM='UNFORMATTED',
     *      FILE='RECUPD.DAT')

c Loop over pairs of symmetries for dipole

      do loop2 = 2, INAST2
       do loop1 = 1,loop2-1
        MJ    = ILSPI(1,loop1)
        MMNP2 = ILSPI(2,loop1)
        MCHAN = ILSPI(3,loop1)
        NJ    = ILSPI(1,loop2)
        MNP2  = ILSPI(2,loop2)
        NCHAN = ILSPI(3,loop2)
        mcfg  = MMNP2-MCHAN*NRANG2
        ncfg  = MNP2 -NCHAN*NRANG2
        jdif  = (abs(MJ)-abs(NJ))/2
        if( MJ*NJ.lt.0 .and. abs(jdif).le.1 ) then
          if(iprint.ge.1)print*,'Dipole transition 2J+1:',mj,nj

c Selection rules ok.  read/write dipole matrix elements

          I = -2
          J = -2

c Continuum-continuum dipole
          allocate(wrkd(ND24,ND24,2))
          do nch = 1, NCHAN
            m2 = 0
            do mch = 1, MCHAN
              if( I.eq.-2 .and. J.eq.-2 ) READ (ifile,end=98) I,J,K,L
c              print*,i,j,k,l
              IF (I.EQ.-1 .AND. J.EQ.-1) GOTO 20
              if( I.eq.nch .and. J.eq.mch ) then
                I = -2
                J = -2
                READ (ifile) ((wrkd(M,N,1),M=1,K),N=1,L)
                READ (ifile) ((wrkd(M,N,2),M=1,K),N=1,L)
                WRITE(kfile) ((wrkd(M,N,1),M=1,K),N=1,L)
                WRITE(kfile) ((wrkd(M,N,2),M=1,K),N=1,L)
              else
                WRITE(kfile) (zero,M=1,NRANG2*NRANG2)
                WRITE(kfile) (zero,M=1,NRANG2*NRANG2)
              endif
              m2 = m2 + NRANG2
            enddo
            if( m2.lt.MMNP2 ) then
              if( I.eq.-2 .and. J.eq.-2 ) READ (ifile,end=98) I,J,K,L
c              print*,i,j,k,l
              IF (I.EQ.-1 .AND. J.EQ.-1) GOTO 20
              if( I.eq.nch .and. J.eq.0 ) then
                if( mcfg.ne.L ) then
                  print *,'STOP: mcfg and L must be the same:',mcfg,L
                  stop
                endif
                I = -2
                J = -2
                READ (ifile) ((wrkd(M,N,1),M=1,K),N=1,L)
                READ (ifile) ((wrkd(M,N,2),M=1,K),N=1,L)
                WRITE(kfile) ((wrkd(M,N,1),M=1,K),N=1,L)
                WRITE(kfile) ((wrkd(M,N,2),M=1,K),N=1,L)
              else
                WRITE(kfile) (zero,M=1,NRANG2*mcfg)
                WRITE(kfile) (zero,M=1,NRANG2*mcfg)
              endif
            endif
            WRITE(kfile) ((ABUTL,M=1,NRANG2),N=1,MCHAN)
            WRITE(kfile) ((ABUTV,M=1,NRANG2),N=1,MCHAN)
          enddo
 
c Bound-continuum dipole

          if( ncfg.gt.0 ) then
            do mch = 1, MCHAN
              if( I.eq.-2 .and. J.eq.-2 ) READ (ifile,end=98) I,J,K,L
c              print*,i,j,k,l
              IF (I.EQ.-1 .AND. J.EQ.-1) GOTO 20
              if( I.eq.0 .and. J.eq.mch ) then
                if( ncfg.ne.K ) then
                  print *,'STOP: ncfg and K must be the same:',ncfg,K
                  stop
                endif
                I = -2
                J = -2
                READ (ifile) ((wrkd(M,N,1),M=1,K),N=1,L)
                READ (ifile) ((wrkd(M,N,2),M=1,K),N=1,L)
                WRITE(kfile) ((wrkd(M,N,1),M=1,K),N=1,L)
                WRITE(kfile) ((wrkd(M,N,2),M=1,K),N=1,L)
              else
                WRITE(kfile) (zero,M=1,ncfg*NRANG2)
                WRITE(kfile) (zero,M=1,ncfg*NRANG2)
              endif
            enddo

c Bound-bound dipole

            if( mcfg.gt.0 ) then
              if( I.eq.-2 .and. J.eq.-2 ) READ (ifile,end=98) I,J,K,L
c              print*,i,j,k,l
              IF (I.EQ.-1 .AND. J.EQ.-1) GOTO 20
              if( I.eq.0 .and. J.eq.0 ) then
                if( mcfg.ne.L .or. ncfg.ne.K ) then
                  print *,'STOP: mcfg,ncfg and L,K=',mcfg,L,ncfg,K
                  stop
                endif
                I = -2
                J = -2
                READ (ifile) ((wrkd(M,N,1),M=1,K),N=1,L)
                READ (ifile) ((wrkd(M,N,2),M=1,K),N=1,L)
              else
                do N=1,mcfg
                  do M=1,ncfg
                    wrkd(M,N,1)=zero
                    wrkd(M,N,2)=zero
                  enddo
                enddo
              endif
              NDIMEN = NRANG2
              NTIMES = (ncfg-1)/NDIMEN + 1
              n2 = 0
              do nt = 1,NTIMES
                n22 = min(n2+NDIMEN,ncfg)
                write(0,*)n22,mcfg
                WRITE(kfile) ((wrkd(M,N,1),M=n2+1,n22),N=1,mcfg)
                WRITE(kfile) ((wrkd(M,N,2),M=n2+1,n22),N=1,mcfg)
                n2 = n22
              enddo
            endif
          endif
             deallocate(wrkd)
c Write zeros for dipole matrix Buttle correction

   20     READ (ifile,end=98) I,J,K,L
c          print*,i,j,k,l
          if( I.ne.-1 .and. J.ne.-1 .and. K.ne.-1 .and. L.ne.-1 ) then
            print*,'WARNING: dipole matrices not finished off:',i,j,k,l
          endif
          n2 = NCHAN*NRANG2
          if( n2.LT.MNP2 ) then
            WRITE(kfile) ((ABUTL,J=1,MNP2-n2),I=1,MCHAN)
            WRITE(kfile) ((ABUTV,J=1,MNP2-n2),I=1,MCHAN)
          endif
          do MCH = 1,MCHAN
            WRITE(kfile) ((BBUTL,J=1,NCHAN),I=1,NRANG2)
            WRITE(kfile) ((BBUTV,J=1,NCHAN),I=1,NRANG2)
          enddo
          m2 = MCHAN*NRANG2
          if( m2.LT.MMNP2 ) then
            WRITE(kfile) ((BBUTL,J=1,NCHAN),I=1,MMNP2-m2)
            WRITE(kfile) ((BBUTV,J=1,NCHAN),I=1,MMNP2-m2)
          endif
          WRITE(kfile) ((BBUTL,J=1,NCHAN),I=1,MCHAN)
          WRITE(kfile) ((BBUTV,J=1,NCHAN),I=1,MCHAN)

c CLEBSCH GORDAN COEFFICIENTS FOR THE POLARIZABILITY CALC

          READ (ifile) MAXM1
          if(MAXM1.gt.MZLR2)then
            print *,'STOP! Increase dimension MZLR2 to:',MAXM1
            stop
          endif

          READ (ifile) (temp(M),M=1,MAXM1)
          WRITE(kfile) MAXM1,(temp(M),M=1,MAXM1)

c ANGULAR INTEGRALS NEEDED FOR OUTER REGION INTEGRATION
c
c
          allocate(AC(NCHAN,MCHAN),BLC(NCHAN,MCHAN),BVC(NCHAN,MCHAN))
c
c
          if(outer) then
            READ (ifile,end=96) ((AC(I,J),J=1,MCHAN),I=1,NCHAN)
            READ (ifile) ((BLC(I,J),J=1,MCHAN),I=1,NCHAN)
            READ (ifile) ((BVC(I,J),J=1,MCHAN),I=1,NCHAN)
            if(iprint.ge.1)print*,'outer region contributions found'
            go to 97
          endif
   96     print*,'outer region contributions not found: zeros inserted'
          do j=1,MZCHF
            do i=1,MZCHF
              AC(i,j)  = zero
              BLC(i,j) = zero
              BVC(i,j) = zero
            enddo
          enddo
   97     WRITE(kfile) ((AC (I,J),J=1,MCHAN),I=1,NCHAN)
          WRITE(kfile) ((BLC(I,J),J=1,MCHAN),I=1,NCHAN)
          WRITE(kfile) ((BVC(I,J),J=1,MCHAN),I=1,NCHAN)
          deallocate(AC,BLC,BVC)
        endif
       enddo
      enddo

   98 if(iprint.ge.1)print*,'RECUPD.DAT created'
      close(kfile)
   99 if(iprint.ge.1.and.nproc.eq.1)print*,'RECUPH.DAT created'
      if(iprint.ge.1.and.nproc.gt.1)print*,'RECUPH.XXX created'
      close(jfile)
      close(ifile)

      stop 'Normal End'
      end

C-----------------------------------------------------------------------

      SUBROUTINE BUTFIT(IMAX,E,F,RA,EMAX,ALPHA,BETA,NBUT,DELTA)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C  FITTING OF BUTTLE CORRECTIONS -- M J SEATON, J.PHYS.B20(1987)L69-72.
C
C  IMAX = NUMBER OF POINTS FOR WHICH CORRECTION CALCULATED.
C  E(I) = ENERGY POINTS
C  F(I) = BUTTLE CORRECTION
C  RA = BOUNDARY RADIUS
C  EMAX = EIGENS(NRANG2,L)
C  ALPHA, BETA = FIT PARAMETERS, NBUT =...
C  DELTA = ACCURACY ACHIEVED
C
      LOGICAL POLE
C
      DIMENSION E(IMAX),F(IMAX)
C-----------------------------------------------------------------------
C
C  INITIALISATIONS
C
      D = RA*RA
      DO 10 I = 1,IMAX
        E(I) = E(I)*D
   10 CONTINUE
      T = 0.5D0
      IF(EMAX.GT.0.0D0)T=T + RA*SQRT(EMAX)/3.141592654D0
      NBUT = T
      DELTA = 1.0D30
      ALPHA = 1.0D0
      BETA = 0.0D0
C
C START ITERATIONS FOR FIT
C
      DO 50 KK = 1,15
C
        X11 = 0.0D0
        X12 = 0.0D0
        X22 = 0.0D0
        Y1 = 0.0D0
        Y2 = 0.0D0
        DELTA0 = DELTA
        DELTA = 0.0D0
C
C    START SUM OVER POINTS I
C
        DO 40 I = 1,IMAX
          U = BETA + E(I)
C
C    CALCULATE FUNCTIONS B(U) AND C(U)
C
          B = 0.0D0
          C = 0.0D0
C    CASE OF U.GT.0.04
          IF (U.GT.0.04D0) THEN
            FK = SQRT(U)
            POLE = .FALSE.
            G = -1.5707963D0
            DO 20 N = 0,NBUT
              G = G + 3.141592654D0
              IF (ABS(FK-G).GT.0.3D0) THEN
                A = 1.0D0/ (U-G*G)
                B = B + A
                C = C - A*A
C
              ELSE
                POLE = .TRUE.
                D1 = FK - G
              ENDIF
C
   20       CONTINUE
            IF (POLE) THEN
              D2 = D1*D1
              D = 0.33333333D0*
     A            D1* (1.0D0+0.066666667D0*
     B            D2* (1.0D0+0.0952381D0*D2))
              A = 1.0D0/ (2.0D0*FK-D1)
              BB = (D+A)/FK
              D = 0.33333333D0* (1.0D0+
     A            D2* (0.2D0+0.031746032D0*D2))
              C = 2.0D0*C + 0.5D0* (D-A*A-BB)/U
              B = 2.0D0*B + BB
C
            ELSE
              T = TAN(FK)
              TK = T/FK
              B = 2.0D0*B + TK
              C = 2.0D0*C + 0.5D0*
     A            (1.0D0+T*T-TK)/U
            ENDIF
C
C    SUM FOR U.LE..04
          ELSE
            G = -1.5707963D0
            DO 30 N = 0,NBUT
              G = G + 3.141592654D0
              A = 1.0D0/ (U-G*G)
              B = B + A
              C = C - A*A
   30       CONTINUE
C
C    CASE OF U.LT..04 AND U.GT.-.04
            IF (U.GT.-0.04D0) THEN
              B = (0.4D0*U+1.0D0)
     A            *U*0.33333333D0
     A            + B*2.0D0 + 1.0D0
              C = ((0.48571429D0*U+0.8D0)*U
     A            +1.0D0)*0.33333333D0
     B            + C*2.0D0
C
C    CASE OF U.LT.-.04
            ELSE
              FK = SQRT(-U)
              T = TANH(FK)
              TK = T/FK
              B = 2.0D0*B + TK
              C = 2.0D0*C + 0.5D0
     A            * (1.0D0-T*T-TK)/U
            ENDIF
C
          ENDIF
C
C    INCREMENT MATRICES X AND Y
          DF = F(I) - ALPHA*B
          DD = ABS(DF)
          IF (DELTA.LT.DD) DELTA = DD
          X11 = X11 + B*B
          X12 = X12 + B*C
          X22 = X22 + C*C
          Y1 = Y1 + DF*B
          Y2 = Y2 + DF*C
C
   40   CONTINUE
C
C  SOLVE EQUATIONS AND INCREMENT ALPHA AND BETA
C
C+++ MODIFICATIONS MADE BY MJS, 22.12.86.
CJZ    X12=ALPHA*X12
        DET = 1.0D0/ (X11*X22-X12*X12)
CJZ    BETA=BETA+DET*(-X12*Y1+X11*Y2)
        BETA = BETA + DET* (-X12*Y1+X11*Y2)/ALPHA
C+++ END MODIFICATIONS
        ALPHA = (X22*Y1-X12*Y2)*DET + ALPHA
C
C    CHECK CONVERGENCE
C
        IF (DELTA.LT.2.0D-4) GOTO 60
C
C  END ITERATIONS
C
   50 CONTINUE
      NBUT = -NBUT
      IF (DELTA.GT.DELTA0) NBUT = 0
C     AS PROCEDURE NOT CONVERGED OR DISTINCTLY DIVERGING - WE'88NOV17.
C
   60 CONTINUE
C
      END

C-----------------------------------------------------------------------

      SUBROUTINE HEAD(jfilee)
      use datamod
      IMPLICIT REAL*8 (A-H,O-Z)

      INCLUDE 'PARAM'

      COMMON/COM1/ENDS(MZNR2,MZLR2),ENAT(MZTAR),
     *            EIGENS(MZNR2,MZLR2),COEFF(3,MZLR2)
      COMMON/COM2/JAT(MZTAR),LPTY(MZTAR), MAXNLG(MZLR2)
      COMMON/COM3/NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,LAM,
     *            NAST,NIX,RA,BSTO,HINT
C
c Write header to standard RECUPH.DAT hamilitonian output file
C
      WRITE(jfilee) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM,IZESP,
     *             JRELOP,JRELOP,JRELOP
      WRITE(jfilee) (MAXNLG(L),L=1,LRANG1), (MAXNLG(L),L=1,LRANG1),
     *             (maxnc,L=1,LRANG1)
C
      DO L = 1,LRANG2
        WRITE(jfilee) (2*EIGENS(N,L),N=1,NRANG2)     !a.u. to Ryd - nrb
        WRITE(jfilee) (ENDS(N,L),N=1,NRANG2)
      ENDDO
C
      WRITE(jfilee) RA,BSTO,HINT,eta,eta,NIX
      WRITE(jfilee) ((COEFF(I,L),I=1,3),L=1,LRANG2)
      WRITE(jfilee) NAST
      WRITE(jfilee) (ENAT(I),I=1,NAST), (JAT(I),I=1,NAST),
     *             (isat,I=1,NAST), (LPTY(I),I=1,NAST)
C
      RETURN
      END
      SUBROUTINE DA2(KEY,IREC,JDISC,LENGTH,ARRAY)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'PARAM'
C-----------------------------------------------------------------------
C
C    TO STORE A LARGE ARRAY IN A DA FILE OF RECORD LENGTH 8*LREC BYTES.
C
C    KEY = 1 FOR READ,
C        = 2 FOR WRITE,
C        = 0 FOR FINDING NUMBER OF DA RECORDS A GIVEN ARRAY TAKES.
C
C    IREC= (ON CALL) POINTER TO FIRST DA RECORD FOR ARRAY,
C        = 0 FOR OPENING DA FILE (BY NAME),
C        =-1 FOR OPENING DA FILE (SCRATCH),
C        = (ON RETURN) POINTER TO NEXT AVAILABLE DA RECORD.
C
C    JDISC = DA FILE UNIT NUMBER.
C
C    ARRAY(LENGTH) = ARRAY TO READ OR WRITE.
C
C-----------------------------------------------------------------------
      PARAMETER (LREC=400)
C
      DIMENSION ARRAY(LENGTH)

C-----------------------------------------------------------------------
C
      IF (IREC.GT.0) GOTO 20
C
      IRECL = 8*LREC
C
C
      if(LREC.ne.LENGTH)then
      write(0,*)'SET LREC in DA2 to NRANG2**2'
      stop
      endif
C
      IF (KEY.EQ.2) THEN
        OPEN (JDISC,STATUS='UNKNOWN',FILE='WCC.DAT',ACCESS='DIRECT',
     A       RECL=IRECL)
      RETURN
      ENDIF
C
      IF (KEY.EQ.3) THEN
        CLOSE(JDISC,status='delete')
      RETURN
      ENDIF
C
C
   20 CONTINUE
      I1 = 1 
      I2 = LENGTH 
C
      IF (KEY.EQ.2) THEN
        WRITE (JDISC,REC=IREC) (ARRAY(I),I=I1,I2)
      ELSE
        READ (JDISC,REC=IREC) (ARRAY(I),I=I1,I2)
      ENDIF
C
C
      END
