C	extrapolation for a "new line" in the extrapolation table, 
C	developed particularly for adintxfi.f (starting 
C	from Ilkka Karasalos "cdadbs.for")
C	ch. 6 of the numerical ODE book by Gear from 1971 is relevant

C	Sven Ivansson 1990



	SUBROUTINE XPOLBUL ( NHYD,LI, IXTYP, HSHARE,EPS, 
     $	HQI,FIHQI,QHQI, JPX,JCOLX, 
     $	IA_TABNEW,IA_VALTERM,IA_TAB,IA_YTAB, CWORK,OKSUB )
C	notation as in adintxfi.for
C	  note that the "new line" is the "i:th line"
C	  hqi() and jcolx will "depend on i"
C	input:      nhyd,li
C	            ixtyp  < "1,-1,0" for "polynomial" ; 
C	                     "2,-2"   for "rational" ; 
C	                     "3,-3"   for "inverse-polynomial" ; 
C	                     "4 & -4" for "rational without x-terms" >
C	            hshare,eps(1:li)
C	            hqi(1:jcolx)  < only hqi(2:jcolx) needed >
C	            fihqi(1:jcolx)  < only fihqi(1:jcolx) needed >
C	                          < only needed when ixtyp=4,-4 , except 
C	                            fihqi(1:2) when "Cxut" is active >
C	            qhqi(1:3)     < only needed when ixtyp=-4 >
C	            jpx           < indicates whether "all extrapolated 
C	                            columns" will be "finally possible" in 
C	                            the extrapolation; if so "jpx=0" else 
C	                            "jpx=1 (final column overrides)" >
C	            jcolx         < extrapolation columns: 2,...,jcolx
C	                            note: "initial value" from column 1 >
C	            ia_tabnew,ia_valterm,ia_tab,ia_ytab  are integer 
C	                  functions providing suitable addresses 
C	                  within cwork(*)
C	input/output: cwork(*)
C	                 input; 
C	                    "tabnew" is "initial value for line i (from 
C	                                                       column 1)"
C	                    "tab(1:jcolx-1)" is for "line i-1"
C	                    "ytab(2:jcolx-1)"  is for "line i-1"
C	                          "ytab(1)"  is for "line i-1" (artificial)
C	                 output; 
C	                    "tabnew" is "dcmplx(delmin,dble(jdel))"     <OBS>
C	                    "valterm" is "estimated value from line i"
C	                    "tab(1:jcolx)" is for "line i"
C	                    "ytab(2:jcolx)"  is for "line i"
C	                           "ytab(2)" is for "line i" (artificial)
C	              oksub
C	                 "true on exit" for "estimate ok"

	IMPLICIT NONE

	INTEGER NHYD,LI, IXTYP, JPX,JCOLX, 
     $	IA_TABNEW,IA_VALTERM,IA_TAB,IA_YTAB
	REAL*8 HSHARE,EPS(*), HQI(*),FIHQI(*),QHQI(*)

	COMPLEX*16 CWORK(*)
	LOGICAL OKSUB

	INTEGER IH,IR, IHR, J, JDEL
	REAL*8 EPSLOK, DELMIN, SLASK
	COMPLEX*16 TABNW,USLASK,VSLASK,CSLASK, BSLASK,B1SLASK, 
     $	YTABNW,YUSLASK,YVSLASK,YCSLASK, FITABNW,FIUSLASK, FITSLASK, 
     $	Q1LOK,Q2LOK, QLOK, W1LOK,W2LOK, VSLASKOLD,YVALOLD,ALFAOLD,
     $	ALFAVAL, ZMOD
	REAL*8 DENOMIN, SAFMINFACT, DLAMCH
	PARAMETER ( SAFMINFACT = 3.0D0 )   ! empirically ok value
	SAVE DENOMIN
	LOGICAL SETFIRST
	DATA SETFIRST /.TRUE./
Cval	logical oksubin
Cerr	integer igang
Cerr	real*8 errlok(100), errsum(100)
Cerr	complex*16 valsum(100)
Cerr	save igang, errsum,valsum
Cerr	data igang /0/

	EXTERNAL IA_TABNEW,IA_VALTERM,IA_TAB,IA_YTAB

Cval	oksubin = oksub

Cerr	if (igang.eq.0) then
Cerr	  if (nhyd*li.gt.100) stop 6061
Cerr	  do ihr=1,nhyd*li
Cerr	    valsum(ihr) = (0.0d0,0.0d0)
Cerr	    errsum(ihr) = 0.0d0
Cerr	  end do
Cerr	  igang = 1
Cerr	end if

	IF (SETFIRST) THEN
	  DENOMIN = SAFMINFACT * DLAMCH('S')
	  SETFIRST = .FALSE.
	END IF

	J = IABS(IXTYP)
	GOTO (1,2,3,4,5) J

 1	CONTINUE    ! ixtyp=1,-1,0
	IF (J.GT.1) STOP 1561
C	polynomial-function extrapolation in this case
	DO IR=1,LI
	  EPSLOK = HSHARE * EPS(IR)
	  IHR = IR
	  DO IH=1,NHYD
C	    extrapolate to i:th row (for the ihr:th integral)
C	    we use now  "tab(j,i)" := "val(j,i)" - "val(j-1,i)"  
C	     (obvious abuse of notation,  "val(0,i)" = 0 )
C	    val(j,i) :=  val(j-1,i)  +  del(j-1,i) / ( hqi(j,i) - 1 )   where 
C	     del(j,i) = cslask(j,i)-tab(j,i-1) = val(j,i)-val(j,i-1)
C	     cslask(j,i) = val(j,i) - val(j-1,i-1)
	    TABNW = CWORK(IA_TABNEW(IHR))  ! val(1,i)
Clog	    write(99,*) 1,tabnw,' ! ',sngl(epslok)
Cxut	    if (ihr.eq.1) then
Cxut	      if (jcolx.eq.2) then
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*) -1,'  ! NEW, first from first x row'
Cxut	        write(99,*) 1/(fihqi(2)-fihqi(1)),
Cxut     $	        cwork(ia_tab(ihr,1)),' &'
Cxut	      end if
Cxut	      write(99,*)
Cxut	      write(99,*) jcolx,' nr of columns'
Cxut	      write(99,*) 1/fihqi(1),tabnw,' &'
Cxut	    end if
	    CSLASK = TABNW       ! val(1,i) ( - val(0,i-1) )
	    USLASK = TABNW       ! val(1,i) ( - val(0,i) )
	    VSLASK = CWORK(IA_TAB(IHR,1))
C	                      ! val(1,i-1) ( - val(0,i-1) )
	    CWORK(IA_TAB(IHR,1)) = USLASK
C	                      ! val(1,i) ( - val(0,i) )
	    USLASK = (CSLASK-VSLASK)/(HQI(2)-1) ! val(2,i) - val(1,i)
	    CSLASK = HQI(2)*USLASK ! val(2,i) - val(1,i-1)
	    TABNW = TABNW + USLASK   ! val(2,i)
	    IF (IXTYP.EQ.1) THEN
	      DELMIN = DMAX1 ( DABS(DREAL(CSLASK)), 
     $	      DABS(DIMAG(CSLASK)) )        ! cdabs ( val(2,i) - val(1,i-1) )
	    ELSE IF (IXTYP.EQ.-1) THEN
	      DELMIN = DMAX1 ( DABS(DREAL(USLASK)), 
     $	      DABS(DIMAG(USLASK)) )        ! cdabs ( val(2,i) - val(1,i) )
	    ELSE
	      DELMIN = DMAX1 ( DABS(DREAL(CSLASK-USLASK)), 
     $	      DABS(DIMAG(CSLASK-USLASK)) ) ! cdabs ( val(1,i) - val(1,i-1) )
	    END IF
