C     sven ivansson 1990, starting from ilkka karasalos "cdabs.for"

C****************************************************************
C
C     ADINTXFI
C	(ADaptive INTegration with eXtrapolation and with FIlon)
C	integral of nhyd*li complex*16 functions (obtained by calling 
C	the subroutines dhcomp and rcomp) of a COMPLEX*16 argument by 
C	an adaptive extrapolation method (see DBA 7.2.2)
C	< Ilkka Karasalos "dcadbs.for" amended 1990 by Sven Ivansson with 
C	  *  option for rational-function Bulirsch-Stoer extrapolation
C	       despite its name the program dcadbs.for only contained 
C	       polynomial-function extrapolation
C	  *  options for Filon integration 
C	       note that "I(h)" will still be a series with the even powers 
C	       of h only, the extrapolation methods that are used (polynomial 
C	       and rational Bulirsch-Stoer) are "tuned" to this case
C	       difference simple/double Filon only for "end point" terms
C	       "trapezoidal non-Filon" :  ifilon = 0
C	       "trapezoidal Filon"     :  ifilon = 1,2
C	       "Simpson Filon"         :  ifilon = -1,-2
C	       "Simpson non-Filon"     :  ifilon = -3
C	       OBS!  for "Simpson" ( i.e.  ifilon < 0 ) , the option 
C	       ixtyp=2,-2 will not make explicit use of the fact that there 
C	       is a "gap" in the asymptotic expansion
C	     - simple (ifilon=+-1)
C	       " f(up) = g(up)*cdexp(gamma*up) "  where g is smooth 
C	       f    is one of the integrands obtained with calls to the 
C	            formal subroutines dhcomp and rcomp
C	            f(up) = dhfakt(up)*rfakt(up)
C	     - double (ifilon=+-2)
C	       " f(up) = gpos(up)*cdexp(gamma*up) + gneg(up)*cdexp(-gamma*up)
C	               =      fpos(up)           +      fneg(up)             "
C	       where gpos and gneg are smooth 
C	       f    is one of the integrands obtained with calls to the 
C	            formal subroutines dhcomp and rcomp
C	            f(up) = dhfakt(up)*rfakt(up)
C	       fpos is one of the integrands obtained with calls to the 
C	            formal subroutines dhcomp and rcomppos
C	            fpos(up) = dhfakt(up)*rfaktpos(up)
C	       fneg is thus obtained as  fneg(up) = f(up) - fpos(up) = 
C	            dhfakt(up) * ( rfakt(up) - rfaktpos(up) )
C	  *  integration over a polygon in the complex "up" plane rather 
C	       than over a line segment on the real "x" axis
C	  *  "matrix integrand" instead of "scalar integrand"
C	  *  it was possible to let the ii variable run from 1 to nst1-1 
C	       (nst2-1) rather than from 1 to 2*nst1-1 (2*nst2-1) and make 
C	       the storage more economical, furthermore the "val2" array is 
C	       not needed
C	  *  the storage requirements for  st2..  have been further 
C	       reduced by one third, no "double storage" is made
C	  *  introduction of an "indicator function" obtained in 
C	       "rlndet & cdet" ( function value  =  dexp(rlndet) * cdet ) 
C	       with calls to the formal subroutine dhcomp, this indicator 
C	       function is used for additional control of the integration 
C	       step with the understanding that the indicator function is 
C	       analytic and that the integrand factors "dhfakt" are in some 
C	       way inversely proportional to the indicator function (by 
C	       analytic factor functions)
C	       hence the integration path must not interfere with the zeroes 
C	       of the indicator function, small steps should be taken close 
C	       to such a zero and large steps (in general) far from these 
C	       zeroes
C	       logic for "early splitting" of subintervals based on "cdet 
C	       argument" differences has been introduced
C	  *  minor modifications, e.g. 
C	       - introduction of the "smax(pre)" logic
C	       - introduction of the "nback" logic
C	       - absolute tolerances ("eps(1:li)") rather than relative ones 
C	         based on current (dynamic) integral estimates               >
C	       - introduction of the "hkvar,lunkvar,nkvar" logic
C
C     FURTHER IDEAS
C         *  use the indicator function more delicately by also looking for 
C              "smooth parts" of each half of the interval and "splitting 
C              early" if such parts are found (Sven)
C         *  use the indicator function more delicately by adaptively updating 
C              control parameters (like the argdiflim:s) guided by the 
C              observed indicator-function variation over previously treated 
C              subintervals
C         *  include an "early split" facility by testing the integrals for 
C              each half of the interval (Ilkka)
C         *  include an "early split" facility by also using a more coarse 
C              "eps(1:li)" in the extrapolation table (Ilkka)
C
C     USEFUL TEST after program modifications, use lunout>=90 (see below)
C	  let dhfakt(ih)(up)=up**(ih-1), ih=1,2,3,..,nhyd
C	    for ifilon=0,-3 let nhyd>0 be arbitrary and 
C	                    let li=1, rfakt(1)(up)=1
C	    for ifilon=+-1  let nhyd=2,3 for ifilon pos.,neg. and 
C	                    let rfakt(ir)(up)=cdexp(gamma(ir)*up) ir=1,2,..,li
C	    for ifilon=+-2  let nhyd=2,3 for ifilon pos.,neg. and 
C	                    let rfakt(ir)(up)=cdexp(gamma(ir)*up)+
C	                      2*cdexp(-gamma(ir)*up), hence rfaktpos(ir)(up)=
C	                      cdexp(gamma(ir)*up), ir=1,2,..,li
C	  all computed results in the j:th extrapolation column for 
C	    " j=(nhyd+1)/2,..,min0(i,jrus)    if  ifilon = 0 "     and 
C	    " j=1,..,min0(i,jrus)             if  ifilon = 1,2 "   and 
C	    " j=1,..,min0(i-2,jrus)           if  ifilon = -1,-2 " and 
C	    " j=(nhyd-1)/2,..,min0(i-2,jrus)  if  ifilon = -3 "   will then be 
C	    exactly correct in each case (note that nhyd=2,3 for ifilon=
C	    1,-1,2,-2) if ixtyp = 1,-1,0 (OBS) but the program does not "know" 
C           that these values are correct until the relevant values in the 
C	    (j+1):st column have been computed (because the error estimates 
C	    are based on differences from the preceding column) (OBS)
C	  use ixtyp = 1,-1,0, a large negative smaxpre and small eps(1:li)
C	  test different hmax:s to obtain different stack depths ("nel"), 
C	    extrapolation depths ("i") and evaluation depths ("ieva")
C	  test the "back" logic by setting nelmaxpre=1,2
C	  test nrus=2,3,4,5,6
C	  try other ixtyp for comparison (the principles for the log output, 
C	    with the "why questions", are not changed)
C
C****************************************************************

	SUBROUTINE ADINTXFI ( NHYD,LI, DHCOMP,RCOMP, 
     $	IFILON,RCOMPPOS,GAMMA, NPU,UPHOERN, 
     $	HMAX,SMAXPRE,EPSH01Q,EPS,NELMAXPRE,NBACKMAX, 
     $	IXTYP,NRUS,IRUS,JRUS, ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, 
     $	HKVAR,LUNKVAR, NKVAR, 
     $	LUNOUT,TEXTLOG, 
     $	NHYDLI,VALUE, EPSMAGN, NFCALLS, 
     $	DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CWORK )  ! & IXGAP,JPX in /XCHANGE/

C	" call dhcomp (up, rlndet,cdet, dhfakt) "  will give rlndet & cdet 
C	     and dhfakt(1:nhyd)(up)
C	" call  rcomp (up,  rfakt) "  will give  rfakt(1:li)(up)
C	ifilon  is 0 for "trapezoidal non-Filon", 1 for "trapezoidal simple 
C	     Filon", 2 for "trapezoidal double Filon", -1 for "Simpson simple 
C	     Filon", -2 for "Simpson double Filon", -3 for "Simpson non-Filon"
C	" call rcomppos (up, rfakt) "  will give  rfaktpos(1:li)(up), which 
C	     is needed when ifilon=+-2
C	gamma(1:li)  are the exponents to be utilized in the Filon method
C	uphoern(1:npu)  are the corners of the integration-path polygon in 
C	     the complex "up" plane
C	hmax  is maximum allowed final "integration step (hcur)" along the 
C	     integration path
C	smaxpre  is a pre-set upper limit for "maximum allowed length of a 
C	     planned path segment (smax(lok))"   
C	     (  smaxpre >= 0.0d0  will give a computed value for smax ,  
C	     smaxpre < 0.0d0  will give  smax = -smaxpre  )
C	epsh01q  is for demanding less accuracy when things get hard, i.e. 
C	     when "h01" becomes small in relation to "stottot"; typically 
C	     epsh01q is set to 0.01d0 (0.0d0 will imply "not active"); the 
C	     output epsmagn is relevant to check the degree of "less demand"
C	eps(1:li)  are absolute error limits, eps(ir) is for the integrals 
C	     involving rfakt(ir)
C	nelmaxpre  is maximum "subinterval-stack depth (nel)", i.e.  number of 
C	     simultaneously active subintervals obtained by trials with 
C	     interval splitting for a planned path segment   
C	     ( nelmaxpre < 0  as well as  nelmaxpre = 0  will give a value, 
C	     for "nelmax", adapted to the current "icworklim0,icworkdim" )
C	nbackmax  is the maximum number of "back" steps with wasted 
C	     evaluations for each planned path segment
C	ixtyp  is 
C	     1,-1,0  for "polynomial-function extrapolation" 
C	     2,-2    for "rational-function Bulirsch-Stoer extrapolation", 
C	             "without using possible a priori knowledge about a gap 
C	             in the asymptotic series (when ifilon<0)" ,      see 
C	             the book by Gear (referred to in DBA), Bulirsch & Stoer 
C	             then recommend using nrus=11 (but nrus=8,7 seems better) 
C	             and jrus=7
C	     3,-3    for "inverse-polynomial-function extrapolation" 
C	     4 & -4  for "rational-function Bulirsch-Stoer extrapolation 
C	             using a priori knowledge about a gap in the asymptotic 
C	             series", only when ifilon<0
C	nrus  is maximum "extrapolation depth (i)", i.e. number of rows in 
C	     the extrapolation table  ( use nrus >= 2 , nrus odd is most 
C	     economical as concerns exploitation of the declared arrays )
C	     nrus =  2, 3 -->  ngroup =  1  , nst1 =    2, nst2 =    3
C	     nrus =  4, 5 -->  ngroup =  2  , nst1 =    4, nst2 =    6
C	     nrus =  6, 7 -->  ngroup =  3  , nst1 =    8, nst2 =   12
C	     nrus =  8, 9 -->  ngroup =  4  , nst1 =   16, nst2 =   24
C	     nrus = 10,11 -->  ngroup =  5  , nst1 =   32, nst2 =   48
C	     nrus = 12,13 -->  ngroup =  6  , nst1 =   64, nst2 =   96
C	     nrus = 14,15 -->  ngroup =  7  , nst1 =  128, nst2 =  192
C	     nrus = 16,17 -->  ngroup =  8  , nst1 =  256, nst2 =  384
C	     nrus = 18,19 -->  ngroup =  9  , nst1 =  512, nst2 =  768
C	     nrus = 20,21 -->  ngroup = 10  , nst1 = 1024, nst2 = 1536
C	     etc.
C	     note that   ngroup = nrus/2   and   nst1/2 = nst2/3
C	     note that   (nelmax-1)*(nst1/2) + nst1  <=  mxst1   
C	        i.e.   nelmax  <=  mxst1/(nst1/2) - 1   
C	        is a necessary condition
C	     set   nrus2 = 1 + 2*ngroup   so   
C	        nrus = nrus2  or  nrus = nrus2 - 1
C	     each subinterval is successively considered as 
C	              1  2 3  *  4 6   8 12   16 24   ...      nst1 nst2
C	        ( i = 1  2 3     4 5   6  7    8  9   ...   nrus2-1 nrus2   )
C	        ( j =     1       2     3       4      .         ngroup     )
C	        ( m =     1       2     4       8      .     2**(ngroup-1)  )
C	                    ( * pair-wise doubling from here )
C	        parts  ( "mhq(i)" )    ,    i=1,2,...,nrus
C	        for group j there are  2*m 3*m  parts where  m = 2**(j-1)
C	irus  is minimum allowed final "extrapolation depth (i)" for each 
C	     subinterval  ( use  2 <= irus <= nrus  if ixgaplok=0  ,  
C	     4 <= irus <= nrus  if ixgaplok=1 )
C	jrus  is maximum number of columns in the extrapolation table  
C	     ( use  2 <= jrus <= nrus  if ixgaplok=0  ,  
C	     2 <= jrus <= nrus-2  if ixgaplok=1 )
C	argdiflim1,argdiflim0,argdiflim2  are certain limits for the maximal 
C            argument differences (in degrees) between consecutive points 
C	     used at a particular extrapolation depth for a particular 
C	     subinterval  ( use  argdiflim1 < argdiflim0 < argdiflim2 ) 
C	     * argdiflim0 concerns the whole subinterval
C	     * argdiflim1,argdiflim2 can be used to bring about an early 
C	       subinterval splitting (which can be economical) (the logic is 
C	       adapted to cases of simple poles, i.e. zeroes of the indicator 
C	       function, close to the integration path)
C	hkvar  is strict "previous hcur" limit for postponing the treatment 
C	     of a certain subinterval (such a subinterval can be postponed 
C	     only when hkvar is positive and when the subinterval is on the 
C	     positive part of the real up axis)
C	lunkvar  is the unit number for the (formatted) file on which the 
C	     "postponed subintervals" are written (only used when hkvar is 
C	     positive)
C	nkvar  is the number of "postponed subintervals" (output) 
C	     note that nkvar:=0 unless hkvar is positive
C	lunout  is the unit number for log output  ( no log if  lunout <= 0 , 
C	     "evaluation density pre-plot" on "lunout+1" if  lunout = 6  or  
C	     lunout >= 80  , extended "self test appropriate for integrands 
C	     giving correct computed results for all i=1,2,..,nrus" if  
C	     lunout >= 90  )
C	textlog  is a character*50 header for the log  (only relevant when  
C	     lunout > 0 )
C	nhydli  is nhyd*li (input)
C	value(1:nhydli)  is input/output (OBS), 
C	     the integral of  dhfakt(ih)*rfakt(ir)  is added to  
C	     value(ihr)  where  ihr = (ih-1)*li + ir , 
C	     ih=1,2,...,nhyd  ir=1,2,...,li
C	epsmagn  sums the "hshare" values (output), if epsh01q is zero it 
C	     will be 1.0d0
C	nfcalls  counts the performed "function calls" (output)
C	work-space:
C	  dhfakt(1:nhyd)  rfakt(1:li)
C	         ( actually copies of "cwork"-entries, used for convenience )
C	  icworklim0,icworkdim
C	  cwork(1:icworkdim)   "containing cwork(icworklim0+1,icworkdim)" 
C	                       ( cwork(1:icworklim0) is not accessed )     
C	     rfaktpos ( 1:li , 0:nelmax )
C	     st1dh ( 1:nhyd , 0:((nelmax+1)*nst1h) )
C	     st1r ( 1:li , 0:((nelmax+1)*nst1h) )
C	     st2dh ( 1:nhyd , 1:((nelmax+1)*nst1h) )
C	     st2r ( 1:li , 1:((nelmax+1)*nst1h) )
C	     ss1 ( 1:nhydli )
C	     ss2 ( 1:nhydli )
C	     ssd ( 1:nhydli )
C	     ss ( 1:nhydli )
C	     tabnew ( 1:nhydli )
C	     valterm ( 1:nhydli )
C	     tab ( 1:nhydli , 1:(min0(nrus,jrus)) )
C	             ( note: first "tab column" irrelevant when ixgaplok=1 )
C	ixgap  (from /xchange/) indicates whether " i=1 & i=3 " are to 
C	     be used for the extrapolation; if so "ixgaplok=0 (no gap)", 
C	     else "ixgaplok=1 (gap)" ;   when  ifilon >= 0  ixgaplok:=ixgap
C	     this is a technical matter due to the fact that 
C	     "Simpson (Filon)" needs an even number of subintervals to 
C	     produce a value and that it was convenient to keep the 
C	     structure for computing the trapezoidal sums etc which are 
C	     still useful building blocks
C	     * with ixgaplok=0:  ( when ifilon>=0 this is the normal choice )
C	             iminx=2 
C	               (i.e. i = 1,  2,  3,  4,  5,  6,   7,...
C	                         H,H/2,H/3,H/4,H/6,H/8,H/12,...
C	                      j = 1,2,3,...,i-1,i )
C	     * with ixgaplok=1:  ( NOTE: when ifilon<0 ixgaplok:=1 !!! )
C	             iminx=4 
C	               ( i.e. i =     2,      4,  5,  6,   7,...
C	                            H/2,    H/4,H/6,H/8,H/12,...
C	                      j =   1,2,3,...,i-2  ; but j=1 for i=2 
C	                 so it can be viewed as  " H --> H/2 " )
C	jpx  (from /xchange/) indicates whether "all extrapolated columns" 
C	     will be "finally possible" in the extrapolation; if so "jpx=0" 
C	     else "jpx=1 (final column overrides)"

	IMPLICIT NONE

	INTEGER NTAB,JTAB, MXHQ, MXST1
	REAL*8 DLIMEXP
	PARAMETER ( NTAB = 17 , JTAB = 10 , MXHQ = 99 , MXST1 = 1024, 
     $  DLIMEXP = 75.0D0 )
