二分法求弹塑性比 elastic plastic ratio using Bi-section

我的Fortran 源程序：
subroutine UWAPIPE_epp(p, force, forcen, dforcee, epp)
!
!! the elastic-plastic portion
use UWAPIPE_module
implicit none

type(propUWAPIPE), intent(in) :: p
real(kind=gp), dimension(3), intent(in) :: force, dforcee
real(kind=gp), dimension(2), intent(in) :: forcen
real(kind=gp), intent(out) :: epp

real(kind=gp), dimension(3) :: forceyd
real(kind=gp) :: a, b, f
integer(kind=8) :: n
a = 0.0_gp
b = 1.0_gp
n = 0
if (p%model == 'ep') then
do
epp = (a+b)/2.0_gp
forceyd = force + epp*dforcee
call UWAPIPE_yield(p, forceyd, p%V0, p%h0, f)
if (abs(f).LT.ftol) then
exit
else if (f.LT.-ftol) then
a = epp
else
b = epp
end if
n = n + 1
if (n.GE.1e6) then
print *, 'error. Can not find the epp. set to 0'
epp = 0.0_gp
exit
end if
end do
else if (p%model == 'bm') then
do
epp = (a+b)/2.0_gp
forceyd = force + epp*dforcee
call UWAPIPE_bubyield(p, forceyd, forcen, p%V0, p%h0, f)
if (abs(f).LT.ftol) then
exit
else if (f.LT.-ftol) then
a = epp
else
b = epp
end if
n = n + 1
if (n.GE.1e6) then
print *, 'error. Can not find the epp. set to 0'
epp = 0.0_gp
exit
end if
end do
end if

end subroutine UWAPIPE_epp