program f90_FEM_HEAT
    !*******************************************************
    !2D unsteady heat conduction analysis
    !*******************************************************
    implicit none
    integer::i,j,k,n,nnn,ia,ja,ne,it
    real(8)::alpc1,alpc2,alpc3,alpc4

    character(len=50)::ftinp         !Input file name for temperature history
    character(len=100)::strcom       !Comment
    integer::NODT                    !Number of nodes
    integer::NELT                    !Number of elements
    integer::MATEL                   !Number of material sets
    integer::KOT                     !Number of nodes with given temperature
    integer::KOC                     !Number of sides with the condition of heat transger boundary
    integer::nod=4                   !Number of nodes per element
    integer::nhen=1                  !Number of degree of freedom per node
    integer::nt                      !Number of degrees of freedom of all nodes (NODT x nhen)
    integer::mm                      !Number of degrees of freedom of reduced FE equation
    integer::ib                      !Band width of reduced matrix
    integer,allocatable::kakom(:,:)  !Element connectivity
    real(8),allocatable::Ak(:)       !Heat conductivity coefficient of element
    real(8),allocatable::Ac(:)       !Specific heat of element
    real(8),allocatable::Arho(:)     !Unit weight of element
    real(8),allocatable::Tk(:)       !maximum temperature rise of element
    real(8),allocatable::Al(:)       !Coefficient of heat release rate of element
    integer,allocatable::matno(:)    !Material set number
    real(8),allocatable::wm(:,:)     !work for material properties
    real(8),allocatable::x(:)        !x-coordinate of node
    real(8),allocatable::y(:)        !y-coordinate of node
    integer,allocatable::nokt(:)     !Node number with given temperature
    integer,allocatable::nekc(:,:)   !Element number and side number with the condition of heat transfer boundary
    real(8),allocatable::alphac(:)   !Heat transfer rate of side of element
    real(8),allocatable::Tinp(:)     !Temperature of given temperature node
    real(8),allocatable::Tcin1(:)    !Temperature of heat transfer boundary node (previous)
    real(8),allocatable::Tcin2(:)    !temperature of heat transfer boundary node (current)
    real(8)::alpc                    !work for heat transfer rate
    real(8)::tc                      !work for temperature of heat transfer boundary
    integer::kchen                   !work for side number of heat transfer boundary
    real(8)::dotq                    !work for heat generation rate
    real(8),allocatable::tck1(:,:)   !Total coefficient matrix (previous)
    real(8),allocatable::tck2(:,:)   !Total coefficient matrix (current)
    real(8),allocatable::rtck2(:,:)  !Total coefficient matrix (current, reduced matrix)
    real(8)::ck1(4,4)                !Element coefficient matrix (previous)
    real(8)::ck2(4,4)                !Element coefficient matrix (current)
    real(8)::ht(4,4)                 !Element heat transfer matrix
    real(8)::fq(4)                   !Element heat generation rate vector
    real(8)::fc(4)                   !Element heat transfer vector
    real(8),allocatable::tempe0(:)   !Initial temperatute
    real(8),allocatable::ftvec(:)    !Heat flux vector
    real(8),allocatable::tempe(:)    !Nodal temperature vector
    real(8),allocatable::reac(:)     !work for nodal temperature
    integer,allocatable::mindex(:)   !Index for treatment of given temperature nodes
    real(8)::delta                   !Time increment
    real(8)::ttime                   !Time
    integer::itmax                   !Number of time history step
    integer::nout                    !Number of temperature output nodes
    integer,allocatable::noout(:)    !Node number of temperature output nodes
    integer::ntout                   !Frequency of output of all node's temperatures
    integer,allocatable::nntout(:)   !Time step number of above
    real(8),allocatable::tempend(:,:)!Nodal temperature for output
    real(8)::fact0=1.0D-30           !0 judgement
    real(8)::elarea1                 !Half area of element
    real(8)::elarea2                 !Half area of element
    integer::io                      !Variable for finding 'EOF'
    character::fnameR*50,fnameW*50   !Input file name, output file name

    character::str*20
    character :: linebuf*1000
    character :: fmt1*200,fmt2*200
    real(4) :: t1,t2

    fmt1="(i5,2(',',e15.7),',',i2,',',e15.7)"               !Node
    fmt2="(i5,4(',',i5),4(',',e15.7),5(',',e15.7),',',i5)"  !Element properties

    call getarg(1,fnameR) !Input file name
    call getarg(2,fnameW) !Output file name

    !To start the measurement of calculation time
    t1=DIFT()
    !*****************************************************
    !Data input
    !*****************************************************
    open(11,file=fnameR,status='old')
    !--------------------------------------------------------------------------------------
    !Comment
    !--------------------------------------------------------------------------------------
    read(11,'(a)') strcom
    !--------------------------------------------------------------------------------------
    !Input file name for temperature history
    !--------------------------------------------------------------------------------------
    read(11,'(a)') ftinp
    !------------------------------------------------------------------------------------------
    !Basic values for analysis
    !------------------------------------------------------------------------------------------
    read(11,*) NODT,NELT,MATEL,KOT,KOC,delta
    !------------------------------------------------------------------------------------------
    !To declare array size
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:4))
    allocate(Ak(1:NELT))
    allocate(Ac(1:NELT))
    allocate(Arho(1:NELT))
    allocate(Tk(1:NELT))
    allocate(Al(1:NELT))
    allocate(matno(1:NELT))
    allocate(wm(1:MATEL,1:5))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(nokt(1:KOT))
    allocate(nekc(1:KOC,1:2))
    allocate(alphac(1:KOC))
    allocate(Tinp(1:KOT))
    allocate(Tcin1(1:KOC))
    allocate(Tcin2(1:KOC))
    allocate(mindex(1:nt))
    allocate(tck1(1:nt,1:nt))
    allocate(tck2(1:nt,1:nt))
    allocate(rtck2(1:nt,1:KOT))
    allocate(tempe0(1:nt))
    allocate(tempe(1:nt))
    allocate(reac(1:nt))
    allocate(ftvec(1:nt))
    !------------------------------------------------------------------------------------------
    !material properties (Ak,Ac,Arho,Tk,Al)
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3),wm(i,4),wm(i,5)
    end do
    !------------------------------------------------------------------------------------------
    !Element-nodes relationship, material set No.
    !------------------------------------------------------------------------------------------
    do ne=1,NELT
        read(11,*) kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4),matno(ne)
        do i=1,MATEL
            if(i==matno(ne))then
                Ak(ne)  =wm(i,1)
                Ac(ne)  =wm(i,2)
                Arho(ne)=wm(i,3)
                Tk(ne)  =wm(i,4)
                Al(ne)  =wm(i,5)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !(x,y) coordinates, initial temperature of node
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),y(i),tempe0(i)
    end do
    !------------------------------------------------------------------------------------------
    !Node number of given temperature nodes
    !------------------------------------------------------------------------------------------
    if(0<KOT)then
        do i=1,KOT
            read(11,*) nokt(i)
        end do
    end if
    !------------------------------------------------------------------------------------------
    !Element number and side number for heat transfer boundary, heat transfer rate
    !------------------------------------------------------------------------------------------
    if(0<KOC)then
        do i=1,KOC
            read(11,*) nekc(i,1),nekc(i,2),alphac(i)
        end do
    end if
    !------------------------------------------------------------------------------------------
    !Number of output node, Node number for output
    !------------------------------------------------------------------------------------------
    read(11,*) nout
    allocate(noout(1:nout))
    read(11,*) (noout(i),i=1,nout)
    !------------------------------------------------------------------------------------------
    !Frequency of outout of all node7s temperature, step number of output
    !------------------------------------------------------------------------------------------
    read(11,*) ntout
    if(0<ntout)then
        allocate(nntout(1:ntout))
        allocate(tempend(1:NODT,1:ntout))
        read(11,*) (nntout(i),i=1,ntout)
    end if
    close(11)

    !To know the number of time step
    open(13,file=ftinp,status='old')
    read(13,'()')!comment
    it=0
    do
        it=it+1
        read(13,'()',iostat=io) !if 'EOF' is found, exit
        if(io<0)exit
    end do
    itmax=it-1
    close(13)

    !*****************************************************
    !To print input data
    !*****************************************************
    open(12, file=fnameW, status='replace')
    write(12,'(a)') trim(strcom)
    write(12,'(a)') 'NODT,NELT,MATEL,KOT,KOC,delta,nout,ntout,itmax'
        write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOT,',',KOC,',',real(delta),',',nout,',',ntout,',',itmax
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,'(a)') '*node characteristics'
    write(12,'(a)') 'node,x,y,fix-T,Tinit'
    do i=1,NODT
        k=0
        if(0<KOT)then
            do j=1,KOT
                if(i==nokt(j))k=1
            end do
        end if
            write(linebuf,FMT=fmt1) i,x(i),y(i),k,tempe0(i)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(12,'(a)') '*element characteristics'
    write(12,'(a)') 'element,node-1,node-2,node-3,node-4,alpc1,alpc2,alpc3,alpc4,Ak,Ac,Arho,Tk,Al,matno'
    do ne=1,NELT
        alpc1=0.0D0
        alpc2=0.0D0
        alpc3=0.0D0
        alpc4=0.0D0
        if(0<KOC)then
            do k=1,KOC
                if(ne==nekc(k,1))then
                    if(nekc(k,2)==kakom(ne,1))alpc1=alphac(k)
                    if(nekc(k,2)==kakom(ne,2))alpc2=alphac(k)
                    if(nekc(k,2)==kakom(ne,3))alpc3=alphac(k)
                    if(nekc(k,2)==kakom(ne,4))alpc4=alphac(k)
                end if
            end do
        end if
            write(linebuf,FMT=fmt2) ne,kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4),&
                alpc1,alpc2,alpc3,alpc4,Ak(ne),Ac(ne),Arho(ne),Tk(ne),Al(ne),matno(ne)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !To verify area of element & to stop if negative area of element is found.
    do ne=1,NELT
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        elarea1=0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        if(elarea1<fact0)then
            write(12,'(a)') 'area of element is less or equal to 0'
            close(12)
            stop
        end if
        i=kakom(ne,1)
        j=kakom(ne,3)
        k=kakom(ne,4)
        elarea2=0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        if(elarea1<fact0)then
            write(12,'(a)') 'area of element is less or equal to 0'
            close(12)
            stop
        end if
    end do
    close(12)

    !*****************************************************
    !Main calculation
    !*****************************************************
    !-----------------------------------------------------
    !Initiarization of array
    !-----------------------------------------------------
    do i=1,nt
        do j=1,nt
            tck1(i,j)=0.0D0
            tck2(i,j)=0.0D0
        end do
    end do
    !-----------------------------------------------------
    !Assembling of total coefficient matrix
    !-----------------------------------------------------
    do ne=1,NELT
        Call STIFF(ne,kakom,x,y,Ak,Ac,Arho,ck1,ck2,delta,NELT,NODT)
        Call ASMAT(ne,nod,nhen,kakom,ck1,tck1,NELT,NODT)!Previous value: -0.5*[K]+1/dt*[C]
        Call ASMAT(ne,nod,nhen,kakom,ck2,tck2,NELT,NODT)!Current value:  +0.5*[K]+1/dt*[C]
        !Treatment of heat transfer boundary
        if(0<KOC)then
            do k=1,KOC
                if(ne==nekc(k,1))then
                    if(nekc(k,2)==kakom(ne,1))kchen=1
                    if(nekc(k,2)==kakom(ne,2))kchen=2
                    if(nekc(k,2)==kakom(ne,3))kchen=3
                    if(nekc(k,2)==kakom(ne,4))kchen=4
                    alpc=alphac(k)
                    Call HTMAT(ne,kchen,kakom,x,y,alpc,ht,NELT,NODT)!Heat transfer matrix
                    do i=1,4
                        do j=1,4
                            ht(i,j)=0.5D0*ht(i,j)
                        end do
                    end do
                    Call ASMAT(ne,4,1,kakom,ht,tck2,NELT,NODT)!current
                    do i=1,4
                        do j=1,4
                            ht(i,j)=-ht(i,j)
                        end do
                    end do
                    Call ASMAT(ne,nod,nhen,kakom,ht,tck1,NELT,NODT)!previous
                end if
            end do
        end if
    end do
    !-----------------------------------------------------
    !Memory of current total coefficient matrix
    !-----------------------------------------------------
    if(0<KOT)then
        do k=1,KOT
            do i=1,nt
                rtck2(i,k)=tck2(i,nokt(k))
            end do
        end do
    end if
    !-----------------------------------------------------
    !Treatment of given temperature nodes
    !-----------------------------------------------------
    do i=1,nt
        mindex(i)=i
    end do
    if(0<KOT)then
        do i=1,KOT
            n=nokt(i)
            mindex(n)=0
        end do
    end if
    mm=0
    do i=1,nt
        if(mindex(i)>0)then
            mm=mm+1
            mindex(mm)=i
        end if
    end do
    do i=1,mm
        ia=mindex(i)
        do j=1,mm
            ja=mindex(j)
            tck2(i,j)=tck2(ia,ja)
        end do
    end do

    !-----------------------------------------------------
    !calculation of band width and memory of band matrix
    !-----------------------------------------------------
    ib=IBAND(mm,tck2,fact0,nt)
    Call BAND(mm,ib,tck2,nt)
    !-----------------------------------------------------
    !To verify diagonal elements and to stop if necessary
    !-----------------------------------------------------
    do i=1,mm
        if(abs(tck2(i,1))<fact0)then
            open(12, file='fnameW.csv', access='append')
            write(12,'(a)') 'diagonal value is less or equal to 0'
            close(12)
            stop
        end if
    end do
    !-----------------------------------------------------
    !Cholesky method
    !-----------------------------------------------------
    Call CHBAND10(mm,ib,tck2,nt)!triangle matrix

    !-----------------------------------------------------
    !Preparation for time history calculation
    !-----------------------------------------------------
    !To set initial temperatures of nodes
    do i=1,nt
        tempe(i)=tempe0(i)
    end do
    !To set initial temperature for heat transfer boundary
    if(0<KOC)then
        do k=1,KOC
            ne=nekc(k,1)
            if(nekc(k,2)==kakom(ne,1))then
                ia=kakom(ne,1)
                ja=kakom(ne,2)
            end if
            if(nekc(k,2)==kakom(ne,2))then
                ia=kakom(ne,2)
                ja=kakom(ne,3)
            end if
            if(nekc(k,2)==kakom(ne,3))then
                ia=kakom(ne,3)
                ja=kakom(ne,4)
            end if
            if(nekc(k,2)==kakom(ne,4))then
                ia=kakom(ne,4)
                ja=kakom(ne,1)
            end if
            Tcin1(k)=0.5D0*(tempe(ia)+tempe(ja))
        end do
    end if

    !To open output file
    open(12, file=fnameW, access='append')
    write(12,'(a)') '*Temperature records of selected nodes'
    !======================
    ttime=0.0D0
    it=0
    !======================
    linebuf="it,time"
    do i=1,nout
        write (str,'(",node-",i5)') noout(i)
        linebuf=trim(linebuf)//str
    end do
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)

    write (linebuf,'(i5)') it
    write (str,'(",",e15.7)') ttime
    linebuf=trim(linebuf)//str
    do i=1,nout
        write (str,'(",",e15.7)') tempe(noout(i))
        linebuf=trim(linebuf)//str
    end do
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    !*****************************************************
    !Calculation of time history (untill end row of 'ftinp')
    !*****************************************************
    !To open the file of temperature time history
    open(13,file=ftinp,status='old')
    read(13,'()')!comment
    nnn=0
    !************************
    do it=1,itmax
    !************************
        ttime=delta*real(it)
        !-----------------------------------------------------
        !To input given temperatures each time
        !-----------------------------------------------------
        !Tepmerature of given temperature nodes, temperatures of nodes with heat transfer boundary
        if(KOT==0.and.KOC==0)then
            read(13,*) i
        end if
        if(0<KOT.and.KOC==0)then
            read(13,*) i,(Tinp(k),k=1,KOT)
        end if
        if(0<KOT.and.0<KOC)then
            read(13,*) i,(Tinp(k),k=1,KOT),(Tcin2(k),k=1,KOC)
        end if
        if(KOT==0.and.0<KOC)then
            read(13,*) i,(Tcin2(k),k=1,KOC)
        end if
        !-----------------------------------------------------
        !To assemble heat flux vector
        !-----------------------------------------------------
        !Previous temperature vector
        do i=1,nt
            reac(i)=0.0D0
            do j=1,nt
                reac(i)=reac(i)+tck1(i,j)*tempe(j)
            end do
        end do
        !Treatment of temperature given nodes
        if(0<KOT)then
            do k=1,KOT
                do i=1,nt
                    reac(i)=reac(i)-Tinp(k)*rtck2(i,k)
                end do
            end do
        end if
        !Heat transfer vector & heat generation rate vector
        do i=1,nt
            ftvec(i)=0.0D0
        end do
        do ne=1,NELT
            !Heat generation rate vector
            dotq=Arho(ne)*Ac(ne)*Tk(ne)*Al(ne)*exp(-Al(ne)*ttime)
            dotq=dotq+Arho(ne)*Ac(ne)*Tk(ne)*Al(ne)*exp(-Al(ne)*(ttime-delta))
            dotq=0.5*dotq
            Call FQVEC(ne,kakom,x,y,dotq,fq,NELT,NODT)
            Call ASVEC(ne,nod,nhen,kakom,fq,ftvec,NELT,NODT)
            !Heat transfer vector
            if(0<KOC)then
                do k=1,KOC
                    if(ne==nekc(k,1))then
                        if(nekc(k,2)==kakom(ne,1))kchen=1
                        if(nekc(k,2)==kakom(ne,2))kchen=2
                        if(nekc(k,2)==kakom(ne,3))kchen=3
                        if(nekc(k,2)==kakom(ne,4))kchen=4
                        alpc=alphac(k)
                        tc=0.5D0*(Tcin1(k)+Tcin2(k))
                        Call HTVEC(ne,kchen,kakom,x,y,alpc,tc,fc,NELT,NODT)
                        Call ASVEC(ne,nod,nhen,kakom,fc,ftvec,NELT,NODT)
                    end if
                end do
            end if
        end do
        !Total heat flux vector
        do i=1,nt
            reac(i)=reac(i)+ftvec(i)
        end do
        !-----------------------------------------------------
        !To reduce heat flux vector
        !-----------------------------------------------------
        do i=1,mm
            ia=mindex(i)
            reac(i)=reac(ia)
        end do
        !-----------------------------------------------------
        !Cholesky method
        !-----------------------------------------------------
        Call CHBAND11(mm,ib,tck2,reac,nt)
        !-----------------------------------------------------
        !Calculation of current temperature
        !-----------------------------------------------------
        do i=1,nt
            tempe(i)=0.0D0
        end do
        do i=1,mm
            ia=mindex(i)
            tempe(ia)=reac(i)
        end do
        do k=1,KOT
            tempe(nokt(k))=Tinp(k)
        end do
        !-----------------------------------------------------
        !Uodate of memory
        !-----------------------------------------------------
        if(0<KOC)then
            do k=1,KOC
                Tcin1(k)=Tcin2(k)
            end do
        end if
        !-----------------------------------------------------
        !To print results
        !-----------------------------------------------------
        write (linebuf,'(i5)') it
        write (str,'(",",e15.7)') ttime
        linebuf=trim(linebuf)//str
        do i=1,nout
            write (str,'(",",e15.7)') tempe(noout(i))
            linebuf=trim(linebuf)//str
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
        !-----------------------------------------------------
        !To save nodal temperatures
        !-----------------------------------------------------
        if(0<ntout)then
            do n=1,ntout
                if(it==nntout(n))then
                    nnn=nnn+1
                    do i=1,NODT
                        tempend(i,nnn)=tempe(i)
                    end do
                end if
            end do
        end if
    !************************
    end do
    !************************
    close(13)!To close the file of temperature time history

    if(0<ntout)then
        !To print nodal temperatures at specified time step
        write(12,'(a)') '*Temperature of all nodes for selected times'
        linebuf="Node,initial"
        do i=1,ntout
            write (str,'(",step-",i5)') nntout(i)
            linebuf=trim(linebuf)//str
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)

        do i=1,NODT
            write (linebuf,'(i5)') i
            write (str,'(",",e15.7)') tempe0(i)
            linebuf=trim(linebuf)//str
            do j=1,ntout
                write (str,'(",",e15.7)') tempend(i,j)
                linebuf=trim(linebuf)//str
            end do
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do

        !To print nodal temperatures to show temperature distribution
        write(12,'(a)') '*Temperature distributions of selected nodes'
        linebuf="Node,x,y,initial"
        do i=1,ntout
            write (str,'(",step-",i5)') nntout(i)
            linebuf=trim(linebuf)//str
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)

        do i=1,nout
            k=noout(i)
            write (linebuf,'(i5)') k
            write (str,'(",",e15.7)') x(k)
            linebuf=trim(linebuf)//str
            write (str,'(",",e15.7)') y(k)
            linebuf=trim(linebuf)//str
            write (str,'(",",e15.7)') tempe0(k)
            linebuf=trim(linebuf)//str
            do j=1,ntout
                write (str,'(",",e15.7)') tempend(k,j)
                linebuf=trim(linebuf)//str
            end do
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    !-----------------------------------------------------
    !To finish the measurement of calculation time
    !-----------------------------------------------------
    t2=DIFT()
    write(12,'("NODT=",i0," nt=",i0," mm=",i0," ib=",i0)') ,NODT,nt,mm,ib
    write(12,'("Calculation time=",f8.3,"(sec)")') t2-t1
    write(linebuf,*) 'Date_time='//trim(CALLDATE())//'_(gfortran)'
    write (12,'(a)') trim(adjustl(linebuf))
    close(12)!To close output file

