
C	to compute modes as functions of depth

C	Sven Ivansson 1992
C	option "nsourcfw<-1" for "beam source array" added may 2020



	PROGRAM RMODPFW

	USE INMODELa
	USE FWYSa
	USE FWPARTa
	USE FWDISa
	USE FWYTa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C       nol <= noldim
	INTEGER NFRDIM, LLDIM
	PARAMETER ( NFRDIM = 4400 , LLDIM = 3500 )
C       nfr <= nfrdim  ,  ll <= lldim

	INCLUDE '../AIRYP.INC'
	INCLUDE '../INPARAM.INC'
	INCLUDE '../DPCFSUB.INC'
	INCLUDE '../INSPEC.INC'
	INCLUDE '../SOURCFW.INC'
	INCLUDE 'FWPARAM.INC'
	INCLUDE 'FWINIT.INC'
	INCLUDE '../single/BEAMARR.INC'  ! to handle beam source array
C	initialization using "common /fwinit/" with "overwrite" (OBS!!!)

	INTEGER*2 N1,N234,N789, NHYD2, I2
	INTEGER MHQ0,MHQNOW,MHQEFF, NFR0, JF,NFR,NFR1,NFR2, IUROT, IUTTYP, 
     $	INOLX, LLMAX, NOLORIG,IFSMNOLORIG,NSWORIG,NSOORIG, 
     $	LL, IPOSMODU,IPOSMODP,IPOS, IH, ID, I, IAB, II,IIA,IIB,IIW, 
     $	IDMAX,ICMAX, ISLASK, ISLASK1,ISLASK2,ISLASK3
	REAL*8 FREQ,FREQMHQ0, HSTEPFSUBPRE, FREQMIN,FREQINT, 
     $  ALADAUER,PREALA,ALA, FREQDOM,HALFBAND, 
     $  DLNDO, DLNDU, DCHINUREMAX,DCHINUIMMAX, SLASK, SLASK1,SLASK2, PI, 
     &	ZMARR,HZARR,XARR,SLWARR,WEIGHTSUM, DSIND, 
     $  OFREQ(nfrdim)
	COMPLEX*16 AZIMILN, UP,UP2, U,U2, DERIVDO, DERIVDU, DUDO, CSLASK, 
     $	CSUM, CSLASK0, QFAKT, 
     $	UPARR(lldim)
	LOGICAL JFCYCLE, FNORMAL, INKLUD
	COMPLEX CESLASK
	CHARACTER*42 FILM, MAINMODU,MAINMODP, FILMODU,FILMODP

	INTEGER,    DIMENSION(:), ALLOCATABLE :: INDABLOK
	REAL*8,     DIMENSION(:,:), ALLOCATABLE :: RHOLOK
	COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: UA2LOK,UB2LOK
	COMPLEX*16, DIMENSION(:), ALLOCATABLE :: UGRARR

	DATA PI /3.141592653589793D0/

 	UPREF = (0.0D0,0.0D0)
	IDET = 2        ! different from -1
	IRKSPEC = 0

	WRITE(*,*)
	WRITE(*,*)
	WRITE(*,*) ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
	WRITE(*,*) ' NOTE! for the strain-energy densities, one ',
     $	'would like to use '
	WRITE(*,*) ' "Lame-paremeters with relaxation neglected", ',
     $	'cf. Hudson ch.9; '
	WRITE(*,*) ' however, this "neglect" is not done here since ',
     $	'the real parts '
	WRITE(*,*) ' of the "relaxation contributions" are not ',
     $	'provided among '
	WRITE(*,*) ' the indata (moreover, they are in general "small")'
	WRITE(*,*) ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
	WRITE(*,*)

	WRITE(*,*)
	WRITE(*,*) 'Give infile FILM:'
	FILM = ' '
	DO WHILE (FILM.EQ.' ') 
	  READ(*,'(A)') FILM
	END DO
	CALL BLANKTOEND (FILM)
	WRITE(*,*) 'Give MHQ0, FREQMHQ0 (Hz), IAZIMI :'
	READ(*,*) MHQ0,FREQMHQ0, IAZIMI
	IF (MHQ0.GE.0) THEN
	  AI4E = .TRUE.
	ELSE
	  MHQ0 = -MHQ0
	  AI4E = .FALSE.
	END IF
	IF (IAZIMI.NE.0) THEN
	  WRITE(*,*) 'Give ALFA1AZIM:'
	  READ(*,*) ALFA1AZIM
	END IF
	WRITE(*,*) 'Give ICUTTYP ( 1="hyperbolic for up" 2="vertical for up" '
	WRITE(*,*) '3="mixture with aimcut" ) :'
	WRITE(*,*) '( "-5 or -15 or -45 or -60" for Muir )'
	WRITE(*,*) '( irrelevant for rigid or free boundaries )'
	WRITE(*,*) '( rk=omeg*u, +- irrelevant if "preala"=1.0 )'
	READ(*,*) ICUTTYP
	IF (ICUTTYP.EQ.3) THEN
	  WRITE(*,*) 'Give AIMCUT:'
	  READ(*,*) 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
	WRITE(*,*) 'Give ISHEETTYP ( 1="++" 2="-+" 3="+-" 4="--" ) for up" :'
	WRITE(*,*) '( irrelevant for rigid or free booundaries )'
	WRITE(*,*) '( "1 or 3" "2 or 4" irrelevant for a fluid half-space )'
	READ(*,*) ISHEETTYP
	WRITE(*,*) 'Give EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,',
     $	'EPSMAGNIFMAX :'
	WRITE(*,*) '( standard is  epsfakt=0.1, hfaktmax=1.05, hoff=6d-16 )'
	READ(*,*) EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,EPSMAGNIFMAX
	WRITE(*,*) 'Give EPSFSUB,HSTEPFSUBPRE :'
	WRITE(*,*) '( hstepfsub will be "hstepfsubpre/freq" , which '
	WRITE(*,*) '  is a "compromise" regarding rkra & rkrb )'
	READ(*,*) EPSFSUB,HSTEPFSUBPRE

	WRITE(*,*)
	WRITE(*,*) 'Give FREQMIN(>0.0),FREQINT (both in Hz) , NFR :'
	WRITE(*,*) '( "nfr1,nfr2 for cycling" if nfr=0 )'
	WRITE(*,*) '( " freqmin<0.0 & freqint=0.0 "  for  ',
     $	'"explicit ofreq(1:nfr) " )'
	WRITE(*,*) '( "freq in mainmodu" if nfr<0 )'
	WRITE(*,*) '< overread with "neg aladauer" below >'
	READ(*,*) FREQMIN,FREQINT, NFR
	IF (NFR.GE.0) THEN
	  IF (NFR.EQ.0) THEN
	    WRITE(*,*) 'Give NFR1,NFR2, NFR :'
	    WRITE(*,*) '( "cycle" for jf<nfr1 & jf>nfr2 )'
	    READ(*,*) NFR1,NFR2, NFR
	    IF (NFR.LE.0) STOP '%%% nfr'
	    JFCYCLE = .TRUE.
	  ELSE
	    JFCYCLE = .FALSE.
	  END IF
	  FNORMAL = .TRUE.
	  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
	    WRITE(*,*) 'give ofreq(1:nfr) increasing'
	    READ(*,*) (OFREQ(JF),JF=1,NFR)
	  END IF
	ELSE
	  JFCYCLE = .FALSE.
	  FNORMAL = .FALSE.
	END IF
	WRITE(*,*) 'Give ALADAUER (s), PREALA :'
	WRITE(*,*) '( "neg aladauer" for "freqdom" & ".J* i/o file names" )'
	WRITE(*,*) '( 0.0<preala<=1.0, preala=1.0 for real "omeg & rk" )'
	WRITE(*,*) '( "preala" is "temporary reduct factor during aladauer" )'
	WRITE(*,*) '( "neg preala" for "freqdom,halfband" below )'
	WRITE(*,*) '( "omeg" will be " 2*pi*freq + i*ala " '
	WRITE(*,*) '   where "ala" = "-ln(preala)/aladauer" )'
	READ(*,*) ALADAUER,PREALA
	IF (PREALA.LT.0.0D0) THEN
	  IF (ALADAUER.GE.0.0D0) STOP '%%% rmodpfw, aladauer & preala'
	END IF
	IF (PREALA.EQ.1.0D0) THEN
	  ALA = 0.0D0
	ELSE
	  ALA = -DLOG(DABS(PREALA))/DABS(ALADAUER)
	END IF
	WRITE(*,*) 'Give IUROT ( 0 for "u as specified", 1 for ',
     $	'"u rotation with -arg(omeg)" :'
	WRITE(*,*) '( u = "up rotated by -arg(omeg)" )'
	WRITE(*,*) '( irrelevant if "preala"=1.0 )'
	READ(*,*) IUROT

	IF (ALADAUER.LT.0.0D0) THEN
	  IF (.NOT.FNORMAL) STOP '%%% rmodpfw, fnormal and aladauer'
	  WRITE(*,*) 'Give FREQDOM (Hz) :'
	  IF (PREALA.LT.0.0D0) WRITE(*,*) '( give halfband as well )'
	  IF (PREALA.GE.0.0D0) THEN
            READ(*,*) FREQDOM
            HALFBAND = 0.0D0
          ELSE
            READ(*,*) FREQDOM,HALFBAND
            IF (FREQDOM.LE.0.0D0) STOP '%%% freqdom,halfband'
            IF (FREQDOM.LE.HALFBAND) STOP '%%% freqdom,halfband'
          END IF
	  IF (FREQDOM.GT.0.0D0) THEN
	    IF (HALFBAND.LE.0.0D0) THEN
	      NFR0 = 0.5D0 + (2*FREQDOM)*DABS(ALADAUER)
	    ELSE
	      NFR0 = 0.5D0 + (FREQDOM+HALFBAND)*DABS(ALADAUER)
	    END IF
	  ELSE
	    WRITE(*,*) 'Give NFR0 (=LX/2, integer power of 2) :'
	    READ(*,*) NFR0
	  END IF
	  NFR = NFR0-1
	  IF (NFR.GT.nfrdim) STOP 60
          SLASK = 1/DABS(ALADAUER)
	  DO JF = 1, NFR
	    OFREQ(JF) = (NFR+1-JF) * SLASK      ! decreasing
	  END DO
	  WRITE(*,*) 'equivalent "freqmin freqint nfr"  ', 
     &	   NFR*SLASK, -NFR*SLASK, NFR
	  IF (FREQDOM.GT.0.0D0) THEN
	    IF (HALFBAND.GT.0.0D0) THEN
	      DO JF = 1, NFR
	        IF ( OFREQ(JF) .LE. (FREQDOM-HALFBAND) ) THEN
	          WRITE(*,*) '!equivalent "freqmin freqint nfr"  ', 
     &	           NFR*SLASK, -(JF-1)*SLASK, JF-1
	          EXIT
	        END IF
	      END DO
	    END IF
	  END IF
	END IF

	WRITE(*,*)
	WRITE(*,*) 'Give "main name" MAINMODU for the nfr "infiles-mode" :'
	WRITE(*,*) '(freq indices as extensions)'
	MAINMODU = ' '
	DO WHILE (MAINMODU.EQ.' ') 
	  READ(*,'(A)') MAINMODU
	END DO
	CALL BLANKTOEND (MAINMODU)
	I = 42
	DO WHILE (MAINMODU(I:I).EQ.' ')
	  I = I-1
	END DO
	IPOSMODU = I+1

	WRITE(*,*)
	WRITE(*,*) 'Give IUTTYP ( 1=displ, 2=stress, 3="energy&power"-dens, '
	WRITE(*,*) '12=displ&stress, 123="all", 0="spec", -1,-2=+"m-coupl", '
	WRITE(*,*) '-3,-4=+"m-couplS", -2,-4="unformatted output") :'
	READ(*,*) IUTTYP
	IF (IUTTYP.EQ.1) THEN
	  IIA = 1
	  IIB = 2
	ELSE IF (IUTTYP.EQ.2) THEN
	  IIA = 3
	  IIB = 4
	ELSE IF (IUTTYP.EQ.3) THEN
	  IIA = 5
	  IIB = 7
	ELSE IF (IUTTYP.EQ.12) THEN
	  IIA = 1
	  IIB = 4
	ELSE IF (IUTTYP.EQ.123) THEN
	  IIA = 1
	  IIB = 7
	ELSE IF (IUTTYP.EQ.0) THEN
	  IF (.NOT.FNORMAL) STOP 41677
	  WRITE(*,*) 'Give KFWTYP (for ys comp; 1,2,3,4) :'
	  READ(*,*) KFWTYP
	  IIA = 100
	  IIB = 100
	ELSE IF ( IUTTYP.LE.-1 .AND. IUTTYP.GE.-4 ) THEN
	  IF (.NOT.FNORMAL) STOP '%%% iuttyp,fnormal'
	  WRITE(*,*) 'Give KFWTYP (for ys comp; 0, 1,2,3,4, 5,6,7; neg for QFIL) :'
	  READ(*,*) KFWTYP
	  IF (KFWTYP.GT.0) THEN
	    IIA = 101
	    IIB = 102
	  ELSE IF (KFWTYP.EQ.0) THEN
	    IIA = 102
	    IIB = 102
	  ELSE
	    IIA = 101
	    IIB = 101
	  END IF
	ELSE
	  STOP 1515
	END IF
	WRITE(*,*) 'Give "main name" MAINMODP for the nfr sets of ',
     $	'"outfiles-mode" :'
	WRITE(*,*) '(freq indices as extensions, plot with plotfkn)'
	READ(*,'(A)') MAINMODP
	CALL BLANKTOEND (MAINMODP)
	I = 42
	DO WHILE (MAINMODP(I:I).EQ.' ')
	  I = I-1
	END DO
	IPOSMODP = I+1

	OPEN(UNIT=20,FILE=FILM,STATUS='OLD')
	ISLASK = IAZIMI
	IAZIMI = 0              ! for this initial "particular" call
	CALL RMODEL ( 20,0, MHQ0 )
	IAZIMI = ISLASK         ! reset iazimi
	CLOSE(20)
	NHYD2 = NHYD
	ALLOCATE ( YS1FW(0:noldim-1), YS2FW(0:noldim-1), 
     &	 YS3FW(0:noldim-1), YS4FW(0:noldim-1) )  ! FWYSa
	ALLOCATE ( EXPTYP_FW(noldim-1), 
     &	 SLNREF_FW(noldim-1), E1X_FW(noldim-1), 
     &	 CP_FW(noldim-1), CQ_FW(noldim-1), STP_FW(noldim-1), 
     &	 STQ_FW(noldim-1), SDP_FW(noldim-1), SDQ_FW(noldim-1), 
     &	 TLNREF_FW(noldim-1), E2T_FW(noldim-1), 
     &	 URA_FW(noldim-1), URB_FW(noldim-1), 
     &	 CPX_FW(noldim-1) )  ! FWPARTa   ! cpx for Airy
	ALLOCATE ( CFAKT_FW(noldim-1), 
     &	 BCFAKT_FW(noldim-1), BCFAKTINV_FW(noldim-1), 
     &	 YA_FW(noldim-1), YB_FW(noldim-1) )  ! FWPARTa (NO equivalence)
	ALLOCATE ( DIS1FW(noldim-1), DIS2FW(noldim-1), 
     &	 DIS3FW(noldim-1), DIS4FW(noldim-1) )  ! FWDISa
	ALLOCATE ( TLNFW(noldim-1), YT1FW(noldim-1), 
     &	 YT2FW(noldim-1), YT3FW(noldim-1), YT4FW(noldim-1), 
     &	 YT5FW(noldim-1), YT6FW(noldim-1), 
     &	 TLNWF(noldim-1), YT1WF(noldim-1), YT2WF(noldim-1), 
     &	 YT3WF(noldim-1), YT4WF(noldim-1), YT5WF(noldim-1), 
     &	 YT6WF(noldim-1) )  ! FWYTa
	ALLOCATE ( TLNOPT(noldim), YT1OPT(noldim),
     &	 YT2OPT(noldim), YT3OPT(noldim), YT4OPT(noldim), 
     $	YT5OPT(noldim), YT6OPT(noldim) )  ! FWYTa (NO equivalence)
	ALLOCATE ( INDABLOK(0:noldim), RHOLOK(1:2,0:noldim), 
     &	 UA2LOK(1:2,0:noldim), UB2LOK(1:2,0:noldim) )
	
	WRITE(*,*)
	WRITE(*,*) '"initial rmodel" nss,nso,nsw,nol ',NSS,NSO,NSW,NOL
	WRITE(*,*) '                         nolhoml ',NOLHOML
	WRITE(*,*) 'Give INOLX (1="nolhoml with liqsixfw,solsixfw") :'
	WRITE(*,*) '[ > 1 for "+llmax" & "inolx=1" ]'
	WRITE(*,*) '[ neg for "-llmax" & "inolx=0" ]'
	READ(*,*) INOLX
	IF ( INOLX.EQ.0 .OR. INOLX.EQ.1 ) THEN
	  LLMAX = 0
	ELSE IF (INOLX.GT.1) THEN
	  LLMAX = INOLX
	  INOLX = 1
	ELSE
	  LLMAX = -INOLX
	  INOLX = 0
	END IF
	WRITE(*,*) 'Give IYSTAB (for ys stabilization; -1="ordinary", ',
     $	'1="x" ) :'
	READ(*,*) IYSTAB

	RW0 = (0.0D0,0.0D0)

	IF (.NOT.FNORMAL) THEN
	  FILMODU = MAINMODU
	  OPEN(UNIT=1,FILE=FILMODU,STATUS='OLD')
	  READ(1,*) ISLASK
	  DO I = 1, ISLASK
	    READ(1,*)
	  END DO
	  READ(1,*) ISLASK, ISLASK1,ISLASK2,SLASK,ISLASK3
	  IF (ISLASK.LT.-NFR) STOP 41561
	  IF ( ISLASK1.NE.ICUTTYP .OR. ISLASK2.NE.ISHEETTYP .OR. 
     $    ISLASK3.NE.IUROT ) STOP 34155
	  IF (ALA.EQ.0.0D0) THEN
	    IF (SLASK.NE.0.0D0) STOP '%%% ala inconsistency'
	  ELSE
	    IF ( DABS((SLASK/ALA)-1) .GT. 0.00001D0 ) 
     $	    STOP '%%% ala inconsistency'
	  END IF
	  DO I = 1, -NFR-1
	    READ(1,*) ISLASK
	    DO ISLASK1 = 1, ISLASK
	      READ(1,*) SLASK, CSLASK
	    END DO
	  END DO
	  READ(1,*) NFR
	  IF ( NFR.GT.nfrdim .OR. NFR.GT.lldim ) STOP 41713
	  DO I = 1, NFR
	    READ(1,*) OFREQ(I), UPARR(I)
	  END DO
	  CLOSE(1)
	END IF

	DO JF = 1, NFR
	  IF (JFCYCLE) THEN
	    IF (JF.LT.NFR1) CYCLE  ! next jf
	    IF (JF.GT.NFR2) CYCLE  ! next jf
	  END IF

	  FREQ = OFREQ(JF)
	  HSTEPFSUB = HSTEPFSUBPRE/FREQ
          OMEG = DCMPLX ( 2*PI*FREQ, ALA )

	  IF (ALADAUER.LT.0.0D0) THEN
	    IF (FREQDOM.GT.0.0D0) THEN
	      IF (HALFBAND.GT.0.0D0) THEN
	        IF ( FREQ .LE. (FREQDOM-HALFBAND) ) THEN
	          GOTO 99
	        ELSE IF ( FREQ .GE. (FREQDOM+HALFBAND) ) THEN
	          GOTO 99
	        END IF
	      END IF
	    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 * MAX0(1,MHQNOW)

	  IF ( IIMHQSOME.NE.0 .OR. IAZIMI.NE.0 ) THEN
	    MHQEFF = MHQNOW
	    CALL RMODEL ( 0,0, MHQEFF )
	  END IF
	  IF (NSOURCFW.EQ.-1) STOP '%%% nsourcfw (-1 not for rmodpfw)'
	  IF (NSOURCFW.LT.-1) THEN
