! copyright (c) 2008 Swedish Defence Research Agency (FOI) 
! may not be copied or distributed without written permission from FOI



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

C	Sven Ivansson  1990






	COMPLEX*32 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*16 AIMCUT
	COMPLEX*32 U0,UPBP,UPBP2,UPREF, UP,UP2

	COMPLEX*32 Z2, MUIR

	IF ( REAL(UPBP).LE.0.0Q0 .OR. AIMAG(UPBP).LT.0.0Q0 ) 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.0Q0,0.0Q0)
	  ELSE
	    URVALUE = SQRT ( UPBP2 - UP2 )
	    IF (AIMAG(URVALUE).LT.0.0Q0) 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.0Q0,0.0Q0)
	  ELSE IF (REAL(UP).EQ.REAL(UPBP)) THEN
	    IF (UP.EQ.UPBP) THEN
	      URVALUE = (0.0Q0,0.0Q0)
	    ELSE IF (AIMAG(UP).GT.AIMAG(UPBP)) THEN
	      IF ( REAL(UPREF) .GT. REAL(UPBP) ) THEN
	        URVALUE = SQRT ( UPBP2 - UP2 )
	        IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      ELSE
	        URVALUE = U0 * SQRT ( UPBP2 - UP2 )
	      END IF
	    ELSE
	      URVALUE = SQRT ( UPBP2 - UP2 )
	      IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    END IF
	  ELSE IF (REAL(UP).EQ.REAL(-UPBP)) THEN
	    IF (UP.EQ.-UPBP) THEN
	      URVALUE = (0.0Q0,0.0Q0)
	    ELSE IF (AIMAG(UP).LT.AIMAG(-UPBP)) THEN
	      IF ( REAL(UPREF) .LT. REAL(-UPBP) ) THEN
	        URVALUE = SQRT ( UPBP2 - UP2 )
	        IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      ELSE
	        URVALUE = U0 * SQRT ( UPBP2 - UP2 )
	      END IF
	    ELSE
	      URVALUE = SQRT ( UPBP2 - UP2 )
	      IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    END IF
	  ELSE
	    IF ( ABS(REAL(UP)) .GT. ABS(REAL(UPBP)) ) THEN
	      URVALUE = SQRT ( UPBP2 - UP2 )
	      IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    ELSE
	      URVALUE = U0 * SQRT ( UPBP2 - UP2 )
	    END IF
	  END IF
	ELSE IF (ICUTTYP.EQ.3) THEN             ! "mixed" cut for up
	  IF ( ABS(AIMAG(UP)) .GT. AIMCUT ) THEN
C	    as for "icuttyp.eq.1"
	    IF (UP2.EQ.UPBP2) THEN
	      URVALUE = (0.0Q0,0.0Q0)
	    ELSE
	      URVALUE = SQRT ( UPBP2 - UP2 )
	      IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE       ! cdsqrt12
	      URVALUE = U0 * URVALUE
	    END IF
	  ELSE
C	    as for "icuttyp.eq.2"
	    IF (UP2.EQ.UPBP2) THEN
	      URVALUE = (0.0Q0,0.0Q0)
	    ELSE IF (REAL(UP).EQ.REAL(UPBP)) THEN
	      IF (UP.EQ.UPBP) THEN
	        URVALUE = (0.0Q0,0.0Q0)
	      ELSE IF (AIMAG(UP).GT.AIMAG(UPBP)) THEN
	        IF ( REAL(UPREF) .GT. REAL(UPBP) ) THEN
	          URVALUE = SQRT ( UPBP2 - UP2 )
	          IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE   ! cdsqrt12
	          URVALUE = U0 * URVALUE
	        ELSE
	          URVALUE = U0 * SQRT ( UPBP2 - UP2 )
	        END IF
	      ELSE
	        URVALUE = SQRT ( UPBP2 - UP2 )
	        IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      END IF
	    ELSE IF (REAL(UP).EQ.REAL(-UPBP)) THEN
	      IF (UP.EQ.-UPBP) THEN
	        URVALUE = (0.0Q0,0.0Q0)
	      ELSE IF (AIMAG(UP).LT.AIMAG(-UPBP)) THEN
	        IF ( REAL(UPREF) .LT. REAL(-UPBP) ) THEN
	          URVALUE = SQRT ( UPBP2 - UP2 )
	          IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE   ! cdsqrt12
	          URVALUE = U0 * URVALUE
	        ELSE
	          URVALUE = U0 * SQRT ( UPBP2 - UP2 )
	        END IF
	      ELSE
	        URVALUE = SQRT ( UPBP2 - UP2 )
	        IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      END IF
	    ELSE
	      IF ( ABS(REAL(UP)) .GT. ABS(REAL(UPBP)) ) THEN
	        URVALUE = SQRT ( UPBP2 - UP2 )
	        IF (AIMAG(URVALUE).LT.0.0Q0) URVALUE = -URVALUE     ! cdsqrt12
	        URVALUE = U0 * URVALUE
	      ELSE
	        URVALUE = U0 * SQRT ( 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*32 FUNCTION MUIR (IDEG, Z2)
C	from Claerbout "Imaging the Earth's Interior" p. 84

	IMPLICIT NONE

	INTEGER IDEG
	COMPLEX*32 Z2

	IF (IDEG.EQ.5) THEN
	  MUIR = (1.0Q0,0.0Q0)
	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)

	USE INMODELa
	USE PREDHASa

	IMPLICIT NONE

	STOP '%%% dhsub, SHf'

	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*16 ZLOK            ! depth in km
	COMPLEX*32 YP(2)       ! ww=yp(1) pp=yp(2)

	COMPLEX*32 DYP(2)

	STOP '%%% dhsub, SHf'

	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*16 ZLOK            ! depth in km
	COMPLEX*32 YT(2)       ! t1=yt(1) t3=yt(2)  OBS

	COMPLEX*32 DYT(2)

	STOP '%%% dhsub, SHf'

	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*16 ZLOK       ! depth in km
	COMPLEX*32 YP(5)  ! p1=yp(1) p2=yp(2) p3=yp(3) p4=yp(4) p6=yp(5)  OBS
C	                                                      p5=-p2=-yp(2)

	COMPLEX*32 DYP(5)

	STOP '%%% dhsub, SHf'

	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*16 ZLOK       ! depth in km
	COMPLEX*32 YT(5)  ! t1=yt(1) t2=yt(2) t3=yt(3) t4=yt(4) t6=yt(5)  OBS
C	                                                      t5=-t2=-yt(2)

	COMPLEX*32 DYT(5)

	STOP '%%% dhsub, SHf'

	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*16 ZLOK            ! depth in km
	COMPLEX*32 YS(4)       ! "omeg*r1"=ys(1) "omeg*r2"=ys(2) 
C	                                r3=ys(3)        r4=ys(4)

	COMPLEX*32 DYS(4)

	STOP '%%% dhsub, SHf'

	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*16 ZLOK            ! depth in km
	COMPLEX*32 YZ(4)       ! "omeg*r1"=yz(1) "omeg*r2"=yz(2) 
C	                                r3=yz(3)        r4=yz(4)

	COMPLEX*32 DYZ(4)

	STOP '%%% dhsub, SHf'

	RETURN
	END