Clog	    write(99,*) 2,uslask,sngl(delmin)
Cxut	    if (ihr.eq.1) write(99,*) 1/fihqi(1),tabnw
	    JDEL = 2
	    CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	    errlok(ihr) = delmin
	    DO J = 3, JCOLX
	      VSLASK = CWORK(IA_TAB(IHR,J-1)) ! val(j-1,i-1) - val(j-2,i-1)
	      CWORK(IA_TAB(IHR,J-1)) = USLASK ! val(j-1,i) - val(j-2,i)
	      USLASK = (CSLASK-VSLASK)/(HQI(J)-1) ! val(j,i) - val(j-1,i)
	      CSLASK = HQI(J)*USLASK  ! val(j,i) - val(j-1,i-1)
	      TABNW = TABNW + USLASK     ! val(j,i)
	      IF (IXTYP.EQ.1) THEN
	        SLASK = DMAX1 ( DABS(DREAL(CSLASK)), 
     $	        DABS(DIMAG(CSLASK)) )        ! cdabs(val(j,i)-val(j-1,i-1))
	      ELSE IF (IXTYP.EQ.-1) THEN
	        SLASK = DMAX1 ( DABS(DREAL(USLASK)), 
     $	        DABS(DIMAG(USLASK)) )        ! cdabs(val(j,i)-val(j-1,i))
	      ELSE
	        SLASK = DMAX1 ( DABS(DREAL(CSLASK-USLASK)), 
     $	        DABS(DIMAG(CSLASK-USLASK)) ) ! cdabs(val(j-1,i)-val(j-1,i-1))
	      END IF
Clog	      write(99,*) j,uslask,sngl(slask)
Cxut	      if (ihr.eq.1) write(99,*) 1/fihqi(1),tabnw
	      IF (SLASK.LT.DELMIN) THEN
	        DELMIN = SLASK
	        JDEL = J
	        CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	        errlok(ihr) = slask
	      END IF   ! picking the most accurate value
	    END DO
	    IF (JPX.NE.0) THEN
	      IF (JCOLX.GE.3) THEN
	        DELMIN = SLASK
	        JDEL = JCOLX
	        CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	        errlok(ihr) = slask
	      END IF   ! last column "overrides"
	    END IF
	    CWORK(IA_TAB(IHR,JCOLX)) = USLASK ! val(jcolx,i) - val(jcolx-1,i)
C	    check if sufficient accuracy in this subinterval
	    IF ( DELMIN .GT. EPSLOK ) OKSUB = .FALSE.
	    CWORK(IA_TABNEW(IHR)) = DCMPLX(DELMIN,DBLE(JDEL))
	    IHR = IHR + LI
	  END DO
	END DO
	GOTO 999

 2	CONTINUE    ! ixtyp=2,-2
C	rational-function Bulirsch-Stoer extrapolation in this case
	DO IR=1,LI
	  EPSLOK = HSHARE * EPS(IR)
	  IHR = IR
	  DO IH=1,NHYD
C	    extrapolate to i:th row (for the ihr:th integral)
C	    we use now  "tab(j,i)" := "val(j,i)" - "val(j-1,i)"  
C	     (obvious abuse of notation,  "val(0,i)" = 0 )
C	    val(j,i) :=  val(j-1,i)  +  del(j-1,i) * cslask(j-1,i) / 
C	           ( hqi(j,i)*tab(j-1,i-1) - cslask(j-1,i) )   where 
C	     del(j,i) = cslask(j,i)-tab(j,i-1) = val(j,i)-val(j,i-1)
C	     cslask(j,i) = val(j,i) - val(j-1,i-1)
	    TABNW = CWORK(IA_TABNEW(IHR))  ! val(1,i)
Clog	    write(99,*) 1,tabnw,' ! ',sngl(epslok)
Cxut	    if (ihr.eq.1) then
Cxut	      if (jcolx.eq.2) then
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*) -1,'  ! NEW, first from first x row'
Cxut	        write(99,*) 1/(fihqi(2)-fihqi(1)),
Cxut     $	        cwork(ia_tab(ihr,1)),' &'
Cxut	      end if
Cxut	      write(99,*)
Cxut	      write(99,*) jcolx,' nr of columns'
Cxut	      write(99,*) 1/fihqi(1),tabnw,' &'
Cxut	    end if
	    CSLASK = TABNW       ! val(1,i) ( - val(0,i-1) )
	    USLASK = TABNW       ! val(1,i) ( - val(0,i) )
	    VSLASK = CWORK(IA_TAB(IHR,1))
C	                      ! val(1,i-1) ( - val(0,i-1) )
	    CWORK(IA_TAB(IHR,1)) = USLASK
C	                      ! val(1,i) ( - val(0,i) )
	    B1SLASK = VSLASK*HQI(2)
	    BSLASK = B1SLASK - CSLASK ! cslask=val(1,i)(-val(0,i-1))
	    IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $      DENOMIN ) THEN
	      BSLASK = (CSLASK-VSLASK)/BSLASK ! val(1,i) - val(1,i-1)
	      USLASK = CSLASK*BSLASK  ! val(2,i) - val(1,i)
	      CSLASK = B1SLASK*BSLASK ! val(2,i) - val(1,i-1)
	    ELSE
	      USLASK = VSLASK    ! cslask = val(2,i) - val(1,i-1)  ok
	    END IF
	    TABNW = TABNW + USLASK   ! val(2,i)
	    IF (IXTYP.EQ.2) THEN
	      DELMIN = DMAX1 ( DABS(DREAL(CSLASK)), 
     $	      DABS(DIMAG(CSLASK)) )    ! cdabs ( val(2,i) - val(1,i-1) )
	    ELSE
	      DELMIN = DMAX1 ( DABS(DREAL(USLASK)), 
     $	      DABS(DIMAG(USLASK)) )    ! cdabs ( val(2,i) - val(1,i) )
	    END IF
Clog	    write(99,*) 2,uslask,sngl(delmin)
Cxut	    if (ihr.eq.1) write(99,*) 1/fihqi(1),tabnw
	    JDEL = 2
	    CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	    errlok(ihr) = delmin
	    DO J = 3, JCOLX
	      VSLASK = CWORK(IA_TAB(IHR,J-1)) ! val(j-1,i-1) - val(j-2,i-1)
	      CWORK(IA_TAB(IHR,J-1)) = USLASK ! val(j-1,i) - val(j-2,i)
	      B1SLASK = VSLASK*HQI(J)
	      BSLASK = B1SLASK-CSLASK      ! cslask=val(j-1,i)-val(j-2,i-1)
	      IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $        DENOMIN ) THEN
	        BSLASK = (CSLASK-VSLASK)/BSLASK ! val(j-1,i) - val(j-1,i-1)
	        USLASK = CSLASK*BSLASK   ! val(j,i) - val(j-1,i)
	        CSLASK = B1SLASK*BSLASK  ! val(j,i) - val(j-1,i-1)
	      ELSE
	        USLASK = VSLASK  ! cslask = val(j,i) - val(j-1,i-1) ok
	      END IF
	      TABNW = TABNW + USLASK     ! val(j,i)
	      IF (IXTYP.EQ.2) THEN
	        SLASK = DMAX1 ( DABS(DREAL(CSLASK)), 
     $	        DABS(DIMAG(CSLASK)) )    ! cdabs ( val(j,i) - val(j-1,i-1) )
	      ELSE
	        SLASK = DMAX1 ( DABS(DREAL(USLASK)), 
     $	        DABS(DIMAG(USLASK)) )    ! cdabs ( val(j,i) - val(j-1,i) )
	      END IF
