
C	additional subroutines (and functions) for "Full Wave" computations 
C	fwderiv and liqderivfw, solderivfw
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 FWDERIV ( UP, DLNDO,DERIVDO, DLNDU,DERIVDU )

	USE INMODELa
	USE FWYSa

	IMPLICIT NONE

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

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

	COMPLEX*16 UP

	REAL*8 DLNDO, DLNDU
	COMPLEX*16 DERIVDO, DERIVDU

	REAL*8 DLNDULOK, PLN, SLN, RH2, SLASK
	COMPLEX*16 DERIVDULOK, 
     $  YP2, YSA1,YSA2,YSA3,YSA4, YSB1,YSB2,YSB3,YSB4, 
     $	U,RK,UP2,U2,RK2, URA,URB, URANEGI,URBNEGI, G0,G0U,G1,G2, 
     $	ALFAPDOWN,BETASDOWN, 
     $	URVALUE

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

	UCOP = U
	U2COP = U2
	RK2GOMEG = RK2/OMEG

	DLNDO = 0.0D0
	DERIVDO = (0.0D0,0.0D0)
	DLNDU = 0.0D0
	DERIVDU = (0.0D0,0.0D0)

	IF (NSS.GT.0) THEN
C	  upper solid
	  CALL SOLDERIVFW ( 1,NSS, OMEG,OMEG2,U,U2,RK,RK2, 
     $	  DLNDO,DERIVDO, DLNDU,DERIVDU )
	END IF

	IF (NSW.GT.NSS) THEN
C	  fluid
	  CALL LIQDERIVFW ( NSS+1,NSW, OMEG,OMEG2,U,U2,RK,RK2, 
     $	  DLNDO,DERIVDO, DLNDU,DERIVDU )
	END IF

	IF ( NOL .GT. (NSW+1) ) THEN
C	  lower solid
	  CALL SOLDERIVFW ( NSW+1,NOL-1, OMEG,OMEG2,U,U2,RK,RK2, 
     $	  DLNDO,DERIVDO, DLNDU,DERIVDU )
	END IF

	IF (RHO(NOL).GT.0.0D0) THEN
C	  lower boundary term
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
C	    solid half-space below
	    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
	    SLN = 0.0D0
	    YSA1 = YS1FW(NOL-1)
	    YSA2 = YS2FW(NOL-1)
	    YSA3 = YS3FW(NOL-1)
	    YSA4 = YS4FW(NOL-1)
	    RH2 = 1/RHO(NOL)                           ! could be stored
	    G0 = TWOVB2(NOL)                           ! local twovb2(nol)
	    G0U = G0*U
	    G1 = G0*U2
	    G2 = 1-G1
	    ALFAPDOWN = 0.5D0 * ( G0U*YSA1 + (G2/URANEGI)*YSA2 + 
     $	    RH2*(YSA4-(U/URANEGI)*YSA3) )
	    BETASDOWN = 0.5D0 * ( (G2/URBNEGI)*YSA1 + G0U*YSA2 + 
     $	    RH2*(YSA3-(U/URBNEGI)*YSA4) )
	    DLNDULOK = 0.0D0
	    DERIVDULOK = RHO(NOL)*U * ( ALFAPDOWN*ALFAPDOWN/URANEGI + 
     $	    BETASDOWN*BETASDOWN/URBNEGI )
	    YSA3 = RH2 * YSA3
	    YSA4 = RH2 * YSA4
	    YSB1 = -RHO(NOL)*BETASDOWN
	    YSB2 = -RHO(NOL)*ALFAPDOWN
	    YSB3 = -YSB1*URBNEGI
	    YSB4 = -YSB2*URANEGI
	    DERIVDULOK = DERIVDULOK - YSB1*YSA4 - YSB2*YSA3 + 
     $	    G0 * ( (-2*U*YSB1+YSB4)*YSA1 + (-2*U*YSB2+YSB3)*YSA2 ) 
C	              ! 2*u could be stored                       ! final 
	  ELSE
