C	subroutines for rpress.f, branch-cut integral(s)

C	Sven Ivansson 1990







	SUBROUTINE DHBRANCHSUM (UP, RLNDET,CDET, DHFAKT)
C	for coincident parts of P-cut & S-cut, still assuming that 
C	the P-cut is below to the left of the S-cut (for "up")

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       ... , nhyd <= nhyddim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 UP

	REAL*8 RLNDET
	COMPLEX*16 CDET, DHFAKT(*)

	INTEGER IH
	COMPLEX*16 DHFAKTQ(nhyddim)

	IPS = 1
	CALL DHBRANCH (UP, RLNDET,CDET, DHFAKTQ)

	IPS = 2
	CALL DHBRANCH (UP, RLNDET,CDET, DHFAKT)

	DO IH = 1, NHYD
	  DHFAKT(IH) = DHFAKT(IH) + DHFAKTQ(IH)
	END DO

	RLNDET = 0.0D0
	CDET = (1.0D0,0.0D0)

	RETURN
	END







	SUBROUTINE DHBRANCH (UP, RLNDET,CDET, DHFAKT)
C	modelled after dhcomp
C	underlying assumptions: 
C	 * up in the closed first quadrant   (  Re(up) >= 0 , Im(up) >= 0  )
C	 * branch-point(s) in the open right half-plane (for "up")
C	 * cut(s) in the closed first quadrant (for "up")
C	 * P-cut below to the left of S-cut (if the latter present) (for "up")

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       ... , nhyd <= nhyddim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INPARAM.INC'
	INCLUDE 'FSUBCOM.INC'
	INCLUDE 'HARKLIQ.INC'
	INCLUDE 'BRANCHPS.INC'
	INCLUDE 'INSPEC.INC'

	COMPLEX*16 UP

	REAL*8 RLNDET
	COMPLEX*16 CDET, DHFAKT(*)

	INTEGER ISIDE, IHARK0, IH
	REAL*8 TLNL,TLNU, RHOCDET
	COMPLEX*16 U,RK,UP2,U2,RK2, YPL1,YPL2, YPU1,YPU2, RW, 
     $	URA,URB, URVALUE, 
     $	DHFAKTQ(nhyddim)

	U = U0*UP
	RK = RK0*UP
	UP2 = UP*UP
	U2 = U*U
	RK2 = RK*RK

	UCOP = U
	U2COP = U2
	RK2GOMEG = RK2/OMEG

	IF (IRKSPEC.NE.0) DLNSPEC = ORMINSPEC * DIMAG(RK)

	DO ISIDE = 1, 2
C	  iside = 1   for the contribution from the lower side of the cut (-)
C	  iside = 2   for the contribution from the upper side of the cut (+)

	  URA = URVALUE(2,0.0D0, U0,UPBPA,UPBPA2,(0.0D0,0.0D0), 
     $	  UP,UP2)  ! with "icuttyp 2", "any aimcut", "upref below to the left"
	  IF (IPS.EQ.1) THEN
	    IF (ISIDE.EQ.2) URA = -URA
	  END IF      ! assumed that "P-cut below to the left of S-cut"
	  IF (ISOLID(NOL).NE.0) THEN
C	    solid half-space below
	    URB = URVALUE(2,0.0D0, U0,UPBPB,UPBPB2,(0.0D0,0.0D0), 
     $	    UP,UP2) ! with "icuttyp 2", "any aimcut", "upref below to the left"
	    IF (IPS.EQ.2) THEN
	      IF (ISIDE.EQ.2) URB = -URB
	    END IF      ! assumed that "P-cut below to the left of S-cut"
	    IHARK0 = 1
	    YTL1 = (1.0D0,0.0D0)
	    YTL2 = (0.0D0,0.0D0)
	    YTL3 = DCMPLX(-DIMAG(URA),DREAL(URA))       !  i*ura
	    YTL4 = DCMPLX(DIMAG(URB),-DREAL(URB))       ! -i*urb
	    YTL6 = -URA*URB
	    TLNL = 0.0D0
	    CALL HARKSI ( NOL-1,NSW+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    IHARK0, IHARKLIQ, RHOCDET, YTL1,YTL2,YTL3,YTL4,YTL6, TLNL, 
     $	    YPL1,YPL2 )
	  ELSE
