
C	additional subroutines "dhcompfw etc" for "Full Wave" computations

C	Sven Ivansson  1992






	SUBROUTINE CDETSUBCIRCFW ( VX,VY, RLNDET,CDET )

	USE INMODELa

	IMPLICIT NONE

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

	INCLUDE '../CIRCDATA.INC'  ! up0circ and radiecirc

	REAL*8 VX,VY, PX,PY

	REAL*8 RLNDET
	COMPLEX*16 CDET

	COMPLEX*16 UP, DHFAKT(nhyddim)

	ENTRY CDETSUBCIRC ( VX,VY, RLNDET,CDET )       ! called from rpressint

	UP = UP0CIRC + (RADIECIRC*DEXP(-VY)) * DCMPLX(DCOS(VX),DSIN(VX))
	GOTO 100

	ENTRY CDETSUBFW ( PX,PY, RLNDET,CDET )
	ENTRY CDETSUB ( PX,PY, RLNDET,CDET )           ! called from rpressint
	UP = DCMPLX(PX,PY)

 100	CALL DHCOMPFW (UP, RLNDET,CDET, DHFAKT)

	RETURN
	END






	SUBROUTINE DHCOMPCIRCFW ( VP, RLNDET,CDET, DHFAKT )
C	note that zsolsifw,solsifw,liqsifw use stored values of cp,cq etc 
C	they are stored in the fwsubytzsol,fwsubytzliq calls to hsolpart,
C	liqpart,airypart

	USE INMODELa
	USE FWYSa
	USE FWPARTa
	USE FWDISa
	USE FWYTa
	USE FWZa

	IMPLICIT NONE

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

	INCLUDE '../INPARAM.INC'
	INCLUDE '../FSUBCOM.INC'
	INCLUDE '../INSPEC.INC'
	INCLUDE '../DDULOK.INC'
	INCLUDE '../CIRCDATA.INC'   ! up0circ and radiecirc
	INCLUDE '../SOURCFW.INC'
	INCLUDE 'FWINTF.INC'
	INCLUDE 'FWNSO.INC'

	COMPLEX*16 VP,UPDUMMY

	REAL*8 RLNDET
	COMPLEX*16 CDET, DHFAKT(*)

	INTEGER ID,IDLOK, IH, ICLOK, IFSMTYPTMP, I, IGANG
	REAL*8 RLNDETLOK, ALN,BLN, SLASK, SLASK1,SLASK2,SLASK3, PI
	COMPLEX*16 DUPGDVP, UP, U,RK,UP2,U2,RK2, CDETLOK, 
     $	URA,URB, URANEGI,URBNEGI, G0,G0U,G1,G2,GG,G2R, 
     $	CSLASK, CSLASK1,CSLASK2,CSLASK3,CSLASK4, URVALUE, 
     &	BFAKT,CFAKT,BCFAKT, YY, AIE,AIDE,BIE,BIDE, C60,C120
	LOGICAL DUPGDVPCORR

	SAVE IGANG, RLNDETLOK, CDETLOK

	DATA PI /3.141592653589793D0/ , IGANG /0/

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

	ENTRY DHCOMPCIRC ( VP, RLNDET,CDET, DHFAKT )   ! called from rpressint

	DUPGDVPCORR = .TRUE.
	DUPGDVP = (RADIECIRC*DEXP(-DIMAG(VP))) * 
     $	DCMPLX(-DSIN(DREAL(VP)),DCOS(DREAL(VP)))   ! i*radiecirc*cdexp(i*vp)
	UP = UP0CIRC + DCMPLX(DIMAG(DUPGDVP),-DREAL(DUPGDVP))
	GOTO 100

	ENTRY DHCOMPFW ( UPDUMMY, RLNDET,CDET, DHFAKT )
	ENTRY DHCOMP ( UPDUMMY, RLNDET,CDET, DHFAKT )  ! called from rpressint
	DUPGDVPCORR = .FALSE.
	UP = UPDUMMY

 100	IF (IDHTYP.NE.0) STOP 41156 ! not implemented
	U = U0*UP
	RK = RK0*UP
	UP2 = UP*UP
	U2 = U*U
	RK2 = RK*RK

	UCOP = U
	U2COP = U2
	RK2GOMEG = RK2/OMEG

	IF (IRKSPEC.NE.0) DLNSPEC = ORMINSPEC * DIMAG(RK)

	IF (IGANG.EQ.0) THEN
	  ALLOCATE ( YS1FW(0:noldim-1), YS2FW(0:noldim-1), 
     &	   YS3FW(0:noldim-1), YS4FW(0:noldim-1) )  ! FWYSa
	  ALLOCATE ( EXPTYP_FW(noldim-1),
     &	   SLNREF_FW(noldim-1), E1X_FW(noldim-1),
     &	   CP_FW(noldim-1), CQ_FW(noldim-1), STP_FW(noldim-1),
     &	   STQ_FW(noldim-1), SDP_FW(noldim-1), SDQ_FW(noldim-1),
     &	   TLNREF_FW(noldim-1), E2T_FW(noldim-1),
     &	   URA_FW(noldim-1), URB_FW(noldim-1),
     &	   CPX_FW(noldim-1) )  ! FWPARTa   ! cpx for Airy
	  ALLOCATE ( CFAKT_FW(noldim-1),
     &	   BCFAKT_FW(noldim-1), BCFAKTINV_FW(noldim-1),
     &	   YA_FW(noldim-1), YB_FW(noldim-1) )  ! FWPARTa (NO equivalence)
	  ALLOCATE ( DIS1FW(noldim-1), DIS2FW(noldim-1), 
     &	   DIS3FW(noldim-1), DIS4FW(noldim-1) )  ! FWDISa
	  ALLOCATE ( TLNFW(noldim-1), YT1FW(noldim-1), 
     &	   YT2FW(noldim-1), YT3FW(noldim-1), YT4FW(noldim-1), 
     &	   YT5FW(noldim-1), YT6FW(noldim-1), 
     &	   TLNWF(noldim-1), YT1WF(noldim-1), YT2WF(noldim-1), 
     &	   YT3WF(noldim-1), YT4WF(noldim-1), YT5WF(noldim-1), 
     &	   YT6WF(noldim-1) )  ! FWYTa
	  ALLOCATE ( Z1FW(0:noldim-1), Z2FW(0:noldim-1), 
     &	   Z3FW(0:noldim-1), Z4FW(0:noldim-1) )  ! FWZa
	  ALLOCATE ( Z1WF(0:noldim-1), Z2WF(0:noldim-1), 
     &	   Z3WF(0:noldim-1), Z4WF(0:noldim-1) )  ! FWZa
	  IF (NSOURCFW.GT.0) THEN
