module defconst
    implicit none
    real(8),parameter::ggg=9.8D0
    real(8),parameter::pi=3.14159265358979323846D0
end module defconst


program f90_FEM_V2DIM
    !*****************************************************************
    !2D Stress Time History Response Analysis
    !*****************************************************************
    use defconst
    implicit none
    integer::i,j,k,l,m,kk,ne,n,ia,ja,iii
    character::strcom*50                !Comment
    integer::NSTRES                     !0: Plane strain, 1: plane stress
    integer::NODT                       !Number of nodes
    integer::NELT                       !Number of elements
    integer::MATEL                      !Number of material sets
    integer::KOX                        !Number of restricted nodes in x-direction
    integer::KOY                        !Number of restricted nodes in y-direction
    integer::NF                         !Number of loaded nodes
    integer,allocatable::kakom(:,:)     !Element connectivity
    integer,allocatable::matno(:)       !Material set number
    real(8),allocatable::wm(:,:)        !work for material properties
    real(8),allocatable::t(:)           !Thickness of element
    real(8),allocatable::Em(:)          !Elastic modulus of element
    real(8),allocatable::po(:)          !Poisson's ratio of element
    real(8),allocatable::gamma(:)       !Unit weight of element
    real(8),allocatable::zetaM(:)       !Constant for damping (proportional term to mass)
    real(8),allocatable::zetaK(:)       !Constant for damping (proportional term to stiffness)
    real(8),allocatable::x(:)           !x-coordinate of node
    real(8),allocatable::y(:)           !y-coordinate od node
    integer,allocatable::nokx(:)        !restricted node number in x-direction
    integer,allocatable::noky(:)        !restricted node number in y-direction
    real(8),allocatable::f(:)           !External force vector
    integer,allocatable::mindex(:)      !Index for treatment of restricted nodes
    real(8),allocatable::ak(:,:)        !Element stiffness matrix
    real(8),allocatable::am(:,:)        !Element mass matrix
    real(8),allocatable::ac(:,:)        !Element damping matrix
    real(8),allocatable::tak(:,:)       !Total stiffness matrix
    real(8),allocatable::tam(:,:)       !Total mass matrix
    real(8),allocatable::tac(:,:)       !Total damping matrix
    real(8),allocatable::tlm(:,:)       !Total left hand side matrix
    real(8),allocatable::acc2(:)        !Current nodal acceleration
    real(8),allocatable::vel2(:)        !Current nodal velocity
    real(8),allocatable::dis2(:)        !Current nodal displacement
    real(8),allocatable::acc1(:)        !Previous nodal acceleration
    real(8),allocatable::vel1(:)        !Previous nodal velocity
    real(8),allocatable::dis1(:)        !Previous nodal displacement
    real(8),allocatable::strsave(:,:)   !Stresses of element
    integer::ndata                      !Number of time history data
    real(8)::dt                         !Time increment
    real(8),allocatable::accinp(:)      !Time history data of accekeration
    real(8)::betaN=1.0D0/4.0D0          !coefficient of Newmark-beta-method (average acceleration method)
    integer::kprnod,kprfre              !Node number and displacement degree number for output
    real(8),allocatable::prdis(:)       !Memory of displacement for output
    real(8),allocatable::prstr(:,:)     !Memory of stresses for output
    real(8)::dismax                     !work for search for Max.displacement
    integer::kpriii                     !Memory of Max.disp. step
    real(8)::sigx,sigy,tauxy,ps1,ps2,ang
    real(8)::xg,yg

    real(8),allocatable::ftvec(:)       !Load vector
    real(8),allocatable::reac(:)        !work array
    real(8),allocatable::wva(:)         !work array
    real(8),allocatable::wvb(:)         !work array
    real(8),allocatable::wk1(:)         !work array
    real(8),allocatable::wk2(:)         !work array
    integer::nt                         !Number of degrees of freedom of all nodes (NODT x nhen)
    integer::nod                        !Number of nodes per element
    integer::nhen=2                     !Number of degree of freedom per node
    integer::mm                         !Number of degrees of freedom of reduced FE equation
    integer::ib                         !Band width of reduced matrix
    real(8)::fact0=1.0D-30              !0 judgement
    character::fnameR*50,fnameW*50,fnameA*50  !Input & output & acceleration input file name
    integer::numnd,numsr                !Number of displacement degrees and stress for output of time history data
    integer,allocatable::nndnum(:,:)    !Node number and displacement degree number for output of time history data
    integer,allocatable::nsrnum(:,:)    !Element number and stress degree number for output of time history data

    character :: linebuf*1000,dummy*200
    character :: fmt1*200,fmt23*200,fmt24*200,fmt3*200,fmt4*200
    real(4) :: t1,t2

    fmt1="(i5,2(',',e15.7),2(',',e15.7),2(',',i2))"             !Node
    fmt23="(i5,3(',',i5),6(',',e15.7),',',i5)"                  !Element
    fmt24="(i5,4(',',i5),6(',',e15.7),',',i5)"                  !Element
    fmt3="(i5,2(',',e15.7),2(',',e15.7))"                       !Displacement
    fmt4="(i5,2(',',e15.7),3(',',e15.7),3(',',e15.7),',',i5)"   !Stress

    call getarg(1,fnameR) !Input file name (structure data)
    call getarg(2,fnameW) !Output file name
    call getarg(3,fnameA) !Input file name for acceleration time history data

    !To start the measurement of calculation time
    t1=DIFT()
    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    !/*--------------------------------------------------------------------------------------*/
    !/*Comment*/
    !/*--------------------------------------------------------------------------------------*/
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
   !1Basic values for analysis
    !------------------------------------------------------------------------------------------
    read(11,*) nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES
    !------------------------------------------------------------------------------------------
    !To declare array size
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:MATEL,1:6))
    allocate(t(1:NELT))
    allocate(Em(1:NELT))
    allocate(po(1:NELT))
    allocate(gamma(1:NELT))
    allocate(zetaM(1:NELT))
    allocate(zetaK(1:NELT))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    if(1<=KOX)allocate(nokx(1:KOX))
    if(1<=KOY)allocate(noky(1:KOY))
    allocate(f(1:nt))
    allocate(mindex(1:nt))
    allocate(ak(1:nod*nhen,1:nod*nhen))
    allocate(am(1:nod*nhen,1:nod*nhen))
    allocate(ac(1:nod*nhen,1:nod*nhen))
    allocate(tak(1:nt,1:nt))
    allocate(tam(1:nt,1:nt))
    allocate(tac(1:nt,1:nt))
    allocate(tlm(1:nt,1:nt))
    allocate(acc2(1:nt))
    allocate(vel2(1:nt))
    allocate(dis2(1:nt))
    allocate(acc1(1:nt))
    allocate(vel1(1:nt))
    allocate(dis1(1:nt))
    allocate(strsave(1:NELT,1:3))
    allocate(ftvec(1:nt))
    allocate(reac(1:nt))
    allocate(wva(1:nt))
    allocate(wvb(1:nt))
    allocate(wk1(1:nt))
    allocate(wk2(1:nt))
    allocate(prdis(1:nt))
    allocate(prstr(1:NELT,1:3))
    !------------------------------------------------------------------------------------------
    !Material properties (t,Em,po,gamma,zetaM,zetaK)
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3),wm(i,4),wm(i,5),wm(i,6)
    end do
    !------------------------------------------------------------------------------------------
    !Element-nodes relationship, material set number
    !------------------------------------------------------------------------------------------
    do ne=1,NELT
        read(11,*) (kakom(ne,j),j=1,nod),matno(ne)
        do i=1,MATEL
            if(i==matno(ne))then
                t(ne)    =wm(i,1)
                Em(ne)   =wm(i,2)
                po(ne)   =wm(i,3)
                gamma(ne)=wm(i,4)
                zetaM(ne)=wm(i,5)
                zetaK(ne)=wm(i,6)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !(x,y) coordinates of node
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),y(i)
    end do
    !------------------------------------------------------------------------------------------
    !Restricted node number
    !------------------------------------------------------------------------------------------
    if(0<KOX)then
        do i=1,KOX
            read(11,*) nokx(i)
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            read(11,*) noky(i)
        end do
    end if
    !------------------------------------------------------------------------------------------
    !Loaded node number, load in x-direction, load in y-direction
    !------------------------------------------------------------------------------------------
    do i=1,nt
        f(i)=0.0D0
    end do
    if(0<NF)then
        do i=1,NF
            read(11,*) n,f(2*n-1),f(2*n)
        end do
    end if
    !------------------------------------------------------------------------------------------
    !Information of nodes for output of time history data
    !------------------------------------------------------------------------------------------
    read(11,*) numnd
    allocate(nndnum(1:numnd,1:2))
    do i=1,numnd
        read(11,*) (nndnum(i,j),j=1,2) !Node number, displacement degree number
    end do
    !------------------------------------------------------------------------------------------
    !Information of element for output of time history data of resultant force
    !------------------------------------------------------------------------------------------
    read(11,*) numsr
    allocate(nsrnum(1:numsr,1:2))
    do i=1,numsr
        read(11,*) (nsrnum(i,j),j=1,2) !Element number, stress degree number
    end do
    !------------------------------------------------------------------------------------------
    !Node number & displacement degree number for output of all displacements and all stresses
    !------------------------------------------------------------------------------------------
    read(11,*) kprnod,kprfre
    close(11)

    !*****************************************************
    !To input acceleration time history data
    !*****************************************************
    open(11,file=fnameA,status='old')
    read(11,'(a)') linebuf
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    allocate(accinp(1:ndata))
    do i=1,ndata
        read(11,*) accinp(i)
    end do
    close(11)


    !*****************************************************
    !To print input data
    !*****************************************************
    open(12, file=fnameW, status='replace')
    write(12,'(a)') trim(strcom)
    write(12,'(a)') 'nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,dt,ndata'
        write(linebuf,*) nod,',',NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',NF,',',NSTRES,',',dt,',',ndata
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,'(a)') '*node characteristics'
    write(12,'(a)') 'node,x,y,fx,fy,fix-x,fix-y'
    do i=1,NODT
        k=0
        l=0
        if(0<KOX)then
            do j=1,KOX
                if(i==nokx(j))k=1
            end do
        end if
        if(0<KOY)then
            do j=1,KOY
                if(i==noky(j))l=1
            end do
        end if
            write(linebuf,FMT=fmt1) i,x(i),y(i),f(2*i-2),f(2*i),k,l
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(12,'(a)') '*element characteristics'
    if(nod==3)then
        write(12,'(a)') 'element,node-1,node-2,node-3,t,E,po,gamma,zetaM,zetaK,matno'
        do ne=1,NELT
                write(linebuf,FMT=fmt23) ne,kakom(ne,1),kakom(ne,2),kakom(ne,3),&
                    t(ne),Em(ne),po(ne),gamma(ne),zetaM(ne),zetaK(ne),matno(ne)
                call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    if(nod==4)then
        write(12,'(a)') 'element,node-1,node-2,node-3,node-4,t,E,po,gamma,zetaM,zetaK,matno'
        do ne=1,NELT
                write(linebuf,FMT=fmt24) ne,kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4),&
                    t(ne),Em(ne),po(ne),gamma(ne),zetaM(ne),zetaK(ne),matno(ne)
                call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    close(12)


    !*****************************************************
    !main calculation
    !*****************************************************
    !-----------------------------------------------------
    !Initialization of matrices
    !-----------------------------------------------------
    do i=1,nt
        do j=1,nt
            tak(i,j)=0.0D0
            tam(i,j)=0.0D0
            tac(i,j)=0.0D0
            tlm(i,j)=0.0D0
        end do
    end do
    !-----------------------------------------------------
    !To assemble matrix
    !-----------------------------------------------------
    if(nod==3)then
        do ne=1,NELT
            !Stiffness matrix
            call AK_PL3(ne,kakom,x,y,NSTRES,Em,po,t,ak,NELT,NODT)
            call ASMAT(ne,nod,nhen,kakom,ak,tak,NELT,NODT)
            !Mass matrix
            call AM_PL3(ne,kakom,x,y,t,gamma,am,NELT,NODT)
            call ASMAT(ne,nod,nhen,kakom,am,tam,NELT,NODT)
            !Damping matrix
            call AC_PL3(ne,zetaM,zetaK,am,ak,ac,NELT)
            call ASMAT(ne,nod,nhen,kakom,ac,tac,NELT,NODT)
        end do
    end if
    if(nod==4)then
        do ne=1,NELT
            !Stiffness matrix
            call AK_PL4(ne,kakom,x,y,NSTRES,Em,po,t,ak,NELT,NODT)
            call ASMAT(ne,nod,nhen,kakom,ak,tak,NELT,NODT)
            !Mass matrix
            call AM_PL4(ne,kakom,x,y,t,gamma,am,NELT,NODT)
            call ASMAT(ne,nod,nhen,kakom,am,tam,NELT,NODT)
            !Damping matrix
            call AC_PL4(ne,zetaM,zetaK,am,ak,ac,NELT)
            call ASMAT(ne,nod,nhen,kakom,ac,tac,NELT,NODT)
        end do
    end if
    !To create left hand side matrix  [M]+dt/2*[C]+betaN*dt*dt*[K]
    do i=1,nt
        do j=1,nt
            tlm(i,j)=1.0D0/betaN/dt/dt*tam(i,j)+1.0D0/2.0D0/betaN/dt*tac(i,j)+tak(i,j)
        end do
    end do
    !-----------------------------------------------------
    !Treatment of restricted nodes
    !-----------------------------------------------------
    do i=1,nt
        mindex(i)=i
    end do
    if(0<KOX)then
        do i=1,KOX
            n=2*nokx(i)-1
            mindex(n)=0
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            n=2*noky(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)
            tlm(i,j)=tlm(ia,ja)
        end do
    end do
    !-----------------------------------------------------
    !To memory band matrix
    !-----------------------------------------------------
    ib=IBAND(mm,tlm,fact0,nt)
    call BAND(mm,ib,tlm,nt)
    !-----------------------------------------------------
    !To verify diagonal elements & stop if necessary
    !-----------------------------------------------------
    do i=1,mm
        if(abs(tlm(i,1))<fact0)stop
    end do
    !-----------------------------------------------------
    !Cholesky method (1)
    !-----------------------------------------------------
    call CHBAND10(mm,ib,tlm,nt)!triangle matrix

    do i=1,nt
        acc1(i)=0.0D0
        vel1(i)=0.0D0
        dis1(i)=0.0D0
        acc2(i)=0.0D0
        vel2(i)=0.0D0
        dis2(i)=0.0D0
    end do
    dismax=0.0D0
    kpriii=0
    !*****************************************************
    !calculation of time history response
    !*****************************************************
    open(12, file=fnameW, access='append')
    write(12,'(a)') '*time history response'
    linebuf='iii,time,accinp'
    do i=1,numnd
        write(dummy,'(",nod",i0,"_",i0,"a,nod",i0,"_",i0,"v,nod",i0,"_",i0,"d")') &
            nndnum(i,1),nndnum(i,2),nndnum(i,1),nndnum(i,2),nndnum(i,1),nndnum(i,2)
        linebuf=trim(linebuf)//dummy
    end do
    do i=1,numsr
        write(dummy,'(",elm",i0,"_",i0)') nsrnum(i,1),nsrnum(i,2)
        linebuf=trim(linebuf)//dummy
    end do
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)

    do iii=1,ndata
        !-----------------------------------------------------
        !Right hand side matrix
        !-----------------------------------------------------
        !(Average acceleration method)
        do i=1,nt
            wva(i)=1.0D0/betaN/dt/dt*dis1(i)+1.0D0/betaN/dt*vel1(i)+(1.0D0/2.0D0/betaN-1.0D0)*acc1(i)
            wvb(i)=1.0D0/2.0D0/betaN/dt*dis1(i)+(1.0D0/2.0D0/betaN-1.0D0)*vel1(i)+(1.0D0/4.0D0/betaN-1.0D0)*dt*acc1(i)
        end do
        !reac=[M]{f}
        !wk1=[M]{wva}
        !wk2=[C]{wvb}
        do i=1,nt
            reac(i)=0.0D0
            wk1(i)=0.0D0
            wk2(i)=0.0D0
            do j=1,nt
                reac(i)=reac(i)+tam(i,j)*f(j)
                wk1(i)=wk1(i)+tam(i,j)*wva(j)
                wk2(i)=wk2(i)+tac(i,j)*wvb(j)
            end do
        end do
        !-----------------------------------------------------
        !To sey load vector {ft}={f}+[M]{wva}+[C]{wvb}
        !-----------------------------------------------------
        do i=1,nt
            ftvec(i)=-accinp(iii)*reac(i)+wk1(i)+wk2(i)
        end do
        !-----------------------------------------------------
        !Cholesky method (2)
        !-----------------------------------------------------
        do i=1,mm
            ia=mindex(i)
            reac(i)=ftvec(ia)
        end do
        call CHBAND11(mm,ib,tlm,reac,nt)!Forward elimination & Backward substitution
        !-----------------------------------------------------
        !Calculation of acc, vel, disp and stress
        !-----------------------------------------------------
        do i=1,nt
            dis2(i)=0.0D0
        end do
        do i=1,mm
            ia=mindex(i)
            dis2(ia)=reac(i)
        end do
        do i=1,nt
            vel2(i)=1.0D0/2.0D0/betaN/dt*dis2(i)-1.0D0/2.0D0/betaN/dt*dis1(i)&
                  -(1.0D0/2.0D0/betaN-1.0D0)*vel1(i)-(1.0D0/4.0D0/betaN-1.0D0)*dt*acc1(i)
            acc2(i)=1.0D0/betaN/dt/dt*dis2(i)-1.0D0/betaN/dt/dt*dis1(i)&
                  -1.0D0/betaN/dt*vel1(i)-(1.0D0/2.0D0/betaN-1.0D0)*acc1(i)
        end do
        if(nod==3)then
            do ne=1,NELT
                call CALST_PL3(ne,kakom,x,y,dis2,NSTRES,Em,po,t,strsave,NELT,NODT,nt)
            end do
        end if
        if(nod==4)then
            do ne=1,NELT
                call CALST_PL4(ne,kakom,x,y,dis2,NSTRES,Em,po,t,strsave,NELT,NODT,nt)
            end do
        end if
        do i=1,nt
            acc1(i)=acc2(i)
            vel1(i)=vel2(i)
            dis1(i)=dis2(i)
        end do
        !To print time history data
        write(dummy,'(i0,",",e15.7,",",e15.7)') iii,dt*dble(iii),accinp(iii)
        linebuf=trim(dummy)
        do i=1,numnd
            k=2*(nndnum(i,1)-1)+nndnum(i,2)
            write(dummy,'(",",e15.7,",",e15.7,",",e15.7)') acc2(k),vel2(k),dis2(k)
            linebuf=trim(linebuf)//dummy
        end do
        do i=1,numsr
            ne=nsrnum(i,1)
            k =nsrnum(i,2)
            sigx =strsave(ne,1)
            sigy =strsave(ne,2)
            tauxy=strsave(ne,3)
            call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
            if(k==1)write(dummy,'(",",e15.7)') sigx
            if(k==2)write(dummy,'(",",e15.7)') sigy
            if(k==3)write(dummy,'(",",e15.7)') tauxy
            if(k==4)write(dummy,'(",",e15.7)') ps1
            if(k==5)write(dummy,'(",",e15.7)') ps2
            linebuf=trim(linebuf)//dummy
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
        !To update memory for output
        if(dismax<abs(dis2(2*(kprnod-1)+kprfre)))then
            dismax=abs(dis2(2*(kprnod-1)+kprfre))
            kpriii=iii
            do i=1,nt
                prdis(i)=dis2(i)
            end do
            do ne=1,NELT
                do j=1,3
                    prstr(ne,j)=strsave(ne,j)
                end do
            end do
        end if
    write(6,*) iii
    !*****************************************************
    end do
    !*****************************************************
    close(12)


    !*****************************************************
    !To print all displacements and all stresses at defined time step
    !*****************************************************
    open(12, file=fnameW, access='append')
    write(12,'(a)') '*displacements and stressresultant at max.displacement'
    write(12,'("kprnod=",i0," kprfre=",i0," kpriii=",i0," time=",e15.7)') ,kprnod,kprfre,kpriii,dt*dble(kpriii)
    !Displacement
    write(12,'(a)') '*displacements'
    write(12,'(a)') 'node,coord-x,coord-y,dis-x,dis-y'
    do i=1,NODT
        write(linebuf,FMT=fmt3) i,x(i),y(i),prdis(2*i-1),prdis(2*i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !Stress
    write(12,'(a)') '*stresses'
    write(12,'(a)') 'element,xg,yg,sigx,sigy,tauxy,ps1,ps2,ang,matno'
    do ne=1,NELT
        sigx=prstr(ne,1)
        sigy=prstr(ne,2)
        tauxy=prstr(ne,3)
        call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
        if(nod==3)call XGYG3M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
        if(nod==4)call XGYG4M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
        write(linebuf,FMT=fmt4) ne,xg,yg,sigx,sigy,tauxy,ps1,ps2,ang,matno(ne)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do

    !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)

