program f90_fem_gfrm
    !***************************************************************************
    !2D Frame Geometrically Nonlineqr Analysis using arc-length increment method
    !***************************************************************************
    implicit none
    integer::i,j,k,l,m,kk,ne,n,ia,ja
    character::strcom*50                !Comment
    integer::NODT                       !Number of nodes
    integer::NELT                       !Number of elements
    integer::MATEL                      !Number of material sets
    integer::KOX                        !Number of fixed nodes in x-direction
    integer::KOY                        !Number of fixed nodes in y-direction
    integer::KOZ                        !Number of fixed 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(:)          !Elastic modulus of element
    real(8),allocatable::AA(:)          !Section area of element
    real(8),allocatable::AI(:)          !Moment of inertia of element
    real(8),allocatable::x(:)           !x-coordinate of node
    real(8),allocatable::y(:)           !y-coordinate of node
    integer,allocatable::nokx(:)        !Fixed node number in x-direction
    integer,allocatable::noky(:)        !Fixed node number in y-direction
    integer,allocatable::nokz(:)        !Fixed node number in rotation
    real(8),allocatable::df(:)          !External force increment
    real(8),allocatable::ff(:)          !Accumulated external force
    integer,allocatable::mindex(:)      !Index for treatment of fixed nodes
    real(8),allocatable::ftvec(:)       !Unbalanced force
    real(8),allocatable::tsm(:,:)       !Total stiffness matrix
    real(8),allocatable::dis(:)         !Displacement increment
    real(8),allocatable::dist(:)        !Accumulated displacement
    real(8),allocatable::reac(:)        !work array
    real(8),allocatable::sresul(:,:)    !Stress resultant
    real(8),allocatable::sm(:,:)        !Element stiffness matrix
    real(8),allocatable::hen(:,:)       !transform matrix
    real(8),allocatable::fevec(:)       !Element load vector
    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::nload,iteration            !Counter of loading and iterative calculation
    integer::nlmax                      !Limit of loading step
    integer::itmax                      !Limit of iteration step
    real(8)::sumd1,sumd2,sumf1,sumf2    !work for judgement of convergence

    !Variables for arc-length method
    real(8),allocatable::dis0(:)        !Displacement increment due to base load increment
    real(8),allocatable::disR(:)        !Displacement increment due to unbalanced force
    real(8),allocatable::ddis(:)        !Displacement increment in iterative step
    real(8)::deltaS0,deltaS             !Arc-length
    real(8)::coefA,Spara                !Scaling parameter
    real(8)::lam,dlam,dlam0             !Load control arameter
    real(8)::sum1,sum2,sum3,sum4        !work for calculation
    real(8)::gamma0

    character::fnameR*50,fnameW*50      !Input & output file name

    character :: linebuf*1000,dummy*100
    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),3(',',e15.7),',',i5)"                        !Element
    fmt3="(i5,2(',',e15.7),3(',',e15.7),3(',',e15.7),3(',',e15.7))"  !Displacement & reaction
    fmt4="(i5,3(',',e15.7),3(',',e15.7))"                            !Stress resultant

    call getarg(1,fnameR) !Input file name
    call getarg(2,fnameW) !Output file name
    call getarg(3,dummy)  !nlmax
    read(dummy,*) nlmax   !To convert character 'dummy' to integer 'nlmax'
    call getarg(4,dummy)  !itmax
    read(dummy,*) itmax   !To convert character 'dummy' to integer 'itmax'
    call getarg(5,dummy)  !coefA
    read(dummy,*) coefA   !To convert character 'dummy' to real number of double presision 'coefA'
    call getarg(6,dummy)
    read(dummy,*) gamma0

    !To start the measurement of calculation time
    t1=DIFT()
    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    !/*--------------------------------------------------------------------------------------*/
    !/*Comment*/
    !/*--------------------------------------------------------------------------------------*/
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
    !Basic values 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:NELT,1:3))
    allocate(Em(1:NELT))
    allocate(AA(1:NELT))
    allocate(AI(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(df(1:nt))
    allocate(ff(1:nt))
    allocate(mindex(1:nt))
    allocate(ftvec(1:nt))
    allocate(tsm(1:nt,1:nt))
    allocate(dis(1:nt))
    allocate(dist(1:nt))
    allocate(reac(1:nt))
    allocate(sresul(1:NELT,1:nod*nhen))
    allocate(sm(1:nod*nhen,1:nod*nhen))
    allocate(hen(1:nod*nhen,1:nod*nhen))
    allocate(fevec(1:nod*nhen))
    allocate(dis0(1:nt))
    allocate(disR(1:nt))
    allocate(ddis(1:nt))
    !------------------------------------------------------------------------------------------
    !Material properties (Em,AA,AI)
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3)
    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)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !(x,y) coordinates of nodes
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),y(i)
    end do
    !------------------------------------------------------------------------------------------
    !Fixed 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 increment in x-dir., load increment in y-dir., moment increment
    !------------------------------------------------------------------------------------------
    do i=1,nt
        df(i)=0.0D0
    end do
    if(0<NF)then
        do i=1,NF
            read(11,*) n,df(3*n-2),df(3*n-1),df(3*n)
        end do
    end if
    close(11)

    !*****************************************************
    !To print input data
    !*****************************************************
    open(12, file=fnameW, status='replace')
    write(12,'(a)') strcom
    write(12,'(a)') 'NODT,NELT,MATEL,KOX,KOY,KOZ,NF,nlmax,itmax'
        write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',KOZ,',',NF,',',nlmax,',',itmax
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,'(a)') '*node characteristics'
    write(12,'(a)') 'node,x,y,dfx,dfy,dfz,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),df(3*i-2),df(3*i-1),df(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,matno'
    do ne=1,NELT
            write(linebuf,FMT=fmt2) ne,kakom(ne,1),kakom(ne,2),Em(ne),AA(ne),AI(ne),matno(ne)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do

    !***************************
    !Scaling parameter
    !***************************
    sum2=0.0D0
    do i=1,nt
        sum2=sum2+df(i)*df(i)
    end do
    Spara=sqrt(coefA/sum2)
    !***************************
    !Initialization
    !***************************
    do i=1,nt
        ff(i)=0.0D0
        dis(i)=0.0D0
        dist(i)=0.0D0
        ddis(i)=0.0D0
    end do
    do ne=1,NELT
        do i=1,6
            sresul(ne,i)=0.0D0
        end do
    end do
    !***************************
    !Treatment of fixed 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 nload=1,nlmax    !Loading
    !=====================================================
            !-----------------------------------------------------
            !Initialization
            !-----------------------------------------------------
            lam=0.0D0
            do i=1,nt
                reac(i)=0.0D0
                do j=1,nt
                    tsm(i,j)=0.0D0
                end do
            end do
            !-----------------------------------------------------
            !To assemble matrix
            !-----------------------------------------------------
            do ne=1,NELT
                !Stiffness matrix
                call FRAMESM(ne,kakom,matno,x,y,dist,Em,AA,AI,sresul,sm,NELT,NODT,MATEL,nt,gamma0)
                call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
            end do
            !-----------------------------------------------------
            !Treatment of fixed nodes
            !-----------------------------------------------------
            do i=1,mm
                ia=mindex(i)
                do j=1,mm
                    ja=mindex(j)
                    tsm(i,j)=tsm(ia,ja)
                end do
            end do
            !-----------------------------------------------------
            !Inverse matrix
            !-----------------------------------------------------
            call MATINV(mm,tsm,nt)
            !dis0
            do i=1,nt
                reac(i)=df(i)
            end do
            call CALDIS(mm,mindex,reac,tsm,nt)
            do i=1,nt
                dis0(i)=reac(i)
            end do
            !--------------------------------------------------------------------------------------
            ! Treatment of 'dlam' and 'deltaS'
            ! + The sign of 'dlam' is defined using inner product of 'ddis' and 'dis0.'
            !   If inner product is positive, the angle is less than or equal to 90 degree.
            !   In this case, the sign of 'dlam' is defined as positive.
            ! + If load-displacement curve becomes horizontal, 'dlam' becomes nearly equal zero.
            !   It means displacement and load increments are zero.
            !   Therefore, when 'dlam' becomes less than or equal to 1/10 times initial value,
            !   arc-length 'deltaS' shall be increased.
            !--------------------------------------------------------------------------------------
            sum1=0.0D0
            sum2=0.0D0
            sum4=0.0D0
            do i=1,nt
                sum1=sum1+dis0(i)*dis0(i)
                sum2=sum2+df(i)*df(i)
                sum4=sum4+ddis(i)*dis0(i)
            end do
            if(nload==1)then
                deltaS0=sqrt(sum1+Spara*Spara*sum2)
                deltaS=deltaS0
                dlam0=sqrt(deltaS*deltaS/(sum1+Spara*Spara*sum2))
            end if
            dlam=sqrt(deltaS*deltaS/(sum1+Spara*Spara*sum2))
            if(abs(dlam)<0.1D0*abs(dlam0))deltaS=deltaS*1.2D0
            if(1.0D0*abs(dlam0)<abs(dlam))deltaS=deltaS
            dlam=sqrt(deltaS*deltaS/(sum1+Spara*Spara*sum2))
            if(sum4<0.0D0)dlam=-dlam
            do i=1,nt
                dis(i)=dlam*dis0(i)
            end do
            do i=1,nt
                ddis(i)=0.0D0 !Displacement increment in irterative step
            end do
        !=====================================================
        do iteration=1,itmax    !Iterative process
        !=====================================================
            lam=lam+dlam
            do i=1,nt
                dist(i)=dist(i)+dis(i)
                ddis(i)=ddis(i)+dis(i)
            end do
            !-----------------------------------------------------
            !calculation of internal stress
            !-----------------------------------------------------
            do i=1,nt
                reac(i)=0.0D0
            end do
            do ne=1,NELT
                call CALSTFRM(ne,kakom,matno,x,y,dis,dist,Em,AA,AI,fevec,sresul,NELT,NODT,MATEL,nt,gamma0)
                call ASVEC(ne,nod,nhen,kakom,fevec,reac,NELT,NODT)
            end do
            !-----------------------------------------------------
            !calculation of unbalanced force
            !-----------------------------------------------------
            sumd1=0.0D0
            sumd2=0.0D0
            sumf1=0.0D0
            sumf2=0.0D0
            do i=1,nt
                ftvec(i)=ff(i)+lam*df(i)-reac(i)
                sumd1=sumd1+abs(dis(i))
                sumd2=sumd2+abs(dist(i))
            end do
            do i=1,mm
                n=mindex(i)
                sumf1=sumf1+abs(ftvec(n))
                sumf2=sumf2+abs(reac(n))
            end do
            !-----------------------------------------------------
            !Judgement of convergence
            !-----------------------------------------------------
            print "('nload=',i0,' iteration=',i0,' ratio=',E15.7,' lam=',E15.7)",nload,iteration,sumf1/sumf2,lam
            if(sumd1/sumd2<1.0D-6.and.sumf1/sumf2<1.0D-3)exit
            !-----------------------------------------------------
            !initialization
            !-----------------------------------------------------
            do i=1,nt
                do j=1,nt
                    tsm(i,j)=0.0D0
                end do
            end do
            !-----------------------------------------------------
            !To assemble matrix
            !-----------------------------------------------------
            do ne=1,NELT
                !Stiffness matrix
                call FRAMESM(ne,kakom,matno,x,y,dist,Em,AA,AI,sresul,sm,NELT,NODT,MATEL,nt,gamma0)
                call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
            end do
            !-----------------------------------------------------
            !Treatment of fixed nodes
            !-----------------------------------------------------
            do i=1,mm
                ia=mindex(i)
                do j=1,mm
                    ja=mindex(j)
                    tsm(i,j)=tsm(ia,ja)
                end do
            end do
            !-----------------------------------------------------
            !Inverse matrix
            !-----------------------------------------------------
            call MATINV(mm,tsm,nt)
            !dis0
            do i=1,nt
                reac(i)=df(i)
            end do
            call CALDIS(mm,mindex,reac,tsm,nt)
            do i=1,nt
                dis0(i)=reac(i)
            end do
            !disR
            do i=1,nt
                reac(i)=ftvec(i)
            end do
            call CALDIS(mm,mindex,reac,tsm,nt)
            do i=1,nt
                disR(i)=reac(i)
            end do
            !dlam
            sum1=0.0D0
            sum2=0.0D0
            sum3=0.0D0
            do i=1,nt
                sum1=sum1+dis0(i)*dis0(i)
                sum2=sum2+df(i)*df(i)
                sum3=sum3+dis0(i)*disR(i)
            end do
            dlam=-sum3/(sum1+Spara*Spara*sum2)
            do i=1,nt
                dis(i)=dlam*dis0(i)+disR(i)
            end do
        !=====================================================
        end do    !Loop of iterative process
        !=====================================================
        !*****************************************************
        !To print the result
        !*****************************************************
        !To define load value
        do i=1,nt
            ff(i)=ff(i)+lam*df(i)
        end do
        write(12,'(a)') '*** New load step ***'
        write(12,'(a)') 'nload,iteration,lam,deltaS,coefA'
        write(linebuf,'(i0,",",i0,",",e15.7,",",e15.7,",",e15.7)') nload,iteration,lam,deltaS,coefA
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
        write(12,'(a)') '*displacements and forces'
        write(12,'(a)') 'node,coord-x,coord-y,dist-x,dist-y,dist-z,ff-x,ff-y,ff-z,ftvec-x,ftvec-y,ftvec-z'
        do i=1,NODT
            write(linebuf,FMT=fmt3) i,x(i),y(i),dist(3*i-2),dist(3*i-1),dist(3*i),&
                ff(3*i-2),ff(3*i-1),ff(3*i),ftvec(3*i-2),ftvec(3*i-1),ftvec(3*i)
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
        write(12,'(a)') '*stress resultants'
        write(12,'(a)') 'element,Ni,Si,Mi,Nj,Sj,Mj'
        do ne=1,NELT
            write(linebuf,FMT=fmt4) ne,sresul(ne,1),sresul(ne,2),sresul(ne,3),sresul(ne,4),sresul(ne,5),sresul(ne,6)
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
        if(itmax<=iteration)exit
    !=====================================================
    end do    !Loop of loading
    !=====================================================

    !To finish the measurement of calculation time
    t2=DIFT()
    write(12,'("NODT=",i0," nt=",i0," mm=",i0)') ,NODT,nt,mm
    write(12,'("Calculation time=",f10.3,"(sec)")') t2-t1
    write(linebuf,*) 'Date_time='//trim(CALLDATE())//' (gfortran)'
    write (12,'(a)') trim(adjustl(linebuf))
    close(12)

