C	subroutine "zerofind" to find zeros of an analytic function

C	Sven Ivansson 1990






	SUBROUTINE ZEROFIND (CDETSUB,CDETSUBCIRC, UPREAL1,UPREAL2,
     $	UPIMAG1,UPIMAG2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,
     $	DRELMAX,DRELGOALFAKT, ARGDIFSTORE, SLIMLEFT, DLOGREDUCFAKT,
     $	DDELT,XYTOL, IOPTION,ICHECK,LUNLEFT,NZERODIM, NZERO,NZLEFT,  
     $	ZREAL,ZIMAG, XYONEMIN,XYONEMAX, 
     $	NELMAX,NPLINMAX,IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK,
     $	IDZHBACK,ICHKBACK)
C	subroutine corntestsub (x1,x2,y1,y2, skiprect)
C	                     for early removal of uninteresting rectangles 
C	                     subroutine corntestsub must be provided, a 
C	                     dummy version is as 
C	                     -----
C	                     subroutine corntestsub (x1,x2,y1,y2, skiprect)
C	                     implicit none
C	                     real*8 x1,x2,y1,y2   ! rectangle coordinates
C	                     logical skiprect
C	                     skiprect = .false.
C	                     return
C	                     end
C	                     -----
C	common /mittchange/ mchange
C	                     mchange is 0 for "dzeroholo start in the 
C	                     middle of the rectangle" and nonzero for 
C	                     "dzeroholo start with contour-integration 
C	                     estimate"
C	                     note that only the "data stored in argnel & 
C	                     alnnel from the argseg calls" are used in the 
C	                     latter case, the data storage is controlled 
C	                     by "argdifstore", hence it may be very sparse 
C	                     when the argument of the analytic function 
C	                     is constant along come line segment (an extended 
C	                     storage in argseg could of course be possible 
C	                     but has not been implemented)
C	cdetsub              subroutine cdetsub (px,py, rlndet,cdet)  gives 
C	                     "the function" as  dexp(rlndet)*cdet  , "the 
C	                     function" must be analytic inside the rectangle 
C	                     and continuous up to the boundary and zero-free 
C	                     on the boundary
C	cdetsubcirc          subroutine cdetsubcirc (vx,vy, rlndet,cdet)  
C	                     as cdetsub but for another independent variable, 
C	                     the arguments are related by "(px,py) = 
C	                     up0circ + radiecirc*cdexp(i*(vx,vy))" with 
C	                     up0circ & radiecirc from "common /circdata/"
C	upreal1,upreal2      horizontal limits of the rectangle
C	upimag1,upimag2      vertical limits of the rectangle
C	ddsupmin             minimal step that is handled (strict limit),  
C	                     always take  ddsupmin >= 0.0d0  ("argseg")
C	ddsupmax             maximal step ("argseg")
C	ddsupfaktmax         maximal factor for step change ("argseg")
C	derivdrelmax         maximum tolerable "relative difference" in the 
C	                     "derivative estimation at up1" for the "safe 
C	                     start", always take  derivdrelmax < 1.0d0  
C	                     ("argseg")
C	drelmax              maximum allowed relative change of the analytic 
C                            function in each step, always take  
C	                     drelmax < 1.0d0  ("argseg")
C	drelgoalfakt         for choosing "next step", always take  
C	                     drelgoalfakt < 1.0d0  ("argseg")
C	argdifstore          lower limit in degrees for storage, always 
C	                     take  argdifstore < 180.0d0  ("argseg")
C	slimleft	     side limit for stopping rectangle division
C	dlogreducfakt        "dlog" of required reduction factor in each step 
C	                     for the absolute value of "the function"  
C	                     ("dzeroholo")
C	ddelt                distance-difference limit for "secant derivative"
C	                     ("dzeroholo")
C	xytol                desired accuracy for the zeros  ("dzeroholo")
C	ioption              0 for "just nzero", 1 for "normal", 2 for "+ log"
C	icheck               1 for additional "argseg" check
C	lunleft              unit number for writing remaining zeros
C	nzerodim             dimensions of zreal,zimag as declared in calling 
C	                     program
C	nzero              in put/output (OBS)  total number of zeros so far 
C	                           (the new ones are added at the end)
C	zreal(1:nzero)     input/output (OBS)  real parts of the zeros
C	zimag(1:nzero)     input/output (OBS)  imag parts of the zeros
C	xyonemin(1:nzero)  input/output (OBS)  estimated accuracy for the 
C	                           zero, obtained from "xytolut"
C	xyonemax(1:nzero)  input/output (OBS)  distance to border of cell, 
C	                           obtained from "x1,x2,y1,y2"
C	nelmax             input/output (OBS)  maximum current value of nel
C	nplinmax           input/output (OBS)  maximum current value of nplin
C	iasfcall           input/output (OBS)  total number of cdetsub calls 
C	                           in connection with and from argseg
C	idzhfcall          input/output (OBS)  total number of cdetsub calls 
C	                           in connection with and from dzeroholo
C	ichkfcall          input/output (OBS)  total number of cdetsubcirc 
C	                           calls for "the check"
C	inelcompr          input/output (OBS)  total number of "nelcompr" calls
C	iasback            input/output (OBS)  total number of "back steps" 
C	                           upon unsuccessful return from argseg
C	idzhback           input/output (OBS)  total number of "division back 
C	                           steps" upon unsuccessful return from 
C	                           dzeroholo
C	ichkback           input/output (OBS)  total number of "back steps" 
C	                           upon "check failure"

	IMPLICIT NONE

	INTEGER NPLINDIM, NELDIM
	PARAMETER ( NPLINDIM = 3000, NELDIM = 30 )
