
C	additional subroutines (and functions) for "Full Wave" computations 
C	special "fw" versions :  liqsifw  solsifw  liqsixfw  solsixfw

C	Sven Ivansson  1992






	SUBROUTINE LIQSIFW ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	TLNFW, YT1FW,YT3FW, Z1FW, DIS2FW,DIS4FW )
C	special version for "full wave"
C	in particular: "parameter change"
C	additional "fw input & output" in "common /fwparam,fwpart,fwys/"
C	propagation of "ys2fw,ys4fw"
C	stabilization with "tlnfw, yt1fw,yt3fw, z1fw"
C	"ngsifw and zliqsifw" must have been called in advance (OBS)

	USE INMODELa
	USE FWYSa
	USE FWPARTa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C	nol <= noldim

	INCLUDE '../DPCFSUB.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE 'FWPARAM.INC'

	INTEGER NL1,NL2,NRIKT
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2
	REAL*8 TLNFW(*)
	COMPLEX*16 YT1FW(*), YT3FW(*), Z1FW(*), DIS2FW(*), DIS4FW(*)

	INTEGER ID, IDQ, IYSTABLOK, IND, ICLOK
	REAL*8 PLN, PLNREF, 
     $	ZENDFSUB, Z0FSUB, SLASK, SLASK1,SLASK3, PI
	COMPLEX*16 YP1,YP2, YT1,YT3, Z1, URA, URANEGI, 
     $	CP,CPX, STP,SDP, WWLOK, 
     $	CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	YP(2)
	LOGICAL EXPTYP, INHOMOG, STABKVAR

	EXTERNAL FSUBLIQ2

	DATA PI /3.141592653589793D0/

	INHOMOG = .TRUE.
	GOTO 100

	ENTRY LIQSIFW0 ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	TLNFW, YT1FW,YT3FW )

	INHOMOG = .FALSE.

 100	CONTINUE

