
C	additional subroutines (and functions) for "Full Wave" computations 
C	special "fw" versions :  fwsubytzsol  fwsubytzliq

C	Sven Ivansson  1994






	SUBROUTINE FWSUBYTZSOL ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, DIS1FW,DIS2FW,DIS3FW,DIS4FW, 
     $	TLN, YT1,YT2,YT3,YT4,YT6, Z1,Z2,Z3,Z4, 
     $	TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, Z1FW,Z2FW,Z3FW,Z4FW )
C	special version for "full wave"
C	in particular: "parameter change"
C	               "always ihark0=0" (so ihark0 not needed)
C	               "no joinlf"
C	               "always allcomp" (so iharkliq & allcomp not needed)
C	               "always dofirst" (so dofirst not needed)
C	               "lupdate" not needed
C	additional "fw input and output" in "common /fwparam,fwpart/"
C	note: -> propagation of "row-vectors from the left" 
C	         in particular: "downward matrices when nrikt=-1"
C	                        "  upward matrices when nrikt=+1"
C	         (however, dpclin & fsubhark5 & fsubzsol4)
C	      -> for yt; homogeneous case
C	         for z; inhomogeneous case

	USE INMODELa
	USE FWPARTa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C       nol <= noldim

	INCLUDE '../DPCFSUB.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE 'FWPARAM.INC'

	INTEGER NL1,NL2,NRIKT
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2, 
     $	DIS1FW(*), DIS2FW(*), DIS3FW(*), DIS4FW(*)
	LOGICAL INIT

	REAL*8 TLN                                      ! input/output
	COMPLEX*16 YT1,YT2,YT3,YT4,YT6, Z1,Z2,Z3,Z4     ! input/output

	REAL*8 TLNFW(*)
	COMPLEX*16 YT1FW(*), YT2FW(*), YT3FW(*), YT4FW(*), 
     $	YT5FW(*), YT6FW(*), 
     $	Z1FW(*), Z2FW(*), Z3FW(*), Z4FW(*)

	INTEGER INDOPT         ! input (0, 1,2,12)
	INTEGER NSOOPT         ! input/output
	REAL*8 QIOPT           ! input/output
	REAL*8 TLNOPT(*)       ! output if indopt.eq.0, else input
	COMPLEX*16 YT1OPT(*), YT2OPT(*), YT3OPT(*), YT4OPT(*), 
     $	YT6OPT(*)              ! output if indopt.eq.0, else input

	INTEGER ID, IDQ, IND
	REAL*8 TLNREF,E2T, ZLN, ZLNREF,E1X, RH1,RH2, 
     $	ZENDFSUB, Z0FSUB, RLNDETLOK, TLNOPTLOK, SLASK
	COMPLEX*16 G0, G0U, 
     $	URA,URB, URAURB,URANEGI,URBNEGI, 
     $	CP,CQ, STP,STQ,SDP,SDQ, 
     $	T1Q,T2Q,T3Q,T4Q,T6Q, S1Q,S2Q, 
     $	CDETLOK, YT1OPTLOK,YT2OPTLOK,YT3OPTLOK,YT4OPTLOK,YT6OPTLOK, 
     $	YT(5), Z(4)
	LOGICAL EXPTYP, SETZLOK, OPTLOK

	EXTERNAL FSUBHARK5, FSUBZSOL4

	SETZLOK = .TRUE.
	OPTLOK = .FALSE.
	GOTO 10

	ENTRY FWSUBYTSOL ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, 
     $	TLN, YT1,YT2,YT3,YT4,YT6, 
     $	TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW )

	SETZLOK = .FALSE.
	OPTLOK = .FALSE.
	GOTO 10

	ENTRY FWOPTYTZSOL ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, DIS1FW,DIS2FW,DIS3FW,DIS4FW, 
     $	TLN, YT1,YT2,YT3,YT4,YT6, Z1,Z2,Z3,Z4, 
     $	TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, Z1FW,Z2FW,Z3FW,Z4FW, 
     $	INDOPT, TLNOPT,YT1OPT,YT2OPT,YT3OPT,YT4OPT,YT6OPT, NSOOPT,QIOPT )
C	note: for indopt.ne.0, same act arg for "ytopt","ytfw" ok (overwrite)

	SETZLOK = .TRUE.
	GOTO 20

	ENTRY FWOPTYTSOL ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, 
     $	TLN, YT1,YT2,YT3,YT4,YT6, 
     $	TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, 
     $	INDOPT, TLNOPT,YT1OPT,YT2OPT,YT3OPT,YT4OPT,YT6OPT, NSOOPT,QIOPT )