C       nplin <= nplindim ,  nel <= neldim

	INCLUDE '../cpm/CIRCDATA.INC'            ! with up0circ and radiecirc

	COMMON /MITTCHANGE/ MCHANGE
	INTEGER MCHANGE
	DATA MCHANGE /1/

	INTEGER IOPTION,ICHECK,LUNLEFT,NZERODIM
	REAL*8 UPREAL1,UPREAL2,UPIMAG1,UPIMAG2, DDSUPMIN,DDSUPMAX,
     $	DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,DRELGOALFAKT, ARGDIFSTORE, 
     $	SLIMLEFT, DLOGREDUCFAKT,DDELT,XYTOL

	INTEGER NZERO,NZLEFT, NELMAX,NPLINMAX,IASFCALL,IDZHFCALL,
     $	ICHKFCALL,INELCOMPR,IASBACK,IDZHBACK,ICHKBACK
	REAL*8 ZREAL(*), ZIMAG(*), XYONEMIN(*), XYONEMAX(*)

	INTEGER NPLIN, IASFAIL,IDZHFAIL, NROOT, NEL,NEL1, 
     $	NPLOK1,NPLOK2,NPLOK3,NPLOK4, INDEX, I, 
     $	NPNEL(4,neldim)
	REAL ARGNEL(0:nplindim,4,neldim), ALNNEL(0:nplindim,4,neldim), 
     $	DDSUPNEL(nplindim,4,neldim)     !!!
	REAL*8 DDSUP1,ARGDIFTEST, X1,X2,Y1,Y2, PX,PY, XYTOLUT, 
     $	XLOK,YLOK, XLOK2,YLOK2, ARGSLASK, SLASK, SLASK1,SLASK2, PI, 
     $	RLNDET,RLNDET0, X1ONE,X2ONE,Y1ONE,Y2ONE, 
     $	ARGLIN(0:nplindim), ALNLIN(0:nplindim), DDSUPLIN(nplindim), 
     $	XYNEL(0:nplindim,4,neldim)
	COMPLEX*16 UP1,UP2, CDET,CDET0, 
     $	UPLIN(0:nplindim)
	LOGICAL ONEROOT, SKIPRECT

        EXTERNAL CDETSUB, CDETSUBCIRC

	DATA PI /3.141592653589793D0/

	CALL CORNTESTSUB (UPREAL1,UPREAL2,UPIMAG1,UPIMAG2, SKIPRECT)
	IF (SKIPRECT) RETURN

	IF (IOPTION.EQ.2) THEN
	  WRITE(99,*)
	  WRITE(99,*) ' --------------------------------------------'
	  WRITE(99,*)
	  WRITE(99,*) ' ZEROFIND log'
	END IF

	ARGDIFTEST = 0.5D0*ARGDIFSTORE + DASIND(DRELMAX)

C	---------------------------------------------------
C	  a "nel" rectangle     !       attributes for
C	with its side indices   !    the "nel" rectangles
C	          2            /
C	   --------------     /            npnel(1:4,1:nel) , numbers of points
C	   !            !    /                    (excluding the initial ones)
C	  3!            !4  !  xynel (0:npnel(,),1:4,1:nel) , x/y (increasing)
C	   !            !   ! argnel (0:npnel(,),1:4,1:nel) , arg(cdet) degrees
C	   !            !   ! alnnel (0:npnel(,),1:4,1:nel) , tlnsu+ln(a(cdet))
C	   --------------  / ddsupnel (1:npnel(,),1:4,1:nel) , "ddsupinf"
C	          1       /
C	                 ! more compact storage could be achieved using 
C	                 ! "pointers n0nel(1:4,1:nel)" to one-dimensional 
C	                 ! arrays xynel,argnel,alnnel,ddsupnel
C	                 ! multiple storage for the same "parent side" for 
C	                 ! different "nel" rectangles could of course be 
C	                 ! avoided by using suitable pointers
C	---------------------------------------------------