C	"  upwards when  nrikt = -1  ,  downwards when  nrikt = +1  "
	IF (NRIKT.LT.0) THEN
	  IDQ = NL1
	ELSE
	  IDQ = NL1-1
	END IF
	PLN = 0.0D0
	YP1 = YS2FW(IDQ)
	YP2 = YS4FW(IDQ)

	DO ID = NL1, NL2, NRIKT

	  IF (IYSTAB.NE.0) THEN
	    STABKVAR = .TRUE.
	    YT1 = YT1FW(ID)
	    YT3 = YT3FW(ID)
	    IF (INHOMOG) Z1 = Z1FW(ID)
	  ELSE
	    STABKVAR = .FALSE.
	  END IF

	  IF (INHOMOG) THEN
	    IF (NRIKT.LT.0) THEN
	      YP1 = YP1 - DIS2FW(ID)
	      YP2 = YP2 - DIS4FW(ID)
	    END IF
	  END IF

	  IF (IFSMTYP(ID).LE.-4) THEN

	    IF ( IFSMTYP(ID).NE.-6 .AND. IFSMTYP(ID).NE.(NRIKT-6) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
C	      stored "ura" (freq-indep) and "call liqpart" (freq-dep)
	      URA = URA_FW(ID)
	      EXPTYP = EXPTYP_FW(ID)
	      PLNREF = TLNREF_FW(ID)
	      CP = CP_FW(ID)
	      STP = STP_FW(ID)
	      SDP = SDP_FW(ID)
	      IF (EXPTYP) URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura
	    END IF
	    YP2 = YP2 / RHO(ID)   ! already here   !"1/rho(id)" could be stored
	    IF (EXPTYP) THEN
C	      now the first propagation step
	      WWLOK = YP1 / URANEGI
	      YP1 = WWLOK + YP2
	      YP2 = -WWLOK + YP2
C	      now the main propagation step, freq-dependent
	      IYSTABLOK = IYSTAB
	      IF (IYSTABLOK.GT.0) THEN
		IF (NRIKT.GT.0) THEN
	          IF ( DMAX1(DABS(DREAL(YT3)),DABS(DIMAG(YT3))) .LE.
     $	            XSTABEPS*DMAX1(DABS(DREAL(YT1)),DABS(DIMAG(YT1))) ) THEN
	            WRITE(*,*) '%%% liqsifw, depths ',ZD(ID-1),ZD(ID)
	            WRITE(*,*) 'u ',U,' modif stab replaced %%%'
	            IYSTABLOK = -1
	            IF (INHOMOG) THEN
	              CALL UNDOXLIQ (ID, YT1,YT3, NRIKT, Z1)
	            ELSE
	              CALL UNDOXLIQ0 (ID, YT1,YT3)
	            END IF
	          END IF
	        ELSE
	          IF ( DMAX1(DABS(DREAL(YT1)),DABS(DIMAG(YT1))) .LE.
     $	            XSTABEPS*DMAX1(DABS(DREAL(YT3)),DABS(DIMAG(YT3))) ) THEN
	            WRITE(*,*) '%%% liqsifw, depths ',ZD(ID-1),ZD(ID)
	            WRITE(*,*) 'u ',U,' modif stab replaced %%%'
	            IYSTABLOK = -1
	            IF (INHOMOG) THEN
	              CALL UNDOXLIQ (ID, YT1,YT3, NRIKT, Z1)
	            ELSE
	              CALL UNDOXLIQ0 (ID, YT1,YT3)
	            END IF
	          END IF
	        END IF
	      END IF
	      IF (IYSTABLOK.LE.0) THEN
	        IF (IYSTABLOK.NE.0) THEN
	          IF (PLNREF.GT.XSPLITLIM) THEN
	            WRITE(*,*) '%%% liqsifw, depths ',ZD(ID-1),ZD(ID)
	            WRITE(*,*) 'quotient u ',PLNREF/XSPLITLIM,U
	            WRITE(*,*) 'split needed, automatic split not ',
     $	            'implemented %%%'
	            STOP 3141
	          END IF
	        END IF
	        PLN = PLN + PLNREF - DLOG2       ! -dlog2 from "first step"
	        IF (NRIKT.GT.0) THEN
	          YP1 = DREAL(SDP)*CP * YP1
	          YP2 = DCONJG(CP) * YP2
	        ELSE
	          YP1 = DCONJG(CP) * YP1
	          YP2 = DREAL(SDP)*CP * YP2
	        END IF
	      ELSE
C	        stabilization included now, using ytfw (and zfw if inhomog)
	        PLN = PLN - PLNREF - DLOG2       ! -dlog2 from "first step"
	        IF (NRIKT.GT.0) THEN
	          YP1 = CP * YP1
	          IF (PLN.NE.0.0D0) THEN
	            IF (YP1.NE.(0.0D0,0.0D0)) THEN
	              SLASK = DEXP(PLN)
	              YP1 = SLASK * YP1
	            END IF
	            PLN = 0.0D0
	          END IF
	          IF (INHOMOG) THEN
		    YP2 = ( Z1 - YT1*YP1 ) / YT3
		  ELSE
		    YP2 = -YT1*YP1 / YT3
	          END IF
	        ELSE
	          YP2 = CP * YP2
	          IF (PLN.NE.0.0D0) THEN
	            IF (YP2.NE.(0.0D0,0.0D0)) THEN
	              SLASK = DEXP(PLN)
	              YP2 = SLASK * YP2
	            END IF
	            PLN = 0.0D0
	          END IF
	          IF (INHOMOG) THEN
	            YP1 = ( Z1 - YT3*YP2 ) / YT1
	          ELSE
	            YP1 = -YT3*YP2 / YT1
	          END IF
	        END IF
	        STABKVAR = .FALSE.
	      END IF
C	      now the next propagation step
	      WWLOK = YP1 - YP2
	      YP2 = RHO(ID) * ( YP1 + YP2 )
	      YP1 = URANEGI*WWLOK
	    ELSE
C	      now the main and only propagation step, freq-dependent
	      IF (IYSTAB.NE.0) THEN
	        IF (PLNREF.GT.XSPLITLIM) THEN
	          WRITE(*,*) '%%% liqsifw, depths ',ZD(ID-1),ZD(ID)
	          WRITE(*,*) 'quotient u ',PLNREF/XSPLITLIM,U
	          WRITE(*,*) 'split needed, automatic split not implemented %%%'
	          STOP 3141
		END IF
	      END IF
	      PLN = PLN + PLNREF
	      IF (NRIKT.GT.0) THEN
	        WWLOK = CP*YP1 + STP*YP2
	        YP2 = RHO(ID) * ( -SDP*YP1 + CP*YP2 )
	        YP1 = WWLOK
	      ELSE
	        WWLOK = CP*YP1 - STP*YP2
	        YP2 = RHO(ID) * ( SDP*YP1 + CP*YP2 )
	        YP1 = WWLOK
	      END IF
	    END IF

	  ELSE IF (IFSMTYP(ID).LE.0) THEN

C	    stored "call airypart"
	    CFAKT = CFAKT_FW(ID)
	    BCFAKT = BCFAKT_FW(ID)
	    BCFAKTINV = BCFAKTINV_FW(ID)
	    YA = YA_FW(ID)
	    YB = YB_FW(ID)
	    PLNREF = TLNREF_FW(ID)
	    CP = CP_FW(ID)
	    CPX = CPX_FW(ID)
	    STP = STP_FW(ID)
	    SDP = SDP_FW(ID)
C	    now propagation
	    YP2 = YP2 / RHO(ID)   ! already here   !"1/rho(id)" could be stored
	    IF (IYSTAB.NE.0) THEN
	      IF (PLNREF.GT.XSPLITLIM) THEN
	        WRITE(*,*) '%%% liqsifw, depths ',ZD(ID-1),ZD(ID)
	        WRITE(*,*) 'quotient ',PLNREF/XSPLITLIM
	        WRITE(*,*) 'split needed, automatic split not implemented %%%'
	        STOP 3142
	      END IF
	    END IF
	    PLN = PLN + PLNREF
	    IF (NRIKT.GT.0) THEN
	      WWLOK = PI * ( CP*YP1 + STP*YP2 )
	      YP2 = (PI*RHO(ID)) * ( -SDP*YP1 + CPX*YP2 )  !"pi*rho" could stor
	      YP1 = WWLOK
	    ELSE
	      WWLOK = PI * ( CPX*YP1 - STP*YP2 )
	      YP2 = (PI*RHO(ID)) * ( SDP*YP1 + CP*YP2 )  !"pi*rho" could stored
	      YP1 = WWLOK
	    END IF

	  ELSE

	    IF (IYSTAB.NE.0) THEN
	      WRITE(*,*) '%%% liqsifw, depths ',ZD(ID-1),ZD(ID)
	      WRITE(*,*) 'stabilization not implemented'
	      STOP 3141
	    END IF
	    IFSMTYPLOK = IFSMTYP(ID)
	    IF (NRIKT.LT.0) THEN
	      ZENDFSUB = ZD(ID-1)
	      Z0FSUB = ZD(ID)
	    ELSE
	      ZENDFSUB = ZD(ID)
	      Z0FSUB = ZD(ID-1)
	    END IF
C	    now propagation
	    YP(1) = YP1
	    YP(2) = YP2
	    CALL DPCLIN ( FSUBLIQ2,2, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	    HFAKTMAX,HOFF,EPSMAGNIFMAX, ZENDFSUB, Z0FSUB,PLN,
     $	    YP,EPSFSUB,HSTEPFSUB, IND )
	    IF (IND.NE.0) THEN
	      WRITE(*,*) ' %% u,id ',UCOP,ID
	      STOP 6770
	    END IF
	    YP1 = YP(1)
	    YP2 = YP(2)

	  END IF

	  IF (PLN.NE.0.0D0) THEN
	    IF ( YP1.NE.(0.0D0,0.0D0) .OR. YP2.NE.(0.0D0,0.0D0) ) THEN
	      SLASK = DEXP(PLN)
	      YP1 = SLASK * YP1
	      YP2 = SLASK * YP2
	    END IF
	    PLN = 0.0D0
	  END IF
	  IF (INHOMOG) THEN
	    IF (NRIKT.GT.0) THEN
	      YP1 = YP1 + DIS2FW(ID)
	      YP2 = YP2 + DIS4FW(ID)
	    END IF
	    IF (STABKVAR) THEN
C	      stabilization now, using ytfw and zfw
	      SLASK1 = DMAX1 ( DABS(DREAL(YT1)), DABS(DIMAG(YT1)) )
	      SLASK3 = DMAX1 ( DABS(DREAL(YT3)), DABS(DIMAG(YT3)) )
	      ICLOK = 1
	      SLASK = SLASK1
	      IF (SLASK3.GT.SLASK) ICLOK = 3
	      GOTO (1,92,3) ICLOK
 92	      STOP 92
 1	      YP1 =  ( Z1 - YT3*YP2 ) / YT1
	      GOTO 9
 3	      YP2 = ( Z1 - YT1*YP1 ) / YT3
 9	      CONTINUE
	    END IF
	  ELSE
	    IF (STABKVAR) THEN
C	      "zero" stabilization now, using ytfw
	      SLASK1 = DMAX1 ( DABS(DREAL(YT1)), DABS(DIMAG(YT1)) )
	      SLASK3 = DMAX1 ( DABS(DREAL(YT3)), DABS(DIMAG(YT3)) )
	      ICLOK = 1
	      SLASK = SLASK1
	      IF (SLASK3.GT.SLASK) ICLOK = 3
	      GOTO (10,920,30) ICLOK
 920	      STOP 920
 10	      YP1 = -YT3*YP2 / YT1
	      GOTO 90
 30	      YP2 = -YT1*YP1 / YT3
 90	      CONTINUE
	    END IF
	  END IF
	  IF (NRIKT.LT.0) THEN
	    IDQ = ID-1
	  ELSE
	    IDQ = ID
	  END IF
	  YS2FW(IDQ) = YP1
	  YS4FW(IDQ) = YP2

	END DO

	RETURN
	END









	SUBROUTINE SOLSIFW ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, 
     $	Z1FW,Z2FW,Z3FW,Z4FW, DIS1FW,DIS2FW,DIS3FW,DIS4FW )
