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



C	subroutine "dhcomp etc" for rpress.f etc
C	Sven Ivansson  1990




	SUBROUTINE CDETSUBCIRC (VX,VY, RLNDET,CDET)

	USE INMODELa

	IMPLICIT NONE

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

	INCLUDE '../CIRCDATA.INC'  ! up0circ and radiecirc

	REAL*16 VX,VY, PX,PY

	REAL*16 RLNDET
	COMPLEX*32 CDET

	COMPLEX*32 UP, QCMPLX
	COMPLEX*32, DIMENSION(:), ALLOCATABLE :: DHFAKT
	LOGICAL DIMDONE
	SAVE DHFAKT, DIMDONE
	DATA  DIMDONE/.FALSE./

	UP = UP0CIRC + 
     $	(RADIECIRC*EXP(-VY)) * QCMPLX(COS(VX),SIN(VX))
	GOTO 100

	ENTRY CDETSUB (PX,PY, RLNDET,CDET)
	UP = QCMPLX(PX,PY)

 100	IF (.NOT.DIMDONE) THEN
	  ALLOCATE ( DHFAKT(nhyddim) )
	  DIMDONE = .TRUE.
	END IF

	CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)

	RETURN
	END





	SUBROUTINE DHCOMPCIRC (VP, RLNDET,CDET, DHFAKT)
C	for complex pressure frequency-response in "(10**15 N)/(km**2)" 
C	integrand factors for the respective depths and distances 
C	( "exp(-i*omeg*t)" with Re(omeg)>0.0 is behind )
C	point source assumed with diagonal moment-tensor matrix with 
C	diagonal elements  " exp(-i*omeg*t) "  "(10**15 N)*(km)"
C	corresponding to the volume-force 
C	 " - exp(-i*omeg*t) * grad(Diracfunction) "  "(10**15 N)/(km**3)"
C	so that the pressure in a punctured disc around the source point is 
C	 " - ( (omeg*omeg*ua*ua)/(4*pi) ) * exp(-i*omeg*(t-r*ua)) / r "
C	              "(10**15 N)/(km**2)"  approximately
C	where r is spherical distance (!) in km, omeg is in Hz, ua in s/km, 
C	and t in s   %% note that "(10**15 N)" = "(10**12 kg)*(km/(s**2))" %%
C	%% note that Im(ua) >= 0.0 when Re(omeg) > 0.0 and that 
C	ua(-omeg) is the complex conjugate of ua(omeg) , of course the 
C	corresponding statements are true for the shear slowness %%

C	idhtyp:   0 for "usual , 1 for "asymptotic" , 2 for "usual-asymptotic"

	USE INMODELa
	USE PREDHASa

	IMPLICIT NONE

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

	INCLUDE '../INPARAM.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE '../HARKLIQ.INC'
	INCLUDE '../CIRCDATA.INC'   ! up0circ and radiecirc
	INCLUDE '../INSPEC.INC'

	COMPLEX*32 VP,UPDUMMY

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

	INTEGER IHARK0, IH
	REAL*16 TLNL,TLNU, RHOCDET
	COMPLEX*32 DUPGDVP, UP, U,RK,UP2,U2,RK2, YPL1,YPL2, YPU1,YPU2, 
     $	RW, URA,URB, URVALUE, QCMPLX
	LOGICAL DUPGDVPCORR

	DUPGDVPCORR = .TRUE.
	DUPGDVP = (RADIECIRC*EXP(-AIMAG(VP))) * 
     $	QCMPLX(-SIN(REAL(VP)),COS(REAL(VP)))   ! i*radiecirc*cdexp(i*vp)
	UP = UP0CIRC + QCMPLX(AIMAG(DUPGDVP),-REAL(DUPGDVP))
	GOTO 100

	ENTRY DHCOMP (UPDUMMY, RLNDET,CDET, DHFAKT)
	DUPGDVPCORR = .FALSE.
	UP = UPDUMMY

C       this "dhcomp" is not for "fw" but 'nsourcfw>0' ok from multi

 100	IF (IDHTYP.NE.1) THEN

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

	  UCOP = U
	  U2COP = U2
	  RK2GOMEG = RK2/OMEG

	  RHOCDET = 1.0Q0

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

	  IF (RHO(NOL).GT.0.0Q0) THEN
C	    half-space present, "branch-cut(s) and sheet(s)" important
	    URA = URVALUE (ICUTTYP,AIMCUT, U0,UPBPA,UPBPA2,UPREF, 
     $	    UP,UP2)
	    IF ( ISHEETTYP.EQ.2 .OR. ISHEETTYP.EQ.4 ) URA = -URA
	    IF (ISOLID(NOL).NE.0) THEN
	      STOP '%%% dhcomp, SHf'
	    ELSE
