
C	subroutine "grderiv" to compute dispersion-function value and 
C	derivatives
C	note: derivatives with respect to "omeg & u" at "up"
C	      ( "omeg & up" would not be convenient when "u0" in "u=u0*up" is 
C	        frequency-dependent )

C	Sven Ivansson  1993



	SUBROUTINE GRDERIVDUSUB ( PX,PY, RLNDET,CDET, 
     $	DLNDU,DERIVDU )
	IMPLICIT NONE
	REAL*16 PX,PY
	REAL*16 RLNDET, DLNDU
	COMPLEX*32 CDET, DERIVDU
	COMPLEX*32 UP, QCMPLX
	UP = QCMPLX(PX,PY)
	CALL GRDERIVDU ( UP, RLNDET,CDET, DLNDU,DERIVDU )
	RETURN
	END

	SUBROUTINE GRVALUESUB ( PX,PY, RLNDET,CDET )
	IMPLICIT NONE
	REAL*16 PX,PY
	REAL*16 RLNDET
	COMPLEX*32 CDET
	COMPLEX*32 UP, QCMPLX
	UP = QCMPLX(PX,PY)
	CALL GRVALUE ( UP, RLNDET,CDET )
	RETURN
	END





	SUBROUTINE GRDERIV ( UP, RLNDET,CDET, 
     $	DLNDO,DERIVDO, DLNDU,DERIVDU )

	USE INMODELa
	USE GRLOCALa

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       nol <= noldim , nhyd <= nhyddim

	INCLUDE 'INPARAM.INC'
	INCLUDE 'FSUBCOM.INC'

	COMPLEX*32 UP

	REAL*16 RLNDET, DLNDO, DLNDU
	COMPLEX*32 CDET, DERIVDO, DERIVDU

	INTEGER IHARK0
	REAL*16 DLNDULOK, TLN, PLN, SLASK
	COMPLEX*32 DERIVDULOK, YT1,YT2,YT3,YT4,YT6, 
     $  YP1,YP2, U,RK,UP2,U2,RK2, URANOL,URBNOL, URA0,URB0,
     $  G0U, YT10PRE,YT20,YT30,YT40,YT60, CSLASK,CSLASK2, URVALUE, QCMPLX
	LOGICAL DIMDONE

	SAVE U,RK,UP2,U2,RK2, URANOL,URBNOL, URA0,URB0,
     &	 G0U, YT10PRE,YT20,YT30,YT40,YT60, DIMDONE
	DATA DIMDONE/.FALSE./

	IGRENTRY = 12
	GOTO 100

	ENTRY GRVALUE ( UP, RLNDET,CDET )
	IGRENTRY = 0
	GOTO 100

	ENTRY GRDERIVDO ( UP, RLNDET,CDET, DLNDO,DERIVDO )
	IGRENTRY = 1
	GOTO 100

	ENTRY GRDERIVDU ( UP, RLNDET,CDET, DLNDU,DERIVDU )
	IGRENTRY = 2

 100	IF (.NOT.DIMDONE) THEN
	  ALLOCATE ( TLNB(0:noldim), 
     &	  YTB1(0:noldim-1),YTB2(0:noldim-1),YTB3(0:noldim-1),
     &	  YTB4(0:noldim-1),YTB6(0:noldim-1), 
     &	  YTB1K(1:noldim-1),YTB2K(1:noldim-1),YTB3K(1:noldim-1),
     &	  YTB4K(1:noldim-1),YTB6K(1:noldim-1), 
     &	  YTB1L(1:noldim-1),YTB2L(1:noldim-1),YTB3L(1:noldim-1),
     &	  YTB4L(1:noldim-1),YTB6L(1:noldim-1), 
     &	  YTB1M(1:noldim-1),YTB2M(1:noldim-1),YTB3M(1:noldim-1),
     &	  YTB4M(1:noldim-1),YTB6M(1:noldim-1), 
     &	  YTB1N(1:noldim),YTB2N(1:noldim),YTB3N(1:noldim),
     &	  YTB4N(1:noldim),YTB6N(1:noldim), 
     &	  EXPTYP_GR(noldim-1), TLNREF_GR(noldim-1), E2T_GR(noldim-1), 
     &	  CP_GR(noldim-1), CQ_GR(noldim-1), 
     &	  STP_GR(noldim-1), STQ_GR(noldim-1), 
     &	  SDP_GR(noldim-1), SDQ_GR(noldim-1), 
     &	  URA_GR(noldim-1), URB_GR(noldim-1), 
     &	  CPX_GR(noldim-1) )  ! GRLOCALa           ! cpx for Airy
	ALLOCATE ( CFAKT_GR(noldim-1), 
     &	 BCFAKT_GR(noldim-1), BCFAKTINV_GR(noldim-1), 
     &	 YA_GR(noldim-1), YB_GR(noldim-1) )  ! GRLOCALa
