Fortran - integral doble

Estoy intentando hacer una integral doble en fortran por el método de simpson pero no se como declarar los limites de integracion.El problema dice asi:
Calcular mediante la formula de Simpson la integral de x*exp(-x*y) siendo la region (0 <_ x <_2, u(x) < _y <_v(x)} y siendo
u(x) = x**2-2*x+2 y v(x) =3*exp(x) (para cada valor de x, c = u(x) y d = v(x)). Tomar n = 200 y utilizar dos subprogramas
funcion para calcular u(x) y v(x).
yo he hecho esto y me da error pero no se donde esta



program prosim
implicit none
integer::n,i
real*8::x,a,b,c,d,s,s1,s2,h
real*8,external::f,g_simp,u,v
n=200
a=0d0
b=2d0
c=u(x)
d=v(x)

h=(b-a)/dble(n)
s1=0d0
do i=1,n-1,2
x=a+dble(i)*h
s1=s1+g_simp(x,c,d,h)
end do

s2=0d0
do i=2,n-2,2
x=a+dble(i)*h
s2=s2+g_simp(x,c,d,h)
end do

s=(h/3d0)*(g_simp(a,c,d,h)+g_simp(b,c,d,h)+4*s1+2*s2)

print*, s
end program prosim



function g_simp(x,c,d,h)
implicit none
integer::j,m
real*8::h,s1,s2,c,d,k,g_simp,y,x
real*8,external::f,

m=2*int((d-c+2*h)/(2*h))

k=(d-c)/dble(m)

s1=0d0
do j=1,m-1,2
y=c+dble(j)*k
s1=s1+f(x,y)
end do

s2=0d0
do j=2,m-2,2
y=c+dble(j)*k
s2=s2+f(x,y)
end do

g_simp=(k/3d0)*(f(x,c)+f(x,d)+4*s1+2*s2)

end function g_simp



function f(x,y)
real*8::f,x,y
f=x*exp((-x)*y)
end function

function u(x)
real*8::u,x
u=x**2-2*x+2
end function

function v(x)
real*8::v,x
v=3*exp(x)
end function

Así, por de pronto, pon el programa al final, después de las funciones. Así compila en gfortran y si lo ejecutas devuelve un valor que ignoro si es correcto.

De todas formas, en el programa principal tienes una llamada a la función 'u' y otra a 'v' pasando como argumento la variable 'x', que en ese momento no parece tener ningún valor asignado.