C	    set all "dis" to zero
	    CALL ZERODIS
	  END IF
	  IGANG = 1
	END IF

	IF (NSOURCFW.GT.0) THEN
C	  no source at "ordinary nso" now (it is only the "matching depth")
	  DO I = 1, NSOURCFW
	    ID = NSOFW(I)
	    DIS2FW(ID) = RW0FW(I) * RK
	    DIS3FW(ID) = RS0FW(I) * RK2
	  END DO
	END IF

C	forward treatment of the upper solid if present, downward
	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)
	      G1 = G0*U2
	      G2 = 1-G1
	      GG = RHO(0)*G0U
	      G2R = RHO(0)*G2
C	      boundary-equation 2x4 matrix is
C	               -g2r      gg*urbnegi  urbnegi       u
C	            g0u*uranegi     -g2       u/rho   uranegi/rho
	      Z1FWANSS = -G2R*DISF2FW - G0U*URANEGI*DISF1FW
	      Z2FWANSS = GG*URBNEGI*DISF2FW + G2*DISF1FW
	      Z3FWANSS = URBNEGI*DISF2FW - U*DISF1FW/RHO(0)
	      Z4FWANSS = U*DISF2FW - URANEGI*DISF1FW/RHO(0)
	      IF (NSW.EQ.0) THEN
C               special 'fix' for this case without 'nss-nsw' fluid
		GOTO 111
	      END IF
	    ELSE
C	      fluid half-space above
	      TLNFWANSS = 0.0D0
	      YT1FWANSS = (0.0D0,0.0D0)
	      YT2FWANSS = (0.0D0,0.0D0)
	      YT3FWANSS = (0.0D0,0.0D0)
	      YT4FWANSS = RHO(0)
	      YT6FWANSS = URANEGI
	      Z1FWANSS = (0.0D0,0.0D0)
	      Z2FWANSS = RHO(0)*DISF2FW
	      Z3FWANSS = -DISF1FW
	      Z4FWANSS = -URANEGI*DISF2FW
	    END IF
	  ELSE IF (RHO(0).EQ.-99.0D0) THEN
	    STOP '%%% dhcompfw, "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)
	    Z1FWANSS = DISF2FW
	    Z2FWANSS = -DISF1FW
	    Z3FWANSS = (0.0D0,0.0D0)
	    Z4FWANSS = (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)
	    Z1FWANSS = (0.0D0,0.0D0)
	    Z2FWANSS = (0.0D0,0.0D0)
	    Z3FWANSS = DISF2FW
	    Z4FWANSS = -DISF1FW
	  END IF
	  IF (NSS.GT.0) THEN