C	    fluid half-space below
	    TLNL = 0.0D0
	    YPL1 = DCMPLX(DIMAG(URA),-DREAL(URA))       !  -i*ura
	    YPL2 = RHO(NOL)
	  END IF

	  IF (RHO(0).GT.0.0D0) STOP 45166   ! upper half-space not implemented
	  IF (NSS.GT.0) THEN
C	    free (or rigid) solid above
	    IHARK0 = -2
	    YTU1 = (0.0D0,0.0D0)
	    YTU2 = (0.0D0,0.0D0)
	    YTU3 = (0.0D0,0.0D0)
	    YTU4 = (0.0D0,0.0D0)
	    YTU6 = (1.0D0,0.0D0)
	    IF (RHO(0).LT.0.0D0) THEN
	      YTU1 = (1.0D0,0.0D0)
	      YTU6 = (0.0D0,0.0D0)
	    END IF
	    TLNU = 0.0D0
	    CALL HARKSI ( 1,NSS,+1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $	    IHARKLIQ, RHOCDET, YTU1,YTU2,YTU3,YTU4,YTU6, TLNU, 
     $	    YPU1,YPU2 )
C	    should "zd(0)=zd(nss)" one will get "ypu1:=-ytu3=0 ypu2:=ytu1=0" 
C	    and a pathological "mode" with "nonvanishing r3 in the infinitely 
C	    thin solid layer(s) on top" (but else "zero") is found
	  ELSE
C	    free (or rigid) fluid above
	    TLNU = 0.0D0
	    YPU1 = (1.0D0,0.0D0)
	    YPU2 = (0.0D0,0.0D0)
	    IF (RHO(0).LT.0.0D0) THEN
	      YPU1 = (0.0D0,0.0D0)
	      YPU2 = (1.0D0,0.0D0)
	    END IF
	  END IF

	  RW = RW0*RK

	  IF (IDET.NE.0) RLNDET = TLNL + TLNU
	  IF (ISIDE.EQ.1) THEN
	    CALL LIQSI ( OMEG,OMEG2,U,U2,RK,RK2, RW0,RW, YPL1,YPL2, 
     $	    YPU1,YPU2, IDET,RHOCDET,RLNDET,CDET, CLIQL,CLIQU, DHFAKTQ )
	  ELSE
	    CALL LIQSI ( OMEG,OMEG2,U,U2,RK,RK2, RW0,RW, YPL1,YPL2, 
     $	    YPU1,YPU2, IDET,RHOCDET,RLNDET,CDET, CLIQL,CLIQU, DHFAKT )
	  END IF

	END DO

	DO IH = 1, NHYD
	  DHFAKT(IH) = DHFAKT(IH) - DHFAKTQ(IH)
	END DO

	RLNDET = 0.0D0
	CDET = (1.0D0,0.0D0)

	RETURN
	END








	SUBROUTINE DHBRANCHSS (SSZ, RLNDET,CDET, DHFAKT)
C	modelled after dhbranch, parameter ss on vertical cut

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       ... , nhyd <= nhyddim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INPARAM.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 SSZ

	REAL*8 RLNDET
	COMPLEX*16 CDET, DHFAKT(*)

	INTEGER IH
	REAL*8 SSLASK
	COMPLEX*16 UP, DUPDSS

	IF ( DIMAG(SSZ).NE.0.0D0 .OR. DREAL(SSZ).LT.0.0D0 ) STOP 560 ! on cut
	SSLASK = SSZ

	IF (SSLASK.NE.0.0D0) THEN
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA + DCMPLX(0.0D0,SSLASK*SSLASK)
	  ELSE
	    UP = UPBPB + DCMPLX(0.0D0,SSLASK*SSLASK)
	  END IF
	  DUPDSS = 2 * DCMPLX(0.0D0,SSLASK)
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	  END IF
	  DUPDSS = (0.0D0,0.0D0)
	END IF

	CALL DHBRANCH (UP, RLNDET,CDET, DHFAKT)

	DO IH = 1, NHYD
	  DHFAKT(IH) = DUPDSS * DHFAKT(IH)
	END DO

	RLNDET = 0.0D0
	CDET = (1.0D0,0.0D0)

	RETURN
	END






	SUBROUTINE RBRANCHSS (SSZ, RFAKT)
