C	subroutine harksi for rpress.f etc

C	Sven Ivansson  1990





	SUBROUTINE HARKSI ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	IHARK0, IHARKLIQ, RHOCDET, YT1,YT2,YT3,YT4,YT6, TLN, YP1,YP2 )
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 minor calculations, through solid layers to a 
C	solid/fluid interface, Cramer's rule is used to solve 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 
C	first two components" so that the shear stress is zero when the 
C	fluid is entered
C	rhocdet is input/output
C	the start values (with de-scaling as given by start value for 
C	tln) for the row-vector  (t1,t2,t3,t4,t5,t6) = 
C	(yt1,yt2,yt3,yt4,-yt2,yt6)  with t5=-t2   (OBS)   
C	are given as input by "yt1,yt2,yt3,yt4,yt6"  but different cases 
C	exist as indicated by input parameter ihark0
C	note that "yt1,yt2,yt3,yt4,yt6" is then used as work-space
C	when iharkliq is nonzero on entry, however, output "yt1,yt2,yt3,yt4,
C	yt6" with de-scaling as given by "output tln" is provided as well
C	tln is de-scaling for "yp1,yp2" (input/output)
C	"yp1,yp2" is output ( yp1:=-yt3 yp2:=yt1 )
C	implementation of Harkrider (1970) p. 1940-1942 & Appendix 1 
C	and Watson (1970) Appendix B appropriately modified
C	note: - for his notation see also Harkrider (1964) p. 651
C	      - conjugation in the end to accomodate "exp(-i*omeg*t)" instead 
C	        of "exp(i*omeg*t)" "trivial" for the Aki (!) vector, note 
C	        the start values however
C	      - scaling of the "start values" is done
C	      - scaling with "e2t" for minors (from harksipart)
C	      - t1,t2,t3,t4,t5,t6 --> t1,t2,t3,t4,-t2,t6 since the second and 
C	        fifth columns of a 6x6 minor matrix B are equal but of 
C	        opposite sign (except that b25=-b22+1 and b55=-b52+1) 
C	        and since the corresponding half-space matrix subdeterminants 
C               are also equal but of opposite sign
C             - to incorporate the "t5" contributions to later matrix 
C	        products the values of "t2" are temporarily doubled, this 
C	        will be ok since the second and fifth rows of a 6x6 minor 
C	        matrix are equal but of opposite sign (except that 
C	        b52=-b22+1 and b55=-b25+1), the average of "b22 and -b52" 
C	        is used as "b22"
C	      - my own (Sven Ivansson) improvements of "fast version of 
C	        Knopoff's method" (Schwab in BSSA 60:5 1970) implemented, 
C	        "propagation steps" corresponding to a decomposition of B as 
C	        a product of simple (sparse) matrices where the frequency 
C	        dependence is isolated in one matrix

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       nol <= noldim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'DPCFSUB.INC'
	INCLUDE 'FSUBCOM.INC'

	INTEGER NL1,NL2,NRIKT, IHARK0, IHARKLIQ
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2

	REAL*8 RHOCDET, TLN
	COMPLEX*16 YT1,YT2,YT3,YT4,YT6

	COMPLEX*16 YP1,YP2

	INTEGER ID, IND
	REAL*8 E2T, RH1,RH2,RH1OLD,RH2OLD,RHQ, TMAX, ZENDFSUB, Z0FSUB, TLNREF
	COMPLEX*16 G0, G0OLD, G0D, G0DU2PM, G0U, 
     $	URA,URB, URAURB,URANEGI,URBNEGI, 
     $	CP,CQ, STP,STQ,SDP,SDQ, 
     $	T1Q,T2Q,T3Q,T4Q,T6Q, CSLASK, 
     $	YT(5)
	LOGICAL ALLCOMP, JOINLF, DOFIRST, LUPDATE, EXPTYP

	EXTERNAL FSUBHARK5

C	minor propagation, the "t row-vector with minors from two rows" is 
C	at the leftmost end, matrix product "from left to right" with the 
C	propagation matrices 
C	"  upwards when  nrikt = -1  ,  downwards when  nrikt = +1  "

	IF (IHARK0.EQ.0) THEN
C	  "ordinary" start vector
	  JOINLF = .FALSE.
	  DOFIRST = .TRUE.
	ELSE
C	  required that "do id = nl1, nl2, nrikt" will not be the empty 
C	  loop but this is not tested here;  one exception: the "empty loop" 
C	  case below when "ihark0=1"
	  IF (IHARK0.EQ.1) THEN
