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



C	subroutine liqsi for rpress.f etc
C	Sven Ivansson  1990





	SUBROUTINE LIQSI ( OMEG,OMEG2,U,U2,RK,RK2, RW0,RW, YPL1,YPL2, 
     $	YPU1,YPU2, IDET,RHOCDET,RLNDET,CDET, CLIQL,CLIQU, DHFAKT )
C	!!! ccc  ", dhomw,DHFAKT)"
C	complex "angular frequency" omeg with Re(omeg) > 0.0 ,  omeg2=omeg**2
C	complex "horizontal slowness" u , u2=u**2
C	complex "wave-number" rk=omeg*u , rk2=rk**2
C	propagator-matrix calculations, through the fluid layers, for the 
C	second and fourth components of the "Aki (!) displacement-stress 
C	vector (Aki & Richards 1980, 7.4.2) but with omeg* for the first 
C	two components", the factor 1/(2*pi) and the - from the Aki 
C	vector basis-functions (rel. the Bessel functions) are also 
C	included here
C	note that the input  (wwl,ppl) = (ypl1,ypl2)  and  (wwu,ppu) = 
C	(ypu1,ypu2)  are subsequently "propagated" and their output 
C	values are for the depth zd(nso)
C	"lower ww,pp" (for the depth zd(nsw)) will be cliql*"input wwl,ppl"
C	"upper ww,pp" (for the depth zd(nss)) will be cliqu*"input wwu,ppu"
C	the source is introduced at a certain depth ( as given by "nso", 
C	" nss <= nso <= nsw " )
C	note that rlndet is input/output (but only if idet.ne.0)

	USE INMODELa

	IMPLICIT NONE

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

	INCLUDE '../DPCFSUB.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE '../INSPEC.INC'
	INCLUDE '../DDULOK.INC'
	INCLUDE '../SOURCFW.INC'
	INCLUDE 'BEAMARR.INC'

	INTEGER IDET
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2, RW0,RW

	COMPLEX*16 YPL1,YPL2, YPU1,YPU2

	REAL*8 RHOCDET,RLNDET
	COMPLEX*16 CDET, CLIQL,CLIQU, DHFAKT(*)    !!! ",dhomw(*),DHFAKT(*)"

	INTEGER ID, IND, ISLASK
	REAL*8 RLNDETLOK, PLNREF, PMAX, ZENDFSUB, Z0FSUB,
     $	DPCLNFSUB, RWTLN, E1P,E1PROD, SLASK, PI, 
     $	PLNDELTA(noldim)
	COMPLEX*16 URA, URANEGI, CP,CPX, STP,SDP, 
     $	WWLOK, YA,YB, CFAKT,BCFAKT,BCFAKTINV, CDETLOK, 
     $	RWTOT,RTTOT, CSLASK1,CSLASK2, CSLASK, 
     $	YPL(2), YPU(2)
	LOGICAL EXPTYP

	EXTERNAL FSUBLIQ2

	DATA PI /3.141592653589793D0/

C	upwards now
C	necessary to store de-scaling exponents, this is done with plndelta
	IF ( YPL1.NE.(0.0D0,0.0D0) .OR. YPL2.NE.(0.0D0,0.0D0) ) THEN
	  DO ID = NSW, NSO+1, -1
	    ISLASK = IHDD(ID)
	    IF (ISLASK.GT.0) THEN
ccc	      dhomw(islask) = ypl1
	      DHFAKT(ISLASK) = YPL2
	    END IF
	    IF (IFSMTYP(ID).LE.-4) THEN
	      IF (IFSMTYP(ID).GE.-5) THEN
C	        warning: not allowed to skip for "first id" in non-empty loop
	        CALL LIQPART ( OMEG,OMEG2,U,U2, UA2(ID), DZ(ID), 
     $	        URA, EXPTYP, PLNREF, CP, STP,SDP )
	        IF (EXPTYP) URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))  ! -i*ura
	      END IF
	      YPL2 = YPL2 / RHO(ID)  !already here  !"1/rho" could be stored
	      IF (EXPTYP) THEN
