
C	subroutine for rpress.f

C	Sven Ivansson 1990





	SUBROUTINE RPRESSINT ( FSOURC,PM0GOR, IUROT, NPUDIM, 
     $	IFILONPRE, INTREST, NPARTS,IBESSARR,NPUARR,
     $	FILMOD,XYONENORM,UPHOERN,AIMEND, UPHOERNLOK, 
     $	NHYDLI, ICWORKLIM0,ICWORKDIM,CPRESSW )
C	cpressw(1:nhydli) := complex pressure frequency-response 
C	in "(10**15 N)/(km**2)" for the respective depths and distances 
C	( "exp(-i*omeg*t)" with Re(omeg)>0.0 is behind )

	USE INMODELa
	USE INRANGEa
	USE PREDHASa
	USE INTWORKa

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       nol <= noldim , nhyd <= nhyddim , li <= lidim
	INTEGER NZERODIM
	REAL*8 TOLDIFF
	PARAMETER ( NZERODIM = 700 , TOLDIFF = 1.0D-14 )
C	nzero <= nzerodim

	INCLUDE 'INPARAM.INC'
	INCLUDE 'INPINT.INC'
	INCLUDE 'INPINTX.INC'
	INCLUDE 'INPINTASX.INC'
	INCLUDE 'INTAIL.INC'
	INCLUDE 'INTWORK.INC'
	INCLUDE 'BRANCHPS.INC'
	INCLUDE 'CIRCDATA.INC'     ! up0circ and radiecirc
	INCLUDE 'INSPEC.INC'
	INCLUDE 'DDULOK.INC'
	INCLUDE 'SOURCFW.INC'

	INTEGER IUROT, NPUDIM, INTREST, IFILONPRE, NPARTS,
     $	IBESSARR(*),NPUARR(*), NHYDLI
	REAL*8 PM0GOR(*), XYONENORM(*), AIMEND(*)
	COMPLEX*16 FSOURC, UPHOERN(NPUDIM,*)
	CHARACTER*32 FILMOD(*)

	COMPLEX*16 UPHOERNLOK(*)

	INTEGER ICWORKLIM0,ICWORKDIM
	COMPLEX*16 CPRESSW(*)

	INTEGER IH,IR, IHR, JL, JP, NPU, IQUARTP, NEU,KEU, NKVAR, 
     $	IFILONLOK, NFCAL, NK, NPLINMAX, IASBACK, 
     $	ICUTTYPOLD, NZERO, IREAD, I, ISLASK,ISLASK1, 
     $	ISHEET(nzerodim)
	REAL*8 EPSEU, EPSMAGNLOK, HSHAREKVAR, AQ,BQ, UPREAL1,UPREAL2,
     $	UPIMAG1,UPIMAG2, YBELOW, AIMSLUT, XYONEMIN,XYONEMAX, 
     $	SLASK, SLASK1,SLASK2, PI, 
     $	EPSCONT(lidim), UPZRE(nzerodim), UPZIM(nzerodim), EPSONE(1)
	COMPLEX*16 UP, CFRES0,CFRES, CSLASK, 
     $	GAMMA(lidim), UPHTAIL(2), UPHKVAR(4), ZMB(2), SSZ(2), 
     $	CPRESSWLOK(nhyddim), RFAKTONE(1)
	LOGICAL PREDONE, DIMDONE

	EXTERNAL DHCOMP,DHCOMPCIRC, RCOMP,RCOMPPOS, R1COMP,R1COMPPOS, 
     $	DHBRANCHSUM, DHBRANCH, DHBRANCHSS, RBRANCHSS,RBRANCHPOSSS, 
     $	CDETSUB, RCOMPONE,RCOMPPOSONE, TERFUN

C       -----
	INTEGER IUPVAL, IZERO1,IZERO2
	REAL*8  UPRE1,UPRE2, UPIM1,UPIM2
	LOGICAL UPVALFIRST
	SAVE IUPVAL, IZERO1,IZERO2, UPRE1,UPRE2, UPIM1,UPIM2, UPVALFIRST
C       -----

	SAVE DIMDONE
	DATA PI /3.141592653589793D0/, UPVALFIRST /.TRUE./, DIMDONE /.FALSE./

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

	IF (LI.GT.lidim) STOP '%%% rpressint, lidim'
	IF (NHYD.GT.nhyddim) STOP '%%% rpressint, nhyddim'
	IF (.NOT.DIMDONE) THEN
	  ALLOCATE ( EPS(lidim), EPSKVAR(lidim), DHFAKT(nhyddim), RFAKT(lidim) )
	  DIMDONE = .TRUE.
	END IF

	OMEG2 = OMEG*OMEG
	IF (IUROT.EQ.0) THEN
	  U0 = (1.0D0,0.0D0)
	  RK0 = OMEG
	ELSE
	  RK0 = CDABS(OMEG)     ! rk0 is now positive real
	  U0 = RK0 / OMEG       ! rk0=omeg*u0, arg(u0)=-arg(omeg)
	END IF
	IF (RHO(NOL).GT.0.0D0) THEN
	  UPBPA = CDSQRT(UA2(NOL)) / U0
	  UPBPA2 = UPBPA * UPBPA
	  IF (ISOLID(NOL).NE.0) THEN
	    UPBPB = CDSQRT(UB2(NOL)) / U0
	    UPBPB2 = UPBPB * UPBPB
	  ELSE
	    UPBPB = (0.0D0,0.0D0)
	    UPBPB2 = (0.0D0,0.0D0)
	  END IF
	ELSE
	  UPBPA = (0.0D0,0.0D0)
	  UPBPA2 = (0.0D0,0.0D0)
	  UPBPB = (0.0D0,0.0D0)
	  UPBPB2 = (0.0D0,0.0D0)
	END IF
	IF (RHO(0).GT.0.0D0) THEN
	  UPBPASURF = CDSQRT(UA2(0)) / U0
	  UPBPA2SURF = UPBPASURF * UPBPASURF
	  IF (ISOLID(0).NE.0) THEN
	    UPBPBSURF = CDSQRT(UB2(0)) / U0
	    UPBPB2SURF = UPBPBSURF * UPBPBSURF
	  ELSE
	    UPBPBSURF = (0.0D0,0.0D0)
	    UPBPB2SURF = (0.0D0,0.0D0)
	  END IF
	ELSE
	  UPBPASURF = (0.0D0,0.0D0)
	  UPBPA2SURF = (0.0D0,0.0D0)
	  UPBPBSURF = (0.0D0,0.0D0)
	  UPBPB2SURF = (0.0D0,0.0D0)
	END IF

	RW0 = -FSOURC*OMEG*RK0  ! to accomodate our "omeg*" modification of 