C	    "unfinished" start vector from homogeneous solid half-space (when 
C	    "id=nl1-nrikt")
	    ID = NL1-NRIKT
	    RH1 = RHO(ID)                          ! local rho(id)
	    RH2 = 1/RH1                            ! could be stored
	    G0 = TWOVB2(ID)                        ! local twovb2(id)
	    IF (ID.EQ.NL2) THEN
C	      "empty loop" case: as "last section" in this "if-block"
	      G0U = G0*U
	      YT6 = YT6 - U2
	      YT2 = -YT6*G0U - U
	      YT1 = 1 + (YT2-U)*G0U
	      YT1 = RH1*YT1
	      YT6 = RH2*YT6
	      JOINLF = .FALSE.
	      DOFIRST = .TRUE.
	    ELSE IF (IFSMTYP(NL1).LE.-4) THEN
	      RH1OLD = RH1
	      RH2OLD = RH2
	      G0OLD = G0
	      RH1 = RHO(NL1)                       ! local rho(nl1)
	      RH2 = 1/RH1                          ! could be stored
	      RHQ = RH1*RH2OLD      ! "extra"      ! could be stored
	      G0 = TWOVB2(NL1)                     ! local twovb2(nl1)
	      G0D = RHQ*G0 - G0OLD  ! "extra"      ! could be stored
C	      now follows "extra" for this case
	      G0DU2PM = (1-RHQ) + G0D*U2           ! "1-rhq" could be stored
C	      now the "rho*Cdelta frm homogeneous solid half-space & Cinvdelta 
C	      from id=nl1 layer" propagation steps jointly recognizing that 
C	      "start yt1" is one and "start yt2" is zero in this case
	      YT6 = RH1OLD*YT6
	      T6Q = YT6 - RH1OLD*U2
	      T1Q = RH1OLD - T6Q*G0D
	      CSLASK = YT6 + T1Q*U2
	      YT1 = T1Q + CSLASK*G0D
	      YT2 = (YT1-RHQ*T1Q)*U
	      YT6 = YT6 - (CSLASK+RHQ*T6Q)*G0DU2PM
	      YT3 = RH1*YT3
	      YT4 = RH1*YT4
	      JOINLF = .FALSE.    ! already done
	      DOFIRST = .FALSE.   ! already done
	      LUPDATE = .TRUE.
	    ELSE
C	      now follows "extra" for this case
	      G0U = G0*U
C	      now the "rho*Cdelta propagation step from the homogeneous solid 
C	      half-space" recognizing that "start yt1" is one and "start yt2" 
C	      is zero in this case
	      YT6 = YT6 - U2
	      YT2 = -YT6*G0U - U
	      YT1 = 1 + (YT2-U)*G0U
	      YT1 = RH1*YT1
	      YT6 = RH2*YT6
C	      in the "empty loop" case with "iharkliq=0", we would only 
C	      have needed to update yt1 (since yt3 is unchanged)
	      JOINLF = .FALSE.
	      DOFIRST = .TRUE.
	    END IF
	  ELSE IF (IHARK0.EQ.-1) THEN
C	    "ordinary" start vector from "rigid"
	    IF (IFSMTYP(NL1).LE.-4) THEN
	      RH1 = RHO(NL1)                       ! local rho(nl1)
	      RH2 = 1/RH1                          ! could be stored
	      G0 = TWOVB2(NL1)                     ! local twovb2(nl1)
C	      now the "Cinvdelta propagation step from id=nl1 layer" 
C	      recognizing that "start yt1" is one and "start yt2,yt3,yt4,yt6" 
C	      are zero in this case
	      YT2 = U
	      YT6 = -U2
C	      actually the next propagation step could also be 
C	      simplified (yt1 is still one and yt3,yt4 are still zero)
	      JOINLF = .FALSE.
	      DOFIRST = .FALSE.   ! already done
	      LUPDATE = .TRUE.    ! "g0u" has not been updated
	    ELSE
	      JOINLF = .FALSE.
	      DOFIRST = .TRUE.
	    END IF
	  ELSE IF (IHARK0.EQ.-2) THEN
C	    "ordinary" start vector from "pressure-release"
	    IF (IFSMTYP(NL1).LE.-4) THEN
	      RH1 = RHO(NL1)                       ! local rho(nl1)
	      RH2 = 1/RH1                          ! could be stored
	      G0 = TWOVB2(NL1)                     ! local twovb2(nl1)