C	    fluid half-space below (below fluid or solid)
	    PLN = 0.0D0
	    YP2 = YS4FW(NOL-1)
	    DLNDULOK = 0.0D0
	    DERIVDULOK = YP2*YP2 * (U/URANEGI) / RHO(NOL)
	  END IF
	  CALL ADDEC ( DLNDU,DERIVDU, DLNDULOK,DERIVDULOK, 
     $	  DLNDU,DERIVDU )
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )
	END IF

	IF (RHO(0).GT.0.0D0) THEN
C	  upper boundary term
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
C	    solid half-space above
	    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
	    SLN = 0.0D0
	    YSA1 = YS1FW(0)
	    YSA2 = YS2FW(0)
	    YSA3 = YS3FW(0)
	    YSA4 = YS4FW(0)
	    RH2 = 1/RHO(0)                             ! could be stored
	    G0 = TWOVB2(0)                             ! local twovb2(0)
	    G0U = G0*U
	    G1 = G0*U2
	    G2 = 1-G1
	    ALFAPDOWN = 0.5D0 * ( G0U*YSA1 + (G2/URANEGI)*YSA2 + 
     $	    RH2*(YSA4-(U/URANEGI)*YSA3) )
	    BETASDOWN = 0.5D0 * ( (G2/URBNEGI)*YSA1 + G0U*YSA2 + 
     $	    RH2*(YSA3-(U/URBNEGI)*YSA4) )
	    DLNDULOK = 0.0D0
	    DERIVDULOK = RHO(0)*U * ( -ALFAPDOWN*ALFAPDOWN/URANEGI - 
     $	    BETASDOWN*BETASDOWN/URBNEGI )
	    YSA3 = RH2 * YSA3
	    YSA4 = RH2 * YSA4
	    YSB1 = -RHO(0)*BETASDOWN
	    YSB2 = -RHO(0)*ALFAPDOWN
	    YSB3 = -YSB1*URBNEGI
	    YSB4 = -YSB2*URANEGI
	    DERIVDULOK = DERIVDULOK + YSB1*YSA4 + YSB2*YSA3 - 
     $	    G0 * ( (-2*U*YSB1+YSB4)*YSA1 + (-2*U*YSB2+YSB3)*YSA2 ) 
C	              ! 2*u could be stored                       ! final 
	  ELSE
C	    fluid half-space above (above fluid or solid)
	    PLN = 0.0D0
	    YP2 = YS4FW(0)
	    DLNDULOK = 0.0D0
	    DERIVDULOK = -YP2*YP2 * (U/URANEGI) / RHO(0)
	  END IF
	  CALL ADDEC ( DLNDU,DERIVDU, DLNDULOK,DERIVDULOK, DLNDU,DERIVDU )
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )
	END IF

	SLASK = DMAX1 ( DABS(DREAL(DERIVDO)), DABS(DIMAG(DERIVDO)) )
	IF (SLASK.GT.0.0D0) THEN
	  SLASK = 1/SLASK
	  DLNDO = DLNDO - DLOG(SLASK)
	  DERIVDO = SLASK*DERIVDO
	END IF
	SLASK = DMAX1 ( DABS(DREAL(DERIVDU)), DABS(DIMAG(DERIVDU)) )
	IF (SLASK.GT.0.0D0) THEN
	  SLASK = 1/SLASK
	  DLNDU = DLNDU - DLOG(SLASK)
	  DERIVDU = SLASK*DERIVDU
	END IF

	RETURN
	END








	SUBROUTINE LIQDERIVFW ( NL1,NL2, OMEG,OMEG2,U,U2,RK,RK2, 
     $	DLNDO,DERIVDO, DLNDU,DERIVDU )