C	note: for indopt.ne.0, same act arg for "ytopt","ytfw" ok (overwrite)

	SETZLOK = .FALSE.
 20	OPTLOK = .TRUE.

	IF ( INDOPT.EQ.1 .OR. INDOPT.EQ.12 ) THEN
	  IDQ = NL1 - NRIKT
	  TLNOPTLOK = TLNOPT(IDQ)
	  YT1OPTLOK = YT1OPT(IDQ)
	  YT2OPTLOK = YT2OPT(IDQ)
	  YT3OPTLOK = YT3OPT(IDQ)
	  YT4OPTLOK = YT4OPT(IDQ)
	  YT6OPTLOK = YT6OPT(IDQ)
C	  cdetlok sign will fit
	  RLNDETLOK = TLN + TLNOPTLOK
	  CDETLOK = YT1*YT6OPTLOK + 2*YT2*YT2OPTLOK + 
     $	  YT3*YT4OPTLOK + YT4*YT3OPTLOK + YT6*YT1OPTLOK     ! yt5=-yt2
	  SLASK = DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	  ( 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)) ) * 
     $	  DMAX1 ( DABS(DREAL(YT1OPTLOK)),DABS(DIMAG(YT1OPTLOK)),
     $	  DABS(DREAL(YT2OPTLOK)),DABS(DIMAG(YT2OPTLOK)),
     $	  DABS(DREAL(YT3OPTLOK)),DABS(DIMAG(YT3OPTLOK)),
     $	  DABS(DREAL(YT4OPTLOK)),DABS(DIMAG(YT4OPTLOK)),
     $	  DABS(DREAL(YT6OPTLOK)),DABS(DIMAG(YT6OPTLOK)) ) )
	  IF (SLASK.LT.QIOPT) THEN
	    IF (NRIKT.GT.0) THEN
	      NSOOPT = NL1 - 1
	    ELSE
	      NSOOPT = NL1
	    END IF
	    QIOPT = SLASK
	  END IF
	END IF

 10	CONTINUE

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  "

	DO ID = NL1, NL2, NRIKT

	  IF (OPTLOK) THEN
	    IF (INDOPT.EQ.0) THEN
	      TLNOPT(ID) = TLN
	      YT1OPT(ID) = YT1
	      YT2OPT(ID) = YT2
	      YT3OPT(ID) = YT3
	      YT4OPT(ID) = YT4
	      YT6OPT(ID) = YT6
	    ELSE
	      IF (ID.NE.NL1) THEN
C	        cdetlok sign will fit
	        RLNDETLOK = TLN + TLNOPTLOK
	        CDETLOK = YT1*YT6OPTLOK + 2*YT2*YT2OPTLOK + 
     $	        YT3*YT4OPTLOK + YT4*YT3OPTLOK + YT6*YT1OPTLOK     ! yt5=-yt2
	        SLASK = 
     $	        DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	        ( 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)) ) * 
     $	        DMAX1 ( DABS(DREAL(YT1OPTLOK)),DABS(DIMAG(YT1OPTLOK)),
     $	        DABS(DREAL(YT2OPTLOK)),DABS(DIMAG(YT2OPTLOK)),
     $	        DABS(DREAL(YT3OPTLOK)),DABS(DIMAG(YT3OPTLOK)),
     $	        DABS(DREAL(YT4OPTLOK)),DABS(DIMAG(YT4OPTLOK)),
     $	        DABS(DREAL(YT6OPTLOK)),DABS(DIMAG(YT6OPTLOK)) ) )
	        IF (SLASK.LT.QIOPT) THEN
	          IF (NRIKT.GT.0) THEN
	            NSOOPT = ID - 1
	          ELSE
	            NSOOPT = ID
	          END IF
	          QIOPT = SLASK
	        END IF
	      END IF
	      TLNOPTLOK = TLNOPT(ID)
	      YT1OPTLOK = YT1OPT(ID)
	      YT2OPTLOK = YT2OPT(ID)
	      YT3OPTLOK = YT3OPT(ID)
	      YT4OPTLOK = YT4OPT(ID)
	      YT6OPTLOK = YT6OPT(ID)
	    END IF
	  END IF

C	  preliminary setting that can be overwritten below
	  TLNFW(ID) = TLN
	  YT1FW(ID) = YT1
	  YT2FW(ID) = YT2
	  YT3FW(ID) = YT3
	  YT4FW(ID) = YT4
	  YT6FW(ID) = YT6

	  IF (SETZLOK) THEN
	    ZLN = TLN
	    Z1FW(ID) = Z1  ! preliminary setting that can be overwritten below
	    Z2FW(ID) = Z2  !   ...
	    Z3FW(ID) = Z3  !   ...
	    Z4FW(ID) = Z4  !   ...
	    IF (NRIKT.LT.0) THEN
	      Z1 = Z1 - YT1*DIS2FW(ID) - 
     $	      YT2*DIS3FW(ID) - YT3*DIS4FW(ID)
	      Z2 = Z2 + YT1*DIS1FW(ID) - 
     $	      YT4*DIS3FW(ID) + YT2*DIS4FW(ID)   ! yt5=-yt2
	      Z3 = Z3 + YT2*DIS1FW(ID) + 
     $	      YT4*DIS2FW(ID) - YT6*DIS4FW(ID)
	      Z4 = Z4 + YT3*DIS1FW(ID) - 
     $	      YT2*DIS2FW(ID) + YT6*DIS3FW(ID)   ! yt5=-yt2
	    END IF
	  END IF