C	modelled after rcomp, parameter ss on vertical cut

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 SSZ

	COMPLEX*16 RFAKT(*)

	REAL*8 SSLASK
	COMPLEX*16 UP

	IF ( DIMAG(SSZ).NE.0.0D0 .OR. DREAL(SSZ).LT.0.0D0 ) STOP 560 ! on cut
	SSLASK = SSZ

	IF (SSLASK.NE.0.0D0) THEN
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA + DCMPLX(0.0D0,SSLASK*SSLASK)
	  ELSE
	    UP = UPBPB + DCMPLX(0.0D0,SSLASK*SSLASK)
	  END IF
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	  END IF
	END IF

	CALL RCOMP (UP, RFAKT)

	RETURN
	END







	SUBROUTINE RBRANCHPOSSS (SSZ, RFAKT)
C	modelled after rcomppos, parameter ss on vertical cut

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 SSZ

	COMPLEX*16 RFAKT(*)

	REAL*8 SSLASK
	COMPLEX*16 UP


	IF ( DIMAG(SSZ).NE.0.0D0 .OR. DREAL(SSZ).LT.0.0D0 ) STOP 560 ! on cut
	SSLASK = SSZ

	IF (SSLASK.NE.0.0D0) THEN
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA + DCMPLX(0.0D0,SSLASK*SSLASK)
	  ELSE
	    UP = UPBPB + DCMPLX(0.0D0,SSLASK*SSLASK)
	  END IF
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	  END IF
	END IF

	CALL RCOMPPOS (UP, RFAKT)

	RETURN
	END







	SUBROUTINE DHBRANCHT (TZ, RLNDET,CDET, DHFAKT)
C	modelled after dhcomp, parameter t on hyperbolic cut

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       ... , nhyd <= nhyddim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INPARAM.INC'
	INCLUDE 'FSUBCOM.INC'
	INCLUDE 'HARKLIQ.INC'
	INCLUDE 'BRANCHPS.INC'
	INCLUDE 'INSPEC.INC'

	COMPLEX*16 TZ

	REAL*8 RLNDET
	COMPLEX*16 CDET, DHFAKT(*)

	INTEGER ISIDE, IHARK0, IH
	REAL*8 TLNL,TLNU, TSLASK, SLASKROT, YSLASK, RHOCDET
	COMPLEX*16 UP, U,RK,UP2,U2,RK2, YPL1,YPL2, YPU1,YPU2, RW, 
     $	URA,URB, DUPDT, 
     $	DHFAKTQ(nhyddim)

	IF ( DIMAG(TZ).NE.0.0D0 .OR. DREAL(TZ).LT.TPS(IPS) ) STOP 560 ! on cut
	TSLASK = TZ
	SLASKROT = DSQRT ( TSLASK*TSLASK + 4*PSPROD(IPS) )

	IF (TSLASK.NE.TPS(IPS)) THEN
	  YSLASK = 0.5D0 * ( TSLASK + SLASKROT )
	  UP = DCMPLX ( PSPROD(IPS)/YSLASK , YSLASK )
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	  END IF   ! should not be used, use the "s-parameter" close to upbp
	END IF
	DUPDT = (-1/SLASKROT) * DCONJG (UP)

	U = U0*UP
	RK = RK0*UP
	UP2 = UP*UP
	U2 = U*U
	RK2 = RK*RK

	UCOP = U
	U2COP = U2
	RK2GOMEG = RK2/OMEG

	IF (IRKSPEC.NE.0) DLNSPEC = ORMINSPEC * DIMAG(RK)

	DO ISIDE = 1, 2
