C	subroutine "zerodivy" to a zero-free "upper part" of a given 
C	axis-parallel rectangle
C	adapted from zerofind.f

C	Sven Ivansson 1990






	SUBROUTINE ZERODIVY (CDETSUB, UPREAL1,UPREAL2,UPIMAG1,UPIMAG2, 
     $	DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,
     $	DRELGOALFAKT, ARGDIFSTORE, 
     $	NPLINMAX,IASFCALL,IASBACK, 
     $	YBELOW)
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	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  drelmax < 1.0d0 
C	                    ("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	nplinmax         input/output (OBS)  maximum current value of nplin
C	iasfcall         input/output (OBS)  total number of cdetsub calls 
C	                         from argseg
C	iasback          input/output (OBS)  total number of "back steps" 
C	                         upon unsuccessful return from argseg
C	ybelow           output  lower y limit of the zero-free rectangle

	IMPLICIT NONE

	INTEGER NPLINDIM
	PARAMETER ( NPLINDIM = 1000 )
C       nplin <= nplindim

	REAL*8 UPREAL1,UPREAL2,UPIMAG1,UPIMAG2, DDSUPMIN,DDSUPMAX,
     $	DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,DRELGOALFAKT, ARGDIFSTORE
	EXTERNAL CDETSUB

	INTEGER NPLINMAX,IASFCALL,IASBACK

	REAL*8 YBELOW

	INTEGER NPLIN, IASFAIL, NROOT, 
     $	NPLOK1,NPLOK2,NPLOK3,NPLOK4, INDEX, I, 
     $	NPNEL(3:4)
	REAL ARGNEL1D,ARGNEL2D, DDSUPNEL11, 
     $	ARGNEL(0:nplindim,3:4), DDSUPNEL(nplindim,3:3)     !!!
	REAL*8 DDSUP1,ARGDIFTEST, Y1, 
     $	YLOK, ARGSLASK, SLASK, RLNDET,RLNDET0, 
     $	ARGLIN(0:nplindim), ALNLIN(0:nplindim), DDSUPLIN(nplindim), 
     $	XYNEL(0:nplindim,3:4)
	COMPLEX*16 UP1,UP2, CDET,CDET0, 
     $	UPLIN(0:nplindim)

	ARGDIFTEST = 0.5D0*ARGDIFSTORE + DASIND(DRELMAX)

C	---------------------------------------------------
C	  a "nel" rectangle     !  attributes for the "successively stripped 
C	with its side indices   !       from below" "nel" rectangle
C	          2            /
C	   --------------     /             npnel(1:4) , numbers of points
C	   !            !    /                     (excluding the initial ones)
C	  3!            !4  !   xynel (0:npnel(,),1:4) , x/y coordinates
C	   !            !   !  argnel (0:npnel(,),1:4) , arg(cdet) degrees
C	   --------------  / ddsupnel (1:npnel(,),1:4) , "ddsupinf"
C	          1       /
C
C	note: storage in general not necessary for "sides 1 and 2" in this 
C	case, for ddsupnel not for "side 4" either
C
C	---------------------------------------------------

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

	UP1 = DCMPLX ( UPREAL1, UPIMAG1 )
	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
	NPLOK1 = NPLIN   !
	ARGNEL1D = ARGLIN(NPLIN) - ARGLIN(0)   !
	DDSUPNEL11 = DDSUPLIN(1)   !
ccc!	npnel(1) = nplin
ccc!	do i = 0, nplin
ccc!	  xynel(i,1) = dreal(uplin(i))
ccc!	end do
ccc!	do i = 0, nplin
ccc!	  argnel(i,1) = arglin(i)
ccc!	end do
ccc!	do i = 1, nplin
ccc!	  ddsupnel(i,1) = ddsuplin(i)
ccc!	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) = NPLIN
	DO I = 0, NPLIN
	  XYNEL(I,4) = DIMAG(UPLIN(I))
	END DO
	DO I = 0, NPLIN
	  ARGNEL(I,4) = ARGLIN(I)
	END DO
ccc!	do i = 1, nplin
ccc!	  ddsupnel(i,4) = ddsuplin(i)
ccc!	end do
	UP1 = DCMPLX ( UPREAL1, UPIMAG1 )
	UP2 = DCMPLX ( UPREAL1, UPIMAG2 )
	DDSUP1 = DDSUPNEL11
	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) = NPLIN
	DO I = 0, NPLIN
	  XYNEL(I,3) = DIMAG(UPLIN(I))
	END DO
	DO I = 0, NPLIN
	  ARGNEL(I,3) = ARGLIN(I)
	END DO
	DO I = 1, NPLIN
	  DDSUPNEL(I,3) = 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
	NPLOK2 = NPLIN   !
	ARGNEL2D = ARGLIN(NPLIN) - ARGLIN(0)   !