C	it is often economical to have mxst1 as an integer power of two 
C	it is economical to have ntab odd

	INCLUDE 'IADR0.INC'

	INTEGER IXGAP,JPX
	COMMON /XCHANGE/ IXGAP,JPX  ! can be changed from 0 as given in "data"
	DATA IXGAP,JPX /0,0/
	INTEGER IXGAPLOK

	INTEGER NHYD,LI, IFILON, NPU, NELMAXPRE,NBACKMAX, 
     $	IXTYP,NRUS,IRUS,JRUS, LUNKVAR, LUNOUT, NHYDLI
	REAL*8 HMAX,SMAXPRE, EPSH01Q, EPS(*), 
     $	ARGDIFLIM1,ARGDIFLIM0,ARGDIFLIM2, HKVAR
	COMPLEX*16 GAMMA(*), UPHOERN(*)
	CHARACTER*50 TEXTLOG
	EXTERNAL DHCOMP,RCOMP, RCOMPPOS

	INTEGER NKVAR, NFCALLS
	REAL*8 EPSMAGN
	COMPLEX*16 VALUE(*)

	INTEGER ICWORKLIM0,ICWORKDIM
	COMPLEX*16 DHFAKT(*),RFAKT(*), CWORK(*)

	INTEGER*2 I2OUT,IEVA2OUT, N2OUT
	INTEGER NELMAX, NGROUP,NRUS2, IRUSP2, NST1,NST2, 
     $	NST1H,NST2H, NST1M1,NST2M1, IH,IR, IHR, JL, NEL,NEL1, NBACK, 
     $	NFCALLS0, NFCALLSTST, IB0,IB,IBQ, ILOOP, I, IEV, ISTACK, 
     $	ISTA1,ISTA2, ISTP1,ISTP2, II,IIG3, JJ,JJG3, ISTEP,IISTEP, 
     $	M21,M22, II21,II22, I21,I22, J,M, IMINX, JCOLX, LUNOP, NPLUNOP, 
     $	NWHY, MHQ3, JLKVAR, ISLASK, 
     $	IA_RFAKTPOS, IA_ST1DH, IA_ST1R, IA_ST2DH, IA_ST2R, IA_SS1, 
     $	IA_SS2, IA_SSD, IA_SS, IA_TABNEW, IA_VALTERM, IA_TAB, IA_YTAB, 
     $	IADR0,IADR1,IADR2,IADRQ, 
     $	IEVA(mxhq), ISTP1EVA(mxhq), ISTP2EVA(mxhq), MHQ(0:ntab)
	REAL*8 H01, HSHARE, HCUR, HMIN1,HMIN2, X0,X1,XM, 
     $	SMAX,SMAXLOK, STOT, STOTTOT, SLEFTLUNOP, SA,SB, ARGLEFT,ARGRIGHT, 
     $	ARGLRMAX, HSHAREKVAR, SLASK, SLASK1,SLASK2, RLNDET, 
     $	XST(0:mxhq), HQI(jtab,2:ntab), 
     $	HSUM(jtab), FIHQI(jtab,2:ntab), QHQI(3,2:ntab), 
     $	ARG1CDET(0:mxst1), ARG2CDET(mxst1)
	COMPLEX C2OUT0,C2OUT1
	COMPLEX*16 UPH0,UPH1, UPKLIV,UPRIKT,UPRIKTH,UPDELTA, BHYPR,
     $	CHYPR,SHYPR, UP, CSLASK, CSLASK1,CSLASK2, USLASK,
     $	VSLASK,BSLASK, UPHKVAR0,UPHKVAR1, CDET
	LOGICAL EVA, NOTMIXED2, OKSUB

	EXTERNAL IA_RFAKTPOS, IA_ST1DH, IA_ST1R, IA_ST2DH, IA_ST2R, 
     $	IA_SS1, IA_SS2, IA_SSD, IA_SS, IA_TABNEW, IA_VALTERM, IA_TAB, 
     $	IA_YTAB

	IF (IFILON.GE.0) THEN
	  IXGAPLOK = IXGAP
	  IF (IABS(IXTYP).GE.4) STOP 3501
	ELSE
	  IXGAPLOK = 1         ! "Simpson (Filon)" needs even nr of parts
	END IF

	IF (IXGAPLOK.EQ.0) THEN
	  IMINX = 2
	  IF (MIN0(NRUS,JRUS).GT.JTAB ) STOP 5165
	ELSE
	  IMINX = 4
	  IF (MIN0(NRUS-2,JRUS).GT.JTAB ) STOP 5165
	END IF

	NGROUP = NRUS/2
	NRUS2 = 1 + 2*NGROUP
	IRUSP2 = IRUS + 2
	IF ( NRUS.LT.IMINX .OR. NRUS2.GT.ntab ) STOP 5164
	IF ( IRUS.LT.IMINX .OR. JRUS.LT.2 ) STOP 5166
	MHQ(0) = 1
	MHQ(1) = 1
	MHQ(2) = 2
	MHQ(3) = 3
	M = 1
	I = 3
	DO J=2,NGROUP
	  M = 2*M
	  I = I+1
	  MHQ(I) = 2*M    ! even i,  2 * (2**(j-1))
	  I = I+1
	  MHQ(I) = 3*M    !  odd i,  3 * (2**(j-1))
	END DO
	NST1 = 2*M
	NST2 = 3*M
	NST1H = M
	NST2H = NST2/2  ! will only be used when nst2 is even
	NST1M1 = NST1 - 1
	NST2M1 = NST2 - 1

	IF (IXGAPLOK.NE.0) THEN
	  MHQ3 = MHQ(3)
	  MHQ(3) = MHQ(2)   ! temporary only, "i=2 present but not i=3"
	END IF
	IF (IFILON.GE.0) THEN
	  DO I = IMINX, NRUS
	    JCOLX = MIN0(I+2-IMINX,JRUS)
	    DO J = 2, JCOLX
	      HQI(J,I)  =  ( DBLE(MHQ(I)) / MHQ(I-J+1) ) ** 2
	    END DO  ! adapted to "I(h)" series with the even powers of h
	  END DO
	ELSE 
	  IF ( IABS(IXTYP).LE.1 .OR. IABS(IXTYP).EQ.3 ) THEN
	    DO I = IMINX, NRUS
	      JCOLX = MIN0(I+2-IMINX,JRUS)
	      SLASK = MHQ(I-1)**2
	      DO J = JCOLX, 3, -1
	        HSUM(J) = HSUM(J-1) + SLASK
	      END DO
	      HSUM(2) = SLASK
	      SLASK2 = MHQ(I)**2
	      DO J = 2, JCOLX
	        SLASK =  ( DBLE(MHQ(I)) / MHQ(I-J+1) ) ** 2
	        HQI(J,I) = SLASK + (SLASK-1)*SLASK2/HSUM(J)
	      END DO  ! "adapted to "I(h)" series with the 
C	                even powers of h except h**2
	    END DO
	  ELSE
	    DO I = IMINX, NRUS
	      JCOLX = MIN0(I+2-IMINX,JRUS)
	      DO J = 2, JCOLX
	        HQI(J,I) = ( DBLE(MHQ(I)) / MHQ(I-J+1) ) ** 2
	      END DO  ! adapted to "I(h)" series with the even powers of h
	    END DO
	    IF (IABS(IXTYP).EQ.4) THEN
	      DO I = IMINX, NRUS
	        JCOLX = MIN0(I+2-IMINX,JRUS)
	        DO J = 1, JCOLX
	          ISLASK = 0
	          DO M = I-J+1, I
	            ISLASK = ISLASK + MHQ(M)**2
	          END DO
	          FIHQI(J,I) = ISLASK
	        END DO
	      END DO
	      DO I = IMINX+3, NRUS
	        QHQI(1,I) = 1.0D0 / (MHQ(I-3)**2)
	        QHQI(2,I) = 1.0D0 / ((MHQ(I-2)*MHQ(I-1))**2)
	        QHQI(3,I) = 1.0D0/(MHQ(I-2)**2) + 
     $	        1.0D0/(MHQ(I-1)**2)
	      END DO
	    END IF
	  END IF
	END IF
	IF (IABS(IXTYP).NE.4) THEN
	  DO I = IMINX, NRUS
	    FIHQI(1,I) = MHQ(I)**2
	    FIHQI(2,I) = FIHQI(1,I) + MHQ(I-1)**2
	  END DO
	END IF  ! only for "cxut" in xpolbul
	IF (IXGAPLOK.NE.0) THEN
	  MHQ(3) = MHQ3   ! reset
	END IF

