C	subroutines and functions for rpress.f etc
C	urvalue  muir  predhcompas  fsubliq2  ... 

C	Sven Ivansson  1990






	COMPLEX*16 FUNCTION URVALUE ( ICUTTYP,AIMCUT, U0,UPBP,UPBP2,
     $	UPREF, UP,UP2 )
C	computation of vertical slowness urvalue
C	upbp is the branch-point for up=u/u0 ("upbpa" or "upbpb" in the 
C	usual calls) with  0 <= arg(upbp) < pi   (the factor "exp(-i*omeg*t)" 
C	is behind)
C	note: the "+" sheet is assumed, meaning that for real and positive 
C	    omeg and u(=u0*up), the wave  exp(i*omeg*(urvalue*z-t))  is 
C	    downgoing or exponentially decreasing when the depth z increases 
C	    (except at the branch-points +-upbp, where "urvalue" is zero and 
C	    up/down waves are meaningless)
C	    correction for "isheettyp" must be done afterwards if the "-" 
C	    sheet is desired (OBS)

	IMPLICIT NONE

	INTEGER ICUTTYP
	REAL*8 AIMCUT
	COMPLEX*16 U0,UPBP,UPBP2,UPREF, UP,UP2

	COMPLEX*16 Z2, MUIR

	IF ( DREAL(UPBP).LE.0.0D0 .OR. DIMAG(UPBP).LT.0.0D0 ) THEN
	  WRITE(*,*) '%%% upbp=',UPBP
	  STOP 5617    ! not implemented yet
	END IF

	IF (ICUTTYP.EQ.1) THEN                  ! "axis-hyperbolic" cut for up
C	  since the "+" sheet is assumed, with this cut and for real and 
C	  positive omeg, the wave  exp(i*omeg*(urvalue*z-t))  will be 
C	  downgoing or exponentially decreasing when the depth z increases 
C	  (except at the branch-point +-upbp, where "urvalue" is zero and 
C	  up/down waves are meaningless)
C	  with the "+" sheet and for real and positive omeg this cut gives 
C	  stable functions in the upper half-plane (the factor 
C	  "exp(-i*omeg*t)" is behind)
	  IF (UP2.EQ.UPBP2) THEN
	    URVALUE = (0.0D0,0.0D0)
	  ELSE
	    URVALUE = CDSQRT ( UPBP2 - UP2 )
	    IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE       ! cdsqrt12
	    URVALUE = U0 * URVALUE
	  END IF
	ELSE IF (ICUTTYP.EQ.2) THEN             ! "vertical" cut for up
	  IF (UP2.EQ.UPBP2) THEN
	    URVALUE = (0.0D0,0.0D0)
	  ELSE IF (DREAL(UP).EQ.DREAL(UPBP)) THEN
	    IF (UP.EQ.UPBP) THEN
	      URVALUE = (0.0D0,0.0D0)
	    ELSE IF (DIMAG(UP).GT.DIMAG(UPBP)) THEN
	      IF ( DREAL(UPREF) .GT. DREAL(UPBP) ) THEN
	        URVALUE = CDSQRT ( UPBP2 - UP2 )
	        IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      ELSE
	        URVALUE = U0 * CDSQRT ( UPBP2 - UP2 )
	      END IF
	    ELSE
	      URVALUE = CDSQRT ( UPBP2 - UP2 )
	      IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    END IF
	  ELSE IF (DREAL(UP).EQ.DREAL(-UPBP)) THEN
	    IF (UP.EQ.-UPBP) THEN
	      URVALUE = (0.0D0,0.0D0)
	    ELSE IF (DIMAG(UP).LT.DIMAG(-UPBP)) THEN
	      IF ( DREAL(UPREF) .LT. DREAL(-UPBP) ) THEN
	        URVALUE = CDSQRT ( UPBP2 - UP2 )
	        IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      ELSE
	        URVALUE = U0 * CDSQRT ( UPBP2 - UP2 )
	      END IF
	    ELSE
	      URVALUE = CDSQRT ( UPBP2 - UP2 )
	      IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    END IF
	  ELSE
	    IF ( DABS(DREAL(UP)) .GT. DABS(DREAL(UPBP)) ) THEN
	      URVALUE = CDSQRT ( UPBP2 - UP2 )
	      IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    ELSE
	      URVALUE = U0 * CDSQRT ( UPBP2 - UP2 )
	    END IF
	  END IF
	ELSE IF (ICUTTYP.EQ.3) THEN             ! "mixed" cut for up
	  IF ( DABS(DIMAG(UP)) .GT. AIMCUT ) THEN