C	---------------------------------------------------
C	initial round

	NEL = 1
	IF (NEL.GT.NELMAX) NELMAX = NEL
	UP1 = DCMPLX ( UPREAL1, UPIMAG1 )
	IASFCALL = IASFCALL + 1
	CALL CDETSUB (UPREAL1,UPIMAG1, RLNDET0,CDET0)
	UP2 = DCMPLX ( UPREAL2, UPIMAG1 )
	DDSUP1 = 0.0D0
	CALL ARGSEG (CDETSUB, UP1,UP2, RLNDET0,CDET0, DDSUP1, 
     $	DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, 
     $	NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET,CDET, 
     $	IASFCALL,IASFAIL)
	IF (NPLIN.GT.NPLINMAX) NPLINMAX = NPLIN
	IF (IASFAIL.NE.0) THEN
          WRITE(*,*) ' %%% ARGSEG FAILURE ON "BIG PART 1" %%% '
	  WRITE(*,*) ' iasfail=',IASFAIL
	  WRITE(*,*) UP1
	  WRITE(*,*) UP2
	  STOP 316
	END IF
	NPNEL(1,1) = NPLIN
	DO I = 0, NPLIN
	  XYNEL(I,1,1) = DREAL(UPLIN(I))
	END DO
	DO I = 0, NPLIN
	  ARGNEL(I,1,1) = ARGLIN(I)
	END DO
	DO I = 0, NPLIN
	  ALNNEL(I,1,1) = ALNLIN(I)
	END DO
	DO I = 1, NPLIN
	  DDSUPNEL(I,1,1) = DDSUPLIN(I)
	END DO
	UP1 = UP2
	UP2 = DCMPLX ( UPREAL2, UPIMAG2 )
	DDSUP1 = DDSUPLIN(NPLIN)
	CALL ARGSEG (CDETSUB, UP1,UP2, RLNDET,CDET, DDSUP1, 
     $	DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, 
     $	NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET,CDET, 
     $	IASFCALL,IASFAIL)
	IF (NPLIN.GT.NPLINMAX) NPLINMAX = NPLIN
	IF (IASFAIL.NE.0) THEN
          WRITE(*,*) ' %%% ARGSEG FAILURE ON "BIG PART 4" %%% '
	  WRITE(*,*) ' iasfail=',IASFAIL
	  WRITE(*,*) UP1
	  WRITE(*,*) UP2
	  STOP 316
	END IF
	NPNEL(4,1) = NPLIN
	DO I = 0, NPLIN
	  XYNEL(I,4,1) = DIMAG(UPLIN(I))
	END DO
	DO I = 0, NPLIN
	  ARGNEL(I,4,1) = ARGLIN(I)
	END DO
	DO I = 0, NPLIN
	  ALNNEL(I,4,1) = ALNLIN(I)
	END DO
	DO I = 1, NPLIN
	  DDSUPNEL(I,4,1) = DDSUPLIN(I)
	END DO
	UP1 = DCMPLX ( UPREAL1, UPIMAG1 )
	UP2 = DCMPLX ( UPREAL1, UPIMAG2 )
	DDSUP1 = DDSUPNEL(1,1,1)
	CALL ARGSEG (CDETSUB, UP1,UP2, RLNDET0,CDET0, DDSUP1, 
     $	DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, 
     $	NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET,CDET, 
     $	IASFCALL,IASFAIL)
	IF (NPLIN.GT.NPLINMAX) NPLINMAX = NPLIN
	IF (IASFAIL.NE.0) THEN
          WRITE(*,*) ' %%% ARGSEG FAILURE ON "BIG PART 3" %%% '
	  WRITE(*,*) ' iasfail=',IASFAIL
	  WRITE(*,*) UP1
	  WRITE(*,*) UP2
	  STOP 316
	END IF
	NPNEL(3,1) = NPLIN
	DO I = 0, NPLIN
	  XYNEL(I,3,1) = DIMAG(UPLIN(I))
	END DO
	DO I = 0, NPLIN
	  ARGNEL(I,3,1) = ARGLIN(I)
	END DO
	DO I = 0, NPLIN
	  ALNNEL(I,3,1) = ALNLIN(I)
	END DO
	DO I = 1, NPLIN
	  DDSUPNEL(I,3,1) = DDSUPLIN(I)
	END DO
	UP1 = UP2
	UP2 = DCMPLX ( UPREAL2, UPIMAG2 )
	DDSUP1 = DDSUPLIN(NPLIN)
	CALL ARGSEG (CDETSUB, UP1,UP2, RLNDET,CDET, DDSUP1, 
     $	DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, 
     $	NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET,CDET, 
     $	IASFCALL,IASFAIL)
	IF (NPLIN.GT.NPLINMAX) NPLINMAX = NPLIN
	IF (IASFAIL.NE.0) THEN
          WRITE(*,*) ' %%% ARGSEG FAILURE ON "BIG PART 2" %%%'
	  WRITE(*,*) ' iasfail=',IASFAIL
	  WRITE(*,*) UP1
	  WRITE(*,*) UP2
	  STOP 316
	END IF
	NPNEL(2,1) = NPLIN
	DO I = 0, NPLIN
	  XYNEL(I,2,1) = DREAL(UPLIN(I))
	END DO
	DO I = 0, NPLIN
	  ARGNEL(I,2,1) = ARGLIN(I)
	END DO
	DO I = 0, NPLIN
	  ALNNEL(I,2,1) = ALNLIN(I)
	END DO
	DO I = 1, NPLIN
	  DDSUPNEL(I,2,1) = DDSUPLIN(I)
	END DO

C	---------------------------------------------------
C	division round

	ONEROOT = .FALSE.

 100	NPLOK1 = NPNEL(1,NEL)
	NPLOK2 = NPNEL(2,NEL)
	NPLOK3 = NPNEL(3,NEL)
	NPLOK4 = NPNEL(4,NEL)
	X1 = XYNEL(0,1,NEL)
	X2 = XYNEL(NPLOK1,1,NEL)
	Y1 = XYNEL(0,3,NEL)
	Y2 = XYNEL(NPLOK3,3,NEL)
	CALL CORNTESTSUB (X1,X2,Y1,Y2, SKIPRECT)
	IF (SKIPRECT) THEN
	  NEL = NEL - 1
	  IF (NEL.EQ.0) RETURN
	  GOTO 100
	END IF
	ARGSLASK = ( + ARGNEL(NPLOK1,1,NEL) - ARGNEL(0,1,NEL) 
     $	             - ARGNEL(NPLOK2,2,NEL) + ARGNEL(0,2,NEL) 
     $	             - ARGNEL(NPLOK3,3,NEL) + ARGNEL(0,3,NEL) 
     $	             + ARGNEL(NPLOK4,4,NEL) - ARGNEL(0,4,NEL) ) / 
     $	           360.0D0
	NROOT = IDNINT (ARGSLASK)
	IF (IOPTION.EQ.2) THEN
	  WRITE(99,*)
	  WRITE(99,*) ' nel,nroot ',NEL,NROOT
	  WRITE(99,*) ' x1,x2 ',XYNEL(0,1,NEL),XYNEL(NPLOK1,1,NEL)
	  WRITE(99,*) ' y1,y2 ',XYNEL(0,3,NEL),XYNEL(NPLOK3,3,NEL)
	  WRITE(99,*) ' "+1" "-2" : ', SNGL(+ARGNEL(NPLOK1,1,NEL)-
     $	  ARGNEL(0,1,NEL)) , SNGL(-ARGNEL(NPLOK2,2,NEL)+ARGNEL(0,2,NEL))
	  WRITE(99,*) ' "-3" "+4" : ', SNGL(-ARGNEL(NPLOK3,3,NEL)+
     $	  ARGNEL(0,3,NEL)) , SNGL(+ARGNEL(NPLOK4,4,NEL)-ARGNEL(0,4,NEL))
	END IF
	IF (NROOT.LT.0) STOP 125
	IF (IOPTION.EQ.0) THEN
	  NZLEFT = NZLEFT + NROOT
	  WRITE(LUNLEFT,*)
	  WRITE(LUNLEFT,*) NROOT,'  (ioption=0)'
	  WRITE(LUNLEFT,*) SNGL(X1),SNGL(X2),SNGL(Y1),SNGL(Y2)
	  RETURN
	END IF

	IF (NROOT.EQ.0) THEN

	  NEL = NEL - 1
	  IF (NEL.EQ.0) RETURN
	  GOTO 100

	ELSE IF (NROOT.EQ.1) THEN

	  IF (.NOT.ONEROOT) THEN
	    X1ONE = X1
	    X2ONE = X2
	    Y1ONE = Y1
	    Y2ONE = Y2
	    ONEROOT = .TRUE.
	  END IF

	  IF (MCHANGE.EQ.0) THEN
