program f90_dsimp
    implicit none
    integer::nx,ny,i,j
    real(8)::a,b,c,d
    real(8)::dx,dy,ss,xa,ya,xb,yb
    character(len=20)::dummy
    real(8),parameter::pi=4.0D0*atan(1.0D0)

    call getarg(1,dummy)
    read(dummy,*) nx
    call getarg(2,dummy)
    read(dummy,*) ny

    a=0.0D0
    b=1.0D0
    c=0.0D0
    d=1.0D0
    dx=(b-a)/dble(nx)
    dy=(d-c)/dble(ny)

    ss=0.0D0
    do i=1,ny
        ya=c+dy*dble(i-1)
        yb=c+dy*dble(i)
        do j=1,nx
            xa=a+dx*dble(j-1)
            xb=a+dx*dble(j)
            ss=ss+DSIMP(xa,ya,xb,yb)
        end do
    end do
    write(6,'(2(i6),2(e20.12))') nx,ny,ss,ss-(exp(1.0D0)-1.0)**2

stop

contains

    real(8) function DSIMP(xa,ya,xb,yb)
        real(8),intent(in)::xa,ya,xb,yb
        real(8)::xc,yc,f1,f2,f3,s1,s2,s3

        xc=0.5D0*(xa+xb)
        yc=0.5D0*(ya+yb)

        f1=FUNC(xa,ya)
        f2=FUNC(xa,yc)
        f3=FUNC(xa,yb)
        s1=(yb-ya)/6.0D0*(f1+4.0D0*f2+f3)

        f1=FUNC(xc,ya)
        f2=FUNC(xc,yc)
        f3=FUNC(xc,yb)
        s2=(yb-ya)/6.0D0*(f1+4.0D0*f2+f3)

        f1=FUNC(xb,ya)
        f2=FUNC(xb,yc)
        f3=FUNC(xb,yb)
        s3=(yb-ya)/6.0D0*(f1+4.0D0*f2+f3)

       DSIMP=(xb-xa)/6.0D0*(s1+4.0D0*s2+s3)

    end function DSIMP


    real(8) function FUNC(x,y)
        real(8),intent(in)::x,y
        FUNC=exp(x+y)
    end function FUNC

end program f90_dsimp