Clog	      write(99,*) j,uslask,sngl(slask)
Cxut	      if (ihr.eq.1) write(99,*) 1/fihqi(1),tabnw
	      IF (SLASK.LT.DELMIN) THEN
	        DELMIN = SLASK
	        JDEL = J
	        CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	        errlok(ihr) = slask
	      END IF   ! picking the most accurate value
	    END DO
	    IF (JPX.NE.0) THEN
	      IF (JCOLX.GE.3) THEN
	        DELMIN = SLASK
	        JDEL = JCOLX
	        CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	        errlok(ihr) = slask
	      END IF   ! last column "overrides"
	    END IF
	    CWORK(IA_TAB(IHR,JCOLX)) = USLASK ! val(jcolx,i) - val(jcolx-1,i)
C	    check if sufficient accuracy in this subinterval
	    IF ( DELMIN .GT. EPSLOK ) OKSUB = .FALSE.
	    CWORK(IA_TABNEW(IHR)) = DCMPLX(DELMIN,DBLE(JDEL))
	    IHR = IHR + LI
	  END DO
	END DO
	GOTO 999

 3	CONTINUE    ! ixtyp=3,-3
C	inverse-polynomial function extrapolation in this case, implemented 
C	by modifying the "rational-function Bulirsch-Stoer extrapolation" 
C	(perhaps it would have been more natural to "modify the polynomial-
C	function extrapolation" using "inverse values")
	DO IR=1,LI
	  EPSLOK = HSHARE * EPS(IR)
	  IHR = IR
	  DO IH=1,NHYD
C	    extrapolate to i:th row (for the ihr:th integral)
C	    we use now  "tab(j,i)" := "val(j,i)" - "val(j-1,i)"  
C	     (obvious abuse of notation,  "val(0,i)" = 0 )
C	    val(j,i) :=  val(j-1,i)  +  del(j-1,i) * val(j-1,i) / 
C	           ( hqi(j,i)*val(j-1,i-1) - val(j-1,i) )   where 
C	     del(j,i) = cslask(j,i)-tab(j,i-1) = val(j,i)-val(j,i-1)
C	     cslask(j,i) = val(j,i) - val(j-1,i-1)
	    TABNW = CWORK(IA_TABNEW(IHR))  ! val(1,i)
Clog	    write(99,*) 1,tabnw,' ! ',sngl(epslok)
Cxut	    if (ihr.eq.1) then
Cxut	      if (jcolx.eq.2) then
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*) -1,'  ! NEW, first from first x row'
Cxut	        write(99,*) 1/(fihqi(2)-fihqi(1)),
Cxut     $	        cwork(ia_tab(ihr,1)),' &'
Cxut	      end if
Cxut	      write(99,*)
Cxut	      write(99,*) jcolx,' nr of columns'
Cxut	      write(99,*) 1/fihqi(1),tabnw,' &'
Cxut	    end if
	    CSLASK = TABNW       ! val(1,i) ( - val(0,i-1) )
	    USLASK = TABNW       ! val(1,i) ( - val(0,i) )
	    VSLASK = CWORK(IA_TAB(IHR,1))
C	                      ! val(1,i-1) ( - val(0,i-1) )
	    CWORK(IA_TAB(IHR,1)) = USLASK
C	                      ! val(1,i) ( - val(0,i) )
	    B1SLASK = VSLASK*HQI(2)    ! vslask=val(1,i-1) here
	    BSLASK = B1SLASK - TABNW     ! tabnw=val(1,i) here
	    IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $      DENOMIN ) THEN
	      BSLASK = (CSLASK-VSLASK)/BSLASK ! val(1,i) - val(1,i-1)
	      USLASK = TABNW*BSLASK  ! val(2,i) - val(1,i)
	      CSLASK = B1SLASK*BSLASK ! val(2,i) - val(1,i-1)
	    ELSE
	      USLASK = VSLASK    ! cslask = val(2,i) - val(1,i-1)  ok
	    END IF
	    TABNW = TABNW + USLASK   ! val(2,i)
	    IF (IXTYP.EQ.2) THEN
	      DELMIN = DMAX1 ( DABS(DREAL(CSLASK)), 
     $	      DABS(DIMAG(CSLASK)) )    ! cdabs ( val(2,i) - val(1,i-1) )
	    ELSE
	      DELMIN = DMAX1 ( DABS(DREAL(USLASK)), 
     $	      DABS(DIMAG(USLASK)) )    ! cdabs ( val(2,i) - val(1,i) )
	    END IF
Clog	    write(99,*) 2,uslask,sngl(delmin)
Cxut	    if (ihr.eq.1) write(99,*) 1/fihqi(1),tabnw
	    JDEL = 2
	    CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	    errlok(ihr) = delmin
	    DO J = 3, JCOLX
	      VSLASK = CWORK(IA_TAB(IHR,J-1)) ! val(j-1,i-1) - val(j-2,i-1)
	      CWORK(IA_TAB(IHR,J-1)) = USLASK ! val(j-1,i) - val(j-2,i)
	      W1LOK = CSLASK-VSLASK   ! val(j-1,i) - val(j-1,i-1)
	      B1SLASK = (TABNW-W1LOK)*HQI(J)  ! val(j-1,i-1)*hqi(j)
	      BSLASK = B1SLASK-TABNW      ! tabnw=val(j-1,i)
	      IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $        DENOMIN ) THEN
	        BSLASK = W1LOK/BSLASK
	        USLASK = TABNW*BSLASK   ! val(j,i) - val(j-1,i)
	        CSLASK = B1SLASK*BSLASK  ! val(j,i) - val(j-1,i-1)
	      ELSE
	        USLASK = VSLASK  ! cslask = val(j,i) - val(j-1,i-1) ok
	      END IF
	      TABNW = TABNW + USLASK     ! val(j,i)
	      IF (IXTYP.EQ.2) THEN
	        SLASK = DMAX1 ( DABS(DREAL(CSLASK)), 
     $	        DABS(DIMAG(CSLASK)) )    ! cdabs ( val(j,i) - val(j-1,i-1) )
	      ELSE
	        SLASK = DMAX1 ( DABS(DREAL(USLASK)), 
     $	        DABS(DIMAG(USLASK)) )    ! cdabs ( val(j,i) - val(j-1,i) )
	      END IF
Clog	      write(99,*) j,uslask,sngl(slask)
Cxut	      if (ihr.eq.1) write(99,*) 1/fihqi(1),tabnw
	      IF (SLASK.LT.DELMIN) THEN
	        DELMIN = SLASK
	        JDEL = J
	        CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	        errlok(ihr) = slask
	      END IF   ! picking the most accurate value
	    END DO
	    IF (JPX.NE.0) THEN
	      IF (JCOLX.GE.3) THEN
	        DELMIN = SLASK
	        JDEL = JCOLX
	        CWORK(IA_VALTERM(IHR)) = TABNW
Cerr	        errlok(ihr) = slask
	      END IF   ! last column "overrides"
	    END IF
	    CWORK(IA_TAB(IHR,JCOLX)) = USLASK ! val(jcolx,i) - val(jcolx-1,i)
C	    check if sufficient accuracy in this subinterval
	    IF ( DELMIN .GT. EPSLOK ) OKSUB = .FALSE.
	    CWORK(IA_TABNEW(IHR)) = DCMPLX(DELMIN,DBLE(JDEL))
	    IHR = IHR + LI
	  END DO
	END DO
	GOTO 999

 4	CONTINUE    ! ixtyp=4,-4

	IF (IXTYP.EQ.4) THEN

C	  rational-function Bulirsch-Stoer extrapolation but "without x terms" 
C	  in this case;  (my,ny) order: (0,0) (0,1) (1,1) (1,2) (2,2) ...
	  DO IR=1,LI
	    EPSLOK = HSHARE * EPS(IR)
	    IHR = IR
	    DO IH=1,NHYD