C	  now yt propagation

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

	    IF ( IFSMTYP(ID).NE.-6 .AND. IFSMTYP(ID).NE.(NRIKT-6) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
	      RH1 = RHO(ID)                        ! local rho(id)
	      RH2 = 1/RH1                          ! could be stored
	      G0 = TWOVB2(ID)                      ! local twovb2(id)
	      G0U = G0*U
	      IF (INIT) THEN
	        CALL HSOLPART ( OMEG,OMEG2,U,U2, UA2(ID),UB2(ID), DZ(ID), 
     $	        URA,URB, EXPTYP, 
     $	        ZLNREF,E1X, CP,CQ, STP,STQ,SDP,SDQ, 
     $	        TLNREF,E2T )                                ! freq-dependent
	        URA_FW(ID) = URA
	        URB_FW(ID) = URB
	        EXPTYP_FW(ID) = EXPTYP
	        SLNREF_FW(ID) = ZLNREF
	        E1X_FW(ID) = E1X
	        CP_FW(ID) = CP
	        CQ_FW(ID) = CQ
	        STP_FW(ID) = STP
	        STQ_FW(ID) = STQ
	        SDP_FW(ID) = SDP
	        SDQ_FW(ID) = SDQ
	        TLNREF_FW(ID) = TLNREF
	        E2T_FW(ID) = E2T
	      ELSE
	        URA = URA_FW(ID)
	        URB = URB_FW(ID)
	        EXPTYP = EXPTYP_FW(ID)
	        ZLNREF = SLNREF_FW(ID)
	        E1X = E1X_FW(ID)
	        CP = CP_FW(ID)
	        CQ = CQ_FW(ID)
	        STP = STP_FW(ID)
	        STQ = STQ_FW(ID)
	        SDP = SDP_FW(ID)
	        SDQ = SDQ_FW(ID)
	        TLNREF = TLNREF_FW(ID)
	        E2T = E2T_FW(ID)
	      END IF
	      IF (EXPTYP) THEN
	        URAURB = URA*URB
	        URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	        URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb
	      END IF
	    ELSE IF (INIT) THEN
	      URA_FW(ID) = URA_FW(ID-NRIKT)
	      URB_FW(ID) = URB_FW(ID-NRIKT)
	      EXPTYP_FW(ID) = EXPTYP_FW(ID-NRIKT)
	      SLNREF_FW(ID) = SLNREF_FW(ID-NRIKT)
	      E1X_FW(ID) = E1X_FW(ID-NRIKT)
	      CP_FW(ID) = CP_FW(ID-NRIKT)
	      CQ_FW(ID) = CQ_FW(ID-NRIKT)
	      STP_FW(ID) = STP_FW(ID-NRIKT)
	      STQ_FW(ID) = STQ_FW(ID-NRIKT)
	      SDP_FW(ID) = SDP_FW(ID-NRIKT)
	      SDQ_FW(ID) = SDQ_FW(ID-NRIKT)
	      TLNREF_FW(ID) = TLNREF_FW(ID-NRIKT)
	      E2T_FW(ID) = E2T_FW(ID-NRIKT)
	    END IF
C	    now the Cinvdelta propagation step, with "special" first, 
C	    neglecting "final rh1 vector-multiplication" which could be done 
C	    simply by adding to tln a stored "dlog(rh1)" <OBS>
	    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
	    IF (EXPTYP) THEN
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
	      IF (IYSTAB.GT.0) THEN
C	        set as if Dinvdelta above
C	        note: "special" above, must be matched for z
	        TLNFW(ID) = TLN
	        YT1FW(ID) = YT1
	        YT2FW(ID) = -2*YT2*URBNEGI
	        YT3FW(ID) = YT3
	        YT4FW(ID) = YT4
	        YT5FW(ID) = -2*YT2*URANEGI   ! before, yt5=-yt2
	        YT6FW(ID) = YT6
	      END IF
C	      now the PEdelta propagation step, freq-dependent
	      TLN = TLN + TLNREF
	      YT2 = E2T*YT2
	      IF (NRIKT.LT.0) THEN
	        YT1 = YT1 * (E2T*E2T)*STP
	        YT3 = YT3 * DREAL(SDQ)*DCONJG(STQ)
	        YT4 = YT4 * DREAL(SDP)*STQ
	        YT6 = YT6 * DCONJG(STP)
	      ELSE
	        YT1 = YT1 * DCONJG(STP)
	        YT3 = YT3 * DREAL(SDP)*STQ
	        YT4 = YT4 * DREAL(SDQ)*DCONJG(STQ)
	        YT6 = YT6 * (E2T*E2T)*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
	      YT6 = T1Q + T4Q
	      TLN = TLN - DLOG4
	      YT1 = -YT1/URAURB  ! could store "1/uraurb=uranegiinv*urbnegiinv"
	      YT2 = 4*YT2
	      YT3 = YT3/URBNEGI  ! could store "1/urbnegi"
	      YT4 = YT4/URANEGI  ! could store "1/uranegi"
	    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
	        YT4 =  T4Q*CP  + T6Q*SDP
	      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
	        YT4 =  T4Q*CP  - T6Q*SDP
	      END IF    ! note that "row vector comes from the left"
	    END IF
C	    now the Cdelta propagation step
	    T2Q = YT2
	    YT2 = YT2+YT2           ! temporary doubling
	    YT6 = YT6 + (YT2-YT1*U)*U
	    T1Q = YT1*U
	    T6Q = YT6*G0U + T1Q
	    YT1 = YT1 + (YT2-T1Q-T6Q)*G0U
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>
	    YT2 = T2Q - T6Q
	    YT1 = RH1*YT1  ! "special"
	    YT6 = RH2*YT6  ! "special"

	  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

	  CALL SCALE5EC ( tspmaxlim1,tspmaxlim2, TLN,YT1,YT2,YT3,YT4,YT6 )

	  IF (SETZLOK) THEN

C	    now z propagation

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

C	      now the Cinv propagation step, with "special" first, 
C	      neglecting "final rh1 vector-multiplication" which could be done 
C	      simply by adding to tln a stored "dlog(rh1)" <OBS>
	      Z1 = RH2 * Z1  ! "special"
	      Z2 = RH2 * Z2  ! "special"
	      Z1 = Z1 - Z4*G0U
	      Z4 = Z4 + Z1*U
	      Z2 = Z2 - Z3*G0U
	      Z3 = Z3 + U*Z2
	      IF (EXPTYP) THEN
C	        now the Dinv propagation step
	        S1Q = Z1*URBNEGI
	        S2Q = Z2*URANEGI
	        Z1 = S1Q + Z3
	        Z2 = S2Q + Z4
	        Z3 = S1Q - Z3
	        Z4 = -S2Q + Z4
	        IF (IYSTAB.GT.0) THEN
C	          note: "special" above, to match that for yt
	          IF (ZLN.NE.TLNFW(ID)) STOP 34156
	          Z1FW(ID) = Z1
	          Z2FW(ID) = Z2
	          Z3FW(ID) = Z3
	          Z4FW(ID) = Z4
	        END IF
C	        now the PE propagation step, freq-dependent
	        ZLN = ZLN + DABS(ZLNREF)
	        IF (NRIKT.LT.0) THEN
	          Z1 = Z1 * DREAL(SDQ)*CQ
	          Z2 = Z2 * DREAL(SDP)*CP
	          Z3 = Z3 * DCONJG(CQ)
	          Z4 = Z4 * DCONJG(CP)
	        ELSE
	          Z1 = Z1 * DCONJG(CQ)
	          Z2 = Z2 * DCONJG(CP)
	          Z3 = Z3 * DREAL(SDQ)*CQ
	          Z4 = Z4 * DREAL(SDP)*CP
	        END IF    ! note that "row vector comes from the left"
	        IF (ZLNREF.GE.0.0D0) THEN
	          Z1 = E1X*Z1
	          Z3 = E1X*Z3
	        ELSE
	          Z2 = E1X*Z2
	          Z4 = E1X*Z4
	        END IF
C	        now the D propagation step
	        ZLN = ZLN - DLOG2
	        S1Q = Z1 + Z3
	        S2Q = Z2 - Z4
	        Z3 = Z1 - Z3
	        Z4 = Z2 + Z4
	        Z1 = S1Q/URBNEGI   ! "1/urbnegi" could be stored
	        Z2 = S2Q/URANEGI   ! "1/uranegi" could be stored
	      ELSE
C	        now the PT propagation step, freq-dependent
	        ZLN = ZLN + DABS(ZLNREF)
	        IF (NRIKT.LT.0) THEN
	          S1Q = Z1*CQ - Z3*SDQ
	          S2Q = Z2*CP - Z4*SDP
	          Z3 = Z1*STQ + Z3*CQ
	          Z4 = Z2*STP + Z4*CP
	        ELSE
	          S1Q = Z1*CQ + Z3*SDQ
	          S2Q = Z2*CP + Z4*SDP
	          Z3 = -Z1*STQ + Z3*CQ
	          Z4 = -Z2*STP + Z4*CP
	        END IF    ! note that "row vector comes from the left"
	        Z1 = S1Q
	        Z2 = S2Q
	        IF (ZLNREF.GE.0.0D0) THEN
	          Z1 = E1X*Z1
	          Z3 = E1X*Z3
	        ELSE
	          Z2 = E1X*Z2
	          Z4 = E1X*Z4
	        END IF
	      END IF
C	      now the C propagation step
	      Z3 = Z3 - Z2*U
	      Z4 = Z4 - Z1*U
	      Z1 = Z1 + Z4*G0U
	      Z2 = Z2 + Z3*G0U
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>
	      Z1 = Z1 * RH1  ! "special"
	      Z2 = Z2 * RH1  ! "special"

	    ELSE

	      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
	      Z(1) = Z1
	      Z(2) = Z2
	      Z(3) = Z3
	      Z(4) = Z4
	      CALL DPCLIN ( FSUBZSOL4,4, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	      HFAKTMAX,HOFF,EPSMAGNIFMAX, ZENDFSUB, Z0FSUB,ZLN,Z,
     $	      EPSFSUB,HSTEPFSUB, IND )
	      IF (IND.NE.0) STOP 6760
	      Z1 = Z(1)
	      Z2 = Z(2)
	      Z3 = Z(3)
	      Z4 = Z(4)

	    END IF

	    IF (ZLN.NE.TLN) THEN
	      SLASK = DEXP(ZLN-TLN)
	      ZLN = TLN
	      Z1 = SLASK * Z1
	      Z2 = SLASK * Z2
	      Z3 = SLASK * Z3
	      Z4 = SLASK * Z4
	    END IF
	    IF (NRIKT.GT.0) THEN
	      Z1 = Z1 + YT1*DIS2FW(ID) + 
     $	      YT2*DIS3FW(ID) + YT3*DIS4FW(ID)
	      Z2 = Z2 - YT1*DIS1FW(ID) + 
     $	      YT4*DIS3FW(ID) - YT2*DIS4FW(ID)   ! yt5=-yt2
	      Z3 = Z3 - YT2*DIS1FW(ID) - 
     $	      YT4*DIS2FW(ID) + YT6*DIS4FW(ID)
	      Z4 = Z4 - YT3*DIS1FW(ID) + 
     $	      YT2*DIS2FW(ID) - YT6*DIS3FW(ID)   ! yt5=-yt2
	    END IF

	  END IF

	END DO

	IF (OPTLOK) THEN
	  IF ( INDOPT.EQ.2  .OR. INDOPT.EQ.12 ) THEN
C	    cdetlok sign will fit
	    RLNDETLOK = TLN + TLNOPTLOK
	    CDETLOK = YT1*YT6OPTLOK + 2*YT2*YT2OPTLOK + 
     $	    YT3*YT4OPTLOK + YT4*YT3OPTLOK + YT6*YT1OPTLOK     ! yt5=-yt2
	    SLASK = DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	    ( 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)) ) * 
     $	    DMAX1 ( DABS(DREAL(YT1OPTLOK)),DABS(DIMAG(YT1OPTLOK)),
     $	    DABS(DREAL(YT2OPTLOK)),DABS(DIMAG(YT2OPTLOK)),
     $	    DABS(DREAL(YT3OPTLOK)),DABS(DIMAG(YT3OPTLOK)),
     $	    DABS(DREAL(YT4OPTLOK)),DABS(DIMAG(YT4OPTLOK)),
     $	    DABS(DREAL(YT6OPTLOK)),DABS(DIMAG(YT6OPTLOK)) ) )
	    IF (SLASK.LT.QIOPT) THEN
	      IF (NRIKT.GT.0) THEN
	        NSOOPT = NL2
	      ELSE
	        NSOOPT = NL2 - 1
	      END IF
	      QIOPT = SLASK
	    END IF
	  END IF
	END IF

	RETURN
	END








	SUBROUTINE FWSUBYTZLIQ ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, DIS2FW,DIS4FW, 
     $	TLN, YT1,YT3, Z1, 
     $	TLNFW, YT1FW,YT3FW, Z1FW )