C	      now follows "extra" for this case
	      G0U = G0*U
C	      now the "Cinvdelta propagation step from id=nl1 layer" recogniz-
C	      ing that "start yt1,yt2,yt3,yt4" are zero and "start yt6" is one 
C	      in this case, neglecting "final rh1*rh1 vector-multiplication" 
C	      which could be done simply by adding to tln a stored 
C	      "dlog(rh1*rh1)" <OBS>
	      RHOCDET = RHOCDET * (RH1*RH1)
	      YT1 = -G0U*G0U
	      T1Q = YT1*U + G0U
	      YT6 = YT6 - (T1Q+G0U)*U
	      YT2 = T1Q
C	      actually the next propagation step could also be 
C	      simplified (yt3,yt4 are still zero)
	      JOINLF = .FALSE.
	      DOFIRST = .FALSE.   ! already done
	      LUPDATE = .FALSE.
	    ELSE
	      JOINLF = .FALSE.
	      DOFIRST = .TRUE.
	    END IF
	  ELSE
	    STOP 4516
	  END IF
	END IF

	DO ID = NL1, NL2, NRIKT

	  IF (ID.NE.NL2) THEN
	    ALLCOMP = .TRUE.
	  ELSE
	    IF (IHARKLIQ.EQ.0) THEN
	      ALLCOMP = .FALSE.
	    ELSE
	      ALLCOMP = .TRUE.
	    END IF
	  END IF

	  IF (IFSMTYP(ID).LE.-4) THEN

C	    with "similar layers in sequence" (which is easily avoided) one 
C	    would need to update certain factors for the matrix elements "only 
C	    the first time", "ifsmtyp" could be used in terms of the condition 
C	    "  if ( ifsmtyp(id).ne.-6 .and. ifsmtyp(id).ne.(nrikt-6) ) then  " 
C	    but the use of "joinlf" makes things complicated (possibly a 
C	    "chain" to follow)
C	    warning if implemented: not allowed to skip for "first id" in 
C	    non-empty loop
	    IF (JOINLF) THEN
	      RH1 = RHO(ID)                        ! local rho(id)
	      RH2 = 1/RH1                          ! could be stored
	      RHQ = RH1*RH2OLD      ! "extra"      ! could be stored
	      G0 = TWOVB2(ID)                      ! local twovb2(id)
	      G0D = RHQ*G0 - G0OLD  ! "extra"      ! could be stored
C	      now follows "extra" for this case
	      G0DU2PM = (1-RHQ) + G0D*U2           ! "1-rhq" could be stored
C	      now the "Cdelta from old layer & Cinvdelta from current layer" 
C	      propagation steps jointly
	      T2Q = RHQ*YT2
	      YT2 = YT2+YT2         ! temporary doubling
	      T6Q = YT6 + (YT2-YT1*U)*U
	      T1Q = YT1 - T6Q*G0D
	      CSLASK = YT6 + T1Q*U2
	      YT1 = T1Q + CSLASK*G0D
	      YT2 = T2Q + (YT1-RHQ*T1Q)*U
	      YT6 = YT6 - (CSLASK+RHQ*T6Q)*G0DU2PM
	      YT3 = RHQ*YT3
	      YT4 = RHQ*YT4
	      LUPDATE = .TRUE.
	    ELSE IF (DOFIRST) THEN
	      RH1 = RHO(ID)                        ! local rho(id)
	      RH2 = 1/RH1                          ! could be stored
	      G0 = TWOVB2(ID)                      ! local twovb2(id)
C	      now follows "extra" for this case
	      G0U = G0*U
C	      now the Cinvdelta propagation step, with "special" first, 
C	      neglecting "final rh1 vector-multiplication" which could be 
C	      done simply by adding to tln a stored "dlog(rh1)" <OBS>
	      RHOCDET = RHOCDET * RH1
	      YT1 = RH2*YT1  ! "special"
	      YT6 = RH1*YT6  ! "special"
	      T2Q = YT2
	      YT2 = YT2+YT2         ! temporary doubling
	      YT1 = YT1 - (YT2+YT6*G0U)*G0U
	      T6Q = YT6*G0U
	      T1Q = YT1*U + T6Q
	      YT6 = YT6 - (YT2+T1Q+T6Q)*U
	      YT2 = T2Q + T1Q
	      LUPDATE = .FALSE.
	    ELSE
