C	routine called from grderiv

C	Sven Ivansson 1993




	SUBROUTINE LIQDERIVA ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2, 
     $	PLN, YP1,YP2, 
     $	DLNDO,DERIVDO, DLNDU,DERIVDU )
C	"downwards"
C	"gr input" in "common /grlocal/"
C	input/output: pln,yp1,yp2, dlndo,derivdo, dlndu,derivdu
C	note: liqderivb 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
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2

	REAL*8 PLN, DLNDO, DLNDU
	COMPLEX*16 YP1,YP2, DERIVDO, DERIVDU

	INTEGER ID
	REAL*8 DLNDOLOK, DLNDULOK, PLNREF, PI
	COMPLEX*16 DERIVDOLOK, DERIVDULOK, YTB3LOK, 
     $	URAURA, URA, URANEGI, 
     $	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/

	ENTRY LIQDERIVADO ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2, 
     $	PLN, YP1,YP2, DLNDO,DERIVDO )

	ENTRY LIQDERIVADU ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2, 
     $	PLN, YP1,YP2, DLNDU,DERIVDU )

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

	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
C	      stored "ura" (freq-indep) and "call liqpart" (freq-dep)
	      URA = URA_GR(ID)
	      EXPTYP = EXPTYP_GR(ID)
	      PLNREF= TLNREF_GR(ID)
	      CP = CP_GR(ID)
	      STP = STP_GR(ID)
	      SDP = SDP_GR(ID)
	      IF (EXPTYP) URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura
	    END IF
	    YP2 = YP2 / RHO(ID)    !already here  !"1/rho" could be stored
	    IF (EXPTYP) THEN
C	      now derivative-term computation 
C	      and column-vector propagation for yp ("from right to left")
	      IF (IGRENTRY.NE.1) DERIVDULOK = (YTB1M(ID)-YTB3M(ID))*YP1/URAURA
	      WWLOK = YP1 / URANEGI
	      YP1 = WWLOK + YP2
	      YP2 = -WWLOK + YP2
	      DERIVDOLOK = -YTB1L(ID)*(DREAL(SDP)*CP)*YP1 + 
     $	      YTB3L(ID)*DCONJG(CP)*YP2
	      IF (IGRENTRY.NE.1) DERIVDULOK = DERIVDULOK + 
     $	      OMEG*DZ(ID)*DERIVDOLOK
	      IF (IGRENTRY.NE.2) DERIVDOLOK = DZ(ID)*URANEGI * 
     $	      DERIVDOLOK                                         ! final
	      PLN = PLN + PLNREF - DLOG2
	      YP1 = DREAL(SDP)*CP * YP1
	      YP2 = DCONJG(CP) * YP2
	      IF (IGRENTRY.NE.1) DERIVDULOK = (U/URANEGI) * 
     $	      ( DERIVDULOK + YTB1(ID)*(YP1-YP2) )                ! final
	      WWLOK = YP1 - YP2
	      YP2 = RHO(ID) * ( YP1 + YP2 )
	      YP1 = URANEGI*WWLOK
	    ELSE
C	      now derivative-term computation 
	      YTB3LOK = YTB3(ID)*RHO(ID)    ! could be stored from liqderivb
	      DCP = DZ(ID)*CP
	      OMEGDCP = OMEG*DCP
	      DSTP = DZ(ID)*STP
	      OMEGDSDP = OMEG*DZ(ID)*SDP       ! "omeg*dz(id)" could be stored
	      IF (IGRENTRY.NE.2) DERIVDOLOK = YTB1(ID) * ( -DSTP*YP1 + 
     $	      URAURA*DCP*YP2 ) - YTB3LOK * ( DCP*YP1 + DSTP*YP2 )      ! final
	      IF (IGRENTRY.NE.1) DERIVDULOK = U * ( 
     $	      YTB1(ID) * ( OMEGDSDP*YP1 - (SDP+OMEGDCP)*YP2 ) + 
     $	      YTB3LOK * ( (OMEGDCP-SDP)*YP1/URAURA + OMEGDSDP*YP2 ) )  ! final
C	      and column-vector propagation for yp ("from right to left")
	      PLN = PLN + PLNREF
	      WWLOK = CP*YP1 + STP*YP2
	      YP2 = RHO(ID) * ( -SDP*YP1 + CP*YP2 )
	      YP1 = WWLOK
	    END IF

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

C	    stored "call airypart"
	    CFAKT = CFAKT_GR(ID)
	    BCFAKT = BCFAKT_GR(ID)
	    BCFAKTINV = BCFAKTINV_GR(ID)
	    YA = YA_GR(ID)
	    YB = YB_GR(ID)
	    PLNREF = TLNREF_GR(ID)
	    CP = CP_GR(ID)
	    CPX = CPX_GR(ID)
	    STP = STP_GR(ID)
	    SDP = SDP_GR(ID)
	    PIOMEG2UTWOCFAKT = PIOMEG2UTWO*CFAKT
C	    now derivative-term computation 
	    YTB3LOK = YTB3(ID)*RHO(ID)      ! could be stored from liqderivb
	    YACP = YA*CP
	    YBCPX = YB*CPX
	    BCFAKTINVSTP = BCFAKTINV*STP
	    BCFAKTSDP = BCFAKT*SDP
	    BCFAKTSDPYA = BCFAKTSDP*YA
	    BCFAKTSDPYB = BCFAKTSDP*YB
	    YP2 = YP2/RHO(ID)   ! already here  ! "1/rho" could be stored
	    IF (IGRENTRY.NE.2) DERIVDOLOK = PIOMEGTHREEINV * ( 
     $	    YTB1(ID) * ( (-BCFAKTINVSTP*YA-BCFAKTSDPYB*YB)*2*YP1 + 
     $	    (-STP+2*BCFAKT*(-YACP*YA+YBCPX*YB))*YP2 ) + 
     $	    YTB3LOK * ( (-SDP-2*BCFAKTINV*(YA*CPX-YB*CP))*YP1 + 
     $	    (BCFAKTSDPYA*YA+BCFAKTINVSTP*YB)*2*YP2 ) )      ! final
	    IF (IGRENTRY.NE.1) DERIVDULOK = PIOMEG2UTWOCFAKT * ( 
     $	    YTB1(ID) * ( (-BCFAKTINVSTP-BCFAKTSDPYB)*YP1 + 
     $	    (BCFAKT*(-YACP+YBCPX))*YP2 ) + 
     $	    YTB3LOK * ( -BCFAKTINV*(CPX-CP)*YP1 + 
     $	    (BCFAKTSDPYA+BCFAKTINVSTP)*YP2 ) )              ! final
C	    and column-vector propagation for yp ("from right to left")
	    PLN = PLN + PLNREF
	    WWLOK = PI * ( CP*YP1 + STP*YP2 )
	    YP2 = (PI*RHO(ID)) * ( -SDP*YP1 + CPX*YP2 )  !"pi*rho" could be sto
	    YP1 = WWLOK

	  ELSE

	    stop 51267                   !!! not implemented yet         !!!

	  END IF

	  DLNDOLOK = PLN + TLNB(ID)
	  DLNDULOK = DLNDOLOK

	  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

	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, PLN,YP1,YP2 )

	END DO

	RETURN
	END