ccc!	npnel(2) = nplin
ccc!	do i = 0, nplin
ccc!	  xynel(i,2) = dreal(uplin(i))
ccc!	end do
ccc!	do i = 0, nplin
ccc!	  argnel(i,2) = arglin(i)
ccc!	end do
ccc!	do i = 1, nplin
ccc!	  ddsupnel(i,2) = ddsuplin(i)
ccc!	end do

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

 100	CONTINUE
ccc!	nplok1 = npnel(1)
ccc!	nplok2 = npnel(2)
	NPLOK3 = NPNEL(3)
	NPLOK4 = NPNEL(4)
	ARGSLASK = ( + ARGNEL1D 
     $	             - ARGNEL2D 
     $	             - ARGNEL(NPLOK3,3) + ARGNEL(0,3) 
     $	             + ARGNEL(NPLOK4,4) - ARGNEL(0,4) ) / 
     $	           360.0D0
	NROOT = IDNINT (ARGSLASK)
	IF (NROOT.EQ.0) THEN
	  YBELOW = XYNEL(0,3)
	  IF (YBELOW.GE.UPIMAG2) STOP 517
	  RETURN
	ELSE
ccc!	  x1 = xynel(0,1)
ccc!	  x2 = xynel(nplok1,1)
	  Y1 = XYNEL(0,3)
ccc!	  y2 = xynel(nplok3,3)
	END IF
	IF (NROOT.LT.0) STOP 125

 110	CONTINUE                                                ! y division

	YLOK = 0.5D0 * (Y1+UPIMAG2)
 20	IF (YLOK.EQ.0.0D0) THEN
	  YLOK = 0.5D0 * (Y1+YLOK)      ! avoid the axes
	END IF

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

	UP1 = DCMPLX (UPREAL1,YLOK)
	CALL CDETSUB (UPREAL1,YLOK, RLNDET,CDET)
	UP2 = DCMPLX (UPREAL2,YLOK)
	DDSUP1 = DDSUPNEL(INDEX+1,3)
	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
	NPLOK1 = NPLIN   !
	ARGNEL1D = ARGLIN(NPLIN) - ARGLIN(0)   !
ccc!	npnel(1) = nplin
ccc!	do i = 0, nplin
ccc!	  xynel(i,1) = dreal(uplin(i))
ccc!	end do
ccc!	do i = 0, nplin
ccc!	  argnel(i,1) = arglin(i)
ccc!	end do
ccc!	do i = 1, nplin
ccc!	  ddsupnel(i,1) = ddsuplin(i)
ccc!	end do

	NPNEL(3) = NPLOK3 - INDEX
	XYNEL(0,3) = YLOK
	DO I = 1, NPNEL(3)
	  XYNEL(I,3) = XYNEL(INDEX+I,3)
	END DO
	SLASK = 0.5D0 * ( ARGNEL(INDEX,3) + ARGNEL(INDEX+1,3) )
	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(UPREAL1),SNGL(UPREAL2),SNGL(Y1),SNGL(UPIMAG2)
	  STOP 103
	END IF
	ARGNEL(0,3) = ARGSLASK
	DO I = 1, NPNEL(3)
	  ARGNEL(I,3) = ARGNEL(INDEX+I,3)
	END DO
	DO I = 1, NPNEL(3)
	  DDSUPNEL(I,3) = DDSUPNEL(INDEX+I,3)
	END DO

	INDEX = NPLOK4/2
	DO WHILE (XYNEL(INDEX,4).LT.YLOK)
	  INDEX = INDEX + 1
	END DO
	DO WHILE (XYNEL(INDEX,4).GT.YLOK)
	  INDEX = INDEX - 1
	END DO
	NPNEL(4) = NPLOK4 - INDEX
	XYNEL(0,4) = YLOK
	DO I = 1, NPNEL(4)
	  XYNEL(I,4) = XYNEL(INDEX+I,4)
	END DO
	SLASK = 0.5D0 * ( ARGNEL(INDEX,4) + ARGNEL(INDEX+1,4) )
	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(UPREAL1),SNGL(UPREAL2),SNGL(Y1),SNGL(UPIMAG2)
	  STOP 103
	END IF
	ARGNEL(0,4) = ARGSLASK
	DO I = 1, NPNEL(4)
	  ARGNEL(I,4) = ARGNEL(INDEX+I,4)
	END DO
ccc!	do i = 1, npnel(4)
ccc!	  ddsupnel(i,4) = ddsupnel(index+i,4)
ccc!	end do

	GOTO 100

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

	END