C	      extrapolate to i:th row (for the ihr:th integral)
C	      we use now  "tab(j,i)" := "val(j,i)" - "val(j-1,i)"  
C	       (obvious abuse of notation,  "val(0,i)" = 0 )
C	      we use now  "ytab(j,i)" := "yval(j,i)" - "yval(j-1,i)"  
C	       (obvious abuse of notation,  "yval(0,i)" = 0 )
C	      val(j,i) :=  val(j-1,i)  +  del(j-1,i) * qqq(j,i)  where 
C	       del(j,i) = cslask(j,i)-tab(j,i-1) = val(j,i)-val(j,i-1)
C	       cslask(j,i) = val(j,i) - val(j-1,i-1)
C	      yval(j,i) :=  yval(j-1,i)  +  ydel(j-1,i) * qqq(j,i)  where 
C	       ydel(j,i) = ycslask(j,i)-ytab(j,i-1) = yval(j,i)-yval(j,i-1)
C	       ycslask(j,i) = yval(j,i) - yval(j-1,i-1)
C	      fival(j,i) :=  fival(j-1,i)  +  fiqqq(j,i)/yval(j,i)
	      TABNW = CWORK(IA_TABNEW(IHR))  ! val(1,i)
	      CSLASK = TABNW       ! val(1,i) ( - val(0,i-1) )
	      USLASK = TABNW       ! val(1,i) ( - val(0,i) )
	      YTABNW = FIHQI(1)  ! yval(1,i)
	      YCSLASK = YTABNW       ! yval(1,i) ( - yval(0,i-1) )
	      YUSLASK = YTABNW       ! yval(1,i) ( - yval(0,i) )
	      FITABNW = TABNW      ! fival(1,i)
Clog	      write(99,*) 1,fitabnw,' ! ',sngl(epslok)
Cxut	      if (ihr.eq.1) then
Cxut	        if (jcolx.eq.2) then
Cxut	          write(99,*)
Cxut	          write(99,*)
Cxut	          write(99,*)
Cxut	          write(99,*) -1,'  ! NEW, first from first x row'
Cxut	          write(99,*) 1/(fihqi(2)-fihqi(1)),
Cxut     $	          cwork(ia_tab(ihr,1)),' &'
Cxut	        end if
Cxut	        write(99,*)
Cxut	        write(99,*) jcolx,' nr of columns'
Cxut	        write(99,*) 1/fihqi(1),fitabnw,' &'
Cxut	      end if
	      FITSLASK = (0.0D0,0.0D0)  ! fival(1,i) - val(1,i)
	      VSLASK = CWORK(IA_TAB(IHR,1))
C	                        ! val(1,i-1) ( - val(0,i-1) )
	      CWORK(IA_TAB(IHR,1)) = USLASK
C	                        ! val(1,i) ( - val(0,i) )
	      YVSLASK = ( FIHQI(2) - FIHQI(1) )
C	                        ! yval(1,i-1) ( - yval(0,i-1) )
	      B1SLASK = VSLASK*HQI(2)    ! val(j-1,i-1)*hqi(j) here
	      BSLASK = B1SLASK - TABNW ! tabnw=val(1,i)
	      IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $        DENOMIN ) THEN
	        Q1LOK = (CSLASK-VSLASK)/BSLASK ! val(1,i) - val(1,i-1)
	        Q2LOK = (YCSLASK-YVSLASK)/BSLASK ! yval(1,i)-yval(1,i-1)
	        USLASK = TABNW*Q1LOK  ! val(2,i) - val(1,i)
	        YUSLASK = TABNW*Q2LOK  ! yval(2,i) - yval(1,i)
	        CSLASK = B1SLASK*Q1LOK ! val(2,i) - val(1,i-1)
	        YCSLASK = B1SLASK*Q2LOK ! yval(2,i) - yval(1,i-1)
	      ELSE
	        USLASK = VSLASK    ! cslask = val(2,i) - val(1,i-1)  ok
	        YUSLASK = YVSLASK  ! ycslask = yval(2,i)-yval(1,i-1)  ok
	      END IF
	      Q2LOK = YTABNW + YUSLASK
	      IF ( DMAX1(DABS(DREAL(Q2LOK)),DABS(DIMAG(Q2LOK))) .GE. 
     $        DENOMIN ) THEN
	        FIUSLASK = ( (FIHQI(2)-YTABNW)*USLASK - 
     $	        YUSLASK*FITSLASK ) / Q2LOK    ! fival(2,i) - fival(1,i)
	      ELSE
	        FIUSLASK = USLASK     ! somewhat arbitrary choice
	      END IF
	      TABNW = TABNW + USLASK   ! val(2,i)
	      YTABNW = Q2LOK            ! yval(2,i)
	      FITABNW = FITABNW + FIUSLASK   ! fival(2,i)
	      FITSLASK = FITSLASK + ( FIUSLASK - USLASK )
C	                        ! fival(2,i) - val(2,i)
	      DELMIN = DMAX1 ( DABS(DREAL(FIUSLASK)), 
     $	      DABS(DIMAG(FIUSLASK)) )   ! cdabs ( fival(2,i)-fival(1,i) )
Clog	      write(99,*) 2,fiuslask,sngl(delmin)
Cxut	      if (ihr.eq.1) then
Cxut	        write(99,*) 1/fihqi(1),fitabnw
Cxut	        write(99,*) '       ',tabnw,'     ! tabnw'
Cxut	        write(99,*) '       ',ytabnw-fihqi(2),
Cxut     $	        '     ! beta'
Cxut	        write(99,*) '       ',ytabnw*fitabnw-
Cxut     $	        fihqi(2)*tabnw,'     ! alfa'
Cxut	      end if
	      JDEL = 2
	      CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	      errlok(ihr) = delmin
	      DO J = 3, JCOLX
	        VSLASK = CWORK(IA_TAB(IHR,J-1))
C	                      ! val(j-1,i-1) - val(j-2,i-1)
	        CWORK(IA_TAB(IHR,J-1)) = USLASK
C	                      ! val(j-1,i) - val(j-2,i)
	        YVSLASK = CWORK(IA_YTAB(IHR,J-1))
C	                      ! yval(j-1,i-1) - yval(j-2,i-1)
	        CWORK(IA_YTAB(IHR,J-1)) = YUSLASK
C	                      ! yval(j-1,i) - val(j-2,i)
	        B1SLASK = VSLASK*HQI(J)
	        BSLASK = B1SLASK - CSLASK ! cslask=val(j-1,i)-val(j-2,i-1)
	        IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $          DENOMIN ) THEN
	          Q1LOK = (CSLASK-VSLASK)/BSLASK ! val(j-1,i)-val(j-1,i-1)
	          Q2LOK = (YCSLASK-YVSLASK)/BSLASK!yval(j-1,i)-yval(j-1,i-1)
	          USLASK = CSLASK*Q1LOK  ! val(j,i) - val(j-1,i)
	          YUSLASK = CSLASK*Q2LOK  ! yval(j,i) - yval(j-1,i)
	          CSLASK = B1SLASK*Q1LOK ! val(j,i) - val(j-1,i-1)
	          YCSLASK = B1SLASK*Q2LOK ! yval(j,i) - yval(j-1,i-1)
	        ELSE
	          USLASK = VSLASK    ! cslask = val(j,i) - val(j-1,i-1)  ok
	          YUSLASK = YVSLASK  ! ycslask = yval(j,i)-yval(j-1,i-1) ok
	        END IF
	        Q2LOK = YTABNW + YUSLASK
	        IF ( DMAX1(DABS(DREAL(Q2LOK)),DABS(DIMAG(Q2LOK))) .GE. 
     $          DENOMIN ) THEN
	          FIUSLASK = ( (FIHQI(J)-YTABNW)*USLASK - 
     $	          YUSLASK*FITSLASK ) / Q2LOK    ! fival(j,i) - fival(j-1,i)
	        ELSE
	          FIUSLASK = USLASK     ! somewhat arbitrary choice
	        END IF
	        TABNW = TABNW + USLASK   ! val(j,i)
	        YTABNW = Q2LOK            ! yval(j,i)
	        FITABNW = FITABNW + FIUSLASK   ! fival(j,i)
	        FITSLASK = FITSLASK + ( FIUSLASK - USLASK )