C	  iside = 1   for the contribution from the lower side of the cut (-)
C	  iside = 2   for the contribution from the upper side of the cut (+)

	  IF (IPS.EQ.1) THEN
	    IF (TSLASK.LE.TPS(1)) THEN
	      URA = (0.0D0,0.0D0)
	    ELSE
	      URA = DSQRT ( HPS(1) + TSLASK*SLASKROT )    ! cdsqrt(upbpa2-up2)
	    END IF
	    IF (ISIDE.EQ.2) URA = -URA
	  ELSE
	    URA = CDSQRT ( UPBPA2 - UP2 )
	    IF (DIMAG(URA).LT.0.0D0) URA = -URA          ! cdsqrt12
	    URA = U0 * URA
	  END IF
	  IF (ISOLID(NOL).NE.0) THEN
C	    solid half-space below
	    IF (IPS.EQ.2) THEN
	      IF (TSLASK.LE.TPS(2)) THEN
	        URB = (0.0D0,0.0D0)
	      ELSE
	        URB = DSQRT ( HPS(2) + TSLASK*SLASKROT )  ! cdsqrt(upbpb2-up2)
	      END IF
	      IF (ISIDE.EQ.2) URB = -URB
	    ELSE
	      URB = CDSQRT ( UPBPB2 - UP2 )
	      IF (DIMAG(URB).LT.0.0D0) URB = -URB          ! cdsqrt12
	      URB = U0 * URB
	    END IF
	    IHARK0 = 1
	    YTL1 = (1.0D0,0.0D0)
	    YTL2 = (0.0D0,0.0D0)
	    YTL3 = DCMPLX(-DIMAG(URA),DREAL(URA))       !  i*ura
	    YTL4 = DCMPLX(DIMAG(URB),-DREAL(URB))       ! -i*urb
	    YTL6 = -URA*URB
	    TLNL = 0.0D0
	    CALL HARKSI ( NOL-1,NSW+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    IHARK0, IHARKLIQ, RHOCDET, YTL1,YTL2,YTL3,YTL4,YTL6, TLNL, 
     $	    YPL1,YPL2 )
	  ELSE
C	    fluid half-space below
	    TLNL = 0.0D0
	    YPL1 = DCMPLX(DIMAG(URA),-DREAL(URA))       !  -i*ura
	    YPL2 = RHO(NOL)
	  END IF

	  IF (RHO(0).GT.0.0D0) STOP 45166   ! upper half-space not implemented
	  IF (NSS.GT.0) THEN
C	    free (or rigid) solid above
	    IHARK0 = -2
	    YTU1 = (0.0D0,0.0D0)
	    YTU2 = (0.0D0,0.0D0)
	    YTU3 = (0.0D0,0.0D0)
	    YTU4 = (0.0D0,0.0D0)
	    YTU6 = (1.0D0,0.0D0)
	    IF (RHO(0).LT.0.0D0) THEN
	      YTU1 = (1.0D0,0.0D0)
	      YTU6 = (0.0D0,0.0D0)
	    END IF
	    TLNU = 0.0D0
	    CALL HARKSI ( 1,NSS,+1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $	    IHARKLIQ, RHOCDET, YTU1,YTU2,YTU3,YTU4,YTU6, TLNU, 
     $	    YPU1,YPU2 )
C	    should "zd(0)=zd(nss)" one will get "ypu1:=-ytu3=0 ypu2:=ytu1=0" 
C	    and a pathological "mode" with "nonvanishing r3 in the infinitely 
C	    thin solid layer(s) on top" (but else "zero") is found
	  ELSE