C	    beam source array, set amsourcfw()
            ZMARR = 0.5D0 * ( ZD(NSOFW(1)) + ZD(NSOFW(-NSOURCFW)) )
            HZARR = 0.5D0 * DABS ( ZD(NSOFW(-NSOURCFW)) - ZD(NSOFW(1)) )
            IF (CUTARR.NE.1.0D0) THEN
              XARR = -DLOG(CUTARR) / (HZARR*HZARR)
            ELSE
              XARR = 0.0D0
            END IF
            SLWARR = DREAL ( CDSQRT(UA2SOURCFW(1)) )
            CSLASK0 = OMEG * DSIND(ANGARR) * SLWARR
            WEIGHTSUM = 0.0D0
	    DO I = 1, -NSOURCFW
              IF ( NSOFW(I).LE.NSS .OR. NSOFW(I).GT.NSW ) STOP '% nsofw()'
              IF ( IFSMTYP(I) .GT. -4 ) STOP '% ifsmtyp(), only homog here'
	      IF (I.GT.1) THEN
                SLASK =  CDABS ( UA2SOURCFW(I)/UA2SOURCFW(1) - 1 )
                IF (SLASK.GT.1.0D-10) STOP '% ua2sourcfw()'
                SLASK =  DABS ( RHOSOURCFW(I)/RHOSOURCFW(1) - 1 )
                IF (SLASK.GT.1.0D-10) STOP '% rhosourcfw()'
	      END IF
	      SLASK = ZD(NSOFW(I)) - ZMARR
              WEIGHTSUM = WEIGHTSUM + DEXP(-XARR*(SLASK*SLASK))
            END DO
            DO I = 1, -NSOURCFW
              SLASK = ZD(NSOFW(I)) - ZMARR
              CSLASK = CSLASK0 * SLASK  ! beamforming phase shift
              CSLASK = DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))     ! i*cslask
              SLASK = DEXP(-XARR*(SLASK*SLASK)) / WEIGHTSUM
              AMSOURCFW(I) = SLASK * CDEXP(CSLASK)
C             note weightsum to ensure "unity total weight"
            END DO
	  END IF

	  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
	  OMEG2 = OMEG*OMEG

	  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

	  IF ( JF.EQ.1 .OR. IIMHQSOME.NE.0 .OR. IAZIMI.NE.0 ) THEN
	    DO ID = 0, NOL-1