ccc	EQUIVALENCE (CQ_GR,CFAKT_GR), 
ccc     &	(STQ_GR,BCFAKT_GR), (SDQ_GR,BCFAKTINV_GR), 
ccc     &	(URA_GR,YA_GR), (URB_GR,YB_GR)  ! GRLOCALa
	  DIMDONE = .TRUE.
	END IF

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

	UCOP = U
	U2COP = U2
	RK2GOMEG = RK2/OMEG

	IF (RHO(NOL).GT.0.0Q0) THEN
C	  half-space present, "branch-cut(s) and sheet(s)" important
	  URANOL = URVALUE (ICUTTYP,AIMCUT, U0,UPBPA,UPBPA2,UPREF, 
     $	  UP,UP2)
	  IF ( ISHEETTYP.EQ.2 .OR. ISHEETTYP.EQ.4 ) URANOL = -URANOL
	  IF (ISOLID(NOL).NE.0) THEN
C	    solid half-space below
	    URBNOL = URVALUE (ICUTTYP,AIMCUT, U0,UPBPB,UPBPB2,UPREF, 
     $	    UP,UP2)
	    IF ( ISHEETTYP.EQ.3 .OR. ISHEETTYP.EQ.4 ) URBNOL = -URBNOL
	    IHARK0 = 1
	    TLNB(NOL) = 0.0Q0
	    YTB1N(NOL) = (1.0Q0,0.0Q0)
	    YTB2N(NOL) = (0.0Q0,0.0Q0)
	    YTB3N(NOL) = QCMPLX(-AIMAG(URANOL),REAL(URANOL))    !  i*ura
	    YTB4N(NOL) = QCMPLX(AIMAG(URBNOL),-REAL(URBNOL))    ! -i*urb
	    YTB6N(NOL) = -URANOL*URBNOL
	    CALL HARKDERIVB ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0 )
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      GOTO 222
	    END IF
	    CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNB(NSW),YTB1(NSW),
     $	    YTB3(NSW) )   ! ytb2(nsw),ytb4(nsw),ytb6(nsw) now irrelevant
	  ELSE
C	    fluid half-space below
	    IF ( NOL .GT. (NSW+1) ) THEN
	      IHARK0 = 0
	      TLNB(NOL-1) = 0.0Q0
	      YTB1(NOL-1) = (0.0Q0,0.0Q0)
	      YTB2(NOL-1) = (0.0Q0,0.0Q0)
	      YTB3(NOL-1) = (0.0Q0,0.0Q0)
	      YTB4(NOL-1) = RHO(NOL)
	      YTB6(NOL-1) = QCMPLX(AIMAG(URANOL),-REAL(URANOL))  ! -i*ura
	      CALL HARKDERIVB ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0 )
	      IF (NSW.EQ.0) THEN
C               special 'fix' for this case without 'nss-nsw' fluid
	        GOTO 222
	      END IF
	      CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNB(NSW),YTB1(NSW),
     $	      YTB3(NSW) )   ! ytb2(nsw),ytb4(nsw),ytb6(nsw) now irrelevant
	    ELSE
	      TLNB(NSW) = 0.0Q0
	      YTB1(NSW) = RHO(NOL)
	      YTB3(NSW) = -QCMPLX(AIMAG(URANOL),-REAL(URANOL))	  ! -(-i*ura)
	    END IF
	  END IF
	ELSE IF (RHO(NOL).LT.0.0Q0) THEN