C	nelmax := "maximum depth of subinterval stack"
	NELMAX = ( ICWORKDIM-ICWORKLIM0 - (NHYD+LI)*(NST1+1) - LI - 
     $	(6+MIN0(NRUS,JRUS))*NHYDLI ) / ( LI + 2*(NHYD+LI)*NST1H )
	IF (NELMAXPRE.GT.0) THEN
	  IF ( NELMAXPRE.GT.MXHQ .OR. ((NELMAXPRE+1)*NST1H).GT.MXST1 
     $	  .OR. NELMAXPRE.GT.NELMAX ) STOP 441
	  NELMAX = NELMAXPRE
	ELSE
	  IF (NELMAXPRE.LT.0) NELMAX = MIN0(-NELMAXPRE,NELMAX)
	  NELMAX = MIN0 ( MXHQ , (MXST1/NST1H)-1 , NELMAX )
	  IF (NELMAX.LT.1) STOP 440
	END IF

	IF (SMAXPRE.GE.0.0D0) THEN
	  IF (NRUS.EQ.NRUS2) THEN
	    SMAX = 0.99D0 * (NST2*HMAX)    !  margin: factor 0.99d0
	  ELSE
	    SMAX = 0.99D0 * (NST1*HMAX)    !  margin: factor 0.99d0
	  END IF
	  IF (SMAXPRE.GT.0.0D0) THEN
	    IF (SMAX.GT.SMAXPRE) SMAX = SMAXPRE
	  END IF
	ELSE
	  SMAX = -SMAXPRE
	END IF
 
	STOTTOT = 0.0D0
	DO JL=2,NPU
	  STOTTOT = STOTTOT + CDABS(UPHOERN(JL)-UPHOERN(JL-1))
	END DO
	EPSMAGN = 0.0D0

	NHYDCOP = NHYD
	LICOP = LI
	NHYDLICOP = NHYDLI
	ISLASK = (NELMAX+1) * NST1H
	IADR0_RFAKTPOS = ICWORKLIM0
	IADR0_ST1DH = IADR0_RFAKTPOS + (1+NELMAX)*LI
	IADR0_ST1R = IADR0_ST1DH + (1+ISLASK)*NHYD
	IADR0_ST2DH = IADR0_ST1R + (1+ISLASK)*LI
	IADR0_ST2R = IADR0_ST2DH + ISLASK*NHYD
	IADR0_SS1 = IADR0_ST2R + ISLASK*LI
	IADR0_SS2 = IADR0_SS1 + NHYDLI
	IADR0_SSD = IADR0_SS2 + NHYDLI
	IADR0_SS = IADR0_SSD + NHYDLI
	IADR0_TABNEW = IADR0_SS + NHYDLI
	IADR0_VALTERM = IADR0_TABNEW + NHYDLI
	IADR0_TAB = IADR0_VALTERM + NHYDLI
	ISLASK = MIN0(NRUS,JRUS)*NHYDLI
	IF (IABS(IXTYP).NE.4) THEN
	  ISLASK = IADR0_TAB + ISLASK
	ELSE
	  IADR0_YTAB = IADR0_TAB + ISLASK
	  ISLASK = IADR0_YTAB + ISLASK
	END IF
	IF ( ISLASK .GT. ICWORKDIM ) THEN
	  WRITE(*,*)
	  WRITE(*,*) '%%% ADINTXFI STOP,  would need  icworkdim >= ',ISLASK
	  STOP 160
	END IF

	IF (LUNOUT.GT.0) THEN
	  IF (LUNOUT.NE.6) OPEN(LUNOUT,FILE='adintxfi.lunout')
	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' ADINTXFI LUNOUT-LOG;    ',TEXTLOG
	  WRITE(LUNOUT,*) ' ifilon,ixgaplok,npu ',IFILON,IXGAPLOK,NPU
	  WRITE(LUNOUT,*) ' hmax smax epsh01q',SNGL(HMAX),SNGL(SMAX),
     $	  SNGL(EPSH01Q)
	  WRITE(LUNOUT,*) ' nelmax,nbackmax ',NELMAX,NBACKMAX
	  WRITE(LUNOUT,*) ' ixtyp nrus irus jrus ',IXTYP,NRUS,IRUS,JRUS
	  WRITE(LUNOUT,*) '     nst1,nst2 ',NST1,NST2
	  WRITE(LUNOUT,*) '     needed icworkdim ',ISLASK
	  WRITE(LUNOUT,*) ' argdiflim1 argdiflim0 argdiflim2 ',
     $	  SNGL(ARGDIFLIM1),SNGL(ARGDIFLIM0),SNGL(ARGDIFLIM2)
	  WRITE(LUNOUT,*) ' hkvar,lunkvar ',SNGL(HKVAR),LUNKVAR
	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' X0   HCUR I ARGDIFMAX   IEVA   ( x1 = ',
     $	  'previous x0 or sb )'
	  WRITE(LUNOUT,*) ' ( i << ieva  indicates bad economy )'
	  WRITE(LUNOUT,*) ' -------------------------'
	  NFCALLS0 = 1
	  NFCALLSTST = 1
	  IF ( LUNOUT.EQ.6 .OR. LUNOUT.GE.80 ) THEN
	    LUNOP = LUNOUT + 1
	    OPEN(LUNOP,FILE='adintxfi.lunop')
	    WRITE(LUNOP,*) 1,'  ! lines to skip'
	    WRITE(LUNOP,*) ' ADINTXFI LUNOP-LOG;    ',TEXTLOG
	    WRITE(LUNOP,*) 1
	    WRITE(LUNOP,*) '     ! nplunop (to be inserted by hand)'
	    NPLUNOP = 0
	  ELSE
	    LUNOP = 0
	  END IF
	  IF (LUNOUT.GE.90) NWHY = 0
	END IF
	SLEFTLUNOP = STOTTOT

	NKVAR = 0

	NFCALLS = 1
	UP = UPHOERN(NPU)
	CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	ARG1CDET(0) = DATAN2D(DIMAG(CDET),DREAL(CDET))
	IADR0 = IA_ST1DH(0,0)
	DO IH=1,NHYD
	  CWORK(IADR0+IH) = DHFAKT(IH)
	END DO
	CALL RCOMP (UP, RFAKT)
	IADR0 = IA_ST1R(0,0)
	DO IR=1,LI
	  CWORK(IADR0+IR) = RFAKT(IR)
	END DO
	IF (IABS(IFILON).EQ.2) THEN
	  IADR1 = IA_RFAKTPOS(1,0)
	  CALL RCOMPPOS ( UP, CWORK(IADR1) )
	END IF

	DO 100 JL = NPU, 2, -1
C	  for the segments of the "integration-path polygon"

	  UPH0 = UPHOERN(JL-1)
	  UPH1 = UPHOERN(JL)
	  UPKLIV = UPH1 - UPH0
	  STOT = CDABS(UPKLIV)
	  IF (STOT.GT.0.0D0) THEN
	    UPRIKT = UPKLIV/STOT
	  ELSE
	    GOTO 100
	  END IF
	  UPRIKTH = 0.5D0 * UPRIKT
	  SMAXLOK = SMAX
	  IF (IFILON.GT.0) THEN
	    DO IR=1,LI
	      SLASK = DABS ( DREAL (GAMMA(IR)*UPRIKT) )
	      IF (SLASK.NE.0.0D0) SMAXLOK = DMIN1 ( SMAXLOK , DLIMEXP/SLASK )
	    END DO      ! to avoid overflow/underflow in hyprest
	  ELSE IF (IFILON.LT.0) THEN
	    IF (IFILON.GT.-3) THEN
	      DO IR=1,LI
	        SLASK = 2 * DABS ( DREAL (GAMMA(IR)*UPRIKT) )
	        IF (SLASK.NE.0.0D0) SMAXLOK = DMIN1 ( SMAXLOK , DLIMEXP/SLASK )
	      END DO      ! to avoid overflow/underflow in hyprest2
	    END IF
	  END IF
	  ISLASK = 1 + STOT/SMAXLOK  ! get equal-length "planned path segments"
	  SMAXLOK = 1.01D0*STOT/ISLASK  ! margin: factor 1.01d0
	  IF (LUNOUT.GT.0) THEN
	    WRITE(LUNOUT,*)
	    I2OUT = JL
	    C2OUT0 = UPH0
	    C2OUT1 = UPH1
	    WRITE(LUNOUT,*) ' jl,uph0&1 ',I2OUT,C2OUT0,C2OUT1
	  END IF
	  SLEFTLUNOP = SLEFTLUNOP - STOT

	  SA = STOT

	  DO 200 WHILE (SA.NE.0.0D0)
C	    for the "planned path segments"

	    NBACK = 0
	    SB = SA
	    SA = DMAX1(0.0D0,SB-SMAXLOK)
	    XST(0) = SA
	    XST(1) = SB
	    IF (LUNOUT.GT.0) WRITE(LUNOUT,*) '   sa,sb ',SNGL(SA),
     $	    SNGL(SB)        !  sa,sb  is a "planned path segment"

	    ARG1CDET(NST1) = ARG1CDET(0)
	    IADR0 = IA_ST1DH(0,0)
	    IADR1 = IA_ST1DH(0,NST1)
	    DO IH=1,NHYD
	      CWORK(IADR1+IH) = CWORK(IADR0+IH)
	    END DO
	    IADR0 = IA_ST1R(0,0)
	    IADR1 = IA_ST1R(0,NST1)
	    DO IR=1,LI
	      CWORK(IADR1+IR) = CWORK(IADR0+IR)
	    END DO

	    NFCALLS = NFCALLS + 1
	    UP = UPH0 + SA*UPRIKT
	    CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	    ARG1CDET(0) = DATAN2D(DIMAG(CDET),DREAL(CDET))
	    IADR0 = IA_ST1DH(0,0)
	    DO IH=1,NHYD
	      CWORK(IADR0+IH) = DHFAKT(IH)
	    END DO
	    CALL RCOMP (UP, RFAKT)
	    IADR0 = IA_ST1R(0,0)
	    DO IR=1,LI
	      CWORK(IADR0+IR) = RFAKT(IR)
	    END DO

	    IF (IABS(IFILON).EQ.2) THEN
	      IADR0 = IA_RFAKTPOS(0,0)
	      IADR1 = IA_RFAKTPOS(0,1)
	      DO IR=1,LI
	        CWORK(IADR1+IR) = CWORK(IADR0+IR)
	      END DO
	      IADR1 = IA_RFAKTPOS(1,0)
	      CALL RCOMPPOS ( UP, CWORK(IADR1) )   !  up = uph0 + sa*uprikt
	    END IF

C	    initialize "subinterval stack" indexed by nel,   attributes are 
C	      xst(0:nel)  rfaktpos(1:li,0:nel)
C	      ieva(1:nel)  istp1eva(1:nel)  istp2eva(1:nel)
C	      arg1cdet(0:nelnum)  st1dh(1:nhyd,0:nelnum)  st1r(1:li,0:nelnum)
C	      arg2cdet(1:nelnum)  st2dh(1:nhyd,1:nelnum)  st2r(1:li,1:nelnum)
C	                                  where  nelnum = ((nel-1)+2) * nst1h
	    NEL = 1
	    IEVA(1) = 1
	    ISTP1EVA(1) = NST1
	    ISTP2EVA(1) = 0
	    IB0 = 0

	    DO 300 WHILE (.TRUE.)
C	      for the "subintervals" in the stack

C	      initialize step-sizes etc
	      X0 = XST(NEL-1)
	      X1 = XST(NEL)
	      H01 = X1 - X0       !  x0,x1  is a "subinterval"
	      HSHARE = EPSH01Q + H01/STOTTOT
	      HMIN1 = H01/NST1
	      HMIN2 = H01/NST2
	      IF ( HMIN1.EQ.0.0D0 .OR. HMIN2.EQ.0.0D0 ) THEN
	        WRITE(*,*)
	        WRITE(*,*) ' %%% ADINTXFI "h01 almost zero" at '
	        WRITE(*,*) '     up = ',UPH0+X0*UPRIKT
	        WRITE(*,*) ' up(sa) = ',UPH0+SA*UPRIKT
	        WRITE(*,*) ' up(sb) = ',UPH0+SB*UPRIKT
	        STOP 770
	      END IF
	      ISTP1 = NST1
	      ISTA1 = NST1H
	      ISTP2 = NST1H   ! i.e. nst2/3
	      ISTA2 = NST1H   ! i.e. nst2/3
	      IB = IB0 + NST1
	      IF (IFILON.EQ.0) THEN
	        IHR = 0
	        DO IH=1,NHYD
	          CSLASK1 = UPRIKTH * CWORK(IA_ST1DH(IH,IB0))
	          CSLASK2 = UPRIKTH * CWORK(IA_ST1DH(IH,IB))
	          DO IR=1,LI
	            IHR = IHR + 1
	            CWORK(IA_SS1(IHR)) = CSLASK1 * CWORK(IA_ST1R(IR,IB0)) + 
     $	            CSLASK2 * CWORK(IA_ST1R(IR,IB))
	            CWORK(IA_SS2(IHR)) = CWORK(IA_SS1(IHR))