C	      indablok(id):
C	              0 "none"
C	             -1 "above" only
C	              1 "below" only
C	             -2 "both" and unequal
C	              2 "both" and equal
C	      ...lok(1,id):
C	             for "above"
C	      ...lok(2,id):
C	           for "below"
	      CALL PARAMZ ( ZD(ID) , INDABLOK(ID), UA2LOK(1,ID),
     $	      UB2LOK(1,ID),RHOLOK(1,ID) )
	      IF (INDABLOK(ID).EQ.0) THEN
		 WRITE(*,*) '%%% rmodpfw, id,indablok=',ID,INDABLOK(ID)
		 STOP 720
	      END IF
	    END DO
	    N1 = 0
	    IF (NSW.GT.0) THEN
	      IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	        DO ID = 0, NSS
                  N1 = N1 + 1
	        END DO
	      END IF
	      DO ID = NSS, NSW
		INKLUD = ISOLID(ID).EQ.0
		IF (ID.EQ.0) INKLUD = RHO(0).GT.0.0D0 .AND. INKLUD
	        IF (INKLUD) N1 = N1 + 1
		INKLUD = ISOLID(ID+1).EQ.0 .AND. INDABLOK(ID).NE.2
		IF ( ID. EQ. (NOL-1) ) INKLUD = RHO(NOL).GT.0.0D0 .AND. INKLUD
	        IF (INKLUD) N1 = N1 + 1
	      END DO
	    END IF
	    IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	      DO ID = NSW, NOL-1
                N1 = N1 + 1
	      END DO
	    END IF
	    N234 = 0
	    IF (NSW.GT.0) THEN
	      IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	        DO ID = 0, NSS
                  N234 = N234 + 1
	        END DO
	      END IF
	      DO ID = NSS, NSW
                N234 = N234 + 1
	      END DO
	    END IF
	    IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	      DO ID = NSW, NOL-1
                N234 = N234 + 1
	      END DO
	    END IF
	    N789 = 0
	    IF (NSW.GT.0) THEN
	      IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	        DO ID = 0, NSS
	          DO IAB = 1, 2
	            IF (IAB.EQ.1) THEN
	              INKLUD = ISOLID(ID).EQ.1
	            ELSE
	              INKLUD = ISOLID(ID+1).EQ.1 .AND. INDABLOK(ID).NE.2
	            END IF
	            IF (INKLUD) N789 = N789 + 1
	          END DO
	        END DO
	      END IF
	      DO ID = NSS, NSW
	        DO IAB = 1, 2
	          IF (IAB.EQ.1) THEN
	            INKLUD = ISOLID(ID).EQ.0
	            IF (ID.EQ.0) INKLUD = RHO(0).GT.0.0D0 .AND. INKLUD
	          ELSE
	            INKLUD = ISOLID(ID+1).EQ.0 .AND. INDABLOK(ID).NE.2
	          IF ( ID .EQ. (NOL-1) ) INKLUD = RHO(NOL).GT.0.0D0 .AND. 
     $              INKLUD
	          END IF
	          IF (INKLUD) N789 = N789 + 1
	        END DO
	      END DO
	    END IF
	    IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	      DO ID = NSW, NOL-1
	        DO IAB = 1, 2
	          IF (IAB.EQ.1) THEN
	            INKLUD = ISOLID(ID).EQ.1
	          ELSE
	            INKLUD = ISOLID(ID+1).EQ.1 .AND. INDABLOK(ID).NE.2
	          END IF
	          IF (INKLUD) N789 = N789 + 1
	        END DO
	      END DO
	    END IF
	  END IF

	  IF (FNORMAL) THEN
	    FILMODU = MAINMODU
	    IPOS = IPOSMODU
	    IF (IPOS.GT.42) STOP 462
	    FILMODU(IPOS:IPOS) = '.'
	    IPOS = IPOS+1
	    IF (ALADAUER.GE.0.0D0) THEN
	      ISLASK = FREQ
	      CALL PUTINTP (ISLASK,IPOS,FILMODU)
	      IF (ISLASK.NE.FREQ) THEN
	        IF (IPOS.GT.42) STOP 460
	        FILMODU(IPOS:IPOS) = 'X'
	        IPOS = IPOS+1
	        CALL PUTINTP (JF,IPOS,FILMODU)
	      END IF
	    ELSE
	      IF (IPOS.GT.42) STOP 460
	      FILMODU(IPOS:IPOS) = 'J'
	      IPOS = IPOS+1
	      CALL PUTINTP (JF,IPOS,FILMODU)
	    END IF
	    OPEN(UNIT=1,FILE=FILMODU,STATUS='OLD')
	    READ(1,*) ISLASK,SLASK
	    IF ( DABS((SLASK/FREQ)-1) .GT. 0.00001D0 ) STOP 41155
	    DO I = 1, ISLASK
	      READ(1,*)
	    END DO
	    READ(1,*) ISLASK, ISLASK1,ISLASK2,SLASK,ISLASK3
	    IF ( ISLASK.LT.1 .OR. ISLASK1.NE.ICUTTYP .OR. 
     $      ISLASK2.NE.ISHEETTYP .OR. ISLASK3.NE.IUROT ) STOP 34155
	    IF (ALA.EQ.0.0D0) THEN
	      IF (SLASK.NE.0.0D0) STOP '%%% ala inconsistency'
	    ELSE
	      IF ( DABS((SLASK/ALA)-1) .GT. 0.00001D0 ) 
     $	      STOP '%%% ala inconsistency'
	    END IF
	    READ(1,*) LL
	    IF (LLMAX.NE.0) LL = MIN0(LL,LLMAX)
	    IF (LL.GT.lldim) STOP '%%% lldim'
	    DO I = 1, LL
	      READ(1,*) SLASK1,SLASK2
	      UPARR(I) = DCMPLX(SLASK1,SLASK2)
	    END DO
	    CLOSE(1)
	  ELSE
	    LL = 1
	  END IF
	  
	  IF ( FNORMAL .OR. JF.EQ.1 ) THEN
	    DO II = IIA, IIB
	      IF ( II.LE.4 .OR. II.GE.100 ) THEN
	        IIW = II
	      ELSE
	        IIW = II+2
	      END IF
	      FILMODP = MAINMODP
	      IPOS = IPOSMODP
	      IF (II.LT.100) THEN
	        IF (IPOS.GT.42) STOP 460
	        FILMODP(IPOS:IPOS) = '_'
	        IPOS = IPOS+1
	        CALL PUTINTP (IIW,IPOS,FILMODP)
	      END IF
	      IF (FNORMAL) THEN
	        IF (IPOS.GT.42) STOP 460
	        FILMODP(IPOS:IPOS) = '.'
	        IPOS = IPOS+1
	        IF (ALADAUER.GE.0.0D0) THEN
	          ISLASK = FREQ
	          CALL PUTINTP (ISLASK,IPOS,FILMODP)
	          IF (ISLASK.NE.FREQ) THEN
	            IF (IPOS.GT.42) STOP 460
	            FILMODP(IPOS:IPOS) = 'X'
	            IPOS = IPOS+1
	            CALL PUTINTP (JF,IPOS,FILMODP)
	          END IF
	        ELSE
	          IF (IPOS.GT.42) STOP 460
	          FILMODP(IPOS:IPOS) = 'J'
	          IPOS = IPOS+1
	          CALL PUTINTP (JF,IPOS,FILMODP)
	        END IF
	      END IF
	      IF (II.EQ.102) THEN
	        IF (IPOS.GT.42) STOP 460
	        FILMODP(IPOS:IPOS) = '-'
	        IPOS = IPOS+1
	      END IF
	      ISLASK = 14
	      IF (IAZIMI.NE.0) ISLASK = ISLASK+1
	      IF (FNORMAL) ISLASK = ISLASK+1
	      IF ( IUTTYP.NE.-2 .AND. IUTTYP.NE.-4 ) THEN
	        OPEN(UNIT=IIW,FILE=FILMODP,STATUS='NEW')
	        IF ( II.LT.100 .OR. II.EQ.102 ) THEN
	          WRITE(IIW,*) ISLASK,SNGL(FREQ),'  ! lines to skip, freq'
	        ELSE
	          WRITE(IIW,*) ISLASK,SNGL(FREQ),IABS(KFWTYP),' !l-to-skip,frq,|kfwtyp|'
	        END IF
	        IF (II.EQ.1) THEN
	          WRITE(IIW,*) ' ( depth , horiz-displ )  from RMODPFW'
	        ELSE IF (II.EQ.2) THEN
	          WRITE(IIW,*) ' ( depth , vert-displ )  from RMODPFW'
	        ELSE IF (II.EQ.3) THEN
	          WRITE(IIW,*) ' ( depth , shear-stress )  from RMODPFW'
	        ELSE IF (II.EQ.4) THEN
	          WRITE(IIW,*) ' ( depth , normal-stress )  from RMODPFW'
	        ELSE IF (II.EQ.5) THEN
	          WRITE(IIW,*) ' ( depth , P energy/power dens )  from RMODPFW'
	        ELSE IF (II.EQ.6) THEN
	          WRITE(IIW,*) ' ( depth , S energy/power dens )  from RMODPFW'
	        ELSE IF (II.EQ.7) THEN
	          WRITE(IIW,*) ' ( depth , kinetic energy/power ',
     $	          'dens )  from RMODPFW'
	        ELSE IF ( II.EQ.100 .OR. II.EQ.101 ) THEN
	          WRITE(IIW,*) ' ( depth , "excited field" )  from RMODPFW'
	        ELSE IF (II.EQ.102) THEN
	          WRITE(IIW,*) ' ( depth , "field" )  from RMODPFW'
	        ELSE
	          STOP '%%% ii problem'
	        END IF
	        WRITE(IIW,*) ' FILM: ',FILM
	        WRITE(IIW,*) ' MHQ0,FREQMHQ0, IAZIMI: ',MHQ0,SNGL(FREQMHQ0), 
     $	         IAZIMI
	        IF (IAZIMI.NE.0) WRITE(IIW,*) '   ALFA1AZIM: ',SNGL(ALFA1AZIM)
	        IF (ICUTTYP.NE.3) THEN
	          WRITE(IIW,*) ' ICUTTYP: ',ICUTTYP
	        ELSE
	          WRITE(IIW,*) ' ICUTTYP,AIMCUT: ',ICUTTYP,SNGL(AIMCUT)
	        END IF
	        WRITE(IIW,*) ' ISHEETTYP: ',ISHEETTYP
	        WRITE(IIW,*) ' EPSGLOBLN,EPSGLOB: ',SNGL(EPSGLOBLN),SNGL(EPSGLOB)
	        WRITE(IIW,*) ' EPSFAKT,HFAKTMAX: ',SNGL(EPSFAKT),SNGL(HFAKTMAX)
	        WRITE(IIW,*) ' HOFF,EPSMAGNIFMAX: ',SNGL(HOFF),SNGL(EPSMAGNIFMAX)
	        WRITE(IIW,*) ' EPSFSUB,HSTEPFSUBPRE: ',SNGL(EPSFSUB),
     $	         SNGL(HSTEPFSUBPRE)
	        IF (FNORMAL) WRITE(IIW,*) ' FREQ: ',FREQ
	        WRITE(IIW,*) ' ALADAUER,PREALA: ',SNGL(ALADAUER),SNGL(PREALA)
	        WRITE(IIW,*) ' IUROT: ',IUROT
	        WRITE(IIW,*) ' FILMODU: ',FILMODU
	        WRITE(IIW,*) ' IUTTYP: ',IUTTYP
	        WRITE(IIW,*) ' INOLX,IYSTAB: ',INOLX,IYSTAB
	        WRITE(IIW,*)
	        IF (FNORMAL) THEN
	          IF (II.LT.100) THEN
	            WRITE(IIW,*) LL,'  !  ll'
	          ELSE IF (II.EQ.100) THEN
	            WRITE(IIW,'(3(X,I7),X,G24.16,X,I7,A)') 
     $	             LL,ICUTTYP,ISHEETTYP,ALA,IUROT,' !  ll, ...'
	          ELSE IF (II.EQ.101) THEN