C	special version for "full wave"
C	in particular: "parameter change"
C	additional "fw input & output" in "common /fwparam,fwpart,fwys/"
C	propagation of "ys1fw,ys2fw,ys3fw,ys4fw"
C	stabilization with "tlnfw, yt1fw,..., z1fw,z2fw,z3fw,z4fw"
C	"harksifw and zsolsifw" must have been called in advance (OBS)

	USE INMODELa
	USE FWYSa
	USE FWPARTa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C	nol <= noldim

	INCLUDE '../DPCFSUB.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE 'FWPARAM.INC'

	INTEGER NL1,NL2,NRIKT
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2
	REAL*8 TLNFW(*)
	COMPLEX*16 YT1FW(*), YT2FW(*), YT3FW(*), YT4FW(*), YT5FW(*),
     $	YT6FW(*), Z1FW(*), Z2FW(*), Z3FW(*), Z4FW(*), 
     $	DIS1FW(*), DIS2FW(*), DIS3FW(*), DIS4FW(*)

	INTEGER ID, IDQ, IYSTABLOK, IND, ICLOK
	REAL*8 SLN, ZENDFSUB, Z0FSUB, SLNREF,E1X, TLNREF, SLASK, 
     $	SLASK1,SLASK2,SLASK3,SLASK4,SLASK6
	COMPLEX*16 YS1,YS2,YS3,YS4, YT1,YT2,YT3,YT4,YT5,YT6, Z1,Z2,Z3,Z4,
     $  G0, G0U, URA,URB, URANEGI,URBNEGI, CP,CQ, STP,STQ,SDP,SDQ, 
     $	S1Q,S2Q, CSLASK, 
     $	YS(4)
	LOGICAL EXPTYP, INHOMOG, STABKVAR

	EXTERNAL FSUBSOL4

	INHOMOG = .TRUE.
	GOTO 100

	ENTRY SOLSIFW0 ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW )

	INHOMOG = .FALSE.

 100	CONTINUE