C	    simple midpoint estimate of "initial (px,py)"
	    PX = 0.5D0 * (X1+X2)
	    PY = 0.5D0 * (Y1+Y2)
	  ELSE
C	    "initial (px,py)" by "integrating z*d(rlndet+ln(abs(cdet)))"
C	    (translation for numerical reasons however)
	    XLOK2 = X1 + X2
	    YLOK2 = Y1 + Y2
	    XLOK = 0.5D0 * XLOK2
	    YLOK = 0.5D0 * YLOK2
	    SLASK1 = 0.0D0
	    SLASK2 = 0.0D0
	    DO I = 1, NPLOK1
	      SLASK = XYNEL(I-1,1,NEL) + XYNEL(I,1,NEL) - XLOK2
	      SLASK1 = SLASK1 + SLASK*(ALNNEL(I,1,NEL)-ALNNEL(I-1,1,NEL))
	      SLASK2 = SLASK2 + SLASK*(ARGNEL(I,1,NEL)-ARGNEL(I-1,1,NEL))
	    END DO
	    DO I = 1, NPLOK2
	      SLASK = XYNEL(I-1,2,NEL) + XYNEL(I,2,NEL) - XLOK2
	      SLASK1 = SLASK1 + SLASK*(ALNNEL(I-1,2,NEL)-ALNNEL(I,2,NEL))
	      SLASK2 = SLASK2 + SLASK*(ARGNEL(I-1,2,NEL)-ARGNEL(I,2,NEL))
	    END DO
	    PX = SLASK2 / 720
	    PY = -SLASK1 / (4*PI)
	    PX = PX + ( (Y1-YLOK)*(ALNNEL(NPLOK1,1,NEL)-
     $	    ALNNEL(0,1,NEL)) + (Y2-YLOK)*(ALNNEL(0,2,NEL)-
     $	    ALNNEL(NPLOK2,2,NEL)) ) / (2*PI)
	    PY = PY + ( (Y1-YLOK)*(ARGNEL(NPLOK1,1,NEL)-
     $	    ARGNEL(0,1,NEL)) + (Y2-YLOK)*(ARGNEL(0,2,NEL)-
     $	    ARGNEL(NPLOK2,2,NEL)) ) / 360
	    SLASK1 = 0.0D0
	    SLASK2 = 0.0D0
	    DO I = 1, NPLOK4
	      SLASK = XYNEL(I-1,4,NEL) + XYNEL(I,4,NEL) - YLOK2
	      SLASK1 = SLASK1 + SLASK*(ALNNEL(I,4,NEL)-ALNNEL(I-1,4,NEL))
	      SLASK2 = SLASK2 + SLASK*(ARGNEL(I,4,NEL)-ARGNEL(I-1,4,NEL))
	    END DO
	    DO I = 1, NPLOK3
	      SLASK = XYNEL(I-1,3,NEL) + XYNEL(I,3,NEL) - YLOK2
	      SLASK1 = SLASK1 + SLASK*(ALNNEL(I-1,3,NEL)-ALNNEL(I,3,NEL))
	      SLASK2 = SLASK2 + SLASK*(ARGNEL(I-1,3,NEL)-ARGNEL(I,3,NEL))
	    END DO
	    PX = PX + SLASK1 / (4*PI)
	    PY = PY + SLASK2 / 720
	    PX = PX + ( (X2-XLOK)*(ARGNEL(NPLOK4,4,NEL)-
     $	    ARGNEL(0,4,NEL)) + (X1-XLOK)*(ARGNEL(0,3,NEL)-
     $	    ARGNEL(NPLOK3,3,NEL)) ) / 360
	    PY = PY - ( (X2-XLOK)*(ALNNEL(NPLOK4,4,NEL)-
     $	    ALNNEL(0,4,NEL)) + (X1-XLOK)*(ALNNEL(0,3,NEL)-
     $	    ALNNEL(NPLOK3,3,NEL)) ) / (2*PI)
	    PX = XLOK + PX
	    PY = YLOK + PY
ccc!!!	    write(111,*)
ccc!!!	    write(111,*) x1,x2,' ! x1,x2'
ccc!!!	    write(111,*) y1,y2,' ! y1,y2'
ccc!!!	    write(111,*) px,py,' ! px,py pre'
	    IF (PX.LE.X1) THEN
	      PX = 0.9D0*X1 + 0.1D0*X2
	    ELSE IF (PX.GE.X2) THEN
	      PX = 0.1D0*X1 + 0.9D0*X2
	    END IF
	    IF (PY.LE.Y1) THEN
	      PY = 0.9D0*Y1 + 0.1D0*Y2
	    ELSE IF (PY.GE.Y2) THEN
	      PY = 0.1D0*Y1 + 0.9D0*Y2
	    END IF
	  END IF
C	  evaluate the analytic function at (px,py)
	  IDZHFCALL = IDZHFCALL + 1
	  CALL CDETSUB (PX,PY, RLNDET,CDET)

	  CALL DZEROHOLO (CDETSUB, X1,X2,Y1,Y2, DLOGREDUCFAKT,DDELT,
     $	  XYTOL, PX,PY,XYTOLUT,RLNDET,CDET, IDZHFCALL,IDZHFAIL)
	  IF (IDZHFAIL.EQ.1) THEN
	    IDZHBACK = IDZHBACK + 1
ccc!!!	    write(111,*) -999,-999,' ! back to 110'
	    GOTO 110
	  ELSE IF (IDZHFAIL.NE.0) THEN
	    STOP 202
	  END IF