C	            trapezoidal sum, mhq(1)=1 subinterval part
	            CWORK(IA_TAB(IHR,1)) = H01 * CWORK(IA_SS1(IHR))  ! should 
C	            ixgaplok be nonzero, this last sentence is unnecessary
	          END DO
	        END DO
	      ELSE IF (IFILON.EQ.1) THEN
	        UPDELTA = H01*UPRIKT
	        DO IR = 1, LI
	          CSLASK = GAMMA(IR)*UPDELTA
	          CALL HYPREST (CSLASK, CHYPR,SHYPR)
	          CHYPR = H01 * CHYPR
	          SHYPR = H01 * SHYPR
	          CSLASK1 = UPRIKTH * CWORK(IA_ST1R(IR,IB0))
	          CSLASK2 = UPRIKTH * CWORK(IA_ST1R(IR,IB))
	          IHR = IR
	          DO IH=1,NHYD
	            USLASK = CSLASK1 * CWORK(IA_ST1DH(IH,IB0))
	            VSLASK = CSLASK2 * CWORK(IA_ST1DH(IH,IB))
	            CWORK(IA_SS1(IHR)) = USLASK + VSLASK
	            CWORK(IA_SS2(IHR)) = CWORK(IA_SS1(IHR))
	            CWORK(IA_SSD(IHR)) = USLASK - VSLASK
C	            "trapez simple Filon" sum, mhq(1)=1 subinterval part
	            CWORK(IA_TAB(IHR,1)) = CHYPR * CWORK(IA_SS1(IHR)) + 
     $	            SHYPR * CWORK(IA_SSD(IHR))                       ! should 
C	            ixgaplok be nonzero, this last sentence is unnecessary
	            IHR = IHR + LI
	          END DO
	        END DO
	      ELSE IF (IFILON.EQ.2) THEN
	        UPDELTA = H01*UPRIKT
	        NEL1 = NEL - 1
	        DO IR = 1, LI
	          CSLASK = GAMMA(IR)*UPDELTA
	          CALL HYPREST (CSLASK, CHYPR,SHYPR)
	          CHYPR = H01 * CHYPR
	          SHYPR = H01 * SHYPR
	          CSLASK1 = UPRIKTH * CWORK(IA_ST1R(IR,IB0))
	          CSLASK2 = UPRIKTH * CWORK(IA_ST1R(IR,IB))
	          BSLASK = UPRIKTH * ( CWORK(IA_ST1R(IR,IB0)) - 
     $	          2*CWORK(IA_RFAKTPOS(IR,NEL1)) )
	          CSLASK = UPRIKTH * ( CWORK(IA_ST1R(IR,IB)) - 
     $	          2*CWORK(IA_RFAKTPOS(IR,NEL)) )
	          IHR = IR
	          DO IH=1,NHYD
	            USLASK = CWORK(IA_ST1DH(IH,IB0))
	            VSLASK = CWORK(IA_ST1DH(IH,IB))
	            CWORK(IA_SS1(IHR)) = CSLASK1*USLASK + CSLASK2*VSLASK
	            CWORK(IA_SS2(IHR)) = CWORK(IA_SS1(IHR))
	            CWORK(IA_SSD(IHR)) = -BSLASK*USLASK + CSLASK*VSLASK
C	            "trapez double Filon" sum, mhq(1)=1 subinterval part
	            CWORK(IA_TAB(IHR,1)) = CHYPR * CWORK(IA_SS1(IHR)) + 
     $	            SHYPR * CWORK(IA_SSD(IHR))                       ! should 
C	            ixgaplok be nonzero, this last sentence is unnecessary
	            IHR = IHR + LI
	          END DO
	        END DO
	      ELSE IF (IFILON.EQ.-1) THEN
	        UPDELTA = H01*UPRIKT
	        DO IR = 1, LI
	          CSLASK1 = UPRIKTH * CWORK(IA_ST1R(IR,IB0))
	          CSLASK2 = UPRIKTH * CWORK(IA_ST1R(IR,IB))
	          IHR = IR
	          DO IH=1,NHYD
	            USLASK = CSLASK1 * CWORK(IA_ST1DH(IH,IB0))
	            VSLASK = CSLASK2 * CWORK(IA_ST1DH(IH,IB))
	            CWORK(IA_SS1(IHR)) = USLASK + VSLASK
	            CWORK(IA_SS2(IHR)) = CWORK(IA_SS1(IHR))
	            CWORK(IA_SSD(IHR)) = USLASK - VSLASK
	            IHR = IHR + LI
	          END DO
	        END DO
	      ELSE IF (IFILON.EQ.-2) THEN
	        UPDELTA = H01*UPRIKT
	        NEL1 = NEL - 1
	        DO IR = 1, LI
	          CSLASK1 = UPRIKTH * CWORK(IA_ST1R(IR,IB0))
	          CSLASK2 = UPRIKTH * CWORK(IA_ST1R(IR,IB))
	          BSLASK = UPRIKTH * ( CWORK(IA_ST1R(IR,IB0)) - 
     $	          2*CWORK(IA_RFAKTPOS(IR,NEL1)) )
	          CSLASK = UPRIKTH * ( CWORK(IA_ST1R(IR,IB)) - 
     $	          2*CWORK(IA_RFAKTPOS(IR,NEL)) )
	          IHR = IR
	          DO IH=1,NHYD
	            USLASK = CWORK(IA_ST1DH(IH,IB0))
	            VSLASK = CWORK(IA_ST1DH(IH,IB))
	            CWORK(IA_SS1(IHR)) = CSLASK1*USLASK + CSLASK2*VSLASK
	            CWORK(IA_SS2(IHR)) = CWORK(IA_SS1(IHR))
	            CWORK(IA_SSD(IHR)) = -BSLASK*USLASK + CSLASK*VSLASK
	            IHR = IHR + LI
	          END DO
	        END DO
	      ELSE
	        UPDELTA = H01*UPRIKT
	        IHR = 0
	        DO IH=1,NHYD
	          CSLASK1 = UPRIKTH * CWORK(IA_ST1DH(IH,IB0))
	          CSLASK2 = UPRIKTH * CWORK(IA_ST1DH(IH,IB))
	          DO IR=1,LI
	            IHR = IHR + 1
	            USLASK = CSLASK1 * CWORK(IA_ST1R(IR,IB0))
	            VSLASK = CSLASK2 * CWORK(IA_ST1R(IR,IB))
	            CWORK(IA_SS1(IHR)) = USLASK + VSLASK
	            CWORK(IA_SS2(IHR)) = CWORK(IA_SS1(IHR))
	            CWORK(IA_SSD(IHR)) = USLASK - VSLASK
	          END DO
	        END DO
	      END IF
	      ISTACK = 1
 
	      DO 400 ILOOP = 2, NRUS
C	        for the "extrapolation rows"
C	        (note that nrus>=2 so the "400 loop" is never ignored)
C	        "i" indicates "extrapolation depth", "ieva" indicates 
C	        available "evaluation depth" for the extrapolation (for 
C	        "ieva"=-1,-2 the minus sign indicates one more available 
C	        evaluation ; whereas for "ieva"=-4 the minus sign indicates 
C	        two, or possibly only one, less available evaluation)
 
	        I = ILOOP

	        IF (IXGAPLOK.NE.0) THEN
	          IF (ILOOP.EQ.3) THEN
	            I = 4  ! "i=3" not necessarily needed
	            ISTACK = 1
	          ELSE IF (ILOOP.EQ.4) THEN
	            I = 3  ! this case was previously left, istack=2 now
	          END IF
	        END IF

	        IADR0 = IA_SS(0)
	        DO IHR=1,NHYDLI
	          CWORK(IADR0+IHR) = (0.0D0,0.0D0)
	        END DO

	        IF (ISTACK.EQ.1) THEN

C	          even i now,   ista1 = nst1/mhq(i)  and  istp1 = 2*ista1
	          HCUR = ISTA1*HMIN1
	          EVA = I.GT.IABS(IEVA(NEL))
	          DO II = ISTA1, NST1M1, ISTP1
	            IF (EVA) THEN
	              NFCALLS = NFCALLS + 1
	              UP = UPH0 + (X0+II*HMIN1)*UPRIKT
	              IB = IB0 + II
	              CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	              ARG1CDET(IB) = DATAN2D(DIMAG(CDET),DREAL(CDET))
	              IADR0 = IA_ST1DH(0,IB)
	              DO IH=1,NHYD
	                CWORK(IADR0+IH) = DHFAKT(IH)
	              END DO
	              CALL RCOMP (UP, RFAKT)
	              IADR0 = IA_ST1R(0,IB)
	              DO IR=1,LI
	                CWORK(IADR0+IR) = RFAKT(IR)
	              END DO
	            ELSE
	              IB = IB0 + II
	              IADR0 = IA_ST1DH(0,IB)
	              DO IH=1,NHYD
	                DHFAKT(IH) = CWORK(IADR0+IH)
	              END DO
	              IADR0 = IA_ST1R(0,IB)
	              DO IR=1,LI
	                RFAKT(IR) = CWORK(IADR0+IR)
	              END DO
	            END IF
	            IADR1 = IA_SS(0)
	            DO IH=1,NHYD
	              CSLASK = UPRIKT * DHFAKT(IH)
	              DO IR=1,LI
	                IADR1 = IADR1 + 1
	                CWORK(IADR1) = CWORK(IADR1) + CSLASK * RFAKT(IR)
	              END DO
	            END DO
	          END DO
	          IF (IFILON.GE.0) THEN
	            IADR0 = IA_SS(0)
	            IADR1 = IA_SS1(0)
	            DO IHR=1,NHYDLI
	              CWORK(IADR1+IHR) = CWORK(IADR1+IHR) + CWORK(IADR0+IHR)
	            END DO
	            IF (IFILON.EQ.0) THEN
	              DO IHR=1,NHYDLI
C	                trapezoidal sum, mhq(i) subinterval parts
	                CWORK(IA_TABNEW(IHR)) = HCUR * CWORK(IA_SS1(IHR))
	              END DO
	            ELSE
	              UPDELTA = HCUR*UPRIKT
	              DO IR = 1, LI
	                CSLASK = GAMMA(IR)*UPDELTA
	                CALL HYPREST (CSLASK, CHYPR,SHYPR)
	                CHYPR = HCUR * CHYPR
	                SHYPR = HCUR * SHYPR
	                IHR = IR
	                DO IH=1,NHYD
C	                  "trapez Filon" sum, mhq(i) subinterval parts
	                  CWORK(IA_TABNEW(IHR)) = 
     $	                  CHYPR * CWORK(IA_SS1(IHR)) + 
     $	                  SHYPR * CWORK(IA_SSD(IHR))
	                  IHR = IHR + LI
	                END DO
	              END DO
	            END IF
	          ELSE
	            IADR0 = IA_SS(0)
	            IADR1 = IA_SS1(0)
	            DO IHR=1,NHYDLI
	              CSLASK = CWORK(IADR0+IHR)
	              CWORK(IADR0+IHR) = CWORK(IADR1+IHR)
	              CWORK(IADR1+IHR) = CWORK(IADR1+IHR) + CSLASK
	            END DO
	            IF (IFILON.GT.-3) THEN
	              UPDELTA = HCUR*UPRIKT
	              DO IR = 1, LI
	                CSLASK = GAMMA(IR)*UPDELTA
	                CALL HYPREST2 (CSLASK, BHYPR,CHYPR,SHYPR)
	                BHYPR = HCUR * BHYPR
	                CHYPR = HCUR * CHYPR
	                SHYPR = HCUR * SHYPR
	                IHR = IR
	                DO IH=1,NHYD
C	                  "Simpson Filon" sum, mhq(i) subinterval parts
	                  CWORK(IA_TABNEW(IHR)) = 
     $	                  BHYPR * CWORK(IA_SS1(IHR)) + 
     $	                  CHYPR * CWORK(IA_SS(IHR)) + 
     $	                  SHYPR * CWORK(IA_SSD(IHR))
	                  IHR = IHR + LI
	                END DO
	              END DO
	            ELSE
	              BHYPR = HCUR * (4.0D0/3)
	              CHYPR = HCUR * (-2.0D0/3)
	              DO IHR=1,NHYDLI
C	                Simpson sum, mhq(i) subinterval parts
	                CWORK(IA_TABNEW(IHR)) = 
     $	                BHYPR * CWORK(IA_SS1(IHR)) + 
     $	                CHYPR * CWORK(IA_SS(IHR))
	              END DO
	            END IF
	          END IF
	          IF (IXGAPLOK.NE.0) THEN
	            IF (I.EQ.2) THEN
	              IADR0 = IA_TABNEW(0)
	              IADR2 = IA_TAB(0,1)
	              DO IHR=1,NHYDLI
	                CWORK(IADR2+IHR) = CWORK(IADR0+IHR)
	              END DO  ! "tab" but not "tabnew" needed in this case
	            END IF
	          END IF
