module defgp
    implicit none
    real(8),parameter::gp1=-0.774596669241483D0,gw1=0.555555555555555D0
    real(8),parameter::gp2= 0.000000000000000D0,gw2=0.888888888888889D0
    real(8),parameter::gp3= 0.774596669241483D0,gw3=0.555555555555555D0
end module defgp

program f90_dgint3
    use defgp
    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+DGINT(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 DGINT(xa,ya,xb,yb)
        real(8),intent(in)::xa,ya,xb,yb
        real(8)::x,y,f1,f2,f3,s

        x=0.5D0*(xb-xa)*gp1+0.5D0*(xa+xb)
        y=0.5D0*(yb-ya)*gp1+0.5D0*(ya+yb); f1=FUNC(x,y)
        y=0.5D0*(yb-ya)*gp2+0.5D0*(ya+yb); f2=FUNC(x,y)
        y=0.5D0*(yb-ya)*gp3+0.5D0*(ya+yb); f3=FUNC(x,y)

        s=gw1*(gw1*f1+gw2*f2+gw3*f3)

        x=0.5D0*(xb-xa)*gp2+0.5D0*(xa+xb)
        y=0.5D0*(yb-ya)*gp1+0.5D0*(ya+yb); f1=FUNC(x,y)
        y=0.5D0*(yb-ya)*gp2+0.5D0*(ya+yb); f2=FUNC(x,y)
        y=0.5D0*(yb-ya)*gp3+0.5D0*(ya+yb); f3=FUNC(x,y)
        s=s+gw2*(gw1*f1+gw2*f2+gw3*f3)

        x=0.5D0*(xb-xa)*gp3+0.5D0*(xa+xb)
        y=0.5D0*(yb-ya)*gp1+0.5D0*(ya+yb); f1=FUNC(x,y)
        y=0.5D0*(yb-ya)*gp2+0.5D0*(ya+yb); f2=FUNC(x,y)
        y=0.5D0*(yb-ya)*gp3+0.5D0*(ya+yb); f3=FUNC(x,y)
        s=s+gw3*(gw1*f1+gw2*f2+gw3*f3)

        DGINT=0.25D0*(xb-xa)*(yb-ya)*s
    end function DGINT


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

end program f90_dgint3
