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

	USE INMODELa
	USE GRLOCALa

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       nol <= noldim

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

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

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

	INTEGER ID
	REAL*16 DLNDOLOK, DLNDULOK, PLNREF, PI
	COMPLEX*32 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, QCMPLX
	COMPLEX*32 CRHO
	LOGICAL EXPTYP

	DATA PI /3.14159265358979323846264338327950280Q0/

	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 = QCMPLX(AIMAG(URA),-REAL(URA))    ! -i*ura
	    END IF
	    CRHO = RHO(ID) * UA2(ID) * REAL(1/UA2(ID))
	    YP2 = YP2/CRHO  ! 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)*(REAL(SDP)*CP)*YP1 + 
     $	      YTB3L(ID)*CONJG(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 = REAL(SDP)*CP * YP1
	      YP2 = CONJG(CP) * YP2
	      IF (IGRENTRY.NE.1) DERIVDULOK = (U/URANEGI) * 
     $	      ( DERIVDULOK + YTB1(ID)*(YP1-YP2) )                ! final
	      WWLOK = YP1 - YP2
	      YP2 = CRHO * ( YP1 + YP2 )
	      YP1 = URANEGI*WWLOK
	    ELSE
C	      now derivative-term computation 
	      YTB3LOK = YTB3(ID)*CRHO    ! 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 = CRHO * ( -SDP*YP1 + CP*YP2 )
	      YP1 = WWLOK
	    END IF

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

	    STOP '%%% liqderiva, SHf'

	  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