C	  rigid below
	  IF ( NOL .GT. (NSW+1) ) THEN
	    IHARK0 = -1
	    TLNB(NOL-1) = 0.0Q0
	    YTB1(NOL-1) = (1.0Q0,0.0Q0)
	    YTB2(NOL-1) = (0.0Q0,0.0Q0)
	    YTB3(NOL-1) = (0.0Q0,0.0Q0)
	    YTB4(NOL-1) = (0.0Q0,0.0Q0)
	    YTB6(NOL-1) = (0.0Q0,0.0Q0)
	    CALL HARKDERIVB ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0 )
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      GOTO 222
	    END IF
	    CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNB(NSW),YTB1(NSW),
     $	    YTB3(NSW) )   ! ytb2(nsw),ytb4(nsw),ytb6(nsw) now irrelevant
	  ELSE
	    TLNB(NSW) = 0.0Q0
	    YTB1(NSW) = (1.0Q0,0.0Q0)
	    YTB3(NSW) = (0.0Q0,0.0Q0)
	  END IF
	ELSE
C	  free below
	  IF ( NOL .GT. (NSW+1) ) THEN
	    IHARK0 = -2
	    TLNB(NOL-1) = 0.0Q0
	    YTB1(NOL-1) = (0.0Q0,0.0Q0)
	    YTB2(NOL-1) = (0.0Q0,0.0Q0)
	    YTB3(NOL-1) = (0.0Q0,0.0Q0)
	    YTB4(NOL-1) = (0.0Q0,0.0Q0)
	    YTB6(NOL-1) = (1.0Q0,0.0Q0)
	    CALL HARKDERIVB ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0 )
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      GOTO 222
	    END IF
	    CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNB(NSW),YTB1(NSW),
     $	    YTB3(NSW) )   ! ytb2(nsw),ytb4(nsw),ytb6(nsw) now irrelevant
	  ELSE
	    TLNB(NSW) = 0.0Q0
	    YTB1(NSW) = (0.0Q0,0.0Q0)
	    YTB3(NSW) = (-1.0Q0,0.0Q0)
	  END IF
	END IF

	IF (NSW.GT.NSS) THEN
	  CALL LIQDERIVB ( NSW,NSS+1, OMEG,OMEG2,U,U2,RK,RK2 )
	ELSE
	  STOP 64430  ! should not occur
	END IF

	IF ( NSS.GT.0 .OR. ISOLID(0).NE.0 ) THEN
	  YTB4(NSS) = -YTB1(NSS)      ! before 'zero:ed'
	  YTB6(NSS) = YTB3(NSS)       ! before 'zero:ed'
	  YTB1(NSS) = (0.0Q0,0.0Q0)
	  YTB2(NSS) = (0.0Q0,0.0Q0)
	  YTB3(NSS) = (0.0Q0,0.0Q0)
	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLNB(NSS),YTB4(NSS),
     $    YTB6(NSS) )
	  IF (NSS.GT.0) THEN
	    IHARK0 = 0
	    CALL HARKDERIVB ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0 )
	  END IF
C	  free or rigid or half-space above
	END IF

 222	RLNDET = TLNB(0)
	IF (RHO(0).GT.0.0Q0) THEN
C	  half-space present, "branch-cut(s) and sheet(s)" important
	  URA0 = URVALUE (ICUTTYP,AIMCUT, U0,UPBPASURF,UPBPA2SURF,UPREF, 
     $	  UP,UP2)
	  IF ( ISHEETTYP.EQ.2 .OR. ISHEETTYP.EQ.4 ) URA0 = -URA0
	  URA0 = -URA0		! up/down viceversa for above and below
	  IF (ISOLID(0).NE.0) THEN
C	    solid half-space above
	    URB0 = URVALUE (ICUTTYP,AIMCUT, U0,UPBPBSURF,UPBPB2SURF,UPREF, 
     $	    UP,UP2)
	    IF ( ISHEETTYP.EQ.3 .OR. ISHEETTYP.EQ.4 ) URB0 = -URB0
	    URB0 = -URB0	! up/down viceversa for above and below
	    G0U = TWOVB2(0)*U
	    YT10PRE = -URA0*URB0 - U2
	    YT20 = -YT10PRE*G0U - U
	    YT30 = QCMPLX(AIMAG(URB0),-REAL(URB0))
	    YT40 = -QCMPLX(AIMAG(URA0),-REAL(URA0))
	    YT60 = ( 1 + (YT20-U)*G0U ) * RHO(0)
	    CDET = YTB1(0)*YT10PRE/RHO(0) + 2*YTB2(0)*YT20 + YTB3(0)*YT30 + 
     $      YTB4(0)*YT40 + YTB6(0)*YT60
	  ELSE