C	          had here  mhq(i) = 2 * 2**((i-2)/2)  subinterval parts
	          IF (EVA) THEN
	            IF (I.GT.4) THEN
	              IEVA(NEL) = I
	            ELSE IF (I.EQ.4) THEN
	              IF (IEVA(NEL).EQ.3) THEN
	                IEVA(NEL) = I
	              ELSE
	                IEVA(NEL) = -4 ! ieva(nel) was 2,-2 here and ixgaplok=1
	              END IF
	            ELSE IF (IEVA(NEL).EQ.1) THEN
	              IEVA(NEL) = I  ! i=2 here
	            ELSE
	              IEVA(NEL) = -2 ! i=2 and ieva(nel) was -1 here
	            END IF
	            ISTP1EVA(NEL) = ISTA1
	          END IF
	          ISTP1 = ISTA1
	          ISTA1 = ISTA1/2
	          ISTACK = 2

	        ELSE

C	          odd i now,   ista2 = nst2/mhq(i)  and  istp2 = 2*ista2
	          HCUR = ISTA2*HMIN2
	          NOTMIXED2 = ISTP2EVA(NEL).GE.0
	          IF (NOTMIXED2) THEN
C	            ieva(nel) >= 2  or ieva(nel) = -4 in this case
	            EVA = I.GT.IEVA(NEL)
	            DO II = ISTA2, NST2M1, ISTP2
	              IIG3 = II/3
	              IF (II.EQ.(3*IIG3)) THEN
	                IB = IB0 + 2*IIG3
	                IADR0 = IA_ST1DH(0,IB)
	                DO IH=1,NHYD
	                  DHFAKT(IH) = CWORK(IADR0+IH)
	                END DO
	                IADR0 = IA_ST1R(0,IB)
	                DO IR=1,LI
	                  RFAKT(IR) = CWORK(IADR0+IR)
	                END DO
	              ELSE IF (EVA) THEN
	                NFCALLS = NFCALLS + 1
	                UP = UPH0 + (X0+II*HMIN2)*UPRIKT
	                IB = IB0 + II - IIG3
	                CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	                ARG2CDET(IB) = DATAN2D(DIMAG(CDET),DREAL(CDET))
	                IADR0 = IA_ST2DH(0,IB)
	                DO IH=1,NHYD
	                  CWORK(IADR0+IH) = DHFAKT(IH)
	                END DO
	                CALL RCOMP (UP, RFAKT)
	                IADR0 = IA_ST2R(0,IB)
	                DO IR=1,LI
	                  CWORK(IADR0+IR) = RFAKT(IR)
	                END DO
	              ELSE
	                IB = IB0 + II - IIG3
	                IADR0 = IA_ST2DH(0,IB)
	                DO IH=1,NHYD
	                  DHFAKT(IH) = CWORK(IADR0+IH)
	                END DO
	                IADR0 = IA_ST2R(0,IB)
	                DO IR=1,LI
	                  RFAKT(IR) = CWORK(IADR0+IR)
	                END DO
	              END IF
	              IADR1 = IA_SS(0)
	              DO IH=1,NHYD
	                CSLASK = UPRIKT * DHFAKT(IH)
	                DO IR=1,LI
	                  IADR1 = IADR1 + 1
	                  CWORK(IADR1) = CWORK(IADR1) + CSLASK * RFAKT(IR)
	                END DO
	              END DO
	            END DO
	          ELSE
C	            i = 3  so  ista2 = istp2 = nst2/3  in this case
	            IF (ISTP2EVA(NEL).EQ.-1) THEN
	              EVA = .TRUE.
	            ELSE
	              EVA = .FALSE.    !  istp2eva(nel) = -2  in this case
	            END IF
	            DO II = ISTA2, NST2M1, ISTP2
	              IIG3 = II/3
	              IF (EVA) THEN
	                NFCALLS = NFCALLS + 1
	                UP = UPH0 + (X0+II*HMIN2)*UPRIKT
	                IB = IB0 + II - IIG3
	                CALL DHCOMP (UP, RLNDET,CDET, DHFAKT)
	                ARG2CDET(IB) = DATAN2D(DIMAG(CDET),DREAL(CDET))
	                IADR0 = IA_ST2DH(0,IB)
	                DO IH=1,NHYD
	                  CWORK(IADR0+IH) = DHFAKT(IH)
	                END DO
	                CALL RCOMP (UP, RFAKT)
	                IADR0 = IA_ST2R(0,IB)
	                DO IR=1,LI
	                  CWORK(IADR0+IR) = RFAKT(IR)
	                END DO
	              ELSE
	                IB = IB0 + II - IIG3
	                IADR0 = IA_ST2DH(0,IB)
	                DO IH=1,NHYD
	                  DHFAKT(IH) = CWORK(IADR0+IH)
	                END DO
	                IADR0 = IA_ST2R(0,IB)
	                DO IR=1,LI
	                  RFAKT(IR) = CWORK(IADR0+IR)
	                END DO
	              END IF
	              IADR1 = IA_SS(0)
	              DO IH=1,NHYD
	                CSLASK = UPRIKT * DHFAKT(IH)
	                DO IR=1,LI
	                  IADR1 = IADR1 + 1
	                  CWORK(IADR1) = CWORK(IADR1) + CSLASK * RFAKT(IR)
	                END DO
	              END DO
	              EVA = .NOT.EVA
	            END DO
                    EVA = .TRUE.   ! another evaluation has been performed
	          END IF
	          IF (IFILON.GE.0) THEN
	            IADR0 = IA_SS(0)
	            IADR2 = IA_SS2(0)
	            DO IHR=1,NHYDLI
	              CWORK(IADR2+IHR) = CWORK(IADR2+IHR) + CWORK(IADR0+IHR)
	            END DO
	            IF (IFILON.EQ.0) THEN
	              DO IHR=1,NHYDLI
C	                trapezoidal sum, mhq(i) subinterval parts
	                CWORK(IA_TABNEW(IHR)) = HCUR * CWORK(IA_SS2(IHR))  ! 
C	                unnecessary when i=3 if ixgaplok=1
	              END DO
	            ELSE
	              UPDELTA = HCUR*UPRIKT
	              DO IR = 1, LI
	                CSLASK = GAMMA(IR)*UPDELTA
	                CALL HYPREST (CSLASK, CHYPR,SHYPR)
	                CHYPR = HCUR * CHYPR
	                SHYPR = HCUR * SHYPR
	                IHR = IR
	                DO IH=1,NHYD
C	                  "trapez Filon" sum, mhq(i) subinterval parts
	                  CWORK(IA_TABNEW(IHR)) = 
     $	                  CHYPR * CWORK(IA_SS2(IHR)) + 
     $	                  SHYPR * CWORK(IA_SSD(IHR))                       ! 
C	                  unnecessary when i=3 if ixgaplok=1
	                  IHR = IHR + LI
	                END DO
	              END DO
	            END IF
	          ELSE IF (I.LT.5) THEN
	            IADR0 = IA_SS(0)
	            IADR2 = IA_SS2(0)
	            DO IHR=1,NHYDLI
	              CWORK(IADR2+IHR) = CWORK(IADR2+IHR) + CWORK(IADR0+IHR)
	            END DO
	          ELSE
	            IADR0 = IA_SS(0)
	            IADR2 = IA_SS2(0)
	            DO IHR=1,NHYDLI
	              CSLASK = CWORK(IADR0+IHR)
	              CWORK(IADR0+IHR) = CWORK(IADR2+IHR)
	              CWORK(IADR2+IHR) = CWORK(IADR2+IHR) + CSLASK
	            END DO
	            IF (IFILON.GT.-3) THEN
	              UPDELTA = HCUR*UPRIKT
	              DO IR = 1, LI
	                CSLASK = GAMMA(IR)*UPDELTA
	                CALL HYPREST2 (CSLASK, BHYPR,CHYPR,SHYPR)
	                BHYPR = HCUR * BHYPR
	                CHYPR = HCUR * CHYPR
	                SHYPR = HCUR * SHYPR
	                IHR = IR
	                DO IH=1,NHYD
C	                  "Simpson Filon" sum, mhq(i) subinterval parts
	                  CWORK(IA_TABNEW(IHR)) = 
     $	                  BHYPR * CWORK(IA_SS2(IHR)) + 
     $	                  CHYPR * CWORK(IA_SS(IHR)) + 
     $	                  SHYPR * CWORK(IA_SSD(IHR))
	                  IHR = IHR + LI
	                END DO
	              END DO
	            ELSE
	              BHYPR = HCUR * (4.0D0/3)
	              CHYPR = HCUR * (-2.0D0/3)
	              DO IHR=1,NHYDLI
C	                Simpson sum, mhq(i) subinterval parts
	                CWORK(IA_TABNEW(IHR)) = 
     $	                BHYPR * CWORK(IA_SS2(IHR)) + 
     $	                CHYPR * CWORK(IA_SS(IHR))
	              END DO
	            END IF
	          END IF
C	          had here  mhq(i) = 3 * 2**((i-3)/2)  subinterval parts
	          IF (IEVA(NEL).LT.I) THEN
	            IEVA(NEL) = I   ! temporarily wrong if ieva(nel) was -4
	            ISTP2EVA(NEL) = ISTA2
	          END IF
C	          have now  ieva(nel) >= 3  and  nst2/3 >= istp2eva(nel) >= 1
	          ISTP2 = ISTA2
	          ISTA2 = ISTA2/2
	          IF (IXGAPLOK.NE.0) THEN
	            IF (I.EQ.3) GOTO 400   ! continue with i=5 & istack=2
	          END IF
	          ISTACK = 1

	        END IF

C	        can now leave the "400 loop" prematurely with "goto 390"
C	        (note that  i >= 2  here )
C	        (this is advantageous if parts of the present subinterval 
C	        exhibit a different behaviour for some integrand, continued 
C	        treatment of the whole subinterval would then be wasteful as 
C	        concerns the number of function evaluations)

	        IF ( (HCUR.LE.HMAX) .AND. (I.GE.IRUS) ) THEN