C	                      ! fival(j,i) - val(j,i)
	        SLASK = DMAX1 ( DABS(DREAL(FIUSLASK)), 
     $	        DABS(DIMAG(FIUSLASK)) )  ! cdabs ( fival(j,i)-fival(j-1,i) )
Clog	        write(99,*) j,fiuslask,sngl(slask)
Cxut	        if (ihr.eq.1) then
Cxut	          write(99,*) 1/fihqi(1),fitabnw
Cxut	          write(99,*) '       ',tabnw,'     ! tabnw'
Cxut	          write(99,*) '       ',ytabnw-fihqi(j),
Cxut     $	          '     ! beta'
Cxut	          write(99,*) '       ',ytabnw*fitabnw-
Cxut     $	          fihqi(j)*tabnw,'     ! alfa'
Cxut	        end if
	        IF (SLASK.LT.DELMIN) THEN
	          DELMIN = SLASK
	          JDEL = J
	          CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	          errlok(ihr) = slask
	        END IF
	      END DO
	      IF (JPX.NE.0) THEN
	        IF (JCOLX.GE.3) THEN
	          DELMIN = SLASK
	          JDEL = JCOLX
	          CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	          errlok(ihr) = slask
	        END IF   ! last column "overrides"
	      END IF
	      CWORK(IA_TAB(IHR,JCOLX)) = USLASK ! val(jcolx,i) - val(jcolx-1,i)
	      CWORK(IA_YTAB(IHR,JCOLX)) = YUSLASK
C	                        ! yval(jcolx,i)-yval(jcolx-1,i)
C	      check if sufficient accuracy in this subinterval
	      IF ( DELMIN .GT. EPSLOK ) OKSUB = .FALSE.
	      CWORK(IA_TABNEW(IHR)) = DCMPLX(DELMIN,DBLE(JDEL))
	      IHR = IHR + LI
	    END DO
	  END DO

	ELSE        ! ixtyp=-4 in this case

C	  rational-function Bulirsch-Stoer extrapolation but "without x terms" 
C	  in this case;  (my,ny) order: (0,0) (0,1) (0,2) (1,2) (2,2) ...
	  DO IR=1,LI
	    EPSLOK = HSHARE * EPS(IR)
	    IHR = IR
	    DO IH=1,NHYD
C	      extrapolate to i:th row (for the ihr:th integral)
C	      we use now  "tab(j,i)" := "val(j,i)" - "val(j-1,i)"  
C	       (obvious abuse of notation,  "val(0,i)" = 0 )
C	      we use now  "ytab(j,i)" := "yval(j,i)" - "yval(j-1,i)"  
C	       (obvious abuse of notation,  "yval(0,i)" = 0 )
C	      val(j,i) :=  val(j-1,i)  +  del(j-1,i) * qqq(j,i)  where 
C	       del(j,i) = cslask(j,i)-tab(j,i-1) = val(j,i)-val(j,i-1)
C	       cslask(j,i) = val(j,i) - val(j-1,i-1)
C	      yval(j,i) :=  yval(j-1,i)  +  ydel(j-1,i) * qqq(j,i)  where 
C	       ydel(j,i) = ycslask(j,i)-ytab(j,i-1) = yval(j,i)-yval(j,i-1)
C	       ycslask(j,i) = yval(j,i) - yval(j-1,i-1)
C	      fival(j,i) :=  fival(j-1,i)  +  fiqqq(j,i)/yval(j,i)
	      TABNW = CWORK(IA_TABNEW(IHR))  ! val(1,i)
	      CSLASK = TABNW       ! val(1,i) ( - val(0,i-1) )
	      USLASK = TABNW       ! val(1,i) ( - val(0,i) )
	      YTABNW = FIHQI(1)  ! yval(1,i)
	      YCSLASK = YTABNW       ! yval(1,i) ( - yval(0,i-1) )
	      YUSLASK = YTABNW       ! yval(1,i) ( - yval(0,i) )
	      FITABNW = TABNW      ! fival(1,i)
Clog	      write(99,*) 1,fitabnw,' ! ',sngl(epslok)
Cxut	      if (ihr.eq.1) then
Cxut	        if (jcolx.eq.2) then
Cxut	          write(99,*)
Cxut	          write(99,*)
Cxut	          write(99,*)
Cxut	          write(99,*) -1,'  ! NEW, first from first x row'
Cxut	          write(99,*) 1/(fihqi(2)-fihqi(1)),
Cxut     $	          cwork(ia_tab(ihr,1)),' &'
Cxut	        end if
Cxut	        write(99,*)
Cxut	        write(99,*) jcolx,' nr of columns'
Cxut	        write(99,*) 1/fihqi(1),fitabnw,' &'
Cxut	      end if
	      FITSLASK = (0.0D0,0.0D0)  ! fival(1,i) - val(1,i)
	      VSLASK = CWORK(IA_TAB(IHR,1))
C	                        ! val(1,i-1) ( - val(0,i-1) )
	      CWORK(IA_TAB(IHR,1)) = USLASK
C	                        ! val(1,i) ( - val(0,i) )
	      YVSLASK = ( FIHQI(2) - FIHQI(1) )
C	                        ! yval(1,i-1) ( - yval(0,i-1) )
	      B1SLASK = VSLASK*HQI(2)    ! val(j-1,i-1)*hqi(j) here
	      BSLASK = B1SLASK - TABNW ! tabnw=val(1,i)
	      IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $        DENOMIN ) THEN
	        Q1LOK = (CSLASK-VSLASK)/BSLASK ! val(1,i) - val(1,i-1)
	        Q2LOK = (YCSLASK-YVSLASK)/BSLASK ! yval(1,i)-yval(1,i-1)
	        USLASK = TABNW*Q1LOK  ! val(2,i) - val(1,i)
	        YUSLASK = TABNW*Q2LOK  ! yval(2,i) - yval(1,i)
	        CSLASK = B1SLASK*Q1LOK ! val(2,i) - val(1,i-1)
	        YCSLASK = B1SLASK*Q2LOK ! yval(2,i) - yval(1,i-1)
	      ELSE
	        USLASK = VSLASK    ! cslask = val(2,i) - val(1,i-1)  ok
	        YUSLASK = YVSLASK  ! ycslask = yval(2,i)-yval(1,i-1)  ok
	      END IF
	      Q2LOK = YTABNW + YUSLASK
	      IF ( DMAX1(DABS(DREAL(Q2LOK)),DABS(DIMAG(Q2LOK))) .GE. 
     $        DENOMIN ) THEN
	        FIUSLASK = ( (FIHQI(2)-YTABNW)*USLASK - 
     $	        YUSLASK*FITSLASK ) / Q2LOK    ! fival(2,i) - fival(1,i)
	      ELSE
	        FIUSLASK = USLASK     ! somewhat arbitrary choice
	      END IF
	      TABNW = TABNW + USLASK   ! val(2,i)
	      YTABNW = Q2LOK            ! yval(2,i)
	      FITABNW = FITABNW + FIUSLASK   ! fival(2,i)
	      FITSLASK = FITSLASK + ( FIUSLASK - USLASK )
