Here is an attempt. I made minimal modernization, just getting rid of the labels and GOTO’s and converting to free-form. There would be a bit more to do…
Actually, most of the goto’s were simulations of IF/THEN/ELSE blocks. Some other ones were more tricky.
The code structure remains non conventional because of the reverse communication constraints. And IMO the main issue with the routine is not the GOTO’s, but the SAVE statement that makes all the local variables static, in order to keep the values between calls. Hence the routine is not thread-safe. At the time the code was written the only alternative would have been to pass all these variables in the argumens. With current Fortran the right way would be to put them in a type
and pass a single argument, and enclose everything in a module
.
I have quickly tested it, but more tests would be needed.
SUBROUTINE dzror(status,x,fx,xlo,xhi,qleft,qhi)
!**********************************************************************
IMPLICIT NONE
! Scalar Arguments
DOUBLE PRECISION fx,x,xhi,xlo,zabstl,zreltl,zxhi,zxlo
INTEGER status
LOGICAL qhi,qleft
! Save statement !!! makes the routine non thread safe !!!
SAVE
! Local Scalars
DOUBLE PRECISION a,abstol,b,c,d,fa,fb,fc,fd,fda,fdb,m,mb,p,q, &
reltol,tol,w,xxhi,xxlo,zx
INTEGER ext
CHARACTER(16) phase
LOGICAL first,qrzero
! Intrinsic Functions
INTRINSIC abs,max,sign
! Statement Functions
DOUBLE PRECISION ftol
! Statement Function definitions
ftol(zx) = 0.5D0*max(abstol,reltol*abs(zx))
! Executable Statements ..
! Fake infinite loop to simulate jumps to specific phases
! phase = 'XXX' ; RETURN
! exits the routine to collect a new function value, and come back at phase "XXX"
! PHASE = 'XXX' [; CYCLE]
! equivalent to jump to the phase "XXX"
! CYCLE can be omitted if "XXX" is located after in the code and if this is the
! last instruction of the current phase
CONTROL_LOOP: DO
IF (phase == 'START') THEN
xlo = xxlo
xhi = xxhi
b = xlo
x = xlo
! GET-FUNCTION-VALUE
status = 1
phase = 'PHASE1'; RETURN
END IF
IF (phase == 'PHASE1') THEN
fb = fx
xlo = xhi
a = xlo
x = xlo
! GET-FUNCTION-VALUE
phase = 'PHASE2'; RETURN
END IF
IF (phase == 'PHASE2') THEN
! Check that F(ZXLO) < 0 < F(ZXHI) or
! F(ZXLO) > 0 > F(ZXHI)
IF (fb.LT.0.0D0 .AND. fx.LT.0.0D0) THEN
qleft = fx .LT. fb
qhi = .FALSE.
status = -1
RETURN ! termination
END IF
IF (fb.GT.0.0D0 .AND. fx.GT.0.0D0) THEN
qleft = fx .GT. fb
qhi = .TRUE.
status = -1
RETURN ! termination
END IF
fa = fx
first = .TRUE.
phase = 'PHASE2_B'
END IF
IF (phase == 'PHASE2_B') THEN
c = a
fc = fa
ext = 0
phase = 'PHASE2_C'
END IF
IF (phase == 'PHASE2_C') THEN
IF (abs(fc).LT.abs(fb)) THEN
IF (abs(fc).LT.abs(fb)) THEN
IF (c.NE.a) THEN
d = a
fd = fa
END IF
a = b
fa = fb
xlo = c
b = xlo
fb = fc
c = a
fc = fa
END IF
END IF
tol = ftol(xlo)
m = (c+b)*.5D0
mb = m - b
IF (.NOT. (abs(mb).GT.tol)) THEN
phase = 'PHASE3_B'
ELSE
IF (ext.GT.3) THEN
w = mb
ELSE
tol = sign(tol,mb)
p = (b-a)*fb
IF (first) THEN
q = fa - fb
first = .FALSE.
ELSE
fdb = (fd-fb)/ (d-b)
fda = (fd-fa)/ (d-a)
p = fda*p
q = fdb*fa - fda*fb
END IF
IF (p.LT.0.0D0) THEN
p = -p
q = -q
END IF
IF (ext.EQ.3) p = p*2.0D0
IF ((p*1.0D0).EQ.0.0D0 .OR. p.LE.(q*tol)) THEN
w = tol
ELSE
IF (p.LT. (mb*q)) THEN
w = p/q
ELSE
w = mb
END IF
CONTINUE
END IF
CONTINUE
END IF
d = a
fd = fa
a = b
fa = fb
b = b + w
xlo = b
x = xlo
! GET-FUNCTION-VALUE
phase = 'PHASE3' ; RETURN
END IF
END IF
IF (phase == 'PHASE3') THEN
fb = fx
IF ((fc*fb).GE.0.0D0) THEN
phase = 'PHASE2_B' ; CYCLE CONTROL_LOOP
END IF
IF (w.EQ.mb) THEN
ext = 0
ELSE
ext = ext + 1
END IF
phase = 'PHASE2_C' ; CYCLE CONTROL_LOOP
END IF
IF (phase == 'PHASE3_B') THEN
xhi = c
qrzero = (fc.GE.0.0D0 .AND. fb.LE.0.0D0) .OR. &
(fc.LT.0.0D0 .AND. fb.GE.0.0D0)
IF (qrzero) THEN
status = 0
ELSE
status = -1
END IF
RETURN ! termination
END IF
EXIT CONTROL_LOOP
END DO CONTROL_LOOP
ENTRY dstzr(zxlo,zxhi,zabstl,zreltl)
!**********************************************************************
! Reset all saved variables to known values
a = 0d0
abstol = 0d0
b = 0d0
c = 0d0
d = 0d0
fa = 0d0
fb = 0d0
fc = 0d0
fd = 0d0
fda = 0d0
fdb = 0d0
m = 0d0
mb = 0d0
p = 0d0
q = 0d0
reltol = 0d0
tol = 0d0
w = 0d0
xxhi = 0d0
xxlo = 0d0
zx = 0d0
ext = 0
phase = 'START'
first = .FALSE.
qrzero = .FALSE.
! Set initial values
xxlo = zxlo
xxhi = zxhi
abstol = zabstl
reltol = zreltl
RETURN
! The following stop statement can never be reached (??)
STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***'
END