C	routine called from grderiv

C	Sven Ivansson 1993



	SUBROUTINE HARKDERIVB ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2, 
     $  IHARK0 )
C	"upwards"
C	"gr input/output" in "common /grlocal/"  (input only "to start")

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       nol <= noldim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'DPCFSUB.INC'
	INCLUDE 'FSUBCOM.INC'
	INCLUDE 'GRLOCAL.INC'

	INTEGER NLNER,NLUPP, IHARK0
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2

	INTEGER ID
	REAL*8 TLN, E2T, RH1,RH2, TLNREF
	COMPLEX*16 YT1,YT2,YT3,YT4,YT6, 
     $	G0, G0U, URA,URB, URAURB,URANEGI,URBNEGI, 
     $	CP,CQ, STP,STQ,SDP,SDQ, 
     $	T1Q,T2Q,T3Q,T4Q,T6Q
	LOGICAL EXPTYP

	IF (IHARK0.EQ.1) THEN
C	  "unfinished" start vector from homogeneous solid half-space (when 
C	  "id=nlner+1")
	  ID = NLNER+1
	  TLN = TLNB(ID)
	  YT1 = YTB1N(ID)
	  YT2 = YTB2N(ID)
	  YT3 = YTB3N(ID)
	  YT4 = YTB4N(ID)
	  YT6 = YTB6N(ID)
	  RH1 = RHO(ID)                          ! local rho(id)
	  RH2 = 1/RH1                            ! could be stored
	  G0 = TWOVB2(ID)                        ! local twovb2(id)
	  G0U = G0*U
C	  now the "rho*Cdelta propagation step from the homogeneous solid 
C	  half-space" recognizing that "start yt1" is one and "start yt2" 
C	  is zero in this case
	  YT6 = YT6 - U2
	  YT2 = -YT6*G0U - U
	  YT1 = 1 + (YT2-U)*G0U
	  YT1 = RH1*YT1
	  YT6 = RH2*YT6
	  ID = NLNER
	  TLNB(ID) = TLN
	  YTB1(ID) = YT1
	  YTB2(ID) = YT2
	  YTB3(ID) = YT3
	  YTB4(ID) = YT4
	  YTB6(ID) = YT6
	ELSE