C	          computation of maximal "cdet argument" differences, in 
C	          degrees and for consecutive point pairs available from "i" 
C	          (rather than "ieva"), for the left and right halves of the 
C	          present subinterval
	          IF ( I .EQ. 2*(I/2)+1 ) THEN
	            ARGLEFT = 0.0D0
	            SLASK1 = ARG1CDET ( IB0 )
	            JJ = 0
	            DO II = ISTP1, NST1H, ISTP1
	              JJ = JJ + ISTP2
	              JJG3 = JJ/3
	              IF (JJ.EQ.(3*JJG3)) THEN
	                JJ = JJ + ISTP2
	                JJG3 = JJ/3
	              END IF
	              SLASK2 = ARG2CDET ( IB0 + JJ - JJG3 )
	              SLASK = DABS(SLASK2-SLASK1)
	              IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	              ARGLEFT = DMAX1(ARGLEFT,SLASK)
	              SLASK1 = SLASK2
	              SLASK2 = ARG1CDET ( IB0 + II )
	              SLASK = DABS(SLASK2-SLASK1)
	              IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	              ARGLEFT = DMAX1(ARGLEFT,SLASK)
	              SLASK1 = SLASK2
	            END DO
	            ARGRIGHT = 0.0D0
	            SLASK1 = ARG1CDET ( IB0 + NST1 )
	            JJ = NST2
	            DO II = NST1-ISTP1, NST1H, -ISTP1
	              JJ = JJ - ISTP2
	              JJG3 = JJ/3
	              IF (JJ.EQ.(3*JJG3)) THEN
	                JJ = JJ - ISTP2
	                JJG3 = JJ/3
	              END IF
	              SLASK2 = ARG2CDET ( IB0 + JJ - (JJ/3) )
	              SLASK = DABS(SLASK2-SLASK1)
	              IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	              ARGRIGHT = DMAX1(ARGRIGHT,SLASK)
	              SLASK1 = SLASK2
	              SLASK2 = ARG1CDET ( IB0 + II )
	              SLASK = DABS(SLASK2-SLASK1)
	              IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	              ARGRIGHT = DMAX1(ARGRIGHT,SLASK)
	              SLASK1 = SLASK2
	            END DO
	          ELSE
	            IF ( IXGAPLOK.EQ.0 .OR. I.NE.4 ) THEN
	              ARGLEFT = 0.0D0
	              SLASK1 = ARG1CDET ( IB0 )
	              JJ = 0
	              ISLASK = 1
	              DO II = ISTP1, NST1H, ISTP1
	                IF (ISLASK.LE.1) THEN
	                  SLASK2 = ARG1CDET ( IB0 + II )
	                  SLASK = DABS(SLASK2-SLASK1)
	                  IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                  ARGLEFT = DMAX1(ARGLEFT,SLASK)
	                  SLASK1 = SLASK2
	                  ISLASK = ISLASK + 1
	                ELSE
	                  JJ = JJ + ISTP2
	                  JJG3 = JJ/3
	                  IF (JJ.EQ.(3*JJG3)) THEN
	                    JJ = JJ + ISTP2
	                    JJG3 = JJ/3
	                  END IF
	                  SLASK2 = ARG2CDET ( IB0 + JJ - JJG3 )
	                  SLASK = DABS(SLASK2-SLASK1)
	                  IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                  ARGLEFT = DMAX1(ARGLEFT,SLASK)
	                  SLASK1 = SLASK2
	                  SLASK2 = ARG1CDET ( IB0 + II )
	                  SLASK = DABS(SLASK2-SLASK1)
	                  IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                  ARGLEFT = DMAX1(ARGLEFT,SLASK)
	                  SLASK1 = SLASK2
	                  ISLASK = ISLASK + 1
	                  IF (ISLASK.EQ.4) ISLASK = 0
	                END IF
	              END DO
	              ARGRIGHT = 0.0D0
	              SLASK1 = ARG1CDET ( IB0 + NST1 )
	              JJ = NST2
	              ISLASK = 1
	              DO II = NST1-ISTP1, NST1H, -ISTP1
	                IF (ISLASK.LE.1) THEN
	                  SLASK2 = ARG1CDET ( IB0 + II )
	                  SLASK = DABS(SLASK2-SLASK1)
	                  IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                  ARGRIGHT = DMAX1(ARGRIGHT,SLASK)
	                  SLASK1 = SLASK2
	                  ISLASK = ISLASK + 1
	                ELSE
	                  JJ = JJ - ISTP2
	                  JJG3 = JJ/3
	                  IF (JJ.EQ.(3*JJG3)) THEN
	                    JJ = JJ - ISTP2
	                    JJG3 = JJ/3
	                  END IF
	                  SLASK2 = ARG2CDET ( IB0 + JJ - JJG3 )
	                  SLASK = DABS(SLASK2-SLASK1)
	                  IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                  ARGRIGHT = DMAX1(ARGRIGHT,SLASK)
	                  SLASK1 = SLASK2
	                  SLASK2 = ARG1CDET ( IB0 + II )
	                  SLASK = DABS(SLASK2-SLASK1)
	                  IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                  ARGRIGHT = DMAX1(ARGRIGHT,SLASK)
	                  SLASK1 = SLASK2
	                  ISLASK = ISLASK + 1
	                  IF (ISLASK.EQ.4) ISLASK = 0
	                END IF
	              END DO
	            ELSE
	              ARGLEFT = 0.0D0
	              SLASK1 = ARG1CDET ( IB0 )
	              DO II = ISTP1, NST1H, ISTP1
	                SLASK2 = ARG1CDET ( IB0 + II )
	                SLASK = DABS(SLASK2-SLASK1)
	                IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                ARGLEFT = DMAX1(ARGLEFT,SLASK)
	                SLASK1 = SLASK2
	              END DO
	              ARGRIGHT = 0.0D0
	              SLASK1 = ARG1CDET ( IB0 + NST1 )
	              DO II = NST1-ISTP1, NST1H, -ISTP1
	                SLASK2 = ARG1CDET ( IB0 + II )
	                SLASK = DABS(SLASK2-SLASK1)
	                IF (SLASK.GT.180.0D0) SLASK = 360.0D0 - SLASK
	                ARGRIGHT = DMAX1(ARGRIGHT,SLASK)
	                SLASK1 = SLASK2
	              END DO
	            END IF
	          END IF
	          ARGLRMAX = DMAX1(ARGLEFT,ARGRIGHT)
	          OKSUB = ARGLRMAX.LE.ARGDIFLIM0
	        ELSE
	          ARGLRMAX = -1.0D0
	          OKSUB = .FALSE.
	        END IF

	        IF (I.GE.IMINX) THEN

	          JCOLX = MIN0(I+2-IMINX,JRUS)