C	    make the call
	    CALL FWSUBYTZSOL ( 1,NSS,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    .TRUE., DIS1FW,DIS2FW,DIS3FW,DIS4FW, 
     $	    TLNFWANSS, YT1FWANSS,YT2FWANSS,YT3FWANSS,YT4FWANSS,YT6FWANSS, 
     $	    Z1FWANSS,Z2FWANSS,Z3FWANSS,Z4FWANSS, 
     $	    TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, 
     $	    Z1FW,Z2FW,Z3FW,Z4FW )
	  END IF
C	  make appropriate settings for "below the nss level" ("29/1 98")
	  TLNFWBNSS = TLNFWANSS
	  YT1FWBNSS = -YT1FWANSS
	  YT3FWBNSS = -YT3FWANSS
	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNFWBNSS,
     $	  YT1FWBNSS,YT3FWBNSS )
	  Z1FWBNSS = -Z1FWANSS + YT2FWANSS*DIS3NSS
	  IF (TLNFWBNSS.NE.TLNFWANSS) Z1FWBNSS = 
     $	  DEXP(TLNFWANSS-TLNFWBNSS) * Z1FWBNSS
	ELSE
C	  now nss.eq.0 and isolid(0).eq.0, set "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
	      TLNFWANSS = 0.0D0
	      YT1FWANSS = (0.0D0,0.0D0)
	      YT2FWANSS = (0.0D0,0.0D0)
	      YT3FWANSS = (0.0D0,0.0D0)
	      YT4FWANSS = RHO(0)
	      YT6FWANSS = URANEGI
	      Z1FWANSS = (0.0D0,0.0D0)
	      Z2FWANSS = RHO(0)*DISF2FW
	      Z3FWANSS = -DISF1FW
	      Z4FWANSS = -URANEGI*DISF2FW
	      GOTO 111
	    END IF
C	    fluid half-space above fluid
	    TLNFWBNSS = 0.0D0
	    YT1FWBNSS = -RHO(0)
	    YT3FWBNSS = URANEGI
	    Z1FWBNSS = DISF1FW
	  ELSE IF (RHO(0).EQ.-99.0D0) THEN
C	    adaptation from avh/multi/liqsi.f (cf articles/ua15bok/GF)
	    IF (NSW.LT.1) STOP '%%% dhcompfw, "multi"'
	    IF ( IFSMTYP(1).LE.-4 .OR. IFSMTYP(1).GE.1 ) 
     &	     STOP '%%% dhcompfw, "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 '%%% dhcompfw, "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 '%%% dhcompfw, "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 '%%% dhcompfw, "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 '%%% dhcompfw, "multi" bfakt'
	    END IF
	    IF ( DISF1FW .NE. (0.0D0,0.0D0) ) STOP '%%% dhcompfw, "multi"'
	    Z1FWBNSS = (0.0D0,0.0D0)
	  ELSE IF (RHO(0).LT.0.0D0) THEN
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      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)
	      Z1FWANSS = DISF2FW
	      Z2FWANSS = -DISF1FW
	      Z3FWANSS = (0.0D0,0.0D0)
	      Z4FWANSS = (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)
	    Z1FWBNSS = DISF1FW
	  ELSE	      
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      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)
	      Z1FWANSS = (0.0D0,0.0D0)
	      Z2FWANSS = (0.0D0,0.0D0)
	      Z3FWANSS = DISF2FW
	      Z4FWANSS = -DISF1FW
	      GOTO 111
	    END IF
C	    free above fluid
	    TLNFWBNSS = 0.0D0
	    YT1FWBNSS = (0.0D0,0.0D0)
	    YT3FWBNSS = (1.0D0,0.0D0)
	    Z1FWBNSS = DISF1FW
	  END IF
	END IF

C	forward treatment of the "fluid, upper part", downward
	IF (NSO.GT.NSS) THEN
	  TLNFWANSO = TLNFWBNSS
	  YT1FWANSO = YT1FWBNSS
	  YT3FWANSO = YT3FWBNSS
	  Z1FWANSO = Z1FWBNSS
