module defparam
    implicit none
    integer,parameter::KBD = 500   !最大境界数
    integer,parameter::KTJ = 10000 !最大節点数
    integer,parameter::KCM = 100   !節点に集まる最大要素数
    integer,parameter::LTE = 100  !多角形の最大辺数
    integer,parameter::KTE = 2 * KTJ + 1
    integer::INPNODT               !指定節点数(入力データ)
    real(8),parameter::pi=3.14159265358979323846D0

end module defparam

!-----------------------------------------------------------------------------------------
program f90_fem_mesh3
!-----------------------------------------------------------------------------------------
    use defparam
    implicit none
    integer::ibex(1:KBD)         =0   !個々の外部境界上の節点数
    integer::ibin(1:KBD)         =0   !個々の内部境界上の節点数
    integer::ibno(1:KBD,1:KTJ)   =0   !外部境界上の節点番号
    integer::mindex(1:KBD+1)     =0   !作業領域
    real(8)::px(1:KTJ+3)         =0.0D0 !節点のx座標
    real(8)::py(1:KTJ+3)         =0.0D0 !節点のy座標
    integer::mtj(1:2*KTJ+1,1:3)  =0   !要素-節点関係
    integer::jac(1:2*KTJ+1,1:3)  =0   !要素-隣接要素関係
    integer::idm(1:2*KTJ+1)      =0   !要素-領域関係

    integer::NEX                 =0   !外部境界数
    integer::NIN                 =0   !内部境界数
    integer::NOB                 =0   !境界上の節点総数
    integer::NIB                 =0   !境界内部の節点総数
    integer::NODE                =0   !節点総数
    integer::NELM                =0   !要素総数
    character(len=50)::fnameR         !入力ファイル名
    character(len=50)::fnameW         !出力ファイル名

    call getarg(1,fnameR)
    call getarg(2,fnameW)
    call INPUT(NEX, NIN, ibex, ibin, ibno, NOB, NIB, px, py, fnameR)
    call MODEL(NEX, NIN, ibex, ibin, ibno, NOB, NIB, mindex, NODE, px, py, NELM, mtj, jac, idm)
    call DATA(NODE, NELM, mtj, jac, idm, px, py, fnameW)

stop

contains
!-----------------------------------------------------------------------------------------
    subroutine INPUT(NEX, NIN, ibex, ibin, ibno, NOB, NIB, px, py, fnameR)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NEX
        integer,intent(inout)::NIN
        integer,intent(inout)::ibex(1:KBD)
        integer,intent(inout)::ibin(1:KBD)
        integer,intent(inout)::ibno(1:KBD,1:KTJ)
        integer,intent(inout)::NOB
        integer,intent(inout)::NIB
        real(8),intent(inout)::px(1:KTJ+3)
        real(8),intent(inout)::py(1:KTJ+3)
        character(len=50)::fnameR

        integer::i,j

        open(11,file=fnameR,status='old')
        read(11,*) INPNODT
        read(11,*) NEX,NIN
        if(0 < NEX)then
            do i = 1,NEX
                read(11,*) ibex(i)
            end do
        end if
        if(0 < NIN)then
            do i = 1,NIN
                read(11,*) ibin(i)
            end do
        end if
        do i = 1,NEX
            read(11,*) (ibno(i,j),j=1,ibex(i))
        end do
        read(11,*) NOB,NIB
        do i = 1,NOB + NIB
            read(11,*) px(i),py(i)
        end do
        close(11)

    end subroutine INPUT
!-----------------------------------------------------------------------------------------
    subroutine MODEL(NEX, NIN, ibex, ibin, ibno, NOB, NIB, mindex, NODE, px, py, NELM, mtj, jac, idm)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NEX
        integer,intent(inout)::NIN
        integer,intent(inout)::ibex(1:KBD)
        integer,intent(inout)::ibin(1:KBD)
        integer,intent(inout)::ibno(1:KBD,1:KTJ)
        integer,intent(inout)::NOB
        integer,intent(inout)::NIB
        integer,intent(inout)::mindex(1:KBD+1)
        integer,intent(inout)::NODE
        real(8),intent(inout)::px(1:KTJ+3)
        real(8),intent(inout)::py(1:KTJ+3)
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::idm(1:2*KTJ+1)

        integer::i,ntp,ia,ib,ic
        real(8)::xmin,xmax,ymin,ymax,rax,ray,dmax
        integer::jnb(1:KTJ)        =0
        integer::nei(1:KTJ,1:KCM)  =0
        integer::ifix(1:KTJ)       =0
        real(8)::angl(1:2*KTJ+1)   =0.0D0
        integer::nedg(1:2*KTJ+1)   =0
        integer::kstack(1:2*KTJ+1) =0
        integer::list(1:KTJ)       =0
        integer::istack(1:KTJ)     =0
        integer::kv(1:2*KTJ+1)     =0
        integer::ibr(1:2*KTJ+1)    =0
        integer::iadres(1:KTJ+3)   =0
        integer::map(1:2*KTJ+1)    =0

        ntp = NOB + NIB
        xmin = px(1)
        xmax = xmin
        ymin = py(1)
        ymax = ymin
        do i = 2,ntp
            if(px(i) < xmin)xmin = px(i)
            if(xmax < px(i))xmax = px(i)
            if(py(i) < ymin)ymin = py(i)
            if(ymax < py(i))ymax = py(i)
        end do
        rax = xmax - xmin
        ray = ymax - ymin
        dmax = rax
        if(dmax < ray)dmax = ray
        do i = 1,ntp
            px(i) = (px(i) - xmin) / dmax
            py(i) = (py(i) - ymin) / dmax
        end do

        NELM = 1
        ia = KTJ + 1
        ib = KTJ + 2
        ic = KTJ + 3
        mtj(1, 1) = ia
        mtj(1, 2) = ib
        mtj(1, 3) = ic
        jac(1, 1) = 0
        jac(1, 2) = 0
        jac(1, 3) = 0

        px(ia) = -1.23D0
        py(ia) = -0.50D0
        px(ib) =  2.23D0
        py(ib) = -0.50D0
        px(ic) =  0.50D0
        py(ic) =  2.50D0

        call ROUGH(NEX, NIN, ibex, ibin, ibno, ntp, NIB, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, &
                  & list, iadres, istack, kv, ibr, map)

        do i = 1,NODE
            ifix(i) = 1
        end do

        call FINE(NODE, px, py, jnb, nei, ifix, NELM, mtj, jac, idm, istack, angl, nedg, kstack, map)
        if(NELM /= 2 * NODE + 1)then
            !Show ERROR mesage and stop the execution
            write(6,*) 'ERROR IN SUB MODEL: stop'
            stop
        end if

        call LAPLAS(NODE, ifix, mtj, px, py, jnb, nei)
        call REMOVE(NEX, mindex, NELM, mtj, jac, idm)
        call CHECK(NELM, mtj, jac)

        do i = 1,NODE
            px(i) = px(i) * dmax + xmin
            py(i) = py(i) * dmax + ymin
        end do

    end subroutine MODEL
