implicit double precision (a-h, o-z)
      parameter(mx=512,my=512)
       dimension ff(0:mx,0:my)
       pi=atan(1.d0)*4.d0
       do 90 k=1,6
       nx=8*2**k
       ny=8*2**k
       dx=1.d0/dble(nx)
       dy=1.d0/dble(ny)
       do 10 i=0,nx 
        x=dble(i)*dx
        do 20 j=0,ny
          y=dble(j)*dy
          ff(i,j)=f(x,y)
   20   continue
   10  continue
       call simp2d(ff,s,dx,dy,nx,ny)
       write(*,*) nx,ny,s,s-4.d0/pi**2
   90  continue
      stop
      end 

      subroutine simp2d(ff,s,dx,dy,nx,ny)
      implicit double precision (a-h, o-z)
      parameter(mx=512,my=512)
       dimension ff(0:mx,0:my),ax(0:mx),ay(0:my)
       do 10 i=0,nx-2,2
        ax(i)=2.d0
        ax(i+1)=4.d0
   10  continue
        ax(0)=1.d0
        ax(nx)=1.d0
       do 20 j=0,ny-2,2
        ay(j)=2.d0
        ay(j+1)=4.d0
   20  continue
        ay(0)=1.d0
        ay(ny)=1.d0
       s=0.d0
       do 30 i=0,nx 
        do 40 j=0, ny
         s=s+ax(i)*ay(j)*ff(i,j)
   40   continue
   30  continue
       s=s*dx*dy/9.d0
       return
       end
      function f(x,y)
      implicit double precision (a-h, o-z)
        pi=atan(1.d0)*4.d0
        f=sin(pi*x)*sin(pi*y)
      return
      end