C	  make the call
	  CALL FWSUBYTZLIQ ( NSS+1,NSO,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	  .TRUE., DIS2FW,DIS4FW, 
     $	  TLNFWANSO, YT1FWANSO,YT3FWANSO, Z1FWANSO, 
     $	  TLNFW, YT1FW,YT3FW, Z1FW )
	ELSE
	  TLNFWANSO = TLNFWBNSS
	  YT1FWANSO = YT1FWBNSS
	  YT3FWANSO = YT3FWBNSS
	  Z1FWANSO = Z1FWBNSS
	END IF

C	forward treatment of the lower solid if present, upward
 111	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)
	      G1 = G0*U2
	      G2 = 1-G1
	      GG = RHO(NOL)*G0U
	      G2R = RHO(NOL)*G2
C	      boundary-equation 2x4 matrix is
C	               -g2r      gg*urbnegi  urbnegi       u
C	            g0u*uranegi     -g2       u/rho   uranegi/rho
	      Z1WFBNSW = -G2R*DISG2FW - G0U*URANEGI*DISG1FW
	      Z2WFBNSW = GG*URBNEGI*DISG2FW + G2*DISG1FW
	      Z3WFBNSW = URBNEGI*DISG2FW - U*DISG1FW/RHO(NOL)
	      Z4WFBNSW = U*DISG2FW - URANEGI*DISG1FW/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
	      Z1WFBNSW = (0.0D0,0.0D0)
	      Z2WFBNSW = RHO(NOL)*DISG2FW
	      Z3WFBNSW = -DISG1FW
	      Z4WFBNSW = -URANEGI*DISG2FW
	    END IF
	  ELSE IF (RHO(NOL).EQ.-99.0D0) THEN
	    STOP '%%% dhcompfw, 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)
	    Z1WFBNSW = DISG2FW
	    Z2WFBNSW = -DISG1FW
	    Z3WFBNSW = (0.0D0,0.0D0)
	    Z4WFBNSW = (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)
	    Z1WFBNSW = (0.0D0,0.0D0)
	    Z2WFBNSW = (0.0D0,0.0D0)
	    Z3WFBNSW = DISG2FW
	    Z4WFBNSW = -DISG1FW
	  END IF
	  IF ( (NSW+1) .LT. NOL ) THEN
C	   make the call
	   CALL FWSUBYTZSOL ( NOL-1,NSW+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	   .TRUE., DIS1FW,DIS2FW,DIS3FW,DIS4FW, 
     $	   TLNWFBNSW, YT1WFBNSW,YT2WFBNSW,YT3WFBNSW,YT4WFBNSW,
     $	   YT6WFBNSW, Z1WFBNSW,Z2WFBNSW,Z3WFBNSW,Z4WFBNSW, 
     $	   TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF, 
     $	   Z1WF,Z2WF,Z3WF,Z4WF )
	  END IF
	  IF (NSW.EQ.0) THEN
C           special 'fix' for this case without 'nss-nsw' fluid
	    GOTO 222
	  END IF
C	  make appropriate settings for "above the nsw level" ("29/1 98")
	  TLNWFANSW = TLNWFBNSW
	  YT1WFANSW = YT1WFBNSW
	  YT3WFANSW = YT3WFBNSW
	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNWFANSW,
     $	  YT1WFANSW,YT3WFANSW )
	  Z1WFANSW = Z1WFBNSW - YT2WFBNSW*DIS3NSW
	  IF (TLNWFANSW.NE.TLNWFBNSW) Z1WFANSW = 
     $	  DEXP(TLNWFBNSW-TLNWFANSW) * Z1WFANSW
	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
	    Z1WFANSW = DISG1FW
	  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)
	    IF ( DISG1FW .NE. (0.0D0,0.0D0) ) STOP '%dhcompfw, loc-react prbl'
	    Z1WFANSW = (0.0D0,0.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)
	    Z1WFANSW = DISG1FW
	  ELSE	      
C	    free below fluid
	    TLNWFANSW = 0.0D0
	    YT1WFANSW = (0.0D0,0.0D0)
	    YT3WFANSW = (-1.0D0,0.0D0)
	    Z1WFANSW = DISG1FW
	  END IF
	END IF

C	forward treatment of the "fluid, lower part", upward
	IF (NSO.LT.NSW) THEN
	  TLNWFBNSO = TLNWFANSW
	  YT1WFBNSO = YT1WFANSW
	  YT3WFBNSO = YT3WFANSW
	  Z1WFBNSO = Z1WFANSW