C	    free (or rigid) fluid above
	    TLNU = 0.0D0
	    YPU1 = (1.0D0,0.0D0)
	    YPU2 = (0.0D0,0.0D0)
	    IF (RHO(0).LT.0.0D0) THEN
	      YPU1 = (0.0D0,0.0D0)
	      YPU2 = (1.0D0,0.0D0)
	    END IF
	  END IF

	  RW = RW0*RK

	  IF (IDET.NE.0) RLNDET = TLNL + TLNU
	  IF (ISIDE.EQ.1) THEN
	    CALL LIQSI ( OMEG,OMEG2,U,U2,RK,RK2, RW0,RW, YPL1,YPL2, 
     $	    YPU1,YPU2, IDET,RHOCDET,RLNDET,CDET, CLIQL,CLIQU, DHFAKTQ )
	  ELSE
	    CALL LIQSI ( OMEG,OMEG2,U,U2,RK,RK2, RW0,RW, YPL1,YPL2, 
     $	    YPU1,YPU2, IDET,RHOCDET,RLNDET,CDET, CLIQL,CLIQU, DHFAKT )
	  END IF

	END DO

	DO IH = 1, NHYD
	  DHFAKT(IH) = DUPDT * ( DHFAKT(IH) - DHFAKTQ(IH) )
	END DO

	RLNDET = 0.0D0
	CDET = (1.0D0,0.0D0)

	RETURN
	END





	SUBROUTINE RBRANCHT (TZ, RFAKT)
C	modelled after rcomp, parameter t on hyperbolic cut

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 TZ

	COMPLEX*16 RFAKT(*)

	REAL*8 TSLASK, SLASKROT, YSLASK
	COMPLEX*16 UP

	IF ( DIMAG(TZ).NE.0.0D0 .OR. DREAL(TZ).LT.TPS(IPS) ) STOP 560 ! on cut
	TSLASK = TZ

	IF (TSLASK.NE.TPS(IPS)) THEN
	  SLASKROT = DSQRT ( TSLASK*TSLASK + 4*PSPROD(IPS) )
	  YSLASK = 0.5D0 * ( TSLASK + SLASKROT )
	  UP = DCMPLX ( PSPROD(IPS)/YSLASK , YSLASK )
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	    UP = UPBPB
	  END IF   ! should not be used, use the "s-parameter" close to upbp
	END IF

	CALL RCOMP (UP, RFAKT)

	RETURN
	END






	SUBROUTINE RBRANCHPOST (TZ, RFAKT)
C	modelled after rcomppos, parameter t on hyperbolic cut

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 TZ

	COMPLEX*16 RFAKT(*)

	REAL*8 TSLASK, SLASKROT, YSLASK
	COMPLEX*16 UP

	IF ( DIMAG(TZ).NE.0.0D0 .OR. DREAL(TZ).LT.TPS(IPS) ) STOP 560 ! on cut
	TSLASK = TZ

	IF (TSLASK.NE.TPS(IPS)) THEN
	  SLASKROT = DSQRT ( TSLASK*TSLASK + 4*PSPROD(IPS) )
	  YSLASK = 0.5D0 * ( TSLASK + SLASKROT )
	  UP = DCMPLX ( PSPROD(IPS)/YSLASK , YSLASK )
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	    UP = UPBPB
	  END IF   ! should not be used, use the "s-parameter" close to upbp
	END IF

	CALL RCOMPPOS (UP, RFAKT)

	RETURN
	END










	SUBROUTINE DHBRANCHS (SZ, RLNDET,CDET, DHFAKT)
