Вычисление определенного интеграла методом трапеций с заданной точностью

program hu
parameter eps=0.001
external fu
intrinsic cos
data m,n,a,b/20,25,0,6.28/
integer max
real s1,s2
 
do i=1,n-1
call in(a,b,n,0,fu,s1,it,mit)
enddo
10 CONTINUE
call in(a,b,n,m,fu,s2,it,mit)
IF (abs(s2-s1).gt.eps) GOTO 10
print*,'iteraziya=',it
print*,'maks it=',mit
end program hu
 
 
subroutine in(a,b,n,m,f1,s1,it,mit)
real s1
n=n+m
dx=(b-a)/n
s1=0
it=0
x1=a
do i=1,n
x2=x1+dx
s1=s1+(f1(x1)+f1(x2))*dx/2
x1=x2
it=i
mit=mit+i
enddo
end
 
function fu(x,f)
real f,x
f=1/(5-3*cos(x))
fu=f
end

program hu
parameter eps=0.001
data m,n,a,b/20,25,0,6.28/ 
intrinsic cos 
integer mit,iter,max,iterat
real s1,s2
 
call in(b,a,n,fu,s1,iter,max)
print*,'int1',s1,iter,max
 
iterat=0
10 CONTINUE
n=n+m
call in(b,a,n,fu,s2,iter,max)
iterat=iterat+1
IF (abs(s2-s1).gt.eps) GOTO 10
print*,'int2',s2,iter,max,iterat,abs(s2-s1)
 
end program hu
 
subroutine in(c,f,er,f1,s1,it,mit)
real s1,f1,c,f
integer it,mit
integer er
dx=(c-f)/er
s1=0
it=0
x1=f
do i=1,er
x2=x1+dx
s1=s1+(1/(5-3*cos(x1))+1/(5-3*cos(x2)))*dx/2
x1=x2
it=i
mit=mit+i
enddo
end

program aaa
real eps,a,b,s,dx,a1,a2
open(9,file='zad.txt',status='old')
read (9,*),a,b,n,m,eps
s1=1
s2=0
a1=a
a2=a
l=10
k=0
do while (abs(s1-s2).gt.eps)
if (k.lt.l) then
s1=0
s2=0
dx1=(b-a)/n
dx2=(b-a)/(n+m)
call bbb (dx1,dx2,a1,a2,n,m,s1,s2)
n=n+m
k=k+1
endif
enddo
print*,s1,s2,k
close(9)
end
subroutine bbb(dxx1,dxx2,aa,aaa,nn,mm,ss1,ss2)
integer nn,mm
real dxx,aa,ss,aaa
f(x)=1/(5-3*cos(x))
do i=1,nn
ss1=ss1+(f(aa)+f(aa+dxx1))*dxx1/2
aa=aa+dxx1
enddo
aa=aaa
do i=1,mm+nn
ss2=ss2+(f(aa)+f(aa+dxx2))*dxx2/2
aa=aa+dxx2
enddo
end subroutine