C	            indicate "cslask,ys4fw() separation" by "-ll"
	            WRITE(IIW,'(3(X,I7),X,G24.16,X,I7,A)') 
     &	             -LL,ICUTTYP,ISHEETTYP,ALA,IUROT,' ! -ll, ...'
	          ELSE IF (II.NE.102) THEN
	            STOP '%%% ii problem'
	          END IF
	        ELSE
	          IF (II.GT.100) STOP '%%% iuttyp,ii,fnormal'
	          WRITE(IIW,*) NFR,'  !  number of modes'
	        END IF
	        WRITE(IIW,*) 
	      ELSE
	        OPEN(UNIT=IIW,FILE=FILMODP,STATUS='NEW',FORM='UNFORMATTED')
	        IF (II.EQ.102) THEN
	          WRITE(IIW) ISLASK,FREQ  ! lines to skip, freq
	        ELSE
	          WRITE(IIW) ISLASK,FREQ,IABS(KFWTYP)  ! l-to-skip,frq,|kfwtyp|
	        END IF
	        WRITE(IIW)  ! ( depth , "field" )  from RMODPFW
	        WRITE(IIW) FILM
	        WRITE(IIW) MHQ0,FREQMHQ0, IAZIMI
	        IF (IAZIMI.NE.0) WRITE(IIW) ALFA1AZIM
	        IF (ICUTTYP.NE.3) THEN
	          WRITE(IIW) ICUTTYP
	        ELSE
	          WRITE(IIW) ICUTTYP,AIMCUT
	        END IF
	        WRITE(IIW) ISHEETTYP
	        WRITE(IIW) EPSGLOBLN,EPSGLOB
	        WRITE(IIW) EPSFAKT,HFAKTMAX
	        WRITE(IIW) HOFF,EPSMAGNIFMAX
	        WRITE(IIW) EPSFSUB,HSTEPFSUBPRE
	        IF (FNORMAL) WRITE(IIW) FREQ
	        WRITE(IIW) ALADAUER,PREALA
	        WRITE(IIW) IUROT
	        WRITE(IIW) FILMODU
	        WRITE(IIW) IUTTYP
	        WRITE(IIW) INOLX,IYSTAB
	        IF (FNORMAL) THEN
	          IF (II.EQ.101) THEN
C	            indicate "cslask,ys4fw() separation" by "-ll"
	            WRITE(IIW) -LL,ICUTTYP,ISHEETTYP,ALA,IUROT  ! -ll, ...
	          ELSE IF (II.NE.102) THEN
	            STOP '%%% ii problem'
	          END IF
	        ELSE
	          STOP '%%% iuttyp,fnormal'
	        END IF
	      END IF
	    END DO  ! ii
	  END IF

	  ALLOCATE ( UGRARR(LL) )
C	  ------------------------------------------------------------------------------------------
	  DO I = 1, LL

	    IF (FNORMAL) THEN
	      UP = UPARR(I)
	    ELSE
	      UP = UPARR(JF)
	    END IF
	    U = U0*UP
	    UP2 = UP*UP
	    U2 = U*U

	    WRITE(*,*)
	    WRITE(*,*) 'NEW MODE  jf,nfr  i,ll are  ',JF,NFR,'  ',I,LL
	    WRITE(*,*) 'freq,up=',FREQ,UP

	    IF (INOLX.NE.0) THEN
	      NOLORIG = NOL
	      NSWORIG = NSW
	      NSOORIG = NSO
	      IFSMNOLORIG = IFSMTYP(NOL)
	      NOL = NOLHOML
	      IF (NOL.GT.noldim) STOP '%%% noldim'
	      NSW = MIN(NSW,NOL-1)
	      IF (NSWORIG.GT.0) THEN
	        IF (NSW.EQ.NSS) THEN
	          NSW = NSS+1  ! keep "nss.lt.nsw except when nsw=0"
	          IF (NSW.EQ.NOL) NOL = NSW+1  ! keep "nsw.lt.nol"
	          IF (DZ(NSW).NE.0.0D0) THEN
	            WRITE(*,*) 'requested "liqsixfw" will not be performed'
	            STOP 341567
	          END IF
	        END IF
	      END IF
	      NSO = MIN(NSO,NSW)  ! not essential, "adaptive nso" in dhcompfw2
	      IFSMTYP(NOL) = -4
	    END IF
	    CALL MODIVAN ( UP, IDMAX,ICMAX )
	    CALL FWDERIV ( UP, DLNDO,DERIVDO, DLNDU,DERIVDU )
	    IF (INOLX.NE.0) THEN
	      NOL = NOLORIG
	      NSW = NSWORIG
	      NSO = NSOORIG
	      IFSMTYP(NOL) = IFSMNOLORIG
	      IF (NOLHOML.LT.NOL) THEN
C	        "track correctly on downgoing exponentials" 
C	        without noise-interference from the upgoing exponentials
	        IF (ISOLID(NOL).NE.0) THEN
                  CALL SOLSIXFW ( NOLHOML,NOL-1, UP,UP2,U,U2 )
	        ELSE
                  CALL LIQSIXFW ( NOLHOML,NOL-1, UP,UP2,U,U2 )
	        END IF
	      END IF
	    END IF

	    IF ( IUTTYP.EQ.1 .OR. IUTTYP.GT.3 ) THEN
	      WRITE(1,*) 
	      WRITE(2,*) 
	      I2 = I
	      CESLASK = UP
	      WRITE(1,*) N1   ,' ! i,freq,up: ',I2,SNGL(FREQ),CESLASK
	      WRITE(2,*) N234 ,' ! i,freq,up: ',I2,SNGL(FREQ),CESLASK
	      IF (NSW.GT.0) THEN
	        IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	          DO ID = 0, NSS
	            CSLASK = YS1FW(ID)/OMEG
	            WRITE(1,*) SNGL(ZD(ID)), CSLASK
	            CSLASK = YS2FW(ID)/OMEG
	            WRITE(2,*) SNGL(ZD(ID)), CSLASK
	          END DO
	        END IF
	        DO ID = NSS, NSW
		  INKLUD = ISOLID(ID).EQ.0
		  IF (ID.EQ.0) INKLUD = RHO(0).GT.0.0D0 .AND. INKLUD
	          IF (INKLUD) THEN
	            CSLASK = (U/RHOLOK(1,ID))*YS4FW(ID)/OMEG
	            WRITE(1,*) SNGL(ZD(ID)), CSLASK
	          END IF
		  INKLUD = ISOLID(ID+1).EQ.0 .AND. INDABLOK(ID).NE.2
		  IF ( ID .EQ. (NOL-1) ) INKLUD = RHO(NOL).GT.0.0D0 .AND. 
     $            INKLUD
	          IF (INKLUD) THEN
	            CSLASK = (U/RHOLOK(2,ID))*YS4FW(ID)/OMEG
	            WRITE(1,*) SNGL(ZD(ID)), CSLASK
	          END IF
	          CSLASK = YS2FW(ID)/OMEG
	          WRITE(2,*) SNGL(ZD(ID)), CSLASK
	        END DO
	      END IF
	      IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	        DO ID = NSW, NOL-1
	          CSLASK = YS1FW(ID)/OMEG
	          WRITE(1,*) SNGL(ZD(ID)), CSLASK
	          CSLASK = YS2FW(ID)/OMEG
	          WRITE(2,*) SNGL(ZD(ID)), CSLASK
	        END DO
	      END IF
	    END IF

	    IF ( IUTTYP.EQ.2 .OR. IUTTYP.GT.3 ) THEN
	      WRITE(3,*) 
	      WRITE(4,*) 
	      I2 = I
	      CESLASK = UP
	      WRITE(3,*) N234,' ! i,freq,up: ',I2,SNGL(FREQ),CESLASK
	      WRITE(4,*) N234,' ! i,freq,up: ',I2,SNGL(FREQ),CESLASK
	      IF (NSW.GT.0) THEN
	        IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	          DO ID = 0, NSS
	            CSLASK = YS3FW(ID)
	            WRITE(3,*) SNGL(ZD(ID)), CSLASK
	            CSLASK = YS4FW(ID)
	            WRITE(4,*) SNGL(ZD(ID)), CSLASK
	          END DO
	        END IF
	        DO ID = NSS, NSW
	          WRITE(3,*) SNGL(ZD(ID)), (0.0,0.0)
	          CSLASK = YS4FW(ID)
	          WRITE(4,*) SNGL(ZD(ID)), CSLASK
	        END DO
	      END IF
	      IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	        DO ID = NSW, NOL-1
	          CSLASK = YS3FW(ID)
	          WRITE(3,*) SNGL(ZD(ID)), CSLASK
	          CSLASK = YS4FW(ID)
	          WRITE(4,*) SNGL(ZD(ID)), CSLASK
	        END DO
	      END IF
	    END IF

	    IF ( IUTTYP.EQ.3 .OR. IUTTYP.EQ.123 ) THEN
	      WRITE(7,*) 
	      WRITE(8,*) 
	      WRITE(9,*) 
	      I2 = I
	      CESLASK = UP
	      WRITE(7,*) N789,' ! i,freq,up: ',I2,SNGL(FREQ),CESLASK
	      WRITE(8,*) N789,' ! i,freq,up: ',I2,SNGL(FREQ),CESLASK
	      WRITE(9,*) N789,' ! i,freq,up: ',I2,SNGL(FREQ),CESLASK