C	        now the first propagation step
	        WWLOK = YPL1 / URANEGI
	        YPL1 = WWLOK + YPL2
	        YPL2 = -WWLOK + YPL2
C	        now the main propagation step, freq-dependent
	        PLNDELTA(ID) = PLNREF - DLOG2    ! -dlog2 from "first step"
	        YPL1 = DCONJG(CP) * YPL1            ! "upwards"
	        YPL2 = DREAL(SDP)*CP * YPL2         ! "upwards"
C	        now the last propagation step
	        WWLOK = YPL1 - YPL2
	        YPL2 = RHO(ID) * ( YPL1 + YPL2 )
	        YPL1 = URANEGI*WWLOK
	      ELSE
C	        now the main and only propagation step, freq-dependent
	        PLNDELTA(ID) = PLNREF
	        WWLOK = CP*YPL1 - STP*YPL2                ! "upwards"
	        YPL2 = RHO(ID) * ( SDP*YPL1 + CP*YPL2 )   ! "upwards"
	        YPL1 = WWLOK
	      END IF
	    ELSE IF (IFSMTYP(ID).LE.0) THEN
	      IF (IFSMTYP(ID).GE.-1) THEN
C	        warning: not allowed to skip for "first id" in non-empty loop
	        ISLASK = 0
	      ELSE
	        ISLASK = -1
	      END IF
	      CALL AIRYPART ( OMEG,OMEG2,U,U2,RK,RK2, 
     $	      UA2(ID),UB2(ID), DZ(ID), ISLASK, 
     $	      CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	      PLNREF, CP,CPX, STP,SDP )
	      PLNDELTA(ID) = PLNREF
	      YPL2 = YPL2 / RHO(ID)               ! "1/rho(id)" could be stored
	      WWLOK = PI * ( CPX*YPL1 - STP*YPL2 )
	      YPL2 = (PI*RHO(ID)) * ( SDP*YPL1 + CP*YPL2 ) !"pi*rho" could stor
	      YPL1 = WWLOK               ! "upwards"
	    ELSE
	      YPL(1) = YPL1
	      YPL(2) = YPL2
	      IFSMTYPLOK = IFSMTYP(ID)
	      ZENDFSUB = ZD(ID-1)
	      Z0FSUB = ZD(ID)
	      DPCLNFSUB = 0.0D0
	      CALL DPCLIN ( FSUBLIQ2,2, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	      HFAKTMAX,HOFF,EPSMAGNIFMAX, ZENDFSUB, Z0FSUB,DPCLNFSUB,
     $	      YPL,EPSFSUB,HSTEPFSUB, IND )
	      PLNDELTA(ID) = DPCLNFSUB
	      IF (IND.NE.0) THEN
	        WRITE(*,*) ' %% u,id ',UCOP,ID
	        STOP 6770
	      END IF
	      YPL1 = YPL(1)
	      YPL2 = YPL(2)
	    END IF
	    PMAX = DMAX1 ( DABS(DREAL(YPL1)), DABS(DIMAG(YPL1)), 
     $	    DABS(DREAL(YPL2)), DABS(DIMAG(YPL2)) )
	    IF ( PMAX.LT.tspmaxlim1 .OR. tspmaxlim2.LT.PMAX ) THEN
	      PMAX = 1/PMAX   ! if pmax was 0.0d0, due to num. canc.
	      PLNDELTA(ID) = PLNDELTA(ID) - DLOG(PMAX)
	      YPL1 = PMAX*YPL1
	      YPL2 = PMAX*YPL2
	    END IF
	  END DO
	ELSE
	  DO ID = NSW, NSO+1, -1
	    ISLASK = IHDD(ID)
	    IF (ISLASK.GT.0) THEN
ccc	      dhomw(islask) = (0.0d0,0.0d0)
	      DHFAKT(ISLASK) = (0.0D0,0.0D0)
	    END IF
	    PLNDELTA(ID) = 0.0D0
	  END DO
	  YPL1 = (0.0D0,0.0D0)
	  YPL2 = (0.0D0,0.0D0)
	END IF