C	  make the call
	  CALL FWSUBYTZLIQ ( NSW,NSO+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	  .TRUE., DIS2FW,DIS4FW, 
     $	  TLNWFBNSO, YT1WFBNSO,YT3WFBNSO, Z1WFBNSO, 
     $	  TLNWF, YT1WF,YT3WF, Z1WF )
	ELSE
	  TLNWFBNSO = TLNWFANSW
	  YT1WFBNSO = YT1WFANSW
	  YT3WFBNSO = YT3WFANSW
	  Z1WFBNSO = Z1WFANSW
	END IF

C	compute the solution at depth "zd(nso)" by matching at that depth
 222	IF (NSW.GT.0) THEN
C	  first - for cdetlok to fit with liqsi, (yp1,yp2) <-> +-(-yt3,yt1)
	  RLNDETLOK = TLNFWANSO + TLNWFBNSO
	  CDETLOK = - ( YT1FWANSO*YT3WFBNSO - YT1WFBNSO*YT3FWANSO )
	  RLNDET = RLNDETLOK
	  CDET = CDETLOK
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, RLNDET,CDET )
	  IF (IDET.EQ.1) RETURN
	  IF (RW0.NE.(0.0D0,0.0D0)) THEN
	    CSLASK2 = Z1WFBNSO*YT3FWANSO - Z1FWANSO*YT3WFBNSO
	    CSLASK4 = Z1FWANSO*YT1WFBNSO - Z1WFBNSO*YT1FWANSO
	    IF (IDET.NE.-1) THEN
	      YS2FW(NSO) = CSLASK2 / CDETLOK
	      YS4FW(NSO) = CSLASK4 / CDETLOK
	    ELSE
	      SLASK = 2*PI * DEXP(RLNDETLOK-DLNDULOK)
	      CSLASK2 = SLASK * CSLASK2 / DERIVDULOK
	      YS2FW(NSO) = DCMPLX(-DIMAG(CSLASK2),DREAL(CSLASK2))   ! 2*pi*i * 
	      CSLASK4 = SLASK * CSLASK4 / DERIVDULOK
	      YS4FW(NSO) = DCMPLX(-DIMAG(CSLASK4),DREAL(CSLASK4))   ! 2*pi*i * 
C	      for residue contribution (positive orientation)
	      SLASK = DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	      ( DMAX1 ( DABS(DREAL(YT1FWANSO)),DABS(DIMAG(YT1FWANSO)),
     $	      DABS(DREAL(YT3FWANSO)),DABS(DIMAG(YT3FWANSO)) ) * 
     $	      DMAX1 ( DABS(DREAL(YT1WFBNSO)),DABS(DIMAG(YT1WFBNSO)),
     $	      DABS(DREAL(YT3WFBNSO)),DABS(DIMAG(YT3WFBNSO)) ) )
	      WRITE(*,*) ' indicators: ',SLASK,
     $	      DEXP(RLNDETLOK-DLNDULOK) * 
     $	      DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	      DMAX1(DABS(DREAL(DERIVDULOK)),DABS(DIMAG(DERIVDULOK)))
	    END IF
	  ELSE
C	    choose either of (yt3fwanso,-yt1fwanso) and (yt3wfbnso,-yt1wfbnso)
	    YS2FW(NSO) = YT3WFBNSO
	    YS4FW(NSO) = -YT1WFBNSO
	  END IF
	ELSE
C         special 'fix' for the case without 'nss-nsw' fluid
	  IF (NSO.NE.0) STOP 34445
	  RLNDETLOK = TLNFWANSS + TLNWFBNSW
	  CDETLOK = YT1FWANSS*YT6WFBNSW + 2*YT2FWANSS*YT2WFBNSW + 
     $    YT3FWANSS*YT4WFBNSW + YT4FWANSS*YT3WFBNSW + YT6FWANSS*YT1WFBNSW
	  RLNDET = RLNDETLOK
	  CDET = CDETLOK
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, RLNDET,CDET )
	  IF (IDET.EQ.1) RETURN
	  IF (RW0.NE.(0.0D0,0.0D0)) THEN