C	    as for "icuttyp.eq.1"
	    IF (UP2.EQ.UPBP2) THEN
	      URVALUE = (0.0D0,0.0D0)
	    ELSE
	      URVALUE = CDSQRT ( UPBP2 - UP2 )
	      IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    END IF
	  ELSE
C	    as for "icuttyp.eq.2"
	    IF (UP2.EQ.UPBP2) THEN
	      URVALUE = (0.0D0,0.0D0)
	    ELSE IF (DREAL(UP).EQ.DREAL(UPBP)) THEN
	      IF (UP.EQ.UPBP) THEN
	        URVALUE = (0.0D0,0.0D0)
	      ELSE IF (DIMAG(UP).GT.DIMAG(UPBP)) THEN
	        IF ( DREAL(UPREF) .GT. DREAL(UPBP) ) THEN
	          URVALUE = CDSQRT ( UPBP2 - UP2 )
	          IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE   ! cdsqrt12
	          URVALUE = U0 * URVALUE
	        ELSE
	          URVALUE = U0 * CDSQRT ( UPBP2 - UP2 )
	        END IF
	      ELSE
	        URVALUE = CDSQRT ( UPBP2 - UP2 )
	        IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      END IF
	    ELSE IF (DREAL(UP).EQ.DREAL(-UPBP)) THEN
	      IF (UP.EQ.-UPBP) THEN
	        URVALUE = (0.0D0,0.0D0)
	      ELSE IF (DIMAG(UP).LT.DIMAG(-UPBP)) THEN
	        IF ( DREAL(UPREF) .LT. DREAL(-UPBP) ) THEN
	          URVALUE = CDSQRT ( UPBP2 - UP2 )
	          IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE   ! cdsqrt12
	          URVALUE = U0 * URVALUE
	        ELSE
	          URVALUE = U0 * CDSQRT ( UPBP2 - UP2 )
	        END IF
	      ELSE
	        URVALUE = CDSQRT ( UPBP2 - UP2 )
	        IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      END IF
	    ELSE
	      IF ( DABS(DREAL(UP)) .GT. DABS(DREAL(UPBP)) ) THEN
	        URVALUE = CDSQRT ( UPBP2 - UP2 )
	        IF (DIMAG(URVALUE).LT.0.0D0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      ELSE
	        URVALUE = U0 * CDSQRT ( UPBP2 - UP2 )
	      END IF
	    END IF
	  END IF
	ELSE                                    ! Muir approx from vertical
	  Z2 = UP2 / UPBP2
	  URVALUE = U0 * UPBP * MUIR(-ICUTTYP,Z2)
	END IF

	RETURN
	END








	COMPLEX*16 FUNCTION MUIR (IDEG, Z2)
C	from Claerbout "Imaging the Earth's Interior" p. 84

	IMPLICIT NONE

	INTEGER IDEG
	COMPLEX*16 Z2

	IF (IDEG.EQ.5) THEN
	  MUIR = (1.0D0,0.0D0)
	ELSE IF (IDEG.EQ.15) THEN
	  MUIR = 1 - Z2/2
	ELSE IF (IDEG.EQ.45) THEN
	  MUIR = 1 - Z2/(2-Z2/2)
	ELSE IF (IDEG.EQ.60) THEN
	  MUIR = 1 - Z2/(2-Z2/(2-Z2/2))
	ELSE
	  STOP 160
	END IF

	RETURN
	END








	SUBROUTINE PREDHCOMPAS
C	to compute exponents and factors for "dhcomp.for with idhtyp=1,2"
C	predhcompas must be called anew if the frequency has been changed
C	note: rho is always assumed to be continuous within each layer
C	      for technical reasons rho is here allowed to be discontinuous 
C	        across the source level (the asymptotics are adapted to this 
C	        case rather than to the underlying "ideal" case for which 
C	        ua2sourc,rhosourc and fsourc,rw0,rw were computed)

	IMPLICIT NONE

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

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INPARAM.INC'
	INCLUDE 'PREDHAS.INC'
	INCLUDE 'SOURCFW.INC'

	INTEGER IH,IH0, ID,IDOLD, I
	REAL*8 RHONSO,RHONSOP1, RHOA,RHOB, RHOAQ,RHOBQ, SLASK
	COMPLEX*16 ASFAKT0, CSLASK, BSLASK

	IF (KFWTYP.NE.4) STOP 3415564   ! not implemented
	IF (NSOURCFW.NE.0) STOP 3415565 ! not implemented

	CSLASK = OMEG * U0
	IH0 = NHYD + 1
	DO IH = 1, NHYD
	  ID = IDHYD(IH)
	  IF ( ID.LT.NSS .OR. ID.GT.NSW ) STOP 3415564   ! not implemented
	  IF (IH0.GT.NHYD) THEN
	    IF (ID.GE.NSO) IH0 = IH
	  END IF
	  ASEXP(IH) = CSLASK * DABS(ZD(ID)-ZD(NSO))
	END DO

	IF (NSS.LT.NSO) THEN
	  IF (IFSMTYP(NSO).LE.0) THEN
	    RHONSO = RHO(NSO)
	  ELSE
	    RHONSO = RHOSOURC
	  END IF
	ELSE
	  IF (IFSMTYP(NSO+1).LE.0) THEN
	    RHONSO = RHO(NSO+1)
	  ELSE
	    RHONSO = RHOSOURC
	  END IF
	END IF
	IF (NSO.LT.NSW) THEN
	  IF (IFSMTYP(NSO+1).LE.0) THEN
	    RHONSOP1 = RHO(NSO+1)
	  ELSE
	    RHONSOP1 = RHOSOURC
	  END IF
	ELSE
	  IF (NSO.EQ.NSS) THEN
	    RHONSOP1 = 1.0D0     ! now nss=nso=nsw, "fictitious fluid layer"
	  ELSE IF (IFSMTYP(NSO).LE.0) THEN
	    RHONSOP1 = RHO(NSO)
	  ELSE
	    RHONSOP1 = RHOSOURC
	  END IF
	END IF
	IF (RHONSO.EQ.RHONSOP1) THEN
	  ASFAKT0 = RW0 * (0.5D0*RHONSO) * OMEG
	ELSE
	  ASFAKT0 = RW0 * ( 1 / ( 1/RHONSO + 1/RHONSOP1 ) ) * OMEG
	END IF

	SLASK = 1.0D0
	RHOA = RHONSOP1
	IF (NSO.LT.NSW) THEN
	  IF (IFSMTYP(NSO+1).LE.0) THEN
	    RHOB = RHOA
	  ELSE
	    CALL FSUBMODEL ( IFSMTYP(NSO+1),ZD(NSO+1), CSLASK,BSLASK,RHOB )
	  END IF
	ELSE
	  IF (NSO.EQ.NSS) THEN
	    RHOB = 1.0D0         ! now nss=nso=nsw, "fictitious fluid layer"
	  ELSE IF (IFSMTYP(NSO).LE.0) THEN
	    RHOB = RHOA
	  ELSE
	    CALL FSUBMODEL ( IFSMTYP(NSO),ZD(NSO), CSLASK,BSLASK,RHOB )
	  END IF
	END IF
	IDOLD = NSO
	DO IH = IH0, NHYD
	  ID = IDHYD(IH)
	  DO I = IDOLD+1, ID
	    IF (I.LT.NSW) THEN
	      IF (IFSMTYP(I+1).LE.0) THEN
	        RHOAQ = RHO(I+1)
	        RHOBQ = RHOAQ
	      ELSE
	        CALL FSUBMODEL ( IFSMTYP(I+1),ZD(I), CSLASK,BSLASK,RHOAQ )
	        CALL FSUBMODEL ( IFSMTYP(I+1),ZD(I+1), CSLASK,BSLASK,RHOBQ )
	      END IF
	      SLASK = SLASK * DSQRT(RHOA/RHOB) * 2*RHOB / (RHOB+RHOAQ)
	      RHOA = RHOAQ
	      RHOB = RHOBQ
	    ELSE IF ( RHO(NSW+1).GT.0.0D0 .AND. ISOLID(NSW+1).EQ.0 ) THEN
	      IF (IFSMTYP(I+1).LE.0) THEN
	        RHOAQ = RHO(I+1)
	        RHOBQ = RHOAQ
	      ELSE
	        CALL FSUBMODEL ( IFSMTYP(I+1),ZD(I), CSLASK,BSLASK,RHOAQ )
	        CALL FSUBMODEL ( IFSMTYP(I+1),ZD(I+1), CSLASK,BSLASK,RHOBQ )
	      END IF
	      SLASK = SLASK * DSQRT(RHOA/RHOB) * 2*RHOB / (RHOB+RHOAQ)
	      RHOA = RHOAQ
	      RHOB = RHOBQ
	    ELSE
	      SLASK = SLASK * DSQRT(RHOA/RHOB)
	    END IF
	  END DO
	  IF (ID.LT.NSW) THEN
	    ASFAKT(IH) = ( SLASK*RHOA/RHONSOP1 ) * ASFAKT0
	  ELSE IF ( RHO(NSW+1).GT.0.0D0 .AND. ISOLID(NSW+1).EQ.0 ) THEN
	    ASFAKT(IH) = ( SLASK*RHOA/RHONSOP1 ) * ASFAKT0
	  ELSE
	    ASFAKT(IH) = ( 2*SLASK*RHOB/RHONSOP1 ) * ASFAKT0
	  END IF
	  IDOLD = ID
	END DO

	SLASK = 1.0D0
	RHOBQ = RHONSO
	IF (NSS.LT.NSO) THEN
	  IF (IFSMTYP(NSO).LE.0) THEN
	    RHOAQ = RHOBQ
	  ELSE
	    CALL FSUBMODEL ( IFSMTYP(NSO),ZD(NSO-1), CSLASK,BSLASK,RHOAQ )
	  END IF
	ELSE
	  IF (IFSMTYP(NSO+1).LE.0) THEN
	    RHOAQ = RHOBQ
	  ELSE
	    CALL FSUBMODEL ( IFSMTYP(NSO+1),ZD(NSO), CSLASK,BSLASK,RHOAQ )
	  END IF
	END IF
	IDOLD = NSO 
	DO IH = IH0-1, 1, -1
	  ID = IDHYD(IH)
	  DO I = IDOLD-1, ID, -1
	    IF (NSS.LT.I) THEN
	      IF (IFSMTYP(I).LE.0) THEN
	        RHOB = RHO(I)
	        RHOA = RHOB
	      ELSE
	        CALL FSUBMODEL ( IFSMTYP(I),ZD(I), CSLASK,BSLASK,RHOB )
	        CALL FSUBMODEL ( IFSMTYP(I),ZD(I-1), CSLASK,BSLASK,RHOA )
	      END IF
	      SLASK = SLASK * DSQRT(RHOBQ/RHOAQ) * 2*RHOAQ / (RHOB+RHOAQ) 
	      RHOBQ = RHOB
	      RHOAQ = RHOA
	    ELSE
	      SLASK = SLASK * DSQRT(RHOBQ/RHOAQ)
	    END IF
	  END DO
	  IF (NSS.LT.ID) THEN
	    ASFAKT(IH) = ( SLASK*RHOB/RHONSO ) * ASFAKT0
	  ELSE IF (ID.EQ.0) THEN
	    IF (RHO(0).NE.0.0D0) STOP 45166   ! only free surface implemented
	    ASFAKT(IH) = (0.0D0,0.0D0)           ! free surface
	  ELSE
	    ASFAKT(IH) = ( 2*SLASK*RHOA/RHONSO ) * ASFAKT0
	  END IF
	  IDOLD = ID
	END DO

	RETURN
	END







	SUBROUTINE FSUBLIQ2 ( N,ZLOK,YP, DYP )
C	for the linear differential-equation system for depth propagation 
C	of  (ww,pp) = (yp(1),yp(2))  in a fluid

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'         ! to get omeg
	INCLUDE 'FSUBCOM.INC'         ! to get rk2gomeg and ifsmtyplok

	INTEGER N              ! not used, but should be 2
	REAL*8 ZLOK            ! depth in km
	COMPLEX*16 YP(2)       ! ww=yp(1) pp=yp(2)

	COMPLEX*16 DYP(2)

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOK,UB2LOK

	CALL FSUBMODEL ( IFSMTYPLOK,ZLOK, UA2LOK,UB2LOK,RHOLOK )
C	ub2 irrelevant (fluid here)

	DYP(1) = ( OMEG*UA2LOK - RK2GOMEG ) * YP(2) / RHOLOK
	DYP(2) = -OMEG*RHOLOK * YP(1)

	RETURN
	END







	SUBROUTINE FSUBNG2 ( N,ZLOK,YT, DYT )
C	for the linear differential-equation system for depth propagation 
C	as  " P(D,z) transposed "  of  (t1,t3) in a fluid
C	P(z2,z1) is the 2x2 propagator in the fluid
C	in relation to fsubliq2 the system matrix is negated and transposed

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'         ! to get omeg
	INCLUDE 'FSUBCOM.INC'         ! to get rk2gomeg and ifsmtyplok

	INTEGER N              ! not used, but should be 2
	REAL*8 ZLOK            ! depth in km
	COMPLEX*16 YT(2)       ! t1=yt(1) t3=yt(2)  OBS

	COMPLEX*16 DYT(2)

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOK,UB2LOK

	CALL FSUBMODEL ( IFSMTYPLOK,ZLOK, UA2LOK,UB2LOK,RHOLOK )
C	ub2 irrelevant (fluid here)

	DYT(1) = OMEG*RHOLOK * YT(2)
	DYT(2) = ( -OMEG*UA2LOK + RK2GOMEG ) * YT(1) / RHOLOK

	RETURN
	END









	SUBROUTINE FSUBHARK5PRE ( N,ZLOK,YP, DYP )
C	for the linear differential-equation system for depth propagation 
C	as  " B(z,H) "  of (t1,t2,t3,t4,t6) in a solid
C	B(z2,z1) is the 6x6 propagator for the 2x2 minors
C	cf. fsubhark5

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'    ! to get omeg
	INCLUDE 'FSUBCOM.INC'    ! to get ucop,u2cop , rk2gomeg and ifsmtyplok

	INTEGER N         ! not used, but should be 5
	REAL*8 ZLOK       ! depth in km
	COMPLEX*16 YP(5)  ! p1=yp(1) p2=yp(2) p3=yp(3) p4=yp(4) p6=yp(5)  OBS
C	                                                      p5=-p2=-yp(2)

	COMPLEX*16 DYP(5)

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOKQ,UB2LOKQ, UA2GUB2, C1,C2, YP2Q

	CALL FSUBMODEL ( IFSMTYPLOK,ZLOK, UA2LOKQ,UB2LOKQ,RHOLOK )
	UA2LOKQ = UA2LOKQ/RHOLOK
	UB2LOKQ = UB2LOKQ/RHOLOK

	UA2GUB2 = UA2LOKQ/UB2LOKQ
	C1 = (1-2*UA2GUB2)*UCOP
	C2 = RHOLOK - 4*(1-UA2GUB2)*U2COP/UB2LOKQ

	YP2Q = YP(2)+YP(2)  ! "temporary doubling" as in harksi  ("p5=-p2")

	DYP(1) = OMEG * ( UA2LOKQ*YP(3) - UB2LOKQ*YP(4) )
	DYP(2) = OMEG * ( C1*YP(3) + UCOP*YP(4) )
	DYP(3) = OMEG * ( -RHOLOK*YP(1) - UCOP*YP2Q + UB2LOKQ*YP(5) )
	DYP(4) = OMEG * ( C2*YP(1) - C1*YP2Q - UA2LOKQ*YP(5) )
	DYP(5) = OMEG * ( -C2*YP(3) + RHOLOK*YP(4) )

	RETURN
	END







	SUBROUTINE FSUBHARK5 ( N,ZLOK,YT, DYT )
C	for the linear differential-equation system for depth propagation 
C	as  " B(D,z) transposed "  of (t1,t2,t3,t4,t6) in a solid
C	B(z2,z1) is the 6x6 propagator for the 2x2 minors
C	in relation to fsubhark5pre the system matrix is negated and transposed

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'    ! to get omeg
	INCLUDE 'FSUBCOM.INC'    ! to get ucop,u2cop , rk2gomeg and ifsmtyplok

	INTEGER N         ! not used, but should be 5
	REAL*8 ZLOK       ! depth in km
	COMPLEX*16 YT(5)  ! t1=yt(1) t2=yt(2) t3=yt(3) t4=yt(4) t6=yt(5)  OBS
C	                                                      t5=-t2=-yt(2)

	COMPLEX*16 DYT(5)

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOKQ,UB2LOKQ, UA2GUB2, C1,C2, YT2Q

	CALL FSUBMODEL ( IFSMTYPLOK,ZLOK, UA2LOKQ,UB2LOKQ,RHOLOK )
	UA2LOKQ = UA2LOKQ/RHOLOK
	UB2LOKQ = UB2LOKQ/RHOLOK

	UA2GUB2 = UA2LOKQ/UB2LOKQ
	C1 = (1-2*UA2GUB2)*UCOP
	C2 = RHOLOK - 4*(1-UA2GUB2)*U2COP/UB2LOKQ

	YT2Q = YT(2)+YT(2)  ! "temporary doubling" as in harksi  ("t5=-t2")

	DYT(1) = OMEG * ( RHOLOK*YT(3) - C2*YT(4) )
	DYT(2) = OMEG * ( UCOP*YT(3) + C1*YT(4) )
	DYT(3) = OMEG * ( -UA2LOKQ*YT(1) - C1*YT2Q + C2*YT(5) )
	DYT(4) = OMEG * ( UB2LOKQ*YT(1) - UCOP*YT2Q - RHOLOK*YT(5) )
	DYT(5) = OMEG * ( -UB2LOKQ*YT(3) + UA2LOKQ*YT(4) )

	RETURN
	END






	SUBROUTINE FSUBSOL4 ( N,ZLOK,YS, DYS )
C	for the linear differential-equation system for depth propagation 
C	of  ( omeg*r1, omeg*r2, r3, r4 )  in a solid

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'    ! to get omeg
	INCLUDE 'FSUBCOM.INC'    ! to get ucop,u2cop , rk2gomeg and ifsmtyplok

	INTEGER N              ! not used, but should be 4
	REAL*8 ZLOK            ! depth in km
	COMPLEX*16 YS(4)       ! "omeg*r1"=ys(1) "omeg*r2"=ys(2) 
C	                                r3=ys(3)        r4=ys(4)

	COMPLEX*16 DYS(4)

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOK,UB2LOK, UA2GUB2, C1U,C2

	CALL FSUBMODEL ( IFSMTYPLOK,ZLOK, UA2LOK,UB2LOK,RHOLOK )

	UA2GUB2 = UA2LOK/UB2LOK
	C1U = ( 1 - 2*UA2GUB2 ) * UCOP
	C2 = (4*RHOLOK) * (1-UA2GUB2) * U2COP / UB2LOK

	DYS(1) = OMEG * ( UCOP*YS(2) + UB2LOK*YS(3)/RHOLOK )
	DYS(2) = OMEG * ( -C1U*YS(1) + UA2LOK*YS(4)/RHOLOK )
	DYS(3) = OMEG * ( (C2-RHOLOK)*YS(1) - C1U*YS(4) )
	DYS(4) = OMEG * ( -RHOLOK*YS(2) - UCOP*YS(3) )

	RETURN
	END







	SUBROUTINE FSUBZSOL4 ( N,ZLOK,YZ, DYZ )
C	for the linear differential-equation system for depth propagation 
C	as  " P(D,z) transposed "  of  ( omeg*r1, omeg*r2, r3, r4 )  in a solid
C	P(z2,z1) is the 4x4 propagator in the solid
C	in relation to fsubsol4 the system matrix is negated and transposed

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'    ! to get omeg
	INCLUDE 'FSUBCOM.INC'    ! to get ucop,u2cop , rk2gomeg and ifsmtyplok

	INTEGER N              ! not used, but should be 4
	REAL*8 ZLOK            ! depth in km
	COMPLEX*16 YZ(4)       ! "omeg*r1"=yz(1) "omeg*r2"=yz(2) 
C	                                r3=yz(3)        r4=yz(4)

	COMPLEX*16 DYZ(4)

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOK,UB2LOK, UA2GUB2, C1U,C2

	CALL FSUBMODEL ( IFSMTYPLOK,ZLOK, UA2LOK,UB2LOK,RHOLOK )

	UA2GUB2 = UA2LOK/UB2LOK
	C1U = ( 1 - 2*UA2GUB2 ) * UCOP
	C2 = (4*RHOLOK) * (1-UA2GUB2) * U2COP / UB2LOK

	DYZ(1) = OMEG * ( C1U*YZ(2) - (C2-RHOLOK)*YZ(3) )
	DYZ(2) = OMEG * ( -UCOP*YZ(1) + RHOLOK*YZ(4) )
	DYZ(3) = OMEG * ( -UB2LOK*YZ(1)/RHOLOK + UCOP*YZ(4) )
	DYZ(4) = OMEG * ( -UA2LOK*YZ(2)/RHOLOK + C1U*YZ(3) )

	RETURN
	END





	SUBROUTINE ROT3 (YA,YB,CFAKT,BCFAKT,BCFAKTINV)
C	rotate +-pi/3 (or 0) , in order to maximize min(alna,alnb) from 
C	" call airy4e (ya, alna,blna, aiea,aidea,biea,bidea) "
C	" call airy4e (yb, alnb,blnb, aieb,aideb,bieb,bideb) "

	IMPLICIT NONE

	COMPLEX*16 YA,YB,CFAKT,BCFAKT,BCFAKTINV

	INTEGER II
	REAL*8 ALN(2)
	COMPLEX*16 Y(2),CZET(2), CSLASK1,CSLASK2, C120,C120N

	DATA C120  /(-0.5D0,+0.8660254037844387)/    ! cdexp(+i*pi/3)
	DATA C120N /(-0.5D0,-0.8660254037844387)/    ! cdexp(-i*pi/3)

	IF ( DIMAG(YA).NE.0.0D0 .OR. DIMAG(YB).NE.0.0D0 ) THEN
	  Y(1) = YA
	  Y(2) = YB
	  CZET(1) = Y(1) * CDSQRT(Y(1))
	  CZET(2) = Y(2) * CDSQRT(Y(2))
	  ALN(1) = DREAL(CZET(1))
	  ALN(2) = DREAL(CZET(2))
	  IF ( DABS(ALN(1)) .GE. DABS(ALN(2)) ) THEN
	    II = 1
	  ELSE
	    II = 2
	  END IF
	  IF ( ALN(II) .LT. 0.0D0 ) THEN
C	    this is when we want to rotate
	    IF ( DIMAG(Y(II)) .GE. 0.0D0 ) THEN
	      CSLASK1 = C120N
	      CSLASK2 = C120
	    ELSE
	      CSLASK1 = C120
	      CSLASK2 = C120N
	    END IF
	    YA = CSLASK1 * YA
	    YB = CSLASK1 * YB
	    CFAKT = CSLASK1 * CFAKT
	    BCFAKT = CSLASK1 * BCFAKT
	    BCFAKTINV = CSLASK2 * BCFAKTINV
	  END IF
	END IF

	RETURN
	END