!-----------------------------------------------------------------------------------------
    subroutine ROUGH(NEX, NIN, ibex, ibin, ibno, ntp, NIB, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, &
                     & list, iadres, istack, kv, ibr, map)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NEX
        integer,intent(in)::NIN
        integer,intent(in)::ibex(1:KBD)
        integer,intent(in)::ibin(1:KBD)
        integer,intent(in)::ibno(1:KBD,1:KTJ)
        integer,intent(in)::ntp
        integer,intent(in)::NIB
        integer,intent(inout)::NODE
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::idm(1:2*KTJ+1)
        integer,intent(inout)::list(1:KTJ)
        integer,intent(inout)::iadres(1:KTJ+3)
        integer,intent(inout)::istack(1:KTJ)
        integer,intent(inout)::kv(1:2*KTJ+1)
        integer,intent(inout)::ibr(1:2*KTJ+1)
        integer,intent(inout)::map(1:2*KTJ+1)

        integer::i,j,k,nb,np,ma,mb,mc,ms,mp,js,jz,loc,it
        real(8)::xs,ys,xa,ya

        do i = 1,NEX
            jz = 0
            nb = i
            np = ibex(nb)
            do j = 1,np
                list(j) = ibno(nb, j)
            end do
            call BOUGEN(jz, np, list, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack, kv, ibr, map, NODE)
            do k = 1,NELM
                if(idm(k) == jz)then
                    ma = iadres(mtj(k, 1))
                    mb = iadres(mtj(k, 2))
                    mc = iadres(mtj(k, 3))
                    ms = ma * mb * mc
                    mp = (mb - ma) * (mc - mb) * (ma - mc)
                    if((ms /= 0).and.(0 < mp))then
                        idm(k) = nb
                    end if
                end if
            end do
        end do

        do i = 1,NIN
            nb = i
            np = ibin(nb)
            do j = 1,np
                list(j) = NODE + j
            end do
            js = list(1)
            xs = px(js)
            ys = py(js)
            call LOCATE(xs, ys, px, py, mtj, jac, NELM, loc)
            jz = idm(loc)
            call BOUGEN(jz, np, list, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack, kv, ibr, map, NODE)

            do k = 1,NELM
                if(idm(k) == jz)then
                    ma = iadres(mtj(k, 1))
                    mb = iadres(mtj(k, 2))
                    mc = iadres(mtj(k, 3))
                    ms = ma * mb * mc
                    mp = (mb - ma) * (mc - mb) * (ma - mc)
                    if((ms /= 0).and.(mp < 0))then
                        idm(k) = 0
                    end if
                end if
            end do
        end do

        do i = 1,KTJ + 3
            iadres(i) = 0
        end do
        do i = 1,NIB
            NODE = NODE + 1
            xa = px(NODE)
            ya = py(NODE)
            call LOCATE(xa, ya, px, py, mtj, jac, NELM, it)
            jz = idm(it)
            call DELAUN(jz, NODE, NODE, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack)
        end do

        if(NODE /= ntp)then
            !Show ERROR mesage and stop the execution
            write(6,*) 'ERROR IN SUB ROUGH: stop'
            stop
        end if

    end subroutine ROUGH
!-----------------------------------------------------------------------------------------
    subroutine BOUGEN(jz, np, list, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack, kv, ibr, map, NODE)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::jz
        integer,intent(in)::np
        integer,intent(in)::list(1:KTJ)
        integer,intent(in)::ntp
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::idm(1:2*KTJ+1)
        integer,intent(inout)::iadres(1:KTJ+3)
        integer,intent(inout)::istack(1:KTJ)
        integer,intent(inout)::kv(1:2*KTJ+1)
        integer,intent(inout)::ibr(1:2*KTJ+1)
        integer,intent(inout)::map(1:2*KTJ+1)
        integer,intent(inout)::NODE

        integer::i,j,k,inp,js,ip,iq,iv

        inp = 0
        do i = 1,KTJ + 3
            iadres(i) = 0
        end do
        do i = 1,np
            iadres(list(i)) = i
        end do
        js = list(1)
        if(NODE < js)then
            inp = inp + 1
            call DELAUN(jz, js, js, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack)
        end if

        loopi: do i = 1,np
            ip = list(mod(i,np) + 1)
            iq = list(i)
            if((NODE < ip).and.(i /= np))then
                inp = inp + 1
                call DELAUN(jz, ip, ip, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack)
            end if
            do j = 1,jnb(iq)
                do k = 1,3
                    if(mtj(nei(iq, j), k) == ip)cycle loopi
                end do
            end do
            call SEARCH(jz, ip, iq, jnb, nei, NELM, mtj, jac, idm, istack, iv, kv, iadres, ibr)
            call POLY(iq, ip, iv, kv, px, py, NELM, mtj, jac, jnb, nei, map)
        end do loopi
        NODE = NODE + inp

    end subroutine BOUGEN
