Fortran中利用梯形积分法怎么解一个内部有积分上限函数的积分

积分函数：

我的解答思路和尝试过的方法

function f3(r,r0,z0,z1,e,L,u,hh,V0,Q,a,r1)
implicit none
real(kind=8)::f3,Q,u,hh,fv,r,r0,z0,z1,Vc,Vn,e,V0,a,L,res,r1,PI,dh2
integer::j,m
real(kind=8),external::f2
m=1000.0
PI=3.1415926
if(r>=r0)then
    Vc=z0*z1*e**2/r
else if(r<r0)then
    Vc=z0*z1*e**2/2.0/r0*(3-(r/r0)**2)
endif
Vn=-V0*(1.0+dcosh(r0/a))/(dcosh(r/a)+dcosh(r0/a))
fV=Vc+Vn+hh**2/2.0/u*(L+0.5)**2/r**2
dh2=(r-r1)/m
res=dh2/2.0*(f2(r,r0,z0,z1,e,L,u,hh,V0,Q,a)+f2(r1,r0,z0,z1,e,L,u,hh,V0,Q,a))
do j=1,m-1
    res=res+dh2*f2(r1+j*dh2,r0,z0,z1,e,L,u,hh,V0,Q,a)
enddo
res=res-PI/4.0
f3=(2*u/hh/hh*abs(Q-fv))**-0.5*dcos(res)**2
end


dh3=(rr(i,2)-rr(i,1))/m
res3=dh3/2.0*(f3(rr(i,2),r0(i),z0,z(i),e,L,u,hh,V0,Q(i),a,rr(i,1))+f3(rr(i,1),r0(i),z0,z(i),e,L,u,hh,V0,Q(i),a,rr(i,1)))
do j=1,m-1
    res3=res3+dh3*f3(rr(i,1)+j*dh2,r0(i),z0,z(i),e,L,u,hh,V0,Q(i),a,rr(i,1),rk)
enddo

写一个梯形积分的公式，嵌套调用。实质就是两层循环的累加。