C	      evaluate  " (lambda+(2.0/3)*my) * e(i,i)*dconjg(e(k,k)) "
C	      evaluate  " (2*my) * edev(i,k)*dconjg(edev(i,k)) "
C	      energy-density as real parts, power-density as imag parts
C	      kinetic as  " rho*omeg*omeg * u(i)*dconjg(u(i)) "
	      IF (NSW.GT.0) THEN
	        IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	          DO ID = 0, NSS
	            DO IAB = 1, 2
	              IF (IAB.EQ.1) THEN
	                INKLUD = ISOLID(ID).EQ.1
	              ELSE
	                INKLUD = ISOLID(ID+1).EQ.1 .AND. INDABLOK(ID).NE.2
	              END IF
	              IF (INKLUD) THEN
	                CSLASK = ( CDABS(UA2LOK(IAB,ID))**2 / 
     $	                RHOLOK(IAB,ID) ) * ( 1/UA2LOK(IAB,ID) - 
     $	                4.0D0/(3*UB2LOK(IAB,ID)) ) * CDABS ( YS4FW(ID) + 
     $	                (2*RHOLOK(IAB,ID)/UB2LOK(IAB,ID))*U*YS1FW(ID) )**2
	                CSLASK = DCMPLX ( 0.25D0*DREAL(CSLASK) , 
     $	                DIMAG(-0.5D0*OMEG*CSLASK) )
	                WRITE(7,*) SNGL(ZD(ID)), CSLASK
	                CSLASK = (2*RHOLOK(IAB,ID))/UB2LOK(IAB,ID)
	                CSLASK = CSLASK * ( CDABS(U2) * 
     $	                CDABS(YS1FW(ID))**2 + (CDABS(UA2LOK(IAB,ID))/
     $	                RHOLOK(IAB,ID))**2 * ( CDABS(YS4FW(ID)-
     $	                (RHOLOK(IAB,ID)/UA2LOK(IAB,ID)-CSLASK)*U*
     $	                YS1FW(ID))**2 - (1.0D0/3)*CDABS(YS4FW(ID)+
     $	                CSLASK*U*YS1FW(ID))**2 ) ) + 
     $	                DCONJG(UB2LOK(IAB,ID))/RHOLOK(IAB,ID) * 
     $	                CDABS(YS3FW(ID))**2
	                CSLASK = DCMPLX ( 0.25D0*DREAL(CSLASK) , 
     $	                DIMAG(-0.5D0*OMEG*CSLASK) )
	                WRITE(8,*) SNGL(ZD(ID)), CSLASK
	                CSLASK = 0.25D0 * RHOLOK(IAB,ID) * 
     $	                ( CDABS(YS1FW(ID))**2 + CDABS(YS2FW(ID))**2 )
	                WRITE(9,*) SNGL(ZD(ID)), CSLASK
	              END IF
	            END DO
	          END DO
	        END IF
	        DO ID = NSS, NSW
	          DO IAB = 1, 2
	            IF (IAB.EQ.1) THEN
	              INKLUD = ISOLID(ID).EQ.0
		      IF (ID.EQ.0) INKLUD = RHO(0).GT.0.0D0 .AND. INKLUD
	            ELSE
	              INKLUD = ISOLID(ID+1).EQ.0 .AND. INDABLOK(ID).NE.2
		      IF ( ID .EQ. (NOL-1) ) INKLUD = RHO(NOL).GT.0.0D0 .AND. 
     $                INKLUD
	            END IF
	            IF (INKLUD) THEN
	              CSLASK = (DCONJG(UA2LOK(IAB,ID))/RHOLOK(IAB,ID)) * 
     $	              CDABS(YS4FW(ID))**2
	              CSLASK = DCMPLX ( 0.25D0*DREAL(CSLASK) , 
     $	              DIMAG(-0.5D0*OMEG*CSLASK) )
	              WRITE(7,*) SNGL(ZD(ID)), CSLASK
	              WRITE(8,*) SNGL(ZD(ID)), (0.0,0.0)
	              CSLASK = (U/RHOLOK(IAB,ID))*YS4FW(ID)
	              CSLASK = 0.25D0 * RHOLOK(IAB,ID) * 
     $	              ( CDABS(CSLASK)**2 + CDABS(YS2FW(ID))**2 )
	              WRITE(9,*) SNGL(ZD(ID)), CSLASK
	            END IF
	          END DO
	        END DO
	      END IF
	      IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	        DO ID = NSW, NOL-1
	          DO IAB = 1, 2
	            IF (IAB.EQ.1) THEN
	              INKLUD = ISOLID(ID).EQ.1
	            ELSE
	              INKLUD = ISOLID(ID+1).EQ.1 .AND. INDABLOK(ID).NE.2
	            END IF
	            IF (INKLUD) THEN
	              CSLASK = ( CDABS(UA2LOK(IAB,ID))**2 / 
     $	              RHOLOK(IAB,ID) ) * ( 1/UA2LOK(IAB,ID) - 
     $	              4.0D0/(3*UB2LOK(IAB,ID)) ) * CDABS ( YS4FW(ID) + 
     $	              (2*RHOLOK(IAB,ID)/UB2LOK(IAB,ID))*U*YS1FW(ID) )**2
	              CSLASK = DCMPLX ( 0.25D0*DREAL(CSLASK) , 
     $	              DIMAG(-0.5D0*OMEG*CSLASK) )
	              WRITE(7,*) SNGL(ZD(ID)), CSLASK
	              CSLASK = (2*RHOLOK(IAB,ID))/UB2LOK(IAB,ID)
	              CSLASK = CSLASK * ( CDABS(U2) * 
     $	              CDABS(YS1FW(ID))**2 + (CDABS(UA2LOK(IAB,ID))/
     $	              RHOLOK(IAB,ID))**2 * ( CDABS(YS4FW(ID)-
     $	              (RHOLOK(IAB,ID)/UA2LOK(IAB,ID)-CSLASK)*U*
     $	              YS1FW(ID))**2 - (1.0D0/3)*CDABS(YS4FW(ID)+
     $	              CSLASK*U*YS1FW(ID))**2 ) ) + 
     $	              DCONJG(UB2LOK(IAB,ID))/RHOLOK(IAB,ID) * 
     $	              CDABS(YS3FW(ID))**2
	              CSLASK = DCMPLX ( 0.25D0*DREAL(CSLASK) , 
     $	              DIMAG(-0.5D0*OMEG*CSLASK) )
	              WRITE(8,*) SNGL(ZD(ID)), CSLASK
	              CSLASK = 0.25D0 * RHOLOK(IAB,ID) * 
     $	              ( CDABS(YS1FW(ID))**2 + CDABS(YS2FW(ID))**2 )
	              WRITE(9,*) SNGL(ZD(ID)), CSLASK
	            END IF
	          END DO
	        END DO
	      END IF
	    END IF

	    IF (IUTTYP.EQ.0) THEN
	      WRITE(100,*) 
	      I2 = I
	      WRITE(100,*) NHYD2,UP,' ! i:',I2
	      IF (NSOURCFW.EQ.0) THEN
C               source at the matching depth
	        IF (NSW.EQ.0) STOP 34443  ! 'nss-nsw' fluid must be present
	        CSLASK = YS4FW(NSO) * UA2SOURC / RHOSOURC
	      ELSE
C               source array
	        CSLASK = (0.0D0,0.0D0)
	        DO II = 1, IABS(NSOURCFW)  ! note stop when nsourcfw.eq.-1
	          ID = NSOFW(II)
	          IF ( NSS.LT.ID .AND. ID.LE.NSW ) THEN
	            CSLASK = CSLASK + 
     $	            AMSOURCFW(II) * (UA2SOURCFW(II)/RHOSOURCFW(II)) * YS4FW(ID)
	          ELSE IF ( ID.LT.NSS .OR. NSW.LT.ID ) THEN
	            CSLASK = CSLASK + 
     $	            AMSOURCFW(II) * (UA2SOURCFW(II)/RHOSOURCFW(II)) * 
     $	            ( YS4FW(ID) + 2*RHOSOURCFW(II)*U*
     $	            YS1FW(ID)/UB2SOURCFW(II) )
		  ELSE
		    STOP 411667
	          END IF
	        END DO
	      END IF
	      IF ( DERIVDU .NE. (0.0D0,0.0D0) ) THEN
	        CSLASK = DEXP(-DLNDU) * CSLASK / DERIVDU
	      ELSE
	        IF ( RHO(NSS).GE.0.0D0 .OR. RHO(NSS).EQ.-99.0D0 ) STOP '%derivdu,rho'
	        IF ( RHO(NSW+1).GE.0.0D0 .OR. RHO(NSW+1).EQ.-99.0D0 ) STOP '%derivdu,rho'
	        IF ( ZD(NSS) .NE. ZD(NSW) ) STOP '%derivdu,zd'
	        CSLASK = (0.0D0,0.0D0)  ! avoid "NaN" for "trivial one-mode medium"
	      END IF
	      IF (KFWTYP.EQ.1) THEN
		CSLASK = CSLASK/OMEG
	        DO IH = 1, NHYD
	          ID = IDHYD(IH)
		  IF ( ID.LE.NSS .AND. (NSS.GT.0.OR.ISOLID(0).NE.0) ) THEN
	            WRITE(100,*) ZD(ID), CSLASK*YS1FW(ID) !r1=ys1/omeg
	          ELSE IF ( ID.GE.NSW .AND. 
     $            ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) ) THEN
	            WRITE(100,*) ZD(ID), CSLASK*YS1FW(ID) !r1=ys1/omeg
		  ELSE
	            WRITE(100,*) ZD(ID), (0.0D0,0.0D0),' !not computed'
	          END IF
	        END DO
	      ELSE IF (KFWTYP.EQ.2) THEN
		CSLASK = CSLASK/OMEG
	        DO IH = 1, NHYD
	          ID = IDHYD(IH)
	          WRITE(100,*) ZD(ID), CSLASK*YS2FW(ID) ! r2 = ys2/omeg
	        END DO
	      ELSE IF (KFWTYP.EQ.3) THEN
	        DO IH = 1, NHYD
	          ID = IDHYD(IH)
	          IF ( ID.LE.NSS .AND. (NSS.GT.0.OR.ISOLID(0).NE.0) ) THEN
	            WRITE(100,*) ZD(ID), CSLASK*YS3FW(ID)  ! r3 = ys3
	          ELSE IF ( ID.GE.NSW .AND. 
     $            ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) ) THEN
	            WRITE(100,*) ZD(ID), CSLASK*YS3FW(ID)  ! r3 = ys3
	          ELSE
	            WRITE(100,*) ZD(ID), (0.0D0,0.0D0) !r3 zero in fluid
	          END IF
	        END DO
	      ELSE IF (KFWTYP.EQ.4) THEN
	        DO IH = 1, NHYD
	          ID = IDHYD(IH)
	          WRITE(100,*) ZD(ID), CSLASK*YS4FW(ID)    ! r4 = ys4
	        END DO
	      ELSE
	        STOP 451617
	      END IF
	    END IF

	    IF ( IUTTYP.EQ.-1 .OR. IUTTYP.EQ.-2 ) THEN
	      IF ( NSS.GT.0 .OR. RHO(0).GT.0.0D0 ) STOP '%%% iuttyp_-1,-2'
	      IF ( RHO(0) .EQ. -99.0D0 ) STOP '%%% iuttyp_-1,-2'  ! not impl
	      IF (NSW.EQ.0) STOP '%%% iuttyp_-1,-2 non-trivial fluid missing'
	      IF ( (NSW+1).NE.NOL .OR. ISOLID(NOL).NE.0 ) STOP '%iuttyp_-1,-2'
	      IF ( KFWTYP.NE.0 .AND. IABS(KFWTYP).NE.4 .AND. IABS(KFWTYP).NE.2 )
     &	       STOP '%iuttyp_-1,-2, kfwtyp'

	      IF (KFWTYP.NE.0) THEN
C	      ---
	      IF (KFWTYP.LT.0) THEN
	        IF (I.EQ.1) THEN
	          OPEN(102,FILE='QFIL',STATUS='OLD')
	          READ(102,*) ISLASK
	          IF (ISLASK.NE.LL) STOP '%%% qfil, ll'
	        END IF
	        READ(102,*) SLASK, ISLASK
	        SLASK1 = DABS(ZD(0)-SLASK)
	        ISLASK1 = 0
	        DO ID = 1, NSW
	          SLASK2 = DABS(ZD(ID)-SLASK)
	          IF (SLASK2.LT.SLASK1) THEN
	            SLASK1 = SLASK2
	            ISLASK1 = ID
	          END IF
	        END DO
	        IF ( SLASK1 .GT. 1.0D-8 ) STOP '%%% qfil, zd() match problem'
	        ID = ISLASK1
	        IF (ISLASK.EQ.2) THEN
		  QFAKT = YS2FW(ID)
	        ELSE IF (ISLASK.EQ.4) THEN
		  QFAKT = YS4FW(ID)
		ELSE
	          STOP '%%% qfil, icmax'
		END IF
	        IF ( QFAKT .EQ. (0.0D0,0.0D0) ) STOP '%%% qfakt'
	        QFAKT = 1/QFAKT
		WRITE(199,*) I, QFAKT
	      ELSE
	        QFAKT = (1.0D0,0.0D0)
	      END IF
	      I2 = I