stop

contains

    subroutine XGYG3M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
        !*************************************************************
        !Calculation of centroid for 3 node element
        !*************************************************************
        integer,intent(in)::ne,kk,NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::xg,yg

        integer::i,j,k

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        xg=(x(i)+x(j)+x(k))/3.0D0
        yg=(y(i)+y(j)+y(k))/3.0D0
    end subroutine XGYG3M


    subroutine XGYG4M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
        !*************************************************************
        !Calculation of centroid for 4 node element
        !*************************************************************
        integer,intent(in)::ne,kk,NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::xg,yg

        integer::i,j,k,l

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        l=kakom(ne,4)
        xg=0.25D0*(x(i)+x(j)+x(k)+x(l))
        yg=0.25D0*(y(i)+y(j)+y(k)+y(l))
    end subroutine XGYG4M


    subroutine DMAT_PL(NSTRES,Eme,poe,d)
        !*************************************************************
        !Stress-strain relationship [De]
        !*************************************************************
        integer,intent(in)::NSTRES
        real(8),intent(in)::Eme,poe
        real(8),intent(out)::d(1:3,1:3)

        integer i,j

        do i=1,3
            do j=1,3
                d(i,j)=0.0D0
            end do
        end do
        select case (NSTRES)
            case (0) !Plane strain
                d(1,1)=Eme*(1.0D0-poe)/((1.0D0+poe)*(1.0D0-2.0D0*poe))
                d(1,2)=Eme*poe/((1.0D0+poe)*(1.0D0-2.0D0*poe))
                d(2,1)=d(1,2)
                d(2,2)=d(1,1)
                d(3,3)=Eme*0.5D0/(1.0D0+poe)
            case (1) !Plane stress
                d(1,1)=Eme/(1.0D0-poe*poe)
                d(1,2)=poe*d(1,1)
                d(2,1)=d(1,2)
                d(2,2)=d(1,1)
                d(3,3)=Eme*0.5D0/(1.0D0+poe)
        end Select
    end subroutine DMAT_PL


    real(8) function CALA_PL3(ne,kakom,x,y,NELT,NODT)
        !*************************************************************
        !Area of 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        integer::i,j,k

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        CALA_PL3=0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))

    end function CALA_PL3


    subroutine BMAT_PL3(ne,kakom,x,y,b,NELT,NODT)
        !*************************************************************
        ![B] matrix for 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::b(1:3,1:6)

        integer::i,j,k
        real(8)::xi,xj,xk,yi,yj,yk,area

        !Strain-displacement matrix [B]
        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        xi=x(i)
        xj=x(j)
        xk=x(k)
        yi=y(i)
        yj=y(j)
        yk=y(k)

        do i=1,3
            do j=1,6
                b(i,j)=0.0D0
            end do
        end do
        b(1,1)=yj-yk
        b(1,3)=yk-yi
        b(1,5)=yi-yj
        b(2,2)=xk-xj
        b(2,4)=xi-xk
        b(2,6)=xj-xi
        b(3,1)=b(2,2)
        b(3,2)=b(1,1)
        b(3,3)=b(2,4)
        b(3,4)=b(1,3)
        b(3,5)=b(2,6)
        b(3,6)=b(1,5)
        do i=1,3
            do j=1,6
                b(i,j)=b(i,j)*0.5D0/area
            end do
        end do

    end subroutine BMAT_PL3


    subroutine AK_PL3(ne,kakom,x,y,NSTRES,Em,po,t,ak,NELT,NODT)
        !*************************************************************
        !Element stiffness matrix for 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,NSTRES
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT)
        real(8),intent(out)::ak(1:6,1:6)

        integer::i,j,k
        real(8)::b(1:3,1:6),d(1:3,1:3),str(1:3,1:6)
        real(8)::area,Eme,poe

        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)
        call BMAT_PL3(ne,kakom,x,y,b,NELT,NODT)
        Eme=Em(ne)
        poe=po(ne)
        call DMAT_PL(NSTRES,Eme,poe,d)
        ![sm]=A*t*[B]T[D][B]
        do i=1,3
            do j=1,6
                str(i,j)=0.0D0
                do k=1,3
                    str(i,j)=str(i,j)+d(i,k)*b(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                ak(i,j)=0.0D0
                do k=1,3
                    ak(i,j)=ak(i,j)+b(k,i)*str(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                ak(i,j)=ak(i,j)*area*t(ne)
            end do
        end do

    end subroutine AK_PL3


    subroutine AM_PL3(ne,kakom,x,y,t,gamma,am,NELT,NODT)
        !*************************************************************
        !Element mass matrix for 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::t(1:NELT),gamma(1:NELT)
        real(8),intent(out)::am(1:6,1:6)

        integer::i,j,k
        real(8)::mmat(1:6,1:6),w(1:6),area

        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)

        mmat(1,1)=0.50D0
        mmat(1,2)=0.00D0
        mmat(1,3)=0.25D0
        mmat(1,4)=0.00D0
        mmat(1,5)=0.25D0
        mmat(1,6)=0.00D0
        mmat(2,1)=0.00D0
        mmat(2,2)=0.50D0
        mmat(2,3)=0.00D0
        mmat(2,4)=0.25D0
        mmat(2,5)=0.00D0
        mmat(2,6)=0.25D0
        mmat(3,1)=0.25D0
        mmat(3,2)=0.00D0
        mmat(3,3)=0.50D0
        mmat(3,4)=0.00D0
        mmat(3,5)=0.25D0
        mmat(3,6)=0.00D0
        mmat(4,1)=0.00D0
        mmat(4,2)=0.25D0
        mmat(4,3)=0.00D0
        mmat(4,4)=0.50D0
        mmat(4,5)=0.00D0
        mmat(4,6)=0.25D0
        mmat(5,1)=0.25D0
        mmat(5,2)=0.00D0
        mmat(5,3)=0.25D0
        mmat(5,4)=0.00D0
        mmat(5,5)=0.50D0
        mmat(5,6)=0.00D0
        mmat(6,1)=0.00D0
        mmat(6,2)=0.25D0
        mmat(6,3)=0.00D0
        mmat(6,4)=0.25D0
        mmat(6,5)=0.00D0
        mmat(6,6)=0.50D0

        do i=1,6
            do j=1,6
                am(i,j)=mmat(i,j)*gamma(ne)*area*t(ne)/3.0D0/ggg
            end do
        end do

    end subroutine AM_PL3


    subroutine AC_PL3(ne,zetaM,zetaK,am,ak,ac,NELT)
        !*************************************************************
        !Element damping matrix for 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT
        real(8),intent(in)::zetaM(1:NELT),zetaK(1:NELT)
        real(8),intent(in)::am(1:6,1:6),ak(1:6,1:6)
        real(8),intent(out)::ac(1:6,1:6)

        integer::i,j

        !Initialization
        do i=1,6
            do j=1,6
                ac(i,j)=0.0
            end do
        end do
        !Element damping matrix
        do i=1,6
            do j=1,6
                ac(i,j)=zetaM(ne)*am(i,j)+zetaK(ne)*ak(i,j)
            end do
        end do
    end subroutine AC_PL3


    subroutine CALST_PL3(ne,kakom,x,y,dist,NSTRES,Em,po,t,strsave,NELT,NODT,nt)
        !*************************************************************
        !Stress calculation for 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,nt,NSTRES
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT),dist(1:nt)
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT)
        real(8),intent(out)::strsave(1:NELT,1:3)

        integer::i,j,ia,k1,k2,k3
        real(8)::b(1:3,1:6),d(1:3,1:3),epsilon(1:3)
        real(8)::wd(1:6),ws(1:3),area,Eme,poe

        !Nodal displacement of element
        do i=1,3
            ia=kakom(ne,i)
            wd(2*i-1)=dist(2*ia-1)
            wd(2*i)=dist(2*ia)
        end do
        !Stresses due to nodal displacement
        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)
        call BMAT_PL3(ne,kakom,x,y,b,NELT,NODT)
        Eme=Em(ne)
        poe=po(ne)
        call DMAT_PL(NSTRES,Eme,poe,d)
        do i=1,3
            epsilon(i)=0.0D0
            do j=1,6
                epsilon(i)=epsilon(i)+b(i,j)*wd(j)
            end do
        end do
        !{sigma}=[D]{epsilon}
        do i=1,3
            ws(i)=0.0D0
            do j=1,3
                ws(i)=ws(i)+d(i,j)*epsilon(j)
            end do
        end do
        !to save calculated stresses
        strsave(ne,1)=ws(1)
        strsave(ne,2)=ws(2)
        strsave(ne,3)=ws(3)

    end subroutine CALST_PL3



    subroutine BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
        !*************************************************************
        ![B] matrix for 4 node element
        !*************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(NELT,4)
        real(8),intent(in)::x(NODT),y(NODT),a,b
        real(8),intent(out)::bm(3,8),detJ

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

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

        !(dN/dadN/db)
        dn1a=-0.25D0*(1.0-b)
        dn2a= 0.25D0*(1.0-b)
        dn3a= 0.25D0*(1.0+b)
        dn4a=-0.25D0*(1.0+b)
        dn1b=-0.25D0*(1.0-a)
        dn2b=-0.25D0*(1.0+a)
        dn3b= 0.25D0*(1.0+a)
        dn4b= 0.25D0*(1.0-a)

        !Jacobi matrix & 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

        !Element of [B] matrix (dN/dx,dN/dy)
        do i=1,3
            do j=1,8
                bm(i,j)=0.0D0
            end do
        end do
        bm(1,1)= J22*dn1a-J12*dn1b
        bm(1,3)= J22*dn2a-J12*dn2b
        bm(1,5)= J22*dn3a-J12*dn3b
        bm(1,7)= J22*dn4a-J12*dn4b
        bm(2,2)=-J21*dn1a+J11*dn1b
        bm(2,4)=-J21*dn2a+J11*dn2b
        bm(2,6)=-J21*dn3a+J11*dn3b
        bm(2,8)=-J21*dn4a+J11*dn4b
        bm(3,1)=-J21*dn1a+J11*dn1b
        bm(3,2)= J22*dn1a-J12*dn1b
        bm(3,3)=-J21*dn2a+J11*dn2b
        bm(3,4)= J22*dn2a-J12*dn2b
        bm(3,5)=-J21*dn3a+J11*dn3b
        bm(3,6)= J22*dn3a-J12*dn3b
        bm(3,7)=-J21*dn4a+J11*dn4b
        bm(3,8)= J22*dn4a-J12*dn4b
        do i=1,3
            do j=1,8
                bm(i,j)=bm(i,j)/detJ
            end do
        end do
    end subroutine BMAT_PL4


    subroutine IPOINT4(kk,a,b)
        !*************************************************************
        !Gauss points
        !*************************************************************
        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 AK_PL4(ne,kakom,x,y,NSTRES,Em,po,t,ak,NELT,NODT)
        !*************************************************************
        !Element stiffness matrix for 4 node element
        !*************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        integer,intent(in)::NSTRES
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT)
        real(8),intent(out)::ak(1:8,1:8)

        integer :: kk,i,j,k
        real(8) :: d(3,3),bm(3,8),str(3,8)
        real(8) :: detJ,a,b

        !Initialization
        do i=1,8
            do j=1,8
                ak(i,j)=0.0D0
            end do
        end do
        !(sm)=(B)T(D,B)*t*detJ
        call DMAT_PL(NSTRES,Em(ne),po(ne),d)
        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
            do i=1,3
                do j=1,8
                    str(i,j)=0.0
                    do k=1,3
                        str(i,j)=str(i,j)+d(i,k)*bm(k,j)
                    end do
                end do
            end do
            do i=1,8
                do j=1,8
                    do k=1,3
                        ak(i,j)=ak(i,j)+bm(k,i)*str(k,j)*t(ne)*detJ
                    end do
                end do
            end do
        end do
    end subroutine AK_PL4


    subroutine AM_PL4(ne,kakom,x,y,t,gamma,am,NELT,NODT)
        !*************************************************************
        !Element mass matrix for 4 node element
        !*************************************************************
        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)::t(1:NELT),gamma(1:NELT)
        real(8),intent(out)::am(1:8,1:8)

        integer :: i,j,k,l,ii,jj,kk,ll
        real(8) :: nm(1:2,1:8),ntn(1:8,1:8)
        real(8) :: detJ,a,b,J11,J12,J21,J22
        real(8) :: dn1a,dn2a,dn3a,dn4a,dn1b,dn2b,dn3b,dn4b

        do ii=1,8
            do jj=1,8
                ntn(ii,jj)=0.0D0
            end do
        end do
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        l=kakom(ne,4)
        do kk=1,4
            call IPOINT4(kk,a,b)
            ![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
            do ii=1,2
                do jj=1,8
                    nm(ii,jj)=0.0D0
                end do
            end do
            nm(1,1)=0.25D0*(1.0D0-a)*(1.0D0-b)
            nm(1,3)=0.25D0*(1.0D0+a)*(1.0D0-b)
            nm(1,5)=0.25D0*(1.0D0+a)*(1.0D0+b)
            nm(1,7)=0.25D0*(1.0D0-a)*(1.0D0+b)
            nm(2,2)=nm(1,1)
            nm(2,4)=nm(1,3)
            nm(2,6)=nm(1,5)
            nm(2,8)=nm(1,7)
            do ii=1,8
                do jj=1,8
                    do ll=1,2
                        ntn(ii,jj)=ntn(ii,jj)+nm(ll,ii)*nm(ll,jj)
                    end do
                end do
            end do
        end do
        do i=1,8
            do j=1,8
                am(i,j)=gamma(ne)*t(ne)/ggg*ntn(i,j)*detJ
            end do
        end do
    end subroutine AM_PL4


    subroutine AC_PL4(ne,zetaM,zetaK,am,ak,ac,NELT)
        !*************************************************************
        !Element damping matrix for 4 node element
        !*************************************************************
        integer,intent(in)::ne,NELT
        real(8),intent(in)::zetaM(1:NELT),zetaK(1:NELT)
        real(8),intent(in)::am(1:8,1:8),ak(1:8,1:8)
        real(8),intent(out)::ac(1:8,1:8)

        integer::i,j

        !Initialization
        do i=1,8
            do j=1,8
                ac(i,j)=0.0
            end do
        end do
        !damping matrix
        do i=1,8
            do j=1,8
                ac(i,j)=zetaM(ne)*am(i,j)+zetaK(ne)*ak(i,j)
            end do
        end do
    end subroutine AC_PL4


    subroutine CALST_PL4(ne,kakom,x,y,dist,NSTRES,Em,po,t,strsave,NELT,NODT,nt)
        !*************************************************************
        !Stress calculation for 4 node element
        !*************************************************************
        integer,intent(in)::NELT,NODT,nt
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT),dist(1:nt)
        integer,intent(in)::NSTRES
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT)
        real(8),intent(out)::strsave(1:NELT,1:3)

        integer :: kk,i,j,ia
        integer :: k1,k2,k3,k4
        real(8) :: detJ,a,b
        real(8) :: d(1:3,1:3),bm(1:3,1:8),epsilon(1:3),stress(1:3),wd(1:8)
        real(8) :: Eme,poe
        real(8) :: sigx,sigy,tauxy

        !Displacement
        do i=1,4
            ia=kakom(ne,i)
            wd(2*i-1)=dist(2*ia-1)
            wd(2*i)=dist(2*ia)
        end do
        k1=kakom(ne,1)
        k2=kakom(ne,2)
        k3=kakom(ne,3)
        k4=kakom(ne,4)

        !Stress in element
        Eme=Em(ne)
        poe=po(ne)
        call DMAT_PL(NSTRES,Eme,poe,d)

        sigx=0.0D0
        sigy=0.0D0
        tauxy=0.0D0
        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
            !{epsilon}=(B){u} 
            do i=1,3
                epsilon(i)=0.0D0
                do j=1,8
                    epsilon(i)=epsilon(i)+bm(i,j)*wd(j)
                end do
            end do
            !{stress}=(D){epsilon} 
            do i=1,3
                stress(i)=0.0D0
                do j=1,3
                    stress(i)=stress(i)+d(i,j)*epsilon(j)
                end do
            end do
            sigx=sigx+stress(1)
            sigy=sigy+stress(2)
            tauxy=tauxy+stress(3)
        end do
        strsave(ne,1)=0.25D0*sigx
        strsave(ne,2)=0.25D0*sigy
        strsave(ne,3)=0.25D0*tauxy

    end subroutine CALST_PL4


    subroutine PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
        !*************************************************************
        !Principal stress calculation
        !*************************************************************
        real(8),intent(in)::sigx,sigy,tauxy
        real(8),intent(out)::ps1,ps2,ang

        ps1=0.5D0*(sigx+sigy)+sqrt(0.25D0*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        ps2=0.5D0*(sigx+sigy)-sqrt(0.25D0*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        if(sigx==sigy)then
            if(tauxy>0.0D0)ang=45.0D0
            if(tauxy<0.0D0)ang=135.0D0
            if(tauxy==0.0D0)ang=0.0D0
        else
            ang=0.5D0*atan(2.0D0*tauxy/(sigx-sigy))
            ang=180.0D0/pi*ang
            if(sigx>sigy.and.tauxy>=0.0D0)ang=ang
            if(sigx>sigy.and.tauxy<0.0D0)ang=ang+180.0D0
            if(sigx<sigy)ang=ang+90.0D0
        end if
    end subroutine PST_PL


    subroutine ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !*************************************************************
        !Assembly of vector
        !*************************************************************
        integer :: NELT,NODT
        integer :: ne,nod,nhen
        integer :: kakom(1:NELT,1:nod)
        real(8) :: fevec(1:nod*nhen)
        real(8) :: 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
        !*************************************************************
        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

        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,mn
        real(8) :: work(1:nt)

        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_V2DIM