ccc!!!	  write(111,*) px,py,' ! px,py post'

	  IF (ICHECK.EQ.1) THEN
	    IF (XYTOLUT.GT.0.0D0) THEN
	      ARGSLASK = 0.0D0
	      UP0CIRC = DCMPLX ( PX, PY )          ! to "common /circdata/"
	      RADIECIRC = XYTOLUT                  ! to "common /circdata/"
	      UP1 = DCMPLX ( 0.0D0, 0.0D0 )
	      ICHKFCALL = ICHKFCALL + 1
	      CALL CDETSUBCIRC (0.0D0, 0.0D0, RLNDET0,CDET0)
	      UP2 = DCMPLX ( 2*PI, 0.0D0 )
	      DDSUP1 = 0.0D0
	      CALL ARGSEG (CDETSUBCIRC, UP1,UP2, RLNDET0,CDET0, DDSUP1, 
     $	      DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	      DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, 
     $	      NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET,CDET, 
     $	      ICHKFCALL,IASFAIL)
	      IF (NPLIN.GT.NPLINMAX) NPLINMAX = NPLIN
	      IF (IASFAIL.EQ.4) THEN
	        STOP 505
	      ELSE IF (IASFAIL.NE.0) THEN
	        ICHKBACK = ICHKBACK + 1
	        GOTO 110
	      END IF
	      ARGSLASK = ARGLIN(NPLIN)-ARGLIN(0)
	      IF ( IDNINT(ARGSLASK/360.0D0) .NE. 1 ) THEN
	        ICHKBACK = ICHKBACK + 1
	        GOTO 110
	      END IF
	    END IF
	  END IF

	  NZERO = NZERO + 1
	  IF (NZERO.GT.nzerodim) STOP 166
	  ZREAL(NZERO) = PX
	  ZIMAG(NZERO) = PY
	  XYONEMIN(NZERO) = XYTOLUT
	  XYONEMAX(NZERO) = DMIN1 ( PX-X1ONE, X2ONE-PX, PY-Y1ONE, Y2ONE-PY )
	  NEL = NEL - 1
	  IF (NEL.EQ.0) RETURN
	  ONEROOT = .FALSE.
	  GOTO 100

	END IF

 110	IF ( (X2-X1) .GE. (Y2-Y1) ) THEN

	  NEL1 = NEL + 1                                          ! x division
	  IF (NEL1.GT.neldim) THEN