stop

Contains

    subroutine NMAT(ne,kakom,x,y,dnk,dnc,detJ,a,b,NELT,NODT)
        !*******************************************************
        ! [N] matrix
        !*******************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::a,b
        real(8),intent(out)::detJ
        real(8),intent(out)::dnk(1:4,1:4),dnc(1:4,1:4)

        integer::i,j,k,l
        real(8)::J11,j12,J21,J22
        real(8)::dn1a,dn2a,dn3a,dn4a
        real(8)::dn1b,dn2b,dn3b,dn4b
        real(8)::dn1x,dn2x,dn3x,dn4x
        real(8)::dn1y,dn2y,dn3y,dn4y
        real(8)::nm1,nm2,nm3,nm4

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

        ![dN/da][dN/db]
        dn1a=-0.25D0*(1.0D0-b)
        dn2a= 0.25D0*(1.0D0-b)
        dn3a= 0.25D0*(1.0D0+b)
        dn4a=-0.25D0*(1.0D0+b)
        dn1b=-0.25D0*(1.0D0-a)
        dn2b=-0.25D0*(1.0D0+a)
        dn3b= 0.25D0*(1.0D0+a)
        dn4b= 0.25D0*(1.0D0-a)
        !Jacobi matrix & 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
        ![dN/dx][dN/dy]
        dn1x=( J22*dn1a-J12*dn1b)/detJ
        dn2x=( J22*dn2a-J12*dn2b)/detJ
        dn3x=( J22*dn3a-J12*dn3b)/detJ
        dn4x=( J22*dn4a-J12*dn4b)/detJ
        dn1y=(-J21*dn1a+J11*dn1b)/detJ
        dn2y=(-J21*dn2a+J11*dn2b)/detJ
        dn3y=(-J21*dn3a+J11*dn3b)/detJ
        dn4y=(-J21*dn4a+J11*dn4b)/detJ
        ![N]
        nm1=0.25D0*(1.0D0-a)*(1.0D0-b)
        nm2=0.25D0*(1.0D0+a)*(1.0D0-b)
        nm3=0.25D0*(1.0D0+a)*(1.0D0+b)
        nm4=0.25D0*(1.0D0-a)*(1.0D0+b)

        dnk(1,1)=dn1x*dn1x+dn1y*dn1y
        dnk(1,2)=dn1x*dn2x+dn1y*dn2y
        dnk(1,3)=dn1x*dn3x+dn1y*dn3y
        dnk(1,4)=dn1x*dn4x+dn1y*dn4y
        dnk(2,1)=dn2x*dn1x+dn2y*dn1y
        dnk(2,2)=dn2x*dn2x+dn2y*dn2y
        dnk(2,3)=dn2x*dn3x+dn2y*dn3y
        dnk(2,4)=dn2x*dn4x+dn2y*dn4y
        dnk(3,1)=dn3x*dn1x+dn3y*dn1y
        dnk(3,2)=dn3x*dn2x+dn3y*dn2y
        dnk(3,3)=dn3x*dn3x+dn3y*dn3y
        dnk(3,4)=dn3x*dn4x+dn3y*dn4y
        dnk(4,1)=dn4x*dn1x+dn4y*dn1y
        dnk(4,2)=dn4x*dn2x+dn4y*dn2y
        dnk(4,3)=dn4x*dn3x+dn4y*dn3y
        dnk(4,4)=dn4x*dn4x+dn4y*dn4y

        dnc(1,1)=nm1*nm1
        dnc(1,2)=nm1*nm2
        dnc(1,3)=nm1*nm3
        dnc(1,4)=nm1*nm4
        dnc(2,1)=nm2*nm1
        dnc(2,2)=nm2*nm2
        dnc(2,3)=nm2*nm3
        dnc(2,4)=nm2*nm4
        dnc(3,1)=nm3*nm1
        dnc(3,2)=nm3*nm2
        dnc(3,3)=nm3*nm3
        dnc(3,4)=nm3*nm4
        dnc(4,1)=nm4*nm1
        dnc(4,2)=nm4*nm2
        dnc(4,3)=nm4*nm3
        dnc(4,4)=nm4*nm4

    end subroutine NMAT


    subroutine IPOINT4(kk,a,b)
        !*******************************************************
        !Gauss points for area integration
        !*******************************************************
        integer,intent(in)::kk
        real(8),intent(out)::a,b

        select case (kk)
            case (1)
                a=-0.5773502692D0
                b=-0.5773502692D0
            case (2)
                a= 0.5773502692D0
                b=-0.5773502692D0
            case (3)
                a= 0.5773502692D0
                b= 0.5773502692D0
            case (4)
                a=-0.5773502692D0
                b= 0.5773502692D0
        end select
    end subroutine IPOINT4


    subroutine STIFF(ne,kakom,x,y,Ak,Ac,Arho,ck1,ck2,delta,NELT,NODT)
        !*******************************************************
        !Element coefficient matrix
        !*******************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::Ak(1:NELT),Ac(1:NELT),Arho(1:NELT)
        real(8),intent(in)::delta
        real(8),intent(out)::ck1(1:4,1:4),ck2(1:4,1:4)

        integer::kk,i,j
        real(8)::dnk(1:4,1:4),dnc(1:4,1:4),ww1(1:4,1:4),ww2(1:4,1:4)
        real(8)::detJ,a,b

        do i=1,4
            do j=1,4
                ww1(i,j)=0.0D0
                ww2(i,j)=0.0D0
                ck1(i,j)=0.0D0
                ck2(i,j)=0.0D0
            end do
        end do
        do kk=1,4
            Call IPOINT4(kk,a,b)
            Call NMAT(ne,kakom,x,y,dnk,dnc,detJ,a,b,NELT,NODT)
            do i=1,4
                do j=1,4
                    ww1(i,j)=ww1(i,j)+Ak(ne)*dnk(i,j)*detJ
                    ww2(i,j)=ww2(i,j)+Arho(ne)*Ac(ne)*dnc(i,j)*detJ
                end do
            end do
        end do
        do i=1,4
            do j=1,4
                ck1(i,j)=-0.5D0*ww1(i,j)+1.0D0/delta*ww2(i,j)
                ck2(i,j)=+0.5D0*ww1(i,j)+1.0D0/delta*ww2(i,j)
            end do
        end do

    end subroutine STIFF


    subroutine LPOINT4(kk,kchen,a,b)
        !*******************************************************
        !Gauss point for line integration
        !*******************************************************
        integer,intent(in)::kk,kchen
        real(8),intent(out)::a,b

        select case(kchen)
            case(1)
                b=-1.0
               if(kk==1)a=-0.5773502692D0
               if(kk==2)a= 0.5773502692D0
            case(2)
               a=1.0
               if(kk==1)b=-0.5773502692D0
               if(kk==2)b= 0.5773502692D0
            case(3)
               b=1.0
               if(kk==1)a=-0.5773502692D0
               if(kk==2)a= 0.5773502692D0
            case(4)
                a=-1.0
                if(kk==1)b=-0.5773502692D0
                if(kk==2)b= 0.5773502692D0
        end select

    end subroutine LPOINT4


    subroutine HTMAT(ne,kchen,kakom,x,y,alpc,ht,NELT,NODT)
        !*******************************************************
        !Element heat transfer matrix
        !*******************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kchen,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::alpc
        real(8),intent(out)::ht(1:4,1:4)

        integer::i,j,kk
        real(8)::a,b,s,vn(1:4)

        select case (kchen)
            case(1)
                i=kakom(ne,1)
                j=kakom(ne,2)
            case(2)
                i=kakom(ne,2)
                j=kakom(ne,3)
            case(3)
                i=kakom(ne,3)
                j=kakom(ne,4)
            case(4)
                i=kakom(ne,4)
                j=kakom(ne,1)
        end select
        s=sqrt((x(j)-x(i))**2.0D0+(y(j)-y(i))**2.0D0)

        do i=1,4
            do j=1,4
                ht(i,j)=0.0D0
            end do
        end do
        do kk=1,2
            Call LPOINT4(kk,kchen,a,b)
            ![N]
            vn(1)=0.25D0*(1.0D0-a)*(1.0D0-b)
            vn(2)=0.25D0*(1.0D0+a)*(1.0D0-b)
            vn(3)=0.25D0*(1.0D0+a)*(1.0D0+b)
            vn(4)=0.25D0*(1.0D0-a)*(1.0D0+b)
            do i=1,4
                do j=1,4
                    ht(i,j)=ht(i,j)+0.5D0*s*alpc*vn(i)*vn(j)
                end do
            end do
        end do

    end subroutine HTMAT


    subroutine FQVEC(ne,kakom,x,y,dotq,fq,NELT,NODT)
        !*******************************************************
        !Heat generation rate vector
        !*******************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::dotq
        real(8),intent(out)::fq(1:4)

        integer::i,j,k,l,kk
        real(8)::detJ,a,b,vn(1:4)
        real(8)::J11,J12,J21,J22
        real(8)::dn1a,dn2a,dn3a,dn4a
        real(8)::dn1b,dn2b,dn3b,dn4b

        do i=1,4
            fq(i)=0.0D0
        end do
        do kk=1,4
            Call IPOINT4(kk,a,b)
            i=kakom(ne,1)
            j=kakom(ne,2)
            k=kakom(ne,3)
            l=kakom(ne,4)
            ![dN/da][dN/db]
            dn1a=-0.25D0*(1.0D0-b)
            dn2a= 0.25D0*(1.0D0-b)
            dn3a= 0.25D0*(1.0D0+b)
            dn4a=-0.25D0*(1.0D0+b)
            dn1b=-0.25D0*(1.0D0-a)
            dn2b=-0.25D0*(1.0D0+a)
            dn3b= 0.25D0*(1.0D0+a)
            dn4b= 0.25D0*(1.0D0-a)
            !Jacobi & 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
            ![N]
            vn(1)=0.25D0*(1.0D0-a)*(1.0D0-b)
            vn(2)=0.25D0*(1.0D0+a)*(1.0D0-b)
            vn(3)=0.25D0*(1.0D0+a)*(1.0D0+b)
            vn(4)=0.25D0*(1.0D0-a)*(1.0D0+b)
            do i=1,4
                fq(i)=fq(i)+dotq*vn(i)*detJ
            end do
        end do

    end subroutine FQVEC


    subroutine HTVEC(ne,kchen,kakom,x,y,alpc,tc,fc,NELT,NODT)
        !*******************************************************
        !Heat transfer vector
        !*******************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kchen,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::alpc,tc
        real(8),intent(out)::fc(1:4)

        integer::i,j,kk
        real(8)::a,b,vn(1:4),s

        select case (kchen)
            case(1)
                i=kakom(ne,1)
                j=kakom(ne,2)
            case(2)
                i=kakom(ne,2)
                j=kakom(ne,3)
            case(3)
                i=kakom(ne,3)
                j=kakom(ne,4)
            case(4)
                i=kakom(ne,4)
                j=kakom(ne,1)
        end select
        s=sqrt((x(j)-x(i))**2.0D0+(y(j)-y(i))**2.0D0)

        do i=1,4
            fc(i)=0.0D0
        end do
        do kk=1,2
            Call LPOINT4(kk,kchen,a,b)
            ![N]
            vn(1)=0.25D0*(1.0D0-a)*(1.0D0-b)
            vn(2)=0.25D0*(1.0D0+a)*(1.0D0-b)
            vn(3)=0.25D0*(1.0D0+a)*(1.0D0+b)
            vn(4)=0.25D0*(1.0D0-a)*(1.0D0+b)
            do i=1,4
                fc(i)=fc(i)+0.5D0*s*alpc*tc*vn(i)
            end do
        end do

    end subroutine HTVEC


    subroutine ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !*******************************************************
        !Assembly of vector
        !*******************************************************
        !(nod: number of nodes of 1 element, nhen: degree of freedom of 1 node)
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,nod,nhen
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::fevec(1:nod*nhen)
        real(8),intent(inout)::ftvec(1:NODT*nhen)

        integer :: i,k,ki,ks,kik,ksk

        do i=1,nod
            ki=(kakom(ne,i)-1)*nhen
            ks=(i-1)*nhen
            do k=1,nhen
                kik=ki+k
                ksk=ks+k
                ftvec(kik)=ftvec(kik)+fevec(ksk)
            end do
        end do
    end subroutine ASVEC


    subroutine ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
        !*******************************************************
        !Assembly of matrix
        !*******************************************************
        !(nod: number of nodes of 1 element, nhen: degree of freedom of 1 node)
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,nod,nhen
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::sm(1:nod*nhen,1:nod*nhen)
        real(8),intent(inout)::tsm(1:NODT*nhen,1:NODT*nhen)

        integer :: i,j,k,l,ki,kj,ks,js,kik,kjl,ksk,jsl

        do i=1,nod
            do j=1,nod
                ki=(kakom(ne,i)-1)*nhen
                kj=(kakom(ne,j)-1)*nhen
                ks=(i-1)*nhen
                js=(j-1)*nhen
                do k=1,nhen
                    do l=1,nhen
                        kik=ki+k
                        kjl=kj+l
                        ksk=ks+k
                        jsl=js+l
                        tsm(kik,kjl)=tsm(kik,kjl)+sm(ksk,jsl)
                    end do
                end do
            end do
        end do
    end subroutine ASMAT


    subroutine CHBAND10(n1,m1,array,nt)
        integer,intent(in)::n1,m1,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        integer :: i,j,k,mm,mn
        real(8) :: z

        array(1,1)=sqrt(array(1,1))
        do j=2,m1
            array(1,j)=array(1,j)/array(1,1)
        end do
        do i=2,n1
            if(n1-i+1>m1)then
                mm=m1
            else
                mm=n1-i+1
            end if
            do j=1,mm
                z=array(i,j)
                if(j/=m1)then
                    if(m1-j+1<i)then
                        mn=m1-j+1
                    else
                        mn=i
                    end if
                    do k=2,mn
                        z=z-array(i-k+1,k)*array(i-k+1,j+k-1)
                    end do
                    if(j>1)then
                        array(i,j)=z/array(i,1)
                    else
                        array(i,1)=sqrt(z)
                    end if
                else
                    array(i,j)=z/array(i,1)
                end if
            end do
        end do
    end subroutine CHBAND10


    subroutine CHBAND11(n1,m1,array,vector,nt)
        integer,intent(in)::n1,m1,nt
        real(8),intent(in)::array(1:nt,1:nt)
        real(8),intent(out)::vector(nt)
        integer :: i,j,mm
        real(8) :: z

        vector(1)=vector(1)/array(1,1)
        do i=2,n1
            z=vector(i)
            if(i>m1)then
                mm=m1
            else
                mm=i
            end if
            do j=2,mm
                z=z-array(i-j+1,j)*vector(i-j+1)
            end do
            vector(i)=z/array(i,1)
        end do

        vector(n1)=vector(n1)/array(n1,1)
        do i=n1-1,1,-1
            z=vector(i)
            do j=2,m1
                if(i+j-1>n1)exit
                z=z-array(i,j)*vector(i+j-1)
            end do
            vector(i)=z/array(i,1)
        end do
    end subroutine CHBAND11


    integer function IBAND(mm,array,fact0,nt)
        integer,intent(in)::mm,nt
        real(8),intent(in)::array(1:nt,1:nt),fact0
        integer :: i,j
        !-----------------------------------------------------
        !Calculation of band
        !-----------------------------------------------------
        IBAND=1
        do i=1,mm-1
            do j=mm,i+1,-1
                if(fact0<=abs(array(i,j)))exit
            end do
            if(IBAND<=j-i+1)IBAND=j-i+1
        end do
        if(mm<IBAND)IBAND=mm
    end function IBAND


    subroutine BAND(mm,ib,array,nt)
        integer,intent(in)::mm,ib,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        integer :: i,j,k,mn
        real(8) :: work(1:nt)
        !-----------------------------------------------------
        !Memory of band matrix
        !-----------------------------------------------------
        do i=2,mm
            if(i+ib-1<mm)then
                mn=ib
            else
                mn=mm-i+1
            end if
            do j=1,ib
                work(j)=0.0D0
            end do
            do k=1,mn
                work(k)=array(i,i+k-1)
            end do
            do j=1,ib
                array(i,j)=work(j)
            end do
        end do
    end subroutine BAND


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


    real(4) function DIFT()
        integer::count,count_rate,count_max
        call system_clock(count,count_rate,count_max)
        DIFT=real(count)/real(count_rate)
    end function DIFT


    character(len=50) function CALLDATE()
        character(len=8)::date
        character(len=10)::time
        character(len=5)::zone
        integer::values(1:8)
        integer::i
        character(len=4)::ye
        character(len=4)::da
        character(len=10)::mo
        character(len=5)::ti
        character(len=10)::str
        character(len=50)::datime

        call DATE_AND_TIME(date,time,zone,values)
        ye=date(1:4)
        mo=date(5:6)
        da=date(7:8)
        ti=time(1:2)//':'//time(3:4)

        if(trim(adjustl(mo))=='01')mo='January'
        if(trim(adjustl(mo))=='02')mo='February'
        if(trim(adjustl(mo))=='03')mo='March'
        if(trim(adjustl(mo))=='04')mo='April'
        if(trim(adjustl(mo))=='05')mo='May'
        if(trim(adjustl(mo))=='06')mo='June'
        if(trim(adjustl(mo))=='07')mo='July'
        if(trim(adjustl(mo))=='08')mo='August'
        if(trim(adjustl(mo))=='09')mo='September'
        if(trim(adjustl(mo))=='10')mo='October'
        if(trim(adjustl(mo))=='11')mo='November'
        if(trim(adjustl(mo))=='12')mo='December'

        CALLDATE=trim(adjustl(da))//' '//trim(adjustl(mo))//' '//trim(adjustl(ye))//' at '//trim(adjustl(ti))
    end function CALLDATE


end program f90_FEM_HEAT