Numerical Methods: Simpsons 1/3 rule for double integral

PLEASE DO NOT TAKE THIS QUESTION IF YOU ARE NOT KNOWLEDGABLE ON THE SUBJECT. I am trying to learn, not get free answers, and would greatly appreacite a thorough explanation of the below code.

Question: We are learning Simpons 1/3 rule for double integrals and our professor gives us some code to look at in Fortan that will aid us in later assignments. I have no idea what he is doing in the last do loop where he says "First integrate w.r.t to x". Can someone please explain this section of code.

subroutine SimpsonDIntRuleForFunc0(a,b,m,c,d,n) IMPLICIT NONE INTEGER :: i,n,m,j REAL :: x(m+1), y(n+1), INTa, a, b, c, d, h, k, fy0, fy1, fy2,SimpFXY h= (b-a)/m k= (d-c)/n x(1) = a DO i= 1, m x(i+1) = x(i)+h ENDDO print*, "x values, m = ",m do i = 1, m print*,(x(i)) end do y(1) = c DO i= 1, n y(i+1) = y(i)+k ENDDO print*, "y values n= ",n do i = 1, n print*,(y(i)) end do INTa = 0 Do i= 1,n-1,2 ! First integrate the function w.r.t. x using Simpson's 1/3 rule. fy0 = 0 fy1 = 0 fy2 = 0 DO j = 1,m-1,2 fy0 = fy0 + SimpFXY(x(j),y(i)) + 4*SimpFXY(x(j+1),y(i)) + SimpFXY(x(j+2),y(i)) fy1 = fy1 + SimpFXY(x(j),y(i+1)) + 4*SimpFXY(x(j+1),y(i+1)) + SimpFXY(x(j+2),y(i+1)) fy2 = fy2 + SimpFXY(x(j),y(i+2)) + 4*SimpFXY(x(j+1),y(i+2)) + SimpFXY(x(j+2),y(i+2)) ENDDO ! Integration step w.r.t. y using Simpson's 1/3 rule. INTa = INTa + fy0 + 4*fy1 + fy2 ENDDO INTa = INTa*h*k/9 print*, 'The value of the integral is: ', INTa END subroutine SimpsonDIntRuleForFunc0 FUNCTION SimpFXY(x,y) IMPLICIT NONE real::x,y,SimpFXY,e e = 2.7182818284590452 SimpFXY = 8 * e**(-x**2 - y**4) END FUNCTION SimpFXY