program simpson2d
   implicit none
   integer,parameter::mx=512,my=512
   double precision,dimension(0:mx,0:my)::ff
   double precision pi,f,x,y,s,dx,dy
   integer nx,ny
   integer i,j,k
   pi=atan(1.0d0)*4.0d0
   do k=1,6
      nx=8*2**k
      ny=8*2**k
      dx=1.0d0/dble(nx)
      dy=1.0d0/dble(ny)
      do i=0,nx
         x=dble(i)*dx
         do j=0,ny
            y=dble(j)*dy
            ff(i,j)=f(x,y)
         end do
      end do
      call simp2d(ff,s,dx,dy,nx,ny)
      write(*,*) nx,ny,s,s-4.0d0/pi**2
   end do
   stop
   end 

subroutine simp2d(ff,s,dx,dy,nx,ny)
   implicit none
   integer,parameter::mx=512,my=512
   double precision,dimension(0:mx,0:my)::ff
   double precision,dimension(0:mx)::ax
   double precision,dimension(0:my)::ay
   double precision s,dx,dy
   integer nx,ny
   integer i,j
   do i=0,nx-2,2
      ax(i)=2.0d0
      ax(i+1)=4.0d0
   end do
   ax(0)=1.0d0
   ax(nx)=1.0d0
   do j=0,ny-2,2
      ay(j)=2.0d0
      ay(j+1)=4.0d0
   end do
   ay(0)=1.0d0
   ay(ny)=1.0d0
   s=0.0d0
   do i=0,nx 
      do j=0, ny
         s=s+ax(i)*ay(j)*ff(i,j)
      end do
   end do
   s=s*dx*dy/9.0d0
   return
   end

function f(x,y)
   implicit none
   double precision pi,f,x,y
   pi=atan(1.0d0)*4.0d0
   f=sin(pi*x)*sin(pi*y)
   return
   end