C	additional "fw input & output" in "common /fwparam,fwpart/"
C	propagation of "tlnfw, yt1fw,yt3fw, z1fw"
C	needed for "zliqsifw" and the stabilization in "liqsifw"
C	note: -> propagation of "row-vector from the left" 
C	         in particular: "downward matrices when nrikt=-1"
C	                        "  upward matrices when nrikt=+1"
C	         (however, dpclin & fsubng2)
C	      -> for yt; homogeneous case
C	         for z; inhomogeneous case

	USE INMODELa
	USE FWPARTa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C       nol <= noldim

	INCLUDE '../DPCFSUB.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE 'FWPARAM.INC'

	INTEGER NL1,NL2,NRIKT
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2, DIS2FW(*), DIS4FW(*)
	LOGICAL INIT

	REAL*8 TLN                 ! input/output
	COMPLEX*16 YT1,YT3, Z1     ! input/output

	REAL*8 TLNFW(*)
	COMPLEX*16 YT1FW(*), YT3FW(*), 
     $	Z1FW(*)

	INTEGER INDOPT                    ! input (0, 1,2,12)
	INTEGER NSOOPT                    ! input/output
	REAL*8 QIOPT                      ! input/output
	REAL*8 TLNOPT(*)                  ! output if indopt.eq.0, else input
	COMPLEX*16 YT1OPT(*), YT3OPT(*)   ! output if indopt.eq.0, else input

	INTEGER ID, IDQ, IND, ISLASK
	REAL*8 TLNREF, ZLN, ZENDFSUB, Z0FSUB, RLNDETLOK, TLNOPTLOK, 
     $	SLASK, PI
	COMPLEX*16 URA, URANEGI, 
     $	CP,CPX, STP,SDP, WWLOK, 
     $	CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	CDETLOK, YT1OPTLOK,YT3OPTLOK, 
     $	YT(2)
	LOGICAL EXPTYP, SETZLOK, OPTLOK

	EXTERNAL FSUBNG2

	DATA PI /3.141592653589793D0/

	SETZLOK = .TRUE.
	OPTLOK = .FALSE.
	GOTO 10

	ENTRY FWSUBYTLIQ ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, 
     $	TLN, YT1,YT3, 
     $	TLNFW, YT1FW,YT3FW )

	SETZLOK = .FALSE.
	OPTLOK = .FALSE.
	GOTO 10

	ENTRY FWOPTYTZLIQ ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, DIS2FW,DIS4FW, 
     $	TLN, YT1,YT3, Z1, 
     $	TLNFW, YT1FW,YT3FW, Z1FW, 
     $	INDOPT, TLNOPT,YT1OPT,YT3OPT, NSOOPT,QIOPT )