!-----------------------------------------------------------------------------------------
    subroutine DELAUN(jz, js, jg, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::jz
        integer,intent(in)::js
        integer,intent(in)::jg
        integer,intent(in)::ntp
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::idm(1:2*KTJ+1)
        integer,intent(in)::iadres(1:KTJ+3)
        integer,intent(inout)::istack(1:KTJ)

        integer::i,itop,maxstk,ip,it,ia,ib,ic,iv1,iv2,iv3,idf,ms,iedge,il,ir,iera,ierb,ierl,iswap
        real(8)::xp,yp

        itop = 0
        maxstk = ntp
        do i = js,jg
            ip = i
            xp = px(ip)
            yp = py(ip)
            call LOCATE(xp, yp, px, py, mtj, jac, NELM, it)
            if(idm(it) /= jz)then
                !Show ERROR mesage and stop the execution
                write(6,*) 'ERROR IN SUB DELAUN: stop'
                stop
            end if

            ia = jac(it, 1)
            ib = jac(it, 2)
            ic = jac(it, 3)
            iv1 = mtj(it, 1)
            iv2 = mtj(it, 2)
            iv3 = mtj(it, 3)
            mtj(it, 1) = ip
            mtj(it, 2) = iv1
            mtj(it, 3) = iv2
            jac(it, 1) = NELM + 2
            jac(it, 2) = ia
            jac(it, 3) = NELM + 1

            NELM = NELM + 1
            idm(NELM) = jz
            mtj(NELM, 1) = ip
            mtj(NELM, 2) = iv2
            mtj(NELM, 3) = iv3
            jac(NELM, 1) = it
            jac(NELM, 2) = ib
            jac(NELM, 3) = NELM + 1
            NELM = NELM + 1
            idm(NELM) = jz
            mtj(NELM, 1) = ip
            mtj(NELM, 2) = iv3
            mtj(NELM, 3) = iv1
            jac(NELM, 1) = NELM - 1
            jac(NELM, 2) = ic
            jac(NELM, 3) = it

            call INCR(iv1, NELM, jnb, nei)
            call INCR(iv2, NELM - 1, jnb, nei)
            if(iv3 <= KTJ)then
                nei(iv3, NEIBOR(iv3, it, jnb, nei)) = NELM - 1
            end if
            call INCR(iv3, NELM, jnb, nei)
            jnb(ip) = 3
            nei(ip, 1) = it
            nei(ip, 2) = NELM - 1
            nei(ip, 3) = NELM
            if(ia /= 0)then
                ms = iadres(iv1) * iadres(iv2)
                idf = abs(iadres(iv1) - iadres(iv2))
                if((idm(ia) == jz).and.((ms == 0).or.(idf /= 1)))then
                    itop = itop + 1
                    istack(itop) = IPUSH(it, maxstk, itop)
                end if
            end if
            if(ib /= 0)then
                call EDGE(ib, it, jac, iedge)
                jac(ib, iedge) = NELM - 1
                ms = iadres(iv2) * iadres(iv3)
                idf = abs(iadres(iv2) - iadres(iv3))
                if((idm(ib) == jz).and.((ms == 0).or.(idf /= 1)))then
                    itop = itop + 1
                    istack(itop) = IPUSH(NELM - 1, maxstk, itop)
                end if
            end if
            if(ic /= 0)then
                call EDGE(ic, it, jac, iedge)
                jac(ic, iedge) = NELM
                ms = iadres(iv3) * iadres(iv1)
                idf = abs(iadres(iv3) - iadres(iv1))
                if((idm(ic) == jz).and.((ms == 0).or.(idf /= 1)))then
                    itop = itop + 1
                    istack(itop) = IPUSH(NELM, maxstk, itop)
                end if
            end if
            do
                if(0 < itop)then
                    il = istack(itop)
                    itop = itop - 1
                    ir = jac(il, 2)
                    call EDGE(ir, il, jac, ierl)
                    iera = mod(ierl,3) + 1
                    ierb = mod(iera,3) + 1
                    iv1 = mtj(ir, ierl)
                    iv2 = mtj(ir, iera)
                    iv3 = mtj(ir, ierb)
                    call SWAP(px(iv1), py(iv1), px(iv2), py(iv2), px(iv3), py(iv3), xp, yp, iswap)
                    if(iswap == 1)then
                        ia = jac(ir, iera)
                        ib = jac(ir, ierb)
                        ic = jac(il, 3)
                        mtj(il, 3) = iv3
                        jac(il, 2) = ia
                        jac(il, 3) = ir
                        mtj(ir, 1) = ip
                        mtj(ir, 2) = iv3
                        mtj(ir, 3) = iv1
                        jac(ir, 1) = il
                        jac(ir, 2) = ib
                        jac(ir, 3) = ic
                        call DECR(iv1, il, jnb, nei)
                        call DECR(iv2, ir, jnb, nei)
                        call INCR(ip, ir, jnb, nei)
                        call INCR(iv3, il, jnb, nei)
                        if(ia /= 0)then
                            call EDGE(ia, ir, jac, iedge)
                            jac(ia, iedge) = il
                            ms = iadres(iv2) * iadres(iv3)
                            idf = abs(iadres(iv2) - iadres(iv3))
                            if((idm(ia) == jz).and.((ms == 0).or.(idf /= 1)))then
                                itop = itop + 1
                                istack(itop) = IPUSH(il, maxstk, itop)
                            end if
                        end if
                        if(ib /= 0)then
                            ms = iadres(iv3) * iadres(iv1)
                            idf = abs(iadres(iv3) - iadres(iv1))
                            if((idm(ib) == jz).and.((ms == 0).or.(idf /= 1)))then
                                itop = itop + 1
                                istack(itop) = IPUSH(ir, maxstk, itop)
                            end if
                        end if
                        if(ic /= 0)then
                            call EDGE(ic, il, jac, iedge)
                            jac(ic, iedge) = ir
                        end if
                    end if
                else
                    exit
                end if
            end do
        end do

    end subroutine DELAUN
!-----------------------------------------------------------------------------------------
    subroutine SEARCH(iz, ip, iq, jnb, nei, NELM, mtj, jac, idm, istack, iv, kv, iadres, ibr)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::iz
        integer,intent(in)::ip
        integer,intent(in)::iq
        integer,intent(in)::jnb(1:KTJ)
        integer,intent(in)::nei(1:KTJ,1:KCM)
        integer,intent(in)::NELM
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        integer,intent(in)::jac(1:2*KTJ+1,1:3)
        integer,intent(in)::idm(1:2*KTJ+1)
        integer,intent(inout)::istack(1:KTJ)
        integer,intent(inout)::iv
        integer,intent(inout)::kv(1:2*KTJ+1)
        integer,intent(in)::iadres(1:KTJ+3)
        integer,intent(inout)::ibr(1:2*KTJ+1)

        integer::i,j,ia,ib,ja,jb,mstk,nbr,idf,ms,ielm,jelm,kelm,mmin,jr

        iv = 0
        mstk = 0
        nbr = 0
        do i = 1,NELM
            kv(i) = 0
            ibr(i) = 0
        end do
        do i = 1,KTJ
            istack(i) = 0
        end do
        do i = 1,jnb(iq)
            nbr = nbr + 1
            ibr(nei(iq, i)) = nbr
        end do
        do i = 1,jnb(iq)
            ielm = nei(iq, i)
            j = IVERT(ielm, iq, mtj)
            ja = mod(j,3) + 1
            jb = mod(ja,3) + 1
            ia = mtj(ielm, ja)
            ib = mtj(ielm, jb)
            idf = abs(iadres(ia) - iadres(ib))
            ms = iadres(ia) * iadres(ib)
            if((idf == 1).and.(ms /= 0))cycle
            jelm = jac(ielm, ja)
            if(jelm == 0)cycle
            if(idm(jelm) /= iz)cycle
            nbr = nbr + 1
            ibr(jelm) = nbr
            if(mtj(jelm, 1) == ip)goto 80
            if(mtj(jelm, 2) == ip)goto 80
            if(mtj(jelm, 3) == ip)goto 80
            mstk = mstk + 1
            istack(mstk) = jelm
        end do
        do
            if(mstk == 0)then
                !Show ERROR mesage and stop the execution
                write(6,*) 'ERROR IN SUB SEARCH: stop'
                stop
            end if
            ielm = istack(1)
            mstk = mstk - 1
            do i = 1,mstk
                istack(i) = istack(i + 1)
            end do
            istack(mstk + 1) = 0
            do j = 1,3
                ja = mod(j,3) + 1
                ia = mtj(ielm, j)
                ib = mtj(ielm, ja)
                idf = abs(iadres(ia) - iadres(ib))
                ms = iadres(ia) * iadres(ib)
                if((idf == 1).and.(ms /= 0))cycle
                jelm = jac(ielm, j)
                if(jelm == 0)cycle
                if(ibr(jelm) /= 0)cycle
                if(idm(jelm) /= iz)cycle
                nbr = nbr + 1
                ibr(jelm) = nbr
                if(mtj(jelm, 1) == ip)goto 80
                if(mtj(jelm, 2) == ip)goto 80
                if(mtj(jelm, 3) == ip)goto 80
                mstk = mstk + 1
                istack(mstk) = jelm
            end do
        end do
80 continue
        do
            iv = iv + 1
            kv(iv) = jelm
            if(ibr(jelm) <= jnb(iq))exit
            mmin = KTE + 1
            do j = 1,3
                jr = jac(jelm, j)
                if(jr == 0)cycle
                if(ibr(jr) == 0)cycle
                ja = mod(j,3) + 1
                ia = mtj(jelm, j)
                ib = mtj(jelm, ja)
                idf = abs(iadres(ia) - iadres(ib))
                ms = iadres(ia) * iadres(ib)
                if((idf == 1).and.(ms /= 0))cycle
                if(ibr(jr) < mmin)then
                    kelm = jr
                    mmin = ibr(jr)
                end if
            end do
            if(mmin == KTE + 1)then
                !Show ERROR mesage and stop the execution
                write(6,*) 'ERROR IN SUB SEARCH: stop'
                stop
            end if
            jelm = kelm
        end do

    end subroutine SEARCH
!-----------------------------------------------------------------------------------------
    subroutine POLY(iq, ip, iv, kv, px, py, NELM, mtj, jac, jnb, nei, map)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::iq
        integer,intent(in)::ip
        integer,intent(in)::iv
        integer,intent(in)::kv(1:2*KTJ+1)
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(in)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)
        integer,intent(inout)::map(1:2*KTJ+1)

        integer::i,j,k,l,ia,ja,ir,jr,iv1,iv2,iedge,npa,npb,nta,ntb
        integer::ielm,jelm,kelm,ivx,ips,ipg,iva,ivb

        integer::iena(1:LTE,1:3)     =0
        integer::ienb(1:LTE,1:3)     =0
        integer::jeea(1:LTE,1:3)     =0
        integer::jeeb(1:LTE,1:3)     =0
        integer::nsra(1:LTE+2)       =0
        integer::nsrb(1:LTE+2)       =0
        integer::ihen(1:2*LTE+1,1:2) =0
        integer::jhen(1:2*LTE+1)     =0
        integer::iad(1:2*LTE+1)      =0
        integer::jstack(1:LTE+2)     =0

        if(iv == 2)then
            call EDGE(kv(1), kv(2), jac, ia)
            call EDGE(kv(2), kv(1), jac, ja)

            ir = jac(kv(1), mod(ia,3) + 1)
            jr = jac(kv(2), mod(ja,3) + 1)

            mtj(kv(1), mod(ia,3) + 1) = iq
            jac(kv(1), ia) = jr
            jac(kv(1), mod(ja,3) + 1) = kv(2)

            mtj(kv(2), mod(ja,3) + 1) = ip
            jac(kv(2), ja) = ir
            jac(kv(2), mod(ja,3) + 1) = kv(1)

            if(ir /= 0)then
                call EDGE(ir, kv(1), jac, iedge)
                jac(ir, iedge) = kv(2)
            end if
            if(jr /= 0)then
                call EDGE(jr, kv(2), jac, iedge)
                jac(jr, iedge) = kv(1)
            end if

            iv1 = mtj(kv(1), ia)
            iv2 = mtj(kv(2), ja)
            call DECR(iv1, kv(2), jnb, nei)
            call DECR(iv2, kv(1), jnb, nei)
            call INCR(iq, kv(1), jnb, nei)
            call INCR(ip, kv(2), jnb, nei)
        else
            npa = 0
            npb = 0
            do i = 1,LTE + 2
                nsra(i) = 0
                nsrb(i) = 0
            end do
            do i = 1,NELM
                map(i) = 0
            end do
            do i = 1,iv
                map(kv(1)) = 1
            end do
            do i = 1,iv
                ielm = kv(i)
                do j = 1,3
                    ivx = mtj(ielm, j)
                    call DECR(ivx, ielm, jnb, nei)
                end do
            end do

            call PICK(iq, ip, iv, kv, mtj, jac, map, npa, npb, nsra, nsrb)
            call SUBDIV(npa, nsra, px, py, nta, iena, jeea, ihen, jhen, iad, jstack)
            call SUBDIV(npb, nsrb, px, py, ntb, ienb, jeeb, ihen, jhen, iad, jstack)

            if(iv /= nta + ntb)then
                !Show ERROR mesage and stop the execution
                write(6,*) 'ERROR IN SUB POLY: stop'
                stop
            end if

            do i = 1,iv
                do j = 1,3
                    jac(kv(i), j) = 0
                end do
            end do
            do i = 1,nta
                ielm = kv(i)
                do j = 1,3
                    mtj(ielm, j) = iena(i, j)
                    if(jeea(i, j) /= 0)then
                        jac(ielm, j) = kv(jeea(i, j))
                    end if
                end do
            end do
            do i = 1,ntb
                ielm = kv(nta + i)
                do j = 1,3
                    mtj(ielm, j) = ienb(i, j)
                    if(jeeb(i, j) /= 0)then
                        jac(ielm, j) = kv(nta + jeeb(i, j))
                    end if
                end do
            end do
            do i = 1,iv
                ielm = kv(i)
                do j = 1,3
                    ivx = mtj(ielm, j)
                    call INCR(ivx, ielm, jnb, nei)
                end do
            end do

            do i = 1,iv
                ielm = kv(i)
                loopj: do j = 1,3
                    jelm = jac(ielm, j)
                    if(jelm /= 0)cycle
                    ips = mtj(ielm, j)
                    ipg = mtj(ielm, mod(j,3) + 1)
                    do k = 1,jnb(ipg)
                        kelm = nei(ipg, k)
                        do l = 1,3
                            iva = mtj(kelm, l)
                            ivb = mtj(kelm, mod(l,3) + 1)
                            if((iva == ipg).and.(ivb == ips))then
                                jac(ielm, j) = kelm
                                jac(kelm, l) = ielm
                                cycle loopj
                            end if
                        end do
                    end do
                end do loopj
            end do
        end if

    end subroutine POLY