C           concatenation
	    CSLASK1 = YT1FWANSS*Z3WFBNSW - YT2FWANSS*Z2WFBNSW + 
     $      Z1FWANSS*YT4WFBNSW + YT4FWANSS*Z1WFBNSW - Z2FWANSS*YT2WFBNSW + 
     $      Z3FWANSS*YT1WFBNSW
	    CSLASK2 = YT1FWANSS*Z4WFBNSW - YT3FWANSS*Z2WFBNSW - 
     $      Z1FWANSS*YT2WFBNSW - YT2FWANSS*Z1WFBNSW - Z2FWANSS*YT3WFBNSW + 
     $      Z4FWANSS*YT1WFBNSW
	    CSLASK3 = YT2FWANSS*Z4WFBNSW - YT3FWANSS*Z3WFBNSW + 
     $      Z1FWANSS*YT6WFBNSW + YT6FWANSS*Z1WFBNSW - Z3FWANSS*YT3WFBNSW + 
     $      Z4FWANSS*YT2WFBNSW
	    CSLASK4 = YT4FWANSS*Z4WFBNSW + YT2FWANSS*Z3WFBNSW + 
     $      Z2FWANSS*YT6WFBNSW + YT6FWANSS*Z2WFBNSW + Z3FWANSS*YT2WFBNSW + 
     $      Z4FWANSS*YT4WFBNSW
	    IF (IDET.NE.-1) THEN
	      YS1FW(NSO) = -CSLASK4 / CDETLOK
	      YS2FW(NSO) = +CSLASK3 / CDETLOK
	      YS3FW(NSO) = -CSLASK2 / CDETLOK
	      YS4FW(NSO) = +CSLASK1 / CDETLOK
	    ELSE
	      SLASK = 2*PI * DEXP(RLNDETLOK-DLNDULOK)
	      CSLASK4 = SLASK * CSLASK4 / DERIVDULOK
	      YS1FW(NSO) = -DCMPLX(-DIMAG(CSLASK4),DREAL(CSLASK4))  ! 2*pi*i * 
	      CSLASK3 = SLASK * CSLASK3 / DERIVDULOK
	      YS2FW(NSO) = +DCMPLX(-DIMAG(CSLASK3),DREAL(CSLASK3))  ! 2*pi*i * 
	      CSLASK2 = SLASK * CSLASK2 / DERIVDULOK
	      YS3FW(NSO) = -DCMPLX(-DIMAG(CSLASK2),DREAL(CSLASK2))  ! 2*pi*i * 
	      CSLASK1 = SLASK * CSLASK1 / DERIVDULOK
	      YS4FW(NSO) = +DCMPLX(-DIMAG(CSLASK1),DREAL(CSLASK1))  ! 2*pi*i * 
C	      for residue contribution (positive orientation)
	      SLASK = DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	      ( DMAX1 ( DABS(DREAL(YT1FWANSS)),DABS(DIMAG(YT1FWANSS)),
     $	      DABS(DREAL(YT2FWANSS)),DABS(DIMAG(YT2FWANSS)),
     $	      DABS(DREAL(YT3FWANSS)),DABS(DIMAG(YT3FWANSS)),
     $	      DABS(DREAL(YT4FWANSS)),DABS(DIMAG(YT4FWANSS)),
     $	      DABS(DREAL(YT6FWANSS)),DABS(DIMAG(YT6FWANSS)) ) * 
     $	      DMAX1 ( DABS(DREAL(YT1WFBNSW)),DABS(DIMAG(YT1WFBNSW)),
     $	      DABS(DREAL(YT2WFBNSW)),DABS(DIMAG(YT2WFBNSW)),
     $	      DABS(DREAL(YT3WFBNSW)),DABS(DIMAG(YT3WFBNSW)),
     $	      DABS(DREAL(YT4WFBNSW)),DABS(DIMAG(YT4WFBNSW)),
     $	      DABS(DREAL(YT6WFBNSW)),DABS(DIMAG(YT6WFBNSW)) ) )
	      WRITE(*,*) ' indicators: ',SLASK,
     $	      DEXP(RLNDETLOK-DLNDULOK) * 
     $	      DMAX1(DABS(DREAL(CDETLOK)),DABS(DIMAG(CDETLOK))) / 
     $	      DMAX1(DABS(DREAL(DERIVDULOK)),DABS(DIMAG(DERIVDULOK)))
	    END IF
	  ELSE
	    STOP 34447    ! should not occur
	  END IF
	END IF
	IF (NHYD.LE.0) RETURN

C	compute the rest of the solution by stabilized back-propagation
	IDLOK = MAX ( NSS, IDHYD(1) )
	IF (NSO.GT.IDLOK) THEN
