module defsp
    implicit none
    real(8),parameter::sp1=-1.0D0
    real(8),parameter::sp2= 0.0D0
    real(8),parameter::sp3= 1.0D0
end module defsp

program f90_adsimp
    use defsp
    implicit none
    integer::NODT,NELT
    integer,allocatable::kakom(:,:)
    real(8),allocatable::x(:),y(:)
    character(len=50)::fnameR
    integer::i,j,k,ne
    real(8)::tarea,earea
    real(8)::value

    call getarg(1,fnameR)
    open(11,file=fnameR,status='old')
        read(11,*) NODT,NELT
        allocate(kakom(1:NELT,1:4))
        allocate(x(1:NODT),y(1:NODT))
        do ne=1,NELT
            read(11,*) (kakom(ne,j),j=1,4)
        end do
        do i=1,NODT
            read(11,*) x(i),y(i)
        end do
    close(11)

    ! Total area of the integrated domain
    tarea=0.0D0
    do ne=1,NELT
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        earea=0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        i=kakom(ne,1)
        j=kakom(ne,3)
        k=kakom(ne,4)
        earea=earea+0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        tarea=tarea+earea
    end do


    value=0.0D0
    do ne=1,NELT
        value=value+ASIMP(ne,kakom,x,y,NELT,NODT)
    end do

    write(6,'(2(i7),2(f15.12))') NODT,NELT,tarea,value

contains

    subroutine CALXYJ(ne,kakom,x,y,a,b,xip,yip,detJ,NELT,NODT)
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(NELT,4)
        real(8),intent(in)::x(NODT),y(NODT),a,b
        real(8),intent(out)::xip,yip,detJ

        integer :: i,j,k,l
        real(8) :: J11,J12,J21,J22
        real(8) :: N1,N2,N3,N4,dn1a,dn2a,dn3a,dn4a,dn1b,dn2b,dn3b,dn4b

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        l=kakom(ne,4)

        !(N and xip, yip)
        N1=0.25D0*(1.0D0-a)*(1.0D0-b)
        N2=0.25D0*(1.0D0+a)*(1.0D0-b)
        N3=0.25D0*(1.0D0+a)*(1.0D0+b)
        N4=0.25D0*(1.0D0-a)*(1.0D0+b)
        xip=N1*x(i)+N2*x(j)+N3*x(k)+N4*x(l)
        yip=N1*y(i)+N2*y(j)+N3*y(k)+N4*y(l)

        !(dN/da,dN/db)
        dn1a=-0.25D0*(1.0-b)
        dn2a= 0.25D0*(1.0-b)
        dn3a= 0.25D0*(1.0+b)
        dn4a=-0.25D0*(1.0+b)
        dn1b=-0.25D0*(1.0-a)
        dn2b=-0.25D0*(1.0+a)
        dn3b= 0.25D0*(1.0+a)
        dn4b= 0.25D0*(1.0-a)
        !Jacobi matrix and det(J)
        J11=dn1a*x(i)+dn2a*x(j)+dn3a*x(k)+dn4a*x(l)
        J12=dn1a*y(i)+dn2a*y(j)+dn3a*y(k)+dn4a*y(l)
        J21=dn1b*x(i)+dn2b*x(j)+dn3b*x(k)+dn4b*x(l)
        J22=dn1b*y(i)+dn2b*y(j)+dn3b*y(k)+dn4b*y(l)
        detJ=J11*J22-J12*J21
    end subroutine CALXYJ


    real(8) function ASIMP(ne,kakom,x,y,NELT,NODT)
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(NELT,4)
        real(8),intent(in)::x(NODT),y(NODT)
        real(8)::xip,yip,detJ,a,b
        real(8)::f1,f2,f3,dj1,dj2,dj3,s1,s2,s3

        a=sp1
        b=sp1;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj1,NELT,NODT);f1=FUNC(xip,yip)
        b=sp2;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj2,NELT,NODT);f2=FUNC(xip,yip)
        b=sp3;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj3,NELT,NODT);f3=FUNC(xip,yip)
        s1=1.0D0/3.0D0*(f1*dj1+4.0D0*f2*dj2+f3*dj3)

        a=sp2
        b=sp1;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj1,NELT,NODT);f1=FUNC(xip,yip)
        b=sp2;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj2,NELT,NODT);f2=FUNC(xip,yip)
        b=sp3;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj3,NELT,NODT);f3=FUNC(xip,yip)
        s2=1.0D0/3.0D0*(f1*dj1+4.0D0*f2*dj2+f3*dj3)

        a=sp3
        b=sp1;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj1,NELT,NODT);f1=FUNC(xip,yip)
        b=sp2;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj2,NELT,NODT);f2=FUNC(xip,yip)
        b=sp3;call CALXYJ(ne,kakom,x,y,a,b,xip,yip,dj3,NELT,NODT);f3=FUNC(xip,yip)
        s3=1.0D0/3.0D0*(f1*dj1+4.0D0*f2*dj2+f3*dj3)

        ASIMP=1.0D0/3.0D0*(s1+4.0D0*s2+s3)
    end function ASIMP


    real(8) function FUNC(x,y)
        real(8),intent(in)::x,y
        real(8)::d
!        FUNC=1.0D0
        d=1.0D0-x*x-y*y
        FUNC=0.0D0
        if(0.0D0<d)FUNC=3.0D0/2.0D0*sqrt(d)
    end function FUNC

end program f90_adsimp