C	    actually the following "nel compression" could be done implicitly 
C	    much earlier by avoiding to "store with a nel" when "nroot=0"
	    INELCOMPR = INELCOMPR + 1
	    CALL NELCOMPR (nplindim, NEL, NPNEL, XYNEL,ARGNEL,ALNNEL,DDSUPNEL)
	    NEL1 = NEL + 1
	    IF (NEL1.GT.neldim) STOP 715
	  END IF
	  XLOK = 0.5D0 * (X1+X2)
 10	  IF (XLOK.EQ.0.0D0) THEN
	    XLOK = 0.5D0 * (X1+XLOK)      ! avoid the axes
	  END IF
	  IF ( XLOK.LE.(X1+SLIMLEFT) .OR. XLOK.GE.(X2-SLIMLEFT) ) THEN 
	    NZLEFT = NZLEFT + NROOT
	    WRITE(LUNLEFT,*)
	    WRITE(LUNLEFT,*) NROOT,'   (slimleft active here)'
	    WRITE(LUNLEFT,*) SNGL(X1),SNGL(X2),SNGL(Y1),SNGL(Y2)
	    NEL = NEL - 1
	    IF (NEL.EQ.0) RETURN
	    ONEROOT = .FALSE.
	    GOTO 100
	  END IF

	  INDEX = NPLOK1/2
	  DO WHILE (XYNEL(INDEX,1,NEL).LT.XLOK)
	    INDEX = INDEX + 1
	  END DO
	  DO WHILE (XYNEL(INDEX,1,NEL).GT.XLOK)
	    INDEX = INDEX - 1
	  END DO

	  UP1 = DCMPLX (XLOK,Y1)
	  IASFCALL = IASFCALL + 1
	  CALL CDETSUB (XLOK,Y1, RLNDET,CDET)
	  UP2 = DCMPLX (XLOK,Y2)
	  DDSUP1 = DDSUPNEL(INDEX+1,1,NEL)
	  CALL ARGSEG (CDETSUB, UP1,UP2, RLNDET,CDET, DDSUP1, 
     $	  DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	  DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, 
     $	  NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET,CDET, 
     $	  IASFCALL,IASFAIL)
	  IF (NPLIN.GT.NPLINMAX) NPLINMAX = NPLIN
	  IF (IASFAIL.NE.0) THEN
	    IASBACK = IASBACK + 1
	    XLOK = 0.5D0 * (X1+XLOK)
	    GOTO 10
	  END IF
	  NPNEL(4,NEL1) = NPLOK4
	  DO I = 0, NPLOK4
	    XYNEL(I,4,NEL1) = XYNEL(I,4,NEL)
	  END DO
	  DO I = 0, NPLOK4
	    ARGNEL(I,4,NEL1) = ARGNEL(I,4,NEL)
	  END DO
	  DO I = 0, NPLOK4
	    ALNNEL(I,4,NEL1) = ALNNEL(I,4,NEL)
	  END DO
	  DO I = 1, NPLOK4
	    DDSUPNEL(I,4,NEL1) = DDSUPNEL(I,4,NEL)
	  END DO
	  NPNEL(3,NEL1) = NPLIN
	  NPNEL(4,NEL) = NPLIN
	  DO I = 0, NPLIN
	    XYNEL(I,3,NEL1) = DIMAG(UPLIN(I))
	  END DO
	  DO I = 0, NPLIN
	    XYNEL(I,4,NEL) = DIMAG(UPLIN(I))
	  END DO
	  DO I = 0, NPLIN
	    ARGNEL(I,3,NEL1) = ARGLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,3,NEL1) = ALNLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,3,NEL1) = ALNLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ARGNEL(I,4,NEL) = ARGLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,4,NEL) = ALNLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,4,NEL) = ALNLIN(I)
	  END DO
	  DO I = 1, NPLIN
	    DDSUPNEL(I,3,NEL1) = DDSUPLIN(I)
	  END DO
	  DO I = 1, NPLIN
	    DDSUPNEL(I,4,NEL) = DDSUPLIN(I)
	  END DO

	  NPNEL(1,NEL1) = NPLOK1 - INDEX
	  NPNEL(1,NEL) = INDEX + 1
	  XYNEL(0,1,NEL1) = XLOK
	  DO I = 1, NPNEL(1,NEL1)
	    XYNEL(I,1,NEL1) = XYNEL(INDEX+I,1,NEL)
	  END DO
	  XYNEL(NPNEL(1,NEL),1,NEL) = XLOK
	  SLASK = 0.5D0 * ( ARGNEL(INDEX,1,NEL) + ARGNEL(INDEX+1,1,NEL) )
	  ARGSLASK = ARGLIN(0) + 360.0D0*IDNINT((SLASK-ARGLIN(0))/360.0D0)
	  IF (DABS(ARGSLASK-SLASK).GE.ARGDIFTEST) THEN
	    WRITE(*,*)
	    WRITE(*,*) NROOT,'  %%% argdiftest failure %%%'
	    WRITE(*,*) SNGL(X1),SNGL(X2),SNGL(Y1),SNGL(Y2)
	    WRITE(*,*) SNGL(ARGSLASK),SNGL(SLASK),
     $	    SNGL(ARGSLASK-SLASK),SNGL(ARGDIFTEST)
	    STOP 103
	  END IF
	  ARGNEL(0,1,NEL1) = ARGSLASK
	  ALNNEL(0,1,NEL1) = ALNLIN(0)
	  DO I = 1, NPNEL(1,NEL1)
	    ARGNEL(I,1,NEL1) = ARGNEL(INDEX+I,1,NEL)
	  END DO
	  DO I = 1, NPNEL(1,NEL1)
	    ALNNEL(I,1,NEL1) = ALNNEL(INDEX+I,1,NEL)
	  END DO
	  ARGNEL(NPNEL(1,NEL),1,NEL) = ARGSLASK
	  ALNNEL(NPNEL(1,NEL),1,NEL) = ALNLIN(0)
	  DO I = 1, NPNEL(1,NEL1)
	    DDSUPNEL(I,1,NEL1) = DDSUPNEL(INDEX+I,1,NEL)
	  END DO

	  INDEX = NPLOK2/2
	  DO WHILE (XYNEL(INDEX,2,NEL).LT.XLOK)
	    INDEX = INDEX + 1
	  END DO
	  DO WHILE (XYNEL(INDEX,2,NEL).GT.XLOK)
	    INDEX = INDEX - 1
	  END DO
	  NPNEL(2,NEL1) = NPLOK2 - INDEX
	  NPNEL(2,NEL) = INDEX + 1
	  XYNEL(0,2,NEL1) = XLOK
	  DO I = 1, NPNEL(2,NEL1)
	    XYNEL(I,2,NEL1) = XYNEL(INDEX+I,2,NEL)
	  END DO
	  XYNEL(NPNEL(2,NEL),2,NEL) = XLOK
	  SLASK = 0.5D0 * ( ARGNEL(INDEX,2,NEL) + ARGNEL(INDEX+1,2,NEL) )
	  ARGSLASK = ARGLIN(NPLIN) + 360.0D0*
     $	  IDNINT((SLASK-ARGLIN(NPLIN))/360.0D0)
	  IF (DABS(ARGSLASK-SLASK).GE.ARGDIFTEST) THEN
	    WRITE(*,*)
	    WRITE(*,*) NROOT,'  %%% argdiftest failure %%%'
	    WRITE(*,*) SNGL(X1),SNGL(X2),SNGL(Y1),SNGL(Y2)
	    WRITE(*,*) SNGL(ARGSLASK),SNGL(SLASK),
     $	    SNGL(ARGSLASK-SLASK),SNGL(ARGDIFTEST)
	    STOP 103
	  END IF
	  ARGNEL(0,2,NEL1) = ARGSLASK
	  ALNNEL(0,2,NEL1) = ALNLIN(NPLIN)
	  DO I = 1, NPNEL(2,NEL1)
	    ARGNEL(I,2,NEL1) = ARGNEL(INDEX+I,2,NEL)
	  END DO
	  DO I = 1, NPNEL(2,NEL1)
	    ALNNEL(I,2,NEL1) = ALNNEL(INDEX+I,2,NEL)
	  END DO
	  ARGNEL(NPNEL(2,NEL),2,NEL) = ARGSLASK
	  ALNNEL(NPNEL(2,NEL),2,NEL) = ALNLIN(NPLIN)
	  DO I = 1, NPNEL(2,NEL1)
	    DDSUPNEL(I,2,NEL1) = DDSUPNEL(INDEX+I,2,NEL)
	  END DO

	ELSE

	  NEL1 = NEL + 1                                          ! y division
	  IF (NEL1.GT.neldim) THEN