C	                        ! fival(2,i) - val(2,i)
	      DELMIN = DMAX1 ( DABS(DREAL(FIUSLASK)), 
     $	      DABS(DIMAG(FIUSLASK)) )   ! cdabs ( fival(2,i)-fival(1,i) )
Clog	      write(99,*) 2,fiuslask,sngl(delmin)
Cxut	      if (ihr.eq.1) then
Cxut	        write(99,*) 1/fihqi(1),fitabnw
Cxut	        write(99,*) '       ',tabnw,'     ! tabnw'
Cxut	        write(99,*) '       ',ytabnw-fihqi(2),
Cxut     $	        '     ! beta'
Cxut	        write(99,*) '       ',ytabnw*fitabnw-
Cxut     $	        fihqi(2)*tabnw,'     ! alfa'
Cxut	      end if
	      JDEL = 2
	      CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	      errlok(ihr) = delmin
	      IF (JCOLX.GT.2) THEN
	        J = 3
	        VSLASK = CWORK(IA_TAB(IHR,2))
C	                        ! val(j-1,i-1) - val(j-2,i-1)
	        CWORK(IA_TAB(IHR,2)) = USLASK
C	                        ! val(j-1,i) - val(j-2,i)
	        YVSLASK = CWORK(IA_YTAB(IHR,2))
C	                        ! yval(j-1,i-1) - yval(j-2,i-1)
	        CWORK(IA_YTAB(IHR,2)) = YUSLASK
C	                        ! yval(j-1,i) - val(j-2,i)
	        W1LOK = CSLASK-VSLASK   ! val(j-1,i) - val(j-1,i-1)
	        B1SLASK = (TABNW-W1LOK)*HQI(J) ! val(j-1,i-1)*hqi(j)
	        BSLASK = B1SLASK - TABNW  ! tabnw=val(j-1,i)
	        W2LOK = YCSLASK-YVSLASK  ! yval(j-1,i)-yval(j-1,i-1)
	        IF (JCOLX.GE.5) YVALOLD = YTABNW - W2LOK  ! yval(2,i-1)
	        IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $          DENOMIN ) THEN
	          Q1LOK = W1LOK/BSLASK
	          Q2LOK = W2LOK/BSLASK
	          USLASK = TABNW*Q1LOK  ! val(j,i) - val(j-1,i)
	          YUSLASK = TABNW*Q2LOK  ! yval(j,i) - yval(j-1,i)
	          CSLASK = B1SLASK*Q1LOK ! val(j,i) - val(j-1,i-1)
	          YCSLASK = B1SLASK*Q2LOK ! yval(j,i) - yval(j-1,i-1)
	        ELSE
	          USLASK = VSLASK    ! cslask = val(j,i) - val(j-1,i-1)  ok
	          YUSLASK = YVSLASK  ! ycslask = yval(j,i)-yval(j-1,i-1)  ok
	        END IF
	        Q2LOK = YTABNW + YUSLASK
	        IF ( DMAX1(DABS(DREAL(Q2LOK)),DABS(DIMAG(Q2LOK))) .GE. 
     $          DENOMIN ) THEN
	          FIUSLASK = ( (FIHQI(J)-YTABNW)*USLASK - 
     $	          YUSLASK*FITSLASK ) / Q2LOK    ! fival(j,i) - fival(j-1,i)
	        ELSE
	          FIUSLASK = USLASK     ! somewhat arbitrary choice
	        END IF
	        TABNW = TABNW + USLASK   ! val(j,i)
	        YTABNW = Q2LOK            ! yval(j,i)
	        FITABNW = FITABNW + FIUSLASK   ! fival(j,i)
	        FITSLASK = FITSLASK + ( FIUSLASK - USLASK )
C	                        ! fival(j,i) - val(j,i)
	        SLASK = DMAX1 ( DABS(DREAL(FIUSLASK)), 
     $	        DABS(DIMAG(FIUSLASK)) )  ! cdabs ( fival(j,i)-fival(j-1,i) )
Clog	        write(99,*) j,fiuslask,sngl(slask)
Cxut	        if (ihr.eq.1) then
Cxut	          write(99,*) 1/fihqi(1),fitabnw
Cxut	          write(99,*) '       ',tabnw,'     ! tabnw'
Cxut	          write(99,*) '       ',ytabnw-fihqi(j),
Cxut     $	          '     ! beta'
Cxut	          write(99,*) '       ',ytabnw*fitabnw-
Cxut     $	          fihqi(j)*tabnw,'     ! alfa'
Cxut	        end if
	        IF (SLASK.LT.DELMIN) THEN
	          DELMIN = SLASK
	          JDEL = J
	          CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	          errlok(ihr) = slask
	        ELSE IF (JPX.NE.0) THEN
	          IF (JCOLX.EQ.J) THEN
	            DELMIN = SLASK
	            JDEL = J
	            CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	            errlok(ihr) = slask
	          END IF   ! last column "overrides"
	        END IF
	        DO J = 4, JCOLX
	          VSLASK = CWORK(IA_TAB(IHR,J-1))
C	                        ! val(j-1,i-1) - val(j-2,i-1)
	          CWORK(IA_TAB(IHR,J-1)) = USLASK
C	                        ! val(j-1,i) - val(j-2,i)
	          YVSLASK = CWORK(IA_YTAB(IHR,J-1))
C	                        ! yval(j-1,i-1) - yval(j-2,i-1)
	          CWORK(IA_YTAB(IHR,J-1)) = YUSLASK
C	                        ! yval(j-1,i) - val(j-2,i)
	          IF (J.NE.5) THEN
	            B1SLASK = VSLASK*HQI(J)
	            BSLASK = B1SLASK - CSLASK ! cslask=val(j-1,i)-val(j-2,i-1)
	            IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $              DENOMIN ) THEN
	              Q1LOK = (CSLASK-VSLASK)/BSLASK ! val(j-1,i)-val(j-1,i-1)
	              Q2LOK =(YCSLASK-YVSLASK)/BSLASK!yval(j-1,i)-yval(j-1,i-1)
	              USLASK = CSLASK*Q1LOK  ! val(j,i) - val(j-1,i)
	              YUSLASK = CSLASK*Q2LOK  ! yval(j,i) - yval(j-1,i)
	              CSLASK = B1SLASK*Q1LOK ! val(j,i) - val(j-1,i-1)
	              YCSLASK = B1SLASK*Q2LOK ! yval(j,i) - yval(j-1,i-1)
	            ELSE
	              USLASK = VSLASK    ! cslask = val(j,i) - val(j-1,i-1)  ok
	              YUSLASK = YVSLASK  ! ycslask = yval(j,i)-yval(j-1,i-1) ok
	            END IF
	          ELSE
	            ALFAOLD = CWORK(IA_YTAB(IHR,1)) ! artificial use (ytab)
	            QLOK = QHQI(1)*VSLASKOLD - 
     $	            (QHQI(2)*YVALOLD-QHQI(3))*(TABNW-CSLASK)
	            W1LOK = CSLASK-VSLASK
	            B1SLASK = ( VSLASKOLD*(TABNW-W1LOK) + 
     $	            QLOK*ALFAOLD ) * HQI(J)
	            QLOK = VSLASKOLD*TABNW + QLOK*ALFAVAL
	            BSLASK = B1SLASK - QLOK
	            IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $              DENOMIN ) THEN
	              Q1LOK = W1LOK/BSLASK
	              Q2LOK = (YCSLASK-YVSLASK)/BSLASK
	              USLASK = QLOK*Q1LOK
	              YUSLASK = QLOK*Q2LOK
	              CSLASK = B1SLASK*Q1LOK
	              YCSLASK = B1SLASK*Q2LOK
	            ELSE
	              USLASK = VSLASK
	              YUSLASK = YVSLASK
	            END IF
	          END IF
	          Q2LOK = YTABNW + YUSLASK
	          IF ( DMAX1(DABS(DREAL(Q2LOK)),DABS(DIMAG(Q2LOK))) .GE. 
     $            DENOMIN ) THEN
	            FIUSLASK = ( (FIHQI(J)-YTABNW)*USLASK - 
     $	            YUSLASK*FITSLASK ) / Q2LOK    ! fival(j,i) - fival(j-1,i)
	          ELSE
	            FIUSLASK = USLASK     ! somewhat arbitrary choice
	          END IF
	          TABNW = TABNW + USLASK   ! val(j,i)
	          YTABNW = Q2LOK            ! yval(j,i)
	          FITABNW = FITABNW + FIUSLASK   ! fival(j,i)
	          FITSLASK = FITSLASK + ( FIUSLASK - USLASK )