!-----------------------------------------------------------------------------------------
    subroutine PICK(iq, ip, iv, kv, mtj, jac, map, npa, npb, nsra, nsrb)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::iq
        integer,intent(in)::ip
        integer,intent(in)::iv
        integer,intent(in)::kv(1:2*KTJ+1)
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        integer,intent(in)::jac(1:2*KTJ+1,1:3)
        integer,intent(in)::map(1:2*KTJ+1)
        integer,intent(inout)::npa
        integer,intent(inout)::npb
        integer,intent(inout)::nsra(1:LTE+2)
        integer,intent(inout)::nsrb(1:LTE+2)

        integer::i,ivx,jvx,jelm

        npa = 1
        nsra(npa) = ip
        ivx = IVERT(kv(1), ip, mtj)
        npa = npa + 1
        nsra(npa) = mtj(kv(1), mod(ivx,3) + 1)
        do i = 2,iv - 1
            jvx = IVERT(kv(i), nsra(npa), mtj)
            jelm = jac(kv(i), jvx)
            if(jelm == 0)then
                npa = npa + 1
                nsra(npa) = mtj(kv(i), mod(jvx,3) + 1)
            else if(map(jelm) == 0)then
                npa = npa + 1
                nsra(npa) = mtj(kv(i), mod(jvx,3) + 1)
            end if
        end do
        npa = npa + 1
        nsra(npa) = iq

        npb = 1
        nsrb(npb) = iq
        ivx = IVERT(kv(iv), iq, mtj)
        npb = npb + 1
        nsrb(npb) = mtj(kv(iv), mod(ivx,3) + 1)
        do i = iv - 1,2,-1
            jvx = IVERT(kv(i), nsrb(npb), mtj)
            jelm = jac(kv(i), jvx)
            if(jelm == 0)then
                npb = npb + 1
                nsrb(npb) = mtj(kv(i), mod(jvx,3) + 1)
            else if(map(jelm) == 0)then
                npb = npb + 1
                nsrb(npb) = mtj(kv(i), mod(jvx,3) + 1)
            end if
        end do
        npb = npb + 1
        nsrb(npb) = ip

        if(iv /= npa + npb - 4)then
            !Show ERROR mesage and stop the execution
            write(6,*) 'ERROR IN SUB PICK: stop'
            stop
        end if

    end subroutine PICK