C	  backward treatment of the "fluid, upper part", 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
	  IF (IDET.NE.-1) THEN
	    CALL LIQSIFW ( NSO,IDLOK+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    TLNFW, YT1FW,YT3FW, Z1FW, DIS2FW,DIS4FW )
	  ELSE
	    CALL LIQSIFW0 ( NSO,IDLOK+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    TLNFW, YT1FW,YT3FW )
	  END IF
	  IFSMTYP(NSO) = IFSMTYPTMP     ! reset
	END IF
	IF ( NSW.GT.0 .AND. NSS.GE.IDHYD(1) ) THEN
	  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
	    IF (IDET.NE.-1) THEN
	      YS3FW(NSS) = DIS3NSS
	    ELSE
	      YS3FW(NSS) = (0.0D0,0.0D0)
	    END IF
	    GOTO (1,2,3) ICLOK
 1	    IF (IDET.NE.-1) THEN
	      YS1FW(NSS) = ( YT4FWANSS*DIS3NSS - 
     $	      YT2FWANSS*YS4FW(NSS) - Z2FWANSS ) / YT1FWANSS         ! yt5=-yt2
	    ELSE
	      YS1FW(NSS) =  -YT2FWANSS*YS4FW(NSS) / YT1FWANSS       ! yt5=-yt2
	    END IF
	    GOTO 9
 2	    IF (IDET.NE.-1) THEN
	      YS1FW(NSS) = ( -YT4FWANSS*YS2FW(NSS) + 
     $	      YT6FWANSS*YS4FW(NSS) - Z3FWANSS ) / YT2FWANSS
	    ELSE
	      YS1FW(NSS) = ( -YT4FWANSS*YS2FW(NSS) + 
     $	      YT6FWANSS*YS4FW(NSS) ) / YT2FWANSS
	    END IF
	    GOTO 9
 3	    IF (IDET.NE.-1) THEN
	      YS1FW(NSS) = ( YT2FWANSS*YS2FW(NSS) - 
     $	      YT6FWANSS*DIS3NSS - Z4FWANSS ) / YT3FWANSS            ! yt5=-yt2
	    ELSE
	      YS1FW(NSS) =  YT2FWANSS*YS2FW(NSS) / YT3FWANSS        ! yt5=-yt2
	    END IF
	  END IF
 9	  IF (NSS.GT.IDHYD(1)) THEN
C	    backward treatment of the upper solid, upward
	    IDLOK = IDHYD(1)
	    IF (IDET.NE.-1) THEN
	      CALL SOLSIFW ( NSS,IDLOK+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW, 
     $	      Z1FW,Z2FW,Z3FW,Z4FW, DIS1FW,DIS2FW,DIS3FW,DIS4FW )
	    ELSE
	      CALL SOLSIFW0 ( NSS,IDLOK+1,-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNFW, YT1FW,YT2FW,YT3FW,YT4FW,YT5FW,YT6FW )
	    END IF
	  END IF
	END IF
	IDLOK = MIN ( NSW, IDHYD(NHYD) )
	IF (NSO.LT.IDLOK) THEN