C	    fluid half-space above
	    CDET = YTB1(0)*QCMPLX(AIMAG(URA0),-REAL(URA0)) + YTB3(0)*RHO(0)
	  END IF
	ELSE IF (RHO(0).EQ.0.0Q0) THEN
C         free above
	  CDET = YTB1(0)
	ELSE IF (ISOLID(1).NE.0) THEN
C	  rigid above solid, rho(0) less than zero
	  CDET = YTB6(0)
	ELSE
C	  rigid above fluid, rho(0) less than zero
	  CDET = YTB3(0)
	END IF
	IF (IGRENTRY.EQ.0) RETURN
	GOTO 200

	ENTRY GRDERIVREST ( UP, DLNDO,DERIVDO, DLNDU,DERIVDU )  ! post grvalue
	IGRENTRY = 12
	GOTO 200

	ENTRY GRDERIVDOREST ( UP, DLNDO,DERIVDO )               ! post grvalue
	IGRENTRY = 1
	GOTO 200

	ENTRY GRDERIVDUREST ( UP, DLNDU,DERIVDU )               ! post grvalue
	IGRENTRY = 2

 200	IF (IGRENTRY.NE.2) THEN
	  DLNDO = 0.0Q0
	  DERIVDO = (0.0Q0,0.0Q0)
	END IF
	IF (IGRENTRY.NE.1) THEN
	  DLNDU = 0.0Q0
	  DERIVDU = (0.0Q0,0.0Q0)
	END IF

	IF (RHO(0).GT.0.0Q0) THEN
C	  half-space present, "branch-cut(s) and sheet(s)" important
C	  ura0 already set
	  IF (ISOLID(0).NE.0) THEN
C           solid half-space above
C           urb0 and yt10pre,yt20,yt30,yt40,yt60 already set
	    IF (IGRENTRY.NE.1) THEN
	      DLNDU = TLNB(0)
	      CSLASK = U * ( URA0/URB0 + URB0/URA0 - 2 )
	      CSLASK2 = -1 - TWOVB2(0)*(YT10PRE+U*CSLASK)
	      DERIVDU = YTB1(0)*CSLASK/RHO(0) + 2*YTB2(0)*CSLASK2 +
     $        U * ( YTB3(0)/YT30 + YTB4(0)/YT40 ) + 
     $        YTB6(0)*TWOVB2(0)*(YT20+U*(CSLASK2-2))*RHO(0)
	      CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )
	    END IF
	    IF ( NSW.EQ.0 .OR. NSS.GT.0 ) THEN
	      TLN = 0.0Q0
	      YT1 = YT10PRE/RHO(0)
	      YT2 = YT20
	      YT3 = YT30
	      YT4 = YT40
	      YT6 = YT60
	      IF (NSW.EQ.0) THEN
C               special 'fix' for this case without 'nss-nsw' fluid
	        GOTO 333
	      END IF
	      IHARK0 = 0
	      IF (IGRENTRY.EQ.12) THEN
	        CALL HARKDERIVA ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $          TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO, DLNDU,DERIVDU )
	      ELSE IF (IGRENTRY.EQ.1) THEN
		CALL HARKDERIVADO ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $          TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO )
	      ELSE
		CALL HARKDERIVADU ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $          TLN,YT1,YT2,YT3,YT4,YT6, DLNDU,DERIVDU )
	      END IF
	      PLN = TLN
	      YP1 = -YT4
	      YP2 = YT6
	    ELSE
	      PLN = 0.0Q0
	      YP1 = -YT40
	      YP2 = YT60
	    END IF
	  ELSE
C           fluid half-space above
	    IF (IGRENTRY.NE.1) THEN
	      DLNDU = TLNB(0)
	      DERIVDU = YTB1(0)*U/QCMPLX(AIMAG(URA0),-REAL(URA0))
	      CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )
	    END IF
	    IF ( NSW.EQ.0 .OR. NSS.GT.0 ) THEN
	      TLN = 0.0Q0
	      YT1 = QCMPLX(AIMAG(URA0),-REAL(URA0))
	      YT2 = (0.0Q0,0.0Q0)
	      YT3 = RHO(0)
	      YT4 = (0.0Q0,0.0Q0)
	      YT6 = (0.0Q0,0.0Q0)
	      IF (NSW.EQ.0) THEN
