C	routine called from grderiv

C	Sven Ivansson 1993



	SUBROUTINE HARKDERIVA ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2, 
     $	IHARK0, TLN, YT1,YT2,YT3,YT4,YT6, 
     $	DLNDO,DERIVDO, DLNDU,DERIVDU )
C	"downwards"
C	"gr input" in "common /grlocal/"
C	input/output: tln,yt1,yt2,yt3,yt4,yt6, dlndo,derivdo, dlndu,derivdu
C	note: harkderivb must have been called in advance

	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

	REAL*8 TLN, DLNDO, DLNDU
	COMPLEX*16 YT1,YT2,YT3,YT4,YT6, DERIVDO, DERIVDU

	INTEGER ID
	REAL*8 DLNDOLOK, DLNDULOK, E2T, RH1,RH2, TMAX, TLNREF
	COMPLEX*16 DERIVDOLOK, DERIVDULOK, 
     $	G0, G0U, G1,G2,G3,G4,G0G0U,G0G3, 
     $	URAURA,URBURB, URAURAINV,URBURBINV, URA,URB, 
     $	URAURB,URANEGI,URBNEGI, URANEGIINV,URBNEGIINV, 
     $	CP,CQ, STP,STQ,SDP,SDQ, 
     $	T1Q,T2Q,T3Q,T4Q,T6Q, DZOMEG
	LOGICAL EXPTYP

	ENTRY HARKDERIVADO ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2, 
     $	IHARK0, TLN, YT1,YT2,YT3,YT4,YT6, DLNDO,DERIVDO )

	ENTRY HARKDERIVADU ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2, 
     $	IHARK0, TLN, YT1,YT2,YT3,YT4,YT6, DLNDU,DERIVDU )

	DO ID = NLUPP, NLNER

	  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       ! could be stored from harkderivb
	      URAURAINV = 1/URAURA        ! could use "-uranegiinv*uranegiinv"
	      URBURB = UB2(ID) - U2       ! could be stored from harkderivb
	      URBURBINV = 1/URBURB        ! could use "-urbnegiinv*urbnegiinv"
C	      stored "ura,urb" (freq-indep) and "call harksipart" (freq-dep)
	      URA = URA_GR(ID)
	      URB = URB_GR(ID)
	      EXPTYP = EXPTYP_GR(ID)
	      TLNREF = TLNREF_GR(ID)
	      E2T = E2T_GR(ID)
	      CP = CP_GR(ID)
	      CQ = CQ_GR(ID)
	      STP = STP_GR(ID)
	      STQ = STQ_GR(ID)
	      SDP = SDP_GR(ID)
	      SDQ = SDQ_GR(ID)
	      IF (EXPTYP) THEN
	        URAURB = URA*URB          ! could be stored from harkderivb
	        URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))   ! -i*ura
	        URBNEGI = DCMPLX(DIMAG(URB),-DREAL(URB))   ! -i*urb
	        URANEGIINV = 1/URANEGI    ! could be stored from harkderivb
	        URBNEGIINV = 1/URBNEGI    ! could be stored from harkderivb
	      END IF
	    END IF
	    DZOMEG = DZ(ID) * OMEG               ! could be stored
	    RH1 = RHO(ID)                        ! local rho(id)
	    RH2 = 1/RH1                          ! could be stored
	    G0 = TWOVB2(ID)                      ! local twovb2(id)
	    G0U = G0*U                  ! could be stored from harkderivb
	    G1 = G0*U2
	    G2 = 1-G1
	    G3 = 1-3*G1
	    G4 = 4*G0U*G2
	    G0G0U = G0*G0U
	    G0G3 = G0*G3