C	note: for indopt.ne.0, same act arg for "ytopt","ytfw" ok (overwrite)

	SETZLOK = .TRUE.
	GOTO 20

	ENTRY FWOPTYTLIQ ( NL1,NL2,NRIKT, OMEG,OMEG2,U,U2,RK,RK2, 
     $	INIT, 
     $	TLN, YT1,YT3, 
     $	TLNFW, YT1FW,YT3FW, 
     $	INDOPT, TLNOPT,YT1OPT,YT3OPT, NSOOPT,QIOPT )
C	note: for indopt.ne.0, same act arg for "ytopt","ytfw" ok (overwrite)

	SETZLOK = .FALSE.
 20	OPTLOK = .TRUE.

	IF ( INDOPT.EQ.1  .OR. INDOPT.EQ.12 ) THEN
	  IDQ = NL1 - NRIKT
	  TLNOPTLOK = TLNOPT(IDQ)
	  YT1OPTLOK = YT1OPT(IDQ)
	  YT3OPTLOK = YT3OPT(IDQ)
C	  first - for cdetlok to fit with liqsi, (yp1,yp2) <-> +-(-yt3,yt1)
	  RLNDETLOK = TLN + TLNOPTLOK
	  CDETLOK = - ( YT1*YT3OPTLOK - YT1OPTLOK*YT3 )
	  SLASK = DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	  ( DMAX1 ( DABS(DREAL(YT1)),DABS(DIMAG(YT1)),
     $	  DABS(DREAL(YT3)),DABS(DIMAG(YT3)) ) * 
     $	  DMAX1 ( DABS(DREAL(YT1OPTLOK)),DABS(DIMAG(YT1OPTLOK)),
     $	  DABS(DREAL(YT3OPTLOK)),DABS(DIMAG(YT3OPTLOK)) ) )
	  IF (SLASK.LT.QIOPT) THEN
	    IF (NRIKT.GT.0) THEN
	      NSOOPT = NL1 - 1
	    ELSE
	      NSOOPT = NL1
	    END IF
	    QIOPT = SLASK
	  END IF
	END IF

 10	CONTINUE

	DO ID = NL1, NL2, NRIKT

	  IF (OPTLOK) THEN
	    IF (INDOPT.EQ.0) THEN
	      TLNOPT(ID) = TLN
	      YT1OPT(ID) = YT1
	      YT3OPT(ID) = YT3
	    ELSE
	      IF (ID.NE.NL1) THEN