C	modelled after dhcomp, parameter s on hyperbolic cut

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       ... , nhyd <= nhyddim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INPARAM.INC'
	INCLUDE 'FSUBCOM.INC'
	INCLUDE 'HARKLIQ.INC'
	INCLUDE 'BRANCHPS.INC'
	INCLUDE 'INSPEC.INC'

	COMPLEX*16 SZ

	REAL*8 RLNDET
	COMPLEX*16 CDET, DHFAKT(*)

	INTEGER ISIDE, IHARK0, IH
	REAL*8 TLNL,TLNU, SSLASK, SQ, SLASKROT, XSLASK, RHOCDET
	COMPLEX*16 UP, U,RK,UP2,U2,RK2, YPL1,YPL2, YPU1,YPU2, RW, 
     $	URA,URB, DUPDS, 
     $	DHFAKTQ(nhyddim)

	IF ( DIMAG(SZ).NE.0.0D0 .OR. DREAL(SZ).LT.0.0D0 ) STOP 560 ! on cut
	SSLASK = SZ

	IF (SSLASK.NE.0.0D0) THEN
	  SQ = HPS(IPS) - SSLASK*SSLASK
	  SLASKROT = DSQRT ( SQ*SQ + 4*PSPROD(IPS)*PSPROD(IPS) )
	  XSLASK = DSQRT ( 0.5D0 * ( SQ + SLASKROT ) )
	  UP = DCMPLX ( XSLASK , PSPROD(IPS)/XSLASK )
	  DUPDS = -(SSLASK/SLASKROT) * DCONJG (UP)
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	  END IF
	  DUPDS = (0.0D0,0.0D0)
	END IF

	U = U0*UP
	RK = RK0*UP
	UP2 = UP*UP
	U2 = U*U
	RK2 = RK*RK

	UCOP = U
	U2COP = U2
	RK2GOMEG = RK2/OMEG

	IF (IRKSPEC.NE.0) DLNSPEC = ORMINSPEC * DIMAG(RK)

	DO ISIDE = 1, 2
C	  iside = 1   for the contribution from the lower side of the cut (-)
C	  iside = 2   for the contribution from the upper side of the cut (+)

	  IF (IPS.EQ.1) THEN
	    URA = DCMPLX (SSLASK,0.0D0)
	    IF (ISIDE.EQ.2) URA = -URA
	  ELSE
	    URA = CDSQRT ( UPBPA2 - UP2 )
	    IF (DIMAG(URA).LT.0.0D0) URA = -URA          ! cdsqrt12
	    URA = U0 * URA
	  END IF
	  IF (ISOLID(NOL).NE.0) THEN
C	    solid half-space below
	    IF (IPS.EQ.2) THEN
	      URB = DCMPLX (SSLASK,0.0D0)
	      IF (ISIDE.EQ.2) URB = -URB
	    ELSE
	      URB = CDSQRT ( UPBPB2 - UP2 )
	      IF (DIMAG(URB).LT.0.0D0) URB = -URB          ! cdsqrt12
	      URB = U0 * URB
	    END IF
	    IHARK0 = 1
	    YTL1 = (1.D0,0.0D0)
	    YTL2 = (0.0D0,0.0D0)
	    YTL3 = DCMPLX(-DIMAG(URA),DREAL(URA))       !  i*ura
	    YTL4 = DCMPLX(DIMAG(URB),-DREAL(URB))       ! -i*urb
	    YTL6 = -URA*URB
	    TLNL = 0.0D0
	    CALL HARKSI ( NOL-1,NSW+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    IHARK0, IHARKLIQ, RHOCDET, YTL1,YTL2,YTL3,YTL4,YTL6, TLNL, 
     $	    YPL1,YPL2 )
	  ELSE
C	    fluid half-space below
	    TLNL = 0.0D0
	    YPL1 = DCMPLX(DIMAG(URA),-DREAL(URA))       !  -i*ura
	    YPL2 = RHO(NOL)
	  END IF

	  IF (RHO(0).GT.0.0D0) STOP 45166   ! upper half-space not implemented
	  IF (NSS.GT.0) THEN
C	    free (or rigid) solid above
	    IHARK0 = -2
	    YTU1 = (0.0D0,0.0D0)
	    YTU2 = (0.0D0,0.0D0)
	    YTU3 = (0.0D0,0.0D0)
	    YTU4 = (0.0D0,0.0D0)
	    YTU6 = (1.0D0,0.0D0)
	    IF (RHO(0).LT.0.0D0) THEN
	      YTU1 = (1.0D0,0.0D0)
	      YTU6 = (0.0D0,0.0D0)
	    END IF
	    TLNU = 0.0D0
	    CALL HARKSI ( 1,NSS,+1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $	    IHARKLIQ, RHOCDET, YTU1,YTU2,YTU3,YTU4,YTU6, TLNU, 
     $	    YPU1,YPU2 )