C	downwards now
C	necessary to store de-scaling exponents, this is done with plndelta
	IF ( YPU1.NE.(0.0D0,0.0D0) .OR. YPU2.NE.(0.0D0,0.0D0) ) THEN
	  ISLASK = IHDD(NSS)
	  IF (ISLASK.GT.0) THEN
ccc	    dhomw(islask) = ypu1
	    DHFAKT(ISLASK) = YPU2
	  END IF
	  DO ID = NSS+1, NSO
	    IF (IFSMTYP(ID).LE.-4) THEN
	      IF ( IFSMTYP(ID).EQ.-7 .OR. IFSMTYP(ID).EQ.-4 ) THEN
C	        warning: not allowed to skip for "first id" in non-empty loop
	        CALL LIQPART ( OMEG,OMEG2,U,U2, UA2(ID), DZ(ID), 
     $	        URA, EXPTYP, PLNREF, CP, STP,SDP )
	        IF (EXPTYP) URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))  ! -i*ura
	      END IF
	      YPU2 = YPU2 / RHO(ID)  !already here   !"1/rho" could be stored
	      IF (EXPTYP) THEN
C	        now the first propagation step
	        WWLOK = YPU1 / URANEGI
	        YPU1 = WWLOK + YPU2
	        YPU2 = -WWLOK + YPU2
C	        now the main propagation step, freq-dependent
	        PLNDELTA(ID) = PLNREF - DLOG2    ! -dlog2 from "first step"
	        YPU1 = DREAL(SDP)*CP * YPU1
	        YPU2 = DCONJG(CP) * YPU2
C	        now the last propagation step
	        WWLOK = YPU1 - YPU2
	        YPU2 = RHO(ID) * ( YPU1 + YPU2 )
	        YPU1 = URANEGI*WWLOK
	      ELSE
C	        now the main and only propagation step, freq-dependent
	        PLNDELTA(ID) = PLNREF
	        WWLOK = CP*YPU1 + STP*YPU2
	        YPU2 = RHO(ID) * ( -SDP*YPU1 + CP*YPU2 )
	        YPU1 = WWLOK
	      END IF
	    ELSE IF (IFSMTYP(ID).LE.0) THEN
	      IF ( IFSMTYP(ID).EQ.-3 .OR. IFSMTYP(ID).EQ.0 ) THEN
C	        warning: not allowed to skip for "first id" in non-empty loop
	        ISLASK = 0
	      ELSE
	        ISLASK = 1
	      END IF
	      CALL AIRYPART ( OMEG,OMEG2,U,U2,RK,RK2, 
     $	      UA2(ID),UB2(ID), DZ(ID), ISLASK, 
     $	      CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	      PLNREF, CP,CPX, STP,SDP )
	      PLNDELTA(ID) = PLNREF
	      YPU2 = YPU2 / RHO(ID)               ! "1/rho(id)" could be stored
	      WWLOK = PI * ( CP*YPU1 + STP*YPU2 )
	      YPU2 = (PI*RHO(ID)) * ( -SDP*YPU1 + CPX*YPU2 ) !"pi*rho" could st
	      YPU1 = WWLOK
	    ELSE
	      YPU(1) = YPU1
	      YPU(2) = YPU2
	      IFSMTYPLOK = IFSMTYP(ID)
	      ZENDFSUB = ZD(ID)
	      Z0FSUB = ZD(ID-1)
	      DPCLNFSUB = 0.0D0
	      CALL DPCLIN ( FSUBLIQ2,2, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	      HFAKTMAX,HOFF,EPSMAGNIFMAX, ZENDFSUB, Z0FSUB,DPCLNFSUB,
     $	      YPU,EPSFSUB,HSTEPFSUB, IND )
	      PLNDELTA(ID) = DPCLNFSUB
	      IF (IND.NE.0) THEN
	        WRITE(*,*) ' %% u,id ',UCOP,ID
	        STOP 6770
	      END IF
	      YPU1 = YPU(1)
	      YPU2 = YPU(2)
	    END IF
	    PMAX = DMAX1 ( DABS(DREAL(YPU1)), DABS(DIMAG(YPU1)), 
     $	    DABS(DREAL(YPU2)), DABS(DIMAG(YPU2)) )
	    IF ( PMAX.LT.tspmaxlim1 .OR. tspmaxlim2.LT.PMAX ) THEN
	      PMAX = 1/PMAX   ! if pmax was 0.0d0, due to num. canc.
	      PLNDELTA(ID) = PLNDELTA(ID) - DLOG(PMAX)
	      YPU1 = PMAX*YPU1
	      YPU2 = PMAX*YPU2
	    END IF
	    ISLASK = IHDD(ID)
	    IF (ISLASK.GT.0) THEN
