
C	additional subroutines for "Full Wave" computations
C	dhivanfw2

C	Sven Ivansson  1992






	SUBROUTINE DHIVANFW2 (UP)

	USE INMODELa
	USE FWYSa
	USE FWPARTa
	USE FWYTa

	IMPLICIT NONE

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

	INCLUDE '../INPARAM.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE '../INSPEC.INC'
	INCLUDE '../DDULOK.INC'
	INCLUDE '../fw/FWPARAM.INC'
	INCLUDE '../fw/FWINTF.INC'
	INCLUDE '../fw/FWNSO.INC'

	COMPLEX*16 UP

	INTEGER NSOORIG, ICLOK, IFSMTYPTMP, ISLASK
	REAL*8 QIOPT, TLNFWNOL, TLNWF0, ALN,BLN, SLASK, 
     &	SLASK1,SLASK2,SLASK3,SLASK4
	COMPLEX*16 U,RK,UP2,U2,RK2, 
     $	URA,URB, URANEGI,URBNEGI, G0,G0U, 
     $	YT1FWNOL,YT2FWNOL,YT3FWNOL,YT4FWNOL,YT6FWNOL, 
     $	YT1WF0,YT2WF0,YT3WF0,YT4WF0,YT6WF0, 
     $	CSLASK, CSLASK1,CSLASK2,CSLASK3,CSLASK4, URVALUE, 
     &	BFAKT,CFAKT,BCFAKT, YY, AIE,AIDE,BIE,BIDE, C60,C120

	DATA C60  /(+0.5D0,+0.8660254037844386D0)/  ! cdexp(i*pi/3)
	DATA C120 /(-0.5D0,+0.8660254037844386D0)/  ! cdexp(i*2*pi/3)

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

	UCOP = U
	U2COP = U2
	RK2GOMEG = RK2/OMEG

C	downward treatment of the upper solid if present
	IF (RHO(0).GT.0.0D0) 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
	  URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura
	  IF (ISOLID(0).NE.0) THEN
	    URB = URVALUE (ICUTTYP,AIMCUT, U0,UPBPBSURF,UPBPB2SURF,UPREF, 
     $      UP,UP2)
	    IF ( ISHEETTYP.EQ.3 .OR. ISHEETTYP.EQ.4 ) URB = -URB
	    URB = -URB      ! up/down viceversa for above and below
	    URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))  ! -i*urb
	  END IF
	END IF
	IF  (NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	  IF (RHO(0).GT.0.0D0) THEN
	    IF (ISOLID(0).NE.0) THEN
C	      solid half-space above
	      G0 = TWOVB2(0)
	      G0U = G0*U
	      CSLASK = -URA*URB - U2
	      CSLASK2 = -CSLASK*G0U - U
	      TLNFWANSS = 0.0D0
	      YT1FWANSS = ( 1 + (CSLASK2-U)*G0U ) * RHO(0)
	      YT2FWANSS = CSLASK2
	      YT3FWANSS = -URANEGI
	      YT4FWANSS = URBNEGI
	      YT6FWANSS = CSLASK / RHO(0)
	      IF (NSW.EQ.0) THEN
C               special 'fix' for this case without 'nss-nsw' fluid
	        TLNFWBNSW = TLNFWANSS
	        YT1FWBNSW = YT1FWANSS
	        YT2FWBNSW = YT2FWANSS
	        YT3FWBNSW = YT3FWANSS
	        YT4FWBNSW = YT4FWANSS
	        YT6FWBNSW = YT6FWANSS
	        GOTO 111
	      END IF
	    ELSE
C	      fluid half-space above upper solid
	      TLNFWANSS = 0.0D0
	      YT1FWANSS = (0.0D0,0.0D0)
	      YT2FWANSS = (0.0D0,0.0D0)
	      YT3FWANSS = (0.0D0,0.0D0)
	      YT4FWANSS = RHO(0)
	      YT6FWANSS = URANEGI
	    END IF
	  ELSE IF (RHO(0).EQ.-99.0D0) THEN
	    STOP '%%% dhivanfw2, "multi"'
	  ELSE IF (RHO(0).LT.0.0D0) THEN
C	    rigid above upper solid
	    TLNFWANSS = 0.0D0
	    YT1FWANSS = (1.0D0,0.0D0)
	    YT2FWANSS = (0.0D0,0.0D0)
	    YT3FWANSS = (0.0D0,0.0D0)
	    YT4FWANSS = (0.0D0,0.0D0)
	    YT6FWANSS = (0.0D0,0.0D0)
	  ELSE
C	    free above upper solid
	    TLNFWANSS = 0.0D0
	    YT1FWANSS = (0.0D0,0.0D0)
	    YT2FWANSS = (0.0D0,0.0D0)
	    YT3FWANSS = (0.0D0,0.0D0)
	    YT4FWANSS = (0.0D0,0.0D0)
	    YT6FWANSS = (1.0D0,0.0D0)
	  END IF
	  IF (NSS.GT.0) THEN
C	    make the call
	    CALL FWOPTYTSOL ( 1,NSS,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    .TRUE., 
     $	    TLNFWANSS,YT1FWANSS,YT2FWANSS,YT3FWANSS,YT4FWANSS,YT6FWANSS, 
     $	    TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, 
     $	    0, TLNOPT,YT1OPT,YT2OPT,YT3OPT,YT4OPT,YT6OPT, NSO,QIOPT )
	  END IF
C	  make appropriate settings for "below the nss level" (see "29/1 98")
	  TLNFWBNSS = TLNFWANSS
	  YT1FWBNSS = -YT1FWANSS
	  YT3FWBNSS = -YT3FWANSS
	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNFWBNSS,
     $	  YT1FWBNSS,YT3FWBNSS )
	ELSE
C	  nss.eq.0 and isolid(0).eq.0, set for "below nss level" ("29/1 98")
	  IF (RHO(0).GT.0.0D0) THEN
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      TLNFWBNSW = 0.0D0
	      YT1FWBNSW = (0.0D0,0.0D0)
	      YT2FWBNSW = (0.0D0,0.0D0)
	      YT3FWBNSW = (0.0D0,0.0D0)
	      YT4FWBNSW = RHO(0)
	      YT6FWBNSW = URANEGI
	      GOTO 111
	    END IF