C	                        ! fival(j,i) - val(j,i)
	          SLASK = DMAX1 ( DABS(DREAL(FIUSLASK)), 
     $	          DABS(DIMAG(FIUSLASK)) )  ! cdabs ( fival(j,i)-fival(j-1,i) )
Clog	          write(99,*) j,fiuslask,sngl(slask)
Cxut	          if (ihr.eq.1) then
Cxut	            write(99,*) 1/fihqi(1),fitabnw
Cxut	            write(99,*) '       ',tabnw,'     ! tabnw'
Cxut	            write(99,*) '       ',ytabnw-fihqi(j),
Cxut     $	            '     ! beta'
Cxut	            write(99,*) '       ',ytabnw*fitabnw-
Cxut     $	            fihqi(j)*tabnw,'     ! alfa'
Cxut	          end if
	          IF (SLASK.LT.DELMIN) THEN
	            DELMIN = SLASK
	            JDEL = J
	            CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	            errlok(ihr) = slask
	          END IF
	          IF (J.EQ.4) THEN
	            VSLASKOLD = VSLASK ! actually needed only for jcolx>=5
	            ALFAVAL = YTABNW*FITABNW - FIHQI(J)*TABNW
	          END IF
	        END DO
	        IF (JPX.NE.0) THEN
	          IF (JCOLX.GE.4) THEN
	            DELMIN = SLASK
	            JDEL = JCOLX
	            CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	            errlok(ihr) = slask
	          END IF
	        END IF   ! last column "overrides"
	        IF (4.LE.JCOLX) 
     $	        CWORK(IA_YTAB(IHR,1)) = ALFAVAL   ! artificial use (ytab)
	      END IF
	      CWORK(IA_TAB(IHR,JCOLX)) = USLASK ! val(jcolx,i) - val(jcolx-1,i)
	      CWORK(IA_YTAB(IHR,JCOLX)) = YUSLASK
C	                        ! yval(jcolx,i)-yval(jcolx-1,i)
C	      check if sufficient accuracy in this subinterval
	      IF ( DELMIN .GT. EPSLOK ) OKSUB = .FALSE.
	      CWORK(IA_TABNEW(IHR)) = DCMPLX(DELMIN,DBLE(JDEL))
	      IHR = IHR + LI
	    END DO
	  END DO

	END IF

	GOTO 999

 5	CONTINUE    ! ixtyp=5
	IF (J.NE.5) STOP 43516
C	rational-function Bulirsch-Stoer extrapolation but "with 
C	denominator x term as when ixtyp=2,-2 and appropriate 
C	numerator x term to cancel the x terms in the Taylor development " 
C	in this case;  (my,ny) order: (0,0) (0,1) (1,1) (1,2) (2,2) ...
	DO IR=1,LI
	  EPSLOK = HSHARE * EPS(IR)
	  IHR = IR
	  DO IH=1,NHYD
C	    extrapolate to i:th row (for the ihr:th integral)
C	    we use now  "tab(j,i)" := "val(j,i)" - "val(j-1,i)"  
C	     (obvious abuse of notation,  "val(0,i)" = 0 )
C	    we use now  "ytab(j,i)" := "yval(j,i)" - "yval(j-1,i)"  
C	     (obvious abuse of notation,  "yval(0,i)" = 0 )
C	    val(j,i) :=  val(j-1,i)  +  del(j-1,i) * qqq(j,i)  where 
C	     del(j,i) = cslask(j,i)-tab(j,i-1) = val(j,i)-val(j,i-1)
C	     cslask(j,i) = val(j,i) - val(j-1,i-1)
C	    yval(j,i) :=  yval(j-1,i)  +  ydel(j-1,i) * qqq(j,i)  where 
C	     ydel(j,i) = ycslask(j,i)-ytab(j,i-1) = yval(j,i)-yval(j,i-1)
C	     ycslask(j,i) = yval(j,i) - yval(j-1,i-1)
C	    fival(j,i) :=  fival(j-1,i)  +  fiqqq(j,i)/yval(j,i)
	    TABNW = CWORK(IA_TABNEW(IHR))  ! val(1,i)
	    CSLASK = TABNW       ! val(1,i) ( - val(0,i-1) )
	    USLASK = TABNW       ! val(1,i) ( - val(0,i) )
	    YTABNW = FIHQI(1)  ! yval(1,i)
	    YCSLASK = YTABNW       ! yval(1,i) ( - yval(0,i-1) )
	    YUSLASK = YTABNW       ! yval(1,i) ( - yval(0,i) )
	    FITABNW = TABNW      ! fival(1,i)
Clog	    write(99,*) 1,fitabnw,' ! ',sngl(epslok)
Cxut	    if (ihr.eq.1) then
Cxut	      if (jcolx.eq.2) then
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*)
Cxut	        write(99,*) -1,'  ! NEW, first from first x row'
Cxut	        write(99,*) 1/(fihqi(2)-fihqi(1)),
Cxut     $	        cwork(ia_tab(ihr,1)),' &'
Cxut	      end if
Cxut	      write(99,*)
Cxut	      write(99,*) jcolx,' nr of columns'
Cxut	      write(99,*) 1/fihqi(1),fitabnw,' &'
Cxut	    end if
	    FITSLASK = (0.0D0,0.0D0)  ! fival(1,i) - val(1,i)
	    VSLASK = CWORK(IA_TAB(IHR,1))
C	                      ! val(1,i-1) ( - val(0,i-1) )
	    CWORK(IA_TAB(IHR,1)) = USLASK
C	                      ! val(1,i) ( - val(0,i) )
	    YVSLASK = ( FIHQI(2) - FIHQI(1) )
C	                      ! yval(1,i-1) ( - yval(0,i-1) )
	    B1SLASK = VSLASK*HQI(2)    ! val(j-1,i-1)*hqi(j) here
	    BSLASK = B1SLASK - TABNW ! tabnw=val(1,i)
	    IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $      DENOMIN ) THEN
	      Q1LOK = (CSLASK-VSLASK)/BSLASK ! val(1,i) - val(1,i-1)
	      Q2LOK = (YCSLASK-YVSLASK)/BSLASK ! yval(1,i)-yval(1,i-1)
	      USLASK = TABNW*Q1LOK  ! val(2,i) - val(1,i)
	      YUSLASK = TABNW*Q2LOK  ! yval(2,i) - yval(1,i)
	      CSLASK = B1SLASK*Q1LOK ! val(2,i) - val(1,i-1)
	      YCSLASK = B1SLASK*Q2LOK ! yval(2,i) - yval(1,i-1)
	    ELSE
	      USLASK = VSLASK    ! cslask = val(2,i) - val(1,i-1)  ok
	      YUSLASK = YVSLASK  ! ycslask = yval(2,i)-yval(1,i-1)  ok
	    END IF
	    YTABNW = YTABNW + YUSLASK    ! yval(2,i)
	    Q2LOK = 2*YTABNW - FIHQI(2)
	    ZMOD = (YUSLASK-Q2LOK)*USLASK
	    IF ( DMAX1(DABS(DREAL(Q2LOK)),DABS(DIMAG(Q2LOK))) .GE. 
     $      DENOMIN ) THEN
	      Q1LOK = FITSLASK
	      FITSLASK = ZMOD/Q2LOK         ! fival(2,i) - val(2,i)
	      Q2LOK = FITSLASK - Q1LOK
	    END IF ! else "fitslask unchanged" (somewhat arbitrary choice)
	    TABNW = TABNW + USLASK   ! val(2,i)
	    FITABNW = TABNW + FITSLASK   ! fival(2,i)
	    FIUSLASK = Q2LOK + USLASK    ! fival(2,i) - fival(1,i)
	    DELMIN = DMAX1 ( DABS(DREAL(FIUSLASK)), 
     $	    DABS(DIMAG(FIUSLASK)) )   ! cdabs ( fival(2,i)-fival(1,i) )