ccc	      dhomw(islask) = ypu1
	      DHFAKT(ISLASK) = YPU2
	    END IF
	  END DO
	ELSE
	  ISLASK = IHDD(NSS)
	  IF (ISLASK.GT.0) THEN
ccc	    dhomw(islask) = (0.0d0,0.0d0)
	    DHFAKT(ISLASK) = (0.0D0,0.0D0)
	  END IF
	  DO ID = NSS+1, NSO
	    ISLASK = IHDD(ID)
	    IF (ISLASK.GT.0) THEN
ccc	      dhomw(islask) = (0.0d0,0.0d0)
	      DHFAKT(ISLASK) = (0.0D0,0.0D0)
	    END IF
	    PLNDELTA(ID) = 0.0D0
	  END DO
	  YPU1 = (0.0D0,0.0D0)
	  YPU2 = (0.0D0,0.0D0)
	END IF

	CDETLOK = YPU1*YPL2 - YPU2*YPL1
	IF (IDET.NE.0) THEN
C	  the de-scaling of cdet that has to be done because of scaling within 
C	  liqsi is of the form  " exp ( sum(id=nss+1,nsw) plndelta(id) ) "
	  RLNDETLOK = RLNDET
	  DO ID = NSS+1, NSW
	    RLNDETLOK = RLNDETLOK + PLNDELTA(ID)
	  END DO
	  RLNDET = RLNDETLOK
	  PMAX = DMAX1 ( DABS(DREAL(CDETLOK)), DABS(DIMAG(CDETLOK)) )
	  IF ( PMAX.LT.tspmaxlim1 .OR. tspmaxlim2.LT.PMAX ) THEN
	    IF (PMAX.NE.0.0D0) THEN
	      RLNDET = RLNDET + DLOG(PMAX)
	      CDET = (RHOCDET/PMAX) * CDETLOK
	    ELSE
	      CDET = (0.0D0,0.0D0)
	    END IF
	  ELSE
	    CDET = RHOCDET * CDETLOK
	  END IF
	  IF (IDET.EQ.1) RETURN
	ELSE
	  CDET = CDETLOK      ! store the "direction" anyway
	END IF

	IF (CDETLOK.EQ.(0.0D0,0.0D0)) THEN
	  IF (RW0.NE.(0.0D0,0.0D0)) THEN
	    WRITE(*,*)
	    WRITE(*,*) ' %%% (LIQSI) NORMAL MODE ENCOUNTERED'
	    WRITE(*,*) ' freq ',DREAL(OMEG)/(2*PI)
	    WRITE(*,*) ' omeg ',OMEG
	    WRITE(*,*) ' rk ',RK
	    WRITE(*,*) ' u ',U
	    WRITE(*,*) ' ypl1 ',YPL1
	    WRITE(*,*) ' ypl2 ',YPL2
	    WRITE(*,*) ' ypu1 ',YPU1
	    WRITE(*,*) ' ypu2 ',YPU2
	    STOP 320
	  END IF
	END IF
	IF (RW0.NE.(0.0D0,0.0D0)) THEN
	  IF (IDET.NE.-1) THEN
C	    the discontinuity ("value below" minus "value above"), in the 
C	    second component of the "Aki (!) displacement-stress vector but 
C	    with omeg* for the first two components" at the source level is 
C	    rw, some relevant factors are also included in rw
	    IF (NSARR.EQ.0) THEN
	      CSLASK = -RW / CDETLOK
	      CSLASK1 = YPU2 * CSLASK
	      CSLASK2 = YPL2 * CSLASK
	      RWTLN = 0.0D0
	    ELSE