!-----------------------------------------------------------------------------------------
    subroutine SUBDIV(npl, nsr, px, py, nte, ien, jee, ihen, jhen, iad, jstack)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::npl
        integer,intent(inout)::nsr(1:LTE+2)
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(inout)::nte
        integer,intent(inout)::ien(1:LTE,1:3)
        integer,intent(inout)::jee(1:LTE,1:3)
        integer,intent(inout)::ihen(1:2*LTE+1,1:2)
        integer,intent(inout)::jhen(1:2*LTE+1)
        integer,intent(inout)::iad(1:2*LTE+1)
        integer,intent(inout)::jstack(1:LTE+2)

        integer::i,j,k,nnpl,nbs,ia,ib,ic,ix
        real(8)::xa,ya,xb,yb,see

        nte = 0
        nnpl = npl
        do i = 1,LTE
            do j = 1,3
                ien(i, j) = 0
                jee(i, j) = 0
            end do
        end do

        do
            if(3 <= nnpl)then
                nbs = 1
                do
                    if(nbs <= nnpl - 1)then
                        ia = nsr(nbs)
                        ib = nsr(nbs + 1)
                        ic = nsr(mod((nbs + 1),nnpl) + 1)
                        xa = px(ib) - px(ia)
                        ya = py(ib) - py(ia)
                        xb = px(ic) - px(ia)
                        yb = py(ic) - py(ia)
                        see = xa * yb - xb * ya
                        if(1.0E-10 < see)then
                            nte = nte + 1
                            ien(nte, 1) = ia
                            ien(nte, 2) = ib
                            ien(nte, 3) = ic
                            nnpl = nnpl - 1
                            do i = nbs + 1,nnpl
                                nsr(i) = nsr(i + 1)
                            end do
                        end if
                        nbs = nbs + 1
                    else
                        exit
                    end if
                end do
            else
                exit
            end if
        end do

        ix = 0
        do i = 1,2 * LTE + 1
            ihen(i, 1) = 0
            ihen(i, 2) = 0
            jhen(i) = 0
            iad(i) = 0
        end do

        do i = 1,nte
            loopj: do j = 1,3
                ia = ien(i, j)
                ib = ien(i, mod(j,3) + 1)
                do k = 1,ix
                    if((ihen(k, 1) == ib).and.(ihen(k, 2) == ia))then
                        jee(i, j) = jhen(k)
                        jee(jhen(k), iad(k)) = i
                        cycle loopj
                    end if
                end do
                ix = ix + 1
                jhen(ix) = i
                iad(ix) = j
                ihen(ix, 1) = ia
                ihen(ix, 2) = ib
            end do loopj
        end do

        call LAWSON(nte, ien, jee, npl, px, py, jstack)

    end subroutine SUBDIV
!-----------------------------------------------------------------------------------------
    subroutine LAWSON(nte, ien, jee, npl, px, py, jstack)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::nte
        integer,intent(inout)::ien(1:LTE,1:3)
        integer,intent(inout)::jee(1:LTE,1:3)
        integer,intent(in)::npl
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(inout)::jstack(1:LTE+2)

        integer::i,j,itop,maxstk,ncount,ielm,il,ir,jl1,jl2,jl3,jr1,jr2,jr3
        integer::iv1,iv2,iv3,iv4,ia,ib,iswap,iedge
        real(8)::xx,yy

        itop = 0
        maxstk = npl
        ncount = 0
        do i = 1,nte
            ielm = i
            itop = itop + 1
            jstack(itop) = IPUSH(ielm, maxstk, itop)
        end do

        do
            if(0 < itop)then
                ncount = ncount + 1
                if(LTE < ncount)then
                    !Show ERROR mesage and stop the execution
                    write(6,*) 'ERROR IN SUB LAWSON: stop'
                    stop
                end if
                il = jstack(itop)
                itop = itop - 1
                do j = 1,3
                    jl1 = j
                    jl2 = mod(jl1,3) + 1
                    jl3 = mod(jl2,3) + 1
                    ir = jee(il, jl1)
                    if(ir == 0)cycle
                    iv1 = ien(il, jl1)
                    iv2 = ien(il, jl2)
                    iv3 = ien(il, jl3)
                    xx = px(ien(il, jl3))
                    yy = py(ien(il, jl3))
                    call EDGE(ir, il, jee, jr1)
                    jr2 = mod(jr1,3) + 1
                    jr3 = mod(jr2,3) + 1
                    iv4 = ien(ir, jr3)
                    call SWAP(px(iv2), py(iv2), px(iv1), py(iv1), px(iv4), py(iv4), xx, yy, iswap)
                    if(iswap == 1)then
                        ia = jee(il, jl2)
                        ib = jee(ir, jr2)
                        ien(il, jl2) = iv4
                        jee(il, jl1) = ib
                        jee(il, jl2) = ir
                        ien(ir, jr2) = iv3
                        jee(ir, jr1) = ia
                        jee(ir, jr2) = il
                        if(ia /= 0)then
                            call EDGE(ia, il, jee, iedge)
                            jee(ia, iedge) = ir
                        end if
                        if(ib /= 0)then
                            call EDGE(ib, ir, jee, iedge)
                            jee(ib, iedge) = il
                        end if
                        itop = itop + 1
                        jstack(itop) = IPUSH(il, maxstk, itop)
                        exit
                    end if
                end do
            else
                exit
            end if
        end do

    end subroutine LAWSON
!-----------------------------------------------------------------------------------------
    subroutine FINE(NODE, px, py, jnb, nei, ifix, NELM, mtj, jac,idm, istack, angl, nedg, kstack, map)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NODE
        real(8),intent(inout)::px(1:KTJ+3)
        real(8),intent(inout)::py(1:KTJ+3)
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)
        integer,intent(inout)::ifix(1:KTJ)
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::idm(1:2*KTJ+1)
        integer,intent(inout)::istack(1:KTJ)
        real(8),intent(inout)::angl(1:2*KTJ+1)
        integer,intent(inout)::nedg(1:2*KTJ+1)
        integer,intent(inout)::kstack(1:2*KTJ+1)
        integer,intent(inout)::map(1:2*KTJ+1)

        integer::i,j,k,l,ii,kk,mm
        real(8)::xx,gx,gy,ar,s,xc,yc,cgrax,cgray,xkp,ykp,drate,almt
        integer::ntotal,neib,ielm,jelm,kelm,lelm,melm,j1,j2,j3,nn,ls,le,lm,ied,jed,iv1,iv2,jedg

        ntotal = INPNODT
        if(ntotal <= NODE)ntotal = NODE
        almt = 1.0D0 / real(ntotal)
20 continue
        almt = 0.5D0 * almt
        neib = 0
!        do i = 1,KTE
!            angl(i) = 0.0
!            nedg(i) = 0
!            kstack(i) = 0
!        end do

        do i = 1,NELM
            if(idm(i) /= 0)then
                ielm = i
                neib = neib + 1
                kstack(neib) = ielm
                call DISTOR(ielm, mtj, px, py, drate, jedg)
                s = AREA(ielm, mtj, px, py)
                angl(ielm) = drate
                nedg(ielm) = jedg
            end if
        end do

        do i = 1,neib - 1
            k = kstack(i)
            xx = angl(k)
            ii = i + 1
            l = i
            do j = ii,neib
                kk = kstack(j)
                if(angl(kk) <= xx)cycle
                xx = angl(kk)
                l = j
            end do
            mm = kstack(i)
            kstack(i) = kstack(l)
            kstack(l) = mm
        end do