C	    fluid half-space above fluid
	    TLNFWBNSS = 0.0D0
	    YT1FWBNSS = -RHO(0)
	    YT3FWBNSS = URANEGI
	  ELSE IF (RHO(0).EQ.-99.0D0) THEN
C	    adaptation from avh/multi/liqsi.f (cf articles/ua15bok/GF)
	    IF (NSW.LT.1) STOP '%%% dhivanfw2, "multi"'
	    IF ( IFSMTYP(1).LE.-4 .OR. IFSMTYP(1).GE.1 ) 
     &	     STOP '%%% dhivanfw2, "multi"'
	    BFAKT = -OMEG2 * UB2(1) / DZ(1)
	    IF (DIMAG(BFAKT).EQ.0.0D0) THEN
	      SLASK = DABS(DREAL(BFAKT))                    ! "dabs" (branch)
	      CFAKT = DCMPLX ( SLASK**(-0.6666666666666667D0) , 0.0D0 )
	      IF (SLASK.EQ.0.0D0) STOP '%%% dhivanfw2, "multi" slask'
	    ELSE
	      IF (DREAL(BFAKT).GE.0.0D0) THEN
	        CFAKT = BFAKT**(-0.6666666666666667D0)
	      ELSE
	        CFAKT = (-BFAKT)**(-0.6666666666666667D0)   ! "-" (branch)
	      END IF
	    END IF
	    BCFAKT = -BFAKT*CFAKT / OMEG
	    YY = CFAKT * ( RK2 - OMEG2*(UA2(1)-0.5D0*UB2(1)) )
	    IF (DABS(DIMAG(BCFAKT)).GE.DSQRT(3.0D0)*DABS(DREAL(BCFAKT))) 
     $	     STOP '%%% dhivanfw2, "multi" bcfakt'
	    IF ( DREAL(BFAKT).GT.0.0D0 ) THEN
	      IF ( ICONTCUT0.NE.0 .OR. DIMAG(BFAKT).GT.0.0D0 .OR. 
     $	      ( DIMAG(BFAKT).EQ.0.0D0 .AND. DIMAG(RK2).LE.0.0D0 ) ) THEN
		YY = YY * C120
	        CALL AIRY4E (YY, ALN,BLN, AIE,AIDE,BIE,BIDE)
	        TLNFWBNSS = -ALN
	        YT1FWBNSS = -2 * DCONJG(C60) * AIE
	        YT3FWBNSS = 2 * C60 * AIDE * BCFAKT / RHO(1)
	      ELSE IF ( DIMAG(BFAKT).EQ.0.0D0 .AND. DIMAG(RK2).GT.0.0D0 ) THEN
		YY = YY * DCONJG(C120)
	        CALL AIRY4E (YY, ALN,BLN, AIE,AIDE,BIE,BIDE)
	        TLNFWBNSS = -ALN
	        YT1FWBNSS = -2 * C60 * AIE
	        YT3FWBNSS = 2 * DCONJG(C60) * AIDE * BCFAKT / RHO(1)
	      ELSE
	        STOP '%%% dhivanfw2, "multi" icontcut0,bfakt,rk2'
	      END IF
	    ELSE IF (DREAL(BFAKT).LT.0.0D0) THEN
	      CALL AIRY4E (YY, ALN,BLN, AIE,AIDE,BIE,BIDE)
	      TLNFWBNSS = -ALN
	      YT1FWBNSS = -AIE
	      YT3FWBNSS = BCFAKT * AIDE / RHO(1)
	    ELSE
	      STOP '%%% dhivanfw2, "multi" bfakt'
	    END IF
	  ELSE IF (RHO(0).LT.0.0D0) THEN
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      TLNFWBNSW = 0.0D0
	      YT1FWBNSW = (1.0D0,0.0D0)
	      YT2FWBNSW = (0.0D0,0.0D0)
	      YT3FWBNSW = (0.0D0,0.0D0)
	      YT4FWBNSW = (0.0D0,0.0D0)
	      YT6FWBNSW = (0.0D0,0.0D0)
	      GOTO 111
	    END IF
C	    rigid above fluid
	    TLNFWBNSS = 0.0D0
	    YT1FWBNSS = (-1.0D0,0.0D0)
	    YT3FWBNSS = (0.0D0,0.0D0)
	  ELSE	      
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      TLNFWBNSW = 0.0D0
	      YT1FWBNSW = (0.0D0,0.0D0)
	      YT2FWBNSW = (0.0D0,0.0D0)
	      YT3FWBNSW = (0.0D0,0.0D0)
	      YT4FWBNSW = (0.0D0,0.0D0)
	      YT6FWBNSW = (1.0D0,0.0D0)
	      GOTO 111
	    END IF
C	    free above fluid
	    TLNFWBNSS = 0.0D0
	    YT1FWBNSS = (0.0D0,0.0D0)
	    YT3FWBNSS = (1.0D0,0.0D0)
	  END IF
	END IF

C	downward treatment of the fluid
	IF (NSW.GT.NSS) THEN
	  TLNFWANSW = TLNFWBNSS
	  YT1FWANSW = YT1FWBNSS
	  YT3FWANSW = YT3FWBNSS
C	  make the call
	  CALL FWOPTYTLIQ ( NSS+1,NSW,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	  .TRUE., 
     $	  TLNFWANSW,YT1FWANSW,YT3FWANSW, 
     $	  TLNFW, YT1FW,YT3FW, 
     $	  0, TLNOPT,YT1OPT,YT3OPT, NSO,QIOPT )
	ELSE
	  STOP 64430  ! should not occur
	END IF
	IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