stop

contains

    subroutine CALDIS(mm,mindex,reac,tsm,nt)
        !*************************************************************
        !Calculation of displacement using inverse matrix
        !*************************************************************
        integer,intent(in)::mm,nt
        integer,intent(in)::mindex(1:nt)
        real(8),intent(in)::tsm(1:nt,1:nt)
        real(8),intent(inout)::reac(1:nt)
        integer::i
        real(8)::work(1:nt)

        do i=1,mm
            ia=mindex(i)
            reac(i)=reac(ia)
        end do
        do i=1,mm
            work(i)=0.0D0
            do j=1,mm
                work(i)=work(i)+tsm(i,j)*reac(j)
            end do
        end do
        do i=1,nt
            reac(i)=0.0D0
        end do
        do i=1,mm
            ia=mindex(i)
            reac(ia)=work(i)
        end do
    end subroutine CALDIS


    subroutine FRAMESM(ne,kakom,matno,x,y,dist,Em,AA,AI,sresul,sm,NELT,NODT,MATEL,nt,gamma0)
        !*************************************************************
        !Element stiffness matrix (linear part + geometric stiffness)
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,MATEL,nt
        integer,intent(in)::kakom(1:NELT,1:2),matno(1:MATEL)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::dist(1:nt)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),AI(1:NELT)
        real(8),intent(in)::sresul(1:NELT,1:6)
        real(8),intent(out)::sm(1:6,1:6)
        real(8),intent(in)::gamma0

        integer::i,j,k
        real(8)::hen(1:6,1:6),sme(1:6,1:6),smg(1:6,1:6),w(1:6,1:6)
        real(8)::AL,P,xx1,yy1,xx2,yy2

        do i=1,6
            do j=1,6
                smg(i,j)=0.0D0
                sm(i,j)=0.0D0
            end do
        end do
        !Axial force
        P=0.50D0*(sresul(ne,4)-sresul(ne,1))
        !Geometric stiffness
        i=kakom(ne,1)
        j=kakom(ne,2)
        xx1=x(i)+dist(3*i-2)
        yy1=y(i)+dist(3*i-1)
        xx2=x(j)+dist(3*j-2)
        yy2=y(j)+dist(3*j-1)
        AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        smg(1,1)=0.0D0
        smg(2,1)=0.0D0
        smg(2,2)=6.0D0/5.0D0/AL
        smg(3,1)=0.0D0
        smg(3,2)=1.0D0/10.0D0
        smg(3,3)=2.0D0*AL/15.0D0
        smg(4,1)=0.0D0
        smg(4,2)=0.0D0
        smg(4,3)=0.0D0
        smg(4,4)=0.0D0
        smg(5,1)=0.0D0
        smg(5,2)=-6.0D0/5.0D0/AL
        smg(5,3)=-1.0D0/10.0D0
        smg(5,4)=0.0D0
        smg(5,5)=6.0D0/5.0D0/AL
        smg(6,1)=0.0D0
        smg(6,2)=1.0D0/10.0D0
        smg(6,3)=-AL/30.0D0
        smg(6,4)=0.0D0
        smg(6,5)=-1.0D0/10.0D0
        smg(6,6)=2.0D0*AL/15.0D0

        smg(1,1)=1.0D0/AL
        smg(4,1)=-1.0D0/AL
        smg(4,4)=1.0D0/AL
        do i=1,6
            do j=i,6
                smg(i,j)=smg(j,i)
            end do
        end do
        !Linear stiffness
        call ELMSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sme)
        if(matno(ne)==2)call MODLSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sresul,sme,gamma0)
        do i=1,6
            do j=1,6
                sm(i,j)=sme(i,j)+P*smg(i,j)
            end do
        end do
        call HENMX(ne,xx1,yy1,xx2,yy2,hen)
        ![w]=[T]T*[k]
        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)*sm(k,j)
                end do
            end do
        end do
        ![sm]=[w]*[T]=[T]T*[k]*[T]: [sm] is stiffness matrix in global coordinate system
        do i=1,6
            do j=1,6
                sm(i,j)=0.0D0
                do k=1,6
                    sm(i,j)=sm(i,j)+w(i,k)*hen(k,j)
                end do
            end do
        end do
    end subroutine FRAMESM


    subroutine ELMSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sme)
        !*************************************************************
        !Element linear stiffness matrix
        !*************************************************************
        integer,intent(in)::ne
        real(8),intent(in)::xx1,yy1,xx2,yy2
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),AI(1:NELT)
        real(8),intent(out)::sme(1:6,1:6)

        integer::i,j
        real(8)::AL,EA,EI

        do i=1,6
            do j=1,6
                sme(i,j)=0.0D0
            end do
        end do

        AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        EA=Em(ne)*AA(ne)
        EI=Em(ne)*AI(ne)
        sme(1,1)=EA/AL
        sme(2,1)=0.0D0
        sme(3,1)=0.0D0
        sme(4,1)=-EA/AL
        sme(5,1)=0.0D0
        sme(6,1)=0.0D0
        sme(2,2)=12.0D0*EI/AL/AL/AL
        sme(3,2)=6.0D0*EI/AL/AL
        sme(4,2)=0.0D0
        sme(5,2)=-12.0D0*EI/AL/AL/AL
        sme(6,2)=6.0D0*EI/AL/AL
        sme(3,3)=4.0D0*EI/AL
        sme(4,3)=0.0D0
        sme(5,3)=-6.0D0*EI/AL/AL
        sme(6,3)=2.0D0*EI/AL
        sme(4,4)=EA/AL
        sme(5,4)=0.0D0
        sme(6,4)=0.0D0
        sme(5,5)=12.0D0*EI/AL/AL/AL
        sme(6,5)=-6.0D0*EI/AL/AL
        sme(6,6)=4.0D0*EI/AL
        do i=1,5
            do j=i+1,6
                sme(i,j)=sme(j,i)
            end do
        end do
    end subroutine ELMSM


    subroutine MODLSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sresul,sme,gamma0)
        !*************************************************************
        !Element linear stiffness matrix
        !*************************************************************
        integer,intent(in)::ne
        real(8),intent(in)::xx1,yy1,xx2,yy2
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),AI(1:NELT)
        real(8),intent(in)::sresul(1:NELT,1:6)
        real(8),intent(inout)::sme(1:6,1:6)
        real(8),intent(in)::gamma0

        integer::i,j
        real(8)::AL,EA,EI
        real(8)::Ka,Kb
        real(8)::g1,g2,A,B,C,D,F,G,H

        AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        EA=Em(ne)*AA(ne)
        EI=Em(ne)*AI(ne)

        if(gamma0==1.0D0)then
            g1=1.0D0
            g2=1.0D0
        else if(gamma0==0.0D0)then
            g1=0.0D0
            g2=0.0D0
        else
            Ka=3.0D0*EI/AL*gamma0/(1.0D0-gamma0)
            Kb=3.0D0*EI/AL*gamma0/(1.0D0-gamma0)
            g1=1.0D0/(1.0D0+3.0D0*EI/AL/Ka)
            g2=1.0D0/(1.0D0+3.0D0*EI/AL/Kb)
        end if
        A=g1+g2+g1*g2
        B=g1*(2.0D0+g2)
        C=3.0D0*g1
        D=3.0D0*g2
        F=3.0D0*g1*g2
        G=g2*(2.0D0+g1)
        H=4.0D0-g1*g2

        sme(2,2)=sme(2,2)*A/H
        sme(2,3)=sme(2,3)*B/H
        sme(2,5)=sme(2,5)*A/H
        sme(2,6)=sme(2,6)*G/H
        sme(3,2)=sme(3,2)*B/H
        sme(3,3)=sme(3,3)*C/H
        sme(3,5)=sme(3,5)*B/H
        sme(3,6)=sme(3,6)*F/H
        sme(5,2)=sme(5,2)*A/H
        sme(5,3)=sme(5,3)*B/H
        sme(5,5)=sme(5,5)*A/H
        sme(5,6)=sme(5,6)*G/H
        sme(6,2)=sme(6,2)*G/H
        sme(6,3)=sme(6,3)*F/H
        sme(6,5)=sme(6,5)*G/H
        sme(6,6)=sme(6,6)*D/H
        
    end subroutine MODLSM


    subroutine CALSTFRM(ne,kakom,matno,x,y,dis,dist,Em,AA,AI,fevec,sresul,NELT,NODT,MATEL,nt,gamma0)
        !*************************************************************
        !Calculation of stress resultant
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,MATEL,nt
        integer,intent(in)::kakom(1:NELT,1:2),matno(1:MATEL)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::dis(1:nt),dist(1:nt)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),AI(1:NELT)
        real(8),intent(out)::fevec(1:6)
        real(8),intent(inout)::sresul(1:NELT,1:6)
        real(8),intent(in)::gamma0

        integer::i,k,ii,jj,ia
        real(8)::hen(1:6,1:6),sme(1:6,1:6)
        real(8)::AL,BL
        real(8)::deltai,deltaj
        real(8)::wd(1:6),ws(1:6)
        real(8)::xx1,yy1,xx2,yy2

        !Elimination of rigid rotation
        call ROTATE(ne,kakom,x,y,dis,dist,deltai,deltaj,NELT,NODT,nt)

        ii=kakom(ne,1)
        jj=kakom(ne,2)
        !Previous element length
        xx1=x(ii)+dist(3*ii-2)-dis(3*ii-2)
        yy1=y(ii)+dist(3*ii-1)-dis(3*ii-1)
        xx2=x(jj)+dist(3*jj-2)-dis(3*jj-2)
        yy2=y(jj)+dist(3*jj-1)-dis(3*jj-1)
        BL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        !Previous linear stiffness matrix
        call ELMSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sme)
        if(matno(ne)==2)call MODLSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sresul,sme,gamma0)
        !Current element length
        xx1=x(ii)+dist(3*ii-2)
        yy1=y(ii)+dist(3*ii-1)
        xx2=x(jj)+dist(3*jj-2)
        yy2=y(jj)+dist(3*jj-1)
        AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        wd(1)=0.0D0
        wd(2)=0.0D0
        wd(3)=deltai
        wd(4)=AL-BL
        wd(5)=0.0D0
        wd(6)=deltaj
        !Internal force increment using previous linear stiffness
        do i=1,6
            ws(i)=0.0D0
            do k=1,6
                ws(i)=ws(i)+sme(i,k)*wd(k)
            end do
        end do
        !To save stress resultant
        do i=1,6
            sresul(ne,i)=sresul(ne,i)+ws(i)
        end do
        !Internal force in global coordinate system after deformation
        call HENMX(ne,xx1,yy1,xx2,yy2,hen)
        do i=1,6
            fevec(i)=0.0D0
            do k=1,6
                fevec(i)=fevec(i)+hen(k,i)*sresul(ne,k)
            end do
        end do
    end subroutine CALSTFRM


    subroutine ROTATE(ne,kakom,x,y,dis,dist,deltai,deltaj,NELT,NODT,nt)
        !*************************************************************
        !Elimination of rigid rotation
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,nt
        integer,intent(in)::kakom(1:NELT,1:2)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::dis(1:nt),dist(1:nt)
        real(8),intent(out)::deltai,deltaj

        integer::i,k,ii,jj,ia
        real(8)::hen(1:6,1:6),wd(1:6),wde(1:6)
        real(8)::xx1,yy1,xx2,yy2
        real(8)::BL,Ti,Tj,theta1i,theta1j,theta2i,theta2j

        ii=kakom(ne,1)
        jj=kakom(ne,2)
        !Accumulated displacement at previous step in global coordinate system
        do i=1,2
            ia=kakom(ne,i)
            wd(3*i-2)=dist(3*ia-2)-dis(3*ia-2)
            wd(3*i-1)=dist(3*ia-1)-dis(3*ia-1)
            wd(3*i)=dist(3*ia)-dis(3*ia)
        end do
        !Initial element length
        xx1=x(ii)
        yy1=y(ii)
        xx2=x(jj)
        yy2=y(jj)
        BL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        call HENMX(ne,xx1,yy1,xx2,yy2,hen)
        do i=1,6
            wde(i)=0.0D0
            do k=1,6
                wde(i)=wde(i)+hen(i,k)*wd(k)
            end do
        end do
        !Rotation angle at previous step in local coordinate system
        Ti=tan(wde(3))
        Tj=tan(wde(6))
        theta1i=((BL+wde(4)-wde(1))*Ti-(wde(5)-wde(2)))/((BL+wde(4)-wde(1))+(wde(5)-wde(2))*Ti)
        theta1j=((BL+wde(4)-wde(1))*Tj-(wde(5)-wde(2)))/((BL+wde(4)-wde(1))+(wde(5)-wde(2))*Tj)
        !Accumulated displacement at current step in global coordinate system
        do i=1,2
            ia=kakom(ne,i)
            wd(3*i-2)=dist(3*ia-2)
            wd(3*i-1)=dist(3*ia-1)
            wd(3*i)=dist(3*ia)
        end do
        do i=1,6
            wde(i)=0.0D0
            do k=1,6
                wde(i)=wde(i)+hen(i,k)*wd(k)
            end do
        end do
        !Rotation angle at current step in local coordinate system
        Ti=tan(wde(3))
        Tj=tan(wde(6))
        theta2i=((BL+wde(4)-wde(1))*Ti-(wde(5)-wde(2)))/((BL+wde(4)-wde(1))+(wde(5)-wde(2))*Ti)
        theta2j=((BL+wde(4)-wde(1))*Tj-(wde(5)-wde(2)))/((BL+wde(4)-wde(1))+(wde(5)-wde(2))*Tj)
        !Rotation angle increment for calculation of internal force
        deltai=theta2i-theta1i
        deltaj=theta2j-theta1j
    end subroutine ROTATE


    subroutine HENMX(ne,xx1,yy1,xx2,yy2,hen)
        !*************************************************************
        !Transform matrix
        !*************************************************************
        integer,intent(in)::ne
        real(8),intent(in)::xx1,yy1,xx2,yy2
        real(8),intent(out)::hen(1:6,1:6)

        integer::i,j
        real(8)::AL,cs,sn

        AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        cs=(xx2-xx1)/AL
        sn=(yy2-yy1)/AL
        do i=1,6
            do j=1,6
                hen(i,j)=0.0D0
            end do
        end do
        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
    end subroutine HENMX


    subroutine ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !*************************************************************
        !assembly of vector
        !*************************************************************
        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
        !*************************************************************
        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 MATINV(n,ary,nt)
        !*************************************
        !Inverse matrix
        !*************************************
        integer,intent(in)::n,nt
        real(8),intent(inout)::ary(1:nt,1:nt)

        integer::i,j,k,l,l1,irow,icol,jrow,jcol
        real(8)::aamax,s,t
        real(8)::factor=1.0D-30
        integer::ipivot(n)
        real(8)::pivot(n)
        integer::nindex(n,2)

        do j=1,n
            ipivot(j)=0
        end do
        do i=1,n
            aamax=0.0D0
            do j=1,n
                if(ipivot(j)-1/=0)then
                    do k=1,n
                        if(0<ipivot(k)-1)return
                        if(ipivot(k)-1<0)then
                            if(abs(aamax)-abs(ary(j,k))<0)then
                                irow=j
                                icol=k
                                aamax=ary(j,k)
                            end if
                        end if
                    end do
                end if
            end do
            ipivot(icol)=ipivot(icol)+1
            if(irow-icol/=0)then
                do l=1,n
                    s=ary(irow,l)
                    ary(irow,l)=ary(icol,l)
                    ary(icol,l)=s
                end do
            end if
            nindex(i,1)=irow
            nindex(i,2)=icol
            pivot(i)=ary(icol,icol)
            if(abs(pivot(i))<factor)then
                write(6,*) '#Singular'
                stop
            end if
            ary(icol,icol)=1.0D0
            do l=1,n
                ary(icol,l)=ary(icol,l)/pivot(i)
            end do
            do l1=1,n
                if(l1-icol/=0)then
                    t=ary(l1,icol)
                    ary(l1,icol)=0.0D0
                    do l=1,n
                        ary(l1,l)=ary(l1,l)-ary(icol,l)*t
                    end do
                end if
            end do
        end do
        do i=1,n
            l=n+1-i
            if(nindex(l,1)-nindex(l,2)/=0)then
                jrow=nindex(l,1)
                jcol=nindex(l,2)
                do k=1,n
                    s=ary(k,jrow)
                    ary(k,jrow)=ary(k,jcol)
                    ary(k,jcol)=s
                end do
            end if
        end do
    end subroutine MATINV


    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_gfrm
