C	although the name is "trapezfi", options for "Simpson Filon" have 
C	been included

C	Sven Ivansson  1990

C****************************************************************
C
C     TRAPEZFI
C	integral of nhyd*li complex*16 functions (obtained by 
C	calling the subroutines dhcomp and rcomp) of a COMPLEX*16 
C	argument by trapezoidal integration
C	note: 
C	  *  options for Filon integration 
C	       difference simple/double Filon only for "end point" terms
C	       "trapezoidal non-Filon" :  ifilon = 0
C	       "trapezoidal Filon"     :  ifilon = 1,2
C	       "Simpson Filon"         :  ifilon = -1,-2
C	       "Simpson non-Filon"     :  ifilon = -3
C	     - simple (ifilon=+-1)
C	       " f(up) = g(up)*cdexp(gamma*up) "  where g is smooth 
C	       f    is one of the integrands obtained with calls to the 
C	            formal subroutines dhcomp and rcomp
C	            f(up) = dhfakt(up)*rfakt(up)
C	     - double (ifilon=+-2)
C	       " f(up) = gpos(up)*cdexp(gamma*up) + gneg(up)*cdexp(-gamma*up)
C	               =      fpos(up)           +      fneg(up)             "
C	       where gpos and gneg are smooth 
C	       f    is one of the integrands obtained with calls to the 
C	            formal subroutines dhcomp and rcomp
C	            f(up) = dhfakt(up)*rfakt(up)
C	       fpos is one of the integrands obtained with calls to the 
C	            formal subroutines dhcomp and rcomppos
C	            fpos(up) = dhfakt(up)*rfaktpos(up)
C	       fneg is thus obtained as  fneg(up) = f(up) - fpos(up) = 
C	            dhfakt(up) * ( rfakt(up) - rfaktpos(up) )
C	  *  integration over a polygon in the complex "up" plane rather 
C	       than over a line segment on the real "x" axis
C	  *  "matrix integrand" instead of "scalar integrand"
C
C****************************************************************

	SUBROUTINE TRAPEZFI ( NHYD,LI, DHCOMP,RCOMP, IFILON,RCOMPPOS,
     $	GAMMA, NPU,UPHOERN, HMAX, NHYDLI,VALUE, NFCALLS, 
     $	DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CWORK )
