C	routine called from grderiv

C	Sven Ivansson 1993




	SUBROUTINE LIQDERIVB ( NLNER,NLUPP, OMEG,OMEG2,U,U2,RK,RK2 )
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
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2

	INTEGER ID, ISLASK
	REAL*8 TLN, TLNREF, PI
	COMPLEX*16 YT1,YT3, URA, URANEGI, 
     $	CP,CPX, STP,SDP, WWLOK, 
     $	CFAKT,BCFAKT,BCFAKTINV, YA,YB
	LOGICAL EXPTYP

	DATA PI /3.141592653589793D0/

	ID = NLNER
	TLN = TLNB(ID)
	YT1 = YTB1(ID)
	YT3 = YTB3(ID)
	CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLN,YT1,YT3 )

	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 LIQPART ( OMEG,OMEG2,U,U2, UA2(ID), DZ(ID), 
     $	      URA, EXPTYP, TLNREF, CP, STP,SDP )
	      URA_GR(ID) = URA
	      EXPTYP_GR(ID) = EXPTYP
	      TLNREF_GR(ID) = TLNREF
	      CP_GR(ID) = CP
	      STP_GR(ID) = STP
	      SDP_GR(ID) = SDP
	      IF (EXPTYP) URANEGI = DCMPLX(DIMAG(URA),-DREAL(URA))    ! -i*ura
	    ELSE
	      URA_GR(ID) = URA_GR(ID+1)
	      EXPTYP_GR(ID) = EXPTYP_GR(ID+1)
	      TLNREF_GR(ID) = TLNREF_GR(ID+1)
	      CP_GR(ID) = CP_GR(ID+1)
	      STP_GR(ID) = STP_GR(ID+1)
	      SDP_GR(ID) = SDP_GR(ID+1)
	    END IF
	    YT3 = YT3*RHO(ID)    ! already here
	    IF (EXPTYP) THEN
C	      now the first propagation step
	      WWLOK = YT1*URANEGI
	      YT1 = WWLOK + YT3
	      YT3 = -WWLOK + YT3
	      YTB1L(ID) = YT1
	      YTB3L(ID) = YT3
C	      now the main propagation step, freq-dependent
	      TLN = TLN + TLNREF - DLOG2       ! -dlog2 from "last step"
	      YT1 = YT1 * DREAL(SDP)*CP
	      YT3 = YT3 * DCONJG(CP)
	      YTB1M(ID) = YT1
	      YTB3M(ID) = YT3
C	      now the last propagation step
	      WWLOK = YT1 - YT3
	      YT3 = ( YT1 + YT3 ) / RHO(ID)      ! "1/rho(id)" could be stored
	      YT1 = WWLOK / URANEGI
	    ELSE
C	      now the main and only propagation step, freq-dependent
	      TLN = TLN + TLNREF
	      WWLOK = YT1*CP - YT3*SDP
	      YT3 = ( YT1*STP + YT3*CP ) / RHO(ID)  ! "1/rho" could be stored
	      YT1 = WWLOK
	    END IF

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

	    IF ( IFSMTYP(ID).NE.-2 .AND. IFSMTYP(ID).NE.(-3) ) THEN
C	      warning: not allowed to skip for "first id" in non-empty loop
	      ISLASK = 0
	    ELSE
	      ISLASK = -1
	    END IF
	    CALL AIRYPART ( OMEG,OMEG2,U,U2,RK,RK2, 
     $	    UA2(ID),UB2(ID), DZ(ID), ISLASK, 
     $	    CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	    TLNREF, CP,CPX, STP,SDP )
	    CFAKT_GR(ID) = CFAKT
	    BCFAKT_GR(ID) = BCFAKT
	    BCFAKTINV_GR(ID) = BCFAKTINV
	    YA_GR(ID) = YA
	    YB_GR(ID) = YB
	    TLNREF_GR(ID) = TLNREF
	    CP_GR(ID) = CP
	    CPX_GR(ID) = CPX
	    STP_GR(ID) = STP
	    SDP_GR(ID) = SDP
C	    now propagation
	    TLN = TLN + TLNREF
	    YT3 = YT3*RHO(ID)
	    WWLOK = PI * ( YT1*CP - YT3*SDP )
	    YT3 = (PI/RHO(ID)) * ( YT1*STP + YT3*CPX )  !"pi/rho" could stored
	    YT1 = WWLOK

	  ELSE

	    stop 51267                  !!! not implemented yet !!!

	  END IF

	  CALL SCALE2EC ( tspmaxlim1,tspmaxlim2, TLN,YT1,YT3 )
	  TLNB(ID-1) = TLN
	  YTB1(ID-1) = YT1
	  YTB3(ID-1) = YT3

	END DO

	RETURN
	END