C	  settings for "below the nsw level" (see "29/1 98")
	  TLNFWBNSW = TLNFWANSW
	  YT1FWBNSW = (0.0D0,0.0D0)
	  YT2FWBNSW = (0.0D0,0.0D0)
	  YT3FWBNSW = (0.0D0,0.0D0)
	  YT4FWBNSW = -YT1FWANSW
	  YT6FWBNSW = YT3FWANSW
	END IF

C	downward treatment of the lower solid if present
 111	IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	  TLNFWNOL = TLNFWBNSW
	  YT1FWNOL = YT1FWBNSW
	  YT2FWNOL = YT2FWBNSW
	  YT3FWNOL = YT3FWBNSW
	  YT4FWNOL = YT4FWBNSW
	  YT6FWNOL = YT6FWBNSW
	  IF ( (NSW+1) .LT. NOL ) THEN
C	    make the call
	    CALL FWOPTYTSOL ( NSW+1,NOL-1,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    .TRUE., 
     $	    TLNFWNOL,YT1FWNOL,YT2FWNOL,YT3FWNOL,YT4FWNOL,YT6FWNOL, 
     $	    TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, 
     $	    0, TLNOPT,YT1OPT,YT2OPT,YT3OPT,YT4OPT,YT6OPT, NSO,QIOPT )
	  END IF
	  TLNOPT(NOL) = TLNFWNOL
	  YT1OPT(NOL) = YT1FWNOL
	  YT2OPT(NOL) = YT2FWNOL
	  YT3OPT(NOL) = YT3FWNOL
	  YT4OPT(NOL) = YT4FWNOL
	  YT6OPT(NOL) = YT6FWNOL
	END IF

C	during the upward treatment below, set "new nso"
	NSOORIG = NSO
	NSO = -1
	QIOPT = 10.0D0

C	upward treatment of the lower solid if present
	IF (RHO(NOL).GT.0.0D0) 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
	  URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura
	  IF (ISOLID(NOL).NE.0) THEN
	    URB = URVALUE (ICUTTYP,AIMCUT, U0,UPBPB,UPBPB2,UPREF, UP,UP2)
	    IF ( ISHEETTYP.EQ.3 .OR. ISHEETTYP.EQ.4 ) URB = -URB
	    URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))  ! -i*urb
	  END IF
	END IF
	IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	  IF (RHO(NOL).GT.0.0D0) THEN
	    IF (ISOLID(NOL).NE.0) THEN
C	      solid half-space below
	      G0 = TWOVB2(NOL)
	      G0U = G0*U
	      CSLASK = -URA*URB - U2
	      CSLASK2 = -CSLASK*G0U - U
	      TLNWFBNSW = 0.0D0
	      YT1WFBNSW = ( 1 + (CSLASK2-U)*G0U ) * RHO(NOL)
	      YT2WFBNSW = CSLASK2
	      YT3WFBNSW = -URANEGI
	      YT4WFBNSW = URBNEGI
	      YT6WFBNSW = CSLASK / RHO(NOL)
	    ELSE
C	      fluid half-space below lower solid
	      TLNWFBNSW = 0.0D0
	      YT1WFBNSW = (0.0D0,0.0D0)
	      YT2WFBNSW = (0.0D0,0.0D0)
	      YT3WFBNSW = (0.0D0,0.0D0)
	      YT4WFBNSW = RHO(NOL)
	      YT6WFBNSW = URANEGI
	    END IF
	  ELSE IF (RHO(NOL).EQ.-99.0D0) THEN
	    STOP '%%% dhivanfw2, local-reaction probl'
	  ELSE IF (RHO(NOL).LT.0.0D0) THEN
C	    rigid below lower solid
	    TLNWFBNSW = 0.0D0
	    YT1WFBNSW = (1.0D0,0.0D0)
	    YT2WFBNSW = (0.0D0,0.0D0)
	    YT3WFBNSW = (0.0D0,0.0D0)
	    YT4WFBNSW = (0.0D0,0.0D0)
	    YT6WFBNSW = (0.0D0,0.0D0)
	  ELSE
C	    free below lower solid
	    TLNWFBNSW = 0.0D0
	    YT1WFBNSW = (0.0D0,0.0D0)
	    YT2WFBNSW = (0.0D0,0.0D0)
	    YT3WFBNSW = (0.0D0,0.0D0)
	    YT4WFBNSW = (0.0D0,0.0D0)
	    YT6WFBNSW = (1.0D0,0.0D0)
	  END IF
	  IF (NSW.EQ.0) THEN
C           special 'fix' for this case without 'nss-nsw' fluid
C	    make the call
	    CALL FWOPTYTSOL ( NOL-1,NSW+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    .FALSE., 
     $	    TLNWFBNSW,YT1WFBNSW,YT2WFBNSW,YT3WFBNSW,YT4WFBNSW,YT6WFBNSW, 
     $	    TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF, 
     $	    12, TLNOPT,YT1OPT,YT2OPT,YT3OPT,YT4OPT,YT6OPT, NSO,QIOPT )
	    YT1WFBNSO = YT1WFBNSW
	    YT2WFBNSO = YT2WFBNSW
	    YT3WFBNSO = YT3WFBNSW
	    YT4WFBNSO = YT4WFBNSW
	    YT6WFBNSO = YT6WFBNSW
	    GOTO 222
	  END IF
	  IF ( (NSW+1) .LT. NOL ) THEN
C	    make the call
	    CALL FWOPTYTSOL ( NOL-1,NSW+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    .FALSE., 
     $	    TLNWFBNSW,YT1WFBNSW,YT2WFBNSW,YT3WFBNSW,YT4WFBNSW,YT6WFBNSW, 
     $	    TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF, 
     $	    1, TLNOPT,YT1OPT,YT2OPT,YT3OPT,YT4OPT,YT6OPT, NSO,QIOPT )
	  END IF
C	  make appropriate settings for "above the nsw level" (see "29/1 98")
	  TLNWFANSW = TLNWFBNSW
	  YT1WFANSW = YT1WFBNSW
	  YT3WFANSW = YT3WFBNSW
	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNWFANSW,
     $	  YT1WFANSW,YT3WFANSW )
	ELSE