C	"fw input" in "common /fwpart,fwys/"
C	input/output: dlndo,derivdo, dlndu,derivdu

	USE INMODELa
	USE FWYSa
	USE FWPARTa

	IMPLICIT NONE

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

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

	INTEGER NL1,NL2
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2

	REAL*8 DLNDO, DLNDU
	COMPLEX*16 DERIVDO, DERIVDU

	INTEGER ID
	REAL*8 DLNDOLOK, DLNDULOK, PLNA, PLNB, PLNREF, PI
	COMPLEX*16 DERIVDOLOK, DERIVDULOK, YPA1,YPA2, YPB1,YPB2, 
     $	URAURA, URA, URANEGI, CSLASK0,CSLASK, COSSINC, 
     $	CP,CPX, STP,SDP, WWLOK, 
     $	CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	DCP, OMEGDCP, DSTP, OMEGDSDP, 
     $	PIOMEGTHREEINV, PIOMEG2UTWO, PIOMEG2UTWOCFAKT, 
     $	YACP,YBCPX, BCFAKTINVSTP, BCFAKTSDP,BCFAKTSDPYA,BCFAKTSDPYB
	LOGICAL EXPTYP

	DATA PI /3.141592653589793D0/

	IF (NL1.GT.NL2) STOP 3456

	PIOMEGTHREEINV = PI / (3*OMEG)  ! could be computed at "higher level"
	PIOMEG2UTWO = PI * 2*OMEG2*U

	DO ID = NL1, NL2

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

	    IF ( IFSMTYP(ID).NE.-6 .AND. IFSMTYP(ID).NE.(-5) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
	      URAURA = UA2(ID) - U2
C	      stored "ura" (freq-indep) and "call liqpart" (freq-dep)
	      URA = URA_FW(ID)
	      EXPTYP = EXPTYP_FW(ID)
	      IF (EXPTYP) THEN
	        URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura
	      ELSE
	        PLNREF = TLNREF_FW(ID)
	        CP = CP_FW(ID)
	        STP = STP_FW(ID)
	        SDP = SDP_FW(ID)
	      END IF
	    END IF
	    IF (EXPTYP) THEN
	      DLNDOLOK = 0.0D0
	      DLNDULOK = DLNDOLOK
C	      first
	      YPA1 = YS2FW(ID)
	      YPA2 = YS4FW(ID)
	      YPA2 = YPA2 / (2*RHO(ID))
	      WWLOK = YPA1 / (2*URANEGI)
	      YPA1 = WWLOK + YPA2
	      YPA2 = -WWLOK + YPA2
	      YPB1 = -YS4FW(ID)   ! obs
	      YPB2 = YS2FW(ID)    ! obs
	      YPB2 = YPB2 * RHO(ID)
	      WWLOK = YPB1*URANEGI
	      YPB1 = WWLOK + YPB2
	      YPB2 = -WWLOK + YPB2
	      DERIVDOLOK = -YPB1*YPA1 + YPB2*YPA2
	      DERIVDULOK = OMEG*DZ(ID)*DERIVDOLOK
	      DERIVDOLOK = DZ(ID)*URANEGI * DERIVDOLOK           ! final
C	      second
	      YPB1 = -YS4FW(ID)   ! obs
	      YPB2 = YS2FW(ID)    ! obs
	      DERIVDULOK = DERIVDULOK + YPB1*(YPA1-YPA2)
C	      third
	      YPA1 = YS2FW(ID-1)
	      YPA2 = YS4FW(ID-1)
	      YPB1 = -YS4FW(ID-1)   ! obs
	      YPB2 = YS2FW(ID-1)    ! obs
	      YPB2 = YPB2 * RHO(ID)
	      WWLOK = YPB1*URANEGI
	      YPB1 = WWLOK + YPB2
	      YPB2 = -WWLOK + YPB2
	      DERIVDULOK = DERIVDULOK + (YPB1-YPB2)*YPA1/(2*URAURA)
	      DERIVDULOK = (U/URANEGI) * DERIVDULOK              ! final
	    ELSE
	      PLNA = 0.0D0
	      YPA1 = YS2FW(ID-1)
	      YPA2 = YS4FW(ID-1)
	      CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, PLNA,YPA1,YPA2 )
	      PLNB = 0.0D0
	      YPB1 = -YS4FW(ID)   ! obs
	      YPB2 = YS2FW(ID)    ! obs
	      CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, PLNB,YPB1,YPB2 )
	      YPA2 = YPA2 / RHO(ID)  !obs,ypa2 changed !"1/rho" could be stored
	      YPB2 = YPB2 * RHO(ID)  ! obs, already here
	      IF (PLNREF.GT.XSPLITLIM) THEN
	        WRITE(*,*) '%%% liqderivfw, depths ',ZD(ID-1),ZD(ID)
	        WRITE(*,*) 'quotient u ',PLNREF/XSPLITLIM,U
	        WRITE(*,*) 'split needed, automatic split not implemented %%%'
	        STOP 3141
	      END IF
	      DLNDOLOK = PLNA + PLNREF + PLNB
	      DLNDULOK = DLNDOLOK