C             source at the matching depth
	      IF (NSOURCFW.EQ.0) THEN
C               source at the matching depth
	        IF (NSW.EQ.0) STOP 34443  ! 'nss-nsw' fluid must be present
	        CSLASK = YS4FW(NSO) * UA2SOURC / RHOSOURC
	      ELSE
C               source array
	        CSLASK = (0.0D0,0.0D0)
	        DO II = 1, IABS(NSOURCFW)  ! note stop when nsourcfw.eq.-1
	          ID = NSOFW(II)
	          IF ( NSS.LT.ID .AND. ID.LE.NSW ) THEN
	            CSLASK = CSLASK + 
     $	            AMSOURCFW(II) * (UA2SOURCFW(II)/RHOSOURCFW(II)) * YS4FW(ID)
	          ELSE IF ( ID.LT.NSS .OR. NSW.LT.ID ) THEN
C	            STOP '% iuttyp=-1,-2; source array outside fluid'
		  ELSE
		    STOP 411667
	          END IF
	        END DO
	      END IF
	      IF ( DERIVDU .NE. (0.0D0,0.0D0) ) THEN
	        CSLASK = DEXP(-DLNDU) * CSLASK / DERIVDU
	      ELSE
	        IF ( RHO(NSS).GE.0.0D0 .OR. RHO(NSS).EQ.-99.0D0 ) STOP '%derivdu,rho'
	        IF ( RHO(NSW+1).GE.0.0D0 .OR. RHO(NSW+1).EQ.-99.0D0 ) STOP '%derivdu,rho'
	        IF ( ZD(NSS) .NE. ZD(NSW) ) STOP '%derivdu,zd'
	        CSLASK = (0.0D0,0.0D0)  ! avoid "NaN" for "trivial one-mode medium"
	      END IF
	      IF (IUTTYP.EQ.-1) THEN
	        WRITE(101,*) 
	        WRITE(101,*) I2, CSLASK/QFAKT, ' ! i'
	        WRITE(101,*) NHYD2,UP,' ! i:',I2
	        IF ( IABS(KFWTYP) .EQ. 4 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101,*) ZD(ID), QFAKT*YS4FW(ID)  !r4 = ys4 "without cslask"
	          END DO
	        ELSE IF ( IABS(KFWTYP) .EQ. 2 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101,*) ZD(ID), QFAKT*YS2FW(ID)/OMEG  ! r2 = ys2/omeg
	          END DO                                       ! "without cslask"
	        ELSE
		  STOP '%%% iuttyp_-1,-2, kfwtyp'
	        END IF
	      ELSE
	        WRITE(101) I, CSLASK/QFAKT, ' ! i'
	        WRITE(101) NHYD,UP,' ! i:',I
	        IF ( IABS(KFWTYP) .EQ. 4 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101) ZD(ID), QFAKT*YS4FW(ID)   ! r4 = ys4 "without cslask"
	          END DO
	        ELSE IF ( IABS(KFWTYP) .EQ. 2 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101) ZD(ID), QFAKT*YS2FW(ID)/OMEG !r2=ys2/omeg utan cslask
	          END DO
	        ELSE
		  STOP '%%% iuttyp_-1,-2, kfwtyp'
	        END IF
	      END IF  ! iuttyp
C	      ---
	      ELSE
C	      ---
	      IF (I.EQ.1) THEN
	        OPEN(101,FILE='QFIL',STATUS='UNKNOWN')
		WRITE(101,*) LL
	      END IF
	      IF ( IDMAX.LT.0 .OR. ICMAX.EQ.0 ) STOP '%%% vanishing mode'
	      WRITE(101,*) ZD(IDMAX), ICMAX
C	      ---
	      END IF  ! kfwtyp
	      
	      IF (KFWTYP.GE.0) THEN
	      I2 = I
	      IF (IUTTYP.EQ.-1) THEN
	        IF (I.EQ.1) THEN 
	          WRITE(102,*) 
	          WRITE(102,*) OMEG,' ! omeg'
	          WRITE(102,*) U0,' ! u0'
	          WRITE(102,*) RK0,' ! rk0'
	          WRITE(102,*) ZD(0),' ! zd(0)'
	          IF ( RHO(NOL).LE.0.0D0 .AND. RHO(NOL).NE.-99.0D0 ) THEN
	            WRITE(102,*) NSW,' ! nsw'
	            ISLASK = NSW
	          ELSE
	            WRITE(102,*) -NSW,' ! -nsw' !halfspace/zlocreactnol indic
	            ISLASK = NOL  ! nsw+1
	          END IF
	          DO ID = 1, ISLASK
	            IF ( IFSMTYP(ID) .GT. -4 ) STOP '%iuttyp_-1,-2, ifsmtyp()'
	            WRITE(102,*) ZD(ID), RHO(ID)
	            WRITE(102,*) UA2(ID)
	          END DO
	          IF ( RHO(NOL) .EQ. -99.0D0 ) THEN
	            IF ( UA2(NOL) .NE. (0.0D0,0.0D0) )
     &	             STOP '%iuttyp_-1,-2, zlocreactnol problem'
	            WRITE(102,*) ZLOCREACTNOL,' ! zlocreactnol'
	          END IF
	          WRITE(102,*) 
	          WRITE(102,'(3(X,I7),X,G24.16,X,I7,A)') 
     $	           LL,ICUTTYP,ISHEETTYP,ALA,IUROT,' !  ll, ...'
	        END IF
	        WRITE(102,*) 
	        WRITE(102,*) I2,UP,' ! i,up'
	        WRITE(102,*) DLNDU,DERIVDU
	        DO ID = 0, NSW
	          WRITE(102,*) YS4FW(ID),YS2FW(ID) !r4=ys4,r2=ys2/omeg
	        END DO
	      ELSE
	        IF (I.EQ.1) THEN 
	          WRITE(102) OMEG    ! omeg
	          WRITE(102) U0      ! u0
	          WRITE(102) RK0     ! rk0
	          WRITE(102) ZD(0)   ! zd(0)
	          IF ( RHO(NOL).LE.0.0D0 .AND. RHO(NOL).NE.-99.0D0 ) THEN
	            WRITE(102) NSW   ! nsw
	            ISLASK = NSW
	          ELSE
	            WRITE(102) -NSW  ! halfspace/zlocreactnol indic
	            ISLASK = NOL  ! nsw+1
	          END IF
	          DO ID = 1, ISLASK
	            IF ( IFSMTYP(ID) .GT. -4 ) STOP '%iuttyp_-1,-2, ifsmtyp()'
	            WRITE(102) ZD(ID), RHO(ID)
	            WRITE(102) UA2(ID)
	          END DO
	          IF ( RHO(NOL) .EQ. -99.0D0 ) THEN
	            IF ( UA2(NOL) .NE. (0.0D0,0.0D0) )
     &	             STOP '%iuttyp_-1,-2, zlocreactnol problem'
	            WRITE(102) ZLOCREACTNOL  ! zlocreactnol
	          END IF
	          WRITE(102) LL,ICUTTYP,ISHEETTYP,ALA,IUROT  !  ll, ...
	        END IF
	        WRITE(102) I,UP  ! i,up  !!! avoid "integer*2" i2
	        WRITE(102) DLNDU,DERIVDU
	        DO ID = 0, NSW
	          WRITE(102) YS4FW(ID),YS2FW(ID)  ! r4=ys4,r2=ys2/omeg
	        END DO
	      END IF  ! iuttyp
	      END IF  ! kfwtyp

	    END IF  ! iuttyp.eq.-1 .or. iuttyp-eq.-2

	    IF ( IUTTYP.EQ.-3 .OR. IUTTYP.EQ.-4 ) THEN
	      IF ( RHO(0) .GT. 0.0D0 ) STOP '%%% iuttyp_-3,-4'
	      IF ( RHO(0) .EQ. -99.0D0 ) STOP '%%% iuttyp_-3,-4'  ! not impl
	      IF ( RHO(NOL) .EQ. -99.0D0 ) STOP '%%% iuttyp_-3,-4'  ! not impl

	      IF (KFWTYP.NE.0) THEN
C	      ---
	      IF (KFWTYP.LT.0) THEN
	        IF (I.EQ.1) THEN
	          OPEN(102,FILE='QFIL',STATUS='OLD')
	          READ(102,*) ISLASK
	          IF (ISLASK.NE.LL) STOP '%%% qfil, ll'
	        END IF
	        READ(102,*) SLASK, ISLASK
	        IF (ISLASK.LT.0) THEN
	          IF (NSS.GT.0) THEN
	            SLASK1 = DABS(ZD(0)-SLASK)
	            ISLASK1 = 0
	            DO ID = 1, NSS
	              SLASK2 = DABS(ZD(ID)-SLASK)
	              IF (SLASK2.LT.SLASK1) THEN
	                SLASK1 = SLASK2
	                ISLASK1 = ID
	              END IF
	            END DO
	          END IF
	          IF ( (NSW+1) .LT. NOL ) THEN
	            IF (NSS.EQ.0) THEN     
	              SLASK1 = DABS(ZD(NSW)-SLASK)
	              ISLASK1 = NSW
	            END IF
	            DO ID = NSW, NOL-1
	              SLASK2 = DABS(ZD(ID)-SLASK)
	              IF (SLASK2.LT.SLASK1) THEN
	                SLASK1 = SLASK2
	                ISLASK1 = ID
	              END IF
	            END DO
	          END IF
	        ELSE IF (ISLASK.GT.0) THEN
	          IF (NSS.EQ.0) THEN
	            SLASK1 = DABS(ZD(0)-SLASK)
	            ISLASK1 = 0
	          ELSE IF ( (NSW+1) .EQ. NOL ) THEN
	            SLASK1 = DABS(ZD(NSW)-SLASK)
	            ISLASK1 = NSW
	          ELSE
	            IF ( (NSS+1) .EQ. NSW ) STOP '%%% qfil, nss,nsw'
	            SLASK1 = DABS(ZD(NSS+1)-SLASK)
	            ISLASK1 = NSS+1
	          END IF
	          DO ID = NSS+1, NSW-1
	            SLASK2 = DABS(ZD(ID)-SLASK)
	            IF (SLASK2.LT.SLASK1) THEN
	              SLASK1 = SLASK2
	              ISLASK1 = ID
	            END IF
	          END DO
	          IF ( (NSW+1) .EQ. NOL ) THEN
	            SLASK2 = DABS(ZD(NSW)-SLASK)
	            IF (SLASK2.LT.SLASK1) THEN
	              SLASK1 = SLASK2
	              ISLASK1 = NSW
	            END IF
	          END IF
	        ELSE
	          STOP '%%% qfil, vanishing mode'
	        END IF
	        IF ( SLASK1 .GT. 1.0D-8 ) STOP '%%% qfil, zd() match problem'
	        ID = ISLASK1
	        IF ( IABS(ISLASK) .EQ. 2 ) THEN
		  QFAKT = YS2FW(ID)
	        ELSE IF ( IABS(ISLASK) .EQ. 4 ) THEN
		  QFAKT = YS4FW(ID)
		ELSE IF (ISLASK.EQ.-1) THEN
		  QFAKT = YS1FW(ID)
		ELSE IF (ISLASK.EQ.-3) THEN
		  QFAKT = YS3FW(ID)
	        ELSE
		  STOP '%%% qfil, icmax'
		END IF
	        IF ( QFAKT .EQ. (0.0D0,0.0D0) ) STOP '%%% qfakt'
	        QFAKT = 1/QFAKT
	      ELSE
	        QFAKT = (1.0D0,0.0D0)
	      END IF
	      I2 = I
	      IF (NSOURCFW.EQ.0) THEN