70 continue
        if(neib == 0)goto 20
        if(NODE < ntotal)then
            ielm = kstack(1)
            ied = nedg(ielm)
            jelm = jac(ielm, ied)
            call EDGE(jelm, ielm, jac, jed)
            iv1 = mtj(ielm, ied)
            iv2 = mtj(ielm, mod(ied,3) + 1)
            NODE = NODE + 1
            px(NODE) = 0.5D0 * (px(iv1) + px(iv2))
            py(NODE) = 0.5D0 * (py(iv1) + py(iv2))

            call REMESH(ielm, ied, jelm, jed, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, istack)

            if(idm(ielm) /= idm(jelm))then
                ifix(NODE) = 1
            else
                gx = 0.0D0
                gy = 0.0D0
                ar = 0.0D0
                do j = 1,jnb(NODE)
                    kelm = nei(NODE, j)
                    j1 = mtj(kelm, 1)
                    j2 = mtj(kelm, 2)
                    j3 = mtj(kelm, 3)
                    s = AREA(kelm, mtj, px, py)
                    xc = (px(j1) + px(j2) + px(j3)) / 3.0D0
                    yc = (py(j1) + py(j2) + py(j3)) / 3.0D0
                    ar = ar + s
                    gx = gx + s * xc
                    gy = gy + s * yc
                end do
                cgrax = gx / ar
                cgray = gy / ar
                xkp = px(NODE)
                ykp = py(NODE)
                px(NODE) = cgrax
                py(NODE) = cgray
                do j = 1,jnb(NODE)
                    kelm = nei(NODE, j)
                    s = AREA(kelm, mtj, px, py)
                    if(s <= 0.1 * almt)then
                        px(NODE) = xkp
                        py(NODE) = ykp
                        exit
                    end if
                end do
            end if

            do i = 1,NELM
                map(i) = 1
            end do
            do i = 1,jnb(NODE)
                lelm = nei(NODE, i)
                map(lelm) = 0
                angl(lelm) = 0.0D0
                nedg(lelm) = 0
            end do
            nn = 0
            do i = 1,neib
                if(map(kstack(i)) == 1)then
                    nn = nn + 1
                    kstack(nn) = kstack(i)
                end if
            end do
            neib = nn

            do i = 1,jnb(NODE)
                melm = nei(NODE, i)
                s = AREA(melm, mtj, px, py)
                if((almt < s).and.(idm(melm) /= 0))then
                    call DISTOR(melm, mtj, px, py, drate, jedg)
                    angl(melm) = drate
                    nedg(melm) = jedg
                    if(angl(kstack(1)) < drate)then
                        do j = 1,neib
                            kstack(neib - j + 2) = kstack(neib - j + 1)
                        end do
                        neib = neib + 1
                        kstack(1) = melm
                    else if(drate < angl(kstack(neib)))then
                        neib = neib + 1
                        kstack(neib) = melm
                    else
                        ls = 1
                        le = neib
                        do
                            if(le - ls /= 1)then
                                lm = int((ls + le) / 2)
                                if(angl(kstack(lm)) < drate)then
                                    le = lm
                                else
                                    ls = lm
                                end if
                            else
                                exit
                            end if
                        end do
                        do k = neib,le,-1
                            kstack(k + 1) = kstack(k)
                        end do
                        neib = neib + 1
                        kstack(le) = melm
                    end if
                end if
            end do
            goto 70
        end if

    end subroutine FINE
!-----------------------------------------------------------------------------------------
    subroutine REMESH(ielm, ied, jelm, jed, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, istack)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::ielm
        integer,intent(in)::ied
        integer,intent(in)::jelm
        integer,intent(in)::jed
        integer,intent(in)::NODE
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::idm(1:2*KTJ+1)
        integer,intent(inout)::istack(1:KTJ)

        integer::itop,maxstk,ip,ibj,jbj,ia,ib,ic,id,iv1,iv2,iv3,iv4,iedge,il,ir,iera,ierb,ierl,iswap
        real(8)::xp,yp

        itop = 0
        maxstk = NODE
        ip = NODE
        xp = px(ip)
        yp = py(ip)
        ibj = mod(ied,3) + 1
        jbj = mod(jed,3) + 1
        ia = jac(ielm, ibj)
        ib = jac(ielm, mod(ibj,3) + 1)
        ic = jac(jelm, jbj)
        id = jac(jelm, mod(jbj,3) + 1)
        iv1 = mtj(ielm, ibj)
        iv2 = mtj(ielm, mod(ibj,3) + 1)
        iv3 = mtj(jelm, jbj)
        iv4 = mtj(jelm, mod(jbj,3) + 1)
        mtj(ielm, 1) = ip
        mtj(ielm, 2) = iv1
        mtj(ielm, 3) = iv2
        jac(ielm, 1) = NELM + 2
        jac(ielm, 2) = ia
        jac(ielm, 3) = NELM + 1

        NELM = NELM + 1
        idm(NELM) = idm(ielm)
        mtj(NELM, 1) = ip
        mtj(NELM, 2) = iv2
        mtj(NELM, 3) = iv3
        jac(NELM, 1) = ielm
        jac(NELM, 2) = ib
        jac(NELM, 3) = jelm

        mtj(jelm, 1) = ip
        mtj(jelm, 2) = iv3
        mtj(jelm, 3) = iv4
        jac(jelm, 1) = NELM
        jac(jelm, 2) = ic
        jac(jelm, 3) = NELM + 1

        NELM = NELM + 1
        idm(NELM) = idm(jelm)
        mtj(NELM, 1) = ip
        mtj(NELM, 2) = iv4
        mtj(NELM, 3) = iv1
        jac(NELM, 1) = jelm
        jac(NELM, 2) = id
        jac(NELM, 3) = ielm
        if(iv1 <= KTJ)then
            nei(iv1, NEIBOR(iv1, jelm, jnb, nei)) = NELM
        end if
        call INCR(iv2, NELM - 1, jnb, nei)
        if(iv3 <= KTJ)then
            nei(iv3, NEIBOR(iv3, ielm, jnb, nei)) = NELM - 1
        end if
        call INCR(iv4, NELM, jnb, nei)
        jnb(ip) = 4
        nei(ip, 1) = ielm
        nei(ip, 2) = NELM - 1
        nei(ip, 3) = jelm
        nei(ip, 4) = NELM

        if(ia /= 0)then
            if(idm(ia) == idm(ielm))then
                itop = itop + 1
                istack(itop) = IPUSH(ielm, maxstk, itop)
            end if
        end if
        if(ib /= 0)then
            call EDGE(ib, ielm, jac, iedge)
            jac(ib, iedge) = NELM - 1
            if(idm(ib) == idm(NELM - 1))then
                itop = itop + 1
                istack(itop) = IPUSH(NELM - 1, maxstk, itop)
            end if
        end if
        if(ic /= 0)then
            if(idm(ic) == idm(jelm))then
                itop = itop + 1
                istack(itop) = IPUSH(jelm, maxstk, itop)
            end if
        end if
        if(id /= 0)then
            call EDGE(id, jelm, jac, iedge)
            jac(id, iedge) = NELM
            if(idm(id) == idm(NELM))then
                itop = itop + 1
                istack(itop) = IPUSH(NELM, maxstk, itop)
            end if
        end if

        do
            if(0 < itop)then
                il = istack(itop)
                itop = itop - 1
                ir = jac(il, 2)
                call EDGE(ir, il, jac, ierl)
                iera = mod(ierl,3) + 1
                ierb = mod(iera,3) + 1
                iv1 = mtj(ir, ierl)
                iv2 = mtj(ir, iera)
                iv3 = mtj(ir, ierb)
                call SWAP(px(iv1), py(iv1), px(iv2), py(iv2), px(iv3), py(iv3), xp, yp, iswap)
                if(iswap == 1)then
                    ia = jac(ir, iera)
                    ib = jac(ir, ierb)
                    ic = jac(il, 3)
                    mtj(il, 3) = iv3
                    jac(il, 2) = ia
                    jac(il, 3) = ir
                    mtj(ir, 1) = ip
                    mtj(ir, 2) = iv3
                    mtj(ir, 3) = iv1
                    jac(ir, 1) = il
                    jac(ir, 2) = ib
                    jac(ir, 3) = ic
                    call DECR(iv1, il, jnb, nei)
                    call DECR(iv2, ir, jnb, nei)
                    call INCR(ip, ir, jnb, nei)
                    call INCR(iv3, il, jnb, nei)
                    if(ia /= 0)then
                        call EDGE(ia, ir, jac, iedge)
                        jac(ia, iedge) = il
                        if(idm(ia) == idm(il))then
                            itop = itop + 1
                            istack(itop) = IPUSH(il, maxstk, itop)
                        end if
                    end if
                    if(ib /= 0)then
                        if(idm(ib) == idm(ir))then
                            itop = itop + 1
                            istack(itop) = IPUSH(ir, maxstk, itop)
                        end if
                    end if
                    if(ic /= 0)then
                        call EDGE(ic, il, jac, iedge)
                        jac(ic, iedge) = ir
                    end if
                end if
            else
                exit
            end if
        end do

    end subroutine REMESH