C	      beam array (see "24/8 07")
	      IF (IDOWNARR.GT.0) THEN
	        SLASK = DZ(NSO+1) / (NSARR-1)
	      ELSE
	        SLASK = DZ(NSO) / (NSARR-1)
	      END IF
	      CALL LIQPART ( OMEG,OMEG2,U,U2, UA2SOURC, SLASK, 
     $	       URA, EXPTYP, PLNREF, CP, STP,SDP )
              IF (PLNREF.EQ.0.0D0) THEN
                E1P = 1.0D0
	      ELSE
	        E1P = DEXP(-PLNREF)
              END IF
	      IF (EXPTYP) THEN
	        URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA)) ! -i*ura
	        YA = 0.5D0 * RW/URANEGI
	        YB = -YA
	      ELSE
	        YA = RW
	        YB = (0.0D0,0.0D0)
	      END IF
	      RWTLN = 0.0D0
	      E1PROD = 1.0D0
	      RWTOT = AMSOURCFW(1) * YA
	      RTTOT = AMSOURCFW(1) * YB
	      DO ID = 2, NSARR
	        RWTLN = RWTLN + PLNREF
                E1PROD = E1PROD*E1P
	        CSLASK = E1PROD * AMSOURCFW(ID)
	        IF (EXPTYP) THEN
C	          now the main propagation step, freq-dependent
	          IF (IDOWNARR.GT.0) THEN
C	            propagating upwards towards matching depth
	            RWTOT = CSLASK*YA + DCONJG(CP) * RWTOT
	            RTTOT = CSLASK*YB + DREAL(SDP)*CP * RTTOT
	          ELSE
C	            propagating downwards towards matching depth
	            RWTOT = CSLASK*YA + DREAL(SDP)*CP * RWTOT
	            RTTOT = CSLASK*YB + DCONJG(CP) * RTTOT
	          END IF
	        ELSE
C	          now the main propagation step, freq-dependent
	          IF (IDOWNARR.GT.0) THEN
C	            propagating upwards towards matching depth
	            WWLOK = CSLASK*YA + CP*RWTOT - STP*RTTOT
	            RTTOT = CSLASK*YB + SDP*RWTOT + CP*RTTOT
	            RWTOT = WWLOK
	          ELSE
C	            propagating downwards towards matching depth
	            WWLOK = CSLASK*YA + CP*RWTOT + STP*RTTOT
	            RTTOT = CSLASK*YB - SDP*RWTOT + CP*RTTOT
	            RWTOT = WWLOK
	          END IF
	        END IF
	      END DO
C	      now the last propagation step
	      IF (EXPTYP) THEN
	        WWLOK = RWTOT - RTTOT
	        RTTOT = RHOSOURC * ( RWTOT + RTTOT )
	        RWTOT = URANEGI*WWLOK
	      ELSE
	        RTTOT = RHOSOURC * RTTOT
	      END IF
	      CSLASK1 = ( -YPU2*RWTOT + YPU1*RTTOT ) / CDETLOK
	      CSLASK2 = ( -YPL2*RWTOT + YPL1*RTTOT ) / CDETLOK
C	      note that rwtln remains to be included now
	    END IF
	  ELSE
	    IF (NSARR.NE.0) STOP '%%% liqsi, beam array not impl here'
	    CSLASK = -RHOCDET*DEXP(RLNDETLOK-DLNDULOK) * RW / DERIVDULOK
	    CSLASK = 2*PI * DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))  ! 2*pi*i * 
C	    for residue contribution (positive orientation)
	    WRITE(*,*) ' indicators: ',
     $	    CDABS(CDETLOK) / (DMAX1(CDABS(YPL1),CDABS(YPL2))*
     $	    DMAX1(CDABS(YPU1),CDABS(YPU2))) , 
     $	    DEXP(RLNDET-DLNDULOK) * CDABS(CDET/DERIVDULOK)
	    CSLASK1 = YPU2 * CSLASK
	    CSLASK2 = YPL2 * CSLASK
	    RWTLN = 0.0D0
	  END IF
	  YPL1 = CSLASK1 * YPL1
	  YPL2 = CSLASK1 * YPL2
	  YPU1 = CSLASK2 * YPU1
	  YPU2 = CSLASK2 * YPU2
	ELSE