C               source at the matching depth
	        IF (NSW.EQ.0) STOP 34443  ! 'nss-nsw' fluid must be present
	        CSLASK = YS4FW(NSO) * UA2SOURC / RHOSOURC
	      ELSE
C               source array
	        CSLASK = (0.0D0,0.0D0)
	        DO II = 1, IABS(NSOURCFW)  ! note stop when nsourcfw.eq.-1
	          ID = NSOFW(II)
	          IF ( NSS.LT.ID .AND. ID.LE.NSW ) THEN
	            CSLASK = CSLASK + 
     $	            AMSOURCFW(II) * (UA2SOURCFW(II)/RHOSOURCFW(II)) * YS4FW(ID)
	          ELSE IF ( ID.LT.NSS .OR. NSW.LT.ID ) THEN
	            CSLASK = CSLASK + 
     $	            AMSOURCFW(II) * (UA2SOURCFW(II)/RHOSOURCFW(II)) * 
     $	            ( YS4FW(ID) + 2*RHOSOURCFW(II)*U*
     $	            YS1FW(ID)/UB2SOURCFW(II) )
		  ELSE
		    STOP 411667
	          END IF
	        END DO
	      END IF
	      IF ( DERIVDU .NE. (0.0D0,0.0D0) ) THEN
	        CSLASK = DEXP(-DLNDU) * CSLASK / DERIVDU
	      ELSE
	        IF ( RHO(NSS).GE.0.0D0 .OR. RHO(NSS).EQ.-99.0D0 ) STOP '%derivdu,rho'
	        IF ( RHO(NSW+1).GE.0.0D0 .OR. RHO(NSW+1).EQ.-99.0D0 ) STOP '%derivdu,rho'
	        IF ( ZD(NSS) .NE. ZD(NSW) ) STOP '%derivdu,zd'
	        CSLASK = (0.0D0,0.0D0)  ! avoid "NaN" for "trivial one-mode medium"
	      END IF
	      IF (IUTTYP.EQ.-3) THEN
	        WRITE(101,*) 
	        WRITE(101,*) I2, CSLASK/QFAKT, ' ! i'
	        WRITE(101,*) NHYD2,UP,' ! i:',I2
	        IF ( IABS(KFWTYP) .EQ. 1 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              WRITE(101,*) ZD(ID), QFAKT*YS1FW(ID)/OMEG  ! r1=ys1/omeg
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              WRITE(101,*) ZD(ID), QFAKT*YS1FW(ID)/OMEG  ! r1=ys1/omeg
		    ELSE
	              WRITE(101,*) ZD(ID), QFAKT*(U*YS4FW(ID)/RHO(ID))/OMEG
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 2 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101,*) ZD(ID), QFAKT*YS2FW(ID)/OMEG  ! r2 = ys2/omeg
	          END DO                                       ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 3 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              WRITE(101,*) ZD(ID), QFAKT*YS3FW(ID)  ! r3 = ys3
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              WRITE(101,*) ZD(ID), QFAKT*YS3FW(ID)  ! r3 = ys3
	            ELSE
	              WRITE(101,*) ZD(ID), (0.0D0,0.0D0) !r3 zero in fluid
	            END IF
	          END DO
	        ELSE IF ( IABS(KFWTYP) .EQ. 4 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101,*) ZD(ID), QFAKT*YS4FW(ID)  !r4 = ys4 "without cslask"
	          END DO
	        ELSE IF ( IABS(KFWTYP) .EQ. 5 ) THEN !se "12/8 19" b1-2,5 (tau_rr part)
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              CSUM = ( 1 - 2*UA2(ID)/UB2(ID) ) * YS4FW(ID)  ! r4 = ys4
	              CSLASK = ( 1 - UA2(ID)/UB2(ID) ) * 4*RHO(ID)/UB2(ID)
	              CSUM = CSUM + CSLASK*(RK0*UP)*YS1FW(ID)/OMEG  ! r1=ys1/omeg
	              WRITE(101,*) ZD(ID), QFAKT*CSUM                     ! "r5"
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              ISLASK = MAX0(ID,NSW+1)
	              CSUM = ( 1 - 2*UA2(ISLASK)/UB2(ISLASK) ) * YS4FW(ID) !r4=ys4
	              CSLASK = ( 1 - UA2(ISLASK)/UB2(ISLASK) ) * 
     &	               4*RHO(ISLASK)/UB2(ISLASK)
	              CSUM = CSUM + CSLASK*(RK0*UP)*YS1FW(ID)/OMEG  ! r1=ys1/omeg
	              WRITE(101,*) ZD(ID), QFAKT*CSUM                     ! "r5"
	            ELSE
	              WRITE(101,*) ZD(ID), QFAKT*YS4FW(ID)                ! r4 = ys4
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 6 ) THEN !se "12/8 19" b1-2,5 (tau_rr part)
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              CSLASK = -2*RHO(ID)/UB2(ID)
	              WRITE(101,*) ZD(ID), QFAKT*CSLASK*YS1FW(ID)/OMEG !r1=ys1/omeg
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              ISLASK = MAX0(ID,NSW+1)
	              CSLASK = -2*RHO(ISLASK)/UB2(ISLASK)
	              WRITE(101,*) ZD(ID), QFAKT*CSLASK*YS1FW(ID)/OMEG !r1=ys1/omeg
	            ELSE
	              WRITE(101,*) ZD(ID), (0.0D0,0.0D0)
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 7 ) THEN !se "12/8 19" b1+baks,5 (bulk stress)
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              CSLASK = (RK0*UP) * 2*RHO(ID)/UB2(ID)
	              CSUM = YS4FW(ID) + CSLASK*YS1FW(ID)/OMEG !r4=ys4,r1=ys1/omeg
	              CSLASK = 1 - 4*UA2(ID)/(3*UB2(ID))
	              WRITE(101,*) ZD(ID), QFAKT*CSLASK*CSUM         ! "r7"
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              ISLASK = MAX0(ID,NSW+1)                  !r4=ys4,r1=ys1/omeg
	              CSLASK = (RK0*UP) * 2*RHO(ISLASK)/UB2(ISLASK)
	              CSUM = YS4FW(ISLASK) + CSLASK*YS1FW(ISLASK)/OMEG
	              CSLASK = 1 - 4*UA2(ISLASK)/(3*UB2(ISLASK))
	              WRITE(101,*) ZD(ID), QFAKT*CSLASK*CSUM         ! "r7"
	            ELSE
	              WRITE(101,*) ZD(ID), QFAKT*YS4FW(ID)           ! r4 = ys4
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE
	          STOP '%%% iuttyp_-3,-4, kfwtyp'
	        END IF
	      ELSE
	        WRITE(101) I, CSLASK/QFAKT, ' ! i'
	        WRITE(101) NHYD,UP,' ! i:',I
	        IF ( IABS(KFWTYP) .EQ. 1 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              WRITE(101) ZD(ID), QFAKT*YS1FW(ID)/OMEG  ! r1=ys1/omeg
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              WRITE(101) ZD(ID), QFAKT*YS1FW(ID)/OMEG  ! r1=ys1/omeg
		    ELSE
	              WRITE(101) ZD(ID), QFAKT*(U*YS4FW(ID)/RHO(ID))/OMEG
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 2 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101) ZD(ID), QFAKT*YS2FW(ID)/OMEG  ! r2 = ys2/omeg
	          END DO                                     ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 3 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              WRITE(101) ZD(ID), QFAKT*YS3FW(ID)  ! r3 = ys3
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              WRITE(101) ZD(ID), QFAKT*YS3FW(ID)  ! r3 = ys3
	            ELSE
	              WRITE(101) ZD(ID), (0.0D0,0.0D0) !r3 zero in fluid
	            END IF
	          END DO
	        ELSE IF ( IABS(KFWTYP) .EQ. 4 ) THEN
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            WRITE(101) ZD(ID), QFAKT*YS4FW(ID)  !r4 = ys4 "without cslask"
	          END DO
	        ELSE IF ( IABS(KFWTYP) .EQ. 5 ) THEN !se "12/8 19" b1-2,5 (tau_rr part)
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              CSUM = ( 1 - 2*UA2(ID)/UB2(ID) ) * YS4FW(ID)  ! r4 = ys4
	              CSLASK = ( 1 - UA2(ID)/UB2(ID) ) * 4*RHO(ID)/UB2(ID)
	              CSUM = CSUM + CSLASK*(RK0*UP)*YS1FW(ID)/OMEG  ! r1=ys1/omeg
	              WRITE(101) ZD(ID), QFAKT*CSUM                     ! "r5"
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              ISLASK = MAX0(ID,NSW+1)
	              CSUM = ( 1 - 2*UA2(ISLASK)/UB2(ISLASK) ) * YS4FW(ID) !r4=ys4
	              CSLASK = ( 1 - UA2(ISLASK)/UB2(ISLASK) ) * 
     &	               4*RHO(ISLASK)/UB2(ISLASK)
	              CSUM = CSUM + CSLASK*(RK0*UP)*YS1FW(ID)/OMEG  ! r1=ys1/omeg
	              WRITE(101) ZD(ID), QFAKT*CSUM                     ! "r5"
	            ELSE
	              WRITE(101) ZD(ID), QFAKT*YS4FW(ID)                ! r4 = ys4
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 6 ) THEN !se "12/8 19" b1-2,5 (tau_rr part)
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              CSLASK = -2*RHO(ID)/UB2(ID)
	              WRITE(101) ZD(ID), QFAKT*CSLASK*YS1FW(ID)/OMEG !r1=ys1/omeg
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              ISLASK = MAX0(ID,NSW+1)
	              CSLASK = -2*RHO(ISLASK)/UB2(ISLASK)
	              WRITE(101) ZD(ID), QFAKT*CSLASK*YS1FW(ID)/OMEG !r1=ys1/omeg
	            ELSE
	              WRITE(101) ZD(ID), (0.0D0,0.0D0)
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE IF ( IABS(KFWTYP) .EQ. 7 ) THEN !se "12/8 19" b1+baks,5 (bulk stress)
	          DO IH = 1, NHYD
	            ID = IDHYD(IH)
	            IF ( ID.LE.NSS .AND. NSS.GT.0 ) THEN
	              CSLASK = (RK0*UP) * 2*RHO(ID)/UB2(ID)
	              CSUM = YS4FW(ID) + CSLASK*YS1FW(ID)/OMEG !r4=ys4,r1=ys1/omeg
	              CSLASK = 1 - 4*UA2(ID)/(3*UB2(ID))
	              WRITE(101) ZD(ID), QFAKT*CSLASK*CSUM         ! "r7"
	            ELSE IF ( ID.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	              ISLASK = MAX0(ID,NSW+1)                  !r4=ys4,r1=ys1/omeg
	              CSLASK = (RK0*UP) * 2*RHO(ISLASK)/UB2(ISLASK)
	              CSUM = YS4FW(ISLASK) + CSLASK*YS1FW(ISLASK)/OMEG
	              CSLASK = 1 - 4*UA2(ISLASK)/(3*UB2(ISLASK))
	              WRITE(101) ZD(ID), QFAKT*CSLASK*CSUM         ! "r7"
	            ELSE
	              WRITE(101) ZD(ID), QFAKT*YS4FW(ID)           ! r4 = ys4
	            END IF
	          END DO                                  ! "without cslask"
	        ELSE
	          STOP '%%% iuttyp_-3,-4, kfwtyp'
	        END IF
	      END IF  ! iuttyp
C	      ---
	      ELSE
C	      ---
	      IF (I.EQ.1) THEN
	        OPEN(101,FILE='QFIL',STATUS='UNKNOWN')
		WRITE(101,*) LL
	      END IF
	      IF ( IDMAX.LT.0 .OR. ICMAX.EQ.0 ) STOP '%%% vanishing mode'
	      IF ( IDMAX.LE.NSS .AND. NSS.GT.0 ) THEN
	        ISLASK = -ICMAX
	      ELSE IF ( IDMAX.GE.NSW .AND. (NSW+1).LT.NOL ) THEN
	        ISLASK = -ICMAX
	      ELSE
	        ISLASK = ICMAX
	      END IF
	      WRITE(101,*) ZD(IDMAX), ISLASK
C	      ---
	      END IF  ! kfwtyp

	      IF (KFWTYP.GE.0) THEN
	      I2 = I
	      DO ID = 1, NOL-1
	        IF ( IFSMTYP(ID) .GT. -4 ) STOP '%%% iuttyp_-3,-4, ifsmtyp()'
	      END DO
	      IF (IUTTYP.EQ.-3) THEN
	        IF (I.EQ.1) THEN 
	          WRITE(102,*) 
	          WRITE(102,*) OMEG,' ! omeg'
	          WRITE(102,*) U0,' ! u0'
	          WRITE(102,*) RK0,' ! rk0'
	          WRITE(102,*) ZD(0),' ! zd(0)'
	          IF ( RHO(NOL) .LE. 0.0D0 ) THEN
	            WRITE(102,*) NSS,NSW,NOL,' ! nss,nsw,nol'
	            ISLASK = NOL-1
	          ELSE
	            WRITE(102,*) NSS,NSW,-NOL,' ! nss,nsw,-nol'  ! half-space indicator
	            ISLASK = NOL
	            IF ( ISOLID(NOL) .EQ. 0 ) UB2(NOL) = (0.0D0,0.0D0) !fluid indicator
	          END IF
	          DO ID = 1, NSS
	            WRITE(102,*) ZD(ID), RHO(ID)
	            WRITE(102,*) UA2(ID), UB2(ID)
	          END DO
	          DO ID = NSS+1, NSW
	            WRITE(102,*) ZD(ID), RHO(ID)
	            WRITE(102,*) UA2(ID)
	          END DO
	          DO ID = NSW+1, ISLASK
	            WRITE(102,*) ZD(ID), RHO(ID)
	            WRITE(102,*) UA2(ID), UB2(ID)
	          END DO
	          WRITE(102,*) 
	          WRITE(102,'(3(X,I7),X,G24.16,X,I7,A)') 
     $	           LL,ICUTTYP,ISHEETTYP,ALA,IUROT,' !  ll, ...'
	        END IF
	        WRITE(102,*) 
	        WRITE(102,*) I2,UP,' ! i,up'
	        WRITE(102,*) DLNDU,DERIVDU
	        IF (NSS.GT.0) THEN
	          DO ID = 0, NSS
	            WRITE(102,*) YS1FW(ID)
	            WRITE(102,*) YS2FW(ID)
	            WRITE(102,*) YS3FW(ID)
	            WRITE(102,*) YS4FW(ID)
	          END DO
	          IF (NSW.EQ.NSS) STOP '%%% 0<nss=nsw'
	        ELSE IF (NSW.GT.0) THEN
	          WRITE(102,*) YS2FW(0)
	          WRITE(102,*) YS4FW(0)
	        END IF
	        DO ID = NSS+1, NSW-1
	          WRITE(102,*) YS2FW(ID)
	          WRITE(102,*) YS4FW(ID)
	        END DO
	        IF ( (NSW+1) .LT. NOL ) THEN
	          DO ID = NSW, NOL-1
	            WRITE(102,*) YS1FW(ID)
	            WRITE(102,*) YS2FW(ID)
	            WRITE(102,*) YS3FW(ID)
	            WRITE(102,*) YS4FW(ID)
	          END DO
	        ELSE
	          IF (NSW.EQ.0) STOP '%%% rmodpfw, nsw=0 & nol<=1'
	          WRITE(102,*) YS2FW(NSW)
	          WRITE(102,*) YS4FW(NSW)
	        END IF
	      ELSE
	        IF (I.EQ.1) THEN 
	          WRITE(102) OMEG           ! omeg
	          WRITE(102) U0             ! u0
	          WRITE(102) RK0            ! rk0
	          WRITE(102) ZD(0)          ! zd(0)
	          IF ( RHO(NOL) .LE. 0.0D0 ) THEN
	            WRITE(102) NSS,NSW,NOL  ! nss,nsw,nol
	            ISLASK = NOL-1
	          ELSE
	            WRITE(102) NSS,NSW,-NOL  ! nss,nsw,-nol  !!!half-space indic
	            ISLASK = NOL
	            IF ( ISOLID(NOL) .EQ. 0 ) UB2(NOL) = (0.0D0,0.0D0) !fluid indic
	          END IF
	          DO ID = 1, NSS
	            WRITE(102) ZD(ID), RHO(ID)
	            WRITE(102) UA2(ID), UB2(ID)
	          END DO
	          DO ID = NSS+1, NSW
	            WRITE(102) ZD(ID), RHO(ID)
	            WRITE(102) UA2(ID)
	          END DO
	          DO ID = NSW+1, ISLASK
	            WRITE(102) ZD(ID), RHO(ID)
	            WRITE(102) UA2(ID), UB2(ID)
	          END DO
	          WRITE(102) LL,ICUTTYP,ISHEETTYP,ALA,IUROT  !  ll, ...
	        END IF
	        WRITE(102) I,UP  ! i,up  !!! avoid "integer*2" i2
	        WRITE(102) DLNDU,DERIVDU
	        IF (NSS.GT.0) THEN
	          DO ID = 0, NSS
	            WRITE(102) YS1FW(ID)
	            WRITE(102) YS2FW(ID)
	            WRITE(102) YS3FW(ID)
	            WRITE(102) YS4FW(ID)
	          END DO
	          IF (NSW.EQ.NSS) STOP '%%% 0<nss=nsw'
	        ELSE IF (NSW.GT.0) THEN
	          WRITE(102) YS2FW(0)
	          WRITE(102) YS4FW(0)
	        END IF
	        DO ID = NSS+1, NSW-1
	          WRITE(102) YS2FW(ID)
	          WRITE(102) YS4FW(ID)
	        END DO
	        IF ( (NSW+1) .LT. NOL ) THEN
	          DO ID = NSW, NOL-1
	            WRITE(102) YS1FW(ID)
	            WRITE(102) YS2FW(ID)
	            WRITE(102) YS3FW(ID)
	            WRITE(102) YS4FW(ID)
	          END DO
	        ELSE
	          IF (NSW.EQ.0) STOP '%%% rmodpfw, nsw=0 & nol<=1'
	          WRITE(102) YS2FW(NSW)
	          WRITE(102) YS4FW(NSW)
	        END IF
	      END IF  ! iuttyp
	      END IF  ! kfwtyp

	    END IF  ! iuttyp.eq.-3 .or. iuttyp-eq.-4

	    IF ( DERIVDU .NE. (0.0D0,0.0D0) ) THEN
	      DUDO = - DEXP(DLNDO-DLNDU) * DERIVDO/DERIVDU
	    ELSE
	      IF ( RHO(NSS).GE.0.0D0 .OR. RHO(NSS).EQ.-99.0D0 ) STOP '%derivdu,rho'
	      IF ( RHO(NSW+1).GE.0.0D0 .OR. RHO(NSW+1).EQ.-99.0D0 ) STOP '%derivdu,rho'
	      IF ( ZD(NSS) .NE. ZD(NSW) ) STOP '%derivdu,zd'
	      IF ( DERIVDO .NE. (0.0D0,0.0D0) ) STOP '%derivdu,derivdo'
	      DUDO = (0.0D0,0.0D0)  ! avoid "NaN" for "trivial one-mode medium"
	    END IF
	    WRITE(*,*) '  du/df = ', 2*PI * DUDO
	    WRITE(*,*) '  group-u = ', U + OMEG*DUDO
	    UGRARR(I) =  U + OMEG*DUDO

	    ISLASK = 0
	    DCHINUREMAX = 0.0D0
	    DCHINUIMMAX = 0.0D0
	    DO ID = 1, NOL-1
	      SLASK1 = DZ(ID) * DABS(DIMAG(URA_FW(ID)))
	      IF (SLASK1.GT.DCHINUREMAX) THEN
	        ISLASK = ID
	        DCHINUREMAX = SLASK1
	      END IF
	      SLASK2 = DZ(ID) * DABS(DREAL(URA_FW(ID)))
	      IF (SLASK2.GT.DCHINUIMMAX) THEN
	        ISLASK = ID
	        DCHINUIMMAX = SLASK2
	      END IF
	    END DO
	    DO ID = 1, NSS
	      SLASK1 = DZ(ID) * DABS(DIMAG(URB_FW(ID)))
	      IF (SLASK1.GT.DCHINUREMAX) THEN
	        ISLASK = ID
	        DCHINUREMAX = SLASK1
	      END IF
	      SLASK2 = DZ(ID) * DABS(DREAL(URB_FW(ID)))
	      IF (SLASK2.GT.DCHINUIMMAX) THEN
	        ISLASK = ID
	        DCHINUIMMAX = SLASK2
	      END IF
	    END DO
	    DO ID = NSW+1, NOL-1
	      SLASK1 = DZ(ID) * DABS(DIMAG(URB_FW(ID)))
	      IF (SLASK1.GT.DCHINUREMAX) THEN
	        ISLASK = ID
	        DCHINUREMAX = SLASK1
	      END IF
	      SLASK2 = DZ(ID) * DABS(DREAL(URB_FW(ID)))
	      IF (SLASK2.GT.DCHINUIMMAX) THEN
	        ISLASK = ID
	        DCHINUIMMAX = SLASK2
	      END IF
	    END DO
	    WRITE(*,*) '  dchinumax, Re & Im ; id :', 
     $	    SNGL(DCHINUREMAX),SNGL(DCHINUIMMAX), ISLASK

	  END DO  ! i=1,ll
C	  ------------------------------------------------------------------------------------------
	  IF ( IUTTYP.NE.-2 .AND. IUTTYP.NE.-4 ) THEN
	    DO II = IIA, IIB
	      WRITE(II,*)
	      WRITE(II,*) LL, ' ! group-u for the modes'
	      DO I = 1, LL
	        WRITE(II,*) UGRARR(I)
	      END DO     
	    END DO     
	  ELSE
	    DO II = IIA, IIB
	      WRITE(II) LL, ' ! group-u for the modes'
	      DO I = 1, LL
	        WRITE(II) UGRARR(I)
	      END DO     
	    END DO     
	  END IF
	  DEALLOCATE ( UGRARR )
	  
	  IF (FNORMAL) THEN
	    DO II = IIA, IIB
	      IF ( II.LE.4 .OR. II.GE.100 ) THEN
	        IIW = II
	      ELSE
	        IIW = II+2
	      END IF
	      CLOSE(IIW)
	    END DO
	  END IF

 99	  CONTINUE

	END DO  ! jf=1,nfr

	STOP
	END