C	"  upwards when  nrikt = -1  ,  downwards when  nrikt = +1  "
	IF (NRIKT.LT.0) THEN
	  IDQ = NL1
	ELSE
	  IDQ = NL1-1
	END IF
	SLN = 0.0D0
	YS1 = YS1FW(IDQ)
	YS2 = YS2FW(IDQ)
	YS3 = YS3FW(IDQ)
	YS4 = YS4FW(IDQ)

	DO ID = NL1, NL2, NRIKT

	  IF (IYSTAB.NE.0) THEN
	    STABKVAR = .TRUE.
	    YT1 = YT1FW(ID)
	    YT2 = YT2FW(ID)
	    YT3 = YT3FW(ID)
	    YT4 = YT4FW(ID)
	    YT5 = YT5FW(ID)
	    YT6 = YT6FW(ID)
	    IF (INHOMOG) THEN
	      Z1 = Z1FW(ID)
	      Z2 = Z2FW(ID)
	      Z3 = Z3FW(ID)
	      Z4 = Z4FW(ID)
	    END IF
	  ELSE
	    STABKVAR = .FALSE.
	  END IF

	  IF (INHOMOG) THEN
	    IF (NRIKT.LT.0) THEN
	      YS1 = YS1 - DIS1FW(ID)
	      YS2 = YS2 - DIS2FW(ID)
	      YS3 = YS3 - DIS3FW(ID)
	      YS4 = YS4 - DIS4FW(ID)
	    END IF
	  END IF

	  IF (IFSMTYP(ID).LE.-4) THEN

	    IF ( IFSMTYP(ID).NE.-6 .AND. IFSMTYP(ID).NE.(NRIKT-6) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
	      G0 = TWOVB2(ID)                      ! local twovb2(id)
	      G0U = G0*U                  ! could be stored from "harksifw"
C	      stored "ura,urb" (freq-indep) and "call (h)solpart" (freq-dep)
	      URA = URA_FW(ID)
	      URB = URB_FW(ID)
	      EXPTYP = EXPTYP_FW(ID)
	      SLNREF = SLNREF_FW(ID)
	      E1X = E1X_FW(ID)
	      CP = CP_FW(ID)
	      CQ = CQ_FW(ID)
	      STP = STP_FW(ID)
	      STQ = STQ_FW(ID)
	      SDP = SDP_FW(ID)
	      SDQ = SDQ_FW(ID)
	      TLNREF = TLNREF_FW(ID)
	      IF (EXPTYP) THEN
	        URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	        URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb
	      END IF
	    END IF
C	    now the C propagation step
	    YS3 = YS3 / RHO(ID)                    !"1/rho(id)" could be stored
	    YS4 = YS4 / RHO(ID)                    !"1/rho(id)" could be stored
	    YS3 = G0U*YS2 + YS3
	    YS2 = YS2 - U*YS3
	    YS4 = G0U*YS1 + YS4
	    YS1 = YS1 - U*YS4
	    IF (EXPTYP) THEN
C	      now the D propagation step
	      SLN = SLN - DLOG2
	      S1Q = YS1/URBNEGI  ! "1/urbnegi" could be stored from "harksifw"
	      S2Q = YS2/URANEGI  ! "1/uranegi" could be stored from "harksifw"
	      YS1 = S1Q + YS3
	      YS2 = S2Q + YS4
	      YS3 = S1Q - YS3
	      YS4 = -S2Q + YS4
C	      now the PE propagation step, freq-dependent
	      IYSTABLOK = IYSTAB
	      IF (IYSTABLOK.GT.0) THEN
	        IF (NRIKT.GT.0) THEN
	          IF ( DMAX1(DABS(DREAL(YT6)),DABS(DIMAG(YT6))) .LE. 
     $		  XSTABEPS*DMAX1(DABS(DREAL(YT2)),DABS(DIMAG(YT2)),
     $	          DABS(DREAL(YT3)),DABS(DIMAG(YT3)),
     $	          DABS(DREAL(YT4)),DABS(DIMAG(YT4)),
     $	          DABS(DREAL(YT5)),DABS(DIMAG(YT5))) ) THEN
	            WRITE(*,*) '%%% solsifw, depths ',ZD(ID-1),ZD(ID)
	            WRITE(*,*) 'u ',U,' modif stab replaced %%%'
	            IYSTABLOK = -1
	            IF (INHOMOG) THEN
	              CALL UNDOXSOL (U, ID, YT1,YT2,YT3,YT4,YT6, 
     $	              NRIKT, Z1,Z2,Z3,Z4)
	            ELSE
	              CALL UNDOXSOL0 (U, ID, YT1,YT2,YT3,YT4,YT6 )
	            END IF
	          END IF
	        ELSE
	          IF ( DMAX1(DABS(DREAL(YT1)),DABS(DIMAG(YT1))) .LE. 
     $	          XSTABEPS*DMAX1(DABS(DREAL(YT2)),DABS(DIMAG(YT2)),
     $	          DABS(DREAL(YT3)),DABS(DIMAG(YT3)),
     $	          DABS(DREAL(YT4)),DABS(DIMAG(YT4)),
     $	          DABS(DREAL(YT5)),DABS(DIMAG(YT5))) ) THEN
	            WRITE(*,*) '%%% solsifw, depths ',ZD(ID-1),ZD(ID)
	            WRITE(*,*) 'u ',U,' modif stab replaced %%%'
	            IYSTABLOK = -1
	            IF (INHOMOG) THEN
	              CALL UNDOXSOL (U, ID, YT1,YT2,YT3,YT4,YT6, 
     $	              NRIKT, Z1,Z2,Z3,Z4)
	            ELSE
	              CALL UNDOXSOL0 (U, ID, YT1,YT2,YT3,YT4,YT6 )
	            END IF
	          END IF
	        END IF
	      END IF
	      IF (IYSTABLOK.LE.0) THEN
	        IF (IYSTABLOK.NE.0) THEN
	          IF (DABS(SLNREF).GT.XSPLITLIM) THEN
	            WRITE(*,*) '%%% solsifw, depths ',ZD(ID-1),ZD(ID)
	            WRITE(*,*) 'quotient u ',DABS(SLNREF)/XSPLITLIM,U
	            WRITE(*,*) 'split needed, automatic split not ',
     $	            'implemented %%%'
	            STOP 3141
	          END IF
	        END IF
	        SLN = SLN + DABS(SLNREF)
	        IF (NRIKT.GT.0) THEN
	          YS1 = DREAL(SDQ)*CQ * YS1
	          YS2 = DREAL(SDP)*CP * YS2
	          YS3 = DCONJG(CQ) * YS3
	          YS4 = DCONJG(CP) * YS4
	        ELSE
	          YS1 = DCONJG(CQ) * YS1
	          YS2 = DCONJG(CP) * YS2
	          YS3 = DREAL(SDQ)*CQ * YS3
	          YS4 = DREAL(SDP)*CP * YS4
	        END IF
	        IF (SLNREF.GE.0.0D0) THEN
	          YS1 = E1X*YS1
	          YS3 = E1X*YS3
	        ELSE
	          YS2 = E1X*YS2
	          YS4 = E1X*YS4
	        END IF
	      ELSE
C	        stabilization included now, using ytfw (and zfw if inhomog)
	        IF (SLNREF.GE.0.0D0) THEN
	          SLN = SLN - (TLNREF-SLNREF)
	        ELSE
	          SLN = SLN - (TLNREF+SLNREF)
	        END IF
	        IF (NRIKT.GT.0) THEN
	          YS1 = CQ * YS1
	          YS2 = CP * YS2
	          IF (SLNREF.GE.0.0D0) THEN
	            YS2 = E1X*YS2
	          ELSE
	            YS1 = E1X*YS1
	          END IF
	          IF (SLN.NE.0.0D0) THEN
	            IF ( YS1.NE.(0.0D0,0.0D0) .OR. YS2.NE.(0.0D0,0.0D0) ) THEN
	              SLASK = DEXP(SLN)
	              YS1 = SLASK * YS1
	              YS2 = SLASK * YS2
	            END IF
	            SLN = 0.0D0
	          END IF
	          CSLASK = 1/YT6
	          IF (INHOMOG) THEN
	            YS3 = CSLASK * ( -Z4 - YT3*YS1 - YT5*YS2 )
	            YS4 = CSLASK * ( Z3 + YT2*YS1 + YT4*YS2 )
	          ELSE
	            YS3 = CSLASK * ( -YT3*YS1 - YT5*YS2 )
	            YS4 = CSLASK * ( YT2*YS1 + YT4*YS2 )
	          END IF
	        ELSE
	          YS3 = CQ * YS3
	          YS4 = CP * YS4
	          IF (SLNREF.GE.0.0D0) THEN
	            YS4 = E1X*YS4
	          ELSE
	            YS3 = E1X*YS3
	          END IF
	          IF (SLN.NE.0.0D0) THEN
	            IF ( YS3.NE.(0.0D0,0.0D0) .OR. YS4.NE.(0.0D0,0.0D0) ) THEN
	              SLASK = DEXP(SLN)
	              YS3 = SLASK * YS3
	              YS4 = SLASK * YS4
	            END IF
	            SLN = 0.0D0
	          END IF
	          CSLASK = 1/YT1
	          IF (INHOMOG) THEN
	            YS1 = CSLASK * ( -Z2 + YT4*YS3 + YT5*YS4 )
	            YS2 = CSLASK * ( Z1 - YT2*YS3 - YT3*YS4 )
	          ELSE
	            YS1 = CSLASK * ( YT4*YS3 + YT5*YS4 )
	            YS2 = CSLASK * ( -YT2*YS3 - YT3*YS4 )
	          END IF
	        END IF
	        STABKVAR = .FALSE.
	      END IF
C	      now the Dinv propagation step
	      S1Q = YS1 + YS3
	      S2Q = YS2 - YS4
	      YS3 = YS1 - YS3
	      YS4 = YS2 + YS4
	      YS1 = URBNEGI*S1Q
	      YS2 = URANEGI*S2Q
	    ELSE
C	      now the PT propagation step, freq-dependent
	      IF (IYSTAB.NE.0) THEN
	        IF (DABS(SLNREF).GT.XSPLITLIM) THEN
	          WRITE(*,*) '%%% solsifw, depths ',ZD(ID-1),ZD(ID)
	          WRITE(*,*) 'quotient u ',DABS(SLNREF)/XSPLITLIM,U
	          WRITE(*,*) 'split needed, automatic split not implemented %%%'
	          STOP 3141
		END IF
	      END IF
	      SLN = SLN + DABS(SLNREF)
	      IF (NRIKT.GT.0) THEN
	        S1Q = CQ*YS1 + STQ*YS3
	        S2Q = CP*YS2 + STP*YS4
	        YS3 = -SDQ*YS1 + CQ*YS3
	        YS4 = -SDP*YS2 + CP*YS4
	      ELSE
	        S1Q = CQ*YS1 - STQ*YS3
	        S2Q = CP*YS2 - STP*YS4
	        YS3 = SDQ*YS1 + CQ*YS3
	        YS4 = SDP*YS2 + CP*YS4
	      END IF
	      YS1 = S1Q
	      YS2 = S2Q
	      IF (SLNREF.GE.0.0D0) THEN
	        YS1 = E1X*YS1
	        YS3 = E1X*YS3
	      ELSE
	        YS2 = E1X*YS2
	        YS4 = E1X*YS4
	      END IF
	    END IF
C	    now the Cinv propagation step
	    YS1 = YS1 + U*YS4
	    YS2 = YS2 + U*YS3
	    YS3 = RHO(ID) * ( YS3 - G0U*YS2 )
	    YS4 = RHO(ID) * ( YS4 - G0U*YS1 )

	  ELSE

	    IF (IYSTAB.NE.0) THEN
	      WRITE(*,*) '%%% solsifw, depths ',ZD(ID-1),ZD(ID)
	      WRITE(*,*) 'stabilization not implemented'
	      STOP 3141
	    END IF
	    IFSMTYPLOK = IFSMTYP(ID)
	    IF (NRIKT.GT.0) THEN
	      ZENDFSUB = ZD(ID)
	      Z0FSUB = ZD(ID-1)
	    ELSE
	      ZENDFSUB = ZD(ID-1)
	      Z0FSUB = ZD(ID)
	    END IF
	    YS(1) = YS1
	    YS(2) = YS2
	    YS(3) = YS3
	    YS(4) = YS4
	    CALL DPCLIN ( FSUBSOL4,4, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	    HFAKTMAX,HOFF,EPSMAGNIFMAX, ZENDFSUB, Z0FSUB,SLN,YS,
     $	    EPSFSUB,HSTEPFSUB, IND )
	    IF (IND.NE.0) STOP 6760
	    YS1 = YS(1)
	    YS2 = YS(2)
	    YS3 = YS(3)
	    YS4 = YS(4)

	  END IF

	  IF (SLN.NE.0.0D0) THEN
	    IF ( YS1.NE.(0.0D0,0.0D0) .OR. YS2.NE.(0.0D0,0.0D0) .OR. 
     $	    YS3.NE.(0.0D0,0.0D0) .OR. YS3.NE.(0.0D0,0.0D0) ) THEN
	      SLASK = DEXP(SLN)
	      YS1 = SLASK * YS1
	      YS2 = SLASK * YS2
	      YS3 = SLASK * YS3
	      YS4 = SLASK * YS4
	    END IF
	    SLN = 0.0D0
	  END IF
	  IF (INHOMOG) THEN
	    IF (NRIKT.GT.0) THEN
	      YS1 = YS1 + DIS1FW(ID)
	      YS2 = YS2 + DIS2FW(ID)
	      YS3 = YS3 + DIS3FW(ID)
	      YS4 = YS4 + DIS4FW(ID)
	    END IF
	    IF (STABKVAR) THEN
C	      stabilization now ("harksifw and zsolsifw" must have been called)
	      SLASK1 = DMAX1 ( DABS(DREAL(YT1)), DABS(DIMAG(YT1)) )
	      SLASK2 = DMAX1 ( DABS(DREAL(YT2)), DABS(DIMAG(YT2)) )
	      SLASK3 = DMAX1 ( DABS(DREAL(YT3)), DABS(DIMAG(YT3)) )
	      SLASK4 = DMAX1 ( DABS(DREAL(YT4)), DABS(DIMAG(YT4)) )
	      SLASK6 = DMAX1 ( DABS(DREAL(YT6)), DABS(DIMAG(YT6)) )
	      ICLOK = 1
	      SLASK = SLASK1
	      IF (SLASK2.GT.SLASK) THEN
	        ICLOK = 2
	        SLASK = SLASK2
	      END IF
	      IF (SLASK3.GT.SLASK) THEN
	        ICLOK = 3
	        SLASK = SLASK3
	      END IF
	      IF (SLASK4.GT.SLASK) THEN
	        ICLOK = 4
	        SLASK = SLASK4
	      END IF
	      IF (SLASK6.GT.SLASK) ICLOK = 6
	      GOTO (1,2,3,4,95,6) ICLOK
 95	      STOP 95
 1	      CSLASK = 1/YT1
	      YS1 = CSLASK * ( -Z2 + YT4*YS3 - YT2*YS4 )    ! yt5=-yt2
	      YS2 = CSLASK * ( Z1 - YT2*YS3 - YT3*YS4 )
	      GOTO 9
 2	      CSLASK = 1/YT2
	      YS1 = CSLASK * ( -Z3 - YT4*YS2 + YT6*YS4 )
	      YS3 = CSLASK * ( Z1 - YT1*YS2 - YT3*YS4 )
	      GOTO 9
 3	      CSLASK = 1/YT3
	      YS1 = CSLASK * ( -Z4 + YT2*YS2 - YT6*YS3 )    ! yt5=-yt2
	      YS4 = CSLASK * ( Z1 - YT1*YS2 - YT2*YS3 )
	      GOTO 9
 4	      CSLASK = 1/YT4
	      YS2 = CSLASK * ( -Z3 - YT2*YS1 + YT6*YS4 )
	      YS3 = CSLASK * ( Z2 + YT1*YS1 + YT2*YS4 )     ! yt5=-yt2
	      GOTO 9
 6	      CSLASK = 1/YT6
	      YS3 = CSLASK * ( -Z4 - YT3*YS1 + YT2*YS2 )    ! yt5=-yt2
	      YS4 = CSLASK * ( Z3 + YT2*YS1 + YT4*YS2 )
 9	      CONTINUE
	    END IF
	  ELSE
	    IF (STABKVAR) THEN
C	      "zero" stabilization now ("harksifw" must have been called)
	      SLASK1 = DMAX1 ( DABS(DREAL(YT1)), DABS(DIMAG(YT1)) )
	      SLASK2 = DMAX1 ( DABS(DREAL(YT2)), DABS(DIMAG(YT2)) )
	      SLASK3 = DMAX1 ( DABS(DREAL(YT3)), DABS(DIMAG(YT3)) )
	      SLASK4 = DMAX1 ( DABS(DREAL(YT4)), DABS(DIMAG(YT4)) )
	      SLASK6 = DMAX1 ( DABS(DREAL(YT6)), DABS(DIMAG(YT6)) )
	      ICLOK = 1
	      SLASK = SLASK1
	      IF (SLASK2.GT.SLASK) THEN
	        ICLOK = 2
	        SLASK = SLASK2
	      END IF
	      IF (SLASK3.GT.SLASK) THEN
	        ICLOK = 3
	        SLASK = SLASK3
	      END IF
	      IF (SLASK4.GT.SLASK) THEN
	        ICLOK = 4
	        SLASK = SLASK4
	      END IF
	      IF (SLASK6.GT.SLASK) ICLOK = 6
	      GOTO (10,20,30,40,950,60) ICLOK
 950	      STOP 950
 10	      CSLASK = 1/YT1
	      YS1 = CSLASK * ( YT4*YS3 - YT2*YS4 ) ! yt5=-yt2
	      YS2 = CSLASK * ( -YT2*YS3 - YT3*YS4 )
	      GOTO 90
 20	      CSLASK = 1/YT2
	      YS1 = CSLASK * ( -YT4*YS2 + YT6*YS4 )
	      YS3 = CSLASK * ( -YT1*YS2 - YT3*YS4 )
	      GOTO 90
 30	      CSLASK = 1/YT3
	      YS1 = CSLASK * ( YT2*YS2 - YT6*YS3 ) ! yt5=-yt2
	      YS4 = CSLASK * ( -YT1*YS2 - YT2*YS3 )
	      GOTO 90
 40	      CSLASK = 1/YT4
	      YS2 = CSLASK * ( -YT2*YS1 + YT6*YS4 )
	      YS3 = CSLASK * ( YT1*YS1 + YT2*YS4 ) ! yt5=-yt2
	      GOTO 90
 60	      CSLASK = 1/YT6
	      YS3 = CSLASK * ( -YT3*YS1 + YT2*YS2 ) ! yt5=-yt2
	      YS4 = CSLASK * ( YT2*YS1 + YT4*YS2 )
 90	      CONTINUE
	    END IF
	  END IF
	  IF (NRIKT.LT.0) THEN
	    IDQ = ID-1
	  ELSE
	    IDQ = ID
	  END IF
	  YS1FW(IDQ) = YS1
	  YS2FW(IDQ) = YS2
	  YS3FW(IDQ) = YS3
	  YS4FW(IDQ) = YS4

	END DO

	RETURN
	END






	SUBROUTINE LIQSIXFW ( NL1,NL2, UP,UP2,U,U2 )
C	complex "horizontal slowness" u , u2=u**2   ( u = u0*up , up2 = up*up )
C	complex "wave-number" rk=omeg*u , rk2=rk**2
C	propagator-matrix calculations, through fluid "half-space equal" 
C	modified from solsixfw simply by letting the shear velocity tend to 0

	USE INMODELa
	USE FWYSa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C	nol <= noldim

	INCLUDE '../INPARAM.INC'

	INTEGER NL1,NL2
	COMPLEX*16 UP,UP2,U,U2

	INTEGER ID
	REAL*8 RH2, THK
	COMPLEX*16 YP1,YP2, 
     $	URA, URANEGI, RKRA, P0, 
     $	ALFAPDOWN, PD2,PD4, 
     $	PDFAKT, CSLASK1, URVALUE

	IF (NL1.GT.NL2) STOP 65714

	YP1 = YS2FW(NL1-1)
	YP2 = YS4FW(NL1-1)

	RH2 = 1/RHO(NOL)                           ! could be stored

C	branch-cut choice now important
	URA = URVALUE (ICUTTYP,AIMCUT, U0,UPBPA,UPBPA2,UPREF, UP,UP2)
	IF ( ISHEETTYP.EQ.2 .OR. ISHEETTYP.EQ.4 ) URA = -URA
	RKRA = OMEG*URA
	URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura

	ALFAPDOWN = 0.5D0 * ( YP1/URANEGI + RH2*YP2 )

	PD2 =  ALFAPDOWN * URANEGI
	PD4 =  ALFAPDOWN * RHO(NOL)

	YP1 = PD2
	YP2 = PD4

C	change of "ys2,ys4" could be used as a 
C	check of the "mismatch" at depth zd(nl1-1)

	DO ID = NL1, NL2

	  THK = ZD(ID) - ZD(NL1-1)
	  P0 = THK*RKRA

	  PDFAKT = DCMPLX ( DCOS(DREAL(P0)), DSIN(DREAL(P0)) )
	  CSLASK1 = DEXP(-DIMAG(P0)) * PDFAKT
	  YP1 = CSLASK1*PD2
	  YP2 = CSLASK1*PD4

	  YS2FW(ID) = YP1
	  YS4FW(ID) = YP2

	END DO

	RETURN
	END






	SUBROUTINE SOLSIXFW ( NL1,NL2, UP,UP2,U,U2 )
C	complex "horizontal slowness" u , u2=u**2   ( u = u0*up , up2 = up*up )
C	complex "wave-number" rk=omeg*u , rk2=rk**2
C	propagator-matrix calculations, through solid "half-space equal" 
C	layers (formulas from Aki (7.54-56) but slightly modified), 
C	for the "Aki (!) displacement-stress vector (Aki & Richards 1980, 
C	7.4.2) but with omeg* for the first two components"
C	the start values for the column-vector are given as input by 
C	(ys1,ys2,ys3,ys4)
C	note that ys1,ys2,ys3,ys4 are then updated successively, here 
C	each time starting from the "reference level" and using "only two 
C	vectors"

	USE INMODELa
	USE FWYSa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C	nol <= noldim

	INCLUDE '../INPARAM.INC'

	INTEGER NL1,NL2
	COMPLEX*16 UP,UP2,U,U2

	INTEGER ID
	REAL*8 RH2, THK
	COMPLEX*16 YS1,YS2,YS3,YS4, 
     $	G0, G0U, G1,G2, GG,G2R, URA,URB, 
     $	URANEGI,URBNEGI, RKRA,RKRB, P0,Q0, 
     $	ALFAPDOWN,BETASDOWN, PD1,PD2,PD3,PD4, SD1,SD2,SD3,SD4, 
     $	PDFAKT,SDFAKT, CSLASK1,CSLASK2, URVALUE

	IF (NL1.GT.NL2) STOP 65714

	YS1 = YS1FW(NL1-1)
	YS2 = YS2FW(NL1-1)
	YS3 = YS3FW(NL1-1)
	YS4 = YS4FW(NL1-1)

	RH2 = 1/RHO(NOL)                           ! could be stored
	G0 = TWOVB2(NOL)                           ! local twovb2(nol)
	G0U = G0*U
	G1 = G0*U2
	G2 = 1-G1
	GG = RHO(NOL)*G0U
	G2R = RHO(NOL)*G2