C	      already done the "Cinvdelta propagation step" here
C	      (lupdate has already been set in this case)
	      DOFIRST = .TRUE.   ! will not be changed any more
	    END IF
	    CALL HARKSIPART ( OMEG,OMEG2,U,U2, 
     $	    UA2(ID),UB2(ID), DZ(ID), 
     $	    URA,URB, EXPTYP, 
     $	    TLNREF,E2T, CP,CQ, STP,STQ,SDP,SDQ )               ! freq-dependent
	    IF (EXPTYP) THEN
C	      freq-independent factors
	      URAURB = URA*URB
	      URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	      URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb
C	      now the Dinvdeltamod propagation step
	      YT1 = -YT1*URAURB
	      YT3 = YT3*URBNEGI
	      YT4 = YT4*URANEGI
	      T1Q = YT1 + YT3
	      T2Q = YT1 - YT3
	      T3Q = YT4 + YT6
	      T4Q = YT4 - YT6
	      YT1 = T1Q - T4Q
	      YT3 = -T2Q + T3Q
	      YT4 = -T1Q - T4Q
	      YT6 = -T2Q - T3Q
C	      now the PEdelta propagation step, freq-dependent
	      TLN = TLN + TLNREF
	      YT2 = E2T*YT2
	      IF (NRIKT.LT.0) THEN
	        YT1 = YT1*STP
	        YT3 = YT3*SDP
	        YT4 = YT4*STQ
	        YT6 = YT6*SDQ
	      ELSE
	        YT1 = YT1*SDQ
	        YT3 = YT3*STQ
	        YT4 = YT4*SDP
	        YT6 = YT6*STP
	      END IF    ! note that "row vector comes from the left"
C	      now the Ddeltamod propagation step
	      T1Q = YT1 + YT3
	      T2Q = YT1 - YT3
	      T3Q = YT4 + YT6
	      T4Q = YT4 - YT6
	      YT1 = T2Q - T3Q
	      YT3 = T1Q - T4Q
	      YT4 = -T2Q - T3Q             ! only needed if "allcomp"
	      YT6 = T1Q + T4Q
	      TLN = TLN - DLOG4
	      YT1 = -YT1/URAURB
	      YT2 = 4*YT2
	      YT3 = YT3/URBNEGI
	      YT4 = YT4/URANEGI            ! only needed if "allcomp"
	    ELSE
C	      now the PTdelta propagation step, freq-dependent
	      TLN = TLN + TLNREF
	      YT2 = E2T*YT2   ! decoupled here
	      IF (NRIKT.LT.0) THEN
	        T1Q =  YT1*CQ  + YT4*SDQ
	        T3Q =  YT3*CQ  - YT6*SDQ
	        T4Q = -YT1*STQ + YT4*CQ
	        T6Q =  YT3*STQ + YT6*CQ
	        YT1 =  T1Q*CP  - T3Q*SDP
	        YT3 =  T1Q*STP + T3Q*CP
	        YT6 = -T4Q*STP + T6Q*CP
	        IF (ALLCOMP) THEN
	          YT4 = T4Q*CP + T6Q*SDP      ! otherwise not needed
	        END IF
	      ELSE
	        T1Q =  YT1*CQ  - YT4*SDQ
	        T3Q =  YT3*CQ  + YT6*SDQ
	        T4Q =  YT1*STQ + YT4*CQ
	        T6Q = -YT3*STQ + YT6*CQ
	        YT1 =  T1Q*CP  + T3Q*SDP
	        YT3 = -T1Q*STP + T3Q*CP
	        YT6 =  T4Q*STP + T6Q*CP
	        IF (ALLCOMP) THEN
	          YT4 = T4Q*CP - T6Q*SDP      ! otherwise not needed
	        END IF
	      END IF    ! note that "row vector comes from the left"
	    END IF
	    IF (ID.EQ.NL2) THEN
	      JOINLF = .FALSE.
	    ELSE
	      IF (IFSMTYP(ID+NRIKT).LE.-4) THEN
	        JOINLF = .TRUE.
	      ELSE
	        JOINLF = .FALSE.
	      END IF
	    END IF
	    IF (JOINLF) THEN
C	      postpone the Cdelta propagation step ("next layer, jointly")
	      RH2OLD = RH2
	      G0OLD = G0
	    ELSE
	      IF (LUPDATE) THEN