C	the second component of the "Aki (!) displacement-stress vector 
C	(Aki & Richards 1980, 7.4.2) but with omeg* for the first two 
C	components" and the - for pressure rather than stress and part of 
C	"dk=omeg*du" (including, when iurot=1, the "du rotation")
	IF (NSOURCFW.GT.0) THEN
	  DO I = 1, NSOURCFW
	    RW0FW(I) = AMSOURCFW(I) * 
     $	    ( -RK0 * (0.5D0/PI)*(UA2SOURCFW(I)/RHOSOURCFW(I)) * OMEG )
	    IF (UB2SOURCFW(I).NE.(0.0D0,0.0D0)) THEN
	      RS0FW(I) = AMSOURCFW(I) * 
     $	      ( RK0 * (1/PI)*(UA2SOURCFW(I)/UB2SOURCFW(I)) )
	    ELSE
	      RS0FW(I) = (0.0D0,0.0D0)
	    END IF
	  END DO
	ELSE IF (NSOURCFW.LT.-1) THEN
	  IF (NSOURCFW.NE.-2) STOP 34441   ! special from 'avh/multi'
	  RW0 = 0.5D0 * RW0
          STOP '%%% rpressint, nsourcfw option from avh/multi not impl'
	END IF

	IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	  DO IR = 1, LI
	    CSLASK = RK0*OR(IR)
	    GAMMA(IR) = DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))     ! i*k*r
	  END DO
	ELSE
	  IFILON = IFILONPRE
	END IF

	PREDONE = .FALSE.

	DO IHR = 1, NHYDLI
	  CPRESSW(IHR) = (0.0D0,0.0D0)
	END DO

	IF (IADINTXFI.EQ.66) THEN
	  IADTAILLOG = 6
	ELSE IF ( IADINTXFI.EQ.77 .OR. IADINTXFI.EQ.88 ) THEN
	  IADTAILLOG = IADINTXFI
	ELSE
	  IADTAILLOG = 0
	END IF

C	-------- tails first (first leftwards, then rightwards) -------------

	IF (ICONTTYP.EQ.1) THEN
C	  use  cpressw ( nhydli+1 : nhydli+nhydli )  for "old values"

	  IF (IADINTXFI.EQ.6) THEN
	    WRITE(*,*)
	    WRITE(*,*) ' TAIL log for iconttyp=1 ("continue")'
	  END IF
	  IF (IADINTXFI.NE.0) THEN
	    DO IR = 1, LI
	      EPS(IR) = PREEPSTERM * PM0GOR(IR)
	    END DO  ! absolute tolerance estimates for whole-space radiation
	  END IF
	  RK0OR = RK0 * DABS(OR(IRCONT))
	  DO IT = 2, 1, -1
	    IF (ITAIL(IT).NE.0) THEN
C	      integral along  dimag(rk0or*up) = constant
	      IBESSTYP = IBESSTAIL(IT)
	      IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	        IF (IBESSTYP.GT.0) THEN
	          IFILON = ISIGN(2,IFILONPRE)   !  "+-2 for J0"
	        ELSE
	          IFILON = ISIGN(1,IFILONPRE)   !  "+-1 for 0.5*H0(1)"
	        END IF
	      END IF
	      IF (DIMAG(RK0OR).NE.0.0D0) THEN
	        IRKREAL = 0
	      ELSE IF (DIMAG(UPRTAIL(IT)).NE.0.0D0) THEN
	        IRKREAL = 0
	      ELSE
	        IRKREAL = 1
	      END IF
	      IF ( IRKREAL.EQ.0 .AND. IBESSTYP.EQ.-2 ) THEN
	        IRKSPEC = 1
	      ELSE
	        IRKSPEC = 0
	      END IF
	      IF (INTREST.EQ.0) THEN
	        IF (ITAIL(IT).EQ.1) THEN
	          IDHTYP = 1
	          IF (.NOT.PREDONE) THEN
	            CALL PREDHCOMPAS
	            PREDONE = .TRUE.
	          END IF
	        ELSE
	          IDHTYP = 0
	        END IF
	      ELSE
	        IDHTYP = 2
	        IF (.NOT.PREDONE) THEN
	          CALL PREDHCOMPAS
	          PREDONE = .TRUE.
	        END IF
	      END IF
	      IF (IT.EQ.1) THEN
	        UPHTAIL(1) = UP0TAIL(1)
	        SLASK = DREAL(RK0OR*UPHTAIL(1))
	        IF (SLASK.LE.0.0D0) STOP 416
	        ISLASK = 1 + SLASK/PI
	        SLASK = (ISLASK-0.25D0)*PI - SLASK
	        IF (SLASK.LT.0.0D0) SLASK = SLASK + PI
	        CSLASK = (SLASK/RK0OR) * UPRTAIL(1)
	        UPHTAIL(2) = UPHTAIL(1) + CSLASK ! Re(rk0*uphtail(2))=m*pi-pi/4
	        TEXTLOG(9:12) = 'righ'
	      ELSE
	        UPHTAIL(2) = UP0TAIL(2)
	        SLASK = DREAL(RK0OR*UPHTAIL(2))
	        IF (SLASK.GE.0.0D0) STOP 426
	        ISLASK = -1 + SLASK/PI
	        SLASK = (ISLASK+0.75D0)*PI - SLASK
	        IF (SLASK.GT.0.0D0) SLASK = SLASK - PI
	        CSLASK = (SLASK/RK0OR) * UPRTAIL(2)
	        UPHTAIL(1) = UPHTAIL(2) + CSLASK ! Re(rk0*uphtail(1))=m*pi-pi/4
	        TEXTLOG(9:12) = 'left'
	      END IF
	      HMAXAS = CDABS(CSLASK) / ANPHPMINTERM
	      IF (IADINTXFI.EQ.0) THEN
	        CALL TRAPEZFI (NHYD,LI, DHCOMP,RCOMP, IFILON,
     $	        RCOMPPOS,GAMMA, 2,UPHTAIL, HMAXAS, NHYDLI,CPRESSW, 
     $	        NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	      ELSE
	        CALL ADINTXFI (NHYD,LI, DHCOMP,RCOMP, IFILON,
     $	        RCOMPPOS,GAMMA, 2,UPHTAIL, HMAXAS,SMAXPRE,0.0D0,EPS,
     $	        NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, ARGDIFLIM1,
     $	        ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, IADTAILLOG,
     $	        TEXTLOG, NHYDLI,CPRESSW, SLASK, NFCAL, DHFAKT,RFAKT, 
     $	        ICWORKLIM0,ICWORKDIM,CPRESSW)
	      END IF
	      IF (IT.EQ.1) THEN
	        NFCALLS(5) = NFCALLS(5) + NFCAL
	      ELSE
	        NFCALLS(6) = NFCALLS(6) + NFCAL
	      END IF
	      UPPID = (PI/RK0OR) * UPRTAIL(IT)
	      UPPIDH = 0.5D0 * UPPID
	      HMAXAS = CDABS(UPPID) / ANPHPMINTERM
	      IQUARTP = 0
	      DO IR = 1, LI
	        EPSCONT(IR) = PREEPSCONT * PM0GOR(IR)
	      END DO
	      DO WHILE (.TRUE.)
	        IQUARTP = IQUARTP + 1
	        IF (IT.EQ.1) THEN
	          UPHTAIL(1) = UPHTAIL(2)
	          UPHTAIL(2) = UPHTAIL(1) + UPPIDH
	        ELSE
	          UPHTAIL(2) = UPHTAIL(1)
	          UPHTAIL(1) = UPHTAIL(2) - UPPIDH
	        END IF
	        DO IR = 1, LI
	          IHR = IR
	          DO IH = 1, NHYD
	            CPRESSW(NHYDLI+IHR) = CPRESSW(IHR)
	            IHR = IHR + LI
	          END DO
	        END DO
	        IF (IADINTXFI.EQ.0) THEN
	          CALL TRAPEZFI (NHYD,LI, DHCOMP,RCOMP, IFILON,
     $	          RCOMPPOS,GAMMA, 2,UPHTAIL, HMAXAS, NHYDLI,CPRESSW, 
     $	          NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	        ELSE
	          CALL ADINTXFI (NHYD,LI, DHCOMP,RCOMP, IFILON,
     $	          RCOMPPOS,GAMMA, 2,UPHTAIL, HMAXAS,SMAXPRE,0.0D0,EPS,
     $	          NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, ARGDIFLIM1,
     $	          ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, IADTAILLOG,
     $	          TEXTLOG, NHYDLI,CPRESSW, SLASK, NFCAL, DHFAKT,RFAKT, 
     $	          ICWORKLIM0,ICWORKDIM,CPRESSW)
	        END IF
	        IF (IT.EQ.1) THEN
	          NFCALLS(5) = NFCALLS(5) + NFCAL
	        ELSE
	          NFCALLS(6) = NFCALLS(6) + NFCAL
	        END IF
	        DO IR = 1, LI
	          SLASK = EPSCONT(IR)
	          IHR = IR
	          DO IH = 1, NHYD
	            IF (CDABS(CPRESSW(IHR)-CPRESSW(NHYDLI+IHR)).GT.
     $	            SLASK) GOTO 10
	            IHR = IHR + LI
	          END DO
	        END DO
	        GOTO 11
 10             CONTINUE
	      END DO
 11	      IF (IADINTXFI.EQ.6) THEN
	        WRITE(*,*) '  it,iquartp ',IT,IQUARTP
	        IF (IT.EQ.1) THEN
	          WRITE(*,*) '      to ',UPHTAIL(2)
	        ELSE
	          WRITE(*,*) '    from ',UPHTAIL(1)
	        END IF
	      END IF
	    END IF
	  END DO  ! end it loop

	ELSE IF (ICONTTYP.EQ.2) THEN

	  IF (IADINTXFI.EQ.6) THEN
	    WRITE(*,*)
	    WRITE(*,*) ' TAIL log for iconttyp=2 ("Euler tails")'
	  END IF
	  NHYD2 = 2*NHYD
	  DO IR = 1, LI
C	    use  cpressw ( nhydli+1 : nhydli+2*nhyd )  for "cos & sin" sums
	    IF (OR(IR).NE.0.0D0) THEN
	      RK0OR = RK0 * DABS(OR(IR))
	    ELSE
	      STOP 6003
	    END IF
	    IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	      CSLASK = RK0*OR(IR)
	      CSLASK = DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))     ! i*k*r
	      GAMMAS(1) = CSLASK
	      GAMMAS(2) = CSLASK
	    END IF
	    DO IT = 2, 1, -1
	      IF (ITAIL(IT).NE.0) THEN