C	  mimic continuity for the vector now, in the following way
	  IF (NSO.GT.0) THEN
	    IF (RHO(NSO+1).NE.0.0D0) THEN
	      CSLASK1 = YPU2
	      CSLASK2 = YPL2
	      YPL1 = CSLASK1 * YPL1
	      YPL2 = CSLASK1 * YPL2
	      YPU1 = CSLASK2 * YPU1
	      YPU2 = CSLASK2 * YPU2
	    ELSE
	      CSLASK1 = (1.0D0,0.0D0)   ! now "nso.eq.nsw" and "free there"
	      CSLASK2 = (1.0D0,0.0D0)
	    END IF
	  ELSE
	    IF (RHO(0).NE.0.0D0) THEN
	      CSLASK1 = YPU2
	      CSLASK2 = YPL2
	      YPL1 = CSLASK1 * YPL1
	      YPL2 = CSLASK1 * YPL2
	      YPU1 = CSLASK2 * YPU1
	      YPU2 = CSLASK2 * YPU2
	    ELSE
	      CSLASK1 = (1.0D0,0.0D0)   ! now "nso.eq.0" and "free there"
	      CSLASK2 = (1.0D0,0.0D0)
	    END IF
	  END IF
	  RWTLN = 0.0D0
	END IF
	IF (NSO.LT.NSW) THEN
	  SLASK = DMAX1 ( DABS(DREAL(CSLASK1)), DABS(DIMAG(CSLASK1)) )
	  IF (SLASK.GT.0.0D0) THEN
	    CSLASK1 = CSLASK1/SLASK
	    IF (IRKSPEC.EQ.0) THEN
	      SLASK = DLOG(SLASK) + RWTLN
	    ELSE
	      SLASK = DLOG(SLASK) + RWTLN - DLNSPEC
	    END IF
	    DO ID = NSO+1, NSW-1
	      SLASK = SLASK - PLNDELTA(ID)
	      ISLASK = IHDD(ID)
	      IF (ISLASK.GT.0) THEN
	        CSLASK = DEXP(SLASK) * CSLASK1
ccc	        dhomw(islask) = cslask * dhomw(islask)
	        DHFAKT(ISLASK) = CSLASK * DHFAKT(ISLASK)
	      END IF
	    END DO
	    SLASK = SLASK - PLNDELTA(NSW)
	    CLIQL = DEXP(SLASK) * CSLASK1
	    ISLASK = IHDD(NSW)
	    IF (ISLASK.GT.0) THEN
ccc	      dhomw(islask) = cliql * dhomw(islask)
	      DHFAKT(ISLASK) = CLIQL * DHFAKT(ISLASK)
	    END IF
	  ELSE
	    DO ID = NSO+1, NSW
	      ISLASK = IHDD(ID)
	      IF (ISLASK.GT.0) THEN
ccc	        dhomw(islask) = (0.0d0,0.0d0)
	        DHFAKT(ISLASK) = (0.0D0,0.0D0)
	      END IF
	    END DO
	    CLIQL = (0.0D0,0.0D0)
	  END IF
	ELSE
	  IF ( RWTLN.EQ.0.0D0 .AND. IRKSPEC.EQ.0 ) THEN
	    CLIQL = CSLASK1
	  ELSE IF (IRKSPEC.EQ.0) THEN
	    CLIQL = DEXP(RWTLN) * CSLASK1
	  ELSE
	    CLIQL = DEXP(RWTLN-DLNSPEC) * CSLASK1
	  END IF
	END IF
	ISLASK = IHDD(NSO)
	IF (ISLASK.GT.0) THEN
	  IF ( RWTLN.EQ.0.0D0 .AND. IRKSPEC.EQ.0 ) THEN
ccc	    dhomw(islask) = ypl1
	    DHFAKT(ISLASK) = YPL2
	  ELSE IF (IRKSPEC.EQ.0) THEN
	    SLASK = DEXP(RWTLN)