C	        first - for cdetlok to fit with liqsi, 
C	        (yp1,yp2) <-> +-(-yt3,yt1)
	        RLNDETLOK = TLN + TLNOPTLOK
	        CDETLOK = - ( YT1*YT3OPTLOK - YT1OPTLOK*YT3 )
	        SLASK = 
     $	        DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	        ( DMAX1 ( DABS(DREAL(YT1)),DABS(DIMAG(YT1)),
     $	        DABS(DREAL(YT3)),DABS(DIMAG(YT3)) ) * 
     $	        DMAX1 ( DABS(DREAL(YT1OPTLOK)),DABS(DIMAG(YT1OPTLOK)),
     $	        DABS(DREAL(YT3OPTLOK)),DABS(DIMAG(YT3OPTLOK)) ) )
	        IF (SLASK.LT.QIOPT) THEN
	          IF (NRIKT.GT.0) THEN
	            NSOOPT = ID - 1
	          ELSE
	            NSOOPT = ID
	          END IF
	          QIOPT = SLASK
	        END IF
	      END IF
	      TLNOPTLOK = TLNOPT(ID)
	      YT1OPTLOK = YT1OPT(ID)
	      YT3OPTLOK = YT3OPT(ID)
	    END IF
	  END IF

C	  preliminary setting that can be overwritten below
	  TLNFW(ID) = TLN
	  YT1FW(ID) = YT1
	  YT3FW(ID) = YT3

	  IF (SETZLOK) THEN
	    ZLN = TLN
	    Z1FW(ID) = Z1   ! preliminary setting that can be overwritten below
	    IF (NRIKT.LT.0) Z1 = Z1 - YT1*DIS2FW(ID) - YT3*DIS4FW(ID)
	  END IF

