
C      ************************************************************
C      *                                                          *
C      * program for calculation of frequency response            *
C      *    or synthetic time traces                              *
C      * horizontally layered earth model                         *
C      *    solid layers on top (possibly none)                   *
C      *    fluid layers together "in between"                    *
C      *    solid layers and a solid half-space below             *
C      *                                                          *
C      * developed starting from Rainer Kind's reflectivity       *
C      *    program for synthetic surface seismograms for a       *
C      *    solid structure                                       *
C      *                                                          *
C      * Sven Ivansson 1990                                       *
C      *                                                          *
C      ************************************************************





	SUBROUTINE RPRESSMAIN ( LUNREAD,LUNWRITE, LUNINSGN, 
     $  LX,LITOP,ICWORKDIM, EX,PM0GOR,EPS0, CR9T1E,CPRESSW, CPLOK,CP )
C       positive luninsgn from rpress, negative luninsgn from rpressfw

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       nol <= noldim , nhyd <= nhyddim , li <= lidim (li:==litop)
	INTEGER NFRDIM, NTAB,JTAB, NPARTSDIM, NPUDIM
	PARAMETER ( NFRDIM = 2048, NTAB = 17, JTAB = 9, NPARTSDIM = 7, 
     $	NPUDIM = 100 )
C	nfr <= nfrdim , nrus0 <= ntab , 
C	min0(nrus0,jrus0) <= jtab , nparts <= npartsdim , 
C	npu(1:nparts) <= npudim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INRANGE.INC'
	INCLUDE 'INPARAM.INC'
	INCLUDE 'DPCFSUB.INC'
	INCLUDE 'INPINT.INC'
	INCLUDE 'INPINTX.INC'
	INCLUDE 'INPINTASX.INC'
	INCLUDE 'HARKLIQ.INC'
	INCLUDE 'INTWORK.INC'
	INCLUDE 'IADR00.INC'
	INCLUDE 'INSPEC.INC'
	INCLUDE 'SOURCFW.INC'     ! to check nsourcfw
        INCLUDE 'fw/FWPARAM.INC'  ! to write iystab

	INTEGER LUNREAD,LUNWRITE, LUNINSGN, LX,LITOP, ICWORKDIM
	REAL*8 EX(lx), PM0GOR(litop),EPS0(litop)
	COMPLEX*16 CR9T1E(litop), CPRESSW(icworkdim)
	COMPLEX CPLOK(lx), CP(lx,litop)

	INTEGER LUNIN, MHQ0,MHQNOW,MHQEFF, IXTYP0,NRUS0,IRUS0,JRUS0, NRUS0PRE, 
     $  I0LOG, NFR0,NFR, IUROT, IDB, INTREST, IFILONPRE, KA,NX,NY, 
     $  NPARTS, NHYDLI, ICWORKLIM0, IPOSAM,IPOSPH, IPOSTR, IPOS, 
     $  ID, IH,IR, IHR, JF, JP,JL, J1,J2, M, J, I, 
     $	IOUT, ISLASK, IXROW, JCOLX, IEXPX, NPTYP, 
     $	IA0_TABNEW,IA0_VALTERM,IA0_TAB,IA0_YTAB, MHQ(ntab), 
     $	IBESSARR(npartsdim), NPUARR(npartsdim), IPOSMOD(npartsdim), 
     $	NFCALLSOLD(6)
	REAL*8 FREQ,FREQMHQ0, HSTEPFSUBPRE, HSTEPFSUB0, PREEPS0, 
     $  FREQMIN,FREQINT, TMIN,DAUER,URED,  FREQDOM, PREALA,ALA,EX1, RINT, 
     $  T0,TMO, OMN, DDT, DDF, DDOM, XMAX,XGLUE,YMAX,YGLUE, 
     $	AM0,PM0FAKT,PM0, PM01M, HMAXDH,HMAXRELWP, HKVARPART, 
     $	PREFREL1,PREFREL2, POSTEPS,POSTEPSTERM,POSTEPSCONT,POSTEPSEU, 
     $	POSTEPS0, SLASK, SLASKD, PI, DATAN2DQ, 
     $	ANGRKTAIL(2), XYONENORM(npartsdim), AIMEND(npartsdim), 
     $	OFREQ(nfrdim), ACOEF(jtab-1,ntab), HQI(jtab,2:ntab), 
     $	FIHQI(jtab),QHQI(3)    ! used as dummies only
	COMPLEX*16 FSOURC, CR9,T1E, CSLASK, CDSQRT12, AZIMILN, 
     $	UPHOERN(npudim,npartsdim), UPHOERNLOK(npudim)
	LOGICAL WTXT, OKSUB
	CHARACTER*32 FILM,BESSFIL, FILO, CR9FIL, FILG, MAINAM,FILAM, 
     $	MAINPH,FILPH, MAINTR,FILTR, MAINMOD(npartsdim), 
     $	FILMOD(npartsdim)

	EXTERNAL IA0_TABNEW,IA0_VALTERM,IA0_TAB,IA0_YTAB

	DATA PI /3.141592653589793D0/

	IF (LUNWRITE.EQ.6) THEN
	  WTXT = .TRUE.
	ELSE
	  WTXT = .FALSE.
	END IF
	LUNIN = IABS(LUNINSGN)
	IF (LUNIN.EQ.0) STOP 34155

	UPREF = (0.0D0,0.0D0)
	ISHEETTYP = 1
	IDET = 0
	IHARKLIQ = 0

	IF (WTXT) WRITE(LUNWRITE,*)
	IF (WTXT) WRITE(LUNWRITE,*) 'Give infile FILM (for the model) :'
	FILM = ' '
	DO WHILE (FILM.EQ.' ')
	  READ(LUNREAD,'(A)') FILM
	END DO
	CALL BLANKTOEND (FILM)
	IF (WTXT) WRITE(LUNWRITE,*) 'Give MHQ0, FREQMHQ0 (Hz), IAZIMI :'
	READ(LUNREAD,*) MHQ0,FREQMHQ0, IAZIMI
	IF (IAZIMI.NE.0) THEN
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give ALFA1AZIM:'
	  READ(LUNREAD,*) ALFA1AZIM
	END IF
	IF (WTXT) WRITE(LUNWRITE,*) 'Give ICUTTYP ( 1="axis-hyperbolic ',
     $  'for up" 2="vertical for up" '
	IF (WTXT) WRITE(LUNWRITE,*) '3="mixture with aimcut" ) :'
	IF (WTXT) WRITE(LUNWRITE,*) '( "-5 or -15 or -45 or -60" for Muir )'
	IF (WTXT) WRITE(LUNWRITE,*) '( irrelevant for rigid or free ',
     $  'boundaries )'
	READ(LUNREAD,*) ICUTTYP
	IF (ICUTTYP.EQ.3) THEN
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give AIMCUT:'
	  READ(LUNREAD,*) AIMCUT
	END IF
	IF ( ICUTTYP.NE.1 .AND. ICUTTYP.NE.2 .AND. ICUTTYP.NE.3 .AND. 
     $	ICUTTYP.NE.-5 .AND. ICUTTYP.NE.-15 .AND. 
     $	ICUTTYP.NE.-45 .AND. ICUTTYP.NE.-60 ) STOP 340
	IF (WTXT) WRITE(LUNWRITE,*) 'Give EPSGLOBLN,EPSGLOB,EPSFAKT,',
     $  'HFAKTMAX,HOFF,EPSMAGNIFMAX :'
	IF (WTXT) WRITE(LUNWRITE,*) '( standard is  epsfakt=0.1, ',
     $  'hfaktmax=1.05, hoff=6d-16 )'
	READ(LUNREAD,*) EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,EPSMAGNIFMAX
	IF (WTXT) WRITE(LUNWRITE,*) 'Give EPSFSUB,HSTEPFSUBPRE :'
	IF (WTXT) WRITE(LUNWRITE,*) '( non-positive epsfsub implies ',
     $  '"dpclin with constant hstep without '
	IF (WTXT) WRITE(LUNWRITE,*) 'internal error control )'
	IF (WTXT) WRITE(LUNWRITE,*) '( except for fsmdzmin, hstepfsub ',
     $  'will be "hstepfsubpre/freq" which '
	IF (WTXT) WRITE(LUNWRITE,*) '  is a "compromise" regarding ',
     $  'rkra & rkrb )'
	READ(LUNREAD,*) EPSFSUB,HSTEPFSUBPRE
	IF (WTXT) WRITE(LUNWRITE,*) 'Give (main) PREEPS0, ',
     $  'IXTYP0,NRUS0PRE,IRUS0,JRUS0 :'
	IF (WTXT) WRITE(LUNWRITE,*) '( preeps0 is a relative tolerance, ',
     $  'so it becomes absolute for dB )'
	IF (WTXT) WRITE(LUNWRITE,*) '( ixtyp0 = "1,-1,0 / 2,-2 / 3,-3" for ',
     $	'"polynomial / rational / inv-poly" )'
	IF (WTXT) WRITE(LUNWRITE,*) '( "rational" not adapted to "gaps" for ',
     $	'"iifsmsome=1 & epsfsub<=0.0d0" )'
	IF (WTXT) WRITE(LUNWRITE,*) '( nrus0pre < 0  for extrapolation log ',
     $  'on units 100-102 )'
	READ(LUNREAD,*) PREEPS0, IXTYP0,NRUS0PRE,IRUS0,JRUS0
	IF (IABS(IXTYP0).GE.3) STOP 5710
	IF (NRUS0PRE.GE.0) THEN
	  I0LOG = 0
	  NRUS0 = NRUS0PRE
	ELSE
	  I0LOG = 1
	  NRUS0 = -NRUS0PRE
	END IF
	IF ( NRUS0.EQ.1 .OR. NRUS0.GT.ntab ) STOP 571
	IF (MIN0(NRUS0,JRUS0).GT.jtab) STOP 572

	IF (WTXT) WRITE(LUNWRITE,*)
	IF (LX.LE.1) THEN
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give FREQMIN(>0.0),FREQINT ',
     $    '(both in Hz) , NFR :'
	  IF (WTXT) WRITE(LUNWRITE,*) '( " freqmin<0.0 & freqint=0.0 " ',
     $    ' for  "explicit ofreq(1:nfr) " )'
	  READ(LUNREAD,*) FREQMIN,FREQINT, NFR
	  IF ( FREQMIN.EQ.0.0D0 .OR. NFR.GT.nfrdim ) STOP 6666
	  IF (FREQMIN.GE.0.0D0) THEN
	    SLASK = FREQINT/NFR
	    DO JF = 1, NFR
	      OFREQ(JF) = FREQMIN + (JF-1)*SLASK
	    END DO
	  ELSE
	    IF (FREQINT.NE.0.0D0) STOP 6660
	    IF (WTXT) WRITE(LUNWRITE,*) 'give ofreq(1:nfr) increasing'
	    READ(LUNREAD,*) (OFREQ(JF),JF=1,NFR)
	  END IF
          OMN = 2*PI*OFREQ(NFR)
	  IUROT = 0
	  ALA = 0.0D0
	ELSE
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give TMIN (s), DAUER (s), URED (s/km) :'
          READ(LUNREAD,*) TMIN,DAUER,URED
          DDT = DAUER/LX
	  DDF = 1/DAUER
          DDOM = 2*PI/DAUER
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give PREALA ( 0.0<preala<=1.0, 1.0 ',
     $    'for real omeg ) :'
	  IF (WTXT) WRITE(LUNWRITE,*) '( "preala" is "temporary reduct ',
     $    'factor during dauer" )'
	  READ(LUNREAD,*) PREALA
	  ALA = -DLOG(PREALA)/DAUER
	  EX1 = DEXP(ALA*DDT)
	  EX(1) = 1.0D0
	  DO I = 2, LX
	    EX(I) = EX1*EX(I-1)
	  END DO
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give IUROT ( 0 for "u as ',
     $    'specified", 1 for "u rotation with -arg(omeg)" ) :'
	  IF (WTXT) WRITE(LUNWRITE,*) '( u = "up rotated by -arg(omeg)" )'
	  IF (WTXT) WRITE(LUNWRITE,*) '( irrelevant if "preala"=1.0 )'
	  READ(LUNREAD,*) IUROT
	  SLASK = ( 0.5D0 + LX/2 ) / (2*DAUER)
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give FREQDOM (Hz) , < ',SNGL(SLASK),' :'
	  IF (WTXT) WRITE(LUNWRITE,*) '( negative for "source-spectrum ',
     $    'from file" )'
          READ(LUNREAD,*) FREQDOM
	  IF (FREQDOM.GT.0.0D0) THEN
	    NFR0 = 0.5D0 + (2*FREQDOM)*DAUER
	    IF (NFR0.GT.(LX/2)) STOP 440
	  ELSE
	    NFR0 = LX/2
	    IF (WTXT) WRITE(LUNWRITE,*) 'give source-spectrum file ',
     $      'CR9FIL ("plotfkn format")'
	    IF (WTXT) WRITE(LUNWRITE,*) '( obtained from "cr9rpress"; '
	    IF (WTXT) WRITE(LUNWRITE,*) '  "dauer in rpress" = "tperiod ',
     $	    'in cr9rpress"'
	    IF (WTXT) WRITE(LUNWRITE,*) '            since ',
     $      '"ddom=2*pi/tperiod" is necessary'
	    IF (WTXT) WRITE(LUNWRITE,*) '  "lx in rpress" = "iperiod in ',
     $	    'cr9rpress"'
	    IF (WTXT) WRITE(LUNWRITE,*) '            since "nfr in ',
     $      'rpress & cr9rpress" must be equal'
	    IF (WTXT) WRITE(LUNWRITE,*) '  Nyquist freq is "(lx/2)/dauer" )'
            READ(LUNREAD,'(A)') CR9FIL
	    CALL BLANKTOEND (CR9FIL)
	    OPEN(UNIT=9,FILE=CR9FIL,STATUS='OLD')
	    READ(9,*) I
	    DO J = 1, I
	      READ(9,*)
	    END DO
	    READ(9,*) I,SLASK
	    IF ( I.NE.1 .OR. SLASK.NE.PREALA ) STOP 15614
	    READ(9,*) I
	    IF ( I .NE. NFR0-1 ) STOP 15614
	  END IF
          OMN = NFR0*DDOM
	  NFR = NFR0-1
	  IF (NFR.GT.nfrdim) STOP 60
          SLASK = 1/DAUER
	  DO JF = 1, NFR
	    OFREQ(JF) = (NFR+1-JF) * SLASK      ! decreasing
	  END DO
	  IF (WTXT) WRITE(LUNWRITE,*) 'equivalent "freqmin freqint nfr"  ', 
     $    NFR*SLASK, -NFR*SLASK, NFR
	  IF (FREQDOM.GT.0.0D0) THEN
	    FREQDOM = 0.5D0 * OMN/(2*PI)
	    TMO = (6.5D0/4) * 0.5D0/FREQDOM    ! experimental shift, causality
	  ELSE
	    TMO = 0.0D0
	    DO JF = 1, NFR
	      READ(9,*) SLASK,CPLOK(NFR+1-JF)
	      IF ( DABS(1-SLASK/OFREQ(NFR+1-JF)) .GT. 0.001D0 ) 
     $	      STOP 5162
	    END DO
	    CLOSE(9)
	  END IF
	END IF

	LI = LITOP
	IF (WTXT) WRITE(LUNWRITE,*) 'Give RMIN,RINT (both in km) :'
	IF (WTXT) WRITE(LUNWRITE,*) '( " rmin<0.0 & rint=0.0 "  for  ',
     $	'"explicit or(1:li) " )'
	READ(LUNREAD,*) RMIN,RINT
	IF (LI.GT.lidim) STOP 5555
	IF (RMIN.GE.0.0D0) THEN
	  DR = RINT/LI
	  DO IR = 1, LI
	    OR(IR) = RMIN + (IR-1)*DR               !  horizontal distance
	  END DO
	ELSE
	  IF (RINT.NE.0.0D0) STOP 5550
	  IF (WTXT) WRITE(LUNWRITE,*) 'give or(1:li) increasing'
	  READ(LUNREAD,*) (OR(IR),IR=1,LI)          !  horizontal distances
	END IF
	IF (OR(1).NE.0.0D0) THEN
	  DO IR = 1, LI
	    RL(IR) = DLOG ( 0.5D0 * OR(IR) )
	  END DO
	  DO IR = 1, LI
	    RS(IR) = 0.5D0 / DSQRT ( PI * OR(IR) )
	  END DO
	END IF
	ORMINSPEC = OR(1)

	IF (WTXT) WRITE(LUNWRITE,*)
	IF (WTXT) WRITE(LUNWRITE,*) 'Give outfile FILO:'
	FILO = ' '
	DO WHILE (FILO.EQ.' ')
	  READ(LUNREAD,'(A)') FILO
	END DO
	CALL BLANKTOEND (FILO)
	OPEN(UNIT=1,FILE=FILO,STATUS='UNKNOWN')

	IF (LX.GT.1) THEN
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give "main name" MAINTR for the nhyd ',
     $	  '"outfiles-trans" :'
	  IF (WTXT) WRITE(LUNWRITE,*) '(interface numbers as extensions, ',
     $    'plot with udp)'
	  MAINTR = ' '
	  DO WHILE (MAINTR.EQ.' ')
	    READ(LUNREAD,'(A)') MAINTR
	  END DO
	  CALL BLANKTOEND (MAINTR)
	  I = 32
	  DO WHILE (MAINTR(I:I).EQ.' ')
	    I = I-1
	  END DO
	  MAINTR((I+1):(I+1)) = '.'
	  IPOSTR = I+2
	  OPEN(UNIT=21,FORM='UNFORMATTED',STATUS='SCRATCH')
	END IF

	IF (WTXT) WRITE(LUNWRITE,*) 'Give outfile FILG (for the ',
     $  'depth/distance pattern):'
	IF (LX.LE.1) THEN
	  FILG = ' '
	  DO WHILE (FILG.EQ.' ')
	    READ(LUNREAD,'(A)') FILG
	  END DO
	ELSE
	  READ(LUNREAD,'(A)') FILG
	END IF
	CALL BLANKTOEND (FILG)
	IF (FILG.NE.' ') THEN
	  OPEN(UNIT=23,FILE=FILG,STATUS='UNKNOWN')
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give IDB (1 for dB, 0 for quotient, ',
     $	  '-1 for sngl, -2 for dble, -3 for "Lame-lambda mult") :'
	  READ(LUNREAD,*) IDB
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give "main name" MAINAM for the nfr ',
     $	  '"outfiles-ampl" :'
	  IF (WTXT) WRITE(LUNWRITE,*) '(freq indices as extensions, plot ',
     $    'with plev2)'
	  READ(LUNREAD,'(A)') MAINAM
	  CALL BLANKTOEND (MAINAM)
	  I = 32
	  DO WHILE (MAINAM(I:I).EQ.' ')
	    I = I-1
	  END DO
	  MAINAM((I+1):(I+1)) = '.'
	  IPOSAM = I+2
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give "main name" MAINPH for the nfr ',
     $	  '"outfiles-phase" :'
	  IF (WTXT) WRITE(LUNWRITE,*) '(freq indices as extensions, plot ',
     $    'with plev2)'
	  READ(LUNREAD,'(A)') MAINPH
	  CALL BLANKTOEND (MAINPH)
	  I = 32
	  DO WHILE (MAINPH(I:I).EQ.' ')
	    I = I-1
	  END DO
	  MAINPH((I+1):(I+1)) = '.'
	  IPOSPH = I+2
	ELSE
	  MAINAM = ' '
	  MAINPH = ' '
	END IF

	IF (WTXT) WRITE(LUNWRITE,*)
	IF (WTXT) WRITE(LUNWRITE,*) 'Give INTREST (nonzero for "simultaneous ',
     $	'asymptotic tail handling") :'
	IF (WTXT) WRITE(LUNWRITE,*) '( note that the "rmodel iteration" ',
     $  'may not give a development '
	IF (WTXT) WRITE(LUNWRITE,*) '  in even powers of "h" if ',
     $  '" (intrest is nonzero) and (iconttyp or '
	IF (WTXT) WRITE(LUNWRITE,*) '  some itail is zero) " )'
	READ(LUNREAD,*) INTREST
	IF (WTXT) WRITE(LUNWRITE,*) 'Give IFILONPRE (for appropriate ',
     $	'"ifilon:s" to be set) :'
	IF (WTXT) WRITE(LUNWRITE,*) '( 0 for "trapez non-Filon" ; ',
     $  '1,2 for "trapez Filon" ; '
	IF (WTXT) WRITE(LUNWRITE,*) '-1,-2 for "Simpson Filon" ; ',
     $  '-3 for "Simpson non-Filon" )'
	READ(LUNREAD,*) IFILONPRE
	CALL BESSRT (LUNREAD,9, BESSFIL,KA,XMAX,NX,XGLUE,YMAX,NY,YGLUE)
	IF (WTXT) WRITE(LUNWRITE,*) 'Give IADINTXFI  ( 0 for "trapezfi", ',
     $  'nonzero for "lunout" to "iadintxfi" ) :'
	IF (WTXT) WRITE(LUNWRITE,*) '( iadintxfi "= 6 or >= 80"  for ',
     $  'integration-step histogram '
	IF (WTXT) WRITE(LUNWRITE,*) '  iadintxfi "= 66, 77 or 88"  ',
     $  'for "detailed tail logs" as well )'
	READ(LUNREAD,*) IADINTXFI
	IF (IADINTXFI.NE.0) THEN
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give (adintxfi) SMAXPRE,',
     $    'NELMAXPRE,NBACKMAX:'
	  READ(LUNREAD,*) SMAXPRE,NELMAXPRE,NBACKMAX
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give (adintxfi) IXTYP,NRUS,IRUS,JRUS:'
	  IF (WTXT) WRITE(LUNWRITE,*) '( ixtyp = "1,-1,0 / 2,-2 4 -4 / ',
     $    '3,-3" for "polyn / rational / inv-poly" )'
	  IF (WTXT) WRITE(LUNWRITE,*) '( "rational cases 2,-2" not ',
     $    'adapted to "gap" for ifilonpre<0 )'
	  READ(LUNREAD,*) IXTYP,NRUS,IRUS,JRUS
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give (adintxfi) ARGDIFLIM1,',
     $    'ARGDIFLIM0,ARGDIFLIM2 (degrees) :'
	  READ(LUNREAD,*) ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2
	  TEXTLOG = ' rpress      filo='//FILO
	END IF

	IF (WTXT) THEN
	  WRITE(LUNWRITE,*)
	  WRITE(LUNWRITE,*) 'The Bessel-function arguments will be ',
     $	  '"omeg*u*r" ;'
	  WRITE(LUNWRITE,*) 'hence, the asymptotic-Bessel ',
     $    'u-periods will be  " 1 / ( (omeg/(2*pi)) * r )"'
	  WRITE(LUNWRITE,*) '   ofreq(1),ofreq(nfr)=',SNGL(OFREQ(1)),
     $	  SNGL(OFREQ(NFR))
	  WRITE(LUNWRITE,*) '                   ala=',ALA
	  WRITE(LUNWRITE,*) '   or(1),or(li)=',SNGL(OR(1)),SNGL(OR(LI)),' .'
	  WRITE(LUNWRITE,*) 'For IBESS... below <unsuitable ',
     $    'with "0.5*H0(1)" close to the origin > :'
	  WRITE(LUNWRITE,*) '(  1="J0-bessra"          2="J0-besbod" '
	  WRITE(LUNWRITE,*) '  -1="0.5*H0(1)-bessra"  ',
     $    '-2="0.5*H0(1)-besbod"  -3="0.5*H0(1)-quick" )'
	  WRITE(LUNWRITE,*) '( however, dbes0rat & dhan0rat ',
     $    'preferred on "real-rk" path parts '
	  WRITE(LUNWRITE,*) '  except for ibess...=-3 )'
	END IF

	IF (WTXT) THEN
	  WRITE(LUNWRITE,*)
	  WRITE(LUNWRITE,*) 'Give NPARTS ( nr of "path-mode-branch parts" ',
     $	  'except "tails" ) :'
	  WRITE(LUNWRITE,*) '( order of treatment later on: "jp=nparts,1,-1" )'
	  WRITE(LUNWRITE,*) '(   npu>0 for "integration polygon" where ',
     $	  '"intrest" relevant )'
	  WRITE(LUNWRITE,*) '(   npu=0 for "modes" )'
	  WRITE(LUNWRITE,*) '(   npu<0 for branch-cut integral(s) '
	  WRITE(LUNWRITE,*) '          -1 for P-branch , -2 for S-branch , ',
     $	  '-12 for both branches )'
	END IF
	READ(LUNREAD,*) NPARTS
	IF (NPARTS.GT.npartsdim) STOP 4315
	NPTYP = 0
	DO JP = 1, NPARTS
	  IF (WTXT) WRITE(LUNWRITE,*) '  part jp=',JP
	  IF (WTXT) WRITE(LUNWRITE,*) '    give IBESSARR(jp),NPUARR(jp)'
	  READ(LUNREAD,*) IBESSARR(JP),NPUARR(JP)
	  IF (OR(1).EQ.0.0D0) THEN
	    IF (IBESSARR(JP).LT.0) STOP 4500
	    IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	      IF (IBESSARR(JP).EQ.1) STOP 4501  ! "rl" needed in rcomppos
	    END IF
	  END IF
	  IF ( IBESSARR(JP).EQ.-3 .AND. RMIN.LE.0.0D0 ) STOP 416
	  IF (NPUARR(JP).GT.npudim) STOP 4316
	  IF (NPUARR(JP).GT.0) THEN
	    IF (WTXT) THEN 
	      WRITE(LUNWRITE,*) '    give complex UPHOERN(1:npu), ',
     $        'all in s/km'
	      WRITE(LUNWRITE,*) '    ( real part +-991 for +-Re(upbpa)'
	      WRITE(LUNWRITE,*) '      real part +-992 for +-Re(upbpb)'
	      WRITE(LUNWRITE,*) '      with "upref" locally from ',
     $	      '"uphoern(1)&(2)" )'
	    END IF
	    READ(LUNREAD,*) (UPHOERN(JL,JP),JL=1,NPUARR(JP))
	    NPTYP = 2
	  ELSE IF (NPUARR(JP).EQ.0) THEN
C	    in the normal case here, "ibessarr(jp).lt.0"
	    IF (NRUS0.NE.0) STOP 473727  ! no "outer extrap" for "mode part"
	    IF (INTREST.NE.0) STOP 3415561 ! not to be used for "mode part"
	    IF (WTXT) WRITE(LUNWRITE,*) '    give "main name" MAINMOD for ',
     $      'mode infile'
	    READ(LUNREAD,'(A)') MAINMOD(JP)
	    CALL BLANKTOEND (MAINMOD(JP))
	    I = 32
	    DO WHILE (MAINMOD(JP)(I:I).EQ.' ')
	      I = I-1
	    END DO
	    MAINMOD(JP)((I+1):(I+1)) = '.'
	    IPOSMOD(JP) = I+2
	    IF (WTXT) WRITE(LUNWRITE,*) '    give XYONENORM (pos. for ',
     $      'radius norm, '
	    IF (WTXT) WRITE(LUNWRITE,*) '    zero for "rmodpfw", neg. ',
     $      'for "dhcomp")'
	    READ(LUNREAD,*) XYONENORM(JP)
	    IF (NPTYP.EQ.0) NPTYP = 1
	  ELSE
C	    in the normal case here, "ibessarr(jp).lt.0"
	    IF (INTREST.NE.0) STOP 3415562 ! not to be used for "bransh part"
	    IF (WTXT) WRITE(LUNWRITE,*) '    give positive AIMEND (in s/km)'
	    READ(LUNREAD,*) AIMEND(JP)
	    IF (NPTYP.EQ.0) NPTYP = 1
	  END IF
	END DO
	IF (NPTYP.NE.0) THEN
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give HMAXDH,HMAXRELWP:'
	  READ(LUNREAD,*) HMAXDH,HMAXRELWP
	  IF (IADINTXFI.NE.0) THEN
	    IF (WTXT) THEN
	      WRITE(LUNWRITE,*) 'Give (main/adintxfi) EPSH01Q,PREEPS:'
	      WRITE(LUNWRITE,*) '( preeps is a relative tolerance, so it ',
     $	      'becomes absolute for dB )'
	      WRITE(LUNWRITE,*) '( preeps is for each "path part" (except ',
     $	      '"tails") separately )'
	    END IF
	    READ(LUNREAD,*) EPSH01Q,PREEPS
	    IF (WTXT) THEN 
	      WRITE(LUNWRITE,*) 'Give (main/adintxfi) HKVARPART,LUNKVAR,',
     $	      'FAKTKVAR:'
	      WRITE(LUNWRITE,*) '( hkvar will be "hkvarpart*hmax" )'
	      WRITE(LUNWRITE,*) '( faktkvar typically 1.0d0 , but negative',
     $	      ' to suppress arg-check )'
	    END IF
	    READ(LUNREAD,*) HKVARPART,LUNKVAR,FAKTKVAR
	    IF (HKVARPART.GT.0.0D0) THEN
	      IF (HKVARPART.GE.1.0D0) STOP 3001
	      IF (FAKTKVAR.GT.0.0D0) THEN
	        IF (WTXT) WRITE(LUNWRITE,*) 'Give DDSUPMIN,DDSUPMAX,',
     $          'DDSUPFAKTMAX:'
	        IF (WTXT) WRITE(LUNWRITE,*) '( typically  0.0d0 0.05d0 ',
     $          '5.0d0 )'
	        READ(LUNREAD,*) DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX
	        IF (WTXT) WRITE(LUNWRITE,*) 'Give DERIVDRELMAX,DRELMAX,',
     $          'DRELGOALFAKT:'
	        IF (WTXT) WRITE(LUNWRITE,*) '( typically  0.12d0 0.4d0 ',
     $          '0.75d0 )'
	        READ(LUNREAD,*) DERIVDRELMAX,DRELMAX,DRELGOALFAKT
	        IF (WTXT) WRITE(LUNWRITE,*) 'Give ARGDIFSTORE ',
     $          '(typically 150.0d0) :'
	        READ(LUNREAD,*) ARGDIFSTORE
	      END IF
	    END IF
	  END IF
	END IF

	IF (WTXT) WRITE(LUNWRITE,*)
	IF (WTXT) WRITE(LUNWRITE,*) 'Give ICONTTYP ( 0="no tails" , ',
     $  '1="continue" , 2="Euler tails" ) :'
	READ(LUNREAD,*) ICONTTYP
	IF (ICONTTYP.NE.0) THEN
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give ITAIL(1) (rightwards) and ',
     $    'ITAIL(2) (leftwards) :'
	  IF (INTREST.EQ.0) THEN
	    IF (WTXT) WRITE(LUNWRITE,*) '( 0 for "no tail", 1 for ',
     $      '"asymp idhtyp tail", 2 for "usual idhtyp tail" )'
	  ELSE
	    IF (WTXT) WRITE(LUNWRITE,*) '( 0 for "no additional tail", ',
     $      '2 for additional tail" )'
	  END IF
	  READ(LUNREAD,*) ITAIL(1),ITAIL(2)
	  IF (ITAIL(1).NE.0) THEN
	    IF ( INTREST.NE.0 .AND. ITAIL(1).NE.2 ) STOP 200
	    IF (WTXT) WRITE(LUNWRITE,*) 'Give  IBESSTAIL(1) , ',
     $      '0 <= ANGRKTAIL(1) < 90 , UP0TAIL(1) :'
	    READ(LUNREAD,*) IBESSTAIL(1),ANGRKTAIL(1),UP0TAIL(1)
	    IF (OR(1).EQ.0.0D0) THEN
	      IF (IBESSTAIL(1).LT.0) STOP 4500
	      IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	        IF (IBESSTAIL(1).EQ.1) STOP 4501  ! "rl" needed in rcomppos
	      END IF
	    END IF
	    IF ( IBESSTAIL(1).EQ.-3 .AND. RMIN.LE.0.0D0 ) STOP 416
	    IF (ANGRKTAIL(1).LT.0.0D0 .OR. 90.0D0.LE.ANGRKTAIL(1)) STOP 3333
	    UPRTAIL(1) = DCMPLX ( 1.0D0 , DTAND(ANGRKTAIL(1)) )
	  END IF
	  IF (ITAIL(2).NE.0) THEN
	    IF ( INTREST.NE.0 .AND. ITAIL(2).NE.2 ) STOP 200
	    IF (WTXT) WRITE(LUNWRITE,*) 'Give  IBESSTAIL(2) , ',
     $      '-90 < ANGRKTAIL(2) <= 0 , UP0TAIL(2) :'
	    READ(LUNREAD,*) IBESSTAIL(2),ANGRKTAIL(2),UP0TAIL(2)
	    IF (OR(1).EQ.0.0D0) THEN
	      IF (IBESSTAIL(2).LT.0) STOP 4500
	      IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
	        IF (IBESSTAIL(2).EQ.1) STOP 4501  ! "rl" needed in rcomppos
	      END IF
	    END IF
	    IF ( IBESSTAIL(2).EQ.-3 .AND. RMIN.LE.0.0D0 ) STOP 416
	    IF (ANGRKTAIL(2).LE.-90.0D0 .OR. 0.0D0.LT.ANGRKTAIL(2)) STOP 3333
	    UPRTAIL(2) = DCMPLX ( 1.0D0 , DTAND(ANGRKTAIL(2)) )
	  END IF
	  IF (WTXT) WRITE(LUNWRITE,*) 'Give ANPHPMINTERM:'
	  READ(LUNREAD,*) ANPHPMINTERM
	  IF (IADINTXFI.NE.0) THEN
C	    for the tails, "epsh01q=hkvar=0.0d0"
	    IF (WTXT) WRITE(LUNWRITE,*) 'Give (adintxfi) PREEPSTERM:'
	    READ(LUNREAD,*) PREEPSTERM
	  END IF
	  IF (ICONTTYP.EQ.1) THEN
	    IF (WTXT) WRITE(LUNWRITE,*) 'Give (cont test) IRCONT,PREEPSCONT:'
	    READ(LUNREAD,*) IRCONT,PREEPSCONT
	    IF ( IRCONT.LT.1 .OR. LI.LT.IRCONT .OR. 
     $	    OR(IRCONT).EQ.0.0D0 ) STOP 5002
	  ELSE
	    IF (WTXT) WRITE(LUNWRITE,*) 'Give (euleracc) KMAXEU,PREEPSEU:'
	    READ(LUNREAD,*) KMAXEU,PREEPSEU
	  END IF
	END IF

	IF (WTXT) WRITE(LUNWRITE,*)
	IF (WTXT) WRITE(LUNWRITE,*)

	WRITE(1,*)
	IF (LUNINSGN.GT.0) THEN
	  WRITE(1,*) ' OUTFILE "FILO" FROM RPRESS'
	  WRITE(1,*) ' --------------------------'
	ELSE
	  WRITE(1,*) ' OUTFILE "FILO" FROM RPRESSFW'
	  WRITE(1,*) ' ----------------------------'
	END IF
	WRITE(1,*)
	WRITE(1,*) ' KFWTYP: ',KFWTYP
	WRITE(1,*) ' IYSTAB: ',IYSTAB
	WRITE(1,*) ' LX,LI,ICWORKDIM: ',LX,LI,ICWORKDIM
	WRITE(1,*)
	WRITE(1,*) ' FILM: ',FILM
	WRITE(1,*) ' MHQ0,FREQMHQ0, IAZIMI: ',MHQ0,SNGL(FREQMHQ0), IAZIMI
	IF (IAZIMI.NE.0) WRITE(1,*) '   ALFA1AZIM: ',SNGL(ALFA1AZIM)
	IF (ICUTTYP.NE.3) THEN
	  WRITE(1,*) ' ICUTTYP: ',ICUTTYP
	ELSE
	  WRITE(1,*) ' ICUTTYP,AIMCUT: ',ICUTTYP,SNGL(AIMCUT)
	END IF
	WRITE(1,*) ' EPSGLOBLN,EPSGLOB: ',SNGL(EPSGLOBLN),SNGL(EPSGLOB)
	WRITE(1,*) ' EPSFAKT,HFAKTMAX: ',SNGL(EPSFAKT),SNGL(HFAKTMAX)
	WRITE(1,*) ' HOFF,EPSMAGNIFMAX: ',SNGL(HOFF),SNGL(EPSMAGNIFMAX)
	WRITE(1,*) ' EPSFSUB,HSTEPFSUBPRE: ',SNGL(EPSFSUB),
     $	SNGL(HSTEPFSUBPRE)
	WRITE(1,*) ' PREEPS0, IXTYP0 : ',SNGL(PREEPS0),IXTYP0
	WRITE(1,*) ' NRUS0PRE,IRUS0,JRUS0 : ',NRUS0PRE,IRUS0,JRUS0

	WRITE(1,*)
	IF (LX.LE.1) THEN
	  WRITE(1,*) ' "continuous wave"'
	  IF (FREQMIN.GT.0.0D0) THEN
	    WRITE(1,*) '   FREQMIN,FREQINT, NFR: ',SNGL(FREQMIN),
     $	    SNGL(FREQINT),'  ',NFR
	  ELSE
	    WRITE(1,*) '   NFR: ',NFR,'   "explicit ofreq(1:nfr)"'
	  END IF
	ELSE
	  WRITE(1,*) ' "transient"    LX =',LX
	  WRITE(1,*) '   TMIN,DAUER:   ',SNGL(TMIN),SNGL(DAUER)
	  WRITE(1,*) '   URED:         ',SNGL(URED)
	  WRITE(1,*) '   PREALA,IUROT: ',SNGL(PREALA),IUROT
	  WRITE(1,*) '   FREQDOM:      ',SNGL(FREQDOM)
	  IF (FREQDOM.LT.0.0D0) WRITE(1,*) '       CR9FIL ',CR9FIL
	  WRITE(1,*) '   MAINTR:       ',MAINTR
	  WRITE(1,*) '   ( NFR,DDF are ',NFR,SNGL(DDF),' )'
	END IF
	WRITE(1,*)
	IF (RMIN.GE.0.0D0) THEN
	  WRITE(1,*) ' RMIN,RINT, LI: ',SNGL(RMIN),SNGL(RINT),'  ',LI
	ELSE
	  WRITE(1,*) ' LI: ',LI,'   "explicit or(1:li)"'
	END IF
	WRITE(1,*)
	WRITE(1,*) ' Out FILO     ',FILO,'         ! this file'
	IF (LX.GT.1) WRITE(1,*) '     MAINTR   ',MAINTR
	WRITE(1,*) ' Out FILG    ',FILG
	IF (FILG.NE.' ') THEN
	  WRITE(1,*) '     MAINAM   ',MAINAM,'      IDB=',IDB
	  WRITE(1,*) '     MAINPH   ',MAINPH,'      IDB=',IDB
	END IF
	WRITE(1,*)
	WRITE(1,*) ' INTREST: ',INTREST
	WRITE(1,*) ' IFILONPRE: ',IFILONPRE
	WRITE(1,*) ' BESSFIL: ','     ',BESSFIL
	IF (BESSFIL.NE.'%') THEN
	  WRITE(1,*) '   KA : ',KA
	  WRITE(1,*) '   XMAX,NX, XGLUE: ',SNGL(XMAX),NX,' ',SNGL(XGLUE)
	  WRITE(1,*) '   YMAX,NY, YGLUE: ',SNGL(YMAX),NY,' ',SNGL(YGLUE)
	END IF
	WRITE(1,*) ' IADINTXFI: ',IADINTXFI
	IF (IADINTXFI.NE.0) THEN
	  WRITE(1,*) '   SMAXPRE: ',SNGL(SMAXPRE)
	  WRITE(1,*) '   NELMAXPRE,NBACKMAX: ',NELMAXPRE,NBACKMAX
	  WRITE(1,*) '   IXTYP,NRUS: ',IXTYP,NRUS
	  WRITE(1,*) '   IRUS,JRUS: ',IRUS,JRUS
	  WRITE(1,*) '   ARGDIFLIM"1,0,2": ',
     $	  SNGL(ARGDIFLIM1),SNGL(ARGDIFLIM0),SNGL(ARGDIFLIM2)
	END IF
	WRITE(1,*)
	WRITE(1,*) ' NPARTS: ',NPARTS
	DO JP = 1, NPARTS
	  WRITE(1,*) '  PART JP=',JP,' :'
	  WRITE(1,*) '    ibessarr(jp),npuarr(jp) ',IBESSARR(JP),NPUARR(JP)
	  IF (NPUARR(JP).EQ.0) THEN
	    WRITE(1,*) '    mainmod,xyonenorm'
	    WRITE(1,*) '    ',MAINMOD(JP),' ',SNGL(XYONENORM(JP))
	  ELSE IF (NPUARR(JP).GT.0) THEN
	    WRITE(1,*) '    uphoern(1:npu)'
	    DO JL = 1, NPUARR(JP)
	      WRITE(1,*) '    ',UPHOERN(JL,JP)
	    END DO
	  ELSE
	    WRITE(1,*) '    aimend ',SNGL(AIMEND(JP))
	  END IF
	END DO
	IF (NPTYP.NE.0) THEN
	  WRITE(1,*) ' HMAXDH,HMAXRELWP: ',SNGL(HMAXDH),SNGL(HMAXRELWP)
	  IF (IADINTXFI.NE.0) THEN
	    WRITE(1,*) ' EPSH01Q,PREEPS: ',SNGL(EPSH01Q),SNGL(PREEPS)
	    IF (NPTYP.EQ.2) THEN
	      WRITE(1,*) ' HKVARPART,LUNKVAR,FAKTKVAR: ',
     $	      SNGL(HKVARPART),LUNKVAR,SNGL(FAKTKVAR)
	      IF (HKVARPART.GT.0.0D0) THEN
	        IF (LUNKVAR.LE.0) STOP 500
	        IF (FAKTKVAR.GT.0.0D0) THEN
	          WRITE(1,*) '   DDSUPMIN,DDSUPMAX: ',SNGL(DDSUPMIN),
     $	          SNGL(DDSUPMAX)
	          WRITE(1,*) '   DDSUPFAKTMAX: ',SNGL(DDSUPFAKTMAX)
	          WRITE(1,*) '   DERIVDRELMAX,DRELMAX: ',
     $	          SNGL(DERIVDRELMAX),SNGL(DRELMAX)
	          WRITE(1,*) '   DRELGOALFAKT, ARGDIFSTORE: ',
     $	          SNGL(DRELGOALFAKT),' ',SNGL(ARGDIFSTORE)
	        END IF
	      END IF
	    END IF
	  END IF
	END IF
	WRITE(1,*)
	WRITE(1,*) ' ICONTTYP: ',ICONTTYP
	IF (ICONTTYP.NE.0) THEN
	  WRITE(1,*) ' ITAIL(1),ITAIL(2) : ',ITAIL(1),ITAIL(2)
	  IF (ITAIL(1).NE.0) THEN
	    WRITE(1,*) '   IBESSTAIL(1),ANGRKTAIL(1): ',
     $	    IBESSTAIL(1),SNGL(ANGRKTAIL(1))
	    WRITE(1,*) '   UP0TAIL(1): ',UP0TAIL(1)
	  END IF
	  IF (ITAIL(2).NE.0) THEN
	    WRITE(1,*) '   IBESSTAIL(2),ANGRKTAIL(2): ',
     $	    IBESSTAIL(2),SNGL(ANGRKTAIL(2))
	    WRITE(1,*) '   UP0TAIL(2): ',UP0TAIL(2)
	  END IF
	  WRITE(1,*) ' ANPHPMINTERM: ',SNGL(ANPHPMINTERM)
	  IF (IADINTXFI.NE.0) WRITE(1,*) ' PREEPSTERM: ',SNGL(PREEPSTERM)
	  IF (ICONTTYP.EQ.1) THEN
	    WRITE(1,*) ' IRCONT,PREEPSCONT: ',IRCONT,SNGL(PREEPSCONT)
	  ELSE
	    WRITE(1,*) ' KMAXEU,PREEPSEU: ',KMAXEU,SNGL(PREEPSEU)
	  END IF
	END IF

	OPEN(UNIT=LUNIN,FILE=FILM,STATUS='OLD')
	ISLASK = IAZIMI
	IAZIMI = 0              ! for this initial "particular" call
	CALL RMODEL ( LUNIN,1, MHQ0 )
	IAZIMI = ISLASK         ! reset iazimi
	IF (LUNINSGN.GT.0) THEN
	  IF (NSOURCFW.NE.0) STOP 41553
	  IF (NSW.EQ.0) THEN
	    WRITE(*,*) 'must have "(interior) fluid" ("nss<nsw") for rpress'
	    STOP 453515
	  END IF
	  IF (NHYD.GT.0) THEN
	    IF ( IDHYD(1).LT.NSS .OR. IDHYD(NHYD).GT.NSW ) THEN
	      WRITE(*,*) 'receivers only in "(interior) fluid" for rpress'
	      STOP 453515
	    END IF
	  END IF
	ELSE
	  IF (NSOURCFW.LE.0) STOP 41554
	  IF (INTREST.NE.0) STOP 3415560 !impl "kfwtyp=4,nso sourc,fluid rec"
	END IF
	CLOSE(LUNIN)
	NHYDLI = NHYD*LI
	WRITE(1,*)
	WRITE(1,*) ' BRANCH POINT AT THE U ORIGIN FOR H0(1)'
	WRITE(1,*) '    the corresponding cut is along the negative rk axis'
	IF (ICONTTYP.EQ.0) THEN
	  ISLASK = 0
	ELSE IF (ICONTTYP.EQ.1) THEN
	  ISLASK = NHYDLI
	ELSE
	  ISLASK = 2*NHYD
	END IF
	IF ( ((IIFSMSOME.NE.0.AND.EPSFSUB.LE.0.0D0).OR.
     $	  IIMHQSOME.NE.0) .AND. NRUS0.NE.0 ) THEN
	  IF ( IRUS0.GT.NRUS0 .OR. JRUS0.LT.2 ) STOP 447
	  NHYDLICOP = NHYDLI
	  IADR00_TABNEW = 0     ! "tabnew(1:nhydlicop)" <--> cpressw(1:nhydli)
	  IADR00_VALTERM = IADR00_TABNEW + NHYDLICOP + ISLASK
	  IADR00_TAB = IADR00_VALTERM + NHYDLICOP
	  ICWORKLIM0 = IADR00_TAB + MIN0(NRUS0,JRUS0)*NHYDLICOP
	ELSE
	  ICWORKLIM0 = NHYDLI + ISLASK    !   cpressw ( 1 : nhydli+islask )
	END IF
	WRITE(1,*)
	WRITE(1,*) ' ICWORKLIM0 = ',ICWORKLIM0
	CLOSE(1)
	IF (ICWORKLIM0.GT.icworkdim) STOP 451

	IF (FILG.NE.' ') THEN
	  IF (LUNINSGN.GT.0) THEN
	    WRITE(23,*) -0.00001,' ! from RPRESS; FILO: ',FILO
	  ELSE
	    WRITE(23,*) -0.00001,' ! from RPRESSFW; FILO: ',FILO
	  END IF
	  WRITE(23,*) NHYD,LI,'    ! nhyd,li (i.e. z=depth,x=range)'
	  WRITE(23,*) (SNGL(ZD(IDHYD(IH))),IH=1,NHYD)
	  WRITE(23,*) (SNGL(OR(IR)),IR=1,LI)
	  CLOSE(23)
	END IF

	MHQ(1) = 1

	IF (NRUS0.NE.0) THEN

C	  update the "Bulirsch Folge"
	  MHQ(2) = 2
	  MHQ(3) = 3
	  IXROW = 4
	  DO WHILE (IXROW.LE.NRUS0)
	    MHQ(IXROW) = 2 * MHQ(IXROW-2)
	    IXROW = IXROW + 2
	  END DO
	  IXROW = 5
	  DO WHILE (IXROW.LE.NRUS0)
	    MHQ(IXROW) = 2 * MHQ(IXROW-2)
	    IXROW = IXROW + 2
	  END DO

C	  computation of "hqi(,)"
	  IF ( IIFSMSOME.EQ.0 .OR. EPSFSUB.GT.0.0D0 ) THEN
	    IF (IIMHQSOME.EQ.0) THEN
	      IF (WTXT) WRITE(LUNWRITE,*) ' %%% outer extrapolation not needed'
	      NRUS0 = 0
	    ELSE
C	      now "I(h)" series with the even powers of h ("ODE trapezoidal")
	      DO I=2,NRUS0            ! i=iminx,nrus
	        JCOLX = MIN0(I,JRUS0)
	        DO J=2,JCOLX
	          HQI(J,I)  =  ( DBLE(MHQ(I)) / MHQ(I-J+1) ) ** 2 
C	                                                  ! "i"/"i-j+1"
	        END DO
	      END DO
	    END IF
	  ELSE IF (IABS(IXTYP0).NE.2) THEN
	    JCOLX = MIN0(NRUS0,JRUS0)
	    DO M=1,JCOLX-1
	      IF (IIMHQSOME.EQ.0) THEN
C	        now "I(h)" series with h**5,h**6,... (Dormand-Prince global 
C	                                              error)
	        IEXPX = 4+M
	      ELSE
C	        now "I(h)" series with h**2,h**4,h**5,h**6,... (combination)
	        IF (M.EQ.1) THEN
	          IEXPX = 2
	        ELSE
	          IEXPX = 2+M
	        END IF
	      END IF
	      DO I=1,NRUS0
	        ACOEF(M,I) = ( DBLE(MHQ(1)) / MHQ(I) ) ** IEXPX
	      END DO
	    END DO
	    DO J=2,JCOLX
	      DO I=NRUS0,J,-1       ! iminx=2
	        SLASK = ACOEF(J-1,I-1) / ACOEF(J-1,I)
	        HQI(J,I) = SLASK
	        SLASKD = 1 / (SLASK-1)
	        DO M=J,JCOLX-1
	          ACOEF(M,I) = SLASKD * ( SLASK*ACOEF(M,I) - ACOEF(M,I-1) )
	        END DO
	      END DO
	    END DO
	  ELSE
C	    now "I(h)" series with the even powers of h, disregarding 
C	    a priori information on "gaps" in the asymptotic series
	    DO I=2,NRUS0            ! i=iminx,nrus
	      JCOLX = MIN0(I,JRUS0)
	      DO J=2,JCOLX
	        HQI(J,I)  =  ( DBLE(MHQ(I)) / MHQ(I-J+1) ) ** 2 
C	                                                  ! "i"/"i-j+1"
	      END DO
	    END DO
	  END IF

	END IF

	DO JF = 1, NFR

	  IF (WTXT) WRITE(LUNWRITE,*) 
	  IF (WTXT) WRITE(LUNWRITE,*) 'jf,nfr:',JF,NFR
	  FREQ = OFREQ(JF)
          OMEG = DCMPLX ( 2*PI*FREQ, ALA )

	  IF (LX.GT.1) THEN

	    IF (FREQDOM.GT.0.0D0) THEN
	      CR9 = 0.5D0 * ( (1.0D0,0.0D0) + CDCOS(PI*OMEG/OMN) )
C	      this is to simulate an explosive point source with diagonal 
C	        moment-tensor with diagonal elements " s(t) " "(10**15 N)*(km)"
C	        where s(t) is the pulse-like signal corresponding to the 
C	        spectrum given by "cr9 for jf=1,..,nfr", tmo, ala 
C	      hence the volume-force will be 
C	        " - s(t) * grad(Diracfunction) "  "(10**15 N)/(km**3)"
C	      for the case "iazimi=0" the pressure in a punctured disc around 
C	        the source point will consequently be (approximatively) 
C	        " ( 1/(4*pi*r) ) * f''(t-Re(r*ua)) " "(10**15 N)/(km**2)" 
C	        where   f(t)  is  "s(t)" convolved with the function 
C	           " Im(-ua/pi) * ( t**2 + (r*Im(ua))**2 )**(-1) * 
C	             ( r*Re(ua**2) + 2*Re(ua)*t ) " 
C	           ( ua*ua * Dirac(t)  when Im(ua) = 0.0 )
C	        and r is spherical distance (!) in km, omeg in Hz, ua in s/km, 
C	        and t in s
C	      %% note that "(10**15 N)" = "(10**12 kg)*(km/(s**2))" %%
C	      %% Im(ua) >= 0.0 ( Re(omeg) > 0.0 within this frequency loop ) %%
C	      %% for the case "iazimi=0" the near-field response is not causal 
C	         when Im(ua) > 0.0 
C	         ( causality requires a frequency dependence of Q, cf. 
C	         Aki ch. 5, as implemented for the case "iazimi nonzero" ) %%
	    ELSE
	      CR9 = CPLOK(JF)
	    END IF

	    IF (CR9.EQ.(0.0D0,0.0D0)) THEN
	      DO IHR = 1, NHYDLI
	        CPRESSW(IHR) = (0.0D0,0.0D0)
	      END DO
	      WRITE(21) (CPRESSW(IHR),IHR=1,NHYDLI)
	      GOTO 99
	    END IF

	  END IF

	  IF (IAZIMI.NE.0) THEN
	    CSLASK = ALFA1AZIM * OMEG
	    FREAZIM = 2 * AZIMILN(CSLASK) / PI
	    FIMAZIM = CSLASK / (1+CSLASK)
	  END IF	    
	  MHQNOW = FREQ/FREQMHQ0
	  MHQNOW = MHQ0 * JMAX0(1,MHQNOW)

	  MHQEFF = MHQNOW
	  CALL RMODEL ( 0,0, MHQEFF )

	  FSOURC = UA2SOURC / (2*PI*RHOSOURC)
C	  the preceding factor is for the representation with a discontinuity 
C	  in the second component of the "Aki (!) displacement-stress 
C	  vector (Aki & Richards 1980, 7.4.2) but with omeg* for the 
C	  first two components" at the source level, the factor 1/(2*pi) 
C	  and the - from the Aki vector basis-functions (rel. the Bessel 
C	  functions) are included here

C	  consider our point source located in a homogeneous acoustic (no 
C	  shear) whole-space with the same medium parameters as at the source 
C	  depth in our original medium, here our point source will require a 
C	  certain power (which is in general different from the power 
C	  required in our original medium)
	  AM0 = DSQRT ( (1.0D0/RHOSOURC) * (1.5D0*DREAL(CDSQRT(
     $	  UA2SOURC))) ) * 2.25D0*CDABS(UA2SOURC)
C	  am0 is the volume-force strength factor for
C	   " - exp(-i*omeg*t) * grad(Diracfunction) "  "(10**15 N)*km"
C	  to provide a point source of the same frequency and the specified 
C	  power but in a homogeneous acoustic (no shear) whole-space with 
C	     P-velocity = 1.5 km/s and no anelastic attenuation (infinite Q)
C	     density =  1.0 "10**12 kg"/(km**3), or 1.0 kg/(dm**3)
	  PM0FAKT = (1/1.5D0)**2 * AM0 / (4*PI*1.0D0)
C	  for real omeg,  pm0 = pm0fakt * dreal(omeg)**2  will be the corres-
C	  ponding peak pressure, in (10**15 N)/(km**2), at a distance of 1 km
	  PM0 = PM0FAKT * DREAL(OMEG)**2  
	  PM01M = 1000*PM0             ! at 1 m rather than 1 km

	  DO I = 1, 6
	    NFCALLS(I) = 0
	  END DO
	  DO IR=1,LI
	    IF (OR(IR).GT.0.0D0) THEN
	      PM0GOR(IR) = PM0/OR(IR)
	    ELSE
	      PM0GOR(IR) = PM0                  ! somewhat arbitrarily
	    END IF
	  END DO
	  IF (LX.GT.1) THEN
	    IF (PREALA.NE.1.0D0) THEN
	      DO IR = 1, LI
                T0 = TMIN + OR(IR)*URED - TMO
	        PM0GOR(IR) = PREALA**(T0/DAUER) * PM0GOR(IR)
	      END DO
	    END IF
	  END IF
	  DO IR=1,LI
	    EPS0(IR) = PREEPS0 * PM0GOR(IR)
	  END DO
	  IF (IIFSMSOME.NE.0) THEN
	    IF (EPSFSUB.GT.0.0D0) THEN
	      HSTEPFSUB = HSTEPFSUBPRE/FREQ
	    ELSE
	      HSTEPFSUB0 = DMIN1 ( HSTEPFSUBPRE/FREQ, FSMDZMIN )  ! to avoid 
C	      extrapolating on "equal hstepfsub"
	    END IF
	  END IF
	  IF (NPTYP.NE.0) THEN
	    IF (OR(LI).GT.0.0D0) THEN
	      HMAX = DMIN1 ( HMAXDH , HMAXRELWP/(FREQ*OR(LI)) )
	    ELSE
	      HMAX = HMAXDH
	    END IF
	    IF (NPTYP.EQ.2) HKVAR = HKVARPART*HMAX
	  END IF

	  DO JP = 1, NPARTS
	    IF (NPUARR(JP).EQ.0) THEN
	      FILMOD(JP) = MAINMOD(JP)
	      IPOS = IPOSMOD(JP)
	      ISLASK = FREQ
	      CALL PUTINTP (ISLASK,IPOS,FILMOD(JP))
	      IF (ISLASK.NE.FREQ) THEN
	        IF (IPOS.GT.32) STOP 460
	        FILMOD(JP)(IPOS:IPOS) = 'X'
	        IPOS = IPOS+1
	        CALL PUTINTP (JF,IPOS,FILMOD(JP))
	      END IF
	    END IF
	  END DO

	  IXROW = 0

	  DO WHILE (.TRUE.)

	    IXROW = IXROW + 1

	    IF (IIFSMSOME.NE.0) THEN
	      IF (EPSFSUB.LE.0.0D0) HSTEPFSUB = HSTEPFSUB0 / MHQ(IXROW)
	    END IF
	    IF (IXROW.GT.1) THEN
	      MHQEFF = MHQNOW * MHQ(IXROW)
	      CALL RMODEL ( 0,0, MHQEFF )
	    END IF

	    DO I = 1, 6
	      NFCALLSOLD(I) = NFCALLS(I)
	    END DO

	    CALL RPRESSINT ( FSOURC,PM0GOR, IUROT, npudim, 
     $	    IFILONPRE, INTREST, NPARTS,IBESSARR,NPUARR,
     $	    FILMOD,XYONENORM,UPHOERN,AIMEND, UPHOERNLOK, 
     $	    NHYDLI, ICWORKLIM0,ICWORKDIM,CPRESSW )

	    IF (NRUS0.EQ.0) THEN
	      POSTEPS0 = -1.0D0
	      GOTO 101
	    END IF

	    IF (IXROW.EQ.1) THEN

	      DO IHR=1,NHYDLI
	        CPRESSW(IA0_TAB(IHR,1)) = CPRESSW(IHR)
	      END DO

	    ELSE

C	      do ihr=1,nhydli
C	        cpressw(ia0_tabnew(ihr)) = cpressw(ihr)
C	      end do                      ! done since  iadr00_tabnew = 0
	      JCOLX = MIN0(IXROW,JRUS0)
	      IF (IXROW.GE.IRUS0) THEN
	        OKSUB = .TRUE.
	      ELSE
	        OKSUB = .FALSE.
	      END IF
	      CALL XPOLBUL ( NHYD,LI, IXTYP0, 1.0D0,EPS0, 
     $	      HQI(1,IXROW),FIHQI,QHQI, 2,JCOLX, 
     $	      IA0_TABNEW,IA0_VALTERM,IA0_TAB,IA0_YTAB, 
     $	      CPRESSW,OKSUB )
	      IF ( WTXT .AND. I0LOG.NE.0 ) THEN
	        DO IOUT = 100, 102
	          OPEN(IOUT,STATUS='UNKNOWN')
	          WRITE(IOUT,*) ' extrapolation log fr rpress(fw); iout=',IOUT
	          WRITE(IOUT,*) '       "100" for " cdabs(valterm) "'
	          WRITE(IOUT,*) '       "101" for " delmin / (pm0/or() )"'
	          WRITE(IOUT,*) '       "102" for " jdel "'
	          WRITE(IOUT,*) ' MHQEFF, FREQ: ',MHQEFF,FREQ
	          WRITE(IOUT,*) ' FILO ',FILO
	          WRITE(IOUT,*)
	          WRITE(IOUT,*) ' BELOW FOR IXROW = ',IXROW
	          WRITE(IOUT,*)
	        END DO
	        WRITE(100,*) ( SNGL ( PM0GOR(IR) ) , IR=1,LI )
	        WRITE(100,*)
	        DO IH = 1, NHYD
	          WRITE(100,*) ( SNGL ( CDABS(CPRESSW(IA0_VALTERM(
     $	          IHR))) ) , IHR=(IH-1)*LI+1,IH*LI )
	        END DO
	        DO IH = 1, NHYD
	          WRITE(101,*) ( SNGL ( DREAL(CPRESSW((IH-1)*LI+IR
     $	          ))/(PM0GOR(IR)) ) , IR=1,LI )      ! "tabnew equivalence", 
C	                                              " delmin / (pm0/or()) "
	        END DO
	        DO IH = 1, NHYD
	          WRITE(102,63) ( IDNINT ( DIMAG(CPRESSW(IHR)) ) , 
     $	          IHR=(IH-1)*LI+1,IH*LI )    ! "tabnew equivalence", "jdel"
 63	          FORMAT(1X,25I3)
	        END DO
	        DO IOUT = 100, 102
	          WRITE(IOUT,*)
	          WRITE(IOUT,133) (NFCALLS(I)-NFCALLSOLD(I),I=1,6),
     $	          ' ! nfcalls now'
 133	          FORMAT(1X,6I10,A)
	          WRITE(IOUT,133) (IDNINT(DBLE(NFCALLS(I))/IXROW),
     $	          I=1,6),' ! nfcalls ave'
	          CLOSE(IOUT)
	        END DO
	      END IF
	      IF ( OKSUB .OR. IXROW.GE.NRUS0 ) THEN
	        POSTEPS0 = 0.0D0
	        DO IR = 1, LI
	          SLASK = PM0GOR(IR)
	          IHR = IR
	          DO IH = 1, NHYD
	            POSTEPS0 = DMAX1 ( DREAL(CPRESSW(IHR))/SLASK, 
     $	            POSTEPS0 ) ! "tabnew equivalence", " delmin / (pm0/or()) "
	            IHR = IHR + LI
	          END DO
	        END DO
	        DO IHR=1,NHYDLI
	          CPRESSW(IHR) = CPRESSW(IA0_VALTERM(IHR))
	        END DO
	        GOTO 101
	      END IF

	    END IF

	  END DO

 101	  PREFREL1 = CDABS(CPRESSW(1)) / (PM0GOR(1))
	  PREFREL2 = PREFREL1
	  DO IR = 1, LI
	    SLASK = PM0GOR(IR)
	    IHR = IR
	    DO IH = 1, NHYD
	      SLASKD = CDABS(CPRESSW(IHR)) / SLASK
	      IF (SLASKD.LT.PREFREL1) PREFREL1 = SLASKD
	      IF (SLASKD.GT.PREFREL2) PREFREL2 = SLASKD
	      IHR = IHR + LI
	    END DO
	  END DO

	  IF (MAINAM.NE.' ') THEN
	    FILAM = MAINAM
	    IPOS = IPOSAM
	    ISLASK = FREQ
	    CALL PUTINTP (ISLASK,IPOS,FILAM)
	    IF (ISLASK.NE.FREQ) THEN
	      IF (IPOS.GT.32) STOP 460
	      FILAM(IPOS:IPOS) = 'X'
	      IPOS = IPOS+1
	      CALL PUTINTP (JF,IPOS,FILAM)
	    END IF
	    OPEN(UNIT=24,FILE=FILAM,STATUS='UNKNOWN')
	    WRITE(24,*) ' 6 ',SNGL(FREQ),PM0,IDB,'  ! nskip,freq,pm0,idb'
	    IF (LUNINSGN.GT.0) THEN
	      WRITE(24,*) ' FILAM fr RPRESS  MHQEFF,FILM=',MHQEFF,' ',FILM
	    ELSE
	      WRITE(24,*) ' FILAM fr RPRESSFW  MHQEFF,FILM=',MHQEFF,' ',FILM
	    END IF
	    IF (IAZIMI.EQ.0) THEN
	      WRITE(24,*) ' KFWTYP, IAZIMI: ',KFWTYP, IAZIMI
	    ELSE
	      WRITE(24,*) ' KFWTYP, IAZIMI,ALFA1AZIM: ',KFWTYP, 
     $	      IAZIMI,SNGL(ALFA1AZIM)
	    END IF
	    WRITE(24,*) ' FILO: ',FILO
	    WRITE(24,*) ' FILG: ',FILG
	    WRITE(24,*) ' Out MAINAM: ',MAINAM,'       IDB=',IDB
	    WRITE(24,*) ' Out MAINPH: ',MAINPH,'       IDB=',IDB
	    WRITE(24,*) NHYD,LI,' ! nhyd,li ; freq = ',SNGL(FREQ)
	    IF (IDB.GT.0) THEN
	      DO IH = 1, NHYD
	        WRITE(24,*) ( SNGL ( 20.0 * DLOG10 ( CDABS(CPRESSW(IHR
     $	        )) / PM01M ) ) , IHR=(IH-1)*LI+1,IH*LI )
	      END DO
	    ELSE IF (IDB.EQ.0) THEN
	      DO IH = 1, NHYD
	        WRITE(24,*) ( SNGL ( CDABS(CPRESSW(IHR)) / PM01M ) , 
     $	        IHR=(IH-1)*LI+1,IH*LI )
	      END DO
	    ELSE IF (IDB.EQ.-1) THEN
	      DO IH = 1, NHYD
	        WRITE(24,*) ( SNGL ( CDABS(CPRESSW(IHR)) ) , 
     $	        IHR=(IH-1)*LI+1,IH*LI )
	      END DO
	    ELSE IF (IDB.EQ.-2) THEN
	      DO IH = 1, NHYD
	        WRITE(24,*) (CDABS(CPRESSW(IHR)),IHR=(IH-1)*LI+1,IH*LI)
	      END DO
	    ELSE
	      SLASK = CDABS(RHOSOURC/UA2SOURC)
	      DO IH = 1, NHYD
	        WRITE(24,*) (SLASK*CDABS(CPRESSW(IHR)),IHR=(IH-1)*LI+1,IH*LI)
	      END DO
	    END IF
	  END IF

	  IF (MAINPH.NE.' ') THEN
	    FILPH = MAINPH
	    IPOS = IPOSPH
	    ISLASK = FREQ
	    CALL PUTINTP (ISLASK,IPOS,FILPH)
	    IF (ISLASK.NE.FREQ) THEN
	      IF (IPOS.GT.32) STOP 460
	      FILPH(IPOS:IPOS) = 'X'
	      IPOS = IPOS+1
	      CALL PUTINTP (JF,IPOS,FILPH)
	    END IF
	    OPEN(UNIT=25,FILE=FILPH,STATUS='UNKNOWN')
	    WRITE(25,*) ' 6 ',SNGL(FREQ),PM0,IDB,'  ! nskip,freq,pm0,idb'
	    IF (LUNINSGN.GT.0) THEN
	      WRITE(25,*) ' FILPH fr RPRESS  MHQEFF,FILM=',MHQEFF,' ',FILM
	    ELSE
	      WRITE(25,*) ' FILPH fr RPRESSFW  MHQEFF,FILM=',MHQEFF,' ',FILM
	    END IF
	    IF (IAZIMI.EQ.0) THEN
	      WRITE(25,*) ' KFWTYP, IAZIMI: ',KFWTYP, IAZIMI
	    ELSE
	      WRITE(25,*) ' KFWTYP, IAZIMI,ALFA1AZIM: ',KFWTYP, 
     $	      IAZIMI,SNGL(ALFA1AZIM)
	    END IF
	    WRITE(25,*) ' FILO: ',FILO
	    WRITE(25,*) ' FILG: ',FILG
	    WRITE(25,*) ' Out MAINAM: ',MAINAM,'       IDB=',IDB
	    WRITE(25,*) ' Out MAINPH: ',MAINPH,'       IDB=',IDB
	    WRITE(25,*) NHYD,LI,' ! nhyd,li ; freq = ',SNGL(FREQ)
	    IF (IDB.GE.-1) THEN
	      DO IH = 1, NHYD
	        WRITE(25,*) ( SNGL ( DATAN2DQ(DIMAG(CPRESSW(IHR)),
     $	        DREAL(CPRESSW(IHR))) ) , IHR=(IH-1)*LI+1,IH*LI )
	      END DO
	    ELSE IF (IDB.EQ.-2) THEN
	      DO IH = 1, NHYD
	        WRITE(25,*) ( DATAN2DQ(DIMAG(CPRESSW(IHR)),
     $	        DREAL(CPRESSW(IHR))) , IHR=(IH-1)*LI+1,IH*LI )
	      END DO
	    ELSE
	      CSLASK = RHOSOURC/UA2SOURC
	      DO IH = 1, NHYD
	        WRITE(25,*) ( DATAN2DQ(DIMAG(CSLASK*CPRESSW(IHR)),
     $	        DREAL(CSLASK*CPRESSW(IHR))) , IHR=(IH-1)*LI+1,IH*LI )
	      END DO
	    END IF
	  END IF

	  DO ISLASK = 1, 3
	    IF (ISLASK.EQ.1) THEN
	      IF (MAINAM.NE.' ') THEN
	        IOUT = 24
	      ELSE
	        IOUT = 0
	      END IF
	    ELSE IF (ISLASK.EQ.2) THEN
	      IF (MAINPH.NE.' ') THEN
	        IOUT = 25
	      ELSE
	        IOUT = 0
	      END IF
	    ELSE
	      IF (WTXT) THEN
	        IOUT = LUNWRITE
	      ELSE
		IOUT = 0
	      END IF
	    END IF
	    IF (IOUT.GT.0) THEN
	      WRITE(IOUT,*) 
	      WRITE(IOUT,133) (NFCALLS(I)-NFCALLSOLD(I),I=1,6),' ! nfcalls now'
	      WRITE(IOUT,133) (IDNINT(DBLE(NFCALLS(I))/IXROW),
     $	      I=1,6),' ! nfcalls ave'
	      WRITE(IOUT,*) 
	      WRITE(IOUT,*) IXROW,NRUS0,'   ! ixrow,nrus0'
	      IF ( NRUS0.NE.0 .AND. IXROW.GT.NRUS0 ) WRITE(IOUT,
     $	      *) '   %%%  WARNING  ACCURACY TEST NOT FULFILLED  %%%'
	      WRITE(IOUT,*) 
	      WRITE(IOUT,*) SNGL(AM0),SNGL(PM0),'   ! am0,pm0'
	      WRITE(IOUT,*)
	      WRITE(IOUT,*) SNGL(PREFREL1),SNGL(PREFREL2),
     $	      '  ! prefrel1,prefrel2'
	      IF ( IADINTXFI.EQ.0 .OR. NPTYP.EQ.0 ) THEN
	        WRITE(IOUT,*) -1.0,-1.0,'  ! supepsmagn,posteps'
	      ELSE IF (PREFREL1.GT.0.0D0) THEN
	        POSTEPS = SUPEPSMAGN*PREEPS / PREFREL1
	        WRITE(IOUT,*) SNGL(SUPEPSMAGN),SNGL(POSTEPS),
     $	        '  ! supepsmagn,posteps'
	      ELSE
	        WRITE(IOUT,*) SNGL(SUPEPSMAGN),-1.0,'  ! supepsmagn,posteps'
	      END IF
	      IF ( IADINTXFI.EQ.0 .OR. ICONTTYP.EQ.0 ) THEN
	        WRITE(IOUT,*) -1.0,'  ! postepsterm'
	      ELSE IF (PREFREL1.GT.0.0D0) THEN
	        POSTEPSTERM = PREEPSTERM / PREFREL1
	        WRITE(IOUT,*) SNGL(POSTEPSTERM),'  ! postepsterm'
	      ELSE
	        WRITE(IOUT,*) -1.0,'  ! postepsterm'
	      END IF
	      IF (ICONTTYP.EQ.0) THEN
	        WRITE(IOUT,*) -1.0,'  ! posteps"cont,eu"'
	      ELSE IF (ICONTTYP.EQ.1) THEN
	        IF (PREFREL1.GT.0.0D0) THEN
	          POSTEPSCONT = PREEPSCONT / PREFREL1
	          WRITE(IOUT,*) SNGL(POSTEPSCONT),'  ! postepscont'
	        ELSE
	          WRITE(IOUT,*) -1.0,'  ! postepscont'
	        END IF
	      ELSE
	        IF (PREFREL1.GT.0.0D0) THEN
	          POSTEPSEU = PREEPSEU / PREFREL1
	          WRITE(IOUT,*) SNGL(POSTEPSEU),'  ! postepseu'
	        ELSE
	          WRITE(IOUT,*) -1.0,'  ! postepseu'
	        END IF
	      END IF
	      WRITE(IOUT,*) 
	      WRITE(IOUT,*) IIFSMSOME,SNGL(FSMDZMIN),IIMHQSOME,
     $	      ' ! iifsmsome,fsmdzmin,iimhqsome'
	      WRITE(IOUT,*) SNGL(POSTEPS0),IXROW,' ! posteps0,ixrow'
	      IF (ISLASK.EQ.1) THEN
	        CLOSE(24)
	      ELSE IF (ISLASK.EQ.2) THEN
	        CLOSE(25)
	      END IF
	    END IF
	  END DO

	  CALL UTSPEC (KFWTYP, JF,FREQ,PM0,IDB, 
     $    MHQEFF,FILM, IAZIMI,ALFA1AZIM, FILO,FILG, RHOSOURC,UA2SOURC, 
     $    NHYD,LI, OR, CPRESSW, 
     $    MAINAM,IPOSAM, MAINPH,IPOSPH)

	  IF (LX.GT.1) THEN

	    IF ( NRUS0.NE.0 .AND. IXROW.GT.NRUS0 ) STOP 713

            DO IR = 1, LI
              T0 = TMIN + OR(IR)*URED - TMO
              T1E = T0*OMEG
	      T1E = CDEXP ( DCMPLX(DIMAG(T1E),-DREAL(T1E)) ) !cdexp(-i*t0*omeg)
	      CR9T1E(IR) = CR9 * T1E
	    END DO
	    IHR = 0
	    DO IH = 1, NHYD
	      DO IR = 1, LI
	        IHR = IHR+1
	        CPRESSW(IHR) = CR9T1E(IR) * CPRESSW(IHR)
	      END DO
	    END DO
	    WRITE(21) (CPRESSW(IHR),IHR=1,NHYDLI)

	  END IF

 99	  CONTINUE

	END DO

	IF (LX.GT.1) THEN

	  IOUT = 1
	  DO ID = 1, NOL
	    IF (IFSMTYP(ID).GE.-3) IOUT = 0
	  END DO
	  IF (IOUT.EQ.1) THEN
	    DO ID = 1, NOL
	      IF (UB2(ID).EQ.(0.0D0,0.0D0)) UB2(ID) = (0.0D0,-1.0D0)
	    END DO         ! for the output only
	  END IF

	  DO IH = 1, NHYD
	    FILTR = MAINTR
	    IPOS = IPOSTR
	    CALL PUTINTP (IH,IPOS,FILTR)
	    OPEN(UNIT=22,FILE=FILTR,FORM='UNFORMATTED',STATUS='UNKNOWN')
	    WRITE(22) LI,LX,DDT,5,SNGL(TMIN),SNGL(DAUER),KFWTYP
	    IF (RMIN.GE.0.0D0) THEN
              WRITE(22) NSS,NSO,NSW,NOL,IDHYD(IH),SNGL(RMIN),
     $	      SNGL(RINT),SNGL(URED)
	    ELSE
              WRITE(22) NSS,NSO,NSW,NOL,IDHYD(IH),-1.0,-1.0,SNGL(URED)
	    END IF
            WRITE(22) (SNGL(DZ(I)),I=1,NOL)
	    IF (IOUT.EQ.1) THEN
              WRITE(22) (SNGL ( 1 / DREAL(CDSQRT(UA2(I))) ), I=1,NOL)
	      WRITE(22) (AMAX1( SNGL ( 1 / DREAL(CDSQRT12(UB2(I))) ), 
     $	      0.0 ), I=1,NOL)
              WRITE(22) (SNGL(RHO(I)),I=1,NOL)
	    ELSE
              WRITE(22) (0.0, I=1,NOL)
              WRITE(22) (0.0, I=1,NOL)
              WRITE(22) (0.0, I=1,NOL)
	    END IF
	    IF (RMIN.LT.0.0D0) WRITE(22) (SNGL(OR(I)),I=1,LI)
	    REWIND 21
	    DO JF = 1, NFR
	      J1 = NFR0-JF+1
	      J2 = LX-NFR+JF
	      READ(21) (CPRESSW(IHR),IHR=1,IH*LI)
	      IHR = (IH-1)*LI
	      DO IR = 1, LI
	        IHR = IHR+1
	        CSLASK = CPRESSW(IHR)
	        CP(J1,IR) = DCONJG(CSLASK)     ! the "conjg" is put here and 
	        CP(J2,IR) = CSLASK   ! not here, "exp(-i*omeg*t)" but "cfftal"
	      END DO
	    END DO
	    J1 = NFR0+1
	    J2 = LX-NFR
	    DO IR = 1, LI
	      CP(1,IR) = (0.0,0.0) !zero pressure in fluid for r>0 when omeg=0
	      DO JF = J1, J2
	        CP(JF,IR) = (0.0,0.0)
	      END DO
	      CALL CFFTAL (LX,DDF,CP(1,IR), 1, 1, 0, 
     $	      ISLASK,SLASKD,CPLOK(1))    ! ok here because of "conjg" as above
	      DO I = 1, LX
	        CPLOK(I) = EX(I)*CPLOK(I)
	      END DO
	      WRITE(22) (REAL(CPLOK(I)),I=1,LX)
	    END DO
	    CLOSE(22)
	  END DO

	  CLOSE(21)

	  IF (IOUT.EQ.1) THEN
	    DO ID = 1,NOL
	      IF (UB2(ID).EQ.(0.0D0,-1.0D0)) UB2(ID) = (0.0D0,0.0D0)
	    END DO
	  END IF

	END IF

	RETURN
	END






	REAL*8 FUNCTION DATAN2DQ (Y,X)
	IMPLICIT NONE
	REAL*8 X,Y
	IF ( X.NE.0.0D0 .OR. Y.NE.0.0D0 ) THEN
	  DATAN2DQ = DATAN2D (Y,X)
	ELSE
	  DATAN2DQ = 9999.9D0
	END IF
	RETURN
	END


 

	INTEGER FUNCTION IA0_TABNEW (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR00.INC'
	IA0_TABNEW = IADR00_TABNEW + IHR
	RETURN
	END

	INTEGER FUNCTION IA0_VALTERM (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR00.INC'
	IA0_VALTERM = IADR00_VALTERM + IHR
	RETURN
	END

	INTEGER FUNCTION IA0_TAB (IHR,J)
	IMPLICIT NONE
	INTEGER IHR,J
	INCLUDE 'IADR00.INC'
	IA0_TAB = IADR00_TAB + (J-1)*NHYDLICOP + IHR
	RETURN
	END

	INTEGER FUNCTION IA0_YTAB (IHR,J)
	IMPLICIT NONE
	INTEGER IHR,J
	INCLUDE 'IADR00.INC'
	IA0_YTAB = IADR00_YTAB + (J-1)*NHYDLICOP + IHR
	RETURN
	END