C	        now follows "extra" for this case (not updated previously)
	        G0U = G0*U
	      END IF
C	      now the Cdelta propagation step
	      T2Q = YT2                    ! only needed if "allcomp"
	      YT2 = YT2+YT2           ! temporary doubling
	      YT6 = YT6 + (YT2-YT1*U)*U
	      T1Q = YT1*U
	      T6Q = YT6*G0U + T1Q
	      YT1 = YT1 + (YT2-T1Q-T6Q)*G0U
	      IF (ALLCOMP) THEN
C	        , with "special" last, neglecting "final rh2 vector-
C	        multiplication" which could be done simply by adding 
C	        to tln a stored "dlog(rh2)" <OBS>
	        RHOCDET = RHOCDET * RH2
	        YT2 = T2Q - T6Q
	        YT1 = RH1*YT1  ! "special"
	        YT6 = RH2*YT6  ! "special"
	      ELSE
C	        ( we only need to update yt1 and yt3 in this case )
	        YT3 = RH2*YT3
	      END IF
	    END IF

	  ELSE

	    YT(1) = YT1
	    YT(2) = YT2
	    YT(3) = YT3
	    YT(4) = YT4
	    YT(5) = YT6                   ! obs
	    IFSMTYPLOK = IFSMTYP(ID)
	    IF (NRIKT.LT.0) THEN
	      ZENDFSUB = ZD(ID-1)
	      Z0FSUB = ZD(ID)
	    ELSE
	      ZENDFSUB = ZD(ID)
	      Z0FSUB = ZD(ID-1)
	    END IF
	    CALL DPCLIN ( FSUBHARK5,5, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	    HFAKTMAX,HOFF,EPSMAGNIFMAX, ZENDFSUB, Z0FSUB,TLN,YT,
     $	    EPSFSUB,HSTEPFSUB, IND )
	    IF (IND.NE.0) STOP 6760
	    YT1 = YT(1)
	    YT2 = YT(2)
	    YT3 = YT(3)
	    YT4 = YT(4)
	    YT6 = YT(5)                   ! obs

	  END IF

	  IF (ALLCOMP) THEN
	    TMAX = DMAX1 ( DABS(DREAL(YT1)), DABS(DIMAG(YT1)), 
     $	    DABS(DREAL(YT2)), DABS(DIMAG(YT2)), 
     $	    DABS(DREAL(YT3)), DABS(DIMAG(YT3)), 
     $	    DABS(DREAL(YT4)), DABS(DIMAG(YT4)), 
     $	    DABS(DREAL(YT6)), DABS(DIMAG(YT6)) )
	    IF ( TMAX.LT.tspmaxlim1 .OR. tspmaxlim2.LT.TMAX ) THEN
	      TMAX = 1/TMAX   ! if tmax was 0.0d0, due to num. canc.
	      TLN = TLN - DLOG(TMAX)
	      YT1 = TMAX*YT1
	      YT2 = TMAX*YT2
	      YT3 = TMAX*YT3
	      YT4 = TMAX*YT4
	      YT6 = TMAX*YT6
	    END IF
	  ELSE
	    TMAX = DMAX1 ( DABS(DREAL(YT1)), DABS(DIMAG(YT1)), 
     $	    DABS(DREAL(YT3)), DABS(DIMAG(YT3)) )
	    IF (TMAX.GT.0.0D0) THEN
	      TMAX = 1/TMAX
	      TLN = TLN - DLOG(TMAX)
	      YT1 = TMAX*YT1
	      YT3 = TMAX*YT3
	    END IF
	  END IF

	END DO

	YP1 = -YT3
	YP2 = YT1

C	(ww,pp) = (yp1,yp2)
C	the "Aki (!) displacement-stress vector (Aki & Richards 1980, 7.4.2) 
C	but with omeg* for the first two components" at the 
C	interface is  (...,ww,0.0,pp)  or , if iharkliq was nonzero on entry,  
C	( t5=-t2, -t3, 0.0, t1 ) = ( -yt2, -yt3, 0.0, yt1 )  on a solid side 
C	of the interface
C	(note that the absolute level is irrelevant, scaling has been done)
C	(the de-scaling exponent is kept track of in input/output "tln" to 
C	make it possible to extract an analytic function of u)
C	(if (yp1,yp2) comes out as the zero vector, a normal mode has been 
C	found that is "quiet in the fluid")

	RETURN
	END