C	      derivative-term computation
	      DCP = DZ(ID)*CP
	      OMEGDCP = OMEG*DCP
	      DSTP = DZ(ID)*STP
	      OMEGDSDP = OMEG*DZ(ID)*SDP       ! "omeg*dz(id)" could be stored
	      DERIVDOLOK = YPB1 * ( -DSTP*YPA1 + URAURA*DCP*YPA2 ) - 
     $	      YPB2 * ( DCP*YPA1 + DSTP*YPA2 )                    ! final
	      IF ( DZ(ID).EQ.0.0D0 .OR. OMEG.EQ.(0.0D0,0.0D0) ) THEN
	        CSLASK = (0.0D0,0.0D0)
	      ELSE IF ( URA .EQ. (0.0D0,0.0D0) ) THEN
	        CSLASK = OMEG * DZ(ID)
	        CSLASK = -CSLASK*CSLASK*CSLASK/3
	      ELSE
	        CSLASK0 = OMEG * DZ(ID)
	        CSLASK = CSLASK0 * URA
	        IF ( MAX(ABS(REAL(CSLASK)),ABS(AIMAG(CSLASK))) .GT.
     &	         0.5D0 ) THEN
	          CSLASK = (OMEGDCP-SDP)/URAURA
	        ELSE
	          CSLASK = CSLASK0*CSLASK0*CSLASK0 * COSSINC(CSLASK)
	          IF (PLNREF.NE.0.0D0) CSLASK = DEXP(-PLNREF) * CSLASK
	        END IF
	      END IF  ! now cslask = (omegdcp-sdp)/uraura
	      DERIVDULOK = U * ( YPB1 * ( OMEGDSDP*YPA1 - (SDP+OMEGDCP)*
     $	      YPA2 ) + YPB2 * ( CSLASK*YPA1 + OMEGDSDP*YPA2 ) )  ! final
	    END IF

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

	    PLNA = 0.0D0
	    YPA1 = YS2FW(ID-1)
	    YPA2 = YS4FW(ID-1)
	    CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, PLNA,YPA1,YPA2 )
	    PLNB = 0.0D0
	    YPB1 = -YS4FW(ID)   ! obs
	    YPB2 = YS2FW(ID)    ! obs
	    CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, PLNB,YPB1,YPB2 )
C	    stored "call airypart"
	    CFAKT = CFAKT_FW(ID)
	    BCFAKT = BCFAKT_FW(ID)
	    BCFAKTINV = BCFAKTINV_FW(ID)
	    YA = YA_FW(ID)
	    YB = YB_FW(ID)
	    PLNREF = TLNREF_FW(ID)
	    CP = CP_FW(ID)
	    CPX = CPX_FW(ID)
	    STP = STP_FW(ID)
	    SDP = SDP_FW(ID)
	    PIOMEG2UTWOCFAKT = PIOMEG2UTWO*CFAKT
	    YPA2 = YPA2 / RHO(ID)  !obs,ypa2 changed  !"1/rho" could be stored
	    YPB2 = YPB2 * RHO(ID)   ! obs, already here
	    IF (PLNREF.GT.XSPLITLIM) THEN
	      WRITE(*,*) '%%% liqderivfw, depths ',ZD(ID-1),ZD(ID)
	      WRITE(*,*) 'quotient ',PLNREF/XSPLITLIM
	      WRITE(*,*) 'split needed, automatic split not implemented %%%'
	      STOP 3142
	    END IF
	    DLNDOLOK = PLNA + PLNREF + PLNB
	    DLNDULOK = DLNDOLOK