ccc	    dhomw(islask) = slask * ypl1
	    DHFAKT(ISLASK) = SLASK * YPL2
	  ELSE
	    SLASK = DEXP(RWTLN-DLNSPEC)
ccc	    dhomw(islask) = slask * ypl1
	    DHFAKT(ISLASK) = SLASK * YPL2
	  END IF
	END IF
	IF (NSO.GT.(NSS+1)) THEN
	  SLASK = DMAX1 ( DABS(DREAL(CSLASK2)), DABS(DIMAG(CSLASK2)) )
	  IF (SLASK.GT.0.0D0) THEN
	    CSLASK2 = CSLASK2/SLASK
	    IF (IRKSPEC.EQ.0) THEN
	      SLASK = DLOG(SLASK) + RWTLN
	    ELSE
	      SLASK = DLOG(SLASK) + RWTLN - DLNSPEC
	    END IF
	    DO ID = NSO-1, NSS+1, -1
	      SLASK = SLASK - PLNDELTA(ID+1)
	      ISLASK = IHDD(ID)
	      IF (ISLASK.GT.0) THEN
	        CSLASK = DEXP(SLASK) * CSLASK2
ccc	        dhomw(islask) = cslask * dhomw(islask)
	        DHFAKT(ISLASK) = CSLASK * DHFAKT(ISLASK)
	      END IF
	    END DO
	    SLASK = SLASK - PLNDELTA(NSS+1)
	    CLIQU = DEXP(SLASK) * CSLASK2
	    ISLASK = IHDD(NSS)
	    IF (ISLASK.GT.0) THEN
ccc	      dhomw(islask) = cliqu * dhomw(islask)
	      DHFAKT(ISLASK) = CLIQU * DHFAKT(ISLASK)
	    END IF
	  ELSE
	    DO ID = NSO-1, NSS, -1
	      ISLASK = IHDD(ID)
	      IF (ISLASK.GT.0) THEN
ccc	        dhomw(islask) = (0.0d0,0.0d0)
	        DHFAKT(ISLASK) = (0.0D0,0.0D0)
	      END IF
	    END DO
	    CLIQU = (0.0D0,0.0D0)
	  END IF
	ELSE IF (NSO.EQ.(NSS+1)) THEN
	  IF ( RWTLN.EQ.0.0D0 .AND. IRKSPEC.EQ.0 ) THEN
	    CLIQU = DEXP(-PLNDELTA(NSS+1)) * CSLASK2
	  ELSE IF (IRKSPEC.EQ.0) THEN
	    CLIQU = DEXP(-PLNDELTA(NSS+1)+RWTLN) * CSLASK2
	  ELSE
	    CLIQU = DEXP(-PLNDELTA(NSS+1)+RWTLN-DLNSPEC) * CSLASK2
	  END IF
	  ISLASK = IHDD(NSS)
	  IF (ISLASK.GT.0) THEN
ccc	    dhomw(islask) = cliqu * dhomw(islask)
	    DHFAKT(ISLASK) = CLIQU * DHFAKT(ISLASK)
	  END IF
	ELSE
	  IF ( RWTLN.EQ.0.0D0 .AND. IRKSPEC.EQ.0 ) THEN
	    CLIQU = CSLASK2
	  ELSE IF (IRKSPEC.EQ.0) THEN
	    CLIQU = DEXP(RWTLN) * CSLASK2
	  ELSE
	    CLIQU = DEXP(RWTLN-DLNSPEC) * CSLASK2
	  END IF
	END IF

C	the dhfakt array refers to the "fourth component" of the "Aki (!) 
C	displacement-stress vector (Aki & Richards 1980, 7.4.2) but with 
C	omeg* for the first two components" with certain integration factors 
C	included (see comments to dhcomp)
C	dhfakt(1:nhyd) concern the total part "directly" (cf. Aki 
C	fig. 7.9)
cccC	the dhomw array is analogous but for the "second component" (so 
cccC	"omeg*" is really involved here)

	RETURN
	END