C	        integral along  dimag(rk0or*up) = constant
	        IBESSTYP = IBESSTAIL(IT)
	        IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	          IF (IBESSTYP.GT.0) THEN
	            IFILON = ISIGN(2,IFILONPRE)   !  "+-2 for J0"
	          ELSE
	            IFILON = ISIGN(1,IFILONPRE)   !  "+-1 for 0.5*H0(1)"
	          END IF
	        END IF
	        IF (DIMAG(RK0OR).NE.0.0D0) THEN
	          IRKREAL = 0
	        ELSE IF (DIMAG(UPRTAIL(IT)).NE.0.0D0) THEN
	          IRKREAL = 0
	        ELSE
	          IRKREAL = 1
	        END IF
	        IF ( IRKREAL.EQ.0 .AND. IBESSTYP.EQ.-2 ) THEN
	          IRKSPEC = 1
	        ELSE
	          IRKSPEC = 0
	        END IF
	        IF (INTREST.EQ.0) THEN
	          IF (ITAIL(IT).EQ.1) THEN
	            IDHTYP = 1
	            IF (.NOT.PREDONE) THEN
	              CALL PREDHCOMPAS
	              PREDONE = .TRUE.
	            END IF
	          ELSE
	            IDHTYP = 0
	          END IF
	        ELSE
	          IDHTYP = 2
	          IF (.NOT.PREDONE) THEN
	            CALL PREDHCOMPAS
	            PREDONE = .TRUE.
	          END IF
	        END IF
	        IF (IT.EQ.1) THEN
	          UPHTAIL(1) = UP0TAIL(1)
	          SLASK = DREAL(RK0OR*UPHTAIL(1))
	          IF (SLASK.LE.0.0D0) STOP 416
	          ISLASK = 1 + SLASK/PI
	          SLASK = (ISLASK-0.25D0)*PI - SLASK
	          IF (SLASK.LT.0.0D0) SLASK = SLASK + PI
	          CSLASK = (SLASK/RK0OR) * UPRTAIL(1)
	          UPHTAIL(2) = UPHTAIL(1) +CSLASK !Re(rk0*uphtail(2))=m*pi-pi/4
	          TEXTLOG(9:12) = 'righ'
	        ELSE
	          UPHTAIL(2) = UP0TAIL(2)
	          SLASK = DREAL(RK0OR*UPHTAIL(2))
	          IF (SLASK.GE.0.0D0) STOP 426
	          ISLASK = -1 + SLASK/PI
	          SLASK = (ISLASK+0.75D0)*PI - SLASK
	          IF (SLASK.GT.0.0D0) SLASK = SLASK - PI
	          CSLASK = (SLASK/RK0OR) * UPRTAIL(2)
	          UPHTAIL(1) = UPHTAIL(2) +CSLASK !Re(rk0*uphtail(1))=m*pi-pi/4
	          TEXTLOG(9:12) = 'left'
	        END IF
	        HMAXAS = CDABS(CSLASK) / ANPHPMINTERM
	        DO IH = 1, NHYD2
	          CPRESSW(NHYDLI+IH) = (0.0D0,0.0D0)
	        END DO
	        IF (IADINTXFI.EQ.0) THEN
	          CALL TRAPEZFI (NHYD,2, DHCOMP,R1COMP, IFILON,
     $	          R1COMPPOS,GAMMAS, 2,UPHTAIL, HMAXAS, NHYD2,
     $	          CPRESSW(NHYDLI+1), NFCAL, DHFAKT,RFAKT, 
     $	          ICWORKLIM0,ICWORKDIM,CPRESSW)
	        ELSE
	          EPSAS(1) = PREEPSTERM * PM0GOR(IR)
	          EPSAS(2) = EPSAS(1)
	          CALL ADINTXFI (NHYD,2, DHCOMP,R1COMP, IFILON,
     $	          R1COMPPOS,GAMMAS, 2,UPHTAIL, HMAXAS,SMAXPRE,0.0D0,
     $	          EPSAS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, 
     $	          ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, 
     $	          IADTAILLOG,TEXTLOG, NHYD2,CPRESSW(NHYDLI+1), SLASK, 
     $	          NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	        END IF
	        IF (IT.EQ.1) THEN
	          NFCALLS(5) = NFCALLS(5) + NFCAL
	        ELSE
	          NFCALLS(6) = NFCALLS(6) + NFCAL
	        END IF
	        IHR = IR
	        DO IH = 1, NHYD
	          ISLASK = 2*IH
	          CPRESSW(IHR) = CPRESSW(IHR) + CPRESSW(NHYDLI+
     $	          ISLASK-1) + CPRESSW(NHYDLI+ISLASK)   ! add cos and sin terms
	          IHR = IHR + LI
	        END DO
	        UPPID = (PI/RK0OR) * UPRTAIL(IT)
	        UPPIDH = 0.5D0 * UPPID
	        HMAXAS = CDABS(UPPID) / ANPHPMINTERM
	        IF (IT.EQ.1) THEN
	          UPH0AS = UPHTAIL(2)
	          UPH1AS = UPH0AS + UPPIDH
	        ELSE
	          UPH0AS = UPHTAIL(1)
	          UPH1AS = UPH0AS - UPPIDH
	        END IF
	        CALL TERFUN ( -1, CPRESSW(NHYDLI+1), 
     $	        ICWORKLIM0,ICWORKDIM,CPRESSW  )
	        IHR = IR
	        DO IH = 1, NHYD
	          CPRESSW(IHR) = CPRESSW(IHR) + CPRESSW(NHYDLI+2*IH) ! add sin
	          IHR = IHR + LI
	        END DO
	        EPSEU = PREEPSEU * PM0GOR(IR)
	        CALL EULERACC ( NHYD2,TERFUN,(-1.0D0,0.0D0), 0,-1, 
     $	        KMAXEU,EPSEU, NEU,KEU,CPRESSW(NHYDLI+1), 
     $	        ICWORKLIM0,ICWORKDIM,CPRESSW  )
	        IF (IADINTXFI.EQ.6) THEN
	          WRITE(*,*) '  ir,it,neu,keu ',IR,IT,NEU,KEU
	          WRITE(*,*) '          uppid ',UPPID
	        END IF
	        IHR = IR
	        DO IH = 1, NHYD
	          ISLASK = NHYDLI + 2*IH
	          CPRESSW(IHR) = CPRESSW(IHR) + CPRESSW(ISLASK-1) + 
     $	          CPRESSW(ISLASK)         ! add cos and sin terms
	          IHR = IHR + LI
	        END DO
	      END IF
	    END DO  ! end it loop
	  END DO  ! end ir loop

	END IF

C	-------- main "field parts" next (in the order "jp=nparts,1,-1") -----

	IF (IADINTXFI.NE.0) THEN
	  DO IR = 1, LI
	    EPS(IR) = PREEPS * PM0GOR(IR)
	  END DO  ! absolute tolerance estimates for whole-space radiation
	  SUPEPSMAGN = 0.0D0
	END IF

	DO JP = NPARTS, 1, -1
	  IBESSTYP = IBESSARR(JP)
	  IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	    IF (IBESSTYP.GT.0) THEN
	      IFILON = ISIGN(2,IFILONPRE)   !  "+-2 for J0"
	    ELSE
	      IFILON = ISIGN(1,IFILONPRE)   !  "+-1 for 0.5*H0(1)"
	    END IF
	  END IF
	  NPU = NPUARR(JP)
	  IF (NPU.GT.0) THEN
C	    path(s) now
	    IF (IADINTXFI.NE.0) TEXTLOG(9:12) = 'main'
	    ISLASK = 0
	    DO JL = 1, NPU
	      IF (DREAL(UPHOERN(JL,JP)).EQ.+991.0D0) THEN
	        UPHOERNLOK(JL) = DCMPLX(DREAL(UPBPA),DIMAG(UPHOERN(JL,JP)))
	        ISLASK = 1
	      ELSE IF (DREAL(UPHOERN(JL,JP)).EQ.-991.0D0) THEN
	        UPHOERNLOK(JL) = DCMPLX(-DREAL(UPBPA),DIMAG(UPHOERN(JL,JP)))
	        ISLASK = 1
	      ELSE IF (DREAL(UPHOERN(JL,JP)).EQ.+992.0D0) THEN
	        UPHOERNLOK(JL) = DCMPLX(DREAL(UPBPB),DIMAG(UPHOERN(JL,JP)))
	        ISLASK = 1
	      ELSE IF (DREAL(UPHOERN(JL,JP)).EQ.-992.0D0) THEN
	        UPHOERNLOK(JL) = DCMPLX(-DREAL(UPBPB),DIMAG(UPHOERN(JL,JP)))
	        ISLASK = 1
	      ELSE
	        UPHOERNLOK(JL) = UPHOERN(JL,JP)
	      END IF
	    END DO
	    IF (ISLASK.NE.0) THEN
	      UPREF = 0.5D0 * ( UPHOERNLOK(1) + UPHOERNLOK(2) )
	      IF ( ICUTTYP.NE.2 .AND. ICUTTYP.NE.3 ) STOP 3471
	    END IF
	    IF (INTREST.EQ.0) THEN
	     IDHTYP = 0
	    ELSE
	      IDHTYP = 2
	      CALL PREDHCOMPAS
	      PREDONE = .TRUE.
	    END IF
	    IF (DIMAG(RK0).NE.0.0D0) THEN
	      IRKREAL = 0
	    ELSE
	      DO JL = 1, NPU
	        IF (DIMAG(UPHOERNLOK(JL)).NE.0.0D0) THEN
	          IRKREAL = 0
	          GOTO 210
	        END IF
	      END DO
	      IRKREAL = 1
	    END IF
 210	    IF ( IRKREAL.EQ.0 .AND. IBESSTYP.EQ.-2 ) THEN
	      IRKSPEC = 1
	    ELSE
	      IRKSPEC = 0
	    END IF
	    IF (IADINTXFI.EQ.0) THEN
	      CALL TRAPEZFI (NHYD,LI, DHCOMP,RCOMP, IFILON,RCOMPPOS,
     $	      GAMMA, NPU,UPHOERNLOK, HMAX, NHYDLI,CPRESSW, 
     $	      NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	      NFCALLS(2) = NFCALLS(2) + NFCAL
	    ELSE
	      IF (HKVAR.GT.0.0D0) THEN
	        OPEN(LUNKVAR,FILE='FOR111.DAT') !!! possible change: "scratch" 
	        WRITE(LUNKVAR,*) TEXTLOG,'   "KVAR" FILE'
	        WRITE(LUNKVAR,*) ' freq,jp ',SNGL(DREAL(OMEG)/(2*PI)),JP
	      END IF
	      CALL ADINTXFI (NHYD,LI, DHCOMP,RCOMP, IFILON,RCOMPPOS,
     $	      GAMMA, NPU,UPHOERNLOK, HMAX,SMAXPRE,EPSH01Q,EPS,
     $	      NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, ARGDIFLIM1,
     $	      ARGDIFLIM0,ARGDIFLIM2, HKVAR,LUNKVAR,NKVAR, IADINTXFI,
     $	      TEXTLOG, NHYDLI,CPRESSW, EPSMAGNLOK, NFCAL, 
     $	      DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	      IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	      NFCALLS(2) = NFCALLS(2) + NFCAL
	      IF (HKVAR.GT.0.0D0) THEN
	        REWIND LUNKVAR
	        READ(LUNKVAR,*)
	        READ(LUNKVAR,*)
	        IF (NKVAR.GT.0) THEN
	          CSLASK = OMEG * U0 * OR(LI)
	          AQ = DREAL(CSLASK)
	          BQ = DIMAG(CSLASK)
	          DO NK = 1, NKVAR
	            READ(LUNKVAR,*) UPREAL1,UPREAL2,HSHAREKVAR  !"1&2" positive
	            UPIMAG1 = -DABS(FAKTKVAR) * DMIN1 ( 
     $	            DABS(UPREAL2-UPREAL1) , 100*HKVAR )
	            SLASK = (-dlimexp-BQ*UPREAL1) / AQ
	            IF (SLASK.GT.UPIMAG1) UPIMAG1 = SLASK
	            SLASK = (-dlimexp-BQ*UPREAL2) / AQ
	            IF (SLASK.GT.UPIMAG1) UPIMAG1 = SLASK
	            UPIMAG2 = 0.0D0   ! " = dimag(uphkvar0) = dimag(uphkvar1) "
	            IF (FAKTKVAR.GT.0.0D0) THEN
	              IDET = 1                 ! observe this important setting
	              CALL ZERODIVY (CDETSUB, UPREAL1,UPREAL2,UPIMAG1,
     $	              UPIMAG2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, 
     $	              DERIVDRELMAX,DRELMAX,DRELGOALFAKT, ARGDIFSTORE, 
     $	              NPLINMAX,NFCALLS(4),IASBACK, 
     $	              YBELOW)
	              IDET = 0                 ! observe this resetting
	            ELSE
	              YBELOW = UPIMAG1
	            END IF
	            YBELOW = 0.6D0 * YBELOW  !get away from possible pole below
	            UPHKVAR(1) = DCMPLX ( UPREAL1 , 0.0D0 )
	            UPHKVAR(2) = DCMPLX ( UPREAL1 , YBELOW )
	            UPHKVAR(3) = DCMPLX ( UPREAL2 , YBELOW )
	            UPHKVAR(4) = DCMPLX ( UPREAL2 , 0.0D0 )
	            DO IR = 1, LI
	              EPSKVAR(IR) = HSHAREKVAR * EPS(IR)
	            END DO
	            IRKREAL = 0                ! observe this important setting
	            IRKSPEC = 0
	            CALL ADINTXFI (NHYD,LI, DHCOMP,RCOMP, 
     $	            IFILON,RCOMPPOS,GAMMA, 4,UPHKVAR, HMAX,SMAXPRE,
     $	            0.0D0,EPSKVAR,NELMAXPRE,NBACKMAX, IXTYP,NRUS,
     $	            IRUS,JRUS, ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 
     $	            0.0D0,0,ISLASK, IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, 
     $	            SLASK, NFCAL, DHFAKT,RFAKT, 
     $	            ICWORKLIM0,ICWORKDIM,CPRESSW)
	            NFCALLS(2) = NFCALLS(2) + NFCAL
	          END DO
	        END IF
	        CLOSE(LUNKVAR)
	      END IF
	    END IF
	    UPREF = (0.0D0,0.0D0)
	  ELSE IF (NPU.EQ.0) THEN
C	    mode(s) now (in normal case here, "ibesstyp.lt.0")
	    IF (IADINTXFI.NE.0) TEXTLOG(9:12) = 'mode'
	    IDHTYP = 0
C           -----
	    IF (UPVALFIRST) THEN
	      OPEN(UNIT=8,FILE='UPVAL.DAT',STATUS='OLD')
	      READ(8,*) IUPVAL
	      IF (IUPVAL.GT.0) THEN
	        READ(8,*) UPRE1,UPRE2, UPIM1,UPIM2
	        WRITE(*,*) '*********************************************'
	        WRITE(*,*) 'UPVAL!!!'
	        WRITE(*,*) UPRE1,UPRE2,' ! upre1,upre2'
	        WRITE(*,*) UPIM1,UPIM2,' ! upim1,upim2'
	        WRITE(*,*) '*********************************************'
	      ELSE IF (IUPVAL.LT.0) THEN
	        READ(8,*) IZERO1,IZERO2
	        WRITE(*,*) '*********************************************'
	        WRITE(*,*) 'UPVAL!!!'
	        WRITE(*,*) IZERO1,IZERO2,' ! izero1,izero2'
	        WRITE(*,*) '*********************************************'
	      END IF
	      CLOSE(8)
	      UPVALFIRST = .FALSE.
	    END IF
C           -----
	    OPEN(UNIT=8,FILE=FILMOD(JP),STATUS='OLD')
	    IF (XYONENORM(JP).NE.0.0D0) THEN
	      READ(8,*) ISLASK, SLASK
	    ELSE
	      READ(8,*) ISLASK, SLASK,I
	      IF (I.NE.KFWTYP) STOP 64060
	    END IF
	    IF ( DABS((SLASK*2*PI/DREAL(OMEG))-1) .GT. 0.00001D0 ) STOP 41155
	    DO I = 1, ISLASK
	      READ(8,*)
	    END DO
	    READ(8,*) IREAD, ISLASK,ISLASK1,SLASK,I
	    IF ( ISLASK.NE.ICUTTYP .OR. ISLASK1.NE.ISHEETTYP .OR. 
     $	    I.NE.IUROT ) STOP 34155
	    IF ( DIMAG(OMEG) .EQ. 0.0D0 ) THEN
	      IF (SLASK.NE.0.0D0) STOP '%%% dimag(omeg) inconsistency'
	    ELSE
	      IF ( DABS((SLASK/DIMAG(OMEG))-1) .GT. 0.00001D0 ) THEN
	        WRITE(*,*) OMEG,' ! omeg'
	        WRITE(*,*) SLASK,' ! ala from filmod(jp)'
	        STOP '%%% dimag(omeg) inconsistency'
	      END IF
	    END IF
	    IF (XYONENORM(JP).GT.0.0D0) THEN
C	      integration around the pole (positive orientation)
C             -----
	       IF (IUPVAL.NE.0) STOP '%%% rpressint, upval not impl'
C             -----
	      ICUTTYPOLD = ICUTTYP
	      ICUTTYP = 2
	      IRKREAL = 0
	      IRKSPEC = 0	! "no spec" in "rcompone,rcompposone"
	      IF (IREAD.LT.2) STOP 41151
	      READ(8,*) NZERO
	      IF (NZERO.GT.nzerodim) STOP 640
	      DO I = 1, NZERO
	        READ(8,*) UPZRE(I),UPZIM(I),ISHEET(I)
	        IF (ISHEET(I).EQ.0) STOP 5761
	      END DO
	      READ(8,*) ISLASK
	      IF (ISLASK.NE.NZERO) STOP 370
	      DO I = 1, NZERO
	        READ(8,*) XYONEMIN,XYONEMAX
	        ISHEETTYP = ISHEET(I)
	        UP = DCMPLX(UPZRE(I),UPZIM(I))
	        CALL RCOMP (UP,RFAKT)
	        SLASK1 = DMIN1 ( 2*XYONEMIN, 0.8D0*XYONEMIN+0.2D0*XYONEMAX )
	        SLASK2 = 0.2D0*XYONEMIN+0.8D0*XYONEMAX
	        IF ( XYONENORM(JP) .LE. SLASK1 ) THEN
	          RADIECIRC = SLASK1
	        ELSE IF ( XYONENORM(JP) .LE. SLASK2 ) THEN
	          RADIECIRC = XYONENORM(JP)
	        ELSE
	          RADIECIRC = SLASK2
	        END IF
	        UP0CIRC = DCMPLX ( UPZRE(I), UPZIM(I) )
	        ZMB(1) = (0.0D0,0.0D0)
	        ZMB(2) = DCMPLX(2*PI,0.0D0)
	        DO IH = 1, NHYD
	          CPRESSWLOK(IH) = (0.0D0,0.0D0)
	        END DO
	        IF (IADINTXFI.EQ.0) THEN
	          CALL TRAPEZFI (NHYD,1, DHCOMPCIRC,RCOMPONE, 
     $	          0,RCOMPPOSONE,GAMMA, 2,ZMB, HMAX, NHYD,CPRESSWLOK, 
     $	          NFCAL, DHFAKT,RFAKTONE, ICWORKLIM0,ICWORKDIM,CPRESSW)
	        ELSE
	          EPSONE(1) = -1.0D0
	          DO IR = 1, LI
	            SLASK = 
     $              DMAX1 ( DABS(DREAL(RFAKT(IR))), DABS(DIMAG(RFAKT(IR))) )
	            IF (SLASK.GT.0.0D0) THEN
	              IF (EPSONE(1).LT.0.0D0) THEN
	                EPSONE(1) = PM0GOR(IR)/SLASK
	              ELSE
	                EPSONE(1) = DMIN1 ( PM0GOR(IR)/SLASK, EPSONE(1) )
	              END IF
	            END IF
	          END DO
	          IF (EPSONE(1).LT.0.0D0) THEN
	            DO IH = 1, NHYD
	              DHFAKT(IH) = (0.0D0,0.0D0) ! arbitrary ("all rfakt" zero)
	            END DO
	            NFCAL = 0
	          ELSE
	            EPSONE(1) = (PREEPS/NZERO) * EPSONE(1)
	            CALL ADINTXFI (NHYD,1, DHCOMPCIRC,RCOMPONE, 
     $	            0,RCOMPPOSONE,GAMMA, 2,ZMB, HMAX,SMAXPRE,
     $	            EPSH01Q,EPSONE,NELMAXPRE,NBACKMAX, IXTYP,NRUS,
     $	            IRUS,JRUS, ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 
     $	            0.0D0,0,NKVAR, IADINTXFI,TEXTLOG, NHYD,CPRESSWLOK, 
     $	            EPSMAGNLOK, NFCAL, DHFAKT,RFAKTONE, 
     $	            ICWORKLIM0,ICWORKDIM,CPRESSW)
	            IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	          END IF
	        END IF
	        NFCALLS(1) = NFCALLS(1) + NFCAL
	        IHR = 0
	        DO IH = 1, NHYD
	          CSLASK = CPRESSWLOK(IH)
	          DO IR = 1, LI
	            IHR = IHR + 1
	            CPRESSW(IHR) = CPRESSW(IHR) + CSLASK*RFAKT(IR)
	          END DO
	        END DO
	      END DO
	      CLOSE(8)
	      ICUTTYP = ICUTTYPOLD
	      ISHEETTYP = 1
	    ELSE IF (XYONENORM(JP).EQ.0.0D0) THEN
C	      residue contributions (positive orientation) from "rmodpfw"
	      IRKSPEC = 0   ! "no spec" in "rmodpfw"
	      CFRES0 = -OMEG2*RK0                            !"2*pi*" here
	      CFRES0 = DCMPLX(-DIMAG(CFRES0),DREAL(CFRES0))  !"i*" rest
	      NZERO = IREAD
	      IF (NZERO.GT.nzerodim) STOP 64067
	      DO I = 1, NZERO
	        READ(8,*) ISLASK,UP
	        IF (ISLASK.NE.NHYD) STOP 1671
C               -----
	        IF (IUPVAL.GT.0) THEN
	          IF ( DREAL(UP).LT.UPRE1 .OR. DREAL(UP).GT.UPRE2 .OR. 
     &	           DIMAG(UP).LT.UPIM1.OR.DIMAG(UP).GT.UPIM2 ) THEN
	            DO IH = 1, NHYD
	              READ(8,*) SLASK, CSLASK
	            END DO
	            CYCLE
	          END IF
	        ELSE IF (IUPVAL.LT.0) THEN
	          IF ( I.LT.IZERO1 .OR. I.GT.IZERO2 ) THEN
	            DO IH = 1, NHYD
	              READ(8,*) SLASK, CSLASK
	            END DO
	            CYCLE
	          END IF
	        END IF
C               -----
	        IF (DIMAG(UP).EQ.0.0D0) THEN
	          IRKREAL = 1
	        ELSE
	          IRKREAL = 0
	        END IF
	        CALL RCOMP (UP,RFAKT)
	        CFRES = CFRES0 * UP
	        IHR = 0
	        DO IH = 1, NHYD
	          READ(8,*) SLASK, CSLASK
	          IF ( DABS(SLASK-ZD(IDHYD(IH))) .GT. TOLDIFF ) STOP 4316
	          CSLASK = CFRES * CSLASK
	          DO IR = 1, LI
	            IHR = IHR + 1
	            CPRESSW(IHR) = CPRESSW(IHR) + CSLASK*RFAKT(IR)
	          END DO
	        END DO
	      END DO
	      CLOSE(8)
	    ELSE
C	      residue contributions (positive orientation) from "dhcomp"
C             -----
	      IF (IUPVAL.NE.0) STOP '%%% rpressint, upval not impl'
C             -----
	      IDET = -1                        ! "residue setting" to liqsi
	      ICUTTYPOLD = ICUTTYP
	      ICUTTYP = 2
	      IF (IREAD.LT.3) STOP 41151
	      READ(8,*) NZERO
	      IF (NZERO.GT.nzerodim) STOP 640
	      DO I = 1, NZERO
	        READ(8,*) UPZRE(I),UPZIM(I),ISHEET(I)
	        IF (ISHEET(I).EQ.0) STOP 5761
	      END DO
	      READ(8,*) ISLASK
	      IF (ISLASK.NE.NZERO) STOP 370
	      DO I = 1, NZERO
	        READ(8,*) XYONEMIN,XYONEMAX
	      END DO
	      READ(8,*) ISLASK
	      IF (ISLASK.NE.NZERO) STOP 370
	      DO I = 1, NZERO
	        READ(8,*) DLNDULOK,DERIVDULOK
	        ISHEETTYP = ISHEET(I)
	        UP = DCMPLX(UPZRE(I),UPZIM(I))
	        IF (DIMAG(UP).EQ.0.0D0) THEN
	          IRKREAL = 1
	        ELSE
	          IRKREAL = 0
	        END IF
	        IF ( IRKREAL.EQ.0 .AND. IBESSTYP.EQ.-2 ) THEN
	          IRKSPEC = 1
	        ELSE
	          IRKSPEC = 0
	        END IF
	        CALL DHCOMP (UP, SLASK,CSLASK, DHFAKT)
	        CALL RCOMP (UP,RFAKT)
	        IHR = 0
	        DO IH = 1, NHYD
	          CSLASK = DHFAKT(IH)
	          DO IR = 1, LI
	            IHR = IHR + 1
	            CPRESSW(IHR) = CPRESSW(IHR) + CSLASK*RFAKT(IR)
	          END DO
	        END DO
	      END DO
	      CLOSE(8)
	      NFCALLS(1) = NFCALLS(1) + NZERO
	      ICUTTYP = ICUTTYPOLD
	      ISHEETTYP = 1
	      IDET = 0                    ! resetting
	    END IF
	  ELSE
C	    branch(es) now (in normal case here, "ibesstyp.lt.0")
	    IF (IADINTXFI.NE.0) TEXTLOG(9:12) = 'brch'
	    IF ( NPU.NE.-1 .AND. NPU.NE.-2 .AND. NPU.NE.-12 ) STOP 5671
	    IF (UPBPA.EQ.(0.0D0,0.0D0)) STOP 4600
	    IF ( UPBPB.EQ.(0.0D0,0.0D0) .AND. NPU.NE.-1 ) STOP 4601
	    IDHTYP = 0
	    IRKREAL = 0
	    IF ( IRKREAL.EQ.0 .AND. IBESSTYP.EQ.-2 ) THEN
	      IRKSPEC = 1
	    ELSE
	      IRKSPEC = 0
	    END IF
	    IF (ICUTTYP.EQ.1) THEN
	      CALL RPBRANCHINT ( FSOURC,PM0GOR, IUROT, IFILONPRE, 
     $	      NPU,0.0D0,AIMEND(JP), NHYDLI, 
     $	      ICWORKLIM0,ICWORKDIM,CPRESSW )
	    ELSE
	      IF (RHO(NOL).GT.0.0D0) THEN
	        IF (DREAL(UPBPA).LE.0.0D0) STOP 6704   ! case not impl
	        IF (ISOLID(NOL).NE.0) THEN
	          IF (DREAL(UPBPB).LE.0.0D0) STOP 6705   ! case not impl
	          IF (DREAL(UPBPA).GE.DREAL(UPBPB)) STOP 540  ! case not impl
	        END IF
	      ELSE
                STOP 6703
	      END IF
	      IF (ICUTTYP.EQ.2) THEN
	        AIMSLUT = AIMEND(JP)
	      ELSE
	        IF (ICUTTYP.NE.3) STOP 7117
	        IF (AIMEND(JP).GE.AIMCUT) THEN
	          IF (AIMEND(JP).GT.AIMCUT) THEN
	            CALL RPBRANCHINT ( FSOURC,PM0GOR, IUROT, IFILONPRE, 
     $	            NPU,AIMCUT,AIMEND(JP), NHYDLI, 
     $	            ICWORKLIM0,ICWORKDIM,CPRESSW )
	          ELSE
C	            otherwise done in rpbranchint
	            IF (RHO(NOL).GT.0.0D0) THEN
	              PSPROD(1) = DREAL(UPBPA) * DIMAG(UPBPA)
	              IF (PSPROD(1).LT.0.0D0) STOP 560  ! case not implemented
	              IF (ISOLID(NOL).NE.0) THEN
	                PSPROD(2) = DREAL(UPBPB) * DIMAG(UPBPB)
	                IF (PSPROD(2).LE.PSPROD(1)) STOP 560 ! case not impl
	              ELSE
	                IF (NPU.NE.-1) STOP 560
	              END IF
	            ELSE
	              STOP 560
	            END IF
	          END IF
	          IF (NPU.EQ.-1) THEN
	            IPS = 1
	            IF (AIMCUT.GT.DIMAG(UPBPA)) THEN
	              ZMB(1) = DCMPLX(DREAL(UPBPA),AIMCUT)
	              ZMB(2) = DCMPLX(PSPROD(1)/AIMCUT,AIMCUT)
	              IF (IADINTXFI.EQ.0) THEN
	                CALL TRAPEZFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX, NHYDLI,CPRESSW, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	              ELSE
	                CALL ADINTXFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX,SMAXPRE,EPSH01Q,
     $	                EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, 
     $	                ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, 
     $	                IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, EPSMAGNLOK, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	                IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	              END IF
	              NFCALLS(3) = NFCALLS(3) + NFCAL
	            END IF
	          ELSE IF (NPU.EQ.-2) THEN
	            IPS = 2
	            IF (AIMCUT.GT.DIMAG(UPBPB)) THEN
	              ZMB(1) = DCMPLX(DREAL(UPBPB),AIMCUT)
	              ZMB(2) = DCMPLX(PSPROD(2)/AIMCUT,AIMCUT)
	              IF (IADINTXFI.EQ.0) THEN
	                CALL TRAPEZFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX, NHYDLI,CPRESSW, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	              ELSE
	                CALL ADINTXFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX,SMAXPRE,EPSH01Q,
     $	                EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, 
     $	                ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, 
     $	                IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, EPSMAGNLOK, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	                IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	              END IF
	              NFCALLS(3) = NFCALLS(3) + NFCAL
	            END IF
	          ELSE
	            IF (NPU.NE.-12) STOP 7105
	            IF (AIMCUT.GT.DIMAG(UPBPA)) THEN
	              IPS = 1
	              SLASK1 = DMIN1(DREAL(UPBPA),PSPROD(2)/AIMCUT)
	              ZMB(1) = DCMPLX(SLASK1,AIMCUT)
	              ZMB(2) = DCMPLX(PSPROD(1)/AIMCUT,AIMCUT)
	              IF (IADINTXFI.EQ.0) THEN
	                CALL TRAPEZFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX, NHYDLI,CPRESSW, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	              ELSE
	                CALL ADINTXFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX,SMAXPRE,EPSH01Q,
     $	                EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, 
     $	                ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, 
     $	                IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, EPSMAGNLOK, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	                IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	              END IF
	              NFCALLS(3) = NFCALLS(3) + NFCAL
	            END IF
	            IF (AIMCUT.GT.DIMAG(UPBPB)) THEN
	              IPS = 2
	              SLASK2 = DMAX1(PSPROD(2)/AIMCUT,DREAL(UPBPA))
	              ZMB(1) = DCMPLX(DREAL(UPBPB),AIMCUT)
	              ZMB(2) = DCMPLX(SLASK2,AIMCUT)
	              IF (IADINTXFI.EQ.0) THEN
	                CALL TRAPEZFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX, NHYDLI,CPRESSW, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	              ELSE
	                CALL ADINTXFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	                RCOMPPOS,GAMMA, 2,ZMB, HMAX,SMAXPRE,EPSH01Q,
     $	                EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, 
     $	                ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, 
     $	                IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, EPSMAGNLOK, 
     $	                NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,
     $	                CPRESSW)
	                IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	              END IF
	              NFCALLS(3) = NFCALLS(3) + NFCAL
	            END IF
	            IF ( AIMCUT.GT.DIMAG(UPBPA) .AND. 
     $	            AIMCUT.GT.DIMAG(UPBPB) .AND. 
     $	            SLASK1.LT.DREAL(UPBPA) ) THEN
	              ZMB(1) = DCMPLX(SLASK1,AIMCUT)
	              ZMB(2) = DCMPLX(DREAL(UPBPA),AIMCUT)
	              IF (IADINTXFI.EQ.0) THEN
	                CALL TRAPEZFI (NHYD,LI, DHBRANCHSUM,RCOMP, 
     $	                IFILON,RCOMPPOS,GAMMA, 2,ZMB, HMAX, NHYDLI,
     $	                CPRESSW, NFCAL, DHFAKT,RFAKT, ICWORKLIM0,
     $	                ICWORKDIM,CPRESSW)
	              ELSE
	                CALL ADINTXFI (NHYD,LI, DHBRANCHSUM,RCOMP, 
     $	                IFILON,RCOMPPOS,GAMMA, 2,ZMB, HMAX,SMAXPRE,
     $	                EPSH01Q,EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,
     $	                IRUS,JRUS, ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 
     $	                0.0D0,0,NKVAR, IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, 
     $	                EPSMAGNLOK, NFCAL, DHFAKT,RFAKT, 
     $	                ICWORKLIM0,ICWORKDIM,CPRESSW)
	                IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	              END IF
	              NFCALLS(3) = NFCALLS(3) + NFCAL
	            END IF
	          END IF
	        END IF
	        AIMSLUT = DMIN1 ( AIMEND(JP), AIMCUT )
	      END IF
	      IF ( NPU.NE.-2 .AND. AIMSLUT.GT.DIMAG(UPBPA) ) THEN
	        IPS = 1
	        IF ( AIMSLUT .GT. DIMAG(UPBPA)+0.01D0 ) THEN
	          ZMB(1) = DCMPLX(DREAL(UPBPA),DIMAG(UPBPA)+0.01D0)
	          ZMB(2) = DCMPLX(DREAL(UPBPA),AIMSLUT)
	          IF (IADINTXFI.EQ.0) THEN
	            CALL TRAPEZFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	            RCOMPPOS,GAMMA, 2,ZMB, HMAX, NHYDLI,CPRESSW, 
     $	            NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	          ELSE
	            CALL ADINTXFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	            RCOMPPOS,GAMMA, 2,ZMB, HMAX,SMAXPRE,EPSH01Q,
     $	            EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, 
     $	            ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, 
     $	            IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, EPSMAGNLOK, 
     $	            NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	            IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	          END IF
	          NFCALLS(3) = NFCALLS(3) + NFCAL
	        END IF
	        IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	          IFILONLOK = 0     ! adapted to small ss
	        ELSE
	          IFILONLOK = IFILONPRE
	        END IF
	        SSZ(1) = (0.0D0,0.0D0)    ! up=upbpa
	        SSZ(2) = DCMPLX ( 
     $	        DMIN1(0.1D0,DSQRT(AIMSLUT-DIMAG(UPBPA))), 0.0D0 )
	        IF (IADINTXFI.EQ.0) THEN
	          CALL TRAPEZFI (NHYD,LI, DHBRANCHSS,RBRANCHSS, 
     $	          IFILONLOK,RBRANCHPOSSS,GAMMA, 2,SSZ, HMAX, 
     $	          NHYDLI,CPRESSW, NFCAL, DHFAKT,RFAKT, 
     $	          ICWORKLIM0,ICWORKDIM,CPRESSW)
	        ELSE
	          CALL ADINTXFI (NHYD,LI, DHBRANCHSS,RBRANCHSS, 
     $	          IFILONLOK,RBRANCHPOSSS,GAMMA, 2,SSZ, HMAX,
     $	          SMAXPRE,EPSH01Q,EPS,NELMAXPRE,NBACKMAX, IXTYP,
     $	          NRUS,IRUS,JRUS, ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 
     $	          0.0D0,0,NKVAR, IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, 
     $	          EPSMAGNLOK, NFCAL, DHFAKT,RFAKT, 
     $	          ICWORKLIM0,ICWORKDIM,CPRESSW)
	          IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	        END IF
	        NFCALLS(3) = NFCALLS(3) + NFCAL
	      END IF
	      IF ( NPU.NE.-1 .AND. AIMSLUT.GT.DIMAG(UPBPB) ) THEN
	        IPS = 2
	        IF ( AIMSLUT .GT. DIMAG(UPBPB)+0.01D0 ) THEN
	          ZMB(1) = DCMPLX(DREAL(UPBPB),DIMAG(UPBPB)+0.01D0)
	          ZMB(2) = DCMPLX(DREAL(UPBPB),AIMSLUT)
	          IF (IADINTXFI.EQ.0) THEN
	            CALL TRAPEZFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	            RCOMPPOS,GAMMA, 2,ZMB, HMAX, NHYDLI,CPRESSW, 
     $	            NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	          ELSE
	            CALL ADINTXFI (NHYD,LI, DHBRANCH,RCOMP, IFILON,
     $	            RCOMPPOS,GAMMA, 2,ZMB, HMAX,SMAXPRE,EPSH01Q,
     $	            EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, 
     $	            ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, 
     $	            IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, EPSMAGNLOK, 
     $	            NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	            IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	          END IF
	          NFCALLS(3) = NFCALLS(3) + NFCAL
	        END IF
	        IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	          IFILONLOK = 0     ! adapted to small ss
	        ELSE
	          IFILONLOK = IFILONPRE
	        END IF
	        SSZ(1) = (0.0D0,0.0D0)    ! up=upbpb
	        SSZ(2) = DCMPLX ( 
     $	        DMIN1(0.1D0,DSQRT(AIMSLUT-DIMAG(UPBPB))), 0.0D0 )
	        IF (IADINTXFI.EQ.0) THEN
	          CALL TRAPEZFI (NHYD,LI, DHBRANCHSS,RBRANCHSS, 
     $	          IFILONLOK,RBRANCHPOSSS,GAMMA, 2,SSZ, HMAX, 
     $	          NHYDLI,CPRESSW, NFCAL, DHFAKT,RFAKT, 
     $	          ICWORKLIM0,ICWORKDIM,CPRESSW)
	        ELSE
	          CALL ADINTXFI (NHYD,LI, DHBRANCHSS,RBRANCHSS, 
     $	          IFILONLOK,RBRANCHPOSSS,GAMMA, 2,SSZ, HMAX,
     $	          SMAXPRE,EPSH01Q,EPS,NELMAXPRE,NBACKMAX, IXTYP,
     $	          NRUS,IRUS,JRUS, ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 
     $	          0.0D0,0,NKVAR, IADINTXFI,TEXTLOG, NHYDLI,CPRESSW, 
     $	          EPSMAGNLOK, NFCAL, DHFAKT,RFAKT, 
     $	          ICWORKLIM0,ICWORKDIM,CPRESSW)
	          IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	        END IF
	        NFCALLS(3) = NFCALLS(3) + NFCAL
	      END IF
	    END IF

	  END IF
	END DO

C	-------- analytic part finally --------------------------------------

	IF (INTREST.NE.0) THEN
C	  add Hankel-transforms of "dhcomp with idhtyp=1" analytically
	  IF (KFWTYP.NE.4) STOP 566144     ! not implemented
	  DO IR = 1, LI
	    IHR = IR
	    CSLASK = (RK0*OR(IR))**2
	    DO IH = 1, NHYD
	      CPRESSW(IHR) = CPRESSW(IHR) + 
     $	      ASFAKT(IH)/CDSQRT(ASEXP(IH)**2+CSLASK)
	      IHR = IHR + LI
	    END DO
	  END DO
	END IF

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

	RETURN
	END