C	    derivative-term computation
	    YACP = YA*CP
	    YBCPX = YB*CPX
	    BCFAKTINVSTP = BCFAKTINV*STP
	    BCFAKTSDP = BCFAKT*SDP
	    BCFAKTSDPYA = BCFAKTSDP*YA
	    BCFAKTSDPYB = BCFAKTSDP*YB
	    DERIVDOLOK = PIOMEGTHREEINV * ( 
     $	    YPB1 * ( (-BCFAKTINVSTP*YA-BCFAKTSDPYB*YB)*2*YPA1 + 
     $	    (-STP+2*BCFAKT*(-YACP*YA+YBCPX*YB))*YPA2 ) + 
     $	    YPB2 * ( (-SDP-2*BCFAKTINV*(YA*CPX-YB*CP))*YPA1 + 
     $	    (BCFAKTSDPYA*YA+BCFAKTINVSTP*YB)*2*YPA2 ) )      ! final
	    DERIVDULOK = PIOMEG2UTWOCFAKT * ( 
     $	    YPB1 * ( (-BCFAKTINVSTP-BCFAKTSDPYB)*YPA1 + 
     $	    (BCFAKT*(-YACP+YBCPX))*YPA2 ) + 
     $	    YPB2 * ( -BCFAKTINV*(CPX-CP)*YPA1 + 
     $	    (BCFAKTSDPYA+BCFAKTINVSTP)*YPA2 ) )              ! final

	  ELSE

	    STOP 51267                !!! not implemented yet         !!!

	  END IF

	  CALL ADDEC ( DLNDO,DERIVDO, DLNDOLOK,DERIVDOLOK, DLNDO,DERIVDO )
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDO,DERIVDO )

	  CALL ADDEC ( DLNDU,DERIVDU, DLNDULOK,DERIVDULOK, DLNDU,DERIVDU )
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )

	END DO

	RETURN
	END







	SUBROUTINE SOLDERIVFW ( NL1,NL2, OMEG,OMEG2,U,U2,RK,RK2, 
     $	DLNDO,DERIVDO, DLNDU,DERIVDU )
C	"fw input" in "common /fwpart,fwys/"
C	input/output: dlndo,derivdo, dlndu,derivdu

	USE INMODELa
	USE FWYSa
	USE FWPARTa

	IMPLICIT NONE

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

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

	INTEGER NL1,NL2
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2

	REAL*8 DLNDO, DLNDU
	COMPLEX*16 DERIVDO, DERIVDU

	INTEGER ID
	REAL*8 DLNDOLOK, DLNDULOK, SLNA, SLNB, SLNREF,E1X
	COMPLEX*16 DERIVDOLOK, DERIVDULOK, 
     $	YSA1,YSA2,YSA3,YSA4, YSB1,YSB2,YSB3,YSB4, 
     $	YSA1K,YSA2K,YSA3K,YSA4K, YSA1N,YSA2N,YSA3N,YSA4N, 
     $	G0, G0U, 
     $	URAURA,URBURB, URA,URB, URANEGI,URBNEGI, 
     $	CP,CQ, STP,STQ,SDP,SDQ, 
     $	DCP,DCQ, OMEGDCP,OMEGDCQ, DSTP,DSTQ, OMEGDSDP,OMEGDSDQ, 
     $	S1Q,S2Q, CSLASK
	LOGICAL EXPTYP

	IF (NL1.GT.NL2) STOP 3456

	DO ID = NL1, NL2

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

	    IF ( IFSMTYP(ID).NE.-6 .AND. IFSMTYP(ID).NE.(-5) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
	      G0 = TWOVB2(ID)                      ! local twovb2(id)
	      G0U = G0*U
	      URAURA = UA2(ID) - U2
	      URBURB = UB2(ID) - U2
C	      stored "ura,urb" (freq-indep) and "call (h)solpart" (freq-dep)
	      URA = URA_FW(ID)
	      URB = URB_FW(ID)
	      EXPTYP = EXPTYP_FW(ID)
	      IF (EXPTYP) THEN
	        URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	        URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb
	      ELSE
	        SLNREF = 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)
	      END IF
	    END IF
	    IF (EXPTYP) THEN
	      DLNDOLOK = 0.0D0
	      DLNDULOK = DLNDOLOK