C	  now yt propagation

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

	    IF ( IFSMTYP(ID).NE.-6 .AND. IFSMTYP(ID).NE.(NRIKT-6) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
	      IF (INIT) THEN
	        CALL LIQPART ( OMEG,OMEG2,U,U2, UA2(ID), DZ(ID), 
     $	        URA, EXPTYP, TLNREF, CP, STP,SDP )
	        URA_FW(ID) = URA
	        EXPTYP_FW(ID) = EXPTYP
	        TLNREF_FW(ID) = TLNREF
	        CP_FW(ID) = CP
	        STP_FW(ID) = STP
	        SDP_FW(ID) = SDP
	      ELSE
	        URA = URA_FW(ID)
	        EXPTYP = EXPTYP_FW(ID)
	        TLNREF = TLNREF_FW(ID)
	        CP = CP_FW(ID)
	        STP = STP_FW(ID)
	        SDP = SDP_FW(ID)
	      END IF
	      IF (EXPTYP) URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura
	    ELSE IF (INIT) THEN
	      URA_FW(ID) = URA_FW(ID-NRIKT)
	      EXPTYP_FW(ID) = EXPTYP_FW(ID-NRIKT)
	      TLNREF_FW(ID) = TLNREF_FW(ID-NRIKT)
	      CP_FW(ID) = CP_FW(ID-NRIKT)
	      STP_FW(ID) = STP_FW(ID-NRIKT)
	      SDP_FW(ID) = SDP_FW(ID-NRIKT)
	    END IF
	    YT3 = YT3*RHO(ID)     ! already here
	    IF (EXPTYP) THEN
C	      now the first propagation step
	      WWLOK = YT1*URANEGI
	      YT1 = WWLOK + YT3
	      YT3 = -WWLOK + YT3
	      IF (IYSTAB.GT.0) THEN
	        TLNFW(ID) = TLN
	        YT1FW(ID) = YT1
	        YT3FW(ID) = YT3
	      END IF
C	      now the main propagation step, freq-dependent
	      TLN = TLN + TLNREF - DLOG2       ! -dlog2 from "last step"
	      IF (NRIKT.LT.0) THEN
	        YT1 = YT1 * DREAL(SDP)*CP
	        YT3 = YT3 * DCONJG(CP)
	      ELSE
	        YT1 = YT1 * DCONJG(CP)
	        YT3 = YT3 * DREAL(SDP)*CP
	      END IF    ! note that "row vector comes from the left"
C	      now the next propagation step
	      WWLOK = YT1 - YT3
	      YT3 = ( YT1 + YT3 ) / RHO(ID)      ! "1/rho(id)" could be stored
	      YT1 = WWLOK / URANEGI
	    ELSE
C	      now the main and only propagation step, freq-dependent
	      TLN = TLN + TLNREF
	      IF (NRIKT.LT.0) THEN
	        WWLOK = YT1*CP - YT3*SDP
	        YT3 = ( YT1*STP + YT3*CP ) / RHO(ID)  !"1/rho" could be stored
	        YT1 = WWLOK
	      ELSE
	        WWLOK = YT1*CP + YT3*SDP
	        YT3 = ( -YT1*STP + YT3*CP ) / RHO(ID) !"1/rho" could be stored
	        YT1 = WWLOK
	      END IF    ! note that "row vector comes from the left"
	    END IF

	  ELSE IF (IFSMTYP(ID).LE.0) THEN

	    IF (INIT) THEN
	      IF ( IFSMTYP(ID).NE.-2 .AND. IFSMTYP(ID).NE.(NRIKT-2) ) THEN
C	        warning: not allowed to skip for "first id" in non-empty loop
	        ISLASK = 0
	      ELSE
	        ISLASK = NRIKT
	      END IF
	      CALL AIRYPART ( OMEG,OMEG2,U,U2,RK,RK2, 
     $	      UA2(ID),UB2(ID), DZ(ID), ISLASK, 
     $	      CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	      TLNREF, CP,CPX, STP,SDP )
	      CFAKT_FW(ID) = CFAKT
	      BCFAKT_FW(ID) = BCFAKT
	      BCFAKTINV_FW(ID) = BCFAKTINV
	      YA_FW(ID) = YA
	      YB_FW(ID) = YB
	      TLNREF_FW(ID) = TLNREF
	      CP_FW(ID) = CP
	      CPX_FW(ID) = CPX
	      STP_FW(ID) = STP
	      SDP_FW(ID) = SDP
	    ELSE
	      CFAKT = CFAKT_FW(ID)
	      BCFAKT = BCFAKT_FW(ID)
	      BCFAKTINV = BCFAKTINV_FW(ID)
	      YA = YA_FW(ID)
	      YB = YB_FW(ID)
	      TLNREF = TLNREF_FW(ID)
	      CP = CP_FW(ID)
	      CPX = CPX_FW(ID)
	      STP = STP_FW(ID)
	      SDP = SDP_FW(ID)
	    END IF
C	    now propagation
	    YT3 = YT3*RHO(ID)     ! already here
	    TLN = TLN + TLNREF
	    IF (NRIKT.LT.0) THEN
	      WWLOK = PI * ( YT1*CP - YT3*SDP )
	      YT3 = (PI/RHO(ID)) * ( YT1*STP + YT3*CPX )  !"pi/rho" could be st
	      YT1 = WWLOK
	    ELSE
	      WWLOK = PI * ( YT1*CPX + YT3*SDP )
	      YT3 = (PI/RHO(ID)) * ( -YT1*STP + YT3*CP )  !"pi/rho" could be st
	      YT1 = WWLOK
	    END IF    ! note that "row vector comes from the left"

	  ELSE

	    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
C	    now propagation
	    YT(1) = YT1
	    YT(2) = YT3
	    CALL DPCLIN ( FSUBNG2,2, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	    HFAKTMAX,HOFF,EPSMAGNIFMAX, ZENDFSUB, Z0FSUB,TLN,
     $	    YT,EPSFSUB,HSTEPFSUB, IND )
	    IF (IND.NE.0) THEN
	      WRITE(*,*) ' %% u,id ',UCOP,ID
	      STOP 6770
	    END IF
	    YT1 = YT(1)
	    YT3 = YT(2)

	  END IF

	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLN,YT1,YT3 )

	  IF (SETZLOK) THEN