C	  backward treatment of the "fluid, lower part", 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
	  IF (IDET.NE.-1) THEN
	    CALL LIQSIFW ( NSO+1,IDLOK,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    TLNWF, YT1WF,YT3WF, Z1WF, DIS2FW,DIS4FW )
	  ELSE
	    CALL LIQSIFW0 ( NSO+1,IDLOK,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    TLNWF, YT1WF,YT3WF )
	  END IF
	  IFSMTYP(NSO+1) = IFSMTYPTMP   ! reset
	END IF
	IF (NSW.LE.IDHYD(NHYD)) THEN
	  IF ( (NSW+1).LT.NOL .OR. ISOLID(NOL).NE.0 ) THEN
	    IF (NSW.GT.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
	      IF (IDET.NE.-1) THEN
	        YS3FW(NSW) = DIS3NSW
	      ELSE
	        YS3FW(NSW) = (0.0D0,0.0D0)
	      END IF
	      GOTO (11,22,33) ICLOK
 11	      IF (IDET.NE.-1) THEN
	        YS1FW(NSW) = ( YT4WFBNSW*DIS3NSW - 
     $	        YT2WFBNSW*YS4FW(NSW) - Z2WFBNSW ) / YT1WFBNSW       ! yt5=-yt2
	      ELSE
	        YS1FW(NSW) =  -YT2WFBNSW*YS4FW(NSW) / YT1WFBNSW     ! yt5=-yt2
	      END IF
	      GOTO 99
 22	      IF (IDET.NE.-1) THEN
	        YS1FW(NSW) = ( -YT4WFBNSW*YS2FW(NSW) + 
     $	        YT6WFBNSW*YS4FW(NSW) - Z3WFBNSW ) / YT2WFBNSW
	      ELSE
	        YS1FW(NSW) = ( -YT4WFBNSW*YS2FW(NSW) + 
     $	        YT6WFBNSW*YS4FW(NSW) ) / YT2WFBNSW
	      END IF
	      GOTO 99
 33	      IF (IDET.NE.-1) THEN
	        YS1FW(NSW) = ( YT2WFBNSW*YS2FW(NSW) - 
     $	        YT6WFBNSW*DIS3NSW - Z4WFBNSW ) / YT3WFBNSW          ! yt5=-yt2
	      ELSE
	        YS1FW(NSW) =  YT2WFBNSW*YS2FW(NSW) / YT3WFBNSW      ! yt5=-yt2
	      END IF
	    ELSE
C             ysfw(nsw) already comptd by special 'fix' for no 'nss-nsw' fluid
	    END IF
	  END IF
 99	  IF (NSW.LT.IDHYD(NHYD)) THEN
C	    backward treatment of the lower solid, downward
	    IDLOK = IDHYD(NHYD)
	    IF (IDET.NE.-1) THEN
	      CALL SOLSIFW ( NSW+1,IDLOK,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF, 
     $	      Z1WF,Z2WF,Z3WF,Z4WF, DIS1FW,DIS2FW,DIS3FW,DIS4FW )
	    ELSE
	      CALL SOLSIFW0 ( NSW+1,IDLOK,1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      TLNWF, YT1WF,YT2WF,YT3WF,YT4WF,YT5WF,YT6WF )
	    END IF
	  END IF
	END IF

	IF (NSOURCFW.GT.0) THEN
C	  set the "dhfakt" array
	  IF (KFWTYP.EQ.1) THEN
C           u (sign-changed)
	    CSLASK = 1/OMEG
	    DO IH = 1, NHYD
	      ID = IDHYD(IH)
	      IF ( ID.LE.NSS .AND. (NSS.GT.0.OR.ISOLID(0).NE.0) ) THEN
	        DHFAKT(IH) = CSLASK * YS1FW(ID)         ! r1 = ys1/omeg
	      ELSE IF ( ID.GE.NSW .AND. 
     $        ((NSW+1).LT.NOL.OR.ISOLID(NOL).NE.0) ) THEN
	        DHFAKT(IH) = CSLASK * YS1FW(ID)         ! r1 = ys1/omeg
	      ELSE
	        DHFAKT(IH) = (0.0D0,0.0D0)  ! ys1fw not computed
	      END IF
	    END DO
	  ELSE IF (KFWTYP.EQ.2) THEN
C           w (sign-changed)
	    CSLASK = 1/OMEG
	    DO IH = 1, NHYD
	      DHFAKT(IH) = CSLASK * YS2FW(IDHYD(IH))    ! r2 = ys2/omeg
	    END DO
	  ELSE IF (KFWTYP.EQ.3) THEN
C           tau_rz (sign-changed)
	    DO IH = 1, NHYD
	      ID = IDHYD(IH)
	      IF ( ID.LE.NSS .AND. (NSS.GT.0.OR.ISOLID(0).NE.0) ) THEN
	        DHFAKT(IH) = YS3FW(ID)                  ! r3 = ys3
	      ELSE IF ( ID.GE.NSW .AND. 
     $        ((NSW+1).LT.NOL.OR.ISOLID(NOL).NE.0) ) THEN
	        DHFAKT(IH) = YS3FW(ID)                  ! r3 = ys3
	      ELSE
	        DHFAKT(IH) = (0.0D0,0.0D0)  ! r3 zero in fluid
	      END IF
	    END DO
	  ELSE IF (KFWTYP.EQ.4) THEN
C           tau_zz (sign-changed)
	    DO IH = 1, NHYD
	      DHFAKT(IH) = YS4FW(IDHYD(IH))             ! r4 = ys4
	    END DO
	  ELSE
	    STOP 4127761
	  END IF
	  IF (DUPGDVPCORR) THEN
	    DO IH = 1, NHYD
	      DHFAKT(IH) = DUPGDVP * DHFAKT(IH)
	    END DO
	  END IF
C	  handling of "irkspec" (obs)
	  IF (IRKSPEC.NE.0) THEN
	    SLASK = DEXP(-DLNSPEC)
	    DO IH = 1, NHYD
	      DHFAKT(IH) = SLASK * DHFAKT(IH)
	    END DO
	  END IF
	ELSE
	  STOP 477700
	END IF

	RETURN
	END

