!******************************************************************************
!                                                                             *
!                                                                             *
!                                                                             *
!               DIFFERENTIAL ALGEBRA PACKAGE OF M. BERZ                       *
!                     ****************************                            *
!                               holy3                                         *
!                                                                             *
!                                                                             *
!                                                                             *
!         VERSION FOR MACHINE IN LINE THAT IS NOT COMMENTED OFF               *
!        TO CREATE DIFFERENT VERSIONS, USE THE PROGRAM 'VERSION'              *
!                                                                             *
!                                                                             *
!                                                                             *
!                                                                             *
!        THIS PACKAGE WAS INITIALLY WRITTEN BY PROF. M. BERZ WHILE AT         *
!        THE LAWRENCE BERKELEY LABORATORY.                                    *
!        IT HAS BEEN EXTENSIVELY MODIFIED BY THE MEMBERS OF THE ESG GROUP.    *
!        THEREFORE PROF. BERZ SHOULD NOT BE HELD RESPONSIBLE FOR ANY BUGS.    *
!                                                                             *
!                  NEW RULES OF THE GAME (EXHAUSTIVE)                         *
!                 **********************************                          *
!                         THERE ARE NONE                                      *
!                                                                             *
!******************************************************************************
!
!
!     THIS FILE CONTAINS ROUTINES TO PERFORM DIFFERENTIAL ALGEBRA (DA)
!     AS AN OPTION, ALSO COMPONENTWISE ALGEBRA (CA) CAN BE PERFORMED.
!     A DESCRIPTION OF THE INTERNAL ARRAYS USED BY THE ROUTINES CAN
!     BE FOUND IN BLOCKDATA DABLD.
!
!
!     SHORT REFERENCE CHART
!     *********************
!
!     THE PARAMETERS USED BELOW HAVE THE FOLLOWING MEANING:
!
!     A,B:                NAME OF INPUT DA VECTORS   (INTEGER)
!     C:                  NAME OF OUTPUT DA VECTOR   (INTEGER)
!     X,Y:                NAME OF INPUT DA MATRIX    (INTEGER(...))
!     Z:                  NAME OF OUTPUT DA MATRIX   (INTEGER(...))
!
!     F:                  NAME OF A DA FUNCTION      (CHARACTER*4)
!     G:                  NAME OF EXTERNAL FUNCTION  (DOUBLE PRECISION)
!     JJ:                 ARRAY OF EXPONENTS         (INTEGER(20))
!     O:                  ORDER                      (INTEGER)
!     N:                  NUMBER OF VARIABLES        (INTEGER)
!     I,J,K:              INTEGER NUMBER             (INTEGER
!     R,RA,RB:            REAL NUMBERS               (DOUBLE PRECISION)
!     H:                  ARRAY OF LENGTH LH         (DOUBLE PRECISION)
!     U:                  OUTPUT UNIT NUMBER         (INTEGER)
!     T:                  COMMENT TEXT               (CHARACTER*10)
!
!
!               SUBROUTINES AND THEIR CALLING PARAMETERS
!               ****************************************
!
!     DAINI(O,N,U):       INITIALIZES CONTROL ARRAYS AND SETS MAX. ORDER O AND
!                         MAX. NUMBER OF VARIABLES N. MUST BE CALLED BEFORE ANY
!                         OTHER DA ROUTINE CAN BE USED.
!
!     DAALL(A,I,T,O,N):   ALLOCATES SPACE FOR I VECTORS A. T: CHARACTER NAME
!     DADAL(A,I):         DEALLOCATES THE I VECTORS A.
!!     DAVAR(A,R,I):       MAKES A INDEPENDENT VARIABLE # I WITH INITIAL VALUE R
!!     DACON(A,R):         SETS A TO CONSTANT R
!     DANOT(O):           SETS NEW TRUNCATION ORDER O FOR DA OPERATIONS
!     DAEPS(R):           SETS NEW ZERO TOLERANCE EPSILON
!
!!     DAPEK(A,JJ,R):      RETURNS COEF R OF MONOMIAL WITH EXPONENTS JJ OF A
!!     DAPOK(A,JJ,R):      SETS COEF OF MONOMIAL WITH EXPONENTS JJ OF A TO R
!
!!     DACOP(A,C):         PERFORMS C = A
!!     DAADD(A,B,C):       PERFORMS C = A + B
!!    DASUB(A,B,C):       PERFORMS C = A - B
!!     DAMUL(A,B,C):       PERFORMS C = A * B
!!     DADIV(A,B,C):       PERFORMS C = A / B
!!     DASQR(A,C):         PERFORMS C = A^2           (SQUARE OF A)
!
!!     DACAD(A,RA,C):      PERFORMS C = A + RA
!!     DACSU(A,RA,C):      PERFORMS C = A - RA
!!     DASUC(A,RA,C):      PERFORMS C = RA - A
!!     DACMU(A,RA,C):      PERFORMS C = A * RA
!!    DACDI(A,RA,C):      PERFORMS C = A / RA
!!     DADIC(A,RA,C):      PERFORMS C = RA / A
!!     DACMA(A,B,RB,C):    PERFORMS C = A + RB*B
!!DAMULIN(A,B,RA,C,D,RB,C):    PERFORMS C = A*B*RA + C*D*RB
!!     DALIN(A,RA,B,RB,C): PERFORMS C = A*RA + B*RB
!!     DAFUN(F,A,C):       PERFORMS C = F(A)          (DA FUNCTION)
!
!!     DAABS(A,R):         PERFORMS R = |A|           (NORM OF A)
!!     DACOM(A,B,R):       PERFORMS R = |A-B|         (NORM OF A-B)
!!     DAPOS(A,C):         PERFORMS C(I) = |A(I)|     (MAKE SIGNS POSITIVE)
!
!!     DACCT(X,I,Y,J,Z,K)  CONCATENATES Z = X o Y;   I,J,K: # OF VECTORS IN X,Y,
!!     DAINV(X,I,Z,K)      INVERTS Z = X^-1;           I,J: # OF VECTORS IN X,Y
!!     DAPIN(X,I,Z,K,JJ)   PARTIALLY INVERTS Z = X^-1; I,J: # OF VECTORS IN X,Y,
!                         JJ: ARRAY; NONZERO ENTRIES DENOTE TO BE INVERTED LINES
!
!!     DADER(I,A,C):       PERFORMS C = DA/DI (DERIV. WITH RESPECT TO VARIABLE I
!!     DAPOI(A,B,C,I):     PERFORMS C = [A,B] (POISSON BRACKET, 2*I: # PHASEVARS
!!     DACFU(A,G,C):       MULTIPLIES COEFFICIENTS WITH FUNCTION G(JJ)
!     DAIMP(H,LH,A):      "IMPORTS" THE ARRAY H WITH LENGTH LH TO DA VAR A
!     DAEXP(H,LH,A):      "EXPORTS" THE DA VAR A TO ARRAY H WITH LENGTH LH
!     DAPRI(A,U):         PRINTS DA VECTOR A TO UNIT U
!     DAREA(A,U):         READS DA VECTOR A FROM UNIT U
!     DADEB(U,T,I):       DEBUGGER, DUMPS TO U. T: MEMO, I=0: RETURN, I=1:STOP
!!     DARAN(A,R,seed):         FILLS A WITH RANDOM NUMBERS. R: FILLFACTOR
!     DANUM(O,N,I):       COMPUTES NUMBER OF MONOMIALS IN N VAR THROUGH ORDER O
!
!
!     ADDITIONAL ROUTINES THE USER DOES NOT NEED TO CALL:
!
!     DAINF: RETURNS INFOS ABOUT A DA VECTOR PREVIOUSLY DECLARED
!     DAPAC: PACKS DA VECTORS
!     DACHK: CHECKS IF DA VECTORS HAVE COMPATIBLE ATTRIBUTES
!     DCODE: TRANSFORMS DIGITS IN A CERTAIN BASE TO A DECIMAL INTEGER
!     NCODE: EXTRACTS DIGITS IN A CERTAIN BASE FROM A DECIMAL INTEGER
!
!
!     FURTHER WISHES
!     **************
!
!     - CHECK DAREA AND DAPRI FOR CA VECTORS
!     - MAKE DAFUN USE DASQR
!
!
!      BLOCKDATA DABLD
!     ***************
!
!
!     PARAMETERS:
!
!     LDA: MAXIMUM NUMBER OF DA-VECTORS;    CAN BE CHANGED QUITE ARBITRARILY
!     LST: LENGTH OF MAIN STORAGE STACK;    CAN BE CHANGED QUITE ARBITRARILY
!     LEA: MAXIMUM NUMBER OF MONOMIALS;     CAN BE INCREASED FOR LARGE NO,NV
!     LIA: DIMENSION OF IA1,IA2;            CAN BE INCREASED FOR LARGE NO,NV
!     LNO: MAXIMUM ORDER;                   CAN BE INCREASED TO ABOUT 1000
!     LNV: MAXIMUM NUMBER OF VARIABLES;     CAN BE INCREASED TO ABOUT 1000
!
!     ALL THE CHANGES IN THE VALUES OF PARAMETERS HAVE TO BE MADE BY GLOBAL
!     SUBSTITUTIONS IN ALL SUBROUTINES.
!
!     DANAME:   NAME OF DA VECTOR
!
!     CC:       STACK OF DOUBLE PRECISON COEFFICIENTS
!     I1:       FIRST CHARACTERISTIC INTEGER (CF DAINI)
!     I2:       SECOND CHARACTERISTIC INTEGER (CF DAINI)
!
!     IE1:      CHARACTERISTIC INTEGER 1 OF UNPACKED REPRESENTATION (CF DAINI)
!     IE2:      CHARACTERISTIC INTEGER 2 OF UNPACKED REPRESENTATION (CF DAINI)
!     IEO:      ORDER OF ENTRY IN UNPACKED REPRESENTATION
!     IA1:      REVERSE TO IE1 (CF DAINI)
!     IA2:      REVERSE TO IE2 (CF DAINI)
!
!     IDANO:    ORDER OF DA VECTOR; IN CA, NUMBER OF COMPONENTS
!     IDANV:    NUMBER OF VARIABLES; IF 0, INDICATES CA VECTOR
!     IDAPO:    FIRST ADDRESS IN STACK
!     IDALM:    NUMBER OF RESERVED STACK POSITIONS
!     IDALL:    NUMBER OF MOMENTARILY REQUIRED STACK POSITIONS
!
!     NDA:      NUMBER OF DA VECTORS MOMENTARILY DEFINED
!     NST:      NUMBER OF STACK POSITIONS MOMENTARILY ALLOCATED
!     NOMAX:    MAXIMUM REQUESTED ORDER  (CF DAINI)
!     NVMAX:    MAXIMUM REQUESTED NUMBER OF VARIABLES (CF DAINI)
!     NMMAX:    MAXIMUM NUMBER OF MONOMIALS FOR NOMAX, NVMAX (CF DAINI)
!     NOCUT:    MOMENTARY TRUNCATION ORDER
!     EPS:      TRUNCATION ACCURACY (CAN BE SET BY USER)
!     EPSMAC:   MANTISSA LENGTH OF MACHINE (PESSIMISTIC ESTIMATE)
!
!-----------------------------------------------------------------------------1
!
      subroutine daini(no,nv,iunit)
      implicit none
      integer i,iall,ibase,ic1,ic2,icmax,io1,io2,iout,iunit,j,jd,jj,jjj,&
     &jjjj,jl,js,k,n,nn,no,nv
!     *****************************
!
!     THIS SUBROUTINE SETS UP THE MAJOR ORDERING AND ADDRESSING ARRAYS IN
!     COMMON BLOCK DAINI. IF IUNIT > 0, THE ARRAYS WILL BE PRINTED TO UNIT
!     NUMBER IUNIT. AN EXAMPLE FOR THE ARRAYS GENERATED BY DAINI CAN BE
!     FOUND AFTER THE ROUTINE.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
!      COMMON / DASCR /  IS(20), RS(20)                                        1
      integer idao,is,iscrri
      double precision rs
      common/dascr/is(100),rs(100),iscrri(100),idao
!-----------------------------------------------------------------------------2

      logical allvec(lda)
      integer nhole
      common /hole/nhole
      common /alloc/ allvec

!
      character aa*10
      dimension n(lnv+1),k(0:lnv),j(lnv),jj(lnv)
!
      write(*, *) "daini:", no, nv
      if(eps.le.0.d0) eps=1.d-38
!      if(EPS.le.0.d0) eps=1.d-90
      epsmac=1.d-7
      if(nv.eq.0) return
      ndamaxi=0
!
      do i=1, lda
        allvec(i) = .false.
      enddo
      nhole=0
!*****************************************

!     INITIALIZING VARIABLES IN COMMON / DA /
!     ***************************************
!
      nda   = 0
      nst   = 0
      nomax = no
      nvmax = nv
      call danum(no,nv,nmmax)
      nocut = no
      lfi   = 0
!
      do i=0,lia
        ia1(i) = 0
        ia2(i) = 0
      enddo
!
      do i=1,100
        is(i) = 0
      enddo
!
      if(nv.gt.lnv.or.no.gt.lno) then
        write(6,*)'ERROR IN SUBROUTINE DAINI, NO, NV = ',no,nv
        call dadeb(31,'ERR DAINI ',1)
      endif
!
      ibase = no+1
      js    = nv/2
      if(float(ibase)**((nv+1)/2).gt.float(lia)) then
        write(6,*)'ERROR, NO = ',no,', NV = ',nv,' TOO LARGE FOR',      &
     &  ' LIA = ',lia
        call dadeb(31,'ERR DAINI ',1)
      endif
!
      icmax = 0
      nn    = 0
      k(0)  = 0
!
      do io2=0,no
!     ***************
!
        n(1)  = io2
        jl    = 0
        jd    = 1
!
  50    jl    = jl + jd
!

!old
!      IF(JL.EQ.0) THEN
!old
!     modified according to Wu Ying
!
        if(jl.le.0) then
!
          goto 100
        elseif(jd.eq.1) then
          j(jl) = 0
        else
          j(jl) = j(jl) + 1
        endif
!
        k(jl)    = k(jl-1)*ibase + j(jl)
        n(jl+1)  = n(jl) - j(jl)
!
        if(j(jl).gt.n(jl)) then
          jd    = -1
          goto 50
        elseif(jl.lt.js) then
          jd    = 1
          goto 50
        else
          j(jl) = n(jl)
          k(jl) = k(jl-1)*ibase + j(jl)

          ic2   = k(jl)
          icmax = max(icmax,ic2)
          k(jl) = 0
!
          ia2(ic2) = nn
!
          do io1=0,no-io2
!        ******************
!
            n(js+1) = io1
            jd      = 1
!
  70        jl      = jl + jd
!
            if(jl.eq.js) then
              goto 80
            elseif(jd.eq.1) then
              j(jl) = 0
            else
              j(jl) = j(jl) + 1
            endif
!
            k(jl)    = k(jl-1)*ibase + j(jl)
            n(jl+1)  = n(jl) - j(jl)
!
            if(j(jl).gt.n(jl)) then
              jd    = -1
              goto 70
            elseif(jl.lt.nv) then
              jd    = 1
              goto 70
            else
              jd    = -1
              j(jl) = n(jl)
              k(jl) = k(jl-1)*ibase + j(jl)
              ic1   = k(jl)
              icmax = max(icmax,ic1)
              nn = nn + 1
!
              ie1(nn) = ic1
              ie2(nn) = ic2
              i1 (nn) = ic1
              i2 (nn) = ic2
              if(ic2.eq.0) ia1(ic1) = nn
              ieo(nn) = io1 + io2
!
              goto 70
            endif
!
   80       continue
          enddo
!
          jd = -1
          goto 50
        endif
!
  100   continue
      enddo
!
      if(nn.gt.lea) then
        write(6,*)'ERROR IN DAINI, NN = ',nn,' EXCEEDS LEA'
        call dadeb(31,'ERR DAINI ',1)
      endif
!
!     ALLOCATING SCRATCH VARIABLES
!     ****************************
!
      iall = 0
      call daall(iall,1,'$$UNPACK$$',nomax,nvmax)
!
      do i=0,nomax
        aa = '$$MUL   $$'
        write(aa(6:10),'(I5)') i
        iall = 0
!      CALL DAALL(IALL,1,AA,I,NVMAX)
        call daall(iall,1,aa,nomax,nvmax)
      enddo
!
      idall(1) = nmmax
!
!     DOUBLE CHECKING ARRAYS IE1,IE2,IA1,IA2
!     **************************************
!
      do i=1,nmmax
!
        jjj = ia1(ie1(i)) + ia2(ie2(i))
        if(jjj.ne.i) then
          write(6,*)'ERROR IN DAINI IN ARRAYS IE1,IE2,IA1,IA2 AT I = ',i
          call dadeb(31,'ERR DAINI ',1)
        endif
!
      enddo
!
      if(iunit.eq.0) return
!
      write(6,*)'ARRAY SETUP DONE, BEGIN PRINTING'
!
      iout = 32
      open(iout,file='DAINI.DAT',status='NEW')
!CRAY OPEN(IOUT,FILE='DAINI',STATUS='UNKNOWN',FORM='FORMATTED')          *CRAY
!CRAY REWIND IOUT                                                        *CRAY
!
      write(iout,'(/A/A/)') ' ARRAYS I1 THROUGH I20, IE1,IE2,IEO',      &
     &' **********************************'
      do i=1,nmmax
        call dancd(ie1(i),ie2(i),jj)
        write(iout,'(1X,I5,2X,4(5I2,1X),3I6)') i,(jj(jjjj),jjjj=1,lnv), &
     &  ie1(i),ie2(i),ieo(i)
      enddo
!
      write(iout,'(/A/A/)') ' ARRAYS IA1,IA2',' **************'
      do i=0,icmax
        write(iout,'(3I10)') i,ia1(i),ia2(i)
      enddo
!
      return
      end

      subroutine daexter
      implicit none
      integer i
!     *****************************
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
      logical allvec(lda)
      integer nhole
      common /hole/nhole
      common /alloc/ allvec

!
      do i=1, lda
        allvec(i)=.false.
      enddo

      return
      end

      subroutine dallsta(ldanow)
      implicit none
      integer i,ldanow
!     *****************************
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
      logical allvec(lda)
      integer nhole
      common /hole/nhole
      common /alloc/ allvec

!
      ldanow=0
      do i=1, lda
        if(allvec(i)) ldanow=ldanow+1
      enddo

      write(6,*) ' ALLOCATED ',ldanow

      return
      end

!
! EXAMPLE: ARRAYS I1 THROUGH I20, IE1,IE2,IEO (NOMAX=3,NVMAX=4)
! *************************************************************
!     I   I1               THROUGH               I20     IE1   IE2   IEO
!     1   0 0 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      0     0     0
!     2   1 0 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      1     0     1
!     3   0 1 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      4     0     1
!     4   2 0 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      2     0     2
!     5   1 1 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      5     0     2
!     6   0 2 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      8     0     2
!     7   3 0 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      3     0     3
!     8   2 1 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      6     0     3
!     9   1 2 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0      9     0     3
!    10   0 3 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0     12     0     3
!    11   0 0 0 0 0  0 0 0 0 0  1 0 0 0 0  0 0 0 0 0      0     1     1
!    12   1 0 0 0 0  0 0 0 0 0  1 0 0 0 0  0 0 0 0 0      1     1     2
!    13   0 1 0 0 0  0 0 0 0 0  1 0 0 0 0  0 0 0 0 0      4     1     2
!    14   2 0 0 0 0  0 0 0 0 0  1 0 0 0 0  0 0 0 0 0      2     1     3
!    15   1 1 0 0 0  0 0 0 0 0  1 0 0 0 0  0 0 0 0 0      5     1     3
!    16   0 2 0 0 0  0 0 0 0 0  1 0 0 0 0  0 0 0 0 0      8     1     3
!    17   0 0 0 0 0  0 0 0 0 0  0 1 0 0 0  0 0 0 0 0      0     4     1
!    18   1 0 0 0 0  0 0 0 0 0  0 1 0 0 0  0 0 0 0 0      1     4     2
!    19   0 1 0 0 0  0 0 0 0 0  0 1 0 0 0  0 0 0 0 0      4     4     2
!    20   2 0 0 0 0  0 0 0 0 0  0 1 0 0 0  0 0 0 0 0      2     4     3
!    21   1 1 0 0 0  0 0 0 0 0  0 1 0 0 0  0 0 0 0 0      5     4     3
!    22   0 2 0 0 0  0 0 0 0 0  0 1 0 0 0  0 0 0 0 0      8     4     3
!    23   0 0 0 0 0  0 0 0 0 0  2 0 0 0 0  0 0 0 0 0      0     2     2
!    24   1 0 0 0 0  0 0 0 0 0  2 0 0 0 0  0 0 0 0 0      1     2     3
!    25   0 1 0 0 0  0 0 0 0 0  2 0 0 0 0  0 0 0 0 0      4     2     3
!    26   0 0 0 0 0  0 0 0 0 0  1 1 0 0 0  0 0 0 0 0      0     5     2
!    27   1 0 0 0 0  0 0 0 0 0  1 1 0 0 0  0 0 0 0 0      1     5     3
!    28   0 1 0 0 0  0 0 0 0 0  1 1 0 0 0  0 0 0 0 0      4     5     3
!    29   0 0 0 0 0  0 0 0 0 0  0 2 0 0 0  0 0 0 0 0      0     8     2
!    30   1 0 0 0 0  0 0 0 0 0  0 2 0 0 0  0 0 0 0 0      1     8     3
!    31   0 1 0 0 0  0 0 0 0 0  0 2 0 0 0  0 0 0 0 0      4     8     3
!    32   0 0 0 0 0  0 0 0 0 0  3 0 0 0 0  0 0 0 0 0      0     3     3
!    33   0 0 0 0 0  0 0 0 0 0  2 1 0 0 0  0 0 0 0 0      0     6     3
!    34   0 0 0 0 0  0 0 0 0 0  1 2 0 0 0  0 0 0 0 0      0     9     3
!    35   0 0 0 0 0  0 0 0 0 0  0 3 0 0 0  0 0 0 0 0      0    12     3
!
!    ARRAYS IA1,IA2
!    **************
!    I        IA1       IA2
!    0         1         0   IE1,IE2 AND IA1,IA2 ALLOW THE EASY COMPUTATION
!    1         2        10   OF THE ADDRESS OF THE PRODUCT OF TWO MONOMIALS.
!    2         4        22   LET IX AND IY BE THE POSITIONS OF THE TWO
!    3         7        31   FACTORS. THEN THE POSITION IZ OF THE PRODUCT OF
!    4         3        16   THE TWO FACTORS IS GIVEN BY
!    5         5        25
!    6         8        32   IZ = IA1(IE1(IX)+IE1(IY)) + IA2(IE2(IX)+IE2(IY))
!    7         0         0
!    8         6        28
!    9         9        33   THE OTHER VARIABLES SET BY DAINI WOULD HAVE THE
!   10         0         0   VALUES
!   11         0         0
!   12        10        34   NOMAX = 3,  NVMAX = 4, NMMAX = 35
!
      subroutine daallno(ic,l,ccc)
      implicit none
      integer i,ind,l,ndanum,no,nv
      double precision x
!     ********************************
!
!     THIS SUBROUTINE ALLOCATES STORAGE FOR A DA VECTOR WITH
!     ORDER NOmax AND NUMBER OF VARIABLES NVmax
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"

      logical allvec(lda)
      integer nhole
      common /hole/nhole
      common /alloc/ allvec

!-----------------------------------------------------------------------------9
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      integer ic(*)
      logical incnda
      character c*10,ccc*10
!
      no=nomax
      nv=nvmax
      ind = 1
      do i=1,l
        if(ic(i).gt.0.and.ic(i).le.nda) then
!         DANAME(IC(I)) = C
!         IF(IDANO(IC(I)).EQ.NO.AND.IDANV(IC(I)).EQ.NV) THEN
        else
          if(nv.ne.0.and.(no.gt.nomax.or.nv.gt.nvmax)) then
            write(6,*)'ERROR IN DAALL, VECTOR ',c,' HAS NO, NV = ',     &
     &      no,nv,' NOMAX, NVMAX = ',nomax,nvmax
            call dadeb(31,'ERR DAALL ',1)
          endif
!
          if(nhole.gt.0) then
            ind=nda
20          if (allvec(ind)) then
              ind = ind - 1
              goto 20
            endif
            incnda = .false.
            nhole=nhole-1
          else
            incnda = .true.
            nda = nda + 1
            ind=nda
            if(nda.gt.lda) then
              write(6,*)'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHA',&
     &        'USTED'
              call dadeb(31,'ERR DAALL ',1)
            endif
          endif

          allvec(ind) = .true.
          ic(i) = ind
!
          if(nv.ne.0) then
            call danum(no,nv,ndanum)
          else
            ndanum = no
          endif

          c = ccc
          if(l.ne.1) write(c(6:10),'(I5)') i
          daname(ind) = c

          if (incnda) then
            if(ind.gt.nomax+2) then
              idano(ind) = nomax
              idanv(ind) = nvmax
              idapo(ind) = nst + 1
              idalm(ind) = nmmax
              idall(ind) = 0
              nst = nst + nmmax
            else
              idano(ind) = no
              idanv(ind) = nv
              idapo(ind) = nst + 1
              idalm(ind) = ndanum
              idall(ind) = 0
              nst = nst + ndanum
            endif
          endif
!
          if(nst.gt.lst) then
            x=-1.d0
            write(6,*)'ERROR IN DAALL, STACK EXHAUSTED '
            write(6,*) ' NST,LST '
            write(6,*)  nst,lst
            write(6,*) ' NDA,NDANUM,NDA*NDANUM '
            write(6,*)  nda,ndanum,nda*ndanum
!            X=DSQRT(X)
            call dadeb(31,'ERR DAALL ',1)
          endif
!
          if(nv.eq.0.or.nomax.eq.1) then
            call daclr(ic(i))
            idall(ic(i)) = idalm(ic(i))
          endif
        endif
      enddo
!
      if(nda.gt.ndamaxi) ndamaxi=nda

      return
      end
      subroutine daall(ic,l,ccc,no,nv)
      implicit none
      integer i,ind,l,ndanum,no,nv
      double precision x
!     ********************************
!
!     THIS SUBROUTINE ALLOCATES STORAGE FOR A DA VECTOR WITH
!     ORDER NO AND NUMBER OF VARIABLES NV
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"

      logical allvec(lda)
      integer nhole
      common /hole/nhole
      common /alloc/ allvec
      integer ndat
      common /alloctot/ ndat

!-----------------------------------------------------------------------------9
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      integer ic(*)
      logical incnda
      character c*10,ccc*10
!
      ind = 1

      do i=1,l
        if(ic(i).gt.0.and.ic(i).le.nda) then
!         DANAME(IC(I)) = C
!         IF(IDANO(IC(I)).EQ.NO.AND.IDANV(IC(I)).EQ.NV) THEN
        else
          if(nv.ne.0.and.(no.gt.nomax.or.nv.gt.nvmax)) then
            write(6,*)'ERROR IN DAALL, VECTOR ',c,' HAS NO, NV = ',     &
     &      no,nv,' NOMAX, NVMAX = ',nomax,nvmax
            call dadeb(31,'ERR DAALL ',1)
          endif
!
          if(nhole.gt.0) then
            ind=nda
20          if (allvec(ind)) then
              ind = ind - 1
              goto 20
            endif
            incnda = .false.
            nhole=nhole-1
          else
            incnda = .true.
            nda = nda + 1
            ind=nda
            if(nda.gt.lda) then
              write(6,*)'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHA',&
     &        'USTED'
              call dadeb(31,'ERR DAALL ',1)
            endif
          endif

          allvec(ind) = .true.

          ic(i) = ind
!
          if(nv.ne.0) then
            call danum(no,nv,ndanum)
          else
            ndanum = no
          endif

          c = ccc
          if(l.ne.1) write(c(6:10),'(I5)') i

          daname(ind) = c

          if (incnda) then
            if(ind.gt.nomax+2) then
              idano(ind) = nomax
              idanv(ind) = nvmax
              idapo(ind) = nst + 1
              idalm(ind) = nmmax
              idall(ind) = 0
              nst = nst + nmmax
            else
              idano(ind) = no
              idanv(ind) = nv
              idapo(ind) = nst + 1
              idalm(ind) = ndanum
              idall(ind) = 0
              nst = nst + ndanum
            endif
          endif
!
          if(nst.gt.lst) then
            x=-1.d0
            write(6,*)'ERROR IN DAALL, STACK EXHAUSTED '
            write(6,*) ' NST,LST '
            write(6,*)  nst,lst
            write(6,*) ' NDA,NDANUM,NDA*NDANUM '
            write(6,*)  nda,ndanum,nda*ndanum
!            X=DSQRT(X)
            call dadeb(31,'ERR DAALL ',1)
          endif
!
!          IF(NV.EQ.0) THEN
          if(nv.eq.0.or.nomax.eq.1) then
            call daclr(ic(i))
            idall(ic(i)) = idalm(ic(i))
          endif
        endif
      enddo
!
      if(nda.gt.ndamaxi) ndamaxi=nda

      return
      end
!
      subroutine dadal(idal,l)
      implicit none
      integer i,l
!     ************************
!
!     THIS SUBROUTINE DEALLOCATES THE VECTORS IDAL
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"

      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
      logical allvec(lda)
      integer nhole
      common /hole/nhole
      common /alloc/ allvec

      integer idal(*)

      do i=l,1,-1
        if(idal(i).le.nomax+2.or.idal(i).gt.nda) then
          write(6,*)'ERROR IN ROUTINE DADAL, IDAL(I),NDA = ',idal(i),nda
          call dadeb(31,'ERR DADAL ',1)
        endif
        if(idal(i).eq.nda) then
!       deallocate
          nst = idapo(nda) - 1
          nda = nda - 1
!        else
!        write(6,'(a10)')daname(i)
!        write(6,*)' etienne',idal(i),nda
!        write(6,*) sqrt(-1.d0)
        else
          nhole=nhole+1
        endif

        allvec(idal(i)) = .false.

!        IDANO(IDAL(I)) = 0
!        IDANV(IDAL(I)) = 0
!        IDAPO(IDAL(I)) = 0
!        IDALM(IDAL(I)) = 0
        idall(idal(i)) = 0

        idal(i) = 0
      enddo

      return
      end


      subroutine davar(ina,ckon,i)
      implicit none
      integer i,ibase,ic1,ic2,illa,ilma,ina,inoa,inva,ipoa
      double precision ckon
!     ****************************
!
!     THIS SUBROUTINE DECLARES THE DA VECTOR
!     AS THE INDEPENDENT VARIABLE NUMBER I.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!
!
      if(i.gt.inva) then
        write(6,*)'ERROR IN DAVAR, I = ',i,' EXCEEDS INVA = ',inva
        call dadeb(31,'ERR DAVAR ',1)
      endif
!
      if(nomax.eq.1) then
        if(i.gt.inva) then
          print*,'ERROR IN DAVAR, I = ',i,' EXCEEDS INVA = ',inva
!           CALL DADEB(31,'ERR DAVAR3',1)
        endif
        call daclr(ina)
        cc(ipoa) = ckon
        cc(ipoa+i) = 1d0
        return
      endif

      ibase = nomax+1
!
      if(i.gt.(nvmax+1)/2) then
        ic1 = 0
        ic2 = ibase**(i-(nvmax+1)/2-1)
      else
        ic1 = ibase**(i-1)
        ic2 = 0
      endif
!
      if(dabs(ckon).gt.eps) then
        idall(ina) = 2
        cc(ipoa) = ckon
        i1(ipoa) = 0
        i2(ipoa) = 0
!
        cc(ipoa+1) = 1.d0
        i1(ipoa+1) = ic1
        i2(ipoa+1) = ic2
      else
        idall(ina) = 1
        cc(ipoa) = 1.d0
        i1(ipoa) = ic1
        i2(ipoa) = ic2
      endif
!
      return
      end
!
      subroutine dacon(ina,ckon)
      implicit none
      integer illa,ilma,ina,inoa,inva,ipoa
      double precision ckon
!     **************************
!
!     THIS SUBROUTINE SETS THE VECTOR C TO THE CONSTANT CKON
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!
!
      if(nomax.eq.1) then
        call daclr(ina)
        cc(ipoa) = ckon
        return
      endif

      idall(ina) = 1
      cc(ipoa) = ckon
      i1(ipoa) = 0
      i2(ipoa) = 0
      if(dabs(ckon).lt.eps) idall(ina) = 0
!
      return
      end
!
      subroutine danot(not)
      implicit none
      integer not
!     *********************
!
!     THIS SUBROUTINE RESETS THE TRUNCATION ORDER NOCUT TO A NEW VALUE
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      if(not.gt.nomax) then
        write(6,*)'ERROR, NOCUT = ',not,' EXCEEDS NOMAX = ',nomax
        call dadeb(31,'ERR DANOT ',1)
      endif
!
      nocut = not
!
      return
      end

      subroutine getdanot(not)
      implicit none
      integer not
!     *********************
!
!     THIS SUBROUTINE RESETS THE TRUNCATION ORDER NOCUT TO A NEW VALUE
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      if(not.gt.nomax) then
        write(6,*)'ERROR, NOCUT = ',not,' EXCEEDS NOMAX = ',nomax
        call dadeb(31,'ERR DANOT ',1)
      endif
!
      not=nocut
!
      return
      end

      subroutine daeps(deps)
      implicit none
      double precision deps
!     **********************
!
!     THIS SUBROUTINE RESETS THE TRUNCATION ORDER NOCUT TO A NEW VALUE
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      if(deps.ge.0.d0) then
        eps = deps
      else
        deps=eps
      endif
!
      return
      end
!
      subroutine dapek(ina,jj,cjj)
      implicit none
      integer i,ibase,ic,ic1,ic2,icu,icz,ii1,ikk,illa,ilma,ina,         &
     &inoa,inva,ipek,ipoa,iu,iz,jj,jj1,mchk
      double precision cjj
!     ****************************
!
!     THIS SUBROUTINE DETERMINES THE COEFFICIENT OF THE ARRAY
!     OF EXPONENTS JJ AND RETURNS IT IN CJJ
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension jj(lnv)
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!
!
      jj1 = 1
      if(inva.eq.0.or.nomax.eq.1) then
        if(inva.ne.0.and.nomax.eq.1) then
          if(illa.ge.2) then
            do i=1,illa - 1
              if(jj(i).eq.1) jj1 = i + 1
            enddo
          else
            jj1 = jj(1) + 1
          endif
        else
          jj1 = jj(1)
        endif
        if(jj1.lt.1.or.jj1.gt.illa) then
          print*,'ERROR IN DAPEK, INDEX OUTSIDE RANGE, JJ(1) = ',jj1
!           CALL DADEB(31,'ERR DAPEK1',1)
        endif
        ipek = ipoa + jj1 - 1
        cjj = cc(ipek)
        return
      endif

      ii1 = (nvmax+1)/2
      ibase = nomax+1
!
!     DETERMINE INDEX TO BE SEARCHED FOR
!     **********************************
!
      call dadcd(jj,ic1,ic2)
!
!ETIENNE
      if(ic1.gt.lia.or.ic2.gt.lia) then
        write(6,*) 'DISASTER IN DAPEK, INA= ',ina
        write(32,*) 'DISASTER IN DAPEK, INA= ',ina
        write(32,*) ic1,ic2
        write(32,*) (jj(ikk),ikk=1,lnv)
      endif
!ETIENNE
      ic = ia1(ic1) + ia2(ic2)
!
!     DETERMINE IF MONOMIAL TO BE POKED CONFORMS WITH INOA, INVA,NOCUT
!     ****************************************************************
!
!      IF(ICO.GT.INOA.OR.ICV.GT.INVA.OR.ICO.GT.NOCUT) THEN
!         CJJ = 0
!         RETURN
!      ENDIF
!
!     DETERMINE IF MONOMIAL IS INSIDE FIRST AND LAST MONOMIALS OF A
!     *************************************************************
!
      iu = ipoa
      iz = ipoa + illa - 1
      icu = ia1(i1(iu))+ia2(i2(iu))
      icz = ia1(i1(iz))+ia2(i2(iz))
!
      if(illa.eq.0) then
        cjj = 0
        return
      elseif(ic.eq.icu) then
        cjj = cc(iu)
        return
      elseif(ic.eq.icz) then
        cjj = cc(iz)
        return
      elseif(ic.lt.icu.or.ic.gt.icz) then
        cjj = 0
        return
      endif
!
!     SEARCHING PROPER MONOMIAL
!     *************************
!
 10   continue
      if(iz-iu.le.1) then
        cjj = 0
        return
      endif
      i = (iu+iz)/2
!
!     if(ia1(i1(i))+ia2(i2(i)) - ic) 20,30,40
      mchk=ia1(i1(i))+ia2(i2(i)) - ic
      if(mchk.lt.0) goto 20
      if(mchk.eq.0) goto 30
      if(mchk.gt.0) goto 40
 20   iu = i
      goto 10
 30   cjj = cc(i)
      return
 40   iz = i
      goto 10
!
      end
!
      subroutine dapok(ina,jj,cjj)
      implicit none
      integer i,ic,ic1,ic2,icu,icz,ii,illa,ilma,ina,inoa,inva,          &
     &ipoa,ipok,iu,iz,jj,jj1,mchk
      double precision cjj
!     ****************************
!
!     THIS SUBROUTINE SETS THE COEFFICIENT OF THE ARRAY
!     OF EXPONENTS JJ TO THE VALUE CJJ
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension jj(lnv)
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!
!
      jj1 = 1
      if(inva.eq.0.or.nomax.eq.1) then
        if(inva.ne.0.and.nomax.eq.1) then
          if(illa.ge.2) then
            do i=1,illa - 1
              if(jj(i).eq.1) jj1 = i + 1
            enddo
          else
            jj1 = jj(1) + 1
          endif
        else
          jj1 = jj(1)
        endif
        if(jj1.lt.1.or.jj1.gt.illa) then
          print*,'ERROR IN DAPOK, INDEX OUTSIDE RANGE, JJ(1) = ',jj1
!           CALL DADEB(31,'ERR DAPOK1',1)
        endif
        ipok = ipoa + jj1 - 1
        cc(ipok) = cjj
        return
      endif

!     DETERMINE INDEX TO BE SEARCHED FOR
!     **********************************
!
      call dadcd(jj,ic1,ic2)
!
      ic = ia1(ic1) + ia2(ic2)
!
!     DETERMINE IF MONOMIAL TO BE POKED CONFORMS WITH INOA, INVA,NOCUT
!     ****************************************************************
!
!      IF(ICO.GT.INOA.OR.ICV.GT.INVA) THEN
!         write(6,*)'ERROR IN DAPOK, MONOMIAL NOT ALLOWED FOR ',A
!         CALL DADEB(31,'ERR DAPOK ',1)
!      ENDIF
!      IF(ICO.GT.NOCUT) RETURN
!
      iu = ipoa
      iz = ipoa + illa - 1
!
!     DETERMINE IF MONOMIAL IS INSIDE FIRST AND LAST MONOMIALS OF A
!     *************************************************************
!
      icu = ia1(i1(iu))+ia2(i2(iu))
      icz = ia1(i1(iz))+ia2(i2(iz))
      if(illa.eq.0) then
        i = ipoa
        goto 100
      elseif(ic.eq.icu) then
        cc(iu) = cjj
        i = iu
        goto 200
      elseif(ic.eq.icz) then
        cc(iz) = cjj
        i = iz
        goto 200
      elseif(ic.lt.icu) then
        i = iu
        goto 100
      elseif(ic.gt.icz) then
        i = iz + 1
        goto 100
      endif
!
!
!     SEARCHING PLACE TO POKE INTO OR BEFORE WHICH TO POKE
!     ****************************************************
!
      iu = ipoa
      iz = ipoa + illa
!
 10   continue
      if(iz-iu.le.1) then
        i = iz
        goto 100
      endif
      i = (iu+iz)/2
!
!      if(ia1(i1(i))+ia2(i2(i)) - ic) 20,30,40
      mchk=ia1(i1(i))+ia2(i2(i)) - ic
      if(mchk.lt.0) goto 20
      if(mchk.eq.0) goto 30
      if(mchk.gt.0) goto 40
 20   iu = i
      goto 10
 30   cc(i) = cjj
      goto 200
 40   iz = i
      goto 10
!
!     INSERTING THE MONOMIAL, MOVING THE REST
!     ***************************************
!
 100  continue
!
      if(dabs(cjj).lt.eps) return
!
      do ii=ipoa+illa,i+1,-1
        cc(ii) = cc(ii-1)
        i1(ii) = i1(ii-1)
        i2(ii) = i2(ii-1)
      enddo
!
      cc(i) = cjj
      i1(i) = ic1
      i2(i) = ic2
!
      idall(ina) = illa + 1
      if(idall(ina).gt.idalm(ina)) then
        write(6,*)'ERROR IN DAPOK '
        call dadeb(31,'ERR DAPOK ',1)
      endif
!
      return
!
!     CASE OF CJJ = 0 WHICH MEANS MOVING THE REST
!     *********************************************
!
 200  continue
      if(dabs(cjj).lt.eps) then
        do ii=i,ipoa+illa-2
          cc(ii) = cc(ii+1)
          i1(ii) = i1(ii+1)
          i2(ii) = i2(ii+1)
        enddo
        idall(ina) = illa - 1
      endif
      return
!
      end
!
      subroutine daclr(inc)
      implicit none
      integer i,illc,ilmc,inc,inoc,invc,ipoc
!     *********************
!
!     THIS SUBROUTINE SETS ALL THE STACK SPACE RESERVED FOR VARIABLE
!     C TO ZERO
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
      do i=ipoc,ipoc+ilmc-1
!
        cc(i) = 0.d0
!
      enddo
!
      return
      end
!
      subroutine dacop(ina,inb)
      implicit none
      integer ia,ib,iif,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,&
     &ipoa,ipob
!     *************************
!
!     THIS SUBROUTINE COPIES THE DA VECTOR A TO THE DA VECTOR B
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INB,INOB,INVB)
!
      ipob = idapo(inb)
      ipoa = idapo(ina)
      illa = idall(ina)
      ib = ipob - 1
!
!      iif = 0
!      if(nomax.eq.1.or.inva.eq.0) iif = 1

      do ia = ipoa,ipoa+illa-1
!
        if(nomax.gt.1) then
          if(ieo(ia1(i1(ia))+ia2(i2(ia))).gt.nocut) goto 100
        endif
        ib = ib + 1
        cc(ib) = cc(ia)
        i1(ib) = i1(ia)
        i2(ib) = i2(ia)
!
 100    continue
      enddo
!
      idall(inb) = ib - ipob + 1
!      if(idall(inb).gt.idalm(inb)) then
!         write(6,*)'ERROR IN DACOP'
!         call dadeb(31,'ERR DACOP ',1)
!      endif
!
      return
      end

!
      subroutine datrashn(idif,ina,inbb)
      implicit none
      integer i,ia,idif,illa,ilma,ina,inb,inbb,inoa,inva,ipoa
      double precision rr
!     *************************
!
!     THIS SUBROUTINE COPIES THE DA VECTOR A TO THE DA VECTOR B
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"

      integer jd(lnv)
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      inb=0

      if(inbb.eq.ina) then
        call daall(inb,1,'$$DAADD $$',inoa,inva)
      else
        inb=inbb
      endif

      call daclr(inb)
!
!
      do ia = ipoa,ipoa+illa-1
!
        if(nomax.ne.1) then
          call dancd(i1(ia),i2(ia),jd)
        else
          do i=1,lnv
            jd(i)=0
          enddo
          if(ia.ne.ipoa) then
            jd(ia-ipoa+1)=1
          endif
        endif

        call dapek(ina,jd,rr)
        jd(idif)=0
        if(dabs(rr).gt.0.d0) call dapok(inb,jd,rr)
!
      enddo
!
!
      if(inbb.eq.ina) then
        call dacop(inb,inbb)
        call dadal(inb,1)
      endif

      return
      end
!
      subroutine daadd(ina,inb,inc)
      implicit none
      integer i,ic,ic1,ic2,icu,icz,ii,illa,ilma,ina,inoa,inva,          &
     &ipoa,ipok,iu,iz,jj,jj1
      integer idaadd,illc,ilmc,inb,inc,inoc,invc,ipoc

      include "TPSALib_prm.f"

      integer ipob
!     *****************************
!
!     THIS SUBROUTINE PERFORMS A DA ADDITION OF THE DA VECTORS A AND B.
!     THE RESULT IS STORED IN C.
!
      if(nomax.eq.1) then
        ipoc = idapo(inc)
        ipoa = idapo(ina)
        ipob = idapo(inb)
!         minv = min(inva,invb,invc)
        do i=0,nvmax
          cc(ipoc+i) = cc(ipoa+i)   + cc(ipob+i)
        enddo
        return
      endif

      if(ina.ne.inc.and.inb.ne.inc) then
        call dalin(ina,+1.d0,inb,+1.d0,inc)
      else
        idaadd = 0
        call daall(idaadd,1,'$$DAADD $$',nomax,nvmax)
        call dalin(ina,+1.d0,inb,+1.d0,idaadd)
        call dacop(idaadd,inc)
        call dadal(idaadd,1)
      endif
!
      return
      end
!
      subroutine dasub(ina,inb,inc)
      implicit none
      integer idasub
      integer i,ic,ic1,ic2,icu,icz,ii,illa,ilma,inoa,inva,ina,          &
     &ipoa,ipok,iu,iz,jj,jj1
      integer idaadd,illc,ilmc,inb,inc,inoc,invc,ipoc

      include "TPSALib_prm.f"

      integer ipob
!     THIS SUBROUTINE PERFORMS A DA SUBTRACTION OF THE DA VECTORS A AND B.
!     THE RESULT IS STORED IN C.
!
      if(nomax.eq.1) then
        ipoc = idapo(inc)
        ipoa = idapo(ina)
        ipob = idapo(inb)
!         minv = min(inva,invb,invc)
        do i=0,nvmax
          cc(ipoc+i) = cc(ipoa+i)   - cc(ipob+i)
        enddo
        return
      endif
      if(ina.ne.inc.and.inb.ne.inc) then
        call dalin(ina,+1.d0,inb,-1.d0,inc)
      else
        idasub = -1
!         call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        call daall(idasub,1,'$$DASUB $$',nomax,nvmax)
        call dalin(ina,+1.d0,inb,-1.d0,idasub)
        call dacop(idasub,inc)
        call dadal(idasub,1)
      endif
!
      return
      end
!
      subroutine damulin(ina,inb,coe1,inc,ind,coe2,ine)
      implicit none
      integer ina,inb,inc,incc,ind,ine,inoc,invc
      double precision coe1,coe2
!     *****************************
!
!     THIS SUBROUTINE PERFORMS A DA MULTIPLICATION OF THE DA VECTORS A AND B.
!     THE RESULT IS STORED IN C. AS TEMPORARY STORAGE, THE STACK SPACE
!     OF THE (NOMAX+2) SCRATCH VARIABLES ALLOCATED BY DAINI IS USED.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!

      call daall(incc,1,'$$DAJUNK$$',inoc,invc)
      call damul(ina,inb,incc)
      call damul(inc,ind,ine)
      call dalin(incc,coe1,ine,coe2,ine )
      call dadal(incc,1)

      return
      end

!
!
! ANFANG UNTERPROGRAMM
      subroutine daexx(ina,inb,inc)
      implicit none
      integer illc,ilmc,ina,inaa,inb,inbb,inc,inoc,invc,ipoc

      include "TPSALib_prm.f"

!     ******************************
!
!     THIS SUBROUTINE EXPONENTIATES INE WITH THE CONSTANT CKON
!
!-----------------------------------------------------------------------------1
      write(6,*) "daexx"
      if(ina.ne.inc.and.inb.ne.inc) then
        call daexxt(ina,inb,inc)
      else
        inaa = 0
        inbb = 0
!         call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        call daall(inaa,1,'$$DAADD $$',nomax,nvmax)
        call daall(inbb,1,'$$DAADD $$',nomax,nvmax)
        call dacop(ina,inaa)
        call dacop(inb,inbb)
        call daexxt(inaa,inbb,inc)
        call dadal(inaa,1)
        call dadal(inbb,1)
      endif

      return
      end

! ANFANG UNTERPROGRAMM
      subroutine daexxt(ina,inb,inc)
      implicit none
      integer idaexx,illa,illb,illc,ilma,ilmb,ilmc,ina,inb,inc,inoa,    &
     &inob,inoc,inva,invb,invc,ipoa,ipob,ipoc
!     ******************************
!
!     THIS SUBROUTINE EXPONENTIATES INA WITH INB
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
!      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!      call dainf(inb,inob,invb,ipob,ilmb,illb)
!      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
      idaexx = 0
      call daall(idaexx,1,'$$DAEXX $$',nomax,nvmax)
      call dafun('LOG   ',ina,inc)
      call damul(inb,inc,idaexx)
      call dafun('EXP   ',idaexx,inc)
      call dadal(idaexx,1)
!
      return
      end

! ANFANG UNTERPROGRAMM
      subroutine dacex(ina,ckon,inb)
      implicit none
      integer illc,ilmc,ina,inb,inc,incc,inoc,invc,ipoc
      double precision ckon

      include "TPSALib_prm.f"

!     ******************************
!
!     THIS SUBROUTINE EXPONENTIATES INE WITH THE CONSTANT CKON
!
!-----------------------------------------------------------------------------1
      write(6,*) "dacex"
      if(ina.eq.inb) then
!        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',nomax,nvmax)
        call dacext(ina,ckon,incc)
        call dacop(incc,inb)
        call dadal(incc,1)
      else
        call dacext(ina,ckon,inb)
      endif

      return
      end
      subroutine dacext(ina,ckon,inb)
      implicit none
      integer idacex,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,   &
     &ipoa,ipob
      double precision ckon
!     ******************************
!
!     THIS SUBROUTINE EXPONENTIATES THE CONSTANT CKON WITH INA
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
!      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
      if(ckon.le.0) then
        print*,'ERROR IN DACEX, CKON NOT POSITIVE'
!        CALL DADEB(31,'ERR DACEX1',1)
      endif
!
      idacex = 0
      call daall(idacex,1,'$$DACEX $$',nomax,nvmax)
      ckon = dlog(ckon)
      call dacmu(ina,ckon,idacex)
      call dafun('EXP   ',idacex,inb)
      call dadal(idacex,1)
!
      return
      end

      subroutine daexc(ina,ckon,inb)
      implicit none
      integer illc,ilmc,ina,inb,inc,incc,inoc,invc,ipoc
      double precision ckon

      include "TPSALib_prm.f"

!     ******************************
!
!     THIS SUBROUTINE EXPONENTIATES INE WITH THE CONSTANT CKON
!
!-----------------------------------------------------------------------------1
!        write(6,*) "daexc"

      if(ina.eq.inb) then
!        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',nomax,nvmax)
        call daexct(ina,ckon,incc)
        call dacop(incc,inb)
        call dadal(incc,1)
      else
        call daexct(ina,ckon,inb)
      endif

      return
      end

      subroutine daexct(ina,ckon,inb)
      implicit none
      integer i,ic,idaexc,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,   &
     &invb,ipoa,ipob
      double precision ckon,xic

      include "TPSALib_prm.f"

!     ******************************
!
!     THIS SUBROUTINE EXPONENTIATES INE WITH THE CONSTANT CKON
!
!-----------------------------------------------------------------------------1
!      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
      idaexc = 0
      call daall(idaexc,1,'$$DAEXC $$',nomax,nvmax)

      if(ckon.lt.0.d0) then
        call dafun('LOG   ',ina,inb)
        call dacmu(inb,ckon,idaexc)
        call dafun('EXP   ',idaexc,inb)
      else
        xic=dabs(ckon-dble(idint(ckon)))
        if(xic.gt.eps) then
          call dafun('LOG   ',ina,inb)
          call dacmu(inb,ckon,idaexc)
          call dafun('EXP   ',idaexc,inb)
        else
          ic=idint(ckon)
          call dacon(idaexc,1.d0)
          do i=1,ic
            call damul(idaexc,ina,idaexc)
          enddo
          call dacop(idaexc,inb)
        endif
      endif
      call dadal(idaexc,1)
!
      return
      end

      subroutine damul(ina,inb,inc)
      implicit none
      integer illc,ilmc,ina,inb,inc,incc,inoc,invc,ipoc
!     *****************************
!
!     THIS SUBROUTINE PERFORMS A DA MULTIPLICATION OF THE DA VECTORS A AND B.
!     THE RESULT IS STORED IN C. AS TEMPORARY STORAGE, THE STACK SPACE
!     OF THE (NOMAX+2) SCRATCH VARIABLES ALLOCATED BY DAINI IS USED.
!
!-----------------------------------------------------------------------------1

      include "TPSALib_prm.f"

      double precision ccipoa,ccipob
      integer ipoa,ipob,i
!
      if(nomax.eq.1) then
        ipoa=idapo(ina)
        ipob=idapo(inb)
        ipoc=idapo(inc)
!         minv = min(inva,invb,invc)
        ccipoa = cc(ipoa)
        ccipob = cc(ipob)
        cc(ipoc) = ccipoa*ccipob
        do i=1,nvmax
          cc(ipoc+i) = ccipoa*cc(ipob+i) + ccipob*cc(ipoa+i)
        enddo
!         do 30 i=ipoc+minv+1,ipoc+invc
!  30     cc(i) = 0d0
        return
      endif

      if(ina.eq.inc.or.inb.eq.inc) then
!        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',nomax,nvmax)
        call damult(ina,inb,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call damult(ina,inb,inc)
      endif

      return
      end

      subroutine damult(ina,inb,inc)
      implicit none
      integer i,i1ia,i2ia,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,ina,   &
     &inb,inc,inoa,inob,inoc,inva,invb,invc,ioffb,ipno,ipoa,ipob,ipoc,  &
     &ipos,minv,noff,noib,nom
      double precision ccia,ccipoa,ccipob
!     *****************************
!
!     THIS SUBROUTINE PERFORMS A DA MULTIPLICATION OF THE DA VECTORS A AND B.
!     THE RESULT IS STORED IN C. AS TEMPORARY STORAGE, THE STACK SPACE
!     OF THE (NOMAX+2) SCRATCH VARIABLES ALLOCATED BY DAINI IS USED.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension ipno(0:lno),noff(0:lno)
!
!
!      CALL DACHK(INA,INOA,INVA, INB,INOB,INVB, INC,INOC,INVC)
!
!     CASE OF FIRST ORDER ONLY
!     ************************


      if(nomax.eq.1) then
        ipoa=idapo(ina)
        ipob=idapo(inb)
        ipoc=idapo(inc)
!         minv = min(inva,invb,invc)
        ccipoa = cc(ipoa)
        ccipob = cc(ipob)
        cc(ipoc) = ccipoa*ccipob
        do i=1,nvmax
          cc(ipoc+i) = ccipoa*cc(ipob+i) + ccipob*cc(ipoa+i)
        enddo
!         do 30 i=ipoc+minv+1,ipoc+invc
!  30     cc(i) = 0d0
        return
      endif
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inb,inob,invb,ipob,ilmb,illb)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)

!     GENERAL CASE
!     ************
!
      do i=0,nomax
        noff(i) = idapo(i+2)
        ipno(i) = 0
      enddo
!
      call daclr(1)
!
!     RE-SORTING THE VECTOR B INTO PIECES THAT ARE OF ONLY ONE ORDER
!     *************************************************************
!
      do ib=ipob,ipob+illb-1
!
        noib = ieo(ia1(i1(ib))+ia2(i2(ib)))
        ipos = ipno(noib) + 1
        ipno(noib) = ipos
        inob = noff(noib) + ipos
!
        cc(inob) = cc(ib)
        i1(inob) = i1(ib)
        i2(inob) = i2(ib)
!
      enddo
!
      do i=0,nomax
        idall(i+2) = ipno(i)
      enddo
!
!     PERFORMING ACTUAL MULTIPLICATION
!     ********************************
!
      nom = min(nocut,inoc)
!
      do ia=ipoa,ipoa+illa-1
!
        i1ia = i1(ia)
        i2ia = i2(ia)
        ccia = cc(ia)
!
        do noib = 0,nom-ieo(ia1(i1(ia))+ia2(i2(ia)))
!
          ioffb = noff(noib)
!
          do ib = ioffb+1,ioffb+ipno(noib)
!
            ic = ia2(i2ia+i2(ib)) + ia1(i1ia + i1(ib))
            cc(ic) = cc(ic) + ccia*cc(ib)
!
          enddo
        enddo
      enddo
!
      call dapac(inc)
!
      return
      end
!
      subroutine dadiv(ina,inb,inc)
      implicit none
      integer idadiv  ,inb
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
!     *************************
!
!     THIS SUBROUTINE SQUARES THE VECTOR A AND STORES THE RESULT IN C.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      double precision ck,ck1
      integer ipoa,ipob,i
!
!     THIS SUBROUTINE PERFORMS A DA DIVISION OF THE DA VECTORS A AND B.
!     THE RESULT IS STORED IN C.
!

      if(nomax.eq.1) then
!         minv = min(inva,invb)
        ipoa = idapo(ina)
        ipob = idapo(inb)
        ipoc = idapo(inc)
        ck=1.d0/cc(ipob)
        ck1=cc(ipoa)*ck
        do i=1,nvmax
          cc(ipoc+i) = (cc(ipoa+i)-cc(ipob+i)*ck1)*ck
        enddo
        cc(ipoc)=ck1
        return
      endif

      idadiv = 0
!      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
      call daall(idadiv,1,'$$DADIV $$',nomax,nvmax)
      call dafun('INV ',inb,idadiv)
      call damul(ina,idadiv,inc)
      call dadal(idadiv,1)
!
      return
      end

!
      subroutine dasqr(ina,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
!     *************************
!
!     THIS SUBROUTINE SQUARES THE VECTOR A AND STORES THE RESULT IN C.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer i,ipoa
      double precision ccipoa
!
!     CASE OF FIRST ORDER ONLY
!     ************************
      if(nomax.eq.1) then
        ipoc = idapo(inc)
        ipoa = idapo(ina)
!         minv = min(inva,invc)
        ccipoa = cc(ipoa)
        cc(ipoc) = ccipoa*ccipoa
        do i=1,nvmax
          cc(ipoc+i) = 2.d0*ccipoa*cc(ipoa+i)
        enddo
        return
      endif

      if(ina.eq.inc) then
!        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',nomax,nvmax)
        call dasqrt(ina,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dasqrt(ina,inc)
      endif

      return
      end

      subroutine dasqrt(ina,inc)
      implicit none
      integer i,i1ia,i2ia,ia,ib,ib1,ic,illa,illc,ilma,ilmc,ina,inc,inoa,&
     &inoc,inva,invc,ioffa,ioffb,ipno,ipoa,ipoc,ipos,                   &
     &minv,noff,noia,noib,nom
      double precision ccia,ccipoa
!     *************************
!
!     THIS SUBROUTINE SQUARES THE VECTOR A AND STORES THE RESULT IN C.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension ipno(0:lno),noff(0:lno)
!
!
!      CALL DACHK(INA,INOA,INVA,'          ',-1,-1,INC,INOC,INVC)
!
!
!     CASE OF FIRST ORDER ONLY
!     ************************
      if(nomax.eq.1) then
        ipoc = idapo(inc)
        ipoa = idapo(ina)
!         minv = min(inva,invc)
        ccipoa = cc(ipoa)
        cc(ipoc) = ccipoa*ccipoa
        do i=1,nvmax
          cc(ipoc+i) = 2d0*ccipoa*cc(ipoa+i)
        enddo
        return
      endif
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!     GENERAL CASE
!     ************
!
      do i=0,nomax
        noff(i) = idapo(i+2)
        ipno(i) = 0
      enddo
!
      call daclr(1)
!
!     RESORTING THE VECTOR A INTO PIECES THAT ARE OF ONLY ONE ORDER
!     *************************************************************
!
      do ia=ipoa,ipoa+illa-1
!
        noia = ieo(ia1(i1(ia))+ia2(i2(ia)))
        ipos = ipno(noia) + 1
        ipno(noia) = ipos
        inoa = noff(noia) + ipos
!
        cc(inoa) = cc(ia)
        i1(inoa) = i1(ia)
        i2(inoa) = i2(ia)
!
      enddo
!
      do i=0,nomax
        idall(i+2) = ipno(i)
      enddo
!
!     PERFORMING ACTUAL MULTIPLICATION
!     ********************************
!
      nom = min(nocut,inoc)
!
      do noia = 0,nom/2
!
        ioffa = noff(noia)
!
        do ia=ioffa+1,ioffa+ipno(noia)
!
          i1ia = i1(ia)
          i2ia = i2(ia)
          ccia = cc(ia)
!
          ic = ia2(i2ia+i2ia) + ia1(i1ia+i1ia)
          cc(ic) = cc(ic) + ccia*ccia
          ccia = ccia + ccia
!
          do noib = noia,nom-noia
!
            ioffb = noff(noib)
            if(noib.eq.noia) then
              ib1 = ia + 1
            else
              ib1 = ioffb + 1
            endif
!
            do ib = ib1,ioffb+ipno(noib)
!
              ic = ia2(i2ia+i2(ib)) + ia1(i1ia + i1(ib))
              cc(ic) = cc(ic) + ccia*cc(ib)
!
            enddo
          enddo
        enddo
      enddo
!
      call dapac(inc)
!
      return
      end
!
      subroutine dacad(ina,ckon,inb)
      implicit none
      integer illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,ipoa,ipob
      double precision ckon,const
!     ******************************
!
!     THIS SUBROUTINE ADDS THE CONSTANT CKON TO THE VECTOR A
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer jj(lnv)
      data jj / lnv*0 /
!
!
!
      call dacop(ina,inb)
      if(nomax.eq.1) then
        cc(idapo(inb)) = cc(idapo(inb)) + ckon
        return
      endif
!
      call dapek(inb,jj,const)
      call dapok(inb,jj,const+ckon)
!
      return
      end
!
      subroutine dacsu(ina,ckon,inb)
      implicit none
      integer illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,ipoa,ipob
      double precision ckon,const
!     ******************************
!
!     THIS SUBROUTINE SUBTRACTS THE CONSTANT CKON FROM THE VECTOR A
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer jj(lnv)
      data jj / lnv*0 /
!
!
      call dacop(ina,inb)
!
      if(nomax.eq.1) then
        cc(idapo(inb)) = cc(idapo(inb)) - ckon
        return
      endif
!
      call dapek(inb,jj,const)
      call dapok(inb,jj,const-ckon)
!
      return
      end
!
      subroutine dasuc(ina,ckon,inb)
      implicit none
      integer illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,ipoa,ipob
      double precision ckon
!     ******************************
!
!     THIS SUBROUTINE SUBTRACTS THE VECTOR INA FROM THE CONSTANT CKON
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer i
!
!      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
      ipob=idapo(inb)
      ipoa=idapo(ina)
      if(nomax.eq.1) then
        cc(ipob) = ckon - cc(ipoa)
        do i=1,nvmax
          cc(ipob+i) =-cc(ipoa+i)
        enddo
        return
      endif

      call dacsu(ina,ckon,inb)
      call dacmu(inb,-1.d0,inb)
!
      return
      end
!
      subroutine dacmu(ina,ckon,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
      double precision ckon
!     ******************************
!
!     THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
!     THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
!     THE DA VECTOR DENOTED WITH THE INTEGER E.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer ipoa,i
!
      if(nomax.eq.1) then
!         minv = min(inva,invb)
        ipoa = idapo(ina)
        ipoc = idapo(inc)
        do i=0,nvmax
          cc(ipoc+i) = cc(ipoa+i) * ckon
        enddo
        return
      endif

      if(ina.eq.inc) then
!        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daallno(incc,1,'$$DAJUNK$$')
        call dacmut(ina,ckon,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dacmut(ina,ckon,inc)
      endif

      return
      end

      subroutine dacmut(ina,ckon,inb)
      implicit none
      integer i,ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,  &
     &ipoa,ipob,minv
      double precision ckon
!     ******************************
!
!     THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
!     THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
!     THE DA VECTOR DENOTED WITH THE INTEGER E.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INB,INOB,INVB)
!
      if(nomax.eq.1) then
!         minv = min(inva,invb)
        ipoa = idapo(ina)
        ipob = idapo(inb)
        do i=0,nvmax
          cc(ipob+i) = cc(ipoa+i) * ckon
        enddo
        return
      endif
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inb,inob,invb,ipob,ilmb,illb)
      if(dabs(ckon).lt.eps) then
        idall(inb) = 0
        return
      endif
!
      ib = ipob - 1
!
      do ia=ipoa,ipoa+illa-1
!
        if(ieo(ia1(i1(ia))+ia2(i2(ia))).gt.nocut) goto 100
        ib = ib + 1
        cc(ib) = cc(ia)*ckon
        i1(ib) = i1(ia)
        i2(ib) = i2(ia)
!
 100    continue
      enddo
!
      idall(inb) = ib-ipob+1
      if(idall(inb).gt.idalm(inb)) then
        write(6,*)'ERROR IN DACMU '
        call dadeb(31,'ERR DACMU ',1)
      endif
!
      return
      end
!
      subroutine dacdi(ina,ckon,inb)
      implicit none
      integer illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,ipoa,ipob
      double precision ckon
!     ******************************
!
!     THIS SUBROUTINE DIVIDES THE VECTOR INA BY THE CONSTANT CKON
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer i
!
!      if(ckon.eq.0.d0) then
!         write(6,*)'ERROR IN DACDI, CKON IS ZERO'
!         call dadeb(31,'ERR DACDI ',1)
!      endif
!
!      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
!
      if(nomax.eq.1) then
!         minv = min(inva,invb)
        ipoa = idapo(ina)
        ipob = idapo(inb)
        do i=0,nvmax
          cc(ipob+i) = cc(ipoa+i)/ ckon
        enddo
        return
      endif

      call dacmu(ina,1.d0/ckon,inb)
!
      return
      end
!
!
      subroutine dadic(ina,ckon,inc)
      implicit none
      integer idadic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,   &
     &ipoa,ipoc
      double precision ckon,zero
!     ******************************
!
!     THIS SUBROUTINE DIVIDES THE CONSTANT CKON BY THE VECTOR INA
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      double precision ck
      integer i
!
!
!
      zero=0.d0
      if(nomax.eq.1) then
!         minv = min(inva,invb)
        ipoa = idapo(ina)
        ipoc = idapo(inc)
        cc(ipoc)=ckon/cc(ipoa)
        ck=cc(ipoc)/cc(ipoa)
        do i=1,nvmax
          cc(ipoc+i) = -cc(ipoa+i)* ck
        enddo
        return
      endif


      if(abs(ckon).lt.eps) then
        call dacon(inc,zero)
        return
      endif

      idadic = 0
      call daall(idadic,1,'$$DADIC $$',nomax,nvmax)

      if(ckon.eq.0.d0) then
        write(6,*)'ERROR IN DACDI and DADIC, CKON IS ZERO'
        call dadeb(31,'ERR DACDI ',1)
      endif
      call dacdi(ina,ckon,idadic)
      call dafun('INV ',idadic,inc)
      call dadal(idadic,1)
!
      return
      end
!
      subroutine dacma(ina,inb,bfac,inc)
      implicit none
      integer idacma,illc,ilmc,ina,inb,inc,inoc,invc,ipoc
      double precision bfac
!     **********************************
!
!     THIS SUBROUTINE PERFORMS THE OPERATIONS C = A + B*BFAC, WHERE A,B,C ARE
!     DA VECTORS AND BFAC IS A DOUBLE PRECISION. A AND C CAN BE IDENTICAL.
!     CAN LATER BE REPLACED BY SOMETHING LIKE DAADD WITH MINOR CHANGES.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer  ipob,ipoa,i
!
      if(nomax.eq.1) then
        ipoc = idapo(inc)
        ipoa = idapo(ina)
        ipob = idapo(inb)
!         minv = min(inva,invb,invc)
        do i=0,nvmax
          cc(ipoc+i) = cc(ipoa+i)   + cc(ipob+i) * bfac
        enddo
!         do 8 i=ipoc+minv+1,ipoc+invc
! 8       cc(i) = 0d0
        return
      endif
!      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
      idacma = 0
      call daall(idacma,1,'$$DACMA $$',nomax,nvmax)
      call dalin(ina,+1.d0,inb,bfac,idacma)
      call dacop(idacma,inc)
      call dadal(idacma,1)
!
      return
      end
!
      subroutine dalin(ina,afac,inb,bfac,inc)
      implicit none
      integer illc,ilmc,ina,inb,inc,incc,inoc,invc,ipoc
      double precision afac,bfac
!     ***************************************
!
!     THIS SUBROUTINE COMPUTES THE LINEAR COMBINATION
!     C = AFAC*A + BFAC*B. IT IS ALSO USED TO ADD AND SUBTRACT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer  ipob,ipoa,i
!
      if(nomax.eq.1) then
        ipoc = idapo(inc)
        ipoa = idapo(ina)
        ipob = idapo(inb)
!         minv = min(inva,invb,invc)
        do i=0,nvmax
          cc(ipoc+i) = cc(ipoa+i) * afac + cc(ipob+i) * bfac
        enddo
!         do 8 i=ipoc+minv+1,ipoc+invc
! 8       cc(i) = 0d0
        return
      endif

      if(ina.eq.inc.or.inb.eq.inc) then
!        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',nomax,nvmax)
        call dalint(ina,afac,inb,bfac,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dalint(ina,afac,inb,bfac,inc)
      endif

      return
      end


      subroutine dalint(ina,afac,inb,bfac,inc)
      implicit none
      integer i,ia,iamax,ib,ibmax,ic,icmax,illa,illb,illc,ilma,ilmb,    &
     &ilmc,ina,inb,inc,inoa,inob,inoc,inva,invb,invc,ipoa,ipob,ipoc,is, &
     &ismax,ismin,ja,jb,minv,mchk
      double precision afac,bfac,ccc,copf
!     ***************************************
!
!     THIS SUBROUTINE COMPUTES THE LINEAR COMBINATION
!     C = AFAC*A + BFAC*B. IT IS ALSO USED TO ADD AND SUBTRACT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!
!      CALL DACHK(INA,INOA,INVA, INB,INOB,INVB, INC,INOC,INVC)
!
      if(nomax.eq.1) then
        ipoc = idapo(inc)
        ipoa = idapo(ina)
        ipob = idapo(inb)
!         minv = min(inva,invb,invc)
        do i=0,nvmax
          cc(ipoc+i) = cc(ipoa+i) * afac + cc(ipob+i) * bfac
        enddo
!         do 8 i=ipoc+minv+1,ipoc+invc
! 8       cc(i) = 0d0
        return
      endif
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inb,inob,invb,ipob,ilmb,illb)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
      ia = ipoa
      ib = ipob
      ic = ipoc - 1
      iamax = ipoa+illa-1
      ibmax = ipob+illb-1
      icmax = ipoc+ilmc-1
      ja = ia1(i1(ia)) + ia2(i2(ia))
      jb = ia1(i1(ib)) + ia2(i2(ib))
!
      if(ia.gt.iamax) then
        ismin = ib
        ismax = ibmax
        copf  = bfac
        goto 50
      endif
      if(ib.gt.ibmax) then
        ismin = ia
        ismax = iamax
        copf  = afac
        goto 50
      endif
!
!     COMPARING
!     *********
!
  10  continue
!      if(ja-jb) 30,20,40
      mchk=ja-jb
      if(mchk.lt.0) goto 30
      if(mchk.eq.0) goto 20
      if(mchk.gt.0) goto 40
!
!     ADDING TWO TERMS
!     ****************
!
  20  continue
      ccc = cc(ia)*afac + cc(ib)*bfac
      if(dabs(ccc).lt.eps) goto 25
      if(ieo(ia1(i1(ia))+ia2(i2(ia))).gt.nocut) goto 25
      ic = ic + 1
      cc(ic) = ccc
      i1(ic) = i1(ia)
      i2(ic) = i2(ia)
  25  continue
      ia = ia + 1
      ib = ib + 1
      if(ia.gt.iamax) then
        ismin = ib
        ismax = ibmax
        copf  = bfac
        goto 50
      endif
      if(ib.gt.ibmax) then
        ismin = ia
        ismax = iamax
        copf  = afac
        goto 50
      endif
      ja = ia1(i1(ia)) + ia2(i2(ia))
      jb = ia1(i1(ib)) + ia2(i2(ib))
      goto 10
!
!     STORING TERM A
!     **************
!
  30  continue
      if(ieo(ia1(i1(ia))+ia2(i2(ia))).gt.nocut) goto 35
      ccc = cc(ia)*afac
      if(dabs(ccc).lt.eps) goto 35
      ic = ic + 1
      cc(ic) = ccc
      i1(ic) = i1(ia)
      i2(ic) = i2(ia)
  35  continue
      ia = ia + 1
      if(ia.gt.iamax) then
        ismin = ib
        ismax = ibmax
        copf  = bfac
        goto 50
      endif
      ja = ia1(i1(ia)) + ia2(i2(ia))
      goto 10
!
!     STORING TERM B
!     **************
!
  40  continue
      if(ieo(ia1(i1(ib))+ia2(i2(ib))).gt.nocut) goto 45
      ccc = cc(ib)*bfac
      if(dabs(ccc).lt.eps) goto 45
      ic = ic + 1
      cc(ic) = ccc
      i1(ic) = i1(ib)
      i2(ic) = i2(ib)
  45  continue
      ib = ib + 1
      if(ib.gt.ibmax) then
        ismin = ia
        ismax = iamax
        copf  = afac
        goto 50
      endif
      jb = ia1(i1(ib)) + ia2(i2(ib))
      goto 10
!
!     COPYING THE REST
!     ****************
!
  50  continue
      do is=ismin,ismax
        if(ieo(ia1(i1(is))+ia2(i2(is))).gt.nocut) goto 60
        ccc = cc(is)*copf
        if(dabs(ccc).lt.eps) goto 60
        ic = ic + 1
        cc(ic) = ccc
        i1(ic) = i1(is)
        i2(ic) = i2(is)
  60    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
!
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DALIN, RESULT HAS TOO MANY TERMS '
        call dadeb(31,'ERR DALIN ',1)
      endif
!
      return
      end
!
      subroutine dafun(cf,ina,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
!     ****************************
!
!     THIS SUBROUTINE COMPUTES THE FUNCTION CF OF THE DA VECTOR A
!     AND STORES THE RESULT IN C.
!     AT PRESENT, SOME FUNCTIONS CAN BE COMPUTED ONLY TO FIFTH ORDER.
!     THIS HAS TO BE FIXED IN THE FUTURE.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      character cf*4

      if(ina.eq.inc) then
!       call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',nomax,nvmax)
        call dafunt(cf,ina,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dafunt(cf,ina,inc)
      endif

      return
      end

      subroutine dafunt(cf,ina,inc)
      implicit none
      integer i,illa,illc,ilma,ilmc,ina,inc,ind,inoa,inoc,inon,inva,    &
     &invc,ipoa,ipoc,ipow,iscr,jj,lfun,no
      double precision a0,a1,a2,a3,a4,a5,ca,e1,                         &
     &e2,ea,era,p,ra,rpi4,sa,scr,                                       &
     &t,xf
!     ****************************
!
!     THIS SUBROUTINE COMPUTES THE FUNCTION CF OF THE DA VECTOR A
!     AND STORES THE RESULT IN C.
!     AT PRESENT, SOME FUNCTIONS CAN BE COMPUTED ONLY TO FIFTH ORDER.
!     THIS HAS TO BE FIXED IN THE FUTURE.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      character cf*4,cfh*4,abcs*26,abcc*26
      dimension xf(0:lno),jj(lnv)
!
      data jj /lnv*0/
      data abcs /'abcdefghijklmnopqrstuvwxyz'/
      data abcc /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
!
      if(cf(1:1).eq.' ') then
        cfh(1:3) = cf(2:4)
        cfh(1:4) = ' '
        cf = cfh
      endif
!
      do i=1,4
        ind = index(abcs,cf(i:i))
        if(ind.ne.0) cf(i:i) = abcc(ind:ind)
      enddo
!      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
!     CASE OF NV = 0 WHICH MEANS COORDINATEWISE OPERATION
!     ***************************************************
!
!     CASE OF NV > 0 WHICH MEANS DIFFERENTIAL ALGEBRAIC OPERATION
!     ***********************************************************
!
      if(cf.eq.'SQR ') then
        call dasqr(ina,inc)
        return
      endif
!
!     ALLOCATE VARIABLES, PICK ZEROTH ORDER TERM
!     ******************************************
!
      ipow = 0
      inon = 0
      iscr = 0
!
      call daall(ipow,1,'$$DAFUN1$$',nomax,nvmax)
      call daall(inon,1,'$$DAFUN2$$',nomax,nvmax)
      call daall(iscr,1,'$$DAFUN3$$',nomax,nvmax)
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!
      call dapek(ina,jj,a0)
!
!      no = min(nocut,inoa,inoc)
      no = min(nocut,nomax)
!
!     BRANCHING TO DIFFERENT FUNCTIONS
!     ********************************
!
      if(cf.eq.'INV ') then
!        1/(A0+P) = 1/A0*(1-(P/A0)+(P/A0)**2-...)
        if(a0.eq.0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        xf(0) = 1.d0/a0
        do i=1,no
          xf(i) = -xf(i-1)/a0
        enddo
!
      elseif(cf.eq.'SQRT') then
!        SQRT(A0+P) = SQRT(A0)*(1+1/2(P/A0)-1/8*(P/A0)**2+...)
        if(a0.le.0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        ra = dsqrt(a0)
        xf(0) = ra
        do i=1,no
          xf(i) = -xf(i-1)/a0/dble(2*i)*dble(2*i-3)
        enddo
!
      elseif(cf.eq.'ISRT') then
!        1/SQRT(A0+P) = 1/SQRT(A0)*(1-1/2(P/A0)+3/8*(P/A0)**2-...)
        if(a0.le.0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        era = 1.d0/dsqrt(a0)
        xf(0) = era
        do i=1,no
          xf(i) = -xf(i-1)/a0/dble(2*i)*dble(2*i-1)
        enddo
!
      elseif(cf.eq.'EXP ') then
!        EXP(A0+P) = EXP(A0)*(1+P+P**2/2!+...)
        ea  = dexp(a0)
        xf(0) = ea
        do i=1,no
          xf(i) = xf(i-1)/dble(i)
        enddo
!
      elseif(cf.eq.'LOG ') then
!        LOG(A0+P) = LOG(A0) + (P/A0) - 1/2*(P/A0)**2 + 1/3*(P/A0)**3 - ...)
        if(a0.le.0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        ea  = dlog(a0)
        xf(0) = ea
        xf(1) = 1.d0/a0
        do i=2,no
          xf(i) = -xf(i-1)/a0/dble(i)*dble(i-1)
        enddo
!
      elseif(cf.eq.'SIN ') then
!        SIN(A0+P) = SIN(A0)*(1-P**2/2!+P**4/4!) + COS(A0)*(P-P**3/3!+P**5/5!)
        sa  = dsin(a0)
        ca  = dcos(a0)
        xf(0) = sa
        xf(1) = ca
        do i=2,no
          xf(i) = -xf(i-2)/dble(i*(i-1))
        enddo
!
      elseif(cf.eq.'COS ') then
!        COS(A0+P) = COS(A0)*(1-P**2/2!+P**4/4!) - SIN(A0)*(P-P**3/3!+P**5/5!)
        sa  = dsin(a0)
        ca  = dcos(a0)
        xf(0) = ca
        xf(1) = -sa
        do i=2,no
          xf(i) = -xf(i-2)/dble(i*(i-1))
        enddo
!
      elseif(cf.eq.'SIRX') then
!        SIN(SQRT(P))/SQRT(P) = 1 - P/3! + P**2/5! - P**3/7! + ...
        if(a0.ne.0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        xf(0)=1.d0
        do i=1,no
          xf(i) = -xf(i-1)/dble(2*i*(2*i+1))
        enddo
!
      elseif(cf.eq.'CORX') then
!        COS(SQRT(P)) = 1 - P/2! + P**2/4! - P**3/6! + ...
        if(a0.ne.0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        xf(0)=1.d0
        do i=1,no
          xf(i) = -xf(i-1)/dble(2*i*(2*i-1))
        enddo
!
      elseif(cf.eq.'SIDX') then
!        SIN(P)/P = 1 - P**2/3! + P**4/5! - P**6/7! + ...
        if(a0.ne.0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        xf(0)=1.d0
        xf(1)=0.d0
        do i=2,no
          xf(i) = -xf(i-2)/dble(i*(i+1))
        enddo
!
      elseif(cf.eq.'TAN ') then
        if(dabs(dcos(a0)).lt.epsmac) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        sa  = dsin(a0)
        ca  = dcos(a0)
        xf(0) = sa/ca
        xf(1) = 1.d0/ca/ca
        xf(2) = 2.d0*sa/ca/ca/ca/2.d0
        xf(3) = (2.d0*ca*ca+6.d0*sa*sa)/ca/ca/ca/ca/6.d0
        xf(4) = (16*sa+8.d0*sa*sa*sa)/ca/ca/ca/ca/ca/24.d0
        xf(5) = (16.d0*ca*ca+24.d0*ca*ca*sa*sa+80.d0*sa*sa+             &
     &  40.d0*sa*sa*sa*sa)/ca/ca/ca/ca/ca/ca/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
          stop
        endif
      elseif(cf.eq.'COT ') then
        if(dabs(dsin(a0)).lt.epsmac) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        sa  = dsin(a0)
        ca  = dcos(a0)
        xf(0) = ca/sa
        xf(1) = -1.d0/sa/sa
        xf(2) = 2.d0*ca/sa/sa/sa/2.d0
        xf(3) = -(2.d0*sa*sa+6.d0*ca*ca)/sa/sa/sa/sa/6.d0
        xf(4) = (16*ca+8.d0*ca*ca*ca)/sa/sa/sa/sa/sa/24.d0
        xf(5) = -(16.d0*sa*sa+24.d0*sa*sa*ca*ca+80.d0*ca*ca+            &
     &  40.d0*ca*ca*ca*ca)/sa/sa/sa/sa/sa/sa/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
          stop
        endif
      elseif(cf.eq.'ASIN') then
        if((1.d0-dabs(a0)).lt.0.d0) then
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          lfun = 0
          return
        endif
        xf(0) = dasin(a0)
        xf(1) = (1.d0-a0*a0)**(-0.5d0)
        xf(2) = a0*xf(1)**3.d0/2.d0
        xf(3) = (1+2.d0*a0*a0)*xf(1)**5.d0/6.d0
        xf(4) = (9.d0*a0+6.d0*a0*a0*a0)*xf(1)**7.d0/24.d0
        xf(5) = (9.d0+72.d0*a0*a0+24.d0*a0*a0*a0*a0)*xf(1)**9.d0/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          stop
        endif
      elseif(cf.eq.'ACOS')then
        if((1.d0-dabs(a0)).lt.0.d0) then
          call dadeb(31,'ERR DAFUN ',1)
          write(*,1000) cf,ina,a0
          lfun = 0
          return
        endif
        xf(0) =  dacos(a0)
        scr =  (1.d0-a0*a0)**(-0.5d0)
        xf(1) =  -scr
        xf(2) = -a0*scr**3.d0/2.d0
        xf(3) = -(1+2.d0*a0*a0)*scr**5.d0/6.d0
        xf(4) = -(9.d0*a0+6.d0*a0*a0*a0)*scr**7.d0/24.d0
        xf(5) = -(9.d0+72.d0*a0*a0+24.d0*a0*a0*a0*a0)*scr**9.d0/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'ATAN') then
!        ATAN(A0+P) = ATAN(A0)+1/(1+A0**2)*P-A0/(1+A0**2)**2*P**2+....)
        xf(0) = datan(a0)
        xf(1) = 1.d0/(1.d0+a0*a0)
        xf(2) = -a0*(xf(1)*xf(1))
        xf(3) = (a0*a0-1.d0/3.d0)*xf(1)**3
        xf(4) = (a0-a0*a0*a0)*xf(1)**4
        xf(5) = (1.d0/5.d0+a0**4-2.d0*a0*a0)*xf(1)**5
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'ACOT') then
        xf(0) = 2.d0*datan(1.d0)-datan(a0)
        scr = 1.d0/(1.d0+a0*a0)
        xf(1) = -scr
        xf(2) = a0*(scr*scr)
        xf(3) = -(a0*a0-1.d0/3.d0)*scr**3
        xf(4) = -(a0-a0*a0*a0)*scr**4
        xf(5) = -(1.d0/5.d0+a0**4-2.d0*a0*a0)*scr**5
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'SINH') then
        sa  = dsinh(a0)
        ca  = dcosh(a0)
        xf(0) = sa
        xf(1) = ca
        xf(2) = sa/2.d0
        xf(3) = ca/6.d0
        xf(4) = sa/24.d0
        xf(5) = ca/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'COSH') then
        sa  = dsinh(a0)
        ca  = dcosh(a0)
        xf(0) = ca
        xf(1) = sa
        xf(2) = ca/2.d0
        xf(3) = sa/6.d0
        xf(4) = ca/24.d0
        xf(5) = sa/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'TANH') then
        sa  = dsinh(a0)
        ca  = dcosh(a0)
        xf(0) = sa/ca
        xf(1) = 1.d0/ca/ca
        xf(2) = -2.d0*sa/ca/ca/ca/2.d0
        xf(3) = (-2.d0*ca*ca+6.d0*sa*sa)/ca/ca/ca/ca/6.d0
        xf(4) = (16*sa-8.d0*sa*sa*sa)/ca/ca/ca/ca/ca/24.d0
        xf(5) = (16.d0*ca*ca-24.d0*ca*ca*sa*sa-80.d0*sa*sa+             &
     &  40.d0*sa*sa*sa*sa)/ca/ca/ca/ca/ca/ca/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'COTH') then
        if(dabs(dsinh(a0)).lt.epsmac) then
          lfun = 0
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          return
        endif
        sa  = dsinh(a0)
        ca  = dcosh(a0)
        xf(0) = ca/sa
        xf(1) = -1.d0/sa/sa
        xf(2) =  2.d0*ca/sa/sa/sa/2.d0
        xf(3) = (2.d0*sa*sa-6.d0*ca*ca)/sa/sa/sa/sa/6.d0
        xf(4) = (16*ca+8.d0*ca*ca*ca)/sa/sa/sa/sa/sa/24.d0
        xf(5) = (16.d0*sa*sa+24.d0*sa*sa*ca*ca-80.d0*ca*ca-             &
     &  40.d0*ca*ca*ca*ca)/sa/sa/sa/sa/sa/sa/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'ASNH') then
        xf(0) = dlog(a0+dsqrt(a0*a0+1.d0))
        xf(1) = (1.d0+a0*a0)**(-0.5d0)
        xf(2) = -a0*xf(1)**3.d0/2.d0
        xf(3) = (2.d0*a0*a0-1.d0)*xf(1)**5.d0/6.d0
        xf(4) = (9.d0*a0-6.d0*a0*a0*a0)*xf(1)**7.d0/24.d0
        xf(5) = (9.d0-72.d0*a0*a0+24.d0*a0*a0*a0*a0)*xf(1)**9.d0/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'ACSH') then
        if((1.d0-a0).ge.0.d0) then
          lfun = 0
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          return
        endif
        xf(0) = dlog(a0+dsqrt(a0*a0-1.d0))
        xf(1) = (a0*a0-1.d0)**(-0.5d0)
        xf(2) = -a0*xf(1)**3.d0/2.d0
        xf(3) = (2.d0*a0*a0+1.d0)*xf(1)**5.d0/6.d0
        xf(4) = (-9.d0*a0-6.d0*a0*a0*a0)*xf(1)**7.d0/24.d0
        xf(5) = (9.d0+72.d0*a0*a0+24.d0*a0*a0*a0*a0)*xf(1)**9.d0/120.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'ATNH') then
        if((dabs(a0)-1.d0).ge.0.d0) then
          lfun = 0
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          return
        endif
        xf(0) =  0.5d0*dlog((1+a0)/(1-a0))
        xf(1) =  1.d0/(1.d0-a0*a0)
        xf(2) =  a0*(xf(1)*xf(1))
        xf(3) = (a0*a0+1.d0/3.d0)*xf(1)**3
        xf(4) = (a0+a0*a0*a0)*xf(1)**4
        xf(5) = (1.d0/5.d0+a0**4+2.d0*a0*a0)*xf(1)**5
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      elseif(cf.eq.'ACTH') then
        if(1.d0-dabs(a0).ge.0.d0) then
          lfun = 0
          write(*,1000) cf,ina,a0
          call dadeb(31,'ERR DAFUN ',1)
          return
        endif
        xf(0) =  0.5d0*dlog((a0+1)/(a0-1))
        scr =  1.d0/(-1.d0+a0*a0)
        xf(1) = -scr
        xf(2) =  a0*(scr*scr)
        xf(3) = (-a0*a0-1.d0/3.d0)*scr**3.d0
        xf(4) = (a0+a0*a0*a0)*scr**4.d0
        xf(5) = (-1.d0/5.d0-a0**4-2.d0*a0*a0)*scr**5.d0
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
!      ELSEIF(CF.EQ.'ABF ') THEN
!
!     DIESE FUNKTION BESCHREIBT DEN FELDABFALL BEI IONENOPTISCHEN ELEMENTEN
!     ABF=1/(1+EXP(A0+X))
!        =1/(1+EXP(A0)*(1-EXP(A0)/(1+EXP(A0))*X+....)
!         XF(0) = 1.D0/(1+DEXP(A0))
!         E1  = DEXP(A0)*X1
!         E1  = DEXP(A0)*Xf(0)
!         E2  = E1 * E1
!         E3  = E2 * E1
!         E4  = E3 * E1
!         E5  = E4 * E1
!         XF(1) = X1*(-E1)
!         XF(2) = X1*(-0.5D0* E1 + E2)
!         XF(3) = X1*(-E1/6.D0 + E2 - E3)
!         XF(4) = X1*(-E1/24.D0 + E2*7.D0/12.D0 - E3*3.D0/2.D0 + E4)
!         XF(5) = X1*(-E1/120.D0 + E2/4.D0 - E3*5.D0/4.D0 +
!     *         E4*2.D0 - E5)
!         IF(NO.GT.5) THEN
!            write(6,*)'ERROR IN DAFUN, ',CF, ' ONLY UP TO NO = 5'
!            CALL DADEB(31,'ERR DAFUN ',1)
!         ENDIF
!      ELSEIF(CF.EQ.'GAUS') THEN
!
!     DIESE FUNKTION BESCHREIBT DIE ENTWICKLUNG VON EXP(-X*X)
!
!         XF(0) = DEXP(-A0*A0)
!         XF(1) = -2.D0*A0*X1
!         XF(2) = (-1.D0+2.D0*A0*A0)*X1
!         XF(3) = (12.D0*A0-8.D0*A0*A0*A0)/6.D0*X1
!         XF(4) = (16.D0*A0*A0*A0*A0-48.D0*A0*A0+12.D0)/24.D0*X1
!         XF(5) = (-32.D0*A0*A0*A0*A0*A0+160.D0*A0*A0*A0-120.D0*A0)/
!     *           120.D0*X1
!         IF(NO.GT.5) THEN
!            write(6,*)'ERROR IN DAFUN, ',CF, ' ONLY UP TO NO = 5'
!            CALL DADEB(31,'ERR DAFUN ',1)
!         ENDIF
      elseif(cf.eq.'ERF ') then
!
!    ERF(X) STELLT DAS INTEGRAL VON 0 BIS X VON [ 2/SQRT(PI) * EXP(-X*X) ]
!    DAR
!
        e1 = dexp(-a0*a0)
        a1 = .254829592d0
        a2 = -.284496736d0
        a3 = 1.421413741d0
        a4 = -1.453152027d0
        a5 = 1.061405429d0
        p  = .3275911d0
        rpi4 = sqrt(datan(1.d0))
        t  = 1.d0/(1.d0+p*a0)
        e2 = 1.d0-t*(a1+t*(a2+t*(a3+t*(a4+t*a5))))*e1
        xf(0)= e2
        xf(1) = e1/rpi4
        xf(2) = -a0*e1/rpi4
        xf(3) = (-2.d0+4.d0*a0*a0)/6.d0*e1/rpi4
        xf(4) = (12.d0*a0-8.d0*a0*a0*a0)/24.d0*e1/rpi4
        xf(5) = (16.d0*a0*a0*a0*a0-48.d0*a0*a0+12.d0)/120.d0*e1/rpi4
        if(no.gt.5) then
          write(6,*)'ERROR IN DAFUN, ',cf, ' ONLY UP TO NO = 5'
          call dadeb(31,'ERR DAFUN ',1)
        endif
      else
        write(6,*)'ERROR, UNSOPPORTED FUNCTION ',cf
      endif
!
      call dacon(inc,xf(0))
      call dacop(ina,inon)
      call dapok(inon,jj,0.d0)
      call dacon(ipow,1.d0)
!
      do i=1,min(no,nocut)
!
        call damul(inon,ipow,iscr)
        call dacop(iscr,ipow)
        call dacma(inc,ipow,xf(i),inc)
!
      enddo
!
 1000 format('ERROR IN DAFUN, ',a4,' DOES NOT EXIST FOR VECTOR ',i10,   &
     &'CONST TERM  = ',e12.5)
!
      call dadal(iscr,1)
      call dadal(inon,1)
      call dadal(ipow,1)
!
      return
      end
!

      subroutine daabs(ina,anorm)
      implicit none
      integer i,illa,ilma,ina,inoa,inva,ipoa
      double precision anorm
!     ***************************
!
!     THIS SUBROUTINE COMPUTES THE NORM OF THE DA VECTOR A
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!
      anorm = 0.d0
      do i=ipoa,ipoa+illa-1
        anorm = anorm + dabs(cc(i))
      enddo
!
      return
      end
!
      subroutine daabs2(ina,anorm)
      implicit none
      integer i,illa,ilma,ina,inoa,inva,ipoa
      double precision anorm
!     ***************************
!
!     THIS SUBROUTINE COMPUTES THE NORM OF THE DA VECTOR A
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!
      anorm = 0.d0
      do i=ipoa,ipoa+illa-1
        anorm = anorm + dabs(cc(i)**2)
      enddo
!
      return
      end
!

      subroutine dacom(ina,inb,dnorm)
      implicit none
      integer idacom,illc,ilmc,ina,inb,inc,inoc,invc,ipoc
      double precision dnorm
!     *******************************
!
!     THIS SUBROUTINE COMPARES TWO DA VECTORS BY RETURNING THE NORM
!     OF THE DIFFERENCE
!
      idacom = 0
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
      call daall(idacom,1,'$$DACOM $$',inoc,invc)
      call dasub(ina,inb,idacom)
      call daabs(idacom,dnorm)
      call dadal(idacom,1)
!
      return
      end
!

      subroutine dapos(ina,inb)
      implicit none
      integer ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,    &
     &ipoa,ipob
!     *************************
!
!     THIS SUBROUTINE MAKES THE SIGNS OF ALL THE COEFFICIENTS OF A POSITIVE
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"

      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INB,INOB,INVB)
!
      ib = ipob - 1
!
      do ia = ipoa,ipoa+illa-1
!
        if(ieo(ia1(i1(ia))+ia2(i2(ia))).gt.nocut) goto 100
        ib     = ib + 1
        cc(ib) = dabs(cc(ia))
        i1(ib) = i1(ia)
        i2(ib) = i2(ia)
!
 100    continue
      enddo
!
      idall(inb) = ib - ipob + 1
      if(idall(inb).gt.idalm(inb)) then
        write(6,*)'ERROR IN DAPOS '
        call dadeb(31,'ERR DAPOS ',1)
      endif
!
      return
      end
!
      subroutine dacct(ma,ia,mb,ib,mc,ic)
      implicit none
      integer i,ia,ib,ic,ij,illc,ilmc,inoc,invc,ipoc
!     ***********************************
!
!     THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC
!     WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC
!     DA VECTORS, RESPECTIVELY.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer mon(lnv),mb(*),mc(*),ma(*)

      if(ma(1).eq.mc(1).or.mb(1).eq.mc(1)) then
        call dainf(mc(1),inoc,invc,ipoc,ilmc,illc)
        do ij=1,ic
          mon(ij)=0
        enddo
        call daall(mon,ic,'$$DAJUNK$$',inoc,invc)
        call dacctt(ma,ia,mb,ib,mon,ic)
        do i=1,ic
          call dacop(mon(i),mc(i))
        enddo
        call dadal(mon,ic)
      else
        call dacctt(ma,ia,mb,ib,mc,ic)
      endif

      return
      end

      subroutine dacctt(mb,ib,mc,ic,ma,ia)
      implicit none
      integer i,ia,ib,ic,iia,iib,iic,illa,illb,illc,ilma,ilmb,ilmc,inoa,&
     &inob,inoc,inva,invb,invc,ipoa,ipob,ipoc,iv,jl,jv
      double precision ccf
!     ***********************************
!
!     THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC
!     WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC
!     DA VECTORS, RESPECTIVELY.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!      INTEGER MON(LNO+1),ICC(LNV),MB(*),MC(*),MA(*)
!ETIENNE
      integer mon(lno+1),icc(lno),mb(*),mc(*),ma(*)
!ETIENNE
!
!     CONSISTENCY CHECKS
!     ******************
!
      iia = ma(1)
      iib = mb(1)
      iic = mc(1)
      call dainf(iia,inoa,inva,ipoa,ilma,illa)
      call dainf(iib,inob,invb,ipob,ilmb,illb)
      call dainf(iic,inoc,invc,ipoc,ilmc,illc)
!
      call damch(ma,ia)
      call damch(mb,ib)
!
      if(ia.ne.ib) then
        write(6,*)'ERROR IN DACCT, IA .NE. IB'
        call dadeb(31,'ERR DACCT1',1)
      elseif(ic.ne.invb) then
        write(6,*)'ERROR IN DACCT, IC.NE.INVB'
        call dadeb(31,'ERR DACCT2',1)
      endif
!
!     ALLOCATING LOCAL VECTORS AND CALLING MTREE
!     ******************************************
!
      do i=1,ib
        icc(i) = 0
      enddo
!
      do i=1,nomax+1
        mon(i) = 0
      enddo
!
      call daall(icc,ib,'$$DACCT $$',nomax,nvmax)
      call daall(mon,nomax+1,'$$DAMON $$',inoc,invc)
!
      call mtree(mb,ib,icc,ib)
!
!     PERFORMING CONCATENATION
!     ************************
!
      do i=1,ia
        call dacon(ma(i),cc(idapo(icc(i))))
      enddo
!
      call dacon(mon(1),1.d0)
!
      do i=1,idall(icc(1))-1
!
        jl = i1(idapo(icc(1))+i)
        jv = i2(idapo(icc(1))+i)
!
        call damul(mon(jl),mc(jv),mon(jl+1))
!
        do iv=1,ia
!
          ccf = cc(idapo(icc(iv))+i)
          if(dabs(ccf).gt.eps) call dacma(ma(iv),mon(jl+1),ccf,ma(iv))
!
        enddo
      enddo
!
      call dadal(mon,nomax+1)
      call dadal(icc,ib)
!
      return
      end
!
      subroutine mtree(mb,ib,mc,ic)
      implicit none
      integer i,ib,ib1,ibi,ic,ic1,ic2,icc,ichk,ii,iib,iic,illb,illc,    &
     &ilmb,ilmc,inob,inoc,invb,invc,ipob,ipoc,j,jl,jnon,nterm,ntermf
      double precision apek,bbijj,chkjj
!     *****************************
!
!     THIS SUBROUTINE IS USED FOR CONCATENATION AND TRACKING OF VECTORS
!     THROUGH A DA MAP. IT COMPUTES THE TREE THAT HAS TO BE TRANSVERSED
!     MB IS THE DA MATRIX WITH IA TERMS. THE OUTPUT MC IS A CA MATRIX WHICH
!     CONTAINS COEFFICIENTS AND CONTROL INTEGERS USED FOR THE TRAVERSAL.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(lnv),jv(0:lno),mb(*),mc(*)
!
!     CONSISTENCY CHECKS
!     ******************
!
      iib = mb(1)
      iic = mc(1)
      call dainf(iib,inob,invb,ipob,ilmb,illb)
      call dainf(iic,inoc,invc,ipoc,ilmc,illc)
!
      call damch(mb,ib)
      call damch(mc,ic)
!
      if(ib.ne.ic) then
        write(6,*)'ERROR IN MTREE, IB .NE. IC'
        call dadeb(31,'ERR MTREE1',1)
      endif
!
!     ALLOCATING LOCAL VECTORS
!     ************************
!
      ichk = 0
      call daall(ichk,1,'$$MTREE $$',nomax,nvmax)
!
!     FIND ALL THE ENTRIES TO BE LOOKED FOR
!     *************************************
!
      call daclr(1)
!
      cc(1) = 1.d0
!
      do i=1,ib
        if(nomax.eq.1) then
          do ib1 = 2,7
            cc(ib1) = 1d0
          enddo
        else
          do ibi = idapo(mb(i)),idapo(mb(i))+idall(mb(i))-1
            icc = ia1(i1(ibi)) + ia2(i2(ibi))
            if(ieo(icc).gt.inob) goto 90
            cc(icc) = 1.d0
   90       continue
          enddo
        endif
      enddo
!
      do ii=1,inob
!
!     SEARCHING FOR FATHER FOR EACH TERM
!
        do i=1,nmmax
          if(cc(i).lt.0.5d0) goto 140
!
          jnon = 0
          call dancd(i1(i),i2(i),jj)
          do j=1,invb
            if(jj(j).eq.0) goto 130
            jnon = j
            jj(j) = jj(j) - 1
            call dadcd(jj,ic1,ic2)
            apek = cc(ia1(ic1)+ia2(ic2))
            jj(j) = jj(j) + 1
            if(apek.ge.0.5d0) goto 140
  130       continue
          enddo
!
          if(jnon.eq.0) goto 140
!
!     TERM IS AN ORPHAN, SO CREATE FOSTER FATHER
!
          jj(jnon) = jj(jnon) - 1
          call dadcd(jj,ic1,ic2)
          cc(ia1(ic1)+ia2(ic2)) = 1.d0
!
  140     continue
        enddo
      enddo
!
      call dapac(ichk)
!ETIENNE      CALL DAPRI(ICHK,32)
!
!     SETTING UP TREE STRUCTURE
!     *************************
!
      ntermf = idall(ichk)
!
!     ZEROTH ORDER TERMS
!     ******************
!
      do i=1,lnv
        jj(i) = 0
      enddo
!
      do i=1,ib
        call dapek(mb(i),jj,bbijj)
        i1(idapo(mc(i))) = 0
        i2(idapo(mc(i))) = 0
        cc(idapo(mc(i))) = bbijj
      enddo
!
      call dapek(ichk,jj,chkjj)
      if(chkjj.gt.0.5d0) then
        call dapok(ichk,jj,-1.d0)
      else
        write(6,*)'ERROR IN MTREE, ZEROTH ORDER TERM OF ICHK IS ZERO'
        call dadeb(31,'ERR MTREE2',1)
      endif
!
      nterm = 1
!
!     HIGHER ORDER TERMS
!     ******************
!
      do jl=1,inob
        jv(jl) = 0
      enddo
!
      jl = 0
      chkjj = 1.d0
!
 200  continue
      if(jl.eq.0.and.chkjj.le.0.5d0) goto 250
      if(jl.lt.inob.and.chkjj.gt.0.5d0) then
        jl = jl + 1
        jj(1) = jj(1) + 1
        jv(jl) = 1
      elseif(jv(jl).eq.invb) then
        jj(jv(jl)) = jj(jv(jl)) - 1
        jv(jl) = 0
        jl = jl - 1
        chkjj = 0.d0
        goto 200
      else
        jj(jv(jl)) = jj(jv(jl)) - 1
        jv(jl) = jv(jl) + 1
        jj(jv(jl)) = jj(jv(jl)) + 1
      endif
!
      call dapek(ichk,jj,chkjj)
!
      if(chkjj.le.0.5d0) goto 200
!
      nterm = nterm + 1
      if(nterm.gt.idalm(mc(1))) then
        write(6,*)'ERROR IN MTREE, NTERM TOO LARGE'
        call dadeb(31,'ERR MTREE3',1)
      endif
!
      call dapok(ichk,jj,-1.d0)
!
!     write(6,*)'JL,JV = ',JL,JV(JL)
      do i=1,ib
        call dapek(mb(i),jj,bbijj)
        i1(idapo(mc(i))+nterm-1) = jl
        i2(idapo(mc(i))+nterm-1) = jv(jl)
        cc(idapo(mc(i))+nterm-1) = bbijj
      enddo
!
      goto 200
!
 250  continue
!
      do i=1,ib
        idall(mc(i)) = nterm
      enddo
!
!     PERFORMING CROSS CHECKS
!     ***********************
!
      if(nterm.ne.ntermf.or.nterm.ne.idall(ichk)) then
        write(6,*)'ERROR IN MTREE, NTERM, NTERMF, IDALL(ICHK) =  '      &
     &  ,nterm,ntermf,idall(ichk)
        call dadeb(31,'ERR MTREE4',1)
      endif
!
      do i=idapo(ichk),idapo(ichk)+nterm-1
        if(dabs(cc(i)+1.d0).gt.epsmac) then
          write(6,*)'ERROR IN MTREE, NOT ALL TERMS IN ICHK ARE -1'
          call dadeb(31,'ERR MTREE5',1)
        endif
      enddo
!
      call dadal(ichk,1)
!
      return
      end
!
      subroutine ppushprint(mc,ic,mf,jc,line)
      implicit none
      integer i,ic,iv,jc,jl,jv,mc,mf
      include "TPSALib_prm.f"
      dimension mc(*)
      character*20 line
      if(mf.le.0) return
      write(mf,*) 0,0,jc+1,0,line
      do i=1,ic
        jc=1+jc
        write(mf,*) jc,jl,jv,cc(idapo(mc(i)))
      enddo
!     xf(i) = cc(idapo(mc(i)))
!      xm(1) = 1.d0
      do i=1,idall(mc(1))-1
        jl = i1(idapo(mc(1))+i)
        jv = i2(idapo(mc(1))+i)
!      xx = xm(jl)*xi(jv)
!      xm(jl+1) = xx
        do iv=1,ic
          jc=1+jc
          write(mf,*) jc,jl,jv,cc(idapo(mc(iv))+i)
!        xf(iv) = xf(iv) + cc(idapo(mc(iv))+i) * xx
        enddo
      enddo
      return
      end
!

      subroutine ppush(mc,ic,xi,xf)
      implicit none
      integer i,ic,iv,jl,jv,mc
      double precision xf,xi,xm,xt,xx
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE MATRIX WHOSE TREE IS STORED IN CA VECTOR MC
!     TO THE COORDINATES IN XI AND STORES THE RESULT IN XF
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension mc(*),xf(*),xi(*),xm(lno+1) ,xt(lno)
!
      do i=1,ic
        xt(i)=xi(i)
      enddo
      do i=1,ic
        xf(i) = cc(idapo(mc(i)))
      enddo
!
      xm(1) = 1.d0
!
      do i=1,idall(mc(1))-1
!
        jl = i1(idapo(mc(1))+i)
        jv = i2(idapo(mc(1))+i)
        xx = xm(jl)*xt(jv)
        xm(jl+1) = xx
!
        do iv=1,ic
          xf(iv) = xf(iv) + cc(idapo(mc(iv))+i) * xx
        enddo
      enddo
!
      return
      end

      subroutine ppush1(mc,xi,xf)
      implicit none
      integer i,iv,jl,jv,mc
      double precision xf,xi,xm,xt,xx
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE MATRIX WHOSE TREE IS STORED IN CA VECTOR MC
!     TO THE COORDINATES IN XI AND STORES THE RESULT IN XF
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension xi(*),xm(lno+1) ,xt(lno)
!
      do i=1,nvmax
        xt(i)=xi(i)
      enddo

      xf = cc(idapo(mc))
!
      xm(1) = 1.d0
!
      do i=1,idall(mc)-1
!
        jl = i1(idapo(mc)+i)
        jv = i2(idapo(mc)+i)
        xx = xm(jl)*xt(jv)
        xm(jl+1) = xx
!
        xf = xf + cc(idapo(mc)+i) * xx
      enddo
!
      return
      end

      subroutine dainv(ma,ia,mb,ib)
      implicit none
      integer i,ia,ib,ij,illb,ilmb,inob,invb,ipob
      double precision x
!     *****************************
!
!     THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND
!     STORES THE RESULT IN MI
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(lnv),ml(lnv),ma(*),mb(*)
!
      dimension x(lnv)
!

      do i=1,lnv
        jj(i)=0
      enddo

      if(ma(1).eq.mb(1)) then
        call dainf(mb(1),inob,invb,ipob,ilmb,illb)
        do i=1,ia
          call dapok(ma(i),jj,0.d0)
        enddo
        do ij=1,ib
          ml(ij)=0
        enddo
        call daall(ml,ib,'$$DAJUNK$$',inob,invb)
        call dainvt(ma,ia,ml,ib)
        do i=1,ib
          call dacop(ml(i),mb(i))
        enddo
        call dadal(ml,ib)
      else
        do i=1,ia
          call dapek(ma(i),jj,x(i))
          call dapok(ma(i),jj,0.d0)
        enddo
        call dainvt(ma,ia,mb,ib)
        do i=1,ia
          call dapok(ma(i),jj,x(i))
        enddo
      endif

      return
      end

      subroutine dainvt(ma,ia,mb,ib)
      implicit none
      integer i,ia,ib,ie,ier,illa,illb,ilma,ilmb,inoa,inob,inva,invb,   &
     &ipoa,ipob,j,k,nocut0
      double precision aa,ai,amjj,amsjj,prod
!     *****************************
!
!     THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND
!     STORES THE RESULT IN MI
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(lnv),ms(lnv),ml(lnv),ma(*),mb(*)
!
      dimension aa(lnv,lnv),ai(lnv,lnv)
!
      call dainf(ma(1),inoa,inva,ipoa,ilma,illa)
      call dainf(mb(1),inob,invb,ipob,ilmb,illb)
!
!     CONSISTENCY CHECKS
!     ******************
!
      call damch(ma,ia)
      call damch(mb,ib)
!etienne
      do ie=1,ib
        call dacon(mb(ie),0.d0)
      enddo
!etienne
!
      if(ia.ne.ib) then
        write(6,*)'ERROR IN DAINV, IA .NE. IB'
        call dadeb(31,'ERR DAINV1',1)
      elseif(ia.ne.inva.or.ib.ne.invb) then
        write(6,*)'ERROR IN DAINV, IA.NE.INVA.OR.IB.NE.INVB'
        call dadeb(31,'ERR DAINV2',1)
      endif
!
!     ALLOCATING LOCAL VECTORS
!     ************************
!
      do i=1,ia
        ms(i) = 0
        ml(i) = 0
      enddo
!
      call daall(ms,ia,'$$INV   $$',inoa,inva)
      call daall(ml,ia,'$$INVL  $$',inoa,inva)
!
!     EXTRACTING LINEAR MATRIX, GENERATING NONLINEAR PART OF A
!     ********************************************************
!
      do i=1,ib
        do j=1,ib
          do k=1,ib
            jj(k) = 0
          enddo
          jj(j) = 1
          call dapek(ma(i),jj,amjj)
          if(dabs(amjj).gt.eps) call dapok(ma(i),jj,0.d0)
          aa(i,j) = amjj
        enddo
        call dacmu(ma(i),-1.d0,ma(i))
      enddo
!
!     INVERTING LINEAR MATRIX, CHECKING RESULT AND STORING IN ML
!     **********************************************************
!
      call matinv(aa,ai,ia,lnv,ier)
!
      if(ier.eq.132) then
        write(6,*)'ERROR IN ROUTINE DAINV'
        call dadeb(31,'ERR DAINV3',1)
      endif
!
      ier = 0
      do i=1,ib
        do j=1,ib
          prod = 0.d0
          do k=1,ib
            jj(k) = 0
            prod = prod + aa(i,k)*ai(k,j)
          enddo
          if(i.eq.j) prod = prod - 1.d0
          if(dabs(prod).gt.100*epsmac) then
            write(6,*)'ERROR IN DAINV, INVERSION DID NOT WORK,I,J,PROD',&
     &      ' = ',                                                      &
     &      i,j,prod,epsmac,eps
            ier = 1
!ETIENNE
            return
!ETIENNE
          endif
          jj(j) = 1
          call dapok(mb(i),jj,ai(i,j))
          call dapok(ml(i),jj,ai(i,j))
        enddo
      enddo
!
      if(ier.eq.1) call dadeb(31,'ERR DAINV4',1)
!
!     ITERATIVELY COMPUTING DIFFERENT PARTS OF THE INVERSE
!     ****************************************************
!
!     MB (OF ORDER I) = A1^-1 o [ E - ANL (NONLINEAR) o MB (OF ORDER I) ]
!
      nocut0 = nocut
!
      do i=2,nocut
!
        nocut = i
!
        call dacct(ma,ia,mb,ib,ms,ia)
        do j=1,ib
          do k=1,ib
            jj(k) = 0
          enddo
          jj(j) = 1
          call dapek(ms(j),jj,amsjj)
          call dapok(ms(j),jj,amsjj+1.d0)
        enddo
!
        call dacct(ml,ia,ms,ia,mb,ib)
!
      enddo
!
      nocut = nocut0
!
!     FLIPPING BACK SIGN OF A, FILLING UP FIRST ORDER PART AGAIN
!     **********************************************************
!
      do i=1,ib
        call dacmu(ma(i),-1.d0,ma(i))
        do j=1,ib
          do k=1,ib
            jj(k) = 0
          enddo
          jj(j) = 1
          call dapok(ma(i),jj,aa(i,j))
        enddo
      enddo
!
      call dadal(ml,ia)
      call dadal(ms,ia)
!
      return
      end
!
      subroutine matinv(a,ai,n,nmx,ier)
      implicit none
      integer i,ier,indx,j,n,nmax,nmx
      double precision a,ai,aw,d
!     *********************************
!
!     THIS SUBROUTINE INVERTS THE MATRIX A AND STORES THE RESULT IN AI
!     INPUT  A   - SAVED
!            N   - ORDER OF MATRIX < 100
!     OUTPUT AI  - A INVERSE
!            IER - 0 NO ERROR
!                  132 ZERO DETERMINANT
!
      parameter (nmax=400)
      dimension a(nmx,nmx),ai(nmx,nmx),aw(nmax,nmax),indx(nmax)

      do i=1,n
        do j=1,n
          aw(i,j) = a(i,j)
          ai(i,j) = 0.0
        enddo
        ai(i,i) = 1.d0
      enddo

      call ludcmp(aw,n,nmax,indx,d,ier)
      if (ier .eq. 132) return
      do j=1,n
        call lubksb(aw,n,nmax,indx,ai(1,j),nmx)
      enddo
!
      return
      end
!
      subroutine ludcmp(a,n,np,indx,d,ier)
      implicit none
      integer i,ier,imax,indx,j,k,n,nmax,np
      double precision a,aamax,d,dum,sum,tiny,vv
!     ************************************
!
!     THIS SUBROUTINE DECOMPOSES A MATRIX INTO LU FORMAT
!     INPUT A: NXN MATRIX - WILL BE OVERWRITTEN BY THE LU DECOMP.
!           NP: PHYSICAL DIMENSION OF A
!           INDX: ROW PERMUTATION VECTOR
!           D: EVEN OR ODD ROW INTERCHANGES
!
!     REFERENCE: NUMERICAL RECIPIES BY PRESS ET AL (CAMBRIDGE) PG. 35
!
      parameter (nmax = 400, tiny = 1.0e-20)
      dimension a(np,np), indx(np), vv(nmax)
      ier=0.
      d=1.d0
      do i=1,n
        aamax=0.d0
        do j=1,n
          if(dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j))
        enddo
        if(aamax.eq.0.d0) then
          ier=132
          return
        endif
        vv(i)=1.d0/aamax
      enddo
      do j=1,n
        if(j.gt.1) then
          do i=1,j-1
            sum=a(i,j)
            if(i.gt.1) then
              do k=1,i-1
                sum=sum-a(i,k)*a(k,j)
              enddo
              a(i,j)=sum
            endif
          enddo
        endif
        aamax=0.d0
        do i=j,n
          sum=a(i,j)
          if (j.gt.1) then
            do k=1,j-1
              sum=sum-a(i,k)*a(k,j)
            enddo
            a(i,j)=sum
          endif
          dum=vv(i)*dabs(sum)
          if(dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
        enddo
        if (j.ne.imax) then
          do k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
          enddo
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(j.ne.n) then
          if(a(j,j).eq.0.d0) a(j,j)=tiny
          dum=1./a(j,j)
          do i=j+1,n
            a(i,j)=a(i,j)*dum
          enddo
        endif
      enddo
      if(a(n,n).eq.0.d0) a(n,n)=tiny
      return
      end
!
      subroutine lubksb(a,n,np,indx,b,nmx)
      implicit none
      integer i,ii,indx,j,ll,n,nmx,np
      double precision a,b,sum
!     ************************************
!
!     THIS SUBROUTINE SOLVES SET OF LINEAR EQUATIONS AX=B,
!     INPUT A: NXN MATRIX IN lu FORM GIVEN BY LUDCMP
!           NP: PHYSICAL DIMENSION OF A
!           INDX: ROW PERMUTATION VECTOR
!           D: EVEN OR ODD ROW INTERCHANGES
!           B: RHS OF LINEAR EQUATION - WILL BE OVERWRITTEN BY X
!
!     REFERENCE: NUMERICAL RECIPIES BY PRESS ET AL (CAMBRIDGE) PG. 36
!
      dimension a(np,np), indx(np), b(nmx)
      ii = 0
      do i=1,n
        ll = indx(i)
        sum = b(ll)
        b(ll) = b(i)
        if(ii.ne.0) then
          do j=ii,i-1
            sum = sum-a(i,j)*b(j)
          enddo
        else if (sum.ne.0.d0) then
          ii = i
        endif
        b(i)=sum
      enddo
      do i=n,1,-1
        sum=b(i)
        if(i.lt.n) then
          do j=i+1,n
            sum = sum-a(i,j)*b(j)
          enddo
        endif

        b(i)=sum/a(i,i)

      enddo
      return
      end
!

      subroutine dapin(ma,ia,mb,ib,jx)
      implicit none
      integer i,ia,ib,ij,illb,ilmb,inob,invb,ipob
      double precision x
!     *****************************
!
!     THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND
!     STORES THE RESULT IN MI
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(lnv),ml(lnv),ma(*),mb(*),jx(*)
!
      dimension x(lnv)
!

      do i=1,lnv
        jj(i)=0
      enddo

      if(ma(1).eq.mb(1)) then
        call dainf(mb(1),inob,invb,ipob,ilmb,illb)
        do i=1,ia
          call dapok(ma(i),jj,0.d0)
        enddo
        do ij=1,ib
          ml(ij)=0
        enddo
        call daall(ml,ib,'$$DAJUNK$$',inob,invb)
        call dapint(ma,ia,ml,ib,jx)
        do i=1,ib
          call dacop(ml(i),mb(i))
        enddo
        call dadal(ml,ib)
      else
        do i=1,ia
          call dapek(ma(i),jj,x(i))
          call dapok(ma(i),jj,0.d0)
        enddo
        call dapint(ma,ia,mb,ib,jx)
        do i=1,ia
          call dapok(ma(i),jj,x(i))
        enddo
      endif

      return
      end

      subroutine dapint(ma,ia,mb,ib,jind)
      implicit none
      integer i,ia,ib,illa,ilma,inoa,inva,ipoa,k
!     **********************************
!
!     THIS SUBROUTINE PERFORMS A PARTIAL INVERSION OF THE ROWS MARKED WITH
!     NONZERO ENTRIES IN JJ OF THE MATRIX A. THE RESULT IS STORED IN B.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(lnv),jind(*),ma(*),mb(*),mn(lnv),mi(lnv),me(lnv)
!
      call dainf(ma(1),inoa,inva,ipoa,ilma,illa)
!

      do i=1,ia
        mn(i) = 0
        mi(i) = 0
        me(i) = 0
      enddo
!
      call daall(mn,ia,'$$PIN1  $$',inoa,inva)
      call daall(mi,ia,'$$PIN2  $$',inoa,inva)
      call daall(me,ia,'$$PIN3  $$',inoa,inva)
!
      do i=1,ia
        do k=1,nvmax
          jj(k) = 0
        enddo
        jj(i) = 1
        call dapok(me(i),jj,1.d0)
      enddo
!
      do i=1,ia
        call dacop(ma(i),mn(i))
        if(jind(i).eq.0) call dacop(me(i),mn(i))
      enddo
!
      call dainv(mn,ia,mi,ia)
!
      do i=1,ia
        if(jind(i).eq.0) call dacop(ma(i),me(i))
      enddo
!
      call dacct(me,ia,mi,ia,mb,ib)
!
      call dadal(me,ia)
      call dadal(mi,ia)
      call dadal(mn,ia)
!
      return
      end
!
      subroutine dader(idif,ina,inc)
      implicit none
      integer idif,illc,ilmc,ina,inc,incc,inoc,invc,ipoc
!     ******************************
!
!     THIS SUBROUTINE COMPUTES THE DERIVATIVE WITH RESPECT TO VARIABLE I
!     OF THE VECTOR A AND STORES THE RESULT IN C.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!

      if(ina.eq.inc) then
        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',inoc,invc)
        call dadert(idif,ina,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dadert(idif,ina,inc)
      endif

      return
      end

      subroutine dadert(idif,ina,inc)
      implicit none
      integer i,ibase,ic,ider1,ider1s,ider2,ider2s,idif,iee,ifac,illa,  &
     &illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj
      double precision rr,x,xdivi
!     ******************************
!
!     THIS SUBROUTINE COMPUTES THE DERIVATIVE WITH RESPECT TO VARIABLE I
!     OF THE VECTOR A AND STORES THE RESULT IN C.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer jd(lnv)
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
!
      if(nomax.eq.1) then
!         PRINT*,'ERROR, DADER CALLED WITH NOMAX = 1'
!        CALL DADEB(31,'ERR DADER1',1)
!         stop 666
        do i=1,lnv
          jd(i)=0
        enddo
        jd(idif)=1
        call dapek(ina,jd,rr)
        call dacon(inc,rr)
        return
      endif
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!
      ibase = nomax + 1
!
      if(idif.gt.(nvmax+1)/2) then
        ider1  = 0
        ider1s = 0
        ider2  = idif-(nvmax+1)/2
        ider2s = 1
        do jj=1,ider2-1
          ider2s = ider2s*ibase
        enddo
        xdivi  = ider2s*ibase
      else
        ider1  = idif
        ider1s = 1
        do jj=1,ider1-1
          ider1s = ider1s*ibase
        enddo
        ider2  = 0
        ider2s = 0
        xdivi  = ider1s*ibase
      endif
!
      ibase = nomax+1
!
      ic = ipoc-1
!
      do i=ipoa,ipoa+illa-1
!
        if(ider1.eq.0) then
          iee = i2(i)
        else
          iee = i1(i)
        endif
!
        x = iee/xdivi
        ifac = int(ibase*(x-int(x+epsmac)+epsmac))
!
        if(ifac.eq.0) goto 100
!
        ic = ic + 1
        cc(ic) = cc(i)*ifac
        i1(ic) = i1(i) - ider1s
        i2(ic) = i2(i) - ider2s
!
 100    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DADER '
        call dadeb(31,'ERR DADER2',1)
      endif
!
      return
      end
!
      subroutine dapoi(ina,inb,inc,n)
      implicit none
      integer i,ina,inb,inc,n
!     *******************************
!
!     THIS SUBROUTINE COMPUTES THE POISSON BRACKET OF THE VECTORS A AND
!     B AND STORES THE RESULT IN C. N IS THE DEGREE OF FREEDOM OF THE SYSTEM.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer is(4)
!
      is(1) = 0
      is(2) = 0
      is(3) = 0
      is(4) = 0
      call daall(is,4,'$$DAPOI $$',nomax,nvmax)
!
!
      do i=1,n
!
        call dader(2*i-1,ina,is(1))
        call dader(2*i,  inb,is(2))
        call damul(is(1),is(2),is(3))
        call daadd(is(4),is(3),is(1))
        call dacop(is(1),is(4))
!
        call dader(2*i,  ina,is(1))
        call dader(2*i-1,inb,is(2))
        call damul(is(1),is(2),is(3))
        call dasub(is(4),is(3),is(1))
        call dacop(is(1),is(4))
!
      enddo

      call dacop(is(4),inc)
!
      call dadal(is,4)
!
      return
      end
!
      subroutine dacfuR(ina,fun,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
      double complex fun
      external fun
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
!     OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
!     RESULT IN C
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!

      if(ina.eq.inc) then
        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',inoc,invc)
        call dacfuRt(ina,fun,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dacfuRt(ina,fun,inc)
      endif

      return
      end

      subroutine dacfuRt(ina,fun,inc)
      implicit none
      integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,  &
     &ipoa,ipoc,j
      double precision cfac,rr
      double COMPLEX fun
      external fun
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
!     OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
!     RESULT IN C
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension j(lnv)
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
!
      if(nomax.eq.1) then
        do i=1,lnv
          j(i)=0
        enddo
        call dapek(ina,j,rr)
        cfac = DREAL(fun(j))
        rr=cfac*rr
        call dapok(inc,j,rr)
        do i=1,lnv
          j(i)=1
          call dapek(ina,j,rr)
          cfac = DREAL(fun(j))
          rr=cfac*rr
          call dapok(inc,j,rr)
          j(i)=0
        enddo
!         PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1'
!        CALL DADEB(31,'ERR DACFU ',1)
!         stop 667
        return
      endif
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!
      ic = ipoc - 1
!
      do ia=ipoa,ipoa+illa-1
!
        call dancd(i1(ia),i2(ia),j)
        cfac = DREAL(fun(j))
!      IF(dABS(CFAC).LT.EPS) GOTO 100
!      IF(dABS(CFAC*CC(IA)).LT.EPS) GOTO 100
        if(dabs(cfac*cc(ia)).lt.eps.or.dabs(cc(ia)).lt.eps) goto 100
!
        ic = ic + 1
        cc(ic) = cc(ia)*cfac
        i1(ic) = i1(ia)
        i2(ic) = i2(ia)
!
 100    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DACFU '
        call dadeb(31,'ERR DACFU ',1)
      endif
!
      return
      end
!
      subroutine dacfu(ina,fun,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
      double precision fun
      external fun
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE PRECISION FUNCTION
!     OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
!     RESULT IN C
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!

      if(ina.eq.inc) then
        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',inoc,invc)
        call dacfut(ina,fun,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dacfut(ina,fun,inc)
      endif

      return
      end




      subroutine dacfuI(ina,fun,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
      double complex fun
      external fun
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
!     OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
!     RESULT IN C
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!

      if(ina.eq.inc) then
        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',inoc,invc)
        call dacfuIt(ina,fun,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call dacfuIt(ina,fun,inc)
      endif

      return
      end

      subroutine dacfuIt(ina,fun,inc)
      implicit none
      integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,  &
     &ipoa,ipoc,j
      double precision cfac,rr
      double COMPLEX fun
      external fun
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION
!     OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
!     RESULT IN C
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension j(lnv)
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
!
      if(nomax.eq.1) then
        do i=1,lnv
          j(i)=0
        enddo
        call dapek(ina,j,rr)
        cfac = DIMAG(fun(j))
        rr=cfac*rr
        call dapok(inc,j,rr)
        do i=1,lnv
          j(i)=1
          call dapek(ina,j,rr)
          cfac = DIMAG(fun(j))
          rr=cfac*rr
          call dapok(inc,j,rr)
          j(i)=0
        enddo
!         PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1'
!        CALL DADEB(31,'ERR DACFU ',1)
!         stop 667
        return
      endif
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!
      ic = ipoc - 1
!
      do ia=ipoa,ipoa+illa-1
!
        call dancd(i1(ia),i2(ia),j)
        cfac = DIMAG(fun(j))
!      IF(dABS(CFAC).LT.EPS) GOTO 100
!      IF(dABS(CFAC*CC(IA)).LT.EPS) GOTO 100
        if(dabs(cfac*cc(ia)).lt.eps.or.dabs(cc(ia)).lt.eps) goto 100
!
        ic = ic + 1
        cc(ic) = cc(ia)*cfac
        i1(ic) = i1(ia)
        i2(ic) = i2(ia)
!
 100    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DACFU '
        call dadeb(31,'ERR DACFU ',1)
      endif
!
      return
      end
!

      subroutine dacfut(ina,fun,inc)
      implicit none
      integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,  &
     &ipoa,ipoc,j
      double precision cfac,fun,rr
      external fun
!     *****************************
!
!     THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE PRECISION FUNCTION
!     OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
!     RESULT IN C
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension j(lnv)
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
!
      if(nomax.eq.1) then
        do i=1,lnv
          j(i)=0
        enddo
        call dapek(ina,j,rr)
        cfac = fun(j)
        rr=cfac*rr
        call dapok(inc,j,rr)
        do i=1,lnv
          j(i)=1
          call dapek(ina,j,rr)
          cfac = fun(j)
          rr=cfac*rr
          call dapok(inc,j,rr)
          j(i)=0
        enddo
!         PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1'
!        CALL DADEB(31,'ERR DACFU ',1)
!         stop 667
        return
      endif
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!
      ic = ipoc - 1
!
      do ia=ipoa,ipoa+illa-1
!
        call dancd(i1(ia),i2(ia),j)
        cfac = fun(j)
!      IF(dABS(CFAC).LT.EPS) GOTO 100
!      IF(dABS(CFAC*CC(IA)).LT.EPS) GOTO 100
        if(dabs(cfac*cc(ia)).lt.eps.or.dabs(cc(ia)).lt.eps) goto 100
!
        ic = ic + 1
        cc(ic) = cc(ia)*cfac
        i1(ic) = i1(ia)
        i2(ic) = i2(ia)
!
 100    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DACFU '
        call dadeb(31,'ERR DACFU ',1)
      endif
!
      return
      end
!

      SUBROUTINE DAIMP(r, ic1, ic2, ina)
      implicit none
      integer           ic1, ic2, ic, ina, i, lh
      double precision  r
*     **************************
*
*     THIS SUBROUTINE "IMPORTS" THE ARRAY H WITH LENGTH LH AND PUTS ITS
*     ENTRIES INTO THE DA VECTOR A
*
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9

      DIMENSION r(nmmax+1), ic1(nmmax), ic2(nmmax)

      lh = nint(r(1))
      call daclr(1)
      ic = idapo(ina)

      do i = 1, lh
        cc(ic) = r(i+1)
        i1(ic) = ic1(i)
        i2(ic) = ic2(i)
        ic = ic + 1
      enddo

      idall(ina) = lh

      RETURN
      END
*
      SUBROUTINE DAEXP(ina, r, ic1, ic2, name)
      implicit none
      character         name*10
      integer           ina, ic1, ic2, ic, i, lh, k, jj
      double precision  r
*     **************************
*
*     THIS SUBROUTINE "EXPORTS" THE DA VACTOR A AND STORES ITS COEFFICIENTS
*     IN THE ARRAY H WITH LENGTH LH
*
      integer j
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
      dimension j(lnv)
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
*
      DIMENSION r(nmmax+1), ic1(nmmax), ic2(nmmax), jj(lnv)

!      call dapri77(ina, 6)
      ic = idapo(ina)
      lh = idall(ina)
      r(1) = lh
      name = daname(ina)

      do k = 1, lnv
         jj(k) = 0
      enddo
      do i = 1, lh
         r(i+1) = cc(ic)
         if (nomax .eq. 1) then
            if (i > 1) jj(i-1) = 1;
            call hash(nomax, lnv, jj, ic1(i), ic2(i))
            if (i > 1) jj(i-1) = 0;
         else
            ic1(i) = i1(ic)
            ic2(i) = i2(ic)
         endif
         ic = ic + 1
      enddo

      RETURN
      END
*


      subroutine dapri(ina,iunit)
      implicit none
      integer i,ii,iii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit,j,k
!     ***************************
!       Frank
!     THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
      dimension j(lnv)
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
!
      if(ina.lt.1.or.ina.gt.nda) then
        print*,'ERROR IN DAPRI, INA = ',ina
!        X = SQRT(-ONE)
!        PRINT*,X
        stop
      endif
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
!
      write(iunit,'(/1X,A,A,I5,A,I5,A,I5/1X,A/)')                       &
     &daname(ina),', NO =',inoa,', NV =',inva,', INA =',ina,            &
     &'***********'//'**********************************'
!
      iout = 0
      ioa = 0

      if(inva.eq.0) then
        write(iunit,'(A)')                                              &
     &  '    I  VALUE  '
        do i = ipoa,ipoa+illa-1
          write(iunit,'(I6,2X,G20.14)') i-ipoa, cc(i)
        enddo
      elseif(nomax.eq.1) then
        if(illa.ne.0) write(iunit,'(A)')                                &
     &  '    I  COEFFICIENT          ORDER   EXPONENTS'
        if(illa.eq.0) write(iunit,'(A)') '   ALL COMPONENTS ZERO '
        do i=1,illa
          do k=1,inva
            j(k)=0
          enddo
          iout=iout+1
          if(i.ne.1) then
            j(i-1)=1
            ioa=1
          endif
          if(abs(cc(ipoa+i-1)).gt.eps) then
            write(iunit,'(I6,2X,E21.14,I5,4X,18(2I2,1X))')              &
     &            iout,cc(ipoa+i-1),ioa,(j(iii),iii=1,nvmax)
!            write(iunit,*) cc(ipoa+i-1)
          endif
        enddo
      else
        if(illa.ne.0) write(iunit,'(A)')                                &
     &  '    I  COEFFICIENT          ORDER   EXPONENTS'
        if(illa.eq.0) write(iunit,'(A)') '   ALL COMPONENTS ZERO '
        do ioa = 0,inoa
          do ii=ipoa,ipoa+illa-1
            if(ieo(ia1(i1(ii))+ia2(i2(ii))).ne.ioa) goto 100
            call dancd(i1(ii),i2(ii),j)
!ETIENNE
            if(abs(cc(ii)).gt.eps) then
!ETIENNE
              iout = iout+1
              write(iunit,'(I6,2X,E21.14,I5,4X,18(2I2,1X))')            &
     &        iout,cc(ii),ioa,(j(iii),iii=1,nvmax)
!ETIENNE
!              write(iunit,*) cc(ii)
            endif
!ETIENNE
!
 100        continue
          enddo
        enddo
!
      endif

      write(iunit,'(A)') '                                      '
!
      return
      end

      subroutine dapri77(ina,iunit)
      implicit none
      integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit,j
      character c10*10,k10*10
!     ***************************
!       Etienne
!     THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      dimension j(lnv)
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      if(iunit.eq.0) return
!
      if(ina.lt.1.or.ina.gt.nda) then
        write(6,*)'ERROR IN DAPRI, INA = ',ina
        stop
      endif
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
!
!      WRITE(IUNIT,*) INA, ' in dapri ', DANAME(INA)
!      WRITE(6,*) INA, ' in dapri ', DANAME(INA)
! 611  WRITE(6,*) ' MORE '
!        READ(5,*) MORE
!        IF(MORE.GT.0) THEN
!        WRITE(6,*) MORE,' ',DANAME(MORE)
!        GOTO 611
!        ENDIF
      write(iunit,'(/1X,A10,A6,I5,A6,I5,A7,I5/1X,A/)')                  &
     &daname(ina),', NO =',inoa,', NV =',inva,', INA =',ina,            &
     &'***********'//'**********************************'
!
      if(illa.ne.0) write(iunit,'(A)')                                  &
     &'    I  COEFFICIENT          ORDER   EXPONENTS'
      if(illa.eq.0) write(iunit,'(A)') '   ALL COMPONENTS ZERO '
!
      c10='      NO ='
      k10='      NV ='

      write(iunit,'(A10,I6,A10,I6)') c10,inoa,k10,inva

      iout = 0
!
!      DO 100 IOA = 0,INOA
      do ioa = 0,nocut
        do ii=ipoa,ipoa+illa-1
          if(nomax.ne.1) then
            if(ieo(ia1(i1(ii))+ia2(i2(ii))).ne.ioa) goto 100
          endif
!ETIENNE
          if(dabs(cc(ii)).gt.eps) then
!ETIENNE

            if(nomax.ne.1) then
              call dancd(i1(ii),i2(ii),j)
              iout = iout+1
            else
              if(ii.eq.ipoa.and.ioa.eq.1) goto 100
              if(ii.gt.ipoa.and.ioa.eq.0) goto 100
              do i=1,lnv
                j(i)=0
              enddo
              if(ii.ne.ipoa) j(ii-ipoa)=1
              iout = iout+1
            endif
!

!      WRITE(IUNIT,*) IOA,CC(II),(J(I),I=1,INVA)
            if(dabs(cc(ii)).gt.eps) then
              if(eps.gt.1.e-37) then
                write(iunit,501) ioa,cc(ii),(j(i),i=1,inva)
              else
                write(iunit,503) ioa,cc(ii),(j(i),i=1,inva)
              endif
            endif
 501        format(' ', i3,1x,g23.16,1x,100(1x,i2))
 503        format(' ', i3,1x,g23.16,1x,100(1x,i2))
 502        format(' ', i5,1x,g23.16,1x,100(1x,i2))

          endif
!ETIENNE
!
 100      continue
        enddo
      enddo
!
      do i=1,lnv
        j(i)=0
      enddo

      if(iout.eq.0) iout=1

      write(iunit,502) -iout,0.d0,(j(i),i=1,inva)
!
      return
      end

      subroutine dashift(ina,inc,ishift)
      implicit none
      integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit,j
      double precision c
!       Frank
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
      dimension j(lnv)
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      integer inb,ishift,ich,ik,jd(lnv),inc
!-----------------------------------------------------------------------------3
!
!

      inb=0
      if(ina.lt.1.or.ina.gt.nda) then
        write(6,*)'ERROR IN DAPRI, INA = ',ina
        stop
      endif
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
      call daall(inb,1,'$$DAJUNK$$',inoa,inva)

!      WRITE(IUNIT,*) INA, ' in dapri ', DANAME(INA)
!      WRITE(6,*) INA, ' in dapri ', DANAME(INA)
! 611  WRITE(6,*) ' MORE '
!        READ(5,*) MORE
!        IF(MORE.GT.0) THEN
!        WRITE(6,*) MORE,' ',DANAME(MORE)
!        GOTO 611
!        ENDIF
      iout = 0
!
!      DO 100 IOA = 0,INOA
      do ioa = 0,nocut
        do ii=ipoa,ipoa+illa-1
          if(nomax.ne.1) then
            if(ieo(ia1(i1(ii))+ia2(i2(ii))).ne.ioa) goto 100
          endif
!ETIENNE
          if(dabs(cc(ii)).gt.eps) then
!ETIENNE

            if(nomax.ne.1) then
              call dancd(i1(ii),i2(ii),j)
              iout = iout+1
            else
              if(ii.eq.ipoa.and.ioa.eq.1) goto 100
              if(ii.gt.ipoa.and.ioa.eq.0) goto 100
              do i=1,lnv
                j(i)=0
              enddo
              if(ii.ne.ipoa) j(ii-ipoa)=1
              iout = iout+1
            endif
!

!      WRITE(IUNIT,*) IOA,CC(II),(J(I),I=1,INVA)
            if(dabs(cc(ii)).gt.eps) then
              if(eps.gt.1.e-37) then
!       write(iunit,501) ioa,cc(ii),(j(i),i=1,inva)
                ich=1
                do ik=1,ishift
                  if(j(ik).ne.0) ich=0
                enddo
                if(ich.eq.1) then
                  do ik=1,lnv
                    jd(ik)=0
                  enddo
                  do ik=ishift+1,lnv
                    jd(ik-ishift)=j(ik)  !%%%%%%Etienne
                  enddo
                endif
                call dapok(inb,jd,cc(ii))
              else
!       write(iunit,503) ioa,cc(ii),(j(i),i=1,inva)
                ich=1
                do ik=1,ishift
                  if(j(ik).ne.0) ich=0
                enddo
                if(ich.eq.1) then
                  do ik=1,lnv
                    jd(ik)=0
                  enddo
                  do ik=ishift+1,lnv
                    jd(ik-ishift)=j(ik)  !%%%%%%Etienne
                  enddo
                endif
                call dapok(inb,jd,cc(ii))
              endif
            endif
 501        format(' ', i3,1x,g23.16,1x,100(1x,i2))
 503        format(' ', i3,1x,g23.16,1x,100(1x,i2))
 502        format(' ', i5,1x,g23.16,1x,100(1x,i2))

          endif
!ETIENNE
!
 100      continue
        enddo
      enddo
!
      do i=1,lnv
        j(i)=0
      enddo

      call dacop(inb,inc)
      call dadal(inb,1)
!
      return
      end

      subroutine darea(ina,iunit)
      implicit none
      integer i,ic,iche,ii,ii1,ii2,iin,illa,ilma,ina,inoa,inva,io,io1,  &
     &ipoa,iunit,iwarin,iwarno,iwarnv,j,nno
      double precision c
!       Frank
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      character c10*10
      dimension j(lnv)
!
      if(ina.lt.1.or.ina.gt.nda) then
        print*,'ERROR IN DAREA, INA = ',ina
!        X = SQRT(-ONE)
!        PRINT*,X
        stop
      endif
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
!
      do i=1,lnv
        j(i) = 0
      enddo
!
      call daclr(1)
!
      ic = 0
!
      iwarno = 0
      iwarnv = 0
      iwarin = 0
!
      read(iunit,'(A10)') c10
      read(iunit,'(18X,I4)') nno
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10

!
!
      iin = 0
!
  10  continue
      iin = iin + 1
      read(iunit,'(I6,2X,G20.14,I5,4X,18(2I2,1X))')                     &
     &ii,c,io,(j(i),i=1,inva)
!
      if(ii.eq.0) goto 20
!ETIENNE
      read(iunit,*) c
!ETIENNE
      if(ii.ne.iin) then
        iwarin = 1
      endif
      io1 = 0
      do i=1,inva
        io1 = io1 + j(i)
      enddo
!
      if(io1.ne.io) then
        iwarnv = 1
        goto 10
      endif
      if(io.gt.inoa) then
!        IF(IWARNO.EQ.0) PRINT*,'WARNING IN DAREA, FILE ',
!    *              'CONTAINS HIGHER ORDERS THAN VECTOR '
        iwarno = 1
        goto 10
      endif
!
      if(nomax.ne.1) then
        ic = ic + 1
        call dadcd(j,ii1,ii2)
        ic = ia1(ii1) + ia2(ii2)
        cc(ic) = c
        goto 10
      else
        iche=0
        do i=1,inva
          if(j(i).eq.1) iche=i
        enddo
        cc(ipoa+iche)=c
        goto 10
      endif

!
  20  continue
!
      if(nomax.ne.1) call dapac(ina)
!
      return
      end
!FF
!
      subroutine darea77(ina,iunit)
      implicit none
      integer i,ic,iche,ii,ii1,ii2,iin,illa,ilma,ina,inoa,inva,ipoa,    &
     &iunit,j,k,nojoh,nvjoh
      double precision c
!     ***************************
!     Etienne
!     THIS SUBROUTINE READS THE DA VECTOR INA FROM UNIT IUNIT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      character c10*10,k10*10
      dimension j(lnv)
!
      if(ina.lt.1.or.ina.gt.nda) then
        write(6,*)'ERROR IN DAPRI, INA = ',ina
        stop
      endif
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
!
      do i=1,lnv
        j(i) = 0
      enddo
!
      call daclr(1)
!
!
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10,I6,A10,I6)') c10,nojoh,k10,nvjoh
!
      iin = 0
!
  10  continue
      iin = iin + 1
      read(iunit,*) ii,c,(j(k),k=1,nvjoh)
      if(ii.lt.0) goto 20

      do i=inva+1,nvjoh
        if(j(i).ne.0) goto 10
      enddo
      iche=0
      do i=1,inva
        iche=iche+j(i)
      enddo
      if(iche.gt.nomax) goto 10
      if(nomax.ne.1) then
        call dadcd(j,ii1,ii2)
        ic = ia1(ii1) + ia2(ii2)
        cc(ic) = c
      else
        iche=0
        do i=1,inva
          if(j(i).eq.1) iche=i
        enddo
        cc(ipoa+iche)=c
      endif
      goto 10
!
  20  continue
!
      if(nomax.ne.1) call dapac(ina)
!
      return
      end

      subroutine dadeb(iunit,c,istop)
      implicit none
      integer istop,iunit
!     *******************************
!
!     THIS SUBROUTINE SERVES AS A DEBUGGING TOOL. IT PRINTS ALL
!     NONZERO INFORMATION IN THE COMMON BLOCKS AND ALL DA  VECTORS.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      character c*10
!
!etienne

      write(6,*) '  ',c
      call flush(6)
!     stop
!      write(6,*) daname(lda+1)
      end
!
!
!
      subroutine danum(no,nv,numda)
      implicit none
      integer i,mm,no,numda,nv
!     *****************************
!
!     THIS SUBROUTINE COMPUTES THE NUMBER OF MONOMIALS OF
!     ORDER NO AND NUMBER OF VARIABLES NV
!
      numda = 1
      mm = max(nv,no)
!
      do i=1,min(nv,no)
        numda = (numda*(mm+i))/i
      enddo
!
      return
      end
!
      subroutine dainf(inc,inoc,invc,ipoc,ilmc,illc)
      implicit none
      integer illc,ilmc,inc,inoc,invc,ipoc
!     **********************************************
!
!     THIS SUBROUTINE SEARCHES THE NUMBER OF DA VECTOR C
!     AND RETURS THE INFORMATION IN COMMON DA
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      if(inc.ge.1.and.inc.le.nda) then
        inoc = idano(inc)
        invc = idanv(inc)
        ipoc = idapo(inc)
        ilmc = idalm(inc)
        illc = idall(inc)
        return
      endif
!
      write(6,*) 'ERROR IN DAINF, DA VECTOR ',inc,' NOT FOUND '
      call dadeb(31,'ERR DAINF ',1)
!
      return
      end
!
      subroutine dapac(inc)
      implicit none
      integer i,ic,illc,ilmc,inc,inoc,invc,ipoc
      double precision ccc
!     ************************
!
!     THIS SUBROUTINE PACKS THE INFORMATION IN THE SCRATCH VECTOR 1
!     INTO THE VECTOR INC.
!     INVERSE IS DAUNP.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
      ipoc = idapo(inc)

!
      ic = ipoc - 1
!
      do i=1,nmmax
        ccc = cc(i)
        if(dabs(ccc).lt.eps) goto 100
        ic = ic + 1
        cc(ic) = ccc
        i1(ic) = ie1(i)
        i2(ic) = ie2(i)
 100    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DAPAC '
        call dadeb(31,'ERR DAPAC ',1)
      endif
!
      return
      end
!
!
      subroutine dachk(ina,inoa,inva, inb,inob,invb, inc,inoc,invc)
      implicit none
      integer ierr,ina,inb,inc,inoa,inob,inoc,inva,invb,invc,invsum,lsw
!     *************************************************************
!
!     THIS SUBROUTINE CHECKS IF THE VECTORS A, B AND C
!     HAVE COMPATIBLE ATTRIBUTES
!
      parameter(lsw=1)
!
      if(lsw.eq.1) return
!
      ierr = 0
!
!     CASE OF A UNARY OPERATION
!     *************************
!
      if(inob.eq.-1.and.invb.eq.-1) then
        invsum = inva + invc
        if(invsum.eq.0) then
          if(inoa.gt.inoc) ierr = 1
        elseif(invsum.eq.1) then
          ierr = 1
        else
          if(inoa.gt.inoc.or.inva.gt.invc) ierr = 1
        endif
        if(ierr.eq.1) then
          write(6,*)'ERROR IN DACHK, ',ina,' AND ',inc,                 &
     &    ' ARE INCOMPATIBLE',inoa,inva,inoc,invc
          call dadeb(31,'ERR DACHK1',1)
        endif
!
!     CASE OF A BINARY OPERATION
!     **************************
!
      else
        invsum = inva + invb + invc
        if(invsum.eq.0) then
          if(inoa.gt.inoc.or.inob.gt.inoc) ierr = 1
        elseif(invsum.eq.1.or.invsum.eq.2) then
          ierr = 1
        else
          if(inoa.gt.inoc.or.inob.gt.inoc.or.                           &
     &    inva.gt.invc.or.invb.gt.invc) ierr = 1
        endif
        if(ierr.eq.1) then
          write(6,*)'ERROR IN DACHK, ',ina,',',inb,' AND ',inc,         &
     &    ' ARE INCOMPATIBLE'
          call dadeb(31,'ERR DACHK2',1)
        endif
      endif
!
      return
      end
!
      subroutine damch(iaa,ia)
      implicit none
      integer i,ia,iaa,illa,ilma,ino1,inoi,inv1,invi,ipoa
!     ************************
!
!     THIS SUBROUTINE CHECKS IF THE IA VECTORS IN THE MATRIX IA HAVE
!     IDENTICAL ATTRIBUTES.
!
      dimension iaa(*)
!
      call dainf(iaa(1),ino1,inv1,ipoa,ilma,illa)
!
      do i=2,ia
        call dainf(iaa(i),inoi,invi,ipoa,ilma,illa)
        if(ino1.ne.inoi.or.inv1.ne.invi) then
          write(6,*)'ERROR IN DAMCH, VECTORS ',iaa(1),' AND ',iaa(i),   &
     &    ' ARE INCOMPATIBLE '
          stop
        endif
      enddo
!
      return
      end
!
      subroutine dadcd(jj,ic1,ic2)
      implicit none
      integer i,ibase,ic1,ic2,isplit
!     ****************************
!
!     THIS SUBROUTINE CODES THE EXPONENTS IN JJ INTO THEIR DA CODES I1,I2.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(lnv)
      ibase = nomax + 1
      isplit = (nvmax+1)/2
      ic1 = 0
      ic2 = 0
!
      do i=nvmax,isplit+1,-1
        ic2 = ic2*ibase + jj(i)
      enddo
!
      do i=isplit,1,-1
        ic1 = ic1*ibase + jj(i)
      enddo
!
      return
      end
!
      subroutine dancd(ic1,ic2,jj)
      implicit none
      integer i,ibase,ic,ic1,ic2,isplit
      double precision x
!     ****************************
!
!     THIS SUBROUTINE ENCODES THE EXPONENTS IN JJ FROM THEIR DA CODES I1,I2.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(*)
      ibase = nomax + 1
      isplit = (nvmax+1)/2
!
      ic = ic1
      do i=1,isplit
        x  = dble(ic)/dble(ibase)
        ic = int(x+epsmac)
        jj(i) = nint(ibase*(x-ic))
      enddo
!
      ic = ic2
      do i=isplit+1,nvmax
        x  = dble(ic)/dble(ibase)
        ic = int(x+epsmac)
        jj(i) = nint(ibase*(x-ic))
      enddo
!
      do i=nvmax+1,lnv
        jj(i) = 0
      enddo
!
      return
      end

!ETIENNE
      subroutine datra(idif,ina,inc)
      implicit none
      integer i,ibase,ic,ider1,ider1s,ider2,ider2s,idif,iee,ifac,illa,  &
     &illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj
      double precision x,xdivi
!     ******************************
!
!     THIS SUBROUTINE COMPUTES THE PSEUDO DERIVATIVE WITH RESPECT TO VARIABLE I
!     OF THE VECTOR A AND STORES THE RESULT IN C.
!
!     dx^n/dx= x^(n-1)
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
!
!       CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!

      if(nomax.eq.1) then
        call dader(idif,ina,inc)
        return
      endif
      ibase = nomax + 1
!
      if(idif.gt.(nvmax+1)/2) then
        ider1  = 0
        ider1s = 0
        ider2  = idif-(nvmax+1)/2
        ider2s = 1
        do jj=1,ider2-1
          ider2s = ider2s*ibase
        enddo
        xdivi  = ider2s*ibase
      else
        ider1  = idif
        ider1s = 1
        do jj=1,ider1-1
          ider1s = ider1s*ibase
        enddo
        ider2  = 0
        ider2s = 0
        xdivi  = ider1s*ibase
      endif
!
      ibase = nomax+1
!
      ic = ipoc-1
!
      do i=ipoa,ipoa+illa-1
!
        if(ider1.eq.0) then
          iee = i2(i)
        else
          iee = i1(i)
        endif
!
        x = iee/xdivi
        ifac = int(ibase*(x-int(x+epsmac)+epsmac))
!
        if(ifac.eq.0) goto 100
!
!etienne      IFAC = INT(IBASE*(X-INT(X)+1.D-8))
!
!etienne      IF(IFAC.EQ.0) GOTO 100
!
        ic = ic + 1
        cc(ic) = cc(i)
        i1(ic) = i1(i) - ider1s
        i2(ic) = i2(i) - ider2s
!
 100    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DADTRA'
        call dadeb(111,'ERR DADTRA',1)
      endif
!
      return
      end

      subroutine etred(no1,nv1,ic1,ic2,no2,nv2,i11,i21)
      implicit none
      integer i,i11,i21,ic,ic1,ic2,no1,no2,nv1,nv2
!     ****************************
!
!     THIS SUBROUTINE CODES THE EXPONENTS IN JJ INTO THEIR DA CODES I1,I2.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(lnv)

      if(nv1.gt.lnv.or.nv2.gt.lnv) then
        write(6,*) ' ERROR IN RECODING '
        stop123
      endif
      call dehash(no1,nv1,ic1,ic2,jj)
      ic=0

      do i=nv1+1,nv2
        if(jj(i).gt.0) then
          i11=-1
          i21=0
          return
        endif
      enddo

      do i=1,nv2
        ic=ic+jj(i)
      enddo

      if(ic.gt.no2) then
        i11=-1
        i21=0
        return
      endif

      call hash(no2,nv2,jj,i11,i21)

!
      return
      end

      subroutine hash(no1,nv1,jj,ic1,ic2)
      implicit none
      integer i,ibase,ic1,ic2,isplit,no1,nv1
!     ****************************
!
!     THIS SUBROUTINE CODES THE EXPONENTS IN JJ INTO THEIR DA CODES I1,I2.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      integer jj(*)

      ibase = no1 + 1
      isplit = (nv1+1)/2
      ic1 = 0
      ic2 = 0
!
      do i=nv1,isplit+1,-1
        ic2 = ic2*ibase + jj(i)
      enddo
!
      do i=isplit,1,-1
        ic1 = ic1*ibase + jj(i)
      enddo
!
      return
      end
!
      subroutine dehash(no1,nv1,ic1,ic2,jj)
      implicit none
      integer i,ibase,ic,ic1,ic2,isplit,no1,nv1
      double precision x
!     ****************************
!
!     THIS SUBROUTINE ENCODES THE EXPONENTS IN JJ FROM THEIR DA CODES I1,I2.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jj(*)

      epsmac=1.e-7
      ibase = no1 + 1
      isplit = (nv1+1)/2
!
      ic = ic1
      do i=1,isplit
        x  = dble(ic)/dble(ibase)
        ic = int(x+epsmac)
        jj(i) = nint(ibase*(x-ic))
      enddo
!
      ic = ic2
      do i=isplit+1,nv1
        x  = dble(ic)/dble(ibase)
        ic = int(x+epsmac)
        jj(i) = nint(ibase*(x-ic))
      enddo
!
      return
      end

      subroutine daswap(j1,j2,inb)
      implicit none
      integer ia,ic,ic1,ic2,illb,ilmb,inb,inob,invb,ipob,j1,j2,jj,k1,k2
!     *************************
!
!     SWAP A DA VECTOR
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension jj(lnv)

      call dainf(inb,inob,invb,ipob,ilmb,illb)

      call daclr(1)
!

      do ia = ipob,ipob+illb-1

        call dehash(nomax,nvmax,i1(ia),i2(ia),jj)
        k1=jj(j1)
        k2=jj(j2)
        jj(j1)=k2
        jj(j2)=k1
        call hash(nomax,nvmax,jj,ic1,ic2)

        ic=ia1(ic1)+ia2(ic2)

        cc(ic) = cc(ia)
        i1(ic) = ic1
        i2(ic) = ic2
!
      enddo
!
!
      call dapac(inb)
      return
      end

      subroutine dagauss(ina,inb,nd2,anorm)
      implicit none
      integer i,ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,  &
     &ipoa,ipob,ja,jb,nd2
      double precision anorm,gau
!     ***************************
!
!     THIS SUBROUTINE COMPUTES THE NORM OF THE DA VECTOR A
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension ja(lnv),jb(lnv)

      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
      anorm = 0.d0

      do ia=ipoa,ipoa+illa-1
        do ib=ipob,ipob+illb-1
          call dancd(i1(ia),i2(ia),ja)
          call dancd(i1(ib),i2(ib),jb)
          gau=1.d0
          do i=1,nd2
            gau= facint(ja(i)+jb(i))*gau
          enddo
          anorm = anorm + cc(ia)*cc(ib)*gau
        enddo
      enddo

!
      return
      end

      subroutine daran(ina,cm,xran)
      implicit none
      integer i,illa,ilma,ina,inoa,inva,ipoa
      double precision bran,cm,xran
!     ************************
!
!     THIS SUBROUTINE FILLS THE DA VECTOR A WITH RANDOM ENTRIES.
!     FOR CM > 0, THE VECTOR IS FILLED WITH REALS,
!     FOR CM < 0, THE VECTOR IS FILLED WITH SINGLE DIGIT INTEGERS
!     ABS(CM) IS THE FILLING FACTOR
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
!
      if(inva.eq.0.or.nomax.eq.1) then
        do i=ipoa,ipoa+ilma-1
          if(cm.gt.0d0) then
            cc(i) = bran(xran)
            if(cc(i).gt.cm) cc(i) = 0d0
          elseif(cm.lt.0d0) then
            cc(i) = int(1+10*bran(xran))
            if(cc(i).gt.-1d1*cm) cc(i) = 0d0
          endif
        enddo
        idall(ina) = idalm(ina)
        return
      endif
!
      if(inoa.ne.nomax.or.inva.ne.nvmax) then
        write(6,*)'ERROR IN DARAN, ONLY VECTORS WITH NO = NOMAX AND'    &
     &  //' NV = NVMAX ALLOWED'
        call dadeb(31,'ERR DARAN1',1)
      endif
!
      call daclr(1)
!
      do i=1,nmmax
        if(cm.gt.0.d0) then
          cc(i) = bran(xran)
          if(cc(i).gt.cm) cc(i) = 0.d0
        elseif(cm.lt.0.d0) then
          cc(i) = int(1+10*bran(xran))
          if(cc(i).gt.-10.d0*cm) cc(i) = 0.d0
        else
          write(6,*)'ERROR IN ROUTINE DARAN'
          call dadeb(31,'ERR DARAN2',1)
        endif
      enddo
!
      call dapac(ina)
!
      return
      end
!
      double precision function bran(xran)
      implicit none
      double precision xran
!     ************************************
!
!     VERY SIMPLE RANDOM NUMBER GENERATOR
!
      xran = xran + 10.d0
      if(xran.gt.1.d4) xran = xran - 9999.12345
      bran = dabs(sin(xran))
      bran = 10*bran
      bran = bran - int(bran)
!      IF(BRAN.LT. .1D0) BRAN = BRAN + .1D0
!
      return
      end
!
!

      subroutine danorm2(ina,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
!     ******************************
!
!     THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
!     THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
!     THE DA VECTOR DENOTED WITH THE INTEGER E.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!

      if(ina.eq.inc) then
        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',inoc,invc)
        call danorm2t(ina,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call danorm2t(ina,inc)
      endif

      return
      end

      subroutine danorm2t(ina,inb)
      implicit none
      integer ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,    &
     &ipoa,ipob
!     ******************************
!
!     THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
!     THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
!     THE DA VECTOR DENOTED WITH THE INTEGER E.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INB,INOB,INVB)
!
!
      ib = ipob - 1
!
      do ia=ipoa,ipoa+illa-1
!
        if(ieo(ia1(i1(ia))+ia2(i2(ia))).gt.nocut) goto 100
        ib = ib + 1
        cc(ib) = cc(ia)**2
        i1(ib) = i1(ia)
        i2(ib) = i2(ia)
!
 100    continue
      enddo
!
      idall(inb) = ib-ipob+1
      if(idall(inb).gt.idalm(inb)) then
        write(6,*)'ERROR IN DANORM'
        call dadeb(31,'ERR DANOR1',1)
      endif
!
      return
      end

      subroutine danormr(ina,inc)
      implicit none
      integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc
!     ******************************
!
!     THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
!     THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
!     THE DA VECTOR DENOTED WITH THE INTEGER E.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!

      if(ina.eq.inc) then
        call dainf(inc,inoc,invc,ipoc,ilmc,illc)
        incc=0
        call daall(incc,1,'$$DAJUNK$$',inoc,invc)
        call danormrt(ina,incc)
        call dacop(incc,inc)
        call dadal(incc,1)
      else
        call danormrt(ina,inc)
      endif

      return
      end

      subroutine danormrt(ina,inb)
      implicit none
      integer ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,    &
     &ipoa,ipob
!     ******************************
!
!     THIS SUBROUTINE MULTIPLIES THE DA VECTOR DENOTED BY THE
!     THE INTEGER A WITH THE CONSTANT C AND STORES THE RESULT IN
!     THE DA VECTOR DENOTED WITH THE INTEGER E.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inb,inob,invb,ipob,ilmb,illb)
!
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INB,INOB,INVB)
!
!
      ib = ipob - 1
!
      do ia=ipoa,ipoa+illa-1
!
        if(ieo(ia1(i1(ia))+ia2(i2(ia))).gt.nocut) goto 100
        ib = ib + 1
        cc(ib) = dsqrt(cc(ia))
        i1(ib) = i1(ia)
        i2(ib) = i2(ia)
!
 100    continue
      enddo
!
      idall(inb) = ib-ipob+1
      if(idall(inb).gt.idalm(inb)) then
        write(6,*)'ERROR IN DANORM '
        call dadeb(31,'ERR DANOR2',1)
      endif
!
      return
      end
      subroutine dakey(c)
      implicit none
      character c*(*)
!

      return
!
      end
! ANFANG UNTERPROGRAMM
      subroutine dapri6(ina,result,ien,i56)
      implicit none
      integer i,i56,ien,ihp,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,   &
     &j
      double precision result
!     *************************************
!
!     THIS SUBROUTINE IS FOR REDUCED STORAGE DA VERSION JULY 91
!     RESULT CONTAINS THE (IEN-1) TH DERIVATIVE TO ENERGY AFTER EXECUTION
!     I56 SAYS WHETHER THE 5TH OR THE 6TH COORDINATE IS THE ENERGY
!     AND MUST HAVE THE VALUE 5 OR 6 ACCORDINGLY
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
      dimension j(lnv)
      result=0.
      if(ina.lt.1.or.ina.gt.nda) then
        print*,'ERROR IN DAPRI6, INA = ',ina
        stop
      endif
      inoa=idano(ina)
      inva=idanv(ina)
      ipoa=idapo(ina)
      ilma=idalm(ina)
      illa=idall(ina)
      iout=0
      if(nomax.eq.1) then
        do i=1,illa
          if(ien.eq.1) then
            if(i-1.ne.0) goto 90
            result=cc(ipoa+i-1)
            return
          endif
          if(ien.eq.2) then
            if(i-1.ne.i56) goto 90
            result=cc(ipoa+i-1)
            return
          endif
 90       continue
        enddo
      else
        do ioa=0,inoa
          do ii=ipoa,ipoa+illa-1
            if(ieo(ia1(i1(ii))+ia2(i2(ii))).ne.ioa) goto 100
            iout=iout+1
            call dancd(i1(ii),i2(ii),j)
            if(i56.eq.6) then
              do ihp=1,5
                if(j(ihp).ne.0) goto 100
              enddo
              if(j(6).eq.(ien-1)) then
                result=cc(ii)
                return
              endif
            else if(i56.eq.5) then
              do ihp=1,4
                if(j(ihp).ne.0) goto 100
              enddo
              if(j(5).eq.(ien-1)) then
                result=cc(ii)
                return
              endif
            else if(i56.eq.4) then
              do ihp=1,3
                if(j(ihp).ne.0) goto 100
              enddo
              if(j(4).eq.(ien-1)) then
                result=cc(ii)
                return
              endif
            endif
 100        continue
          enddo
        enddo
      endif
      return
      end
! ANFANG UNTERPROGRAMM
      subroutine darea6(ina,zfeld,i56)
      implicit none
      integer i,i56,ic,ii1,ii2,iin,illa,ilma,ina,inoa,inva,io,io1,ip,   &
     &ipoa,iwarin,iwarno,iwarnv,j
      double precision zfeld
!     *************************************
!
!     THIS SUBROUTINE IS FOR REDUCED STORAGE DA VERSION JULY 91
!     I56 SAYS WHETHER THE 5TH OR THE 6TH COORDINATE IS THE ENERGY
!     AND MUST HAVE THE VALUE 5 OR 6 ACCORDINGLY
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      character daname(lda)*10
      common / daname / daname
      dimension zfeld(100)
!-----------------------------------------------------------------------------3
      dimension j(lnv)
      if(ina.lt.1.or.ina.gt.nda) then
        print*,'ERROR IN DAREA6, INA = ',ina
        stop
      endif
      inoa=idano(ina)
      inva=idanv(ina)
      ipoa=idapo(ina)
      ilma=idalm(ina)
      illa=idall(ina)
      do i=1,lnv
        j(i)=0
      enddo
      call daclr(1)
      ic=0
      iwarno=0
      iwarnv=0
      iwarin=0
      iin=0
      if(nomax.eq.1) then
        do i=1,illa
          if (i-1.eq.0) then
            cc(ipoa+i-1)=zfeld(1)
          else if (i-1.eq.i56) then
            cc(ipoa+i-1)=zfeld(2)
          endif
        enddo
        return
      endif
      do ip=1,inva
        j(ip)=0
      enddo
      io=0
  10  continue
      iin=iin+1
      io1=0
      do i=1,inva
        io1=io1+j(i)
      enddo
      if(io1.ne.io) then
        if(iwarnv.eq.0) print*,'WARNING IN DAREA6, FILE ',              &
     &  'CONTAINS MORE VARIABLES THAN VECTOR'
        iwarnv = 1
        goto 10
      endif
      if(io.gt.inoa) then
        iwarno = 1
        goto 10
      endif
      ic = ic + 1
      call dadcd(j,ii1,ii2)
      ic = ia1(ii1) + ia2(ii2)
      cc(ic) = zfeld(io+1)
      j(i56)=j(i56)+1
      io=io+1
      if (io.gt.inoa) goto 20
      goto 10
  20  continue
      call dapac(ina)
      return
      end
! ANFANG FUNKTION
      double precision function dare(ina)
      implicit none
      integer ii,illa,ilma,ina,inoa,inva,ioa,ipoa,j,jj
!     ***********************************
!     NEW VERSION OF DARE, AUGUST 1992
!     SUPPOSED TO TREAT THE 0TH COMPONENT ACCURATELY
!
!     30.10 1997 E.Mcintosh & F.Schmidt
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      dimension j(lnv)
!-----------------------------------------------------------------------------9
!
!      CALL DAINF(INA,INOA,INVA,IPOA,ILMA,ILLA)
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)

!FRS 30.10.1997
      if(nomax.eq.1) then
        dare = cc(ipoa)
        return
      endif
!FRS 30.10.1997
!FRS March 1997
!      IF(NOMAX.EQ.1) goto 110
!FRS March 1997

      ioa = 0
      do ii=ipoa,ipoa+illa-1
        if(ieo(ia1(i1(ii))+ia2(i2(ii))).ne.ioa) goto 100
        call dancd(i1(ii),i2(ii),j)
        do jj=1,inva
          if(j(jj).ne.0) goto 100
        enddo
! 110    continue
        dare = cc(ipoa)
        return
 100    continue
      enddo
      dare = 0d0
      return
      end
! ANFANG UNTERPROGRAMM
      subroutine daprimax(ina,iunit)
      implicit none
      integer i,ii,iii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit,j
!     ***************************
!
!     THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!-----------------------------------------------------------------------------9
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      dimension j(lnv)
!
      if(ina.lt.1.or.ina.gt.nda) then
        print*,'ERROR IN DAPRI, INA = ',ina
        stop
      endif
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
!
      iout = 0
      ioa = 0

      if(inva.eq.0) then
        do i = ipoa,ipoa+illa-1
          write(iunit,'(I6,2X,G20.14)') i-ipoa, cc(i)
        enddo
      elseif(nomax.eq.1) then
        do i=1,illa
          iout=iout+1
          if(i.ne.1) then
            j(i-1)=1
            ioa=1
          endif
          write(iunit,'(I6,2X,G20.14,I5,4X,18(2I2,1X))')                &
     &    iout,cc(ipoa+i-1),ioa,(j(iii),iii=1,nvmax)
          write(iunit,*) cc(ipoa+i-1)
        enddo
      else
        iout = 0
        do ioa = 0,inoa
          do ii=ipoa,ipoa+illa-1
            if(ieo(ia1(i1(ii))+ia2(i2(ii))).ne.ioa) goto 100
            call dancd(i1(ii),i2(ii),j)
!ETIENNE
            if(abs(cc(ii)).gt.eps) then
!ETIENNE
              iout = iout+1
              write(iunit,'(I6,2X,G20.14,I5,4X,18(2I2,1X))')            &
     &        iout,cc(ii),ioa,(j(iii),iii=1,nvmax)
!ETIENNE
              write(iunit,*) cc(ii)
            endif
!ETIENNE
!
 100        continue
          enddo
        enddo
      endif
!

!     WRITE(IUNIT,'(A)') '                                      '
!
      return
      end
!FF

!  unknown stuff
      subroutine damono(ina,jd,cfac,istart,inc)
      implicit none
      integer ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,    &
     &ipoa,ipoc,istart,jd
      double precision cfac
!     *****************************
!
!     THIS SUBROUTINE RETURNS THE MONOMIALS ONE BY ONE
!     OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE
!     RESULT IN C
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      dimension jd(*)
!
      if(ina.eq.inc) then
        write(6,*) ' USE DIFFERENT POWER SERIES IN DAMONO '
        stop999
      endif
      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
      if(istart.eq.0) then
        istart=illa
        return
      endif
!
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!
      ic = ipoc - 1
!
      ia=ipoa+istart-1
!
      call dancd(i1(ia),i2(ia),jd)

!
      ic = ic + 1
      cc(ic) = cc(ia)
      cfac=cc(ia)
      i1(ic) = i1(ia)
      i2(ic) = i2(ia)
!
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DAMONO'
        call dadeb(31,'ERR DAMONO',1)
      endif
!
      return
      end
!
!

      subroutine dacycle(ina,ipresent,value,j,illa)
      implicit none
      integer i,ii,illa,ilma,ina,inoa,inva,iout,ipoa,ipresent,j
      double precision value
!     ***************************
!
!     THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      dimension j(lnv)
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
!
      if(ina.lt.1.or.ina.gt.nda) then
        write(6,*)'ERROR IN DAPRI, INA = ',ina
        stop
      endif
!
      if(ina.eq.0) then
        value=0.d0
        illa=0
        do i=1,lnv
          j(i)=0
        enddo
        return
      endif

      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
      iout = 0
      ipresent=1+ipresent
      if(ipresent.gt.illa) then
        ipresent=1
      endif
      ii=ipresent+ipoa-1
      call dancd(i1(ii),i2(ii),j)
      value=cc(ii)
      return

      end
      subroutine daorder(ina,iunit,jx,invo,nchop)
      implicit none
      integer i,ic,ii,ii1,ii2,iin,illa,ilma,ina,inoa,inva,invo,io,io1,  &
     &ipoa,iunit,iwarin,iwarno,iwarnv,j,jh,jt,jx,nchop
      double precision c
!     ***************************
!
!     THIS SUBROUTINE READS THE DA VECTOR INA FROM UNIT IUNIT.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
      character daname(lda)*10
      common / daname / daname
!-----------------------------------------------------------------------------3
!
      character c10*10
      dimension j(lnv),jx(lnv),jt(lnv)
!
      if(ina.lt.1.or.ina.gt.nda) then
        write(6,*)'ERROR IN DAPRI, INA = ',ina
        stop
      endif
!
      inoa = idano(ina)
      inva = idanv(ina)
      ipoa = idapo(ina)
      ilma = idalm(ina)
      illa = idall(ina)
!
      do i=1,lnv
        jt(i)=0
        j(i) = 0
      enddo
!
      call daclr(1)
!
      ic = 0
!
      iwarno = 0
      iwarnv = 0
      iwarin = 0
!
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
      read(iunit,'(A10)') c10
!
      iin = 0
!
  10  continue
      iin = iin + 1
      read(iunit,'(I6,2X,G20.14,I5,4X,18(2I2,1X))')                     &
     &ii,c,io,(jt(i),i=1,invo)
!
      if(ii.eq.0) goto 20
!etienne
      read(iunit,*) c
      do jh=1,invo
        j(jh)=jt(jx(jh))
      enddo
      do jh=nchop+1,inva
        j(jh)=0
      enddo
!etienne
      io1 = 0
      do i=1,inva
        io1 = io1 + j(i)
      enddo
!
      ic = ic + 1
      call dadcd(j,ii1,ii2)
      ic = ia1(ii1) + ia2(ii2)
      cc(ic) = c
      goto 10
!
  20  continue
!
      call dapac(ina)
!
      return
      end
!
!ETIENNE
      subroutine datrash(idif,ina,inc)
      implicit none
      integer i,ibase,ic,ider1,ider1s,ider2,ider2s,idif,ikil1,ikil2,    &
     &illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj
      double precision xdivi
!     ******************************
!
!     THIS SUBROUTINE COMPUTES THE DERIVATIVE WITH RESPECT TO VARIABLE I
!     OF THE VECTOR A AND STORES THE RESULT IN C.
!
!-----------------------------------------------------------------------------1
      include "TPSALib_prm.f"
!
      integer jx(lnv)

!      call daclr(1)
!      call dacop(ina,1)
!      call dapac(ina)

      call dainf(ina,inoa,inva,ipoa,ilma,illa)
      call dainf(inc,inoc,invc,ipoc,ilmc,illc)
!
!
!      CALL DACHK(INA,INOA,INVA, '          ',-1,-1,INC,INOC,INVC)
!
      ibase = nomax + 1
!
      if(idif.gt.(nvmax+1)/2) then
        ider1  = 0
        ider1s = 0
        ider2  = idif-(nvmax+1)/2
        ider2s = 1
        do jj=1,ider2-1
          ider2s = ider2s*ibase
        enddo
        xdivi  = ider2s*ibase
      else
        ider1  = idif
        ider1s = 1
        do jj=1,ider1-1
          ider1s = ider1s*ibase
        enddo
        ider2  = 0
        ider2s = 0
        xdivi  = ider1s*ibase
      endif
!
      ibase = nomax+1
!
      ic = ipoc-1
!
      do i=ipoa,ipoa+illa-1
!
        call dancd(i1(i),i2(i),jx)

        ikil1=0
        ikil2=0
        if(idif.gt.(nvmax+1)/2) then
          ikil2=jx(idif)
        else
          ikil1=jx(idif)
        endif
!
!      X = IEE/XDIVI
!etienne      IFAC = INT(IBASE*(X-INT(X)+1.D-8))
!
!etienne      IF(IFAC.EQ.0) GOTO 100
!
        ic = ic + 1
        cc(ic) = cc(i)
        i1(ic) = i1(i) - ikil1*ider1s
        i2(ic) = i2(i) - ikil2*ider2s
!
! 100    continue
      enddo
!
      idall(inc) = ic - ipoc + 1
      if(idall(inc).gt.idalm(inc)) then
        write(6,*)'ERROR IN DATRASH '
        call dadeb(111,'ERR DATRAS',1)
      endif
!
      return
      end
      integer function  mypause(i)
      implicit none
! Replaces obsolescent feature pause
      integer i
      !write (*,'(A,i6)',ADVANCE='NO') ' PAUSE: ',i
      write (*,'(A,i6)') ' PAUSE: ',i
      read(*,*)
      mypause=i
      end function mypause