Clog	          write(99,*) ' ------- x0,x1,i ',sngl(x0),sngl(x1),i
Cval	          do j = 1, jcolx
Cval	            write(50+j,*) i,mhq(i),' -------'
Cval	          end do
CvalC	          if so, rename to "adintxfiVAL" and 
CvalC	                 call "xpolbulVAL" instead of "xpolbul" below and 
CvalC	                 increase the "ntab,jtab,mxst1" parameters as needed
	          CALL XPOLBUL ( NHYD,LI, IXTYP, HSHARE,EPS, 
     $	          HQI(1,I),FIHQI(1,I),QHQI(1,I), JPX,JCOLX, 
     $	          IA_TABNEW,IA_VALTERM,IA_TAB,IA_YTAB, CWORK,OKSUB )

	        END IF

	        IF (OKSUB) THEN

	          IF (LUNOUT.GT.0) THEN
	            I2OUT = I            ! i2out = 2,3,4,5,..,nrus here
	            IEVA2OUT = IEVA(NEL) !ieva2out=-4, -2, 2,3,4,5,..,nrus here
	            WRITE(LUNOUT,*) '     ', SNGL(X0), SNGL(HCUR), 
     $	            I2OUT, SNGL(DMAX1(ARGLEFT,ARGRIGHT)), IEVA2OUT
	            IF ( IEVA2OUT .EQ. 2*(IEVA2OUT/2)+1 ) THEN
	              ISLASK = 4*(MHQ(IEVA2OUT)/3)
	            ELSE IF (IEVA2OUT.GT.2) THEN
	              ISLASK = 3*(MHQ(IEVA2OUT)/2)
	            ELSE  IF (IEVA2OUT.EQ.2) THEN
	              ISLASK = 2
	            ELSE  IF (IEVA2OUT.EQ.-2) THEN
	              ISLASK = 2 + 1
	            ELSE
	              IF (ISTP2EVA(NEL).EQ.0) THEN
	                ISLASK = 4     ! ieva2out=-4 here
	              ELSE
	                ISLASK = 4 + 1 ! ieva2out=-4 & istp2eva(nel)=-1,-2 here
	              END IF
	            END IF
	            IF (LUNOP.GT.0) THEN
	              SLASK = ISLASK / H01        ! evaluation density
	              NPLUNOP = NPLUNOP + 2
	              WRITE(LUNOP,*) SLEFTLUNOP+X1 , SLASK
	              WRITE(LUNOP,*) SLEFTLUNOP+X0 , SLASK
	            END IF
	            IF ( I2OUT .EQ. 2*(I2OUT/2)+1 ) THEN
	              NFCALLS0 = NFCALLS0 + 4*(MHQ(I2OUT)/3)
	            ELSE IF (I2OUT.GT.4) THEN
	              NFCALLS0 = NFCALLS0 + 3*(MHQ(I2OUT)/2)
	            ELSE IF (I2OUT.EQ.4) THEN
	              IF (IXGAPLOK.EQ.0) THEN
	                NFCALLS0 = NFCALLS0 + 3*(MHQ(I2OUT)/2)
	              ELSE
	                NFCALLS0 = NFCALLS0 + 3*(MHQ(I2OUT)/2) - 2
	              END IF
	            ELSE
	              NFCALLS0 = NFCALLS0 + 2       ! i2out=2 here
	            END IF
	            NFCALLSTST = NFCALLSTST + ISLASK
	            IF ( DABS(H01-MHQ(I)*HCUR) .GT. 0.001*H01 ) STOP 666
	            IF (LUNOUT.GE.90) THEN
	              IF (IFILON.NE.-3) THEN
	                ISLASK = 1 + (NHYD+1)/2  ! nhyd=2,3 if ifilon=+-1,+-2
	              ELSE
	                ISLASK = 1 + (NHYD-1)/2
	              END IF
	              IF ( I.GT.IRUS .AND. JRUS.GE.ISLASK ) THEN
	                IF ( H01 .LE. MHQ(I-1)*HMAX ) THEN
	                  IF (IXGAPLOK.EQ.0) THEN
	                    IF ( I .GT. ISLASK ) THEN
	                      WRITE(LUNOUT,*) '             WHY NOT ',
     $	                      'SUCCESSFUL WITH I=I-1 ?'
	                      NWHY = NWHY + 1
	                    END IF
	                  ELSE
	                    IF ( (I-2) .GT. ISLASK ) THEN
	                      WRITE(LUNOUT,*) '             WHY NOT ',
     $	                      'SUCCESSFUL WITH I=I-1 ?'
	                      NWHY = NWHY + 1
	                    END IF
	                  END IF
	                END IF
	              END IF
	              IF ( X0.NE.SA .OR. X1.NE.SB ) THEN
	                IF ( 2*H01 .LE. MHQ(NRUS)*HMAX ) THEN
	                  WRITE(LUNOUT,*) '             WHY NOT ',
     $	                  'SUCCESSFUL WITH THE PARENT SUBINTERVAL ?'
	                  NWHY = NWHY + 1
	                END IF
	              END IF
	            END IF
	          END IF
	          IADR0 = IA_VALTERM(0)
	          DO IHR=1,NHYDLI
	            VALUE(IHR) = VALUE(IHR) + CWORK(IADR0+IHR)
	          END DO
	          EPSMAGN = EPSMAGN + HSHARE
	          NEL = NEL - 1      ! pop "subinterval stack"
	          IF (NEL.GE.1) THEN
	            GOTO 310
	          ELSE
	            GOTO 200
	          END IF

	        ELSE IF (HCUR.LT.HKVAR) THEN

	          IF (DIMAG(UPH0).EQ.0.0D0) THEN
	            IF (DIMAG(UPH1).EQ.0.0D0) THEN
	              CSLASK = UPH0 + X0*UPRIKT
	              IF (DREAL(CSLASK).GT.0.0D0) THEN
	                CSLASK1 = UPH0 + X1*UPRIKT
	                IF (DREAL(CSLASK1).GT.0.0D0) THEN
	                  IF (NKVAR.EQ.0) THEN
	                    JLKVAR = JL
	                    NKVAR = NKVAR + 1
	                    UPHKVAR0 = CSLASK
	                    UPHKVAR1 = CSLASK1
	                    HSHAREKVAR = HSHARE
	                  ELSE IF (JL.EQ.JLKVAR) THEN
	                    IF ( CSLASK1 .EQ. UPHKVAR0 ) THEN
	                      UPHKVAR0 = CSLASK
	                      HSHAREKVAR = HSHAREKVAR + HSHARE
	                    ELSE
	                      WRITE(LUNKVAR,*) DREAL(UPHKVAR0),
     $	                      DREAL(UPHKVAR1),HSHAREKVAR
	                      NKVAR = NKVAR + 1
	                      UPHKVAR0 = CSLASK
	                      UPHKVAR1 = CSLASK1
	                      HSHAREKVAR = HSHARE
	                    END IF
	                  ELSE
	                    WRITE(LUNKVAR,*) DREAL(UPHKVAR0),
     $	                    DREAL(UPHKVAR1),HSHAREKVAR
	                    JLKVAR = JL
	                    NKVAR = NKVAR + 1
	                    UPHKVAR0 = CSLASK
	                    UPHKVAR1 = CSLASK1
	                    HSHAREKVAR = HSHARE
	                  END IF
	                  IF (LUNOUT.GT.0) THEN
	                    I2OUT = I            ! i2out = 2,3,4,5,..,nrus
	                    IEVA2OUT = IEVA(NEL) !ieva2out=-4 ,-2, 2,3,..,nrus
	                    WRITE(LUNOUT,*) '   % ', SNGL(X0), 
     $	                    SNGL(HCUR), I2OUT, SNGL(DMAX1(ARGLEFT,
     $	                    ARGRIGHT)), IEVA2OUT,' %KVAR'
	                    IF ( IEVA2OUT .EQ. 2*(IEVA2OUT/2)+1 ) THEN
	                      ISLASK = 4*(MHQ(IEVA2OUT)/3)
	                    ELSE IF (IEVA2OUT.GT.2) THEN
	                      ISLASK = 3*(MHQ(IEVA2OUT)/2)
	                    ELSE  IF (IEVA2OUT.EQ.2) THEN
	                      ISLASK = 2
	                    ELSE  IF (IEVA2OUT.EQ.-2) THEN
	                      ISLASK = 2 + 1
	                    ELSE
	                      IF (ISTP2EVA(NEL).EQ.0) THEN
	                        ISLASK = 4     ! ieva2out=-4 here
	                      ELSE
	                        ISLASK = 4 + 1 ! ieva2out=-4 & istp2eva()=-1,-2
	                      END IF
	                    END IF
	                    N2OUT = ISLASK - 1
	                    WRITE(LUNOUT,*) '                      %% ',
     $	                    N2OUT,' interior [Beval %%'
	                    IF (LUNOP.GT.0) THEN
	                      SLASK = ISLASK / H01        ! evaluation density
	                      NPLUNOP = NPLUNOP + 2
	                      WRITE(LUNOP,*) SLEFTLUNOP+X1 , -SLASK
	                      WRITE(LUNOP,*) SLEFTLUNOP+X0 , -SLASK
	                    END IF
	                    NFCALLS0 = NFCALLS0 + 1
	                    NFCALLSTST = NFCALLSTST + ISLASK
	                    IF ( DABS(H01-MHQ(I)*HCUR) .GT. 0.001*H01 ) 
     $                      STOP 666
	                  END IF
	                  EPSMAGN = EPSMAGN + HSHARE   ! "included in advance"
	                  NEL = NEL - 1      ! pop "subinterval stack"
	                  IF (NEL.GE.1) THEN
	                    GOTO 310
	                  ELSE
	                    GOTO 200
	                  END IF
	                END IF
	              END IF
	            END IF
	          END IF

	        END IF
 
	        IF (ARGLRMAX.GE.0.0D0) THEN

	          IF ( IXGAPLOK.EQ.0 .OR. I.NE.4 ) THEN
	            ISLASK = (NRUS-I+1)/2
	          ELSE
	            ISLASK = (NRUS-I+2)/2
	          END IF
	          IF ( ARGLRMAX .GT. MHQ(2*ISLASK)*ARGDIFLIM0 ) 
     $	          GOTO 390 ! "early split" (because of the final "argdiflim0 
C	              requirement") , note that  mhq(2*islask) = 2**islask

	          IF (I.GE.IRUSP2) THEN
	            IF (ARGLRMAX.GT.ARGDIFLIM2) GOTO 390        ! "early split"
	            IF (ARGLEFT.GT.ARGDIFLIM0) THEN
	              IF (ARGRIGHT.LT.ARGDIFLIM1) GOTO 390   ! "early split"
	            END IF
	            IF (ARGRIGHT.GT.ARGDIFLIM0) THEN
	              IF (ARGLEFT.LT.ARGDIFLIM1) GOTO 390    ! "early split"
	            END IF
	          END IF

	        END IF

 400	      CONTINUE
 
C	      after completing the "400 loop" :   ieva(nel) = nrus >= 2 ,   
C	      istp1eva(nel) = 1 ,   istp2eva(nel) = "1,2,0" if nrus 
C	      is "odd,even>2,2"
C	      if  ixgaplok=1 & nrus=4 , however, ieva(nel)=-4 & 
C	      istp1eva(nel)=1 & istp2eva(nel)=0

 390	      CONTINUE
C	      now  ieva(nel) = -4, -2, 2,3,...,nrus  so  iabs(ieva(nel)) >= 2  
C	        istp2eva(nel) < 0  is possible only when  ieva(nel) = -2,-4
C	      values at  x0 + ii*(h01/nst1)  for  
C	        " ii = 0,1,2,...,nst1  &  ii a multiple of istp1eva(nel) " 
C	        are now stored in  "st1..(..,ib)"  where ib=ib0+ii
C	      values at  x0 + ii*(h01/nst2)  for  
C	        " ii = 1,2, 4,5, 7,8, ...,nst2-1  &  ii a positive multiple 
C	          of istp2eva(nel) " 
C	        are now stored in  "st2..(..,ib)"  where ib=ib0+ii-(ii/3) 
C	              ( so the ib:s are among ib0+1,ib0+2,ib0+3,...,ib0+nst1 )
C	        ( if  istp2eva(nel) = 0  no storage here )
C	        if  istp2eva(nel) = -1  storage for  ii = 2*nst2/3
C	        if  istp2eva(nel) = -2  storage for  ii = nst2/3

C	      test for "back"
	      IF (NEL.EQ.NELMAX) THEN
	        IF (NBACK.EQ.NBACKMAX) THEN
	          IF (LUNOUT.GT.0) THEN
	            WRITE(LUNOUT,*)
	            WRITE(LUNOUT,*) '     %%   noconv ;  nfcalls ',
     $	            NFCALLS,'   %%'
	          END IF
	          STOP 777
	        ELSE
	          NBACK = NBACK + 1    
	          SB = X1
	          IF (NEL.GT.1) THEN
C	            deepest subinterval in the stack will be wasted
	            SA = XST(1)
	          ELSE
C	            left half of the subinterval in the stack will be wasted
	            SA = 0.5D0 * ( X0 + X1 )
	          END IF
	          IF (LUNOUT.GT.0) THEN
	            IEV = IEVA(1)   ! iev = -4, -2,-1, 1,2,3,4,5,..,nrus here
C	                            ( iev=-1,1 can only occur when  nel > 1  )
C	            in general when  nel > 1 :   iev = nrus-2  if  nrus >= 5
C	                                         iev = -2      if  nrus  = 4
C	                                         iev = -1      if  nrus  = 3
C	                                         iev =  1      if  nrus  = 2
C	                       however, if  ixgaplok = 1  iev=2 if nrus=4
C	            in general when  nel = 1 :   iev = nrus
C	                       however, if  ixgaplok = 1  iev=-4 if nrus=4
C	            (this is unless the "400 loop" was left prematurely)
	            IF (NEL.GT.1) THEN
	              IF ( IABS(IEV) .EQ. 2*(IABS(IEV)/2)+1 ) THEN
	                IF (IEV.GT.1) THEN
	                  N2OUT = 4*(MHQ(IEV)/3)
	                ELSE IF (IEV.EQ.1) THEN
	                  N2OUT = 1
	                ELSE
	                  N2OUT = 1 + 1     ! iev=-1 here
	                END IF
	              ELSE IF (IEV.GT.2) THEN
	                N2OUT = 3*(MHQ(IEV)/2)
	              ELSE  IF (IEV.EQ.2) THEN
	                N2OUT = 2
	              ELSE IF (IEV.EQ.-2) THEN
	                N2OUT = 2 + 1
	              ELSE IF (ISTP2EVA(1).EQ.0) THEN
	                N2OUT = 4       ! iev=-4 here
	              ELSE
	                N2OUT = 4 + 1   ! iev=-4 & istp2eva(1)=-1,-2 here
	              END IF
	            ELSE
	              IF ( IEV .EQ. 2*(IEV/2)+1 ) THEN
	                N2OUT = 2*(MHQ(IEV)/3)
	              ELSE IF (IEV.GT.2) THEN
	                N2OUT = 3*(MHQ(IEV)/4)
	              ELSE  IF (IEV.EQ.2) THEN
	                N2OUT = 1
	              ELSE IF (IEV.EQ.-2) THEN
	                IF (ISTP2EVA(1).EQ.-1) THEN
	                  N2OUT = 1
	                ELSE
	                  N2OUT = 1 + 1   !  istp2eva(1)=-2  here
	                END IF
	              ELSE
	                IF (ISTP2EVA(1).NE.-2) THEN
	                  N2OUT = 2      !  iev=-4, istp2eva(1)=0,-1  here
	                ELSE
	                  N2OUT = 2 + 1  !  iev=-4 here
	                END IF
	              END IF
	            END IF
	            WRITE(LUNOUT,*) '   sa,sb ',SNGL(SA),SNGL(SB),
     $	            '   %% back, ',N2OUT,' wasted eval %%'
	            NFCALLSTST = NFCALLSTST + N2OUT
	          END IF
	          IF (NEL.GT.1) THEN
C	            shift the  xst, ieva,istp1eva,istp2eva  arrays
	            DO J=1,NEL
	              XST(J-1) = XST(J)
	            END DO
	            DO J=2,NEL
	              IEVA(J-1) = IEVA(J)
	              ISTP1EVA(J-1) = ISTP1EVA(J)
	              ISTP2EVA(J-1) = ISTP2EVA(J)
	            END DO
	            IF (IABS(IFILON).EQ.2) THEN
C	              shift the rfaktpos arrays backwards
	              IADR0 = IA_RFAKTPOS(0,0)
	              ISLASK = IADR0 + NEL*LI
	              DO M = IADR0+1, ISLASK
	                CWORK(M) = CWORK(M+LI)
	              END DO
	            END IF
	          END IF
C	          shift stored "cdet arguments" and "evaluations" nst1/2 steps
	          ISLASK = NEL*NST1H
	          DO IBQ = 0, ISLASK
	            IB = IBQ + NST1H
	            ARG1CDET(IBQ) = ARG1CDET(IB)
	            IADR0 = IA_ST1DH(0,IB)
	            IADRQ = IA_ST1DH(0,IBQ)
	            DO IH=1,NHYD
	              CWORK(IADRQ+IH) = CWORK(IADR0+IH)
	            END DO
	            IADR0 = IA_ST1R(0,IB)
	            IADRQ = IA_ST1R(0,IBQ)
	            DO IR=1,LI
	              CWORK(IADRQ+IR) = CWORK(IADR0+IR)
	            END DO
	          END DO
	          DO IBQ = 1, ISLASK
	            IB = IBQ + NST1H
	            ARG2CDET(IBQ) = ARG2CDET(IB)
	            IADR0 = IA_ST2DH(0,IB)
	            IADRQ = IA_ST2DH(0,IBQ)
	            DO IH=1,NHYD
	              CWORK(IADRQ+IH) = CWORK(IADR0+IH)
	            END DO
	            IADR0 = IA_ST2R(0,IB)
	            IADRQ = IA_ST2R(0,IBQ)
	            DO IR=1,LI
	              CWORK(IADRQ+IR) = CWORK(IADR0+IR)
	            END DO
	          END DO
C	          can now get  nel = 0  temporarily
	          NEL = NEL - 1      ! pop "subinterval stack"
	        END IF
	      END IF

C	      subinterval splitting
	      NEL1 = NEL + 1
	      XM = 0.5D0 * ( X0 + X1 )
	      XST(NEL) = XM   ! the "old" subinterval becomes its left half
	      XST(NEL1) = X1  ! the "new" subinterval, right half of the "old"
	      IF (NEL.GT.0) THEN
	        IEV = IEVA(NEL)    !  iev here:  -4, -2, 2, 3,4, 5,6,7,..,nrus
C	        unless the "400 loop" was left prematurely,  iev = nrus (or 
C	        -4 if "ixgaplok=1 & nrus=4")
	        ISTP1EVA(NEL1) = ISTP1EVA(NEL)
	        ISTP2 = ISTP2EVA(NEL)
	        IF (IEV.GE.5) THEN
	          IEV = IEV - 2
	          IEVA(NEL) = IEV
	          IEVA(NEL1) = IEV
	          ISTP2EVA(NEL1) = ISTP2  ! istp2 here: (nst2/3)/2,...,8,4,2,1
	        ELSE IF ( IEV.EQ.3 .OR. IEV.EQ.4 ) THEN
	          IEV = 2 - IEV  ! so  3 -> -1  4 -> -2  ,  istp2 = nst2/3 here
	          IEVA(NEL) = IEV
	          IEVA(NEL1) = IEV
	          ISTP2EVA(NEL) = -1
	          ISTP2EVA(NEL1) = -2
	        ELSE IF (IEV.EQ.2) THEN
	          IEVA(NEL) = 1
	          IEVA(NEL1) = 1
	          ISTP2EVA(NEL1) = 0   !  istp2 = 0  here
	        ELSE IF (IEV.EQ.-2) THEN
	          IF (ISTP2.EQ.-1) THEN
	            IEVA(NEL) = 1
	            IEVA(NEL1) = -1
	            ISTP2EVA(NEL) = 0
	            ISTP2EVA(NEL1) = -2
	          ELSE
	            IEVA(NEL) = -1 !  ( istp2 = -2  in this case )
	            IEVA(NEL1) = 1
	            ISTP2EVA(NEL) = -1
	            ISTP2EVA(NEL1) = 0
	          END IF
	        ELSE
	          IF (ISTP2.EQ.0) THEN
	            IEVA(NEL) = 2  !  ( iev = -4  in this case )
	            IEVA(NEL1) = 2
	            ISTP2EVA(NEL) = 0
	            ISTP2EVA(NEL1) = 0
	          ELSE IF (ISTP2.EQ.-1) THEN
	            IEVA(NEL) = 2  !  ( iev = -4  in this case )
	            IEVA(NEL1) = -2
	            ISTP2EVA(NEL) = 0
	            ISTP2EVA(NEL1) = -2
	          ELSE
	            IEVA(NEL) = -2  !  ( iev=-4 and istp2=-2 in this case )
	            IEVA(NEL1) = 2
	            ISTP2EVA(NEL) = -1
	            ISTP2EVA(NEL1) = 0
	          END IF
	        END IF
	        IF (IABS(IFILON).EQ.2) THEN
	          IADR0 = IA_RFAKTPOS(0,NEL)
	          IADR1 = IA_RFAKTPOS(0,NEL1)
	          DO IR=1,LI
	            CWORK(IADR1+IR) = CWORK(IADR0+IR)
	          END DO
	        END IF
	      ELSE
C	        now  nel = 0  temporarily  ( from "back" with  nelmax = 1 )
	        IEV = IEVA(1)     !  iev here:  -4, -2, 2, 3,4, 5,6,7,..,nrus
C	        unless the "400 loop" was left prematurely,  iev = nrus (or 
C	        -4 if "ixgaplok=1 & nrus=4")
	        ISTP2 = ISTP2EVA(1)
	        IF (IEV.GE.5) THEN
	          IEVA(1) = IEV - 2       ! istp2 here: (nst2/3)/2,...,8,4,2,1
	        ELSE IF ( IEV.EQ.3 .OR. IEV.EQ.4 ) THEN
	          IEVA(1) = 2 - IEV ! so  3 -> -1  4 -> -2 ,  istp2=nst2/3 here
	          ISTP2EVA(1) = -2
	        ELSE IF (IEV.EQ.2) THEN
	          IEVA(1) = 1           !  istp2 = 0  here
	        ELSE IF (IEV.EQ.-2) THEN
	          IF (ISTP2.EQ.-1) THEN
	            IEVA(1) = -1
	            ISTP2EVA(1) = -2
	          ELSE
	            IEVA(1) = 1   !  ( istp2 = -2  in this case )
	            ISTP2EVA(1) = 0
	          END IF
	        ELSE
	          IF (ISTP2.EQ.0) THEN
	            IEVA(1) = 2  !  ( iev = -4  in this case )
	            ISTP2EVA(1) = 0
	          ELSE IF (ISTP2.EQ.-1) THEN
	            IEVA(1) = -2  !  ( iev = -4  in this case )
	            ISTP2EVA(1) = -2
	          ELSE
	            IEVA(1) = 2  !  ( iev=-4 and istp2=-2 in this case )
	            ISTP2EVA(1) = 0
	          END IF
	        END IF
	      END IF
	      IF (IABS(IFILON).EQ.2) THEN
	        UP = UPH0 + XM*UPRIKT
	        IADR0 = IA_RFAKTPOS(1,NEL)
	        CALL RCOMPPOS ( UP, CWORK(IADR0) )
	      END IF
	      NEL = NEL1       ! push "subinterval stack"

C	      "set ib0 and expand" for the top subinterval in the stack
 310	      IB0 = (NEL-1)*NST1H
C	      ( even multiples, at least some, of the new HMIN1 )
	      ISTP1 = ISTP1EVA(NEL)
	      DO I = NST1H, 1, -ISTP1
	        IB = IB0 + I
	        IBQ = IB + I
C	        ib:s  among ib0+nst1h,ib0+nst1h-1,...,ib0+2,ib0+1
C	        ibq:s among ib0+nst1 ,ib0+nst1-2 ,...,ib0+4,ib0+2
	        ARG1CDET(IBQ) = ARG1CDET(IB)
	        IADR0 = IA_ST1DH(0,IB)
	        IADRQ = IA_ST1DH(0,IBQ)
	        DO IH=1,NHYD
	          CWORK(IADRQ+IH) = CWORK(IADR0+IH)
	        END DO
	        IADR0 = IA_ST1R(0,IB)
	        IADRQ = IA_ST1R(0,IBQ)
	        DO IR=1,LI
	          CWORK(IADRQ+IR) = CWORK(IADR0+IR)
	        END DO
	      END DO
	      ISTP1EVA(NEL) = 2*ISTP1
C	      ( even non-three multiples, at least some, of the new HMIN2 )
	      ISTP2 = ISTP2EVA(NEL)
	      IF (ISTP2.GT.0) THEN
C	        now  ieva(nel) >= 3  so  nrus >= 3  and  nst2 is even
C	        istp2 here: nst2/3,...,8,4,2,1
	        IISTEP = 3*ISTP2
	        ISTEP = 2*ISTP2
	        IF ( ISTP2 .EQ. 3*(ISTP2/3)+1 ) THEN
	          M21 = 1
	          M22 = 0
	        ELSE
	          M21 = 0
	          M22 = 1
	        END IF
	        II21 = NST2H - ISTP2
	        I21 = II21 - (II21/3)
	        II22 = NST2H - ISTEP
	        I22 = II22 - (II22/3)
	        DO WHILE (II21.GT.0)
	          IB = IB0 + I21
	          IBQ = IB + I21 - M21
	          ARG2CDET(IBQ) = ARG2CDET(IB)
	          IADR0 = IA_ST2DH(0,IB)
	          IADRQ = IA_ST2DH(0,IBQ)
	          DO IH=1,NHYD
	            CWORK(IADRQ+IH) = CWORK(IADR0+IH)
	          END DO
	          IADR0 = IA_ST2R(0,IB)
	          IADRQ = IA_ST2R(0,IBQ)
	          DO IR=1,LI
	            CWORK(IADRQ+IR) = CWORK(IADR0+IR)
	          END DO
	          II21 = II21 - IISTEP
	          I21 = I21 - ISTEP
	          IF (II22.LE.0) GOTO 22
	          IB = IB0 + I22
	          IBQ = IB + I22 - M22
	          ARG2CDET(IBQ) = ARG2CDET(IB)
	          IADR0 = IA_ST2DH(0,IB)
	          IADRQ = IA_ST2DH(0,IBQ)
	          DO IH=1,NHYD
	            CWORK(IADRQ+IH) = CWORK(IADR0+IH)
	          END DO
	          IADR0 = IA_ST2R(0,IB)
	          IADRQ = IA_ST2R(0,IBQ)
	          DO IR=1,LI
	            CWORK(IADRQ+IR) = CWORK(IADR0+IR)
	          END DO
	          II22 = II22 - IISTEP
	          I22 = I22 - ISTEP
	        END DO
  22	        ISTP2EVA(NEL) = 2*ISTP2
	      ELSE IF (ISTP2.EQ.-1) THEN
C	        now  nrus >= 3  and  ieva(nel) = -1,-2, -4
	        IB = IB0 + NST1H - (NST1H/3)  ! ib0+ii-(ii/3) with ii=nst2/3
	        IBQ = IB0 + NST1 - (NST1/3)   ! ib0+ii-(ii/3) with ii=2*nst2/3
	        ARG2CDET(IBQ) = ARG2CDET(IB)
	        IADR0 = IA_ST2DH(0,IB)
	        IADRQ = IA_ST2DH(0,IBQ)
	        DO IH=1,NHYD
	          CWORK(IADRQ+IH) = CWORK(IADR0+IH)
	        END DO
	        IADR0 = IA_ST2R(0,IB)
	        IADRQ = IA_ST2R(0,IBQ)
	        DO IR=1,LI
	          CWORK(IADRQ+IR) = CWORK(IADR0+IR)
	        END DO
	      ELSE IF (ISTP2.EQ.-2) THEN
	        IF (NST2.GE.6) THEN
C	          ( with  nst2 = 3  nothing needs to be done here, one 
C	          would get  ibq = ib (=ib0+1) )
C	          now   nrus >= 5  and  ieva(nel) = -1,-2, -4
	          IB = IB0 + NST1H - (NST1/3)    ! (ib0-nst1h)+ii-(ii/3) with 
C	          ii=2*nst2/3
	          IBQ = IB0 + NST1H - (NST1H/3)  ! ib0+ii-(ii/3) with ii=nst2/3
	          ARG2CDET(IBQ) = ARG2CDET(IB)
	          IADR0 = IA_ST2DH(0,IB)
	          IADRQ = IA_ST2DH(0,IBQ)
	          DO IH=1,NHYD
	            CWORK(IADRQ+IH) = CWORK(IADR0+IH)
	          END DO
	          IADR0 = IA_ST2R(0,IB)
	          IADRQ = IA_ST2R(0,IBQ)
	          DO IR=1,LI
	            CWORK(IADRQ+IR) = CWORK(IADR0+IR)
	          END DO
	        END IF
	      END IF

 300	    CONTINUE

 200	  CONTINUE

 100	CONTINUE

	IF (LUNOUT.GT.0) THEN
	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' -------------------------'
	  WRITE(LUNOUT,*) ' for a subinterval :  ( length=mhq(i)*hcur ',
     $	  ', mhq(i)="parts" )'
	  WRITE(LUNOUT,*) '   first row "i(eva)", second row "parts", ',
     $	  'third row "evaluations"'
	  WRITE(LUNOUT,*) '     (1) 2 3 4 5  6  7  8  9 10 11 12  13 ',
     $	  ' 14  15  16  17 ..'
	  WRITE(LUNOUT,*) '     (1) 2 3 4 6  8 12 16 24 32 48 64  96 ',
     $	  '128 192 256 384 ..'
	  WRITE(LUNOUT,*) '     (1) 2 4 6 8 12 16 24 32 48 64 96 128 ',
     $	  '192 256 384 512 ..'
	  WRITE(LUNOUT,*) '   i(eva)>=3 odd,even  -->  4/3,3/2 evalua',
     $	  'tions per part'
	  WRITE(LUNOUT,*) '   "-" means 1 more and 2 (or 1) less ',
     $	  'eval for ieva=-1,-2 and -4 , resp'
	  WRITE(LUNOUT,*) ' in total: one starting-point evaluation also'
	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' nfcalls0,nfcalls ',NFCALLS0,NFCALLS
	  WRITE(LUNOUT,*) '       ( nfcalls0 << nfcalls  indicates bad ',
     $	  'economy )' !remedy: reduce "back waste", improve "early split" logic
C	  nfcalls0 is based on "i" rather than "ieva" and does not include 
C	  "back waste"
	  IF (NFCALLS.NE.NFCALLSTST) THEN
	    WRITE(LUNOUT,*) '       NFCALLSTST=',NFCALLSTST
	    STOP 555
	  END IF
	  IF (LUNOP.GT.0) THEN
	    WRITE(LUNOUT,*) ' lunop,nplunop ',LUNOP,NPLUNOP
	    WRITE(LUNOP,*)
	    WRITE(LUNOP,*) ' nplunop = ',NPLUNOP
	    CLOSE(LUNOP)
	  END IF
	  IF (LUNOUT.GE.90) THEN
	    IF (NWHY.EQ.0) THEN
	      WRITE(LUNOUT,*) ' extended "self test" performed, ok'
	    ELSE
	      WRITE(LUNOUT,*) ' extended "self test" performed,  ',
     $	      'NWHY =',NWHY
	    END IF
	  END IF
	  WRITE(LUNOUT,*) ' epsmagn ',EPSMAGN
	  IF (HKVAR.GT.0.0D0) WRITE(LUNOUT,*) ' hkvar,lunkvar, nkvar ',
     $	  HKVAR,LUNKVAR,' ',NKVAR
	  IF (NKVAR.GT.0) WRITE(LUNOUT,*) '        %%%  nkvar > 0  %%%'
	  WRITE(LUNOUT,*) ' ok return from adintxfi'
	  WRITE(LUNOUT,*)
	  IF (LUNOUT.NE.6) CLOSE(LUNOUT)
	END IF

	IF (HKVAR.GT.0.0D0) THEN
	  IF (NKVAR.GT.0) WRITE(LUNKVAR,*) DREAL(UPHKVAR0),
     $	  DREAL(UPHKVAR1),HSHAREKVAR
	END IF	  

	RETURN
 
	END







	INTEGER FUNCTION IA_RFAKTPOS (IR,J)
	IMPLICIT NONE
	INTEGER IR,J
	INCLUDE 'IADR0.INC'
	IA_RFAKTPOS = IADR0_RFAKTPOS + J*LICOP + IR
	RETURN
	END

	INTEGER FUNCTION IA_ST1DH (IH,IB)
	IMPLICIT NONE
	INTEGER IH,IB
	INCLUDE 'IADR0.INC'
	IA_ST1DH = IADR0_ST1DH + IB*NHYDCOP + IH
	RETURN
	END

	INTEGER FUNCTION IA_ST1R (IR,IB)
	IMPLICIT NONE
	INTEGER IR,IB
	INCLUDE 'IADR0.INC'
	IA_ST1R = IADR0_ST1R + IB*LICOP + IR
	RETURN
	END

	INTEGER FUNCTION IA_ST2DH (IH,IB)
	IMPLICIT NONE
	INTEGER IH,IB
	INCLUDE 'IADR0.INC'
	IA_ST2DH = IADR0_ST2DH + (IB-1)*NHYDCOP + IH
	RETURN
	END

	INTEGER FUNCTION IA_ST2R (IR,IB)
	IMPLICIT NONE
	INTEGER IR,IB
	INCLUDE 'IADR0.INC'
	IA_ST2R = IADR0_ST2R + (IB-1)*LICOP + IR
	RETURN
	END

	INTEGER FUNCTION IA_SS1 (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR0.INC'
	IA_SS1 = IADR0_SS1 + IHR
	RETURN
	END

	INTEGER FUNCTION IA_SS2 (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR0.INC'
	IA_SS2 = IADR0_SS2 + IHR
	RETURN
	END

	INTEGER FUNCTION IA_SSD (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR0.INC'
	IA_SSD = IADR0_SSD + IHR
	RETURN
	END

	INTEGER FUNCTION IA_SS (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR0.INC'
	IA_SS = IADR0_SS + IHR
	RETURN
	END

	INTEGER FUNCTION IA_TABNEW (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR0.INC'
	IA_TABNEW = IADR0_TABNEW + IHR
	RETURN
	END

	INTEGER FUNCTION IA_VALTERM (IHR)
	IMPLICIT NONE
	INTEGER IHR
	INCLUDE 'IADR0.INC'
	IA_VALTERM = IADR0_VALTERM + IHR
	RETURN
	END

	INTEGER FUNCTION IA_TAB (IHR,J)
	IMPLICIT NONE
	INTEGER IHR,J
	INCLUDE 'IADR0.INC'
	IA_TAB = IADR0_TAB + (J-1)*NHYDLICOP + IHR
	RETURN
	END

	INTEGER FUNCTION IA_YTAB (IHR,J)
	IMPLICIT NONE
	INTEGER IHR,J
	INCLUDE 'IADR0.INC'
	IA_YTAB = IADR0_YTAB + (J-1)*NHYDLICOP + IHR
	RETURN
	END