C	  (nol-1).eq.nsw and isolid(nol).eq.0, set "above nsw levl" ("29/1 98")
	  IF (RHO(NOL).GT.0.0D0) THEN
C	    fluid half-space below fluid
	    TLNWFANSW = 0.0D0
	    YT1WFANSW = RHO(NOL)
	    YT3WFANSW = -URANEGI
	  ELSE IF (RHO(NOL).EQ.-99.0D0) THEN
C	    local reaction below fluid
	    TLNWFANSW = 0.0D0
	    YT1WFANSW = ZLOCREACTNOL  ! se pierce ovan s 107
	    YT3WFANSW = (0.0D0,+1.0D0)
	  ELSE IF (RHO(NOL).LT.0.0D0) THEN
C	    rigid below fluid
	    TLNWFANSW = 0.0D0
	    YT1WFANSW = (1.0D0,0.0D0)
	    YT3WFANSW = (0.0D0,0.0D0)
	  ELSE	      
C	    free below fluid
	    TLNWFANSW = 0.0D0
	    YT1WFANSW = (0.0D0,0.0D0)
	    YT3WFANSW = (-1.0D0,0.0D0)
	  END IF
	END IF

C	upward treatment of the fluid
	IF (NSS.LT.NSW) THEN
	  TLNWFBNSS = TLNWFANSW
	  YT1WFBNSS = YT1WFANSW
	  YT3WFBNSS = YT3WFANSW
	  SLASK = TLNOPT(NSW+1)
	  CSLASK1 = YT1OPT(NSW+1)
	  CSLASK3 = YT3OPT(NSW+1)
	  TLNOPT(NSW+1) = TLNFWANSW
	  YT1OPT(NSW+1) = YT1FWANSW
	  YT3OPT(NSW+1) = YT3FWANSW
C	  make the call
	  CALL FWOPTYTLIQ ( NSW,NSS+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	  .FALSE., 
     $	  TLNWFBNSS,YT1WFBNSS,YT3WFBNSS, 
     $	  TLNWF, YT1WF,YT3WF, 
     $	  12, TLNOPT,YT1OPT,YT3OPT, NSO,QIOPT )
	  TLNOPT(NSW+1) = SLASK
	  YT1OPT(NSW+1) = CSLASK1
	  YT3OPT(NSW+1) = CSLASK3
	ELSE
	  STOP 64430  ! should not occur
	END IF
	IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
C	  settings for "above the nss level" (see "29/1 98")
	  TLNWFANSS = TLNWFBNSS
	  YT1WFANSS = (0.0D0,0.0D0)
	  YT2WFANSS = (0.0D0,0.0D0)
	  YT3WFANSS = (0.0D0,0.0D0)
	  YT4WFANSS = -YT1WFBNSS
	  YT6WFANSS = YT3WFBNSS
	END IF

C	upward treatment of the upper solid if present
	IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	  TLNWF0 = TLNWFANSS
	  YT1WF0 = YT1WFANSS
	  YT2WF0 = YT2WFANSS
	  YT3WF0 = YT3WFANSS
	  YT4WF0 = YT4WFANSS
	  YT6WF0 = YT6WFANSS
	  IF (NSS.GT.0) THEN
C	    make the call
	    CALL FWOPTYTSOL ( NSS,1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    .FALSE., 
     $	    TLNWF0,YT1WF0,YT2WF0,YT3WF0,YT4WF0,YT6WF0,  
     $	    TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF, 
     $	    2, TLNOPT,YT1OPT,YT2OPT,YT3OPT,YT4OPT,YT6OPT, NSO,QIOPT )
	  END IF
	END IF

 222	WRITE(*,*) ' nso,qiopt: ',NSO,QIOPT
	WRITE(*,*) ' ( nss,nsw,nol are ',NSS,NSW,NOL,' )'