C	    derivative-term computation  ( ytbn(id)*Cdelta'*yta )
C	    propagation with Cdelta for yta (except for "rh2 vector-mult")
	    YT1 = RH1*YT1
	    YT6 = RH2*YT6
	    T2Q = YT2
	    YT2 = YT2+YT2         ! temporary doubling
	    IF (IGRENTRY.NE.1) DERIVDULOK = 
     $	    YTB1N(ID)*(-G4*YT1-G3*YT2-2*U*YT6) + 
     $	    2*YTB2N(ID)*(G0G3*YT1-2*G0U*YT2+YT6) + 
     $	    YTB6N(ID)*(-2*G0G0U*YT1-G0*YT2)
	    YT6 = YT6 - G0U*(YT2+G0U*YT1)
	    T1Q = G0U*YT1
	    T6Q = U*YT6 + T1Q
	    YT1 = YT1 - U*(YT2+T6Q+T1Q)
	    YT2 = T2Q + T6Q
	    IF (EXPTYP) THEN
C	      derivative-term computation  ( ytbm(id)*Ddeltamod'*yta )
C	      propagation with Ddeltamod for yta
	      YT1 = -YT1/URAURB       ! could store "1/uraurb" from harkderivb
	      YT2 = 4*YT2
	      YT3 = YT3*URBNEGIINV
	      YT4 = YT4*URANEGIINV
	      IF (IGRENTRY.NE.1) DERIVDULOK = DERIVDULOK + U * 
     $	      ( YTB1M(ID)*(URAURAINV+URBURBINV)*YT1 + 
     $	      YTB3M(ID)*URBURBINV*YT3 + YTB4M(ID)*URAURAINV*YT4 )
	      T1Q = YT3 + YT6
	      T2Q = -YT3 + YT6
	      T3Q = YT1 - YT4
	      T4Q = -YT1 - YT4
	      YT1 = T1Q + T3Q
	      YT3 = T1Q - T3Q
	      YT4 = T4Q + T2Q
	      YT6 = T4Q - T2Q
C	      derivative-term computation  ( ytbl(id)*PEdelta'*yta )
C	      propagation with PEdelta for yta
	      YT2 = E2T*YT2
	      YT1 = STP*YT1
	      YT3 = SDP*YT3
	      YT4 = STQ*YT4
	      YT6 = SDQ*YT6
	      T1Q = -YTB1L(ID)*YT1 + YTB6L(ID)*YT6  ! igrentry.ne.0
	      T2Q = YTB3L(ID)*YT3 - YTB4L(ID)*YT4   ! igrentry.ne.0
	      IF (IGRENTRY.NE.2) DERIVDOLOK = DZ(ID) * ( 
     $	      (URANEGI+URBNEGI)*T1Q + (URANEGI-URBNEGI)*T2Q )     ! final
	      IF (IGRENTRY.NE.1) DERIVDULOK = DERIVDULOK + DZOMEG*U * ( 
     $	      (URANEGIINV+URBNEGIINV)*T1Q + (URANEGIINV-URBNEGIINV)*T2Q )
C	      derivative-term computation  ( ytbk(id)*Dinvdeltamod'*yta )
C	      propagation with Dinvdeltamod for yta
	      T1Q = YT1 - YT4
	      T2Q = -YT1 - YT4
	      T3Q = -YT3 - YT6
	      T4Q = YT3 - YT6
	      YT1 = T1Q + T3Q
	      YT3 = T1Q - T3Q
	      YT4 = T4Q + T2Q
	      YT6 = T4Q - T2Q
	      YT1 = -URAURB*YT1
	      YT3 = URBNEGI*YT3
	      YT4 = URANEGI*YT4
	      IF (IGRENTRY.NE.1) DERIVDULOK = DERIVDULOK - U * 
     $	      ( YTB1K(ID)*(URAURAINV+URBURBINV)*YT1 + 
     $	      YTB3K(ID)*URBURBINV*YT3 + YTB4K(ID)*URAURAINV*YT4 )
	    ELSE
C	      < PTdelta = PTdeltaQ * PTdeltaP according to Woodhouse >
C	      derivative-term computation  ( ytbl(id)*PTdeltaP'*yta )
	      T1Q =  YTB1L(ID)*YT1 + YTB3L(ID)*YT3 + YTB4L(ID)*YT4 + 
     $	      YTB6L(ID)*YT6
	      T2Q =  YTB1L(ID)*YT3 - YTB4L(ID)*YT6
	      T3Q = -YTB3L(ID)*YT1 + YTB6L(ID)*YT4
	      IF (IGRENTRY.NE.2) DERIVDOLOK = -STP*T1Q + CP*(URAURA*T2Q+T3Q)
	      IF (IGRENTRY.NE.1) DERIVDULOK = DERIVDULOK + 
     $	      U * ( SDP * ( DZOMEG*T1Q - T2Q + URAURAINV*T3Q ) - 
     $	      DZOMEG*CP * ( T2Q + URAURAINV*T3Q ) )
C	      propagation with PTdeltaP for yta
	      YT2 = E2T*YT2   ! decoupled here, somewhat in advance
	      T1Q =    CP*YT1 + STP*YT3
	      YT3 =  -SDP*YT1 +  CP*YT3
	      T4Q =    CP*YT4 - STP*YT6
	      YT6 =   SDP*YT4 +  CP*YT6
	      YT1 = T1Q
	      YT4 = T4Q
C	      derivative-term computation  ( ytbk(id)*PTdeltaQ'*yta )
	      T1Q =  YTB1K(ID)*YT1 + YTB3K(ID)*YT3 + YTB4K(ID)*YT4 + 
     $	      YTB6K(ID)*YT6
	      T2Q = -YTB1K(ID)*YT4 + YTB3K(ID)*YT6
	      T3Q =  YTB4K(ID)*YT1 - YTB6K(ID)*YT3
	      IF (IGRENTRY.NE.2) DERIVDOLOK = DZ(ID) * ( DERIVDOLOK - 
     $	      STQ*T1Q + CQ*(URBURB*T2Q+T3Q) )                 ! final
	      IF (IGRENTRY.NE.1) DERIVDULOK = DERIVDULOK + 
     $	      U * ( SDQ * ( DZOMEG*T1Q - T2Q + URBURBINV*T3Q ) - 
     $	      DZOMEG*CQ * ( T2Q + URBURBINV*T3Q ) )            ! final
C	      propagation with PTdeltaQ for yta
	      T1Q =   CQ*YT1 - STQ*YT4
	      T3Q =   CQ*YT3 + STQ*YT6
	      YT4 =  SDQ*YT1 +  CQ*YT4
	      YT6 = -SDQ*YT3 +  CQ*YT6
	      YT1 = T1Q
	      YT3 = T3Q
	    END IF
C	    derivative-term computation  ( ytb(id)*Cinvdelta'*yta )
C	    propagation with Cinvdelta for yta (except for "rh1 vector-mult")
	    T2Q = YT2
	    YT2 = YT2+YT2           ! temporary doubling
C	    ! "rh2*ytb1" and "rh1*ytb6" could be stored from harkderivb
	    IF (IGRENTRY.NE.1) DERIVDULOK = DERIVDULOK + 
     $	    RH2*YTB1(ID)*(YT2-2*U*YT6) + 
     $	    2*YTB2(ID)*(-G0*YT1-2*G0U*YT2-G3*YT6) + 
     $	    RH1*YTB6(ID)*(-2*G0G0U*YT1+G0G3*YT2-G4*YT6)          ! final
	    YT1 = YT1 + U*(YT2-U*YT6)
	    T6Q = U*YT6
	    T1Q = G0U*YT1 + T6Q
	    YT6 = YT6 + G0U*(YT2-T6Q-T1Q)
	    YT2 = T2Q - T1Q
	    YT1 = RH2*YT1
	    YT6 = RH1*YT6
C	    set the appropriate exponentials
	    IF (EXPTYP) THEN
	      TLN = TLN + TLNREF - DLOG4
	    ELSE
	      TLN = TLN + TLNREF
	    END IF
	    DLNDOLOK = TLNB(ID) + TLN
	    DLNDULOK = DLNDOLOK

	  ELSE

	    stop 51267                  !!! not implemented yet !!!

	  END IF

	  TMAX = DMAX1 ( DABS(DREAL(YT1)), DABS(DIMAG(YT1)), 
     $	  DABS(DREAL(YT2)), DABS(DIMAG(YT2)), 
     $	  DABS(DREAL(YT3)), DABS(DIMAG(YT3)), 
     $	  DABS(DREAL(YT4)), DABS(DIMAG(YT4)), 
     $	  DABS(DREAL(YT6)), DABS(DIMAG(YT6)) )
	  IF ( TMAX.LT.tspmaxlim1 .OR. tspmaxlim2.LT.TMAX ) THEN
	    TMAX = 1/TMAX   ! if tmax was 0.0d0, due to num. canc.
	    TLN = TLN - DLOG(TMAX)
	    YT1 = TMAX*YT1
	    YT2 = TMAX*YT2
	    YT3 = TMAX*YT3
	    YT4 = TMAX*YT4
	    YT6 = TMAX*YT6
	  END IF

	  IF (IGRENTRY.NE.2) THEN
	    CALL ADDEC ( DLNDO,DERIVDO, DLNDOLOK,DERIVDOLOK, DLNDO,DERIVDO )
	    CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDO,DERIVDO )
	  END IF

	  IF (IGRENTRY.NE.1) THEN
	    CALL ADDEC ( DLNDU,DERIVDU, DLNDULOK,DERIVDULOK, DLNDU,DERIVDU )
	    CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )
	  END IF

	END DO

	IF (IHARK0.EQ.1) THEN

	  ID = NLNER+1
	  RH1 = RHO(ID)                        ! local rho(id)
	  RH2 = 1/RH1                          ! could be stored
	  G0 = TWOVB2(ID)                      ! local twovb2(id)
	  G0U = G0*U
	  G1 = G0*U2
	  G2 = 1-G1
	  G3 = 1-3*G1
	  G4 = 4*G0U*G2
	  G0G0U = G0*G0U
	  G0G3 = G0*G3
C	  derivative-term computation  ( ytbn(id)*(rho*Cdelta)'*yta )
C	  recognizing that "ytbn1(id)" is one and "ytbn2(id)" is zero in 
C	  this case
C	  propagation with "rho*Cdelta" for yta
	  YT1 = RH1*YT1
	  YT6 = RH2*YT6
	  T2Q = YT2
	  YT2 = YT2+YT2         ! temporary doubling
	  IF (IGRENTRY.NE.1) THEN
	    DLNDULOK = TLNB(ID) + TLN
	    DERIVDULOK = (-G4*YT1-G3*YT2-2*U*YT6) + YTB6N(ID)*
     $	    (-2*G0G0U*YT1-G0*YT2)     ! + 2*ytb2n(id)*(g0g3*yt1-2*g0u*yt2+yt6)
	    CALL ADDEC ( DLNDU,DERIVDU, DLNDULOK,DERIVDULOK, DLNDU,DERIVDU )
	    CALL SCALEEC ( tspmaxlim1,tspmaxlim2, DLNDU,DERIVDU )
	  END IF
	  YT6 = YT6 - G0U*(YT2+G0U*YT1)
	  T1Q = G0U*YT1
	  T6Q = U*YT6 + T1Q
	  YT1 = YT1 - U*(YT2+T6Q+T1Q)
	  YT2 = T2Q + T6Q

	END IF

	RETURN
	END