C	" call dhcomp (up, rlndet,cdet, dhfakt) "  will give rlndet & cdet 
C	     and dhfakt(1:nhyd)(up)
C	" call  rcomp (up,  rfakt) "  will give  rfakt(1:li)(up)
C	ifilon  is 0 for "trapezoidal non-Filon", 1 for "trapezoidal simple 
C	     Filon", 2 for "trapezoidal double Filon", -1 for "Simpson simple 
C	     Filon", -2 for "Simpson double Filon", -3 for "Simpson non-Filon"
C	" call rcomppos (up, rfakt) "  will give  rfaktpos(1:li)(up), which 
C	     is needed when ifilon=+-2
C	gamma(1:li)  are the exponents to be utilized in the Filon method
C	uphoern(1:npu)  are the corners of the integration-path polygon in 
C	     the complex "up" plane
C	hmax  is the maximum allowed integration step along the integration 
C	     path
C	nhydli  is nhyd*li (input)
C	value(1:nhydli)  is input/output (OBS), 
C	     the integral of  dhfakt(ih)*rfakt(ir)  is added to  
C	     value(ihr)  where  ihr = (ih-1)*li + ir , 
C	     ih=1,2,...,nhyd  ir=1,2,...,li
C	nfcalls  is output which counts the "function calls"
C	work-space:
C	  dhfakt(1:nhyd)
C	  rfakt(1:li)
C	  icworklim0,icworkdim
C	  cwork(1:icworkdim)   "containing cwork(icworklim0+1,icworkdim)" 
C	                       < cwork(1:icworklim0) is not accessed here >
C	     rfaktpos(1:li)
C	     bhyprst(1:li)
C	     chyprst(1:li)
C	     stdh(1:nhyd)
C	     str(1:li)
C	     strpos(1:li)

	IMPLICIT NONE

	INTEGER NHYD,LI, IFILON, NPU, NHYDLI
	REAL*8 HMAX
	COMPLEX*16 GAMMA(*), UPHOERN(*)
	EXTERNAL DHCOMP,RCOMP, RCOMPPOS

	INTEGER NFCALLS
	COMPLEX*16 VALUE(*)

	INTEGER ICWORKLIM0,ICWORKDIM
	COMPLEX*16 DHFAKT(*),RFAKT(*), CWORK(*)

	INTEGER LL, I, JL, IH,IR, IHR, 
     $	IADR0_RFAKTPOS, IADR0_BHYPRST, IADR0_CHYPRST, IADR0_STDH, 
     $	IADR0_STR, IADR0_STRPOS
	REAL*8 HCUR, SLASK, RLNDET
	COMPLEX*16 UPH0, UPKLIV,UPRIKT,UPDELTA, UP, CDF, 
     $	BHYPR,CHYPR,SHYPR, XGP,XGN,XXX,YYY, CSLASK, CDET, dhfaktsum
	LOGICAL CNOW

	IADR0_RFAKTPOS = ICWORKLIM0
	IADR0_BHYPRST = IADR0_RFAKTPOS + LI
	IADR0_CHYPRST = IADR0_BHYPRST + LI
	IADR0_STDH = IADR0_CHYPRST + LI
	IADR0_STR = IADR0_STDH + NHYD
	IADR0_STRPOS = IADR0_STR + LI
	IF ( ( IADR0_STRPOS + LI ) .GT. ICWORKDIM ) STOP 160

	NFCALLS = 1
	UP = UPHOERN(NPU)
	CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	DO IH=1,NHYD
	  CWORK(IADR0_STDH+IH) = DHFAKT(IH)
	END DO
	CALL RCOMP (UP, RFAKT)
	DO IR=1,LI
	  CWORK(IADR0_STR+IR) = RFAKT(IR)
	END DO
	IF (IABS(IFILON).EQ.2) CALL RCOMPPOS ( UP, 
     $	CWORK(IADR0_STRPOS+1) )

	DO JL = NPU, 2, -1

	  UPH0 = UPHOERN(JL-1)
	  UPKLIV = UPHOERN(JL) - UPH0
	  SLASK = CDABS(UPKLIV)
	  IF (SLASK.GT.0.0D0) THEN
	    UPRIKT = UPKLIV/SLASK
	  ELSE
	    UPRIKT = (0.0D0,0.0D0)
	  END IF
	  LL = 1 + SLASK/HMAX
	  IF (IFILON.LT.0) THEN
	    IF ( LL .NE. 2*(LL/2) ) LL = LL + 1  ! even number of subintervals
	  END IF
	  HCUR = SLASK/LL
	  UPDELTA = HCUR*UPRIKT

	  if (npu.ne.2) stop 341441 
	  open(UNIT=44,FILE='trapezfi.xspec',STATUS='NEW')
	  open(UNIT=45,FILE='trapezfi.yspec',STATUS='NEW')
	  open(UNIT=46,FILE='trapezfi.zspec',STATUS='NEW')
	  write(44,*) 'up complex', LL+1,' 1 '
	  write(45,*) 'dhfakt complex', LL+1,NHYD
	  write(46,*) 'dhfaktsum complex', LL+1,1
	  write(44,*) dreal(up),dimag(up)
	  write(45,*) (dreal(dhfakt(ih)),dimag(dhfakt(ih)),ih=1,nhyd)
	  dhfaktsum = (0.0d0,0.0d0)
	  do ih = 1, nhyd
	    dhfaktsum = dhfaktsum + dhfakt(ih)
	  end do
	  write(46,*) dreal(dhfaktsum),dimag(dhfaktsum)

	  NFCALLS = NFCALLS + 1
	  CALL DHCOMP (UPH0, RLNDET,CDET, DHFAKT)
	  CALL RCOMP (UPH0, RFAKT)

	  IF (IFILON.EQ.0) THEN
	    CSLASK = 0.5D0*UPDELTA
	    IHR = 0
	    DO IH=1,NHYD
	      CDF = CSLASK * CWORK(IADR0_STDH+IH)
	      DO IR=1,LI
	        IHR = IHR+1
	        VALUE(IHR) = VALUE(IHR) + CDF * CWORK(IADR0_STR+IR)
	      END DO
	    END DO
	    IHR = 0
	    DO IH=1,NHYD
	      CDF = CSLASK*DHFAKT(IH)
	      DO IR=1,LI
	        IHR = IHR+1
	        VALUE(IHR) = VALUE(IHR) + CDF*RFAKT(IR)
	      END DO
	    END DO
	  ELSE IF (IFILON.EQ.1) THEN
	    DO IR=1,LI
	      CSLASK = GAMMA(IR)*UPDELTA
	      CALL HYPREST (CSLASK, CHYPR,SHYPR)
	      CHYPR = UPDELTA * CHYPR
	      SHYPR = UPDELTA * SHYPR
	      CWORK(IADR0_CHYPRST+IR) = CHYPR
	      XGP = 0.5D0 * (CHYPR+SHYPR)
	      XGN = 0.5D0 * (CHYPR-SHYPR)
	      IHR = IR
	      CDF = XGN * CWORK(IADR0_STR+IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + CWORK(IADR0_STDH+IH) * CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGP*RFAKT(IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + DHFAKT(IH)*CDF
	        IHR = IHR + LI
	      END DO
	    END DO
	  ELSE IF (IFILON.EQ.2) THEN
	    CALL RCOMPPOS ( UPH0, CWORK(IADR0_RFAKTPOS+1) )
	    DO IR=1,LI
	      CSLASK = GAMMA(IR)*UPDELTA
	      CALL HYPREST (CSLASK, CHYPR,SHYPR)
	      CHYPR = UPDELTA * CHYPR
	      SHYPR = UPDELTA * SHYPR
	      CWORK(IADR0_CHYPRST+IR) = CHYPR
	      XGP = 0.5D0 * (CHYPR+SHYPR)
	      XGN = 0.5D0 * (CHYPR-SHYPR)
	      XXX = RFAKT(IR) - CWORK(IADR0_RFAKTPOS+IR)
	      YYY = CWORK(IADR0_STR+IR) - CWORK(IADR0_STRPOS+IR)
	      IHR = IR
	      CDF = XGN * CWORK(IADR0_STRPOS+IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + CWORK(IADR0_STDH+IH) * CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGP*YYY
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + CWORK(IADR0_STDH+IH) * CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGP * CWORK(IADR0_RFAKTPOS+IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + DHFAKT(IH)*CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGN*XXX
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + DHFAKT(IH)*CDF
	        IHR = IHR + LI
	      END DO
	    END DO
	  ELSE IF (IFILON.EQ.-1) THEN
	    DO IR=1,LI
	      CSLASK = GAMMA(IR)*UPDELTA
	      CALL HYPREST2 (CSLASK, BHYPR,CHYPR,SHYPR)
	      CHYPR = UPDELTA * (BHYPR+CHYPR)
	      BHYPR = UPDELTA * BHYPR
	      SHYPR = UPDELTA * SHYPR
	      CWORK(IADR0_BHYPRST+IR) = BHYPR
	      CWORK(IADR0_CHYPRST+IR) = CHYPR
	      XGP = 0.5D0 * (CHYPR+SHYPR)
	      XGN = 0.5D0 * (CHYPR-SHYPR)
	      IHR = IR
	      CDF = XGN * CWORK(IADR0_STR+IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + CWORK(IADR0_STDH+IH) * CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGP*RFAKT(IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + DHFAKT(IH)*CDF
	        IHR = IHR + LI
	      END DO
	    END DO
	  ELSE IF (IFILON.EQ.-2) THEN
	    CALL RCOMPPOS ( UPH0, CWORK(IADR0_RFAKTPOS+1) )
	    DO IR=1,LI
	      CSLASK = GAMMA(IR)*UPDELTA
	      CALL HYPREST2 (CSLASK, BHYPR,CHYPR,SHYPR)
	      CHYPR = UPDELTA * (BHYPR+CHYPR)
	      BHYPR = UPDELTA * BHYPR
	      SHYPR = UPDELTA * SHYPR
	      CWORK(IADR0_BHYPRST+IR) = BHYPR
	      CWORK(IADR0_CHYPRST+IR) = CHYPR
	      XGP = 0.5D0 * (CHYPR+SHYPR)
	      XGN = 0.5D0 * (CHYPR-SHYPR)
	      XXX = RFAKT(IR) - CWORK(IADR0_RFAKTPOS+IR)
	      YYY = CWORK(IADR0_STR+IR) - CWORK(IADR0_STRPOS+IR)
	      IHR = IR
	      CDF = XGN * CWORK(IADR0_STRPOS+IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + CWORK(IADR0_STDH+IH) * CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGP*YYY
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + CWORK(IADR0_STDH+IH) * CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGP * CWORK(IADR0_RFAKTPOS+IR)
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + DHFAKT(IH)*CDF
	        IHR = IHR + LI
	      END DO
	      IHR = IR
	      CDF = XGN*XXX
	      DO IH=1,NHYD
	        VALUE(IHR) = VALUE(IHR) + DHFAKT(IH)*CDF
	        IHR = IHR + LI
	      END DO
	    END DO
	  ELSE IF (IFILON.EQ.-3) THEN
	    BHYPR = UPDELTA * (4.0D0/3)
	    CHYPR = UPDELTA * (2.0D0/3)
	    XGP = 0.5D0 * CHYPR
	    XGN = 0.5D0 * CHYPR
	    IHR = 0
	    DO IH=1,NHYD
	      CDF = XGN * CWORK(IADR0_STDH+IH)
	      DO IR=1,LI
	        IHR = IHR+1
	        VALUE(IHR) = VALUE(IHR) + CDF * CWORK(IADR0_STR+IR)
	      END DO
	    END DO
	    IHR = 0
	    DO IH=1,NHYD
	      CDF = XGP*DHFAKT(IH)
	      DO IR=1,LI
	        IHR = IHR+1
	        VALUE(IHR) = VALUE(IHR) + CDF*RFAKT(IR)
	      END DO
	    END DO
	  ELSE
	    STOP 670
	  END IF

	  DO IH=1,NHYD
	    CWORK(IADR0_STDH+IH) = DHFAKT(IH)
	  END DO
	  DO IR=1,LI
	    CWORK(IADR0_STR+IR) = RFAKT(IR)
	  END DO
	  IF (IABS(IFILON).EQ.2) THEN
	    DO IR=1,LI
	      CWORK(IADR0_STRPOS+IR) = CWORK(IADR0_RFAKTPOS+IR)
	    END DO
	  END IF

	  NFCALLS = NFCALLS + LL-1

	  IF (IFILON.EQ.0) THEN

	    DO I = LL-1, 1, -1

	      UP = UPH0 + (I*HCUR)*UPRIKT
	      CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	      write(44,*) dreal(up),dimag(up)
	      write(45,*) (dreal(dhfakt(ih)),dimag(dhfakt(ih)),ih=1,nhyd)
	      dhfaktsum = (0.0d0,0.0d0)
	      do ih = 1, nhyd
	        dhfaktsum = dhfaktsum + dhfakt(ih)
	      end do
	      write(46,*) dreal(dhfaktsum),dimag(dhfaktsum)
	      CALL RCOMP (UP, RFAKT)
	      IHR = 0
	      DO IH=1,NHYD
	        CDF = UPDELTA*DHFAKT(IH)
	        DO IR=1,LI
	          IHR = IHR+1
	          VALUE(IHR) = VALUE(IHR) + CDF*RFAKT(IR)
	        END DO
	      END DO

	    END DO

	  ELSE IF (IFILON.GT.0) THEN

	    DO I = LL-1, 1, -1

	      UP = UPH0 + (I*HCUR)*UPRIKT
	      CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	      write(44,*) dreal(up),dimag(up)
	      write(45,*) (dreal(dhfakt(ih)),dimag(dhfakt(ih)),ih=1,nhyd)
	      dhfaktsum = (0.0d0,0.0d0)
	      do ih = 1, nhyd
	        dhfaktsum = dhfaktsum + dhfakt(ih)
	      end do
	      write(46,*) dreal(dhfaktsum),dimag(dhfaktsum)
	      CALL RCOMP (UP, RFAKT)
	      DO IR=1,LI
	        IHR = IR
	        CDF = CWORK(IADR0_CHYPRST+IR) * RFAKT(IR)
	        DO IH=1,NHYD
	          VALUE(IHR) = VALUE(IHR) + CDF*DHFAKT(IH)
	          IHR = IHR + LI
	        END DO
	      END DO

	    END DO

	  ELSE IF (IFILON.GT.-3) THEN

	    CNOW = .FALSE.

	    DO I = LL-1, 1, -1

	      UP = UPH0 + (I*HCUR)*UPRIKT
	      CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	      write(44,*) dreal(up),dimag(up)
	      write(45,*) (dreal(dhfakt(ih)),dimag(dhfakt(ih)),ih=1,nhyd)
	      dhfaktsum = (0.0d0,0.0d0)
	      do ih = 1, nhyd
	        dhfaktsum = dhfaktsum + dhfakt(ih)
	      end do
	      write(46,*) dreal(dhfaktsum),dimag(dhfaktsum)
	      CALL RCOMP (UP, RFAKT)
	      IF (CNOW) THEN
	        DO IR=1,LI
	          IHR = IR
	          CDF = CWORK(IADR0_CHYPRST+IR) * RFAKT(IR)
	          DO IH=1,NHYD
	            VALUE(IHR) = VALUE(IHR) + CDF*DHFAKT(IH)
	            IHR = IHR + LI
	          END DO
	        END DO
	      ELSE
	        DO IR=1,LI
	          IHR = IR
	          CDF = CWORK(IADR0_BHYPRST+IR) * RFAKT(IR)
	          DO IH=1,NHYD
	            VALUE(IHR) = VALUE(IHR) + CDF*DHFAKT(IH)
	            IHR = IHR + LI
	          END DO
	        END DO
	      END IF

	      CNOW = .NOT.CNOW

	    END DO

	  ELSE

	    CNOW = .FALSE.

	    DO I = LL-1, 1, -1

	      UP = UPH0 + (I*HCUR)*UPRIKT
	      CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	      write(44,*) dreal(up),dimag(up)
	      write(45,*) (dreal(dhfakt(ih)),dimag(dhfakt(ih)),ih=1,nhyd)
	      dhfaktsum = (0.0d0,0.0d0)
	      do ih = 1, nhyd
	        dhfaktsum = dhfaktsum + dhfakt(ih)
	      end do
	      write(46,*) dreal(dhfaktsum),dimag(dhfaktsum)
	      CALL RCOMP (UP, RFAKT)
	      IF (CNOW) THEN
	        IHR=0
	        DO IH=1,NHYD
	          CDF = CHYPR*DHFAKT(IH)
	          DO IR=1,LI
	            IHR=IHR+1
	            VALUE(IHR) = VALUE(IHR) + CDF*RFAKT(IR)
	          END DO
	        END DO
	      ELSE
	        IHR=0
	        DO IH=1,NHYD
	          CDF = BHYPR*DHFAKT(IH)
	          DO IR=1,LI
	            IHR=IHR+1
	            VALUE(IHR) = VALUE(IHR) + CDF*RFAKT(IR)
	          END DO
	        END DO
	      END IF

	      CNOW = .NOT.CNOW

	    END DO

	  END IF

	END DO

ccc	NFCALLS = NFCALLS + 1
	CALL DHCOMP (UPH0, RLNDET,CDET, DHFAKT)
	write(44,*) dreal(uph0),dimag(uph0)
	write(45,*) (dreal(dhfakt(ih)),dimag(dhfakt(ih)),ih=1,nhyd)
	dhfaktsum = (0.0d0,0.0d0)
	do ih = 1, nhyd
	  dhfaktsum = dhfaktsum + dhfakt(ih)
	end do
	write(46,*) dreal(dhfaktsum),dimag(dhfaktsum)
	close(44)
	close(45)
	close(46)

	RETURN
	END