C	compute the solution at depth zd(nso)
	IF ( NSW.GT.0 .AND. NSO.EQ.NSW ) THEN
	  IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	    YS3FW(NSO) = (0.0D0,0.0D0)
	    CSLASK1 = YT4WFBNSW*YT6FWBNSW - YT6WFBNSW*YT4FWBNSW
	    CSLASK2 = -YT3WFBNSW*YT6FWBNSW
	    CSLASK3 = (0.0D0,0.0D0)
	    CSLASK4 = YT1WFBNSW*YT4FWBNSW
	    SLASK1 = DMAX1(DABS(DREAL(CSLASK1)),DABS(DIMAG(CSLASK1)))
	    SLASK2 = DMAX1(DABS(DREAL(CSLASK2)),DABS(DIMAG(CSLASK2)))
	    SLASK3 = 0.0D0
	    SLASK4 = DMAX1(DABS(DREAL(CSLASK4)),DABS(DIMAG(CSLASK4)))
	    ICLOK = 1
	    SLASK = SLASK1
	    IF (SLASK2.GT.SLASK) THEN
	      ICLOK = 2
	      SLASK = SLASK2
	    END IF
	    IF (SLASK4.GT.SLASK) ICLOK = 4
	    GOTO (11,12,990,14) ICLOK
 990	    STOP 3417
 11	    YS1FW(NSO) = CDSQRT(CSLASK1)
	    CSLASK = 1 / YS1FW(NSO)
	    YS2FW(NSO) = CSLASK * (-YT2WFBNSW*YT6FWBNSW)
	    YS4FW(NSO) = CSLASK * (-YT2WFBNSW*YT4FWBNSW)
	    GOTO 999
 12	    YS2FW(NSO) = CDSQRT(CSLASK2)
	    CSLASK = 1 / YS2FW(NSO)
	    YS1FW(NSO) = CSLASK * (-YT2WFBNSW*YT6FWBNSW)
	    YS4FW(NSO) = 0.5D0 * CSLASK * 
     $	    ( YT1WFBNSW*YT6FWBNSW - YT3WFBNSW*YT4FWBNSW )
	    GOTO 999
 14	    YS4FW(NSO) = CDSQRT(CSLASK4)
	    CSLASK = 1 / YS4FW(NSO)
	    YS1FW(NSO) = CSLASK * (-YT2WFBNSW*YT4FWBNSW)
	    YS2FW(NSO) = 0.5D0 * CSLASK * 
     $	    ( YT1WFBNSW*YT6FWBNSW - YT3WFBNSW*YT4FWBNSW )
 999	    SLASK = 0.0D0
	    CALL SCALE3EC ( 0.0D0,0.0D0, SLASK,YS1FW(NSO),YS2FW(NSO),
     $	    YS4FW(NSO) )
	  ELSE
	    CSLASK1 = -YT3WFANSW*YT3FWANSW
	    CSLASK2 = -YT1WFANSW*YT1FWANSW
	    SLASK1 = DMAX1(DABS(DREAL(CSLASK1)),DABS(DIMAG(CSLASK1)))
	    SLASK2 = DMAX1(DABS(DREAL(CSLASK2)),DABS(DIMAG(CSLASK2)))
	    IF (SLASK1.GT.SLASK2) THEN
	      YS2FW(NSO) = CDSQRT(CSLASK1)
	      YS4FW(NSO) = 0.5D0 * ( YT1WFANSW*YT3FWANSW + 
     $	      YT3WFANSW*YT1FWANSW ) / YS2FW(NSO)
	    ELSE
	      YS4FW(NSO) = CDSQRT(CSLASK2)
	      YS2FW(NSO) = 0.5D0 * ( YT1WFANSW*YT3FWANSW + 
     $	      YT3WFANSW*YT1FWANSW ) / YS4FW(NSO)
	    END IF
	    SLASK = 0.0D0
	    CALL SCALE2EC ( 0.0D0,0.0D0, SLASK,YS2FW(NSO),YS4FW(NSO) )
	  END IF
	ELSE IF (NSW.GT.0 .AND. NSO.EQ.NSS ) THEN
	  IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	    YS3FW(NSO) = (0.0D0,0.0D0)
	    CSLASK1 = YT4WFANSS*YT6FWANSS - YT6WFANSS*YT4FWANSS
	    CSLASK2 = YT6WFANSS*YT3FWANSS
	    CSLASK4 = -YT4WFANSS*YT1FWANSS
	    SLASK1 = DMAX1(DABS(DREAL(CSLASK1)),DABS(DIMAG(CSLASK1)))
	    SLASK2 = DMAX1(DABS(DREAL(CSLASK2)),DABS(DIMAG(CSLASK2)))
	    SLASK4 = DMAX1(DABS(DREAL(CSLASK4)),DABS(DIMAG(CSLASK4)))
	    ICLOK = 1
	    SLASK = SLASK1
	    IF (SLASK2.GT.SLASK) THEN
	      ICLOK = 2
	      SLASK = SLASK2
	    END IF
	    IF (SLASK4.GT.SLASK) ICLOK = 4
	    GOTO (21,22,880,24) ICLOK
 880	    STOP 3417
 21	    YS1FW(NSO) = CDSQRT(CSLASK1)
	    CSLASK = 1 / YS1FW(NSO)
	    YS2FW(NSO) = CSLASK * YT6WFANSS*YT2FWANSS
	    YS4FW(NSO) = CSLASK * YT4WFANSS*YT2FWANSS
	    GOTO 888
 22	    YS2FW(NSO) = CDSQRT(CSLASK2)
	    CSLASK = 1 / YS2FW(NSO)
	    YS1FW(NSO) = CSLASK * YT6WFANSS*YT2FWANSS
	    YS4FW(NSO) = 0.5D0 * CSLASK * 
     $	    ( YT4WFANSS*YT3FWANSS - YT6WFANSS*YT1FWANSS )
	    GOTO 888
 24	    YS4FW(NSO) = CDSQRT(CSLASK4)
	    CSLASK = 1 / YS4FW(NSO)
	    YS1FW(NSO) = CSLASK * YT4WFANSS*YT2FWANSS
	    YS2FW(NSO) = 0.5D0 * CSLASK * 
     $	    ( YT4WFANSS*YT3FWANSS - YT6WFANSS*YT1FWANSS )
 888	    SLASK = 0.0D0
	    CALL SCALE3EC ( 0.0D0,0.0D0, SLASK,YS1FW(NSO),YS2FW(NSO),
     $	    YS4FW(NSO) )
	  ELSE
	    CSLASK1 = -YT3WFBNSS*YT3FWBNSS
	    CSLASK2 = -YT1WFBNSS*YT1FWBNSS
	    SLASK1 = DMAX1(DABS(DREAL(CSLASK1)),DABS(DIMAG(CSLASK1)))
	    SLASK2 = DMAX1(DABS(DREAL(CSLASK2)),DABS(DIMAG(CSLASK2)))
	    IF (SLASK1.GT.SLASK2) THEN
	      YS2FW(NSO) = CDSQRT(CSLASK1)
	      YS4FW(NSO) = 0.5D0 * ( YT1WFBNSS*YT3FWBNSS + 
     $	      YT3WFBNSS*YT1FWBNSS ) / YS2FW(NSO)
	    ELSE
	      YS4FW(NSO) = CDSQRT(CSLASK2)
	      YS2FW(NSO) = 0.5D0 * ( YT1WFBNSS*YT3FWBNSS + 
     $	      YT3WFBNSS*YT1FWBNSS ) / YS4FW(NSO)
	    END IF
	    SLASK = 0.0D0
	    CALL SCALE2EC ( 0.0D0,0.0D0, SLASK,YS2FW(NSO),YS4FW(NSO) )
	  END IF
	ELSE IF ( NSS.LT.NSO .AND. NSO.LT.NSW ) THEN
	  YT1FWANSO = YT1FW(NSO+1)
	  YT3FWANSO = YT3FW(NSO+1)
	  IF (IYSTAB.GT.0) THEN
	    IF (IFSMTYP(NSO+1).LE.-4) THEN
	      IF (EXPTYP_FW(NSO+1)) CALL UNDOXLIQ0 (NSO+1, YT1FWANSO,YT3FWANSO)
	    END IF
	  END IF
	  YT1WFBNSO = YT1WF(NSO)
	  YT3WFBNSO = YT3WF(NSO)
	  IF (IYSTAB.GT.0) THEN
	    IF (IFSMTYP(NSO).LE.-4) THEN
	      IF (EXPTYP_FW(NSO)) CALL UNDOXLIQ0 (NSO, YT1WFBNSO,YT3WFBNSO)
	    END IF
	  END IF
	  CSLASK1 = -YT3WFBNSO*YT3FWANSO
	  CSLASK2 = -YT1WFBNSO*YT1FWANSO
	  SLASK1 = DMAX1(DABS(DREAL(CSLASK1)),DABS(DIMAG(CSLASK1)))
	  SLASK2 = DMAX1(DABS(DREAL(CSLASK2)),DABS(DIMAG(CSLASK2)))
	  IF (SLASK1.GT.SLASK2) THEN
	    YS2FW(NSO) = CDSQRT(CSLASK1)
	    YS4FW(NSO) = 0.5D0 * ( YT1WFBNSO*YT3FWANSO + 
     $	    YT3WFBNSO*YT1FWANSO ) / YS2FW(NSO)
	  ELSE
	    YS4FW(NSO) = CDSQRT(CSLASK2)
	    YS2FW(NSO) = 0.5D0 * ( YT1WFBNSO*YT3FWANSO + 
     $	    YT3WFBNSO*YT1FWANSO ) / YS4FW(NSO)
	  END IF
	  SLASK = 0.0D0
	  CALL SCALE2EC ( 0.0D0,0.0D0, SLASK,YS2FW(NSO),YS4FW(NSO) )
	ELSE
