
C	additional subroutines (and functions) for "Full Wave" computations 
C	special "fw" versions :  undoxsol(0)  undoxliq

C	Sven Ivansson  1994






	SUBROUTINE UNDOXSOL ( U, ID, YT1,YT2,YT3,YT4,YT6, NRIKT, Z1,Z2,Z3,Z4 )

	USE INMODELa
	USE FWPARTa
	USE FWDISa

	IMPLICIT NONE

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

	INTEGER ID, NRIKT
	COMPLEX*16 U

	COMPLEX*16 YT1,YT2,YT3,YT4,YT6, Z1,Z2,Z3,Z4

	REAL*8 RH1,RH2
	COMPLEX*16 G0, G0U, 
     $	URA,URB, URAURB,URANEGI,URBNEGI, 
     $	T1Q,T2Q,T3Q,T4Q,T6Q, S1Q,S2Q
	LOGICAL SETZ, EXPTYP

	SETZ = .TRUE.
	GOTO 100

	ENTRY UNDOXSOL0 ( U, ID, YT1,YT2,YT3,YT4,YT6 )

	SETZ = .FALSE.

 100	CONTINUE

	IF (IFSMTYP(ID).GT.-4) STOP 34156

	RH1 = RHO(ID)                        ! local rho(id)
	RH2 = 1/RH1                          ! could be stored
	G0 = TWOVB2(ID)                      ! local twovb2(id)
	G0U = G0*U
	URA = URA_FW(ID)
	URB = URB_FW(ID)
	EXPTYP = EXPTYP_FW(ID)
	IF (.NOT.EXPTYP) STOP 341555
	URAURB = URA*URB
	URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb

C	now yt propagation

C	now the Ddelta 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
	YT1 = -YT1 / (4*URAURB)
	YT2 = -YT2 / (2*URBNEGI)
	YT3 = YT3 / (4*URBNEGI)
	YT4 = YT4 / (4*URANEGI)
	YT6 = YT6 / 4
C	now the "rh1*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
	YT2 = T2Q - T6Q
	YT1 = RH1*YT1
	YT6 = RH2*YT6

	IF (SETZ) THEN

C	  now z propagation

C	  now the D propagation step
	  S1Q = Z1 + Z3
	  S2Q = Z2 - Z4
	  Z3 = ( Z1 - Z3 ) / 2
	  Z4 = ( Z2 + Z4 ) / 2
	  Z1 = S1Q / (2*URBNEGI)
	  Z2 = S2Q / (2*URANEGI)
C	  now the "rh1*C" propagation step
	  Z3 = Z3 - Z2*U
	  Z4 = Z4 - Z1*U
	  Z1 = Z1 + Z4*G0U
	  Z2 = Z2 + Z3*G0U
	  Z1 = Z1 * RH1
	  Z2 = Z2 * RH1

	  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

	RETURN
	END









	SUBROUTINE UNDOXLIQ ( ID, YT1,YT3, NRIKT, Z1 )

	USE INMODELa
	USE FWPARTa
	USE FWDISa

	IMPLICIT NONE

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

	INTEGER ID, NRIKT

	COMPLEX*16 YT1,YT3, Z1

	COMPLEX*16 URA, URANEGI, WWLOK
	LOGICAL SETZ, EXPTYP

	SETZ = .TRUE.
	GOTO 100

	ENTRY UNDOXLIQ0 ( ID, YT1,YT3 )

	SETZ = .FALSE.

 100	CONTINUE

	IF (IFSMTYP(ID).GT.-4) STOP 344451

	URA = URA_FW(ID)
	EXPTYP = EXPTYP_FW(ID)
	IF (.NOT.EXPTYP) STOP 345155
	URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura

C	now yt propagation

C	now the next propagation step
	WWLOK = YT1 - YT3
	YT3 = ( YT1 + YT3 ) / (2*RHO(ID))
	YT1 = WWLOK / (2*URANEGI)


	IF (SETZ) THEN

C	  now z propagation

	  IF (NRIKT.GT.0) Z1 = Z1 + YT1*DIS2FW(ID) + YT3*DIS4FW(ID)

	END IF

	RETURN
	END