C	      first
	      YSA1 = YS1FW(ID)
	      YSA2 = YS2FW(ID)
	      YSA3 = YS3FW(ID)
	      YSA4 = YS4FW(ID)
	      YSA3 = YSA3 / RHO(ID)
	      YSA4 = YSA4 / RHO(ID)
	      YSA3 = G0U*YSA2 + YSA3
	      YSA2 = YSA2 - U*YSA3
	      YSA4 = G0U*YSA1 + YSA4
	      YSA1 = YSA1 - U*YSA4
	      YSA1 = YSA1 / URBNEGI
	      YSA2 = YSA2 / URANEGI
	      S1Q = YSA1 + YSA3
	      S2Q = YSA2 + YSA4
	      YSA3 = ( YSA1 - YSA3 ) / 2
	      YSA4 = ( -YSA2 + YSA4 ) / 2
	      YSA1 = S1Q / 2
	      YSA2 = S2Q / 2
	      YSB1 = -YS3FW(ID)   ! obs
	      YSB2 = -YS4FW(ID)   ! obs
	      YSB3 = YS1FW(ID)    ! obs
	      YSB4 = YS2FW(ID)    ! obs
	      YSB3 = YSB3 * RHO(ID)
	      YSB4 = YSB4 * RHO(ID)
	      YSB1 = YSB1 - YSB4*G0U
	      YSB2 = YSB2 - YSB3*G0U
	      YSB3 = YSB3 + YSB2*U
	      YSB4 = YSB4 + YSB1*U
	      S1Q = YSB1*URBNEGI
	      S2Q = YSB2*URANEGI
	      YSB1 = S1Q + YSB3
	      YSB2 = S2Q + YSB4
	      YSB3 = S1Q - YSB3
	      YSB4 = -S2Q + YSB4
	      S1Q = -YSB1*YSA1 + YSB3*YSA3
	      S2Q = -YSB2*YSA2 + YSB4*YSA4
	      DERIVDOLOK = DZ(ID) * ( URBNEGI*S1Q + URANEGI*S2Q )   ! final
	      DERIVDULOK = DZ(ID)*OMEG*U * ( S1Q/URBNEGI + 
     $	      S2Q/URANEGI )   ! "dz(id)*omeg" could be stored
C	      second
	      YSB1 = -YS3FW(ID)   ! obs
	      YSB2 = -YS4FW(ID)   ! obs
	      YSB3 = YS1FW(ID)    ! obs
	      YSB4 = YS2FW(ID)    ! obs
	      YSB3 = YSB3 * RHO(ID)
	      YSB4 = YSB4 * RHO(ID)
	      S1Q = YSA1 + YSA3
	      S2Q = YSA2 - YSA4
	      YSA3K = YSA1 - YSA3
	      YSA4K = YSA2 + YSA4
	      YSA1K = URBNEGI*S1Q
	      YSA2K = URANEGI*S2Q
	      DERIVDULOK = DERIVDULOK + YSB1*YSA4K + YSB2*YSA3K - 
     $	      G0 * ( YSB3*(YSA2K+2*U*YSA3K) + 
     $	      YSB4*(YSA1K+2*U*YSA4K) ) !2*u could be stored
	      YSB1 = YSB1 - YSB4*G0U
	      YSB2 = YSB2 - YSB3*G0U
	      YSB3 = YSB3 + YSB2*U
	      YSB4 = YSB4 + YSB1*U
	      S1Q = YSB1*U/URBNEGI
	      S2Q = YSB2*U/URANEGI
	      DERIVDULOK = DERIVDULOK + S1Q*(YSA1+YSA3) + S2Q*(YSA2-YSA4)