C               special 'fix' for this case without 'nss-nsw' fluid
	        GOTO 333
	      END IF
	      IHARK0 = 0
	      IF (IGRENTRY.EQ.12) THEN
		CALL HARKDERIVA ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $          TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO, DLNDU,DERIVDU )
	      ELSE IF (IGRENTRY.EQ.1) THEN
		CALL HARKDERIVADO ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $          TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO )
	      ELSE
		CALL HARKDERIVADU ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $          TLN,YT1,YT2,YT3,YT4,YT6, DLNDU,DERIVDU )
	      END IF
	      PLN = TLN
	      YP1 = -YT4
	      YP2 = YT6
	    ELSE
	      PLN = 0.0Q0
	      YP1 = QCMPLX(AIMAG(URA0),-REAL(URA0))
	      YP2 = RHO(0)
	    END IF
	  END IF
	ELSE
C         free (or rigid) above
	  IF ( NSW.EQ.0 .OR. NSS.GT.0 ) THEN
	    TLN = 0.0Q0
	    YT1 = (1.0Q0,0.0Q0)
	    YT2 = (0.0Q0,0.0Q0)
	    YT3 = (0.0Q0,0.0Q0)
	    YT4 = (0.0Q0,0.0Q0)
	    YT6 = (0.0Q0,0.0Q0)
	    IF (RHO(0).LT.0.0Q0) THEN
	      YT1 = (0.0Q0,0.0Q0)
	      YT6 = (1.0Q0,0.0Q0)
	    END IF
	    IF (NSW.EQ.0) THEN
C             special 'fix' for this case without 'nss-nsw' fluid
	      GOTO 333
	    END IF
	    IHARK0 = 0
	    IF (IGRENTRY.EQ.12) THEN
	      CALL HARKDERIVA ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $        TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO, DLNDU,DERIVDU )
	    ELSE IF (IGRENTRY.EQ.1) THEN
	      CALL HARKDERIVADO ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $        TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO )
	    ELSE
	      CALL HARKDERIVADU ( NSS,1, OMEG,OMEG2,U,U2,RK,RK2, IHARK0, 
     $        TLN,YT1,YT2,YT3,YT4,YT6, DLNDU,DERIVDU )
	    END IF
	    PLN = TLN
	    YP1 = -YT4
	    YP2 = YT6
	  ELSE
	    PLN = 0.0Q0
	    YP1 = (1.0Q0,0.0Q0)
	    YP2 = (0.0Q0,0.0Q0)
	    IF (RHO(0).LT.0.0Q0) THEN
	      YP1 = (0.0Q0,0.0Q0)
	      YP2 = (1.0Q0,0.0Q0)
	    END IF
	  END IF
	END IF

	IF (NSW.GT.NSS) THEN
	  IF (IGRENTRY.EQ.12) THEN
	    CALL LIQDERIVA ( NSW,NSS+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    PLN,YP1,YP2, DLNDO,DERIVDO, DLNDU,DERIVDU )
	  ELSE IF (IGRENTRY.EQ.1) THEN
	    CALL LIQDERIVADO ( NSW,NSS+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    PLN,YP1,YP2, DLNDO,DERIVDO )
	  ELSE
	    CALL LIQDERIVADU ( NSW,NSS+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    PLN,YP1,YP2, DLNDU,DERIVDU )
	  END IF
	ELSE
	  STOP 64430  ! should not occur
	END IF

 333	IF (RHO(NOL).GT.0.0Q0) THEN
C	  half-space present, "branch-cut(s) and sheet(s)" important
C	  uranol already set
	  IF (ISOLID(NOL).NE.0) THEN
C	    solid half-space below
C	    urbnol already set
	    IF (NSW.GT.0) THEN
	      TLN = PLN
	      YT1 = YP1
	      YT2 = (0.0Q0,0.0Q0)
	      YT3 = YP2
	      YT4 = (0.0Q0,0.0Q0)
	      YT6 = (0.0Q0,0.0Q0)
	    ELSE