!-----------------------------------------------------------------------------------------
    subroutine LAPLAS(NODE, ifix, mtj, px, py, jnb, nei)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NODE
        integer,intent(in)::ifix(1:KTJ)
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        real(8),intent(inout)::px(1:KTJ+3)
        real(8),intent(inout)::py(1:KTJ+3)
        integer,intent(in)::jnb(1:KTJ)
        integer,intent(in)::nei(1:KTJ,1:KCM)

        integer::itera = 5
        integer::it,i,j,ielm,j1,j2,j3
        real(8)::gx,gy,ar,s,xc,yc,cgrax,cgray

        do it = 1,itera
            do i = 1,NODE
                if(ifix(i) == 0)then
                    gx = 0.0D0
                    gy = 0.0D0
                    ar = 0.0D0
                    do j = 1,jnb(i)
                        ielm = nei(i, j)
                        j1 = mtj(ielm, 1)
                        j2 = mtj(ielm, 2)
                        j3 = mtj(ielm, 3)
                        s = AREA(ielm, mtj, px, py)
                        xc = (px(j1) + px(j2) + px(j3)) / 3.0D0
                        yc = (py(j1) + py(j2) + py(j3)) / 3.0D0
                        ar = ar + s
                        gx = gx + s * xc
                        gy = gy + s * yc
                    end do
                    cgrax = gx / ar
                    cgray = gy / ar
                    px(i) = cgrax
                    py(i) = cgray
                end if
            end do
        end do

    end subroutine LAPLAS
!-----------------------------------------------------------------------------------------
    subroutine REMOVE(NEX, mindex, NELM, mtj, jac, idm)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NEX
        integer,intent(inout)::mindex(1:KBD+1)
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:2*KTJ+1,1:3)
        integer,intent(inout)::jac(1:2*KTJ+1,1:3)
        integer,intent(inout)::idm(1:2*KTJ+1)

        integer::i,j,k,l,iz,inelm,ielm,jelm,kelm,lelm,ikp,kedg,ledg
        integer::mkp(1:3),jkp(1:3)

        ielm = 0
        mindex(1) = 1
        do i = 1,NEX
            iz = i
            inelm = 0
            do j = mindex(i),NELM
                if(idm(j) == iz)then
                    ielm = ielm + 1
                    jelm = j
                    inelm = inelm + 1
                    if(ielm /= jelm)then
                        do k = 1,3
                            mkp(k) = mtj(ielm, k)
                            jkp(k) = jac(ielm, k)
                        end do
                        ikp = idm(ielm)
                        do k = 1,3
                            kelm = jac(jelm, k)
                            if(kelm /= 0)then
                                call EDGE(kelm, jelm, jac, kedg)
                                jac(kelm, kedg) = ielm + KTE + 1
                            end if
                        end do
                        do l = 1,3
                            lelm = jkp(l)
                            if(lelm /= 0)then
                                call EDGE(lelm, ielm, jac, ledg)
                                jac(lelm, ledg) = jelm + KTE + 1
                            end if
                        end do
                        do k = 1,3
                            jac(ielm, k) = mod(jac(ielm, k),(KTE + 1))
                            jac(jelm, k) = mod(jac(jelm, k),(KTE + 1))
                            kelm = jac(ielm, k)
                            lelm = jac(jelm, k)
                            do l = 1,3
                                if(0<kelm)jac(kelm, l) = mod(jac(kelm, l),(KTE + 1))
                                if(0<lelm)jac(lelm, l) = mod(jac(lelm, l),(KTE + 1))
                            end do
                        end do
                        do k = 1,3
                            jkp(k) = jac(ielm, k)
                        end do
                        do k = 1,3
                            mtj(ielm, k) = mtj(jelm, k)
                            jac(ielm, k) = jac(jelm, k)
                            mtj(jelm, k) = mkp(k)
                            jac(jelm, k) = jkp(k)
                        end do
                        idm(ielm) = idm(jelm)
                        idm(jelm) = ikp
                    end if
                end if
            end do
            mindex(i + 1) = mindex(i) + inelm
        end do

        do i = 1,ielm
            do j = 1,3
                if(ielm < jac(i, j))then
                    jac(i, j) = 0
                end if
            end do
        end do
        do i = ielm + 1,NELM
            do j = 1,3
                mtj(i, j) = 0
                jac(i, j) = 0
            end do
        end do

        NELM = ielm

    end subroutine REMOVE
!-----------------------------------------------------------------------------------------
    subroutine CHECK(NELM, mtj, jac)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NELM
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        integer,intent(in)::jac(1:2*KTJ+1,1:3)

        integer::i,j,k,ielm,jelm,kelm,ia,ib,ja,jb

        do i = 1,NELM
            loopj: do j = 1,3
                ielm = i
                ia = mtj(i, j)
                ib = mtj(i, mod(j,3) + 1)
                jelm = jac(i, j)
                if(jelm == 0)cycle
                do k = 1,3
                    kelm = jac(jelm, k)
                    if(kelm == ielm)then
                        ja = mtj(jelm, k)
                        jb = mtj(jelm, mod(k,3) + 1)
                        if((ia == jb).and.(ib == ja))cycle loopj
                        !Show ERROR mesage and stop the execution
                        write(6,*) 'ERROR IN SUB CHECK: stop'
                        stop
                    end if
                end do
                !Show ERROR mesage and stop the execution
                write(6,*) 'ERROR IN SUB CHECK: stop'
                stop
            end do loopj
        end do

    end subroutine CHECK
!-----------------------------------------------------------------------------------------
    subroutine DATA(NODE, NELM, mtj, jac, idm, px, py, fnameW)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NODE
        integer,intent(in)::NELM
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        integer,intent(in)::jac(1:2*KTJ+1,1:3)
        integer,intent(in)::idm(1:2*KTJ+1)
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        character(len=50),intent(in)::fnameW

        integer::i,j
        real(8)::aa

        open(12, file=fnameW, status='replace')
        write(12,'(2(i5))') NODE,NELM
        do i = 1,NELM
            write(12,'(4(i5))') (mtj(i,j),j=1,3),idm(i)
        end do
        do i = 1,NODE
            write(12,'(2(e15.7),i3)') px(i),py(i)
        end do
        ! Find irregular elements
        do i=1,NELM
            aa=AREA(i, mtj, px, py)
            if(aa<=0.0)write(12,'("Irregular element: No, area=",i5,", ",e15.7)'),i,aa
        end do

        close(12)

    end subroutine DATA
!-----------------------------------------------------------------------------------------
    subroutine INCR(n, l, jnb, nei)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::n
        integer,intent(in)::l
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)

        if(n <= KTJ)then
            jnb(n) = jnb(n) + 1
            if(KCM < jnb(n))then
                !Show ERROR mesage and stop the execution
                write(6,*) 'ERROR IN SUB INCR: stop'
                stop
            end if
            nei(n, jnb(n)) = l
        end if

    end subroutine INCR