C	      third
	      YSA1 = YS1FW(ID-1)
	      YSA2 = YS2FW(ID-1)
	      YSA3 = YS3FW(ID-1)
	      YSA4 = YS4FW(ID-1)
	      YSA3 = YSA3 / RHO(ID)
	      YSA4 = YSA4 / RHO(ID)
	      YSA3K = G0U*YSA2 + YSA3
	      YSA2K = YSA2 - U*YSA3K
	      YSA4K = G0U*YSA1 + YSA4
	      YSA1K = YSA1 - U*YSA4K
	      YSB1 = -YS3FW(ID-1)   ! obs
	      YSB2 = -YS4FW(ID-1)   ! obs
	      YSB3 = YS1FW(ID-1)    ! obs
	      YSB4 = YS2FW(ID-1)    ! obs
	      YSB3 = YSB3 * RHO(ID)
	      YSB4 = YSB4 * RHO(ID)
	      YSB1 = YSB1 - YSB4*G0U
	      YSB2 = YSB2 - YSB3*G0U
	      YSB3 = YSB3 + YSB2*U
	      YSB4 = YSB4 + YSB1*U
	      S1Q = YSB1*URBNEGI
	      S2Q = YSB2*URANEGI
	      YSB1 = S1Q + YSB3
	      YSB2 = S2Q + YSB4
	      YSB3 = S1Q - YSB3
	      YSB4 = -S2Q + YSB4
	      YSA1K = YSA1K / URBNEGI
	      YSA2K = YSA2K / URANEGI
	      DERIVDULOK = DERIVDULOK + (U/2) * ( (YSB1+YSB3)*YSA1K/
     $	      URBURB + (YSB2-YSB4)*YSA2K/URAURA )
	      S1Q = YSB1 + YSB3
	      S2Q = YSB2 - YSB4
	      YSB3 = ( YSB1 - YSB3 ) / 2
	      YSB4 = ( YSB2 + YSB4 ) / 2
	      YSB1 = S1Q / (2*URBNEGI)
	      YSB2 = S2Q / (2*URANEGI)
	      DERIVDULOK = DERIVDULOK - YSB1*YSA4 - YSB2*YSA3 + 
     $	      G0 * ( (-2*U*YSB1+YSB4)*YSA1 + (-2*U*YSB2+YSB3)*YSA2 ) 
C	                ! 2*u could be stored                       ! final 
	    ELSE
	      SLNA = 0.0D0
	      YSA1 = YS1FW(ID-1)
	      YSA2 = YS2FW(ID-1)
	      YSA3 = YS3FW(ID-1)
	      YSA4 = YS4FW(ID-1)
	      CALL SCALE4EC ( tspmaxlim1,tspmaxlim2, SLNA,YSA1,YSA2,YSA3,YSA4 )
	      SLNB = 0.0D0
	      YSB1 = -YS3FW(ID)   ! obs
	      YSB2 = -YS4FW(ID)   ! obs
	      YSB3 = YS1FW(ID)    ! obs
	      YSB4 = YS2FW(ID)    ! obs
	      CALL SCALE4EC ( tspmaxlim1,tspmaxlim2, SLNB,YSB1,YSB2,YSB3,YSB4 )
C	      start column-vector propagation for ysa ("from right to left")
	      YSA3 = YSA3 / RHO(ID)          !"1/rho" could be stored
	      YSA4 = YSA4 / RHO(ID)          !"1/rho" could be stored
	      YSA3K = G0U*YSA2 + YSA3
	      YSA2K = YSA2 - U*YSA3K
	      YSA4K = G0U*YSA1 + YSA4
	      YSA1K = YSA1 - U*YSA4K
	      YSB3 = YSB3 * RHO(ID)   ! obs, already here
	      YSB4 = YSB4 * RHO(ID)   ! obs, already here
	      IF (DABS(SLNREF).GT.XSPLITLIM) THEN
	        WRITE(*,*) '%%% solderivfw, depths ',ZD(ID-1),ZD(ID)
	        WRITE(*,*) 'quotient u ',DABS(SLNREF)/XSPLITLIM,U
	        WRITE(*,*) 'split needed, automatic split not implemented %%%'
	        STOP 3141
	      END IF
	      DLNDOLOK = SLNA + DABS(SLNREF) + SLNB
	      DLNDULOK = DLNDOLOK
C	      now column-vector propagation for ysa ("from right to left")
	      S1Q = CQ*YSA1K + STQ*YSA3K
	      S2Q = CP*YSA2K + STP*YSA4K
	      YSA3N = -SDQ*YSA1K + CQ*YSA3K
	      YSA4N = -SDP*YSA2K + CP*YSA4K
	      YSA1N = S1Q
	      YSA2N = S2Q
	      IF (SLNREF.GE.0.0D0) THEN
	        YSA1N = E1X*YSA1N
	        YSA3N = E1X*YSA3N
	      ELSE
	        YSA2N = E1X*YSA2N
	        YSA4N = E1X*YSA4N
	      END IF	        
