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


program f90_FEM_VFRM
    !*****************************************************************
    !2D Frame 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::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::KOZ                        !Number of restricted nodes in rotation
    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::Em(:)          !Elasyic mudulus of element
    real(8),allocatable::AA(:)          !Section area of element
    real(8),allocatable::AI(:)          !Moment of inertia 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 of node
    integer,allocatable::nokx(:)        !Restricted node number in x-direction
    integer,allocatable::noky(:)        !Restricted node number in y-direction
    integer,allocatable::nokz(:)        !Restricted node number in rotation
    real(8),allocatable::f(:)           !External load 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::hen(:,:)       !Transform 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::sresul(:,:)    !Resultant force of element
    integer::ndata                      !Number of time history data
    real(8)::dt                         !Time increment
    real(8),allocatable::accinp(:)      !Acceleration time histroy data
    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::prsresul(:,:)  !Memory of stress resultant for output
    real(8)::dismax                     !work for search of Max.displacement
    integer::kpriii                     !Memory for Max.disp step

    real(8),allocatable::ftvec(:)       !Total 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=2                      !Number of nodes per element
    integer::nhen=3                     !Number of degrees 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 time history input data file name
    integer::numnd,numsr                !Number of displacement degree and stress resultant 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 resultant degree number for output of time history data

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

    fmt1="(i5,2(',',e15.7),3(',',e15.7),3(',',i2))"    !Node
    fmt2="(i5,2(',',i5),6(',',e15.7),',',i5)"          !Element
    fmt3="(i5,2(',',e15.7),3(',',e15.7))"              !Displacement, reaction
    fmt4="(i5,3(',',e15.7),3(',',e15.7))"              !Stress resultant

    call getarg(1,fnameR) !Input data file (structure data)
    call getarg(2,fnameW) !Output data file
    call getarg(3,fnameA) !Acceleration time history data input file

    !To start the measurement of calculation time
    t1=DIFT()
    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    !/*--------------------------------------------------------------------------------------*/
    !/*Comment*/
    !/*--------------------------------------------------------------------------------------*/
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
    !Basic data for analysis
    !------------------------------------------------------------------------------------------
    read(11,*) NODT,NELT,MATEL,KOX,KOY,KOZ,NF
    !------------------------------------------------------------------------------------------
    !To declare array size
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:MATEL,1:6))
    allocate(Em(1:NELT))
    allocate(AA(1:NELT))
    allocate(AI(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))
    if(1<=KOZ)allocate(nokz(1:KOZ))
    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(hen(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(sresul(1:NELT,1:nod*nhen))
    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(prsresul(1:NELT,1:nod*nhen))
    !------------------------------------------------------------------------------------------
    !Material properties (Em,AA,AI,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 (node-1, node-2, matno)
    !------------------------------------------------------------------------------------------
    do ne=1,NELT
        read(11,*) kakom(ne,1),kakom(ne,2),matno(ne)
        do i=1,MATEL
            if(i==matno(ne))then
                Em(ne)    =wm(i,1)
                AA(ne)   =wm(i,2)
                AI(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
    !------------------------------------------------------------------------------------------
    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
    if(0<KOZ)then
        do i=1,KOZ
            read(11,*) nokz(i)
        end do
    end if
    !------------------------------------------------------------------------------------------
    !Loaded node number, load in x-direction, load in y-direction, moment
    !------------------------------------------------------------------------------------------
    do i=1,nt
        f(i)=0.0D0
    end do
    if(0<NF)then
        do i=1,NF
            read(11,*) n,f(3*n-2),f(3*n-1),f(3*n)
        end do
    end if
    !------------------------------------------------------------------------------------------
    !Information of node 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 stress resultant
    !------------------------------------------------------------------------------------------
    read(11,*) numsr
    allocate(nsrnum(1:numsr,1:2))
    do i=1,numsr
        read(11,*) (nsrnum(i,j),j=1,2) !Element number, stress resultant degree number
    end do
    !------------------------------------------------------------------------------------------
    !Node number & displacement degree number for output of all displacements and all stress resultants
    !------------------------------------------------------------------------------------------
    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)') 'NODT,NELT,MATEL,KOX,KOY,KOZ,NF,dt,ndata'
        write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',KOZ,',',NF,',',dt,',',ndata
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,'(a)') '*node characteristics'
    write(12,'(a)') 'node,x,y,fx,fy,fz,fix-x,fix-y,fix-z'
    do i=1,NODT
        k=0
        l=0
        m=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
        if(0<KOZ)then
            do j=1,KOZ
                if(i==nokz(j))m=1
            end do
        end if
            write(linebuf,FMT=fmt1) i,x(i),y(i),f(3*i-2),f(3*i-1),f(3*i),&
                k,l,m
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(12,'(a)') '*element characteristics'
    write(12,'(a)') 'element,node-1,node-2,E,A,I,gamma,zetaM,zetaK,matno'
    do ne=1,NELT
            write(linebuf,FMT=fmt2) ne,kakom(ne,1),kakom(ne,2),Em(ne),AA(ne),AI(ne),gamma(ne),zetaM(ne),zetaK(ne),matno(ne)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    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 matrices
    !-----------------------------------------------------
    do ne=1,NELT
        !Stiffness matrix
        call AK_FRM(ne,kakom,x,y,Em,AA,AI,ak,hen,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,ak,tak,NELT,NODT)
        !Mass matrix
        call AM_FRM(ne,kakom,x,y,AA,gamma,am,hen,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,am,tam,NELT,NODT)
        !damping matrix
        call AC_FRM(ne,zetaM,zetaK,am,ak,ac,NELT)
        call ASMAT(ne,nod,nhen,kakom,ac,tac,NELT,NODT)
    end do
  !Total 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=3*nokx(i)-2
            mindex(n)=0
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            n=3*noky(i)-1
            mindex(n)=0
        end do
    end if
    if(0<KOZ)then
        do i=1,KOZ
            n=3*nokz(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 & to 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)!triangular 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 vector
        !-----------------------------------------------------
        !(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
        !-----------------------------------------------------
        !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 stresses
        !-----------------------------------------------------
        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
        do ne=1,NELT
            call AK_FRM(ne,kakom,x,y,Em,AA,AI,ak,hen,NELT,NODT)
            call CALST_FRM(ne,kakom,dis2,ak,hen,reac,sresul,Em,AA,NELT,NODT,nt)
        end do
        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=3*(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
            write(dummy,'(",",e15.7)') sresul(nsrnum(i,1),nsrnum(i,2))
            linebuf=trim(linebuf)//dummy
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
        !To update memory for output
        if(dismax<abs(dis2(3*(kprnod-1)+kprfre)))then
            dismax=abs(dis2(3*(kprnod-1)+kprfre))
            kpriii=iii
            do i=1,nt
                prdis(i)=dis2(i)
            end do
            do ne=1,NELT
                do j=1,6
                    prsresul(ne,j)=sresul(ne,j)
                end do
            end do
        end if
    !*****************************************************
    end do
    !*****************************************************
    close(12)

    !*****************************************************
    !To print all displacements and all stress resultants at defined 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,dis-z'
    do i=1,NODT
        write(linebuf,FMT=fmt3) i,x(i),y(i),prdis(3*i-2),prdis(3*i-1),prdis(3*i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !Stress resultant
    write(12,'(a)') '*stress resultants'
    write(12,'(a)') 'element,Ni,Si,Mi,Nj,Sj,Mj'
    do ne=1,NELT
        write(linebuf,FMT=fmt4) ne,prsresul(ne,1),prsresul(ne,2),prsresul(ne,3),&
                                   prsresul(ne,4),prsresul(ne,5),prsresul(ne,6)
        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 AK_FRM(ne,kakom,x,y,Em,AA,AI,ak,hen,NELT,NODT)
        !*************************************************************
        !Element stiffness matrix
        !*************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:2)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),AI(1:NELT)
        real(8),intent(out)::ak(1:6,1:6),hen(1:6,1:6)
        integer :: i,j
        real(8) :: AL,EA,EI,cs,sn,w(1:6,1:6)

        !Initialization
        do i=1,6
            do j=1,6
                hen(i,j)=0.0D0
                ak(i,j)=0.0D0
            end do
        end do
        !Transform matrix
        i=kakom(ne,1)
        j=kakom(ne,2)
        AL=sqrt((x(j)-x(i))**2.0D0+(y(j)-y(i))**2.0D0)
        cs=(x(j)-x(i))/AL
        sn=(y(j)-y(i))/AL
        hen(1,1)= cs
        hen(1,2)= sn
        hen(2,1)=-sn
        hen(2,2)= cs
        hen(3,3)=1.0D0
        do i=1,3
            do j=1,3
                hen(i+3,j+3)=hen(i,j)
            end do
        end do
        !Element stiffness matrix
        EA=Em(ne)*AA(ne)
        EI=Em(ne)*AI(ne)
        ak(1,1)=EA/AL
        ak(2,1)=0.0D0
        ak(3,1)=0.0D0
        ak(4,1)=-EA/AL
        ak(5,1)=0.0D0
        ak(6,1)=0.0D0
        ak(2,2)=12.0D0*EI/AL/AL/AL
        ak(3,2)=6.0D0*EI/AL/AL
        ak(4,2)=0.0D0
        ak(5,2)=-12.0D0*EI/AL/AL/AL
        ak(6,2)=6.0D0*EI/AL/AL
        ak(3,3)=4.0D0*EI/AL
        ak(4,3)=0.0D0
        ak(5,3)=-6.0D0*EI/AL/AL
        ak(6,3)=2.0D0*EI/AL
        ak(4,4)=EA/AL
        ak(5,4)=0.0D0
        ak(6,4)=0.0D0
        ak(5,5)=12.0D0*EI/AL/AL/AL
        ak(6,5)=-6.0D0*EI/AL/AL
        ak(6,6)=4.0D0*EI/AL
        do i=1,5
            do j=i+1,6
                ak(i,j)=ak(j,i)
            end do
        end do
        !To transform local coordinates into global coordinates
        do i=1,6
            do j=1,6
                w(i,j)=0.0D0
                do k=1,6
                    w(i,j)=w(i,j)+hen(k,i)*ak(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                ak(i,j)=0.0D0
                do k=1,6
                    ak(i,j)=ak(i,j)+w(i,k)*hen(k,j)
                end do
            end do
        end do

    end subroutine AK_FRM


    subroutine AM_FRM(ne,kakom,x,y,AA,gamma,am,hen,NELT,NODT)
        !*************************************************************
        !Element mass matrix
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:2)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::AA(1:NELT),gamma(1:NELT)
        real(8),intent(in)::hen(1:6,1:6)
        real(8),intent(out)::am(1:6,1:6)

        integer::i,j
        real(8)::AL,w(1:6,1:6)

        !Initialization
        do i=1,6
            do j=1,6
                am(i,j)=0.0
            end do
        end do
        !Element length
        i=kakom(ne,1)
        j=kakom(ne,2)
        AL=sqrt((x(j)-x(i))**2.0D0+(y(j)-y(i))**2.0D0)
        !Element mass matrix
        am(1,1)=1.0D0/3.0D0
        am(2,1)=0.0D0
        am(3,1)=0.0D0
        am(4,1)=1.0D0/6.0D0
        am(5,1)=0.0D0
        am(6,1)=0.0D0
        am(2,2)=13.0D0/35.0D0
        am(3,2)=11.0D0*AL/210.0D0
        am(4,2)=0.0D0
        am(5,2)=9.0D0/70.0D0
        am(6,2)=-13.0D0*AL/420.0D0
        am(3,3)=AL*AL/105.0D0
        am(4,3)=0.0D0
        am(5,3)=13.0D0*AL/420.0D0
        am(6,3)=-AL*AL/140.0D0
        am(4,4)=1.0D0/3.0D0
        am(5,4)=0.0D0
        am(6,4)=0.0D0
        am(5,5)=13.0D0/35.0D0
        am(6,5)=-11.0D0*AL/210.0D0
        am(6,6)=AL*AL/105.0D0
        do i=1,5
            do j=i+1,6
                am(i,j)=am(j,i)
            end do
        end do
        do i=1,6
            do j=1,6
                am(i,j)=gamma(ne)*AA(ne)*AL*am(i,j)/ggg
            end do
        end do
        !To transform local coordinates into global coordinates
        do i=1,6
            do j=1,6
                w(i,j)=0.0D0
                do k=1,6
                    w(i,j)=w(i,j)+hen(k,i)*am(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                am(i,j)=0.0D0
                do k=1,6
                    am(i,j)=am(i,j)+w(i,k)*hen(k,j)
                end do
            end do
        end do

    end subroutine AM_FRM


    subroutine AC_FRM(ne,zetaM,zetaK,am,ak,ac,NELT)
        !*************************************************************
        !Element damping matrix
        !*************************************************************
        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_FRM


    subroutine CALST_FRM(ne,kakom,dis,ak,hen,reac,sresul,Em,AA,NELT,NODT,nt)
        !*************************************************************
        !Stress resultant of element
        !*************************************************************
        integer,intent(in)::NELT,NODT,nt
        integer,intent(in)::ne,kakom(1:NELT,1:2)
        real(8),intent(in)::dis(1:nt),ak(1:6,1:6),hen(1:6,1:6)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT)
        real(8),intent(inout)::reac(1:nt)
        real(8),intent(out)::sresul(1:NELT,1:6)
        integer :: i,j,k,ia,ir,iw
        real(8) :: wd(1:6),ws(1:6),tempe

        !Displacement in global coordinate system
        do i=1,2
            ia=kakom(ne,i)
            wd(3*i-2)=dis(3*ia-2)
            wd(3*i-1)=dis(3*ia-1)
            wd(3*i)=dis(3*ia)
        end do
        !Nodal forces in global coordinate system
        do i=1,6
            ws(i)=0.0D0
            do k=1,6
                ws(i)=ws(i)+ak(i,k)*wd(k)
            end do
        end do
        !Stress resultant in local coordinate system
        !(To save stress resultant in array '[tsm]')
        do i=1,6
            sresul(ne,i)=0.0D0
            do k=1,6
                sresul(ne,i)=sresul(ne,i)+hen(i,k)*ws(k)
            end do
        end do
        !Stress resultant in global coordinate system
        do i=1,6
            ws(i)=0.0D0
            do j=1,6
                ws(i)=ws(i)+hen(j,i)*sresul(ne,j)
            end do
        end do
        !To set reaction
        do i=1,2
            ia=kakom(ne,i)
            do j=1,3
                ir=3*(ia-1)+j
                iw=3*(i-1)+j
                reac(ir)=reac(ir)+ws(iw)
            end do
        end do
    end subroutine CALST_FRM




    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_VFRM