下面两个程序，第一个是计算辛普森积分的程序，第二个是超松弛迭代法的程序，现在第一个程序要利用第二个程序的结果，请问如何实现。
1
subroutine qsimp(func,a,b,s)
integer jmax
real a,b,func,s,eps
external func
parameter(eps=1.e-6,jmax=20)
!uses trapzd
integer j
real os, ost, st
ost=-1.e30
os=-1.e30
do j=1,jmax
call trapzd(func,a,b,st,j)
s=(4.*st-ost)/3.
if(abs(s-os)<eps*abs(os))return
if(s==0..and.os==0..and.j>6)return
os=s
ost=st
enddo
pause 'too many steps in qsimp'
end subroutine qsimp

subroutine trapzd(func,a,b,s,n)
integer n
real a,b,s,func
external func
integer it,j
real del,sum,tnm,x
if(n==1)then
s=0.5*(b-a)*(func(a)+func(b))
else
it=2**(n-2)
tnm=it
del=(b-a)/tnm
x=a+0.5*del
sum=0.
do j=1,it
   sum=sum+func(x)
   x=x+del
   enddo
   s=0.5*(s+(b-a)*sum/tnm)
endif
end subroutine trapzd


program d3r3
!drive for routine qsimp
parameter(pio2=1.5707963)
external func,fint
a=0.0
b=pio2
write(*,'(1x,a)')&
           'integral of func computed with qsimp'
write(*,'(1x,a,f10.6)') 'actual value of integral is ',&
       fint(b)-fint(a)
call qsimp(func,a,b,s)
write(*,'(1x,a,f10.6)')&
      'result from routine qsimp is',s
end program
function func(x)
 func=(x**2)*(x**2-2.0)*sin(x)
end function func
function fint(x)
! integral of func
fint=4.0*x*((x**2)-7.0)*sin(x)&
        -((x**4)-14.0*(x**2)+28.0)*cos(x)
end function fint

2

subroutine ssor(a,n,b,x,eps,om,ii)
parameter(imax=10000)
real a(n,n),b(n),x(n)
integer i,j,ii
real r,rx
do i=1,n
   r=1/a(i,i)
   b(i)=b(i)*r
   do j=1,n
      a(i,j)=a(i,j)*r
    end do
end do
do ii=1,imax
   rx=0.0
   do i=1,n
      r=b(i)
      do j=1,n
         r=r-a(i,j)*x(j)
      enddo
      if (abs(r)>rx)rx=abs(r)
         x(i)=x(i)+om*r
      enddo
      if(om*rx<=eps) return 
 enddo
pause 'too many iterations'
end subroutine ssor

 program dir12
 ! drive program for routine ssor
 parameter(n=99,eps=1.e-7,om=1.3)
 parameter(v=0.5,u=0.03141593)
 dimension a(n,n),c(n,n),r(n),b(n),x(n)
a(1,1)=-(1+v*cos(0.5*u))**3-(1+v*cos(1.5*u))**3
a(1,2)=(1+v*cos(1.5*u))**3
a(99,98)=(1+v*cos(197*0.5*u))**3
a(99,99)=-(1+v*cos(197*0.5*u))**3-(1+v*cos(199*0.5*u))**3
do i=2,n-1
 a(i,i)=-(1+v*cos((i-0.5)*u))**3-(1+v*cos((i+0.5)*u))**3
 a(i,i+1)=(1+v*cos((i+0.5)*u))**3
 a(i,i-1)=(1+v*cos((i-0.5)*u))**3
 enddo
 do i=1,n
b(i)=3*u*v*(cos((i+0.5)*u)-cos((i-0.5)*u))
 enddo


 do i=1,n
    do j=1,n
    c(i,j)=a(i,j)
    enddo
enddo
do l=1,n
   r(l)=b(l)
   x(l)=0.0
enddo
call ssor(c,n,r,x,eps,om,ii)
write(*,*) 'solution vector'
write(*,'(1x,99f12.9)') (x(l),l=1,n)
!test results with original matrix
write(*,*) 'right-hand side vector'
write(*,'(1x,99f12.9)') (b(l),l=1,n)
write(*,*) 'result of matrix applied to sol''n vector'
do l=1,n
   b(l)=0.0
   do j=1,n
     b(l)=b(l)+a(l,j)*x(j)
    enddo
enddo
write(*,'(1x,99f12.9)') (b(l),l=1,n)
write(*,*) 'iteration times=',ii
end