C	  we have now either nsw.eq.0 or nso.lt.nss or nsw.lt.nso
	  IF ( NSW.EQ.0 .AND. NSO.EQ.0 ) THEN
C           special 'fix' for this case without 'nss-nsw' fluid
	    YT1FWANSO = YT1FWBNSW
	    YT2FWANSO = YT2FWBNSW
	    YT3FWANSO = YT3FWBNSW
	    YT4FWANSO = YT4FWBNSW
	    YT6FWANSO = YT6FWBNSW
	  ELSE IF ( NSO .LT. (NOL-1) ) THEN
	    YT1FWANSO = YT1FW(NSO+1)
	    YT2FWANSO = YT2FW(NSO+1)
	    YT3FWANSO = YT3FW(NSO+1)
	    YT4FWANSO = YT4FW(NSO+1)
	    YT6FWANSO = YT6FW(NSO+1)
	    IF (IYSTAB.GT.0) THEN
	      IF (IFSMTYP(NSO+1).LE.-4) THEN
	        IF (EXPTYP_FW(NSO+1)) CALL UNDOXSOL0 (U, NSO+1, 
     $	        YT1FWANSO,YT2FWANSO,YT3FWANSO,YT4FWANSO,YT6FWANSO)
	      END IF
	    END IF
	  ELSE
	    YT1FWANSO = YT1FWNOL
	    YT2FWANSO = YT2FWNOL
	    YT3FWANSO = YT3FWNOL
	    YT4FWANSO = YT4FWNOL
	    YT6FWANSO = YT6FWNOL
	  END IF
	  IF (NSO.GT.0) THEN
	    YT1WFBNSO = YT1WF(NSO)
	    YT2WFBNSO = YT2WF(NSO)
	    YT3WFBNSO = YT3WF(NSO)
	    YT4WFBNSO = YT4WF(NSO)
	    YT6WFBNSO = YT6WF(NSO)
	    IF (IYSTAB.GT.0) THEN
	      IF (IFSMTYP(NSO).LE.-4) THEN
	        IF (EXPTYP_FW(NSO)) CALL UNDOXSOL0 (U, NSO, 
     $	        YT1WFBNSO,YT2WFBNSO,YT3WFBNSO,YT4WFBNSO,YT6WFBNSO)
	      END IF
	    END IF
	  ELSE IF (NSW.GT.0) THEN
	    YT1WFBNSO = YT1WF0
	    YT2WFBNSO = YT2WF0
	    YT3WFBNSO = YT3WF0
	    YT4WFBNSO = YT4WF0
	    YT6WFBNSO = YT6WF0
	  ELSE
