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 ) 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 'INMODEL.INC' INCLUDE 'INRANGE.INC' INCLUDE 'INPARAM.INC' INCLUDE 'INPINT.INC' INCLUDE 'INPINTX.INC' INCLUDE 'INPINTASX.INC' INCLUDE 'PREDHAS.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 EXTERNAL DHCOMP,DHCOMPCIRC, RCOMP,RCOMPPOS, R1COMP,R1COMPPOS, $ DHBRANCHSUM, DHBRANCH, DHBRANCHSS, RBRANCHSS,RBRANCHPOSSS, $ CDETSUB, RCOMPONE,RCOMPPOSONE DATA PI /3.141592653589793D0/ C --------------------------------------------------------------------- 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.0) THEN IF (NSOURCFW.NE.-2) STOP 34441 ! special for 'multi' RW0 = 0.5D0 * RW0 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 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. $ SLASK.NE.DIMAG(OMEG) .OR. I.NE.IUROT ) STOP 34155 IF (XYONENORM(JP).GT.0.0D0) THEN C integration around the pole (positive orientation) 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 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" 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,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,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,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,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,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,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,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,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,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 SUBROUTINE TERFUN ( N, TERM, ICWORKLIM0,ICWORKDIM,CWORK ) C must be called initially (outside euleracc) with n=-1 IMPLICIT NONE INCLUDE 'RPPARAMETER.INC' C nol <= noldim , nhyd <= nhyddim , li <= lidim , INCLUDE 'INMODEL.INC' INCLUDE 'INRANGE.INC' INCLUDE 'INPINT.INC' INCLUDE 'INTAIL.INC' INCLUDE 'INTWORK.INC' INTEGER N COMPLEX*16 TERM(*) INTEGER ICWORKLIM0,ICWORKDIM COMPLEX*16 CWORK(*) INTEGER IH, NKVAR, NFCAL, ISLASK REAL*8 SLASK COMPLEX*16 CSLASK, $ OLDCOS(nhyddim), UPHTAIL(2) EXTERNAL DHCOMP, R1COMP,R1COMPPOS SAVE OLDCOS IF (N.GE.0) THEN DO IH = 1, NHYD ISLASK = 2*IH TERM(ISLASK-1) = OLDCOS(IH) TERM(ISLASK) = (0.0D0,0.0D0) END DO IF (IT.EQ.1) THEN UPHTAIL(1) = UPH1AS + N*UPPID UPHTAIL(2) = UPHTAIL(1) + UPPIDH TEXTLOG(9:12) = 'righ' ELSE UPHTAIL(2) = UPH1AS - N*UPPID UPHTAIL(1) = UPHTAIL(2) - UPPIDH TEXTLOG(9:12) = 'left' END IF IF (IADINTXFI.EQ.0) THEN CALL TRAPEZFI (NHYD,2, DHCOMP,R1COMP, IFILON, $ R1COMPPOS,GAMMAS, 2,UPHTAIL, HMAXAS, NHYD2,TERM, $ NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CWORK) ELSE 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,TERM, SLASK, NFCAL, DHFAKT,RFAKT, $ ICWORKLIM0,ICWORKDIM,CWORK) END IF IF (IT.EQ.1) THEN NFCALLS(5) = NFCALLS(5) + NFCAL ELSE NFCALLS(6) = NFCALLS(6) + NFCAL END IF DO IH = 1, NHYD ISLASK = 2*IH-1 OLDCOS(IH) = TERM(ISLASK) ! temporary storage TERM(ISLASK) = (0.0D0,0.0D0) END DO ELSE IF (N.NE.-1) STOP 741 DO IH = 1, NHYD2 TERM(IH) = (0.0D0,0.0D0) END DO END IF IF (IT.EQ.1) THEN UPHTAIL(1) = UPH1AS + (N+0.5D0)*UPPID ! uph0as if n=-1 UPHTAIL(2) = UPHTAIL(1) + UPPIDH ! uph1as if n=-1 TEXTLOG(9:12) = 'righ' ELSE UPHTAIL(2) = UPH1AS - (N+0.5D0)*UPPID ! uph0as if n=-1 UPHTAIL(1) = UPHTAIL(2) - UPPIDH ! uph1as if n=-1 TEXTLOG(9:12) = 'left' END IF IF (IADINTXFI.EQ.0) THEN CALL TRAPEZFI (NHYD,2, DHCOMP,R1COMP, IFILON, $ R1COMPPOS,GAMMAS, 2,UPHTAIL, HMAXAS, NHYD2,TERM, $ NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CWORK) ELSE 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,TERM, SLASK, NFCAL, DHFAKT,RFAKT, $ ICWORKLIM0,ICWORKDIM,CWORK) END IF IF (IT.EQ.1) THEN NFCALLS(5) = NFCALLS(5) + NFCAL ELSE NFCALLS(6) = NFCALLS(6) + NFCAL END IF IF (N.GE.0) THEN DO IH = 1, NHYD ISLASK = 2*IH-1 CSLASK = TERM(ISLASK) TERM(ISLASK) = OLDCOS(IH) OLDCOS(IH) = CSLASK END DO ELSE DO IH = 1, NHYD OLDCOS(IH) = TERM(2*IH-1) END DO END IF RETURN END