!-----------------------------------------------------------------------------------------
    subroutine DECR(n, l, jnb, nei)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::n
        integer,intent(in)::l
        integer,intent(inout)::jnb(1:KTJ)
        integer,intent(inout)::nei(1:KTJ,1:KCM)

        integer::i,inb

        if(n <= KTJ)then
            inb = NEIBOR(n, l, jnb, nei)
            do i = inb,jnb(n) - 1
                nei(n, i) = nei(n, i + 1)
            end do
            nei(n, jnb(n)) = 0
            jnb(n) = jnb(n) - 1
            if(jnb(n) == 0)then
                !Show ERROR mesage and stop the execution
                write(6,*) 'ERROR IN SUB DECR: stop'
                stop
            end if
        end if

    end subroutine DECR
!-----------------------------------------------------------------------------------------
    subroutine DISTOR(n, mtj, px, py, drate, jedg)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::n
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        real(8),intent(inout)::drate
        integer,intent(inout)::jedg

        real(8)::dst(1:3)
        real(8)::xa,ya,xb,yb,xc,yc,ang1,ang2,ang3

        xa = px(mtj(n, 1))
        ya = py(mtj(n, 1))
        xb = px(mtj(n, 2))
        yb = py(mtj(n, 2))
        xc = px(mtj(n, 3))
        yc = py(mtj(n, 3))
        ang1 = THETA(xc, yc, xa, ya, xb, yb)
        ang2 = THETA(xa, ya, xb, yb, xc, yc)
        ang3 = pi - ang1 - ang2
        dst(1) = ang1 - pi / 3.0D0
        dst(2) = ang2 - pi / 3.0D0
        dst(3) = ang3 - pi / 3.0D0
        drate = abs(dst(1)) + abs(dst(2)) + abs(dst(3))
        jedg = 1
        if(dst(jedg) < dst(2))jedg = 2
        if(dst(jedg) < dst(3))jedg = 3

    end subroutine DISTOR
!-----------------------------------------------------------------------------------------
    subroutine LOCATE(xp, yp, px, py, mtj, jac, NELM, itri)
!-----------------------------------------------------------------------------------------
        real(8),intent(in)::xp
        real(8),intent(in)::yp
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        integer,intent(in)::jac(1:2*KTJ+1,1:3)
        integer,intent(in)::NELM
        integer,intent(inout)::itri

        integer::i,ia,ib

        itri = NELM
10 continue
        do i = 1,3
            ia = mtj(itri, i)
            ib = mtj(itri, mod(i,3) + 1)
            if((px(ia) - xp) * (py(ib) - yp) < (py(ia) - yp) * (px(ib) - xp))then
                itri = jac(itri, i)
                goto 10
            end if
        end do

    end subroutine LOCATE
!-----------------------------------------------------------------------------------------
    subroutine SWAP(x1, y1, x2, y2, x3, y3, xp, yp, iswap)
!-----------------------------------------------------------------------------------------
        real(8),intent(in)::x1
        real(8),intent(in)::y1
        real(8),intent(in)::x2
        real(8),intent(in)::y2
        real(8),intent(in)::x3
        real(8),intent(in)::y3
        real(8),intent(in)::xp
        real(8),intent(in)::yp
        integer,intent(inout)::iswap

        real(8)::x13,y13,x23,y23,x1p,y1p,x2p,y2p,cosa,cosb,sina,sinb

        x13 = x1 - x3
        y13 = y1 - y3
        x23 = x2 - x3
        y23 = y2 - y3
        x1p = x1 - xp
        y1p = y1 - yp
        x2p = x2 - xp
        y2p = y2 - yp
        cosa = x13 * x23 + y13 * y23
        cosb = x2p * x1p + y1p * y2p
        if((0.0D0 <= cosa).and.(0.0D0 <= cosb))then
            iswap = 0
        else if((cosa < 0.0D0).and.(cosb < 0.0D0))then
            iswap = 1
        else
            sina = x13 * y23 - x23 * y13
            sinb = x2p * y1p - x1p * y2p
            if((sina * cosb + sinb * cosa) < 0.0D0)then
                iswap = 1
            else
                iswap = 0
            end if
        end if

    end subroutine SWAP
!-----------------------------------------------------------------------------------------
    subroutine EDGE(l, k, work, iedge)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::l
        integer,intent(in)::k
        integer,intent(in)::work(:,:)
        integer,intent(inout)::iedge

        integer::i

        iedge = -1
        do i = 1,3
            if(work(l, i) == k)then
                iedge = i
                exit
            end if
        end do
        !Show ERROR mesage and stop the execution
        if(iedge == -1)then
            write(6,*) 'ERROR IN SUB EDGE: stop'
            stop
        end if

    end subroutine EDGE
!-----------------------------------------------------------------------------------------
    integer function IPUSH(item, maxstk, itop)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::item
        integer,intent(in)::maxstk
        integer,intent(in)::itop

        if(maxstk < itop)then
            !Show ERROR mesage and stop the execution
            write(6,*) 'ERROR IN FUNCTION IPUSH: stop'
            stop
        else
            IPUSH = item
        end if

    end function IPUSH
!-----------------------------------------------------------------------------------------
    integer function NEIBOR(n, l, jnb, nei)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::n
        integer,intent(in)::l
        integer,intent(in)::jnb(1:KTJ)
        integer,intent(in)::nei(1:KTJ,1:KCM)

        integer::i

        NEIBOR = -1
        do i = 1,jnb(n)
            if(nei(n, i) == l)then
                NEIBOR = i
                exit
            end if
        end do
        !Show ERROR mesage and stop the execution
        if(NEIBOR == -1)then
            write(6,*) 'ERROR IN FUNCTION NEIBOR: stop'
            stop
        end if

    end function NEIBOR
!-----------------------------------------------------------------------------------------
    integer function IVERT(l, k, mtj)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::l
        integer,intent(in)::k
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)

        integer::i

        IVERT = -1
        do i = 1,3
            if(mtj(l, i) == k)then
                IVERT = i
                exit
            end if
        end do
        !Show ERROR mesage and stop the execution
        if(IVERT == -1)then
            write(6,*) 'ERROR IN FUNCTION IVERT: stop'
            stop
        end if

    end function IVERT
!-----------------------------------------------------------------------------------------
    real(8) function AREA(ielm, mtj, px, py)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::ielm
        integer,intent(in)::mtj(1:2*KTJ+1,1:3)
        real(8),intent(in)::px(1:KTJ+3)
        real(8),intent(in)::py(1:KTJ+3)

        integer::j1,j2,j3

        j1 = mtj(ielm, 1)
        j2 = mtj(ielm, 2)
        j3 = mtj(ielm, 3)
        AREA=0.5*(px(j1)*py(j2)+px(j2)*py(j3)+px(j3)*py(j1)-px(j1)*py(j3)-px(j2)*py(j1)-px(j3)*py(j2))

    end function AREA
!-----------------------------------------------------------------------------------------
    real(8) function THETA(x0, y0, x1, y1, x2, y2)
!-----------------------------------------------------------------------------------------
        real(8),intent(in)::x0
        real(8),intent(in)::y0
        real(8),intent(in)::x1
        real(8),intent(in)::y1
        real(8),intent(in)::x2
        real(8),intent(in)::y2

        real(8)::err,xa,ya,xb,yb,prdin,prdex

        err = 1.0E-10
        xa = x1 - x0
        ya = y1 - y0
        xb = x2 - x0
        yb = y2 - y0
        prdin = xa * xb + ya * yb
        prdex = xa * yb - xb * ya
        if(abs(prdin) < err)then
            THETA = pi / 2.0D0
        else
            THETA = atan(prdex / prdin)
            if(THETA < 0.0D0)then
                THETA = THETA + pi
            end if
            if(pi < THETA)then
                THETA = THETA - pi
            end if
        end if

    end function THETA
!-----------------------------------------------------------------------------------------
end program f90_fem_mesh3
!-----------------------------------------------------------------------------------------