C	      now row-vector propagation for ypb ("from left to right")
C	      and derivative-term computation
	      DERIVDULOK = YSB1*YSA4N + YSB2*YSA3N - G0 * ( YSB3*
     $	      (YSA2N+2*U*YSA3N) + YSB4*(YSA1N+2*U*YSA4N) ) !2*u could be stored
	      YSB1 = YSB1 - YSB4*G0U
	      YSB2 = YSB2 - YSB3*G0U
	      YSB3 = YSB3 + YSB2*U
	      YSB4 = YSB4 + YSB1*U
	      DCQ = DZ(ID)*CQ
	      OMEGDCQ = OMEG*DCQ
	      DSTQ = DZ(ID)*STQ
	      OMEGDSDQ = OMEG*DZ(ID)*SDQ       ! "omeg*dz(id)" could be stored
	      DERIVDOLOK = YSB1 * ( -DSTQ*YSA1K + URBURB*DCQ*YSA3K ) - 
     $	      YSB3 * ( DCQ*YSA1K + DSTQ*YSA3K )
	      IF (SLNREF.GE.0.0D0) DERIVDOLOK = E1X*DERIVDOLOK
	      CSLASK = YSB1 * ( OMEGDSDQ*YSA1K - (SDQ+OMEGDCQ)*YSA3K ) + 
     $	      YSB3 * ( (OMEGDCQ-SDQ)*YSA1K/URBURB + OMEGDSDQ*YSA3K )
	      IF (SLNREF.GE.0.0D0) CSLASK = E1X*CSLASK
	      DERIVDULOK = DERIVDULOK + U*CSLASK
	      DCP = DZ(ID)*CP
	      OMEGDCP = OMEG*DCP
	      DSTP = DZ(ID)*STP
	      OMEGDSDP = OMEG*DZ(ID)*SDP       ! "omeg*dz(id)" could be stored
	      CSLASK = YSB2 * ( -DSTP*YSA2K + URAURA*DCP*YSA4K ) - 
     $	      YSB4 * ( DCP*YSA2K + DSTP*YSA4K )
	      IF (SLNREF.LT.0.0D0) CSLASK = E1X*CSLASK
	      DERIVDOLOK = DERIVDOLOK + CSLASK                   ! final
	      CSLASK = YSB2 * ( OMEGDSDP*YSA2K - (SDP+OMEGDCP)*YSA4K ) + 
     $	      YSB4 * ( (OMEGDCP-SDP)*YSA2K/URAURA + OMEGDSDP*YSA4K )
	      IF (SLNREF.LT.0.0D0) CSLASK = E1X*CSLASK
	      DERIVDULOK = DERIVDULOK + U*CSLASK
	      S1Q = YSB1*CQ - YSB3*SDQ
	      S2Q = YSB2*CP - YSB4*SDP
	      YSB3 = YSB1*STQ + YSB3*CQ
	      YSB4 = YSB2*STP + YSB4*CP
	      YSB1 = S1Q
	      YSB2 = S2Q
	      IF (SLNREF.GE.0.0D0) THEN
	        YSB1 = E1X*YSB1
	        YSB3 = E1X*YSB3
	      ELSE
	        YSB2 = E1X*YSB2
	        YSB4 = E1X*YSB4
	      END IF	        
	      DERIVDULOK = DERIVDULOK - YSB1*YSA4 - YSB2*YSA3 + 
     $	      G0 * ( (-2*U*YSB1+YSB4)*YSA1 + (-2*U*YSB2+YSB3)*YSA2 ) 
C	                ! 2*u could be stored                       ! final 
	    END IF

	  ELSE

	    STOP 51267                 !!! not implemented yet         !!!

	  END IF

	  CALL ADDEC ( DLNDO,DERIVDO, DLNDOLOK,DERIVDOLOK, DLNDO,DERIVDO )
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDO,DERIVDO )

	  CALL ADDEC ( DLNDU,DERIVDU, DLNDULOK,DERIVDULOK, DLNDU,DERIVDU )
	  CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )

	END DO

	RETURN
	END