C	      fluid half-space below
	      IF ( NOL .GT. (NSW+1) ) THEN
	        STOP '%%% dhcomp, SHf'
	      ELSE
	        TLNL = 0.0Q0
	        YPL1 = QCMPLX(AIMAG(URA),-REAL(URA)) !  -i*ura
	        YPL2 = RHO(NOL) * UA2(NOL) * REAL(UA2(NOL))
	      END IF
	    END IF
	  ELSE IF (RHO(NOL).EQ.-99.0Q0) THEN
	    STOP '%%% dhcomp, SHf'
	  ELSE IF (RHO(NOL).LT.0.0Q0) THEN
C	    rigid below
	    IF ( NOL .GT. (NSW+1) ) THEN
	      STOP '%%% dhcomp, SHf'
	    ELSE
	      TLNL = 0.0Q0
	      YPL1 = (0.0Q0,0.0Q0)
	      YPL2 = (1.0Q0,0.0Q0)
	    END IF
	  ELSE
C	    free below
	    IF ( NOL .GT. (NSW+1) ) THEN
	      STOP '%%% dhcomp, SHf'
	    ELSE
	      TLNL = 0.0Q0
	      YPL1 = (1.0Q0,0.0Q0)
	      YPL2 = (0.0Q0,0.0Q0)
	    END IF
	  END IF

	  IF (RHO(0).GT.0.0Q0) THEN
C	    half-space present, "branch-cut(s) and sheet(s)" important
	    URA = URVALUE (ICUTTYP,AIMCUT, U0,UPBPASURF,UPBPA2SURF,UPREF, 
     $	    UP,UP2)
	    IF ( ISHEETTYP.EQ.2 .OR. ISHEETTYP.EQ.4 ) URA = -URA
	    URA = -URA         ! up/down viceversa for above and below
	    IF (ISOLID(0).NE.0) THEN
	      STOP '%%% dhcomp, SHf'
	    ELSE
C	      fluid half-space above
	      IF (NSS.GT.0) THEN
	        STOP '%%% dhcomp, SHf'
	      ELSE
		TLNU = 0.0Q0
		YPU1 = QCMPLX(AIMAG(URA),-REAL(URA)) !  -i*ura
		YPU2 = RHO(0) * UA2(0) * REAL(UA2(0))
	        IF (NSW.EQ.0) THEN
	          STOP '%%% dhcomp, SHf'
	        END IF
	      END IF
	    END IF
	  ELSE IF (RHO(0).EQ.-99.0Q0) THEN
C	    adaptation from avh/multi/liqsi.f (cf articles/ua15bok/GF)
	    IF ( NSS.NE.0 .OR. NSW.LT.1 ) STOP '%%% dhcomp, "multi"'
	    IF ( IFSMTYP(1).LE.-4 .OR. IFSMTYP(1).GE.1 )
     &	     STOP '%%% dhcomp, "multi"'
	    STOP '%dhcomp, "multi" not yet impl for quadruple prec'
	  ELSE IF (RHO(0).LT.0.0Q0) THEN
C	    rigid above
	    IF (NSS.GT.0) THEN
	      STOP '%%% dhcomp, SHf'
	    ELSE IF (NSW.EQ.0) THEN
	      STOP '%%% dhcomp, SHf'
	    ELSE
	      TLNU = 0.0Q0
	      YPU1 = (0.0Q0,0.0Q0)
	      YPU2 = (1.0Q0,0.0Q0)
	    END IF
	  ELSE
C	    free above
	    IF (NSS.GT.0) THEN
	      STOP '%%% dhcomp, SHf'
	    ELSE IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      IF (IDET.NE.1) STOP 3337
	      RLNDET = TLNL
	      CDET = YPL2
	      RETURN
	    ELSE
	      TLNU = 0.0Q0
	      YPU1 = (1.0Q0,0.0Q0)
	      YPU2 = (0.0Q0,0.0Q0)
	    END IF
	  END IF

	  RW = RW0*RK    ! the appropriate discontinuity ("value below" minus 
C	  "value above") in the second component of the "Aki (!) displacement-
C	  stress vector (Aki & Richards 1980, 7.4.2) but with omeg* for the 
C	  first two components" at the source level is  -fsourc*omeg*rk  if 
C	  the factor 1/(2*pi) and the - from the Aki vector basis-functions 
C	  (rel. the Bessel functions) and the - for pressure rather than 
C	  stress are included along with the additional factor rk in the 
C	  integrand
C	  (usually we include in rw0 the factor rk0 for the change of 
C	  variable from k to up)
C	  with (r1,r2,r3,r4) as in Aki 7.4.2 our vector here is then  
C	       rk0 * (k/(2*pi)) * ( omeg*r1, omeg*r2, r3, r4 ) 
C	  and our (ww,pp) = (yp1,yp2) corresponds to the (second,fourth) 
C	  part 
C	  the part  -(k/(2*pi))*r4  is denoted "sigma" %%

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