C	     ytwfbnso already computed by special 'fix' for no 'nss-nsw' fluid
	  END IF
	  CSLASK1 = YT4WFBNSO*YT6FWANSO - YT6WFBNSO*YT4FWANSO
	  CSLASK2 = YT6WFBNSO*YT3FWANSO - YT3WFBNSO*YT6FWANSO
	  CSLASK3 = YT3WFBNSO*YT1FWANSO - YT1WFBNSO*YT3FWANSO
	  CSLASK4 = YT1WFBNSO*YT4FWANSO - YT4WFBNSO*YT1FWANSO
	  SLASK1 = DMAX1(DABS(DREAL(CSLASK1)),DABS(DIMAG(CSLASK1)))
	  SLASK2 = DMAX1(DABS(DREAL(CSLASK2)),DABS(DIMAG(CSLASK2)))
	  SLASK3 = DMAX1(DABS(DREAL(CSLASK3)),DABS(DIMAG(CSLASK3)))
	  SLASK4 = DMAX1(DABS(DREAL(CSLASK4)),DABS(DIMAG(CSLASK4)))
	  ICLOK = 1
	  SLASK = SLASK1
	  IF (SLASK2.GT.SLASK) THEN
	    ICLOK = 2
	    SLASK = SLASK2
	  END IF
	  IF (SLASK3.GT.SLASK) THEN
	    ICLOK = 3
	    SLASK = SLASK3
	  END IF
	  IF (SLASK4.GT.SLASK) ICLOK = 4
	  GOTO (31,32,33,34) ICLOK
 31	  YS1FW(NSO) = CDSQRT(CSLASK1)
	  CSLASK = 1 / YS1FW(NSO)
	  YS2FW(NSO) = CSLASK * ( YT6WFBNSO*YT2FWANSO - YT2WFBNSO*YT6FWANSO )
	  YS3FW(NSO) = 0.5D0 * CSLASK * 
     $	  ( YT1WFBNSO*YT6FWANSO + YT3WFBNSO*YT4FWANSO - 
     $	  YT4WFBNSO*YT3FWANSO - YT6WFBNSO*YT1FWANSO )
	  YS4FW(NSO) = CSLASK * ( YT4WFBNSO*YT2FWANSO - YT2WFBNSO*YT4FWANSO )
	  GOTO 777
 32	  YS2FW(NSO) = CDSQRT(CSLASK2)
	  CSLASK = 1 / YS2FW(NSO)
	  YS1FW(NSO) = CSLASK * ( YT6WFBNSO*YT2FWANSO - YT2WFBNSO*YT6FWANSO )
	  YS3FW(NSO) = CSLASK * ( YT2WFBNSO*YT3FWANSO - YT3WFBNSO*YT2FWANSO )
	  YS4FW(NSO) = 0.5D0 * CSLASK * 
     $	  ( YT1WFBNSO*YT6FWANSO + YT4WFBNSO*YT3FWANSO - 
     $	  YT3WFBNSO*YT4FWANSO - YT6WFBNSO*YT1FWANSO )
	  GOTO 777
 33	  YS3FW(NSO) = CDSQRT(CSLASK3)
	  CSLASK = 1 / YS3FW(NSO)
	  YS1FW(NSO) = 0.5D0 * CSLASK * 
     $	  ( YT1WFBNSO*YT6FWANSO + YT3WFBNSO*YT4FWANSO - 
     $	  YT4WFBNSO*YT3FWANSO - YT6WFBNSO*YT1FWANSO )
	  YS2FW(NSO) = CSLASK * ( YT2WFBNSO*YT3FWANSO - YT3WFBNSO*YT2FWANSO )
	  YS4FW(NSO) = CSLASK * ( YT1WFBNSO*YT2FWANSO - YT2WFBNSO*YT1FWANSO )
	  GOTO 777
 34	  YS4FW(NSO) = CDSQRT(CSLASK4)
	  CSLASK = 1 / YS4FW(NSO)
	  YS1FW(NSO) = CSLASK * ( YT4WFBNSO*YT2FWANSO - YT2WFBNSO*YT4FWANSO )
	  YS2FW(NSO) = 0.5D0 * CSLASK * 
     $	  ( YT1WFBNSO*YT6FWANSO + YT4WFBNSO*YT3FWANSO - 
     $	  YT3WFBNSO*YT4FWANSO - YT6WFBNSO*YT1FWANSO )
	  YS3FW(NSO) = CSLASK * ( YT1WFBNSO*YT2FWANSO - YT2WFBNSO*YT1FWANSO )
 777	  SLASK = 0.0D0
	  CALL SCALE4EC ( 0.0D0,0.0D0, SLASK,YS1FW(NSO),YS2FW(NSO),
     $	  YS3FW(NSO),YS4FW(NSO) )
	END IF

C	compute the rest of the solution by stabilized back-propagation
	IF (NSO.GT.0) THEN
	  IF (NSO.GT.NSS) THEN
	    IF (NSO.GT.NSW) THEN
C	      backward treatment of the lower solid, upward
	      IFSMTYPTMP = IFSMTYP(NSO)
	      IF (IFSMTYPTMP.EQ.-7) THEN
	        IFSMTYP(NSO) = -4
	      ELSE IF (IFSMTYPTMP.EQ.-6) THEN
	        IFSMTYP(NSO) = -5
	      ELSE IF (IFSMTYPTMP.EQ.-3) THEN
	        IFSMTYP(NSO) = 0
	      ELSE IF (IFSMTYPTMP.EQ.-2) THEN
	        IFSMTYP(NSO) = -1
	      END IF
	      CALL SOLSIFW0 ( NSO,NSW+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW )
	      IFSMTYP(NSO) = IFSMTYPTMP     ! reset
	    END IF
	    IF (NSW.GT.NSS) THEN
C	      backward treatment of the "fluid, upper part", upward
	      IF (NSO.GT.NSW) THEN
	        ISLASK = NSW
	      ELSE
	        ISLASK = NSO
	      END IF
	      IFSMTYPTMP = IFSMTYP(ISLASK)
	      IF (IFSMTYPTMP.EQ.-7) THEN
	        IFSMTYP(ISLASK) = -4
	      ELSE IF (IFSMTYPTMP.EQ.-6) THEN
	        IFSMTYP(ISLASK) = -5
	      ELSE IF (IFSMTYPTMP.EQ.-3) THEN
	        IFSMTYP(ISLASK) = 0
	      ELSE IF (IFSMTYPTMP.EQ.-2) THEN
	        IFSMTYP(ISLASK) = -1
	      END IF
	      CALL LIQSIFW0 ( ISLASK,NSS+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNFW, YT1FW,YT3FW )
	      IFSMTYP(ISLASK) = IFSMTYPTMP     ! reset
	      IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	        SLASK1 = DMAX1 ( DABS(DREAL(YT1FWANSS)), 
     $	        DABS(DIMAG(YT1FWANSS)) )
	        SLASK2 = DMAX1 ( DABS(DREAL(YT2FWANSS)), 
     $	        DABS(DIMAG(YT2FWANSS)) )
	        SLASK3 = DMAX1 ( DABS(DREAL(YT3FWANSS)), 
     $	        DABS(DIMAG(YT3FWANSS)) )
	        ICLOK = 1
	        SLASK = SLASK1
	        IF (SLASK2.GT.SLASK) THEN
	          ICLOK = 2
	          SLASK = SLASK2
	        END IF
	        IF (SLASK3.GT.SLASK) ICLOK = 3
	        YS3FW(NSS) = (0.0D0,0.0D0)
	        GOTO (41,42,43) ICLOK
 41	        YS1FW(NSS) =  -YT2FWANSS*YS4FW(NSS) / YT1FWANSS    ! yt5=-yt2
	        GOTO 666
 42	        YS1FW(NSS) = ( -YT4FWANSS*YS2FW(NSS) + 
     $	        YT6FWANSS*YS4FW(NSS) ) / YT2FWANSS
	        GOTO 666
 43	        YS1FW(NSS) =  YT2FWANSS*YS2FW(NSS) / YT3FWANSS     ! yt5=-yt2
 666	        CONTINUE
	      END IF
	    END IF
	  END IF
	  IF (NSW.GT.0) THEN
	    IF (NSS.GT.0) THEN