C	    now z propagation
	    IF (IFSMTYP(ID).LE.-4) THEN
	      IF (EXPTYP) THEN
	        IF (IYSTAB.GT.0) THEN
	          IF (ZLN.NE.TLNFW(ID)) STOP 34156
	          Z1FW(ID) = Z1
	        END IF
	      END IF
	    END IF
	    IF (ZLN.NE.TLN) THEN
	      SLASK = DEXP(ZLN-TLN)
	      ZLN = TLN
	      Z1 = SLASK * Z1
	    END IF
	    IF (NRIKT.GT.0) Z1 = Z1 + YT1*DIS2FW(ID) + YT3*DIS4FW(ID)
	  END IF

	END DO

	IF (OPTLOK) THEN
	  IF ( INDOPT.EQ.2 .OR. INDOPT.EQ.12 ) THEN
C	    first - for cdetlok to fit with liqsi, (yp1,yp2) <-> +-(-yt3,yt1)
	    RLNDETLOK = TLN + TLNOPTLOK
	    CDETLOK = - ( YT1*YT3OPTLOK - YT1OPTLOK*YT3 )
	    SLASK = DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	    ( DMAX1 ( DABS(DREAL(YT1)),DABS(DIMAG(YT1)),
     $	    DABS(DREAL(YT3)),DABS(DIMAG(YT3)) ) * 
     $	    DMAX1 ( DABS(DREAL(YT1OPTLOK)),DABS(DIMAG(YT1OPTLOK)),
     $	    DABS(DREAL(YT3OPTLOK)),DABS(DIMAG(YT3OPTLOK)) ) )
	    IF (SLASK.LT.QIOPT) THEN
	      IF (NRIKT.GT.0) THEN
	        NSOOPT = NL2
	      ELSE
	        NSOOPT = NL2 - 1
	      END IF
	      QIOPT = SLASK
	    END IF
	  END IF
	END IF

	RETURN
	END