C	  note that if "idet.ne.0" then  " dexp(rlndet) * cdet "  will be 
C	  an analytic function of up
C	  - with free and/or rigid boundaries this analytic function is entire
C	  - for a half-space boundary the function is analytic except at the 
C	    branch-cuts for "ura(surf)" (and "urb(surf)") for the half-space(s)

	ELSE

C	  asymptotic expressions valid for large dabs(dreal(rk))
C	  rlndet,cdet only given "default" values
C	  predhcompas must have been called in advance
	  RK = RK0 * UP
	  IF (IRKSPEC.EQ.0) THEN
	    IF (REAL(RK).GE.0.0Q0) THEN
	      DO IH = 1, NHYD
	        DHFAKT(IH) = ASFAKT(IH) * EXP(-UP*ASEXP(IH))
	      END DO
	    ELSE
	      DO IH = 1, NHYD
	        DHFAKT(IH) = - ASFAKT(IH) * EXP(UP*ASEXP(IH))
	      END DO
	    END IF
	  ELSE
	    DLNSPEC = ORMINSPEC * AIMAG(RK)
	    IF (REAL(RK).GE.0.0Q0) THEN
	      DO IH = 1, NHYD
	        DHFAKT(IH) = ASFAKT(IH) * EXP(-UP*ASEXP(IH)-DLNSPEC)
	      END DO
	    ELSE
	      DO IH = 1, NHYD
	        DHFAKT(IH) = - ASFAKT(IH) * EXP(UP*ASEXP(IH)-DLNSPEC)
	      END DO
	    END IF
	  END IF
	  RLNDET = 0.0Q0
	  CDET = (1.0Q0,0.0Q0)

	END IF

	IF (IDHTYP.EQ.2) THEN
C	  reduction with asymptotic expressions valid for large dabs(dreal(rk))
C	  predhcompas must have been called in advance
	  IF (IRKSPEC.EQ.0) THEN
	    IF (REAL(RK).GE.0.0Q0) THEN
	      DO IH = 1, NHYD
	        DHFAKT(IH) = DHFAKT(IH) - 
     $	        ASFAKT(IH)*EXP(-UP*ASEXP(IH))
	      END DO
	    ELSE
	      DO IH = 1, NHYD
	        DHFAKT(IH) = DHFAKT(IH) + 
     $	        ASFAKT(IH)*EXP(UP*ASEXP(IH))
	      END DO
	    END IF
	  ELSE
	    DLNSPEC = ORMINSPEC * AIMAG(RK)
	    IF (REAL(RK).GE.0.0Q0) THEN
	      DO IH = 1, NHYD
	        DHFAKT(IH) = DHFAKT(IH) - 
     $	        ASFAKT(IH)*EXP(-UP*ASEXP(IH)-DLNSPEC)
	      END DO
	    ELSE
	      DO IH = 1, NHYD
	        DHFAKT(IH) = DHFAKT(IH) + 
     $	        ASFAKT(IH)*EXP(UP*ASEXP(IH)-DLNSPEC)
	      END DO
	    END IF
	  END IF
	END IF

	IF (DUPGDVPCORR) THEN
	  DO IH = 1, NHYD
	    DHFAKT(IH) = DUPGDVP * DHFAKT(IH)
	  END DO
	END IF

C	pressure = "integral(a,+infinity)" ( -sigma(k)*BH(k*r)*(dk/dup)*dup )
C	  for each depth and each horizontal distance where either 
C	"a=0.0 and BH(z)=J0(z)" or "a=-infinity and BH(z)=0.5*H0(1)(z)"
C	%%  note that our "dhfakt" corresponds to " -sigma * dk/dup " and 
C	that our "rfakt" (see the rcomp routine) corresponds to " BH(k*r) " 
C	since " du/dup = u0 " and " dk/dup = rk0 " %%
C	%%  however, if (dupgdvpcorr) <i.e. "dhcompcirc"> our "dhfakt" 
C	corresponds to " -sigma * dk/dup * dup/dvp = -sigma * dk/dvp "
C	%% if tapering should be desired, suitable factors would be  
C	" 0.5q0 * ( 1 + cos((1-s)*pi) ) "  and   
C	" 0.5q0 * ( 1 + cos(s*pi) )     "  
C	with s given by  " up = (1-s)*uph1 + s*uph2 "  %%

	RETURN
	END