Clog	    write(99,*) 2,fiuslask,sngl(delmin)
Cxut	    if (ihr.eq.1) then
Cxut	      write(99,*) 1/fihqi(1),fitabnw
Cxut	      write(99,*) '       ',tabnw,'     ! tabnw'
Cxut	      write(99,*) '       ',ytabnw-fihqi(2),
Cxut     $	      '     ! beta'
Cxut	      write(99,*) '       ',ytabnw*fitabnw-
Cxut     $	      fihqi(2)*tabnw,'     ! alfa'
Cxut	    end if
	    JDEL = 2
	    CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	    errlok(ihr) = delmin
	    DO J = 3, JCOLX
	      VSLASK = CWORK(IA_TAB(IHR,J-1))
C	                    ! val(j-1,i-1) - val(j-2,i-1)
	      CWORK(IA_TAB(IHR,J-1)) = USLASK
C	                    ! val(j-1,i) - val(j-2,i)
	      YVSLASK = CWORK(IA_YTAB(IHR,J-1))
C	                    ! yval(j-1,i-1) - yval(j-2,i-1)
	      CWORK(IA_YTAB(IHR,J-1)) = YUSLASK
C	                    ! yval(j-1,i) - val(j-2,i)
	      B1SLASK = VSLASK*HQI(J)
	      BSLASK = B1SLASK - CSLASK ! cslask=val(j-1,i)-val(j-2,i-1)
	      IF ( DMAX1(DABS(DREAL(BSLASK)),DABS(DIMAG(BSLASK))) .GE. 
     $        DENOMIN ) THEN
	        Q1LOK = (CSLASK-VSLASK)/BSLASK ! val(j-1,i)-val(j-1,i-1)
	        Q2LOK = (YCSLASK-YVSLASK)/BSLASK!yval(j-1,i)-yval(j-1,i-1)
	        USLASK = CSLASK*Q1LOK  ! val(j,i) - val(j-1,i)
	        YUSLASK = CSLASK*Q2LOK  ! yval(j,i) - yval(j-1,i)
	        CSLASK = B1SLASK*Q1LOK ! val(j,i) - val(j-1,i-1)
	        YCSLASK = B1SLASK*Q2LOK ! yval(j,i) - yval(j-1,i-1)
	      ELSE
	        USLASK = VSLASK    ! cslask = val(j,i) - val(j-1,i-1)  ok
	        YUSLASK = YVSLASK  ! ycslask = yval(j,i)-yval(j-1,i-1) ok
	      END IF
	      YTABNW = YTABNW + YUSLASK    ! yval(j,i)
	      Q2LOK = 2*YTABNW - FIHQI(J)
	      ZMOD = ZMOD + (YUSLASK-Q2LOK)*USLASK
	      IF ( DMAX1(DABS(DREAL(Q2LOK)),DABS(DIMAG(Q2LOK))) .GE. 
     $        DENOMIN ) THEN
	        Q1LOK = FITSLASK
	        FITSLASK = ZMOD/Q2LOK         ! fival(j,i) - val(j,i)
	        Q2LOK = FITSLASK - Q1LOK
	      END IF ! else "fitslask unchanged" (somewhat arbitrary choice)
	      TABNW = TABNW + USLASK   ! val(j,i)
	      FITABNW = TABNW + FITSLASK   ! fival(j,i)
	      FIUSLASK = Q2LOK + USLASK    ! fival(j,i) - fival(j-1,i)
	      SLASK = DMAX1 ( DABS(DREAL(FIUSLASK)), 
     $	      DABS(DIMAG(FIUSLASK)) )  ! cdabs ( fival(j,i)-fival(j-1,i) )
Clog	      write(99,*) j,fiuslask,sngl(slask)
Cxut	      if (ihr.eq.1) then
Cxut	        write(99,*) 1/fihqi(1),fitabnw
Cxut	        write(99,*) '       ',tabnw,'     ! tabnw'
Cxut	        write(99,*) '       ',ytabnw-fihqi(j),
Cxut     $	        '     ! beta'
Cxut	        write(99,*) '       ',ytabnw*fitabnw-
Cxut     $	        fihqi(j)*tabnw,'     ! alfa'
Cxut	      end if
	      IF (SLASK.LT.DELMIN) THEN
	        DELMIN = SLASK
	        JDEL = J
	        CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	        errlok(ihr) = slask
	      END IF
	    END DO
	    IF (JPX.NE.0) THEN
	      IF (JCOLX.GE.3) THEN
	        DELMIN = SLASK
	        JDEL = JCOLX
	        CWORK(IA_VALTERM(IHR)) = FITABNW
Cerr	        errlok(ihr) = slask
	      END IF   ! last column "overrides"
	    END IF
	    CWORK(IA_TAB(IHR,JCOLX)) = USLASK ! val(jcolx,i) - val(jcolx-1,i)
	    CWORK(IA_YTAB(IHR,JCOLX)) = YUSLASK
C	                      ! yval(jcolx,i)-yval(jcolx-1,i)
C	    check if sufficient accuracy in this subinterval
	    IF ( DELMIN .GT. EPSLOK ) OKSUB = .FALSE.
	    CWORK(IA_TABNEW(IHR)) = DCMPLX(DELMIN,DBLE(JDEL))
	    IHR = IHR + LI
	  END DO
	END DO
	GOTO 999

 999	CONTINUE

Cval	do ihr=1,nhyd*li         ! change nhyd*li if "not everything" desired
Cval	  j = 1
Cval	  cslask = cwork(ia_tab(ihr,j))
Cval	  write(50+j,*) ihr,cslask
Cval	  do j = 2, jcolx
Cval	    cslask = cslask + cwork(ia_tab(ihr,j))
Cval	    write(50+j,*) ihr,cslask
Cval	  end do
Cval	end do
Cval	oksub = oksubin

Cerr	if (oksub) then
Cerr	  do ihr=1,nhyd*li
Cerr	    valsum(ihr) = valsum(ihr) + cwork(ia_valterm(ihr))
Cerr	    errsum(ihr) = errsum(ihr) + errlok(ihr)
Cerr	  end do
Cerr	end if

	RETURN

Cerr	entry cerrout ( nhyd,li )
CerrC	<<< insert "call cerrout(nhyd,li)" afterwards in a "main" program >>>
Cerr	open(98,status='new')
Cerr	write(98,*)
Cerr	write(98,*) nhyd*li,'   ! nhydli'
Cerr	write(98,*)
Cerr	do ihr=1,nhyd*li
Cerr	  write(98,*) ihr,'   ! ihr'
Cerr	  write(98,*) sngl(errsum(ihr)),' ! errsum(ihr)'
Cerr	  write(98,*) sngl(cdabs(valsum(ihr))),' ! cdabs(valsum(ihr))'
Cerr	end do
Cerr	return

	END