C	      backward treatment of the upper solid, upward
	      IF (NSO.GT.NSS) THEN
	        ISLASK = NSS
	      ELSE
	        ISLASK = NSO
	      END IF
	      IFSMTYPTMP = IFSMTYP(ISLASK)
	      IF (IFSMTYPTMP.EQ.-7) THEN
	        IFSMTYP(ISLASK) = -4
	      ELSE IF (IFSMTYPTMP.EQ.-6) THEN
	        IFSMTYP(ISLASK) = -5
	      ELSE IF (IFSMTYPTMP.EQ.-3) THEN
	        IFSMTYP(ISLASK) = 0
	      ELSE IF (IFSMTYPTMP.EQ.-2) THEN
	        IFSMTYP(ISLASK) = -1
	      END IF
	      CALL SOLSIFW0 ( ISLASK,1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW )
	      IFSMTYP(ISLASK) = IFSMTYPTMP     ! reset
	    END IF
	  END IF
	END IF
	IF ( NSO .LT. (NOL-1) ) THEN
	  IF (NSO.LT.NSW) THEN
	    IF (NSO.LT.NSS) THEN
C	      backward treatment of the upper solid, downward
	      IFSMTYPTMP = IFSMTYP(NSO+1)
	      IF (IFSMTYPTMP.EQ.-6) THEN
	        IFSMTYP(NSO+1) = -7
	      ELSE IF (IFSMTYPTMP.EQ.-5) THEN
	        IFSMTYP(NSO+1) = -4
	      ELSE IF (IFSMTYPTMP.EQ.-2) THEN
	        IFSMTYP(NSO+1) = -3
	      ELSE IF (IFSMTYPTMP.EQ.-1) THEN
	        IFSMTYP(NSO+1) = 0
	      END IF
	      CALL SOLSIFW0 ( NSO+1,NSS,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF )
	      IFSMTYP(NSO+1) = IFSMTYPTMP     ! reset
	    END IF
	    IF (NSW.GT.NSS) THEN
C	      backward treatment of the "fluid, lower part", downward
	      IF (NSO.LT.NSS) THEN
	        ISLASK = NSS
	      ELSE
	        ISLASK = NSO
	      END IF
	      IFSMTYPTMP = IFSMTYP(ISLASK+1)
	      IF (IFSMTYPTMP.EQ.-6) THEN
	        IFSMTYP(ISLASK+1) = -7
	      ELSE IF (IFSMTYPTMP.EQ.-5) THEN
	        IFSMTYP(ISLASK+1) = -4
	      ELSE IF (IFSMTYPTMP.EQ.-2) THEN
	        IFSMTYP(ISLASK+1) = -3
	      ELSE IF (IFSMTYPTMP.EQ.-1) THEN
	        IFSMTYP(ISLASK+1) = 0
	      END IF
	      CALL LIQSIFW0 ( ISLASK+1,NSW,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNWF, YT1WF,YT3WF )
	      IFSMTYP(ISLASK+1) = IFSMTYPTMP     ! reset
	    END IF
	    IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	      SLASK1 = DMAX1 ( DABS(DREAL(YT1WFBNSW)), 
     $	      DABS(DIMAG(YT1WFBNSW)) )
	      SLASK2 = DMAX1 ( DABS(DREAL(YT2WFBNSW)), 
     $	      DABS(DIMAG(YT2WFBNSW)) )
	      SLASK3 = DMAX1 ( DABS(DREAL(YT3WFBNSW)), 
     $	      DABS(DIMAG(YT3WFBNSW)) )
	      ICLOK = 1
	      SLASK = SLASK1
	      IF (SLASK2.GT.SLASK) THEN
	        ICLOK = 2
	        SLASK = SLASK2
	      END IF
	      IF (SLASK3.GT.SLASK) ICLOK = 3
	      YS3FW(NSW) = (0.0D0,0.0D0)
	      GOTO (51,52,53) ICLOK
 51	      YS1FW(NSW) =  -YT2WFBNSW*YS4FW(NSW) / YT1WFBNSW      ! yt5=-yt2
	      GOTO 555
 52	      YS1FW(NSW) = ( -YT4WFBNSW*YS2FW(NSW) + 
     $	      YT6WFBNSW*YS4FW(NSW) ) / YT2WFBNSW
	      GOTO 555
 53	      YS1FW(NSW) =  YT2WFBNSW*YS2FW(NSW) / YT3WFBNSW       ! yt5=-yt2
 555	      CONTINUE
	    END IF
	  END IF
	  IF ( (NSW+1) .LT. NOL ) THEN
C	    backward treatment of the lower solid, downward
	    IF (NSO.LT.NSW) THEN
	      ISLASK = NSW
	    ELSE
	      ISLASK = NSO
	    END IF
	    IFSMTYPTMP = IFSMTYP(ISLASK+1)
	    IF (IFSMTYPTMP.EQ.-6) THEN
	      IFSMTYP(ISLASK+1) = -7
	    ELSE IF (IFSMTYPTMP.EQ.-5) THEN
	      IFSMTYP(ISLASK+1) = -4
	    ELSE IF (IFSMTYPTMP.EQ.-2) THEN
	      IFSMTYP(ISLASK+1) = -3
	    ELSE IF (IFSMTYPTMP.EQ.-1) THEN
	      IFSMTYP(ISLASK+1) = 0
	    END IF
	    CALL SOLSIFW0 ( ISLASK+1,NOL-1,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF )
	    IFSMTYP(ISLASK+1) = IFSMTYPTMP     ! reset
	  END IF
	END IF

	NSO = NSOORIG                     ! reset

	RETURN
	END