C	      "tln yt" already computed by special 'fix' for no 'nss-nsw' fluid
	    END IF
	    IHARK0 = 1
	    IF (IGRENTRY.EQ.12) THEN
	      CALL HARKDERIVA ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, 
     $	      DLNDO,DERIVDO, DLNDU,DERIVDU )
	    ELSE IF (IGRENTRY.EQ.1) THEN
	      CALL HARKDERIVADO ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO )
	    ELSE
	      CALL HARKDERIVADU ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	      IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, DLNDU,DERIVDU )
	    END IF
	    IF (IGRENTRY.NE.1) THEN
	      DLNDULOK = TLN
	      CSLASK = -YT3/URANOL + YT4/URBNOL
	      CSLASK = QCMPLX ( -AIMAG(CSLASK), REAL(CSLASK) )
	      DERIVDULOK = U * ( CSLASK + (URANOL/URBNOL+URBNOL/URANOL)*YT6 )
	    END IF
	  ELSE
C	    fluid half-space below
	    IF ( NOL .GT. (NSW+1) ) THEN
	      IF (NSW.GT.0) THEN
	        TLN = PLN
	        YT1 = YP1
	        YT2 = (0.0Q0,0.0Q0)
	        YT3 = YP2
	        YT4 = (0.0Q0,0.0Q0)
	        YT6 = (0.0Q0,0.0Q0)
	      ELSE
C	        "tln yt" already cmptd by special 'fix' for no 'nss-nsw' fluid
	      END IF
	      IHARK0 = 0
	      IF (IGRENTRY.EQ.12) THEN
	        CALL HARKDERIVA ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	        IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, 
     $	        DLNDO,DERIVDO, DLNDU,DERIVDU )
	      ELSE IF (IGRENTRY.EQ.1) THEN
	        CALL HARKDERIVADO ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	        IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO )
	      ELSE
	        CALL HARKDERIVADU ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	        IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, DLNDU,DERIVDU )
	      END IF
	      IF (IGRENTRY.NE.1) THEN
	        DLNDULOK = TLN
		DERIVDULOK = U * YT6 / QCMPLX(AIMAG(URANOL),-REAL(URANOL))
	      END IF
	    ELSE
	      IF (IGRENTRY.NE.1) THEN
	        DLNDULOK = PLN
		DERIVDULOK = -U * YP2 / QCMPLX(AIMAG(URANOL),-REAL(URANOL))
	      END IF
	    END IF
	  END IF
	  IF (IGRENTRY.NE.1) THEN
	    CALL ADDEC ( DLNDU,DERIVDU, DLNDULOK,DERIVDULOK, DLNDU,DERIVDU )
	    CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )
	  END IF
	ELSE IF ( NOL .GT. (NSW+1) ) THEN
C         rigid or free below
	  IF (NSW.GT.0) THEN
	    TLN = PLN
	    YT1 = YP1
	    YT2 = (0.0Q0,0.0Q0)
	    YT3 = YP2
	    YT4 = (0.0Q0,0.0Q0)
	    YT6 = (0.0Q0,0.0Q0)
	  ELSE
C	    "tln yt" already computed by special 'fix' for no 'nss-nsw' fluid
	  END IF
	  IF (RHO(NOL).LT.0.0Q0) THEN
	    IHARK0 = -1
	  ELSE
	    IHARK0 = -2
	  END IF
	  IF (IGRENTRY.EQ.12) THEN
	    CALL HARKDERIVA ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, 
     $	    DLNDO,DERIVDO, DLNDU,DERIVDU )
	  ELSE IF (IGRENTRY.EQ.1) THEN
	    CALL HARKDERIVADO ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO )
	  ELSE
	    CALL HARKDERIVADU ( NOL-1,NSW+1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	    IHARK0, TLN,YT1,YT2,YT3,YT4,YT6, DLNDU,DERIVDU )
	  END IF
	END IF

	IF (IGRENTRY.NE.2) THEN
	  SLASK = MAX ( ABS(REAL(DERIVDO)), ABS(AIMAG(DERIVDO)) )
	  IF (SLASK.GT.0.0Q0) THEN
	    SLASK = 1/SLASK
	    DLNDO = DLNDO - LOG(SLASK)
	    DERIVDO = SLASK*DERIVDO
	  END IF
	END IF
	IF (IGRENTRY.NE.1) THEN
	  SLASK = MAX ( ABS(REAL(DERIVDU)), ABS(AIMAG(DERIVDU)) )
	  IF (SLASK.GT.0.0Q0) THEN
	    SLASK = 1/SLASK
	    DLNDU = DLNDU - LOG(SLASK)
	    DERIVDU = SLASK*DERIVDU
	  END IF
	END IF

	RETURN
	END