C	    should "zd(0)=zd(nss)" one will get "ypu1:=-ytu3=0 ypu2:=ytu1=0" 
C	    and a pathological "mode" with "nonvanishing r3 in the infinitely 
C	    thin solid layer(s) on top" (but else "zero") is found
	  ELSE
C	    free (or rigid) fluid above
	    TLNU = 0.0D0
	    YPU1 = (1.0D0,0.0D0)
	    YPU2 = (0.0D0,0.0D0)
	    IF (RHO(0).LT.0.0D0) THEN
	      YPU1 = (0.0D0,0.0D0)
	      YPU2 = (1.0D0,0.0D0)
	    END IF
	  END IF

	  RW = RW0*RK

	  IF (IDET.NE.0) RLNDET = TLNL + TLNU
	  IF (ISIDE.EQ.1) THEN
	    CALL LIQSI ( OMEG,OMEG2,U,U2,RK,RK2, RW0,RW, YPL1,YPL2, 
     $	    YPU1,YPU2, IDET,RHOCDET,RLNDET,CDET, CLIQL,CLIQU, DHFAKTQ )
	  ELSE
	    CALL LIQSI ( OMEG,OMEG2,U,U2,RK,RK2, RW0,RW, YPL1,YPL2, 
     $	    YPU1,YPU2, IDET,RHOCDET,RLNDET,CDET, CLIQL,CLIQU, DHFAKT )
	  END IF

	END DO

	DO IH = 1, NHYD
	  DHFAKT(IH) = DUPDS * ( DHFAKT(IH) - DHFAKTQ(IH) )
	END DO

	RLNDET = 0.0D0
	CDET = (1.0D0,0.0D0)

	RETURN
	END






	SUBROUTINE RBRANCHS (SZ, RFAKT)
C	modelled after rcomp, parameter s on hyperbolic cut

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 SZ

	COMPLEX*16 RFAKT(*)

	REAL*8 SSLASK, SQ, SLASKROT, XSLASK
	COMPLEX*16 UP

	IF ( DIMAG(SZ).NE.0.0D0 .OR. DREAL(SZ).LT.0.0D0 ) STOP 560 ! on cut
	SSLASK = SZ

	IF (SSLASK.NE.0.0D0) THEN
	  SQ = HPS(IPS) - SSLASK*SSLASK
	  SLASKROT = DSQRT ( SQ*SQ + 4*PSPROD(IPS)*PSPROD(IPS) )
	  XSLASK = DSQRT ( 0.5D0 * ( SQ + SLASKROT ) )
	  UP = DCMPLX ( XSLASK , PSPROD(IPS)/XSLASK )
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	  END IF
	END IF

	CALL RCOMP (UP, RFAKT)

	RETURN
	END






	SUBROUTINE RBRANCHPOSS (SZ, RFAKT)
C	modelled after rcomppos, parameter s on hyperbolic cut

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'
	INCLUDE 'BRANCHPS.INC'

	COMPLEX*16 SZ

	COMPLEX*16 RFAKT(*)

	REAL*8 SSLASK, SQ, SLASKROT, XSLASK
	COMPLEX*16 UP

	IF ( DIMAG(SZ).NE.0.0D0 .OR. DREAL(SZ).LT.0.0D0 ) STOP 560 ! on cut
	SSLASK = SZ

	IF (SSLASK.NE.0.0D0) THEN
	  SQ = HPS(IPS) - SSLASK*SSLASK
	  SLASKROT = DSQRT ( SQ*SQ + 4*PSPROD(IPS)*PSPROD(IPS) )
	  XSLASK = DSQRT ( 0.5D0 * ( SQ + SLASKROT ) )
	  UP = DCMPLX ( XSLASK , PSPROD(IPS)/XSLASK )
	ELSE
	  IF (IPS.EQ.1) THEN
	    UP = UPBPA
	  ELSE
	    UP = UPBPB
	  END IF
	END IF

	CALL RCOMPPOS (UP, RFAKT)

	RETURN
	END