C	branch-cut choice now important
	URA = URVALUE (ICUTTYP,AIMCUT, U0,UPBPA,UPBPA2,UPREF, UP,UP2)
	IF ( ISHEETTYP.EQ.2 .OR. ISHEETTYP.EQ.4 ) URA = -URA
	RKRA = OMEG*URA

C	branch-cut choice now important
	URB = URVALUE (ICUTTYP,AIMCUT, U0,UPBPB,UPBPB2,UPREF, UP,UP2)
	IF ( ISHEETTYP.EQ.3 .OR. ISHEETTYP.EQ.4 ) URB = -URB
	RKRB = OMEG*URB

	URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb

	ALFAPDOWN = 0.5D0 * ( G0U*YS1 + (G2/URANEGI)*YS2 + 
     $	RH2*(YS4-(U/URANEGI)*YS3) )
	BETASDOWN = 0.5D0 * ( (G2/URBNEGI)*YS1 + G0U*YS2 + 
     $	RH2*(YS3-(U/URBNEGI)*YS4) )

	PD1 =  ALFAPDOWN * U
	PD2 =  ALFAPDOWN * URANEGI
	PD3 = -ALFAPDOWN * (GG*URANEGI)
	PD4 =  ALFAPDOWN * G2R
	SD1 =  BETASDOWN * URBNEGI
	SD2 =  BETASDOWN * U
	SD3 =  BETASDOWN * G2R
	SD4 = -BETASDOWN * (GG*URBNEGI)

	YS1 = PD1 + SD1
	YS2 = PD2 + SD2
	YS3 = PD3 + SD3
	YS4 = PD4 + SD4

C	change of "ys1,ys2,ys3,ys4" could be used as a 
C	check of the "mismatch" at depth zd(nl1-1)

	DO ID = NL1, NL2

	  THK = ZD(ID) - ZD(NL1-1)
	  P0 = THK*RKRA
	  Q0 = THK*RKRB

	  PDFAKT = DCMPLX ( DCOS(DREAL(P0)), DSIN(DREAL(P0)) )
	  SDFAKT = DCMPLX ( DCOS(DREAL(Q0)), DSIN(DREAL(Q0)) )
	  CSLASK1 = DEXP(-DIMAG(P0)) * PDFAKT
	  CSLASK2 = DEXP(-DIMAG(Q0)) * SDFAKT
	  YS1 = CSLASK1*PD1 + CSLASK2*SD1
	  YS2 = CSLASK1*PD2 + CSLASK2*SD2
	  YS3 = CSLASK1*PD3 + CSLASK2*SD3
	  YS4 = CSLASK1*PD4 + CSLASK2*SD4

	  YS1FW(ID) = YS1
	  YS2FW(ID) = YS2
	  YS3FW(ID) = YS3
	  YS4FW(ID) = YS4

	END DO

	RETURN
	END