C	  "ordinary" start vector  ( ihark0 = 0, -1, -2 )
	  ID = NLNER
	  TLN = TLNB(ID)
	  YT1 = YTB1(ID)
	  YT2 = YTB2(ID)
	  YT3 = YTB3(ID)
	  YT4 = YTB4(ID)
	  YT6 = YTB6(ID)
	END IF

	DO ID = NLNER, NLUPP, -1

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

	    IF ( IFSMTYP(ID).NE.-6 .AND. IFSMTYP(ID).NE.(-7) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
	      CALL HARKSIPART ( OMEG,OMEG2,U,U2, 
     $	      UA2(ID),UB2(ID), DZ(ID), 
     $	      URA,URB, EXPTYP, 
     $	      TLNREF,E2T, CP,CQ, STP,STQ,SDP,SDQ )            ! freq-dependent
	      URA_GR(ID) = URA
	      URB_GR(ID) = URB
	      EXPTYP_GR(ID) = EXPTYP
	      TLNREF_GR(ID) = TLNREF
	      E2T_GR(ID) = E2T
	      CP_GR(ID) = CP
	      CQ_GR(ID) = CQ
	      STP_GR(ID) = STP
	      STQ_GR(ID) = STQ
	      SDP_GR(ID) = SDP
	      SDQ_GR(ID) = SDQ
	      IF (EXPTYP) THEN
	        URAURB = URA*URB
	        URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	        URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb
	      END IF
	    ELSE
	      URA_GR(ID) = URA_GR(ID+1)
	      URB_GR(ID) = URB_GR(ID+1)
	      EXPTYP_GR(ID) = EXPTYP_GR(ID+1)
	      TLNREF_GR(ID) = TLNREF_GR(ID+1)
	      E2T_GR(ID) = E2T_GR(ID+1)
	      CP_GR(ID) = CP_GR(ID+1)
	      CQ_GR(ID) = CQ_GR(ID+1)
	      STP_GR(ID) = STP_GR(ID+1)
	      STQ_GR(ID) = STQ_GR(ID+1)
	      SDP_GR(ID) = SDP_GR(ID+1)
	      SDQ_GR(ID) = SDQ_GR(ID+1)
	    END IF
	    RH1 = RHO(ID)                        ! local rho(id)
	    RH2 = 1/RH1                          ! could be stored
	    G0 = TWOVB2(ID)                      ! local twovb2(id)
	    G0U = G0*U
C	    now the Cinvdelta propagation step (except for "rh1 vector mult")
	    YT1 = RH2*YT1
	    YT6 = RH1*YT6
	    T2Q = YT2
	    YT2 = YT2+YT2         ! temporary doubling
	    YT1 = YT1 - (YT2+YT6*G0U)*G0U
	    T6Q = YT6*G0U
	    T1Q = YT1*U + T6Q
	    YT6 = YT6 - (YT2+T1Q+T6Q)*U
	    YT2 = T2Q + T1Q
	    YTB1K(ID) = YT1
	    YTB2K(ID) = YT2
	    YTB3K(ID) = YT3
	    YTB4K(ID) = YT4
	    YTB6K(ID) = YT6
	    IF (EXPTYP) THEN
C	      now the Dinvdeltamod propagation step
	      YT1 = -YT1*URAURB
	      YT3 = YT3*URBNEGI
	      YT4 = YT4*URANEGI
	      T1Q = YT1 + YT3
	      T2Q = YT1 - YT3
	      T3Q = YT4 + YT6
	      T4Q = YT4 - YT6
	      YT1 = T1Q - T4Q
	      YT3 = -T2Q + T3Q
	      YT4 = -T1Q - T4Q
	      YT6 = -T2Q - T3Q
	      YTB1L(ID) = YT1
	      YTB2L(ID) = YT2
	      YTB3L(ID) = YT3
	      YTB4L(ID) = YT4
	      YTB6L(ID) = YT6
C	      now the PEdelta propagation step, freq-dependent
	      TLN = TLN + TLNREF
	      YT2 = E2T*YT2
	      YT1 = YT1*STP
	      YT3 = YT3*SDP
	      YT4 = YT4*STQ
	      YT6 = YT6*SDQ
C	      now the Ddeltamod propagation step
	      T1Q = YT1 + YT3
	      T2Q = YT1 - YT3
	      T3Q = YT4 + YT6
	      T4Q = YT4 - YT6
	      YT1 = T2Q - T3Q
	      YT3 = T1Q - T4Q
	      YT4 = -T2Q - T3Q
	      YT6 = T1Q + T4Q
	      YTB1M(ID) = YT1
	      YTB2M(ID) = YT2
	      YTB3M(ID) = YT3
	      YTB4M(ID) = YT4
	      YTB6M(ID) = YT6
	      TLN = TLN - DLOG4
	      YT1 = -YT1/URAURB  ! "1/uraurb=-uranegiinv*urbnegiinv" to hdera
	      YT2 = 4*YT2
	      YT3 = YT3/URBNEGI     ! could store "1/urbnegi" to harkderiva
	      YT4 = YT4/URANEGI     ! could store "1/uranegi" to harkderiva
	    ELSE
C	      now the PTdelta propagation step, freq-dependent
	      TLN = TLN + TLNREF
	      YT2 = E2T*YT2   ! decoupled here
	      T1Q =  YT1*CQ  + YT4*SDQ
	      T3Q =  YT3*CQ  - YT6*SDQ
	      T4Q = -YT1*STQ + YT4*CQ
	      T6Q =  YT3*STQ + YT6*CQ
	      YTB1L(ID) = T1Q
	      YTB2L(ID) = YT2      ! somewhat in advance
	      YTB3L(ID) = T3Q
	      YTB4L(ID) = T4Q
	      YTB6L(ID) = T6Q
	      YT1 =  T1Q*CP  - T3Q*SDP
	      YT3 =  T1Q*STP + T3Q*CP
	      YT6 = -T4Q*STP + T6Q*CP
	      YT4 =  T4Q*CP  + T6Q*SDP
	    END IF
	    YTB1N(ID) = YT1
	    YTB2N(ID) = YT2
	    YTB3N(ID) = YT3
	    YTB4N(ID) = YT4
	    YTB6N(ID) = YT6
C	    now the Cdelta propagation step (except for "rh2 vector mult")
	    T2Q = YT2
	    YT2 = YT2+YT2           ! temporary doubling
	    YT6 = YT6 + (YT2-YT1*U)*U
	    T1Q = YT1*U
	    T6Q = YT6*G0U + T1Q
	    YT1 = YT1 + (YT2-T1Q-T6Q)*G0U
	    YT2 = T2Q - T6Q
	    YT1 = RH1*YT1
	    YT6 = RH2*YT6

	  ELSE

	    stop 51267                  !!! not implemented yet !!!

	  END IF

	  IF (ID.GT.NLUPP) CALL SCALE5EC (tspmaxlim1,tspmaxlim2, 
     $	  TLN,YT1,YT2,YT3,YT4,YT6 )

	  TLNB(ID-1) = TLN     ! obs, for ytb & ytbk & ytbl & ytbm & ytbn
	  YTB1(ID-1) = YT1
	  YTB2(ID-1) = YT2
	  YTB3(ID-1) = YT3
	  YTB4(ID-1) = YT4
	  YTB6(ID-1) = YT6

	END DO

	RETURN
	END