C	    actually the following "nel compression" could be done implicitly 
C	    much earlier by avoiding to "store with a nel" when "nroot=0"
	    INELCOMPR = INELCOMPR + 1
	    CALL NELCOMPR (nplindim, NEL, NPNEL, XYNEL,ARGNEL,ALNNEL,DDSUPNEL)
	    NEL1 = NEL + 1
	    IF (NEL1.GT.neldim) STOP 715
	  END IF
	  YLOK = 0.5D0 * (Y1+Y2)
 20	  IF (YLOK.EQ.0.0D0) THEN
	    YLOK = 0.5D0 * (Y1+YLOK)      ! avoid the axes
	  END IF
	  IF ( YLOK.LE.(Y1+SLIMLEFT) .OR. YLOK.GE.(Y2-SLIMLEFT) ) THEN 
	    NZLEFT = NZLEFT + NROOT
	    WRITE(LUNLEFT,*)
	    WRITE(LUNLEFT,*) NROOT,'   (slimleft active here)'
	    WRITE(LUNLEFT,*) SNGL(X1),SNGL(X2),SNGL(Y1),SNGL(Y2)
	    NEL = NEL - 1
	    IF (NEL.EQ.0) RETURN
	    ONEROOT = .FALSE.
	    GOTO 100
	  END IF

	  INDEX = NPLOK3/2
	  DO WHILE (XYNEL(INDEX,3,NEL).LT.YLOK)
	    INDEX = INDEX + 1
	  END DO
	  DO WHILE (XYNEL(INDEX,3,NEL).GT.YLOK)
	    INDEX = INDEX - 1
	  END DO

	  UP1 = DCMPLX (X1,YLOK)
	  IASFCALL = IASFCALL + 1
	  CALL CDETSUB (X1,YLOK, RLNDET,CDET)
	  UP2 = DCMPLX (X2,YLOK)
	  DDSUP1 = DDSUPNEL(INDEX+1,3,NEL)
	  CALL ARGSEG (CDETSUB, UP1,UP2, RLNDET,CDET, DDSUP1, 
     $	  DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	  DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, 
     $	  NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET,CDET, 
     $	  IASFCALL,IASFAIL)
	  IF (NPLIN.GT.NPLINMAX) NPLINMAX = NPLIN
	  IF (IASFAIL.NE.0) THEN
	    IASBACK = IASBACK + 1
	    YLOK = 0.5D0 * (Y1+YLOK)
	    GOTO 20
	  END IF
	  NPNEL(2,NEL1) = NPLOK2
	  DO I = 0, NPLOK2
	    XYNEL(I,2,NEL1) = XYNEL(I,2,NEL)
	  END DO
	  DO I = 0, NPLOK2
	    ARGNEL(I,2,NEL1) = ARGNEL(I,2,NEL)
	  END DO
	  DO I = 0, NPLOK2
	    ALNNEL(I,2,NEL1) = ALNNEL(I,2,NEL)
	  END DO
	  DO I = 1, NPLOK2
	    DDSUPNEL(I,2,NEL1) = DDSUPNEL(I,2,NEL)
	  END DO
	  NPNEL(1,NEL1) = NPLIN
	  NPNEL(2,NEL) = NPLIN
	  DO I = 0, NPLIN
	    XYNEL(I,1,NEL1) = DREAL(UPLIN(I))
	  END DO
	  DO I = 0, NPLIN
	    XYNEL(I,2,NEL) = DREAL(UPLIN(I))
	  END DO
	  DO I = 0, NPLIN
	    ARGNEL(I,1,NEL1) = ARGLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,1,NEL1) = ALNLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,1,NEL1) = ALNLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ARGNEL(I,2,NEL) = ARGLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,2,NEL) = ALNLIN(I)
	  END DO
	  DO I = 0, NPLIN
	    ALNNEL(I,2,NEL) = ALNLIN(I)
	  END DO
	  DO I = 1, NPLIN
	    DDSUPNEL(I,1,NEL1) = DDSUPLIN(I)
	  END DO
	  DO I = 1, NPLIN
	    DDSUPNEL(I,2,NEL) = DDSUPLIN(I)
	  END DO

	  NPNEL(3,NEL1) = NPLOK3 - INDEX
	  NPNEL(3,NEL) = INDEX + 1
	  XYNEL(0,3,NEL1) = YLOK
	  DO I = 1, NPNEL(3,NEL1)
	    XYNEL(I,3,NEL1) = XYNEL(INDEX+I,3,NEL)
	  END DO
	  XYNEL(NPNEL(3,NEL),3,NEL) = YLOK
	  SLASK = 0.5D0 * ( ARGNEL(INDEX,3,NEL) + ARGNEL(INDEX+1,3,NEL) )
	  ARGSLASK = ARGLIN(0) + 360.0D0*IDNINT((SLASK-ARGLIN(0))/360.0D0)
	  IF (DABS(ARGSLASK-SLASK).GE.ARGDIFTEST) THEN
	    WRITE(*,*)
	    WRITE(*,*) NROOT,'  %%% argdiftest failure %%%'
	    WRITE(*,*) SNGL(X1),SNGL(X2),SNGL(Y1),SNGL(Y2)
	    WRITE(*,*) SNGL(ARGSLASK),SNGL(SLASK),
     $	    SNGL(ARGSLASK-SLASK),SNGL(ARGDIFTEST)
	    STOP 103
	  END IF
	  ARGNEL(0,3,NEL1) = ARGSLASK
	  ALNNEL(0,3,NEL1) = ALNLIN(0)
	  DO I = 1, NPNEL(3,NEL1)
	    ARGNEL(I,3,NEL1) = ARGNEL(INDEX+I,3,NEL)
	  END DO
	  DO I = 1, NPNEL(3,NEL1)
	    ALNNEL(I,3,NEL1) = ALNNEL(INDEX+I,3,NEL)
	  END DO
	  ARGNEL(NPNEL(3,NEL),3,NEL) = ARGSLASK
	  ALNNEL(NPNEL(3,NEL),3,NEL) = ALNLIN(0)
	  DO I = 1, NPNEL(3,NEL1)
	    DDSUPNEL(I,3,NEL1) = DDSUPNEL(INDEX+I,3,NEL)
	  END DO

	  INDEX = NPLOK4/2
	  DO WHILE (XYNEL(INDEX,4,NEL).LT.YLOK)
	    INDEX = INDEX + 1
	  END DO
	  DO WHILE (XYNEL(INDEX,4,NEL).GT.YLOK)
	    INDEX = INDEX - 1
	  END DO
	  NPNEL(4,NEL1) = NPLOK4 - INDEX
	  NPNEL(4,NEL) = INDEX + 1
	  XYNEL(0,4,NEL1) = YLOK
	  DO I = 1, NPNEL(4,NEL1)
	    XYNEL(I,4,NEL1) = XYNEL(INDEX+I,4,NEL)
	  END DO
	  XYNEL(NPNEL(4,NEL),4,NEL) = YLOK
	  SLASK = 0.5D0 * ( ARGNEL(INDEX,4,NEL) + ARGNEL(INDEX+1,4,NEL) )
	  ARGSLASK = ARGLIN(NPLIN) + 360.0D0*
     $	  IDNINT((SLASK-ARGLIN(NPLIN))/360.0D0)
	  IF (DABS(ARGSLASK-SLASK).GE.ARGDIFTEST) THEN
	    WRITE(*,*)
	    WRITE(*,*) NROOT,'  %%% argdiftest failure %%%'
	    WRITE(*,*) SNGL(X1),SNGL(X2),SNGL(Y1),SNGL(Y2)
	    WRITE(*,*) SNGL(ARGSLASK),SNGL(SLASK),
     $	    SNGL(ARGSLASK-SLASK),SNGL(ARGDIFTEST)
	    STOP 103
	  END IF
	  ARGNEL(0,4,NEL1) = ARGSLASK
	  ALNNEL(0,4,NEL1) = ALNLIN(NPLIN)
	  DO I = 1, NPNEL(4,NEL1)
	    ARGNEL(I,4,NEL1) = ARGNEL(INDEX+I,4,NEL)
	  END DO
	  DO I = 1, NPNEL(4,NEL1)
	    ALNNEL(I,4,NEL1) = ALNNEL(INDEX+I,4,NEL)
	  END DO
	  ARGNEL(NPNEL(4,NEL),4,NEL) = ARGSLASK
	  ALNNEL(NPNEL(4,NEL),4,NEL) = ALNLIN(NPLIN)
	  DO I = 1, NPNEL(4,NEL1)
	    DDSUPNEL(I,4,NEL1) = DDSUPNEL(INDEX+I,4,NEL)
	  END DO

	END IF

	NEL = NEL1
	IF (NEL.GT.NELMAX) NELMAX = NEL
	GOTO 100

C	---------------------------------------------------

	END








	SUBROUTINE NELCOMPR (NPLINDIM, NEL, NPNEL, XYNEL,ARGNEL,
     $	ALNNEL,DDSUPNEL)

	IMPLICIT NONE

	INTEGER NPLINDIM

	INTEGER NEL, NPNEL(4,*)
	REAL ARGNEL(0:NPLINDIM,4,*), ALNNEL(0:NPLINDIM,4,*), 
     $	DDSUPNEL(NPLINDIM,4,*)          !!!
	REAL*8 XYNEL(0:NPLINDIM,4,*) 

	INTEGER NELFREE,NELTEST, NPLOK1,NPLOK2,NPLOK3,NPLOK4, NROOT, I
	REAL*8 SLASK

	NELFREE = 0

	DO NELTEST = 1, NEL

	  NPLOK1 = NPNEL(1,NELTEST)
	  NPLOK2 = NPNEL(2,NELTEST)
	  NPLOK3 = NPNEL(3,NELTEST)
	  NPLOK4 = NPNEL(4,NELTEST)
	  SLASK = ( + ARGNEL(NPLOK1,1,NELTEST) - ARGNEL(0,1,NELTEST) 
     $	            - ARGNEL(NPLOK2,2,NELTEST) + ARGNEL(0,2,NELTEST) 
     $	            - ARGNEL(NPLOK3,3,NELTEST) + ARGNEL(0,3,NELTEST) 
     $	            + ARGNEL(NPLOK4,4,NELTEST) - ARGNEL(0,4,NELTEST) ) / 
     $	          360.0D0
	  NROOT = IDNINT (SLASK)

	  IF (NROOT.NE.0) THEN

	    IF (NELFREE.GT.0) THEN

	      NPNEL(1,NELFREE) = NPNEL(1,NELTEST)
	      NPNEL(2,NELFREE) = NPNEL(2,NELTEST)
	      NPNEL(3,NELFREE) = NPNEL(3,NELTEST)
	      NPNEL(4,NELFREE) = NPNEL(4,NELTEST)

	      DO I = 0, NPLOK1
	        XYNEL(I,1,NELFREE) = XYNEL(I,1,NELTEST)
	      END DO
	      DO I = 0, NPLOK2
	        XYNEL(I,2,NELFREE) = XYNEL(I,2,NELTEST)
	      END DO
	      DO I = 0, NPLOK3
	        XYNEL(I,3,NELFREE) = XYNEL(I,3,NELTEST)
	      END DO
	      DO I = 0, NPLOK4
	        XYNEL(I,4,NELFREE) = XYNEL(I,4,NELTEST)
	      END DO

	      DO I = 0, NPLOK1
	        ARGNEL(I,1,NELFREE) = ARGNEL(I,1,NELTEST)
	      END DO
	      DO I = 0, NPLOK1
	        ALNNEL(I,1,NELFREE) = ALNNEL(I,1,NELTEST)
	      END DO
	      DO I = 0, NPLOK2
	        ARGNEL(I,2,NELFREE) = ARGNEL(I,2,NELTEST)
	      END DO
	      DO I = 0, NPLOK2
	        ALNNEL(I,2,NELFREE) = ALNNEL(I,2,NELTEST)
	      END DO
	      DO I = 0, NPLOK3
	        ALNNEL(I,3,NELFREE) = ALNNEL(I,3,NELTEST)
	      END DO
	      DO I = 0, NPLOK4
	        ARGNEL(I,4,NELFREE) = ARGNEL(I,4,NELTEST)
	      END DO
	      DO I = 0, NPLOK4
	        ALNNEL(I,4,NELFREE) = ALNNEL(I,4,NELTEST)
	      END DO

	      DO I = 1, NPLOK1
	        DDSUPNEL(I,1,NELFREE) = DDSUPNEL(I,1,NELTEST)
	      END DO
	      DO I = 1, NPLOK2
	        DDSUPNEL(I,2,NELFREE) = DDSUPNEL(I,2,NELTEST)
	      END DO
	      DO I = 1, NPLOK3
	        DDSUPNEL(I,3,NELFREE) = DDSUPNEL(I,3,NELTEST)
	      END DO
	      DO I = 1, NPLOK4
	        DDSUPNEL(I,4,NELFREE) = DDSUPNEL(I,4,NELTEST)
	      END DO

	      NELFREE = NELFREE + 1

	    END IF

	  ELSE

	    IF (NELFREE.EQ.0) NELFREE = NELTEST

	  END IF

	END DO

	IF (NELFREE.NE.0) NEL = NELFREE - 1

	RETURN
	END