program f90_FRAME
    !平面骨組構造解析
    implicit none
    integer::i,j,k,l,m,kk,ne,n,ia,ja
    character::strcom*50                !書出用コメント
    integer::NODT                       !節点総数
    integer::NELT                       !要素総数
    integer::MATEL                      !材料種類数
    integer::KOX                        !x方向変位既知節点数
    integer::KOY                        !y方向変位既知節点数
    integer::KOZ                        !回転変位既知節点数
    integer::NF                         !載荷節点数
    integer,allocatable::kakom(:,:)     !要素構成接点番号
    integer,allocatable::matno(:)       !材料種別No
    real(8),allocatable::wm(:,:)        !材料物性作業用
    real(8),allocatable::Em(:)          !要素弾性係数
    real(8),allocatable::AA(:)          !要素断面積
    real(8),allocatable::AI(:)          !要素断面２次モーメント
    real(8),allocatable::gamma(:)       !要素単位体積重量(質量ではない)
    real(8),allocatable::gkh(:)         !要素の水平方向物体力係数(重力加速度に対する比率)
    real(8),allocatable::gkv(:)         !要素の鉛直方向物体力係数(重力加速度に対する比率)
    real(8),allocatable::alpha(:)       !要素線膨張係数
    real(8),allocatable::x(:)           !節点x座標
    real(8),allocatable::y(:)           !節点y座標
    real(8),allocatable::deltaT(:)      !節点温度変化量(＋：温度上昇)
    integer,allocatable::nokx(:)        !x方向変位既知節点番号
    integer,allocatable::noky(:)        !y方向変位既知節点番号
    integer,allocatable::nokz(:)        !回転変位既知節点番号
    real(8),allocatable::f(:)           !全体外力ベクトル
    integer,allocatable::index(:)       !既知節点変位処理のためのインデックス
    real(8),allocatable::ftvec(:)       !全体荷重ベクトル
    real(8),allocatable::ftemp(:)      !全体温度ひずみ相当荷重ベクトル
    real(8),allocatable::fmass(:)       !全体物体力ベクトル
    real(8),allocatable::fdis(:)        !既知節点変位変
    real(8),allocatable::rdis(:)        !既知節点変位による荷重成分格納
    real(8),allocatable::tsm(:,:)       !全体剛性マトリックス
    real(8),allocatable::dis(:)         !全体節点変位増分
    real(8),allocatable::reac(:)        !作業用配列
    real(8),allocatable::sm(:,:)        !要素剛性マトリックス
    real(8),allocatable::hen(:,:)       !座標変換マトリックス
    real(8),allocatable::fevec(:)       !要素荷重ベクトル
    integer::nt                         !全自由度（総節点数×２[１節点自由度]）
    integer::nod=2                      !1要素節点数
    integer::nhen=3                     !1節点自由度
    integer::mm                         !既知節点変位処理後自由度（逆行列の寸法）
    integer::ib                         !バンド幅
    real(8)::fact0=1.0D-30              !0判定

    character :: linebuf*1000
    character :: fmt1*200,fmt2*200,fmt3*200,fmt4*200
    real(4) :: t1,t2
    integer :: iv(1:8)
    character date*10,time1*20,time2*20,zone*10

    fmt1="(i5,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7,&
            ',',i2,',',i2,',',i2,',',e15.7,',',e15.7,',',e15.7,',',e15.7)"    !節点入力値
    fmt2="(i5,',',i5,',',i5,',',e15.7,',',e15.7,',',e15.7,&
            ',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',i5)"                  !要素特性
    fmt3="(i5,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7,&
            ',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7)"     !節点変位・反力出力
    fmt4="(i5,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7)"   !断面力出力

    !時間計測開始
    call DATE_AND_TIME(date,time1,zone,iv)
    read(time1,*) t1
    !***************************
    !data input
    !***************************
    open(11,file='fnameR.csv',status='old')
    !/*--------------------------------------------------------------------------------------*/
    !/*書出用コメント入力*/
    !/*--------------------------------------------------------------------------------------*/
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
    !節点数，要素数，材料種類数,X方向変位既知節点数，Y方向変位既知節点数，回転変位既知節点数，載荷点数
    !------------------------------------------------------------------------------------------
    read(11,*) NODT,NELT,MATEL,KOX,KOY,KOZ,NF
    !------------------------------------------------------------------------------------------
    !配列寸法宣言
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:7))
    allocate(Em(1:NELT))
    allocate(AA(1:NELT))
    allocate(AI(1:NELT))
    allocate(gamma(1:NELT))
    allocate(gkh(1:NELT))
    allocate(gkv(1:NELT))
    allocate(alpha(1:NELT))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(deltaT(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(index(1:nt))
    allocate(ftvec(1:nt))
    allocate(ftemp(1:nt))
    allocate(fmass(1:nt))
    allocate(fdis(1:nt))
    allocate(rdis(1:nt))
    allocate(tsm(1:nt,1:nt))
    allocate(dis(1:nt))
    allocate(reac(1:nt))
    allocate(sm(1:nod*nhen,1:nod*nhen))
    allocate(hen(1:nod*nhen,1:nod*nhen))
    allocate(fevec(1:nod*nhen))
    !------------------------------------------------------------------------------------------
    !材料特性（Em,AA,AI,gamma,gkh,gkv,alpha）
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3),wm(i,4),wm(i,5),wm(i,6),wm(i,7)
    end do
    !------------------------------------------------------------------------------------------
    !要素構成節点と材料種別（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)
                gkh(ne)  =wm(i,5)
                gkv(ne)  =wm(i,6)
                alpha(ne)=wm(i,7)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !節点座標（x,y），温度変化量
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),y(i),deltaT(i)
    end do
    !------------------------------------------------------------------------------------------
    !変位既知節点番号，既知節点変位量
    !------------------------------------------------------------------------------------------
    do i=1,nt
        rdis(i)=0.0D0
    end do
    if(0<KOX)then
        do i=1,KOX
            read(11,*) nokx(i),rdis(3*nokx(i)-2)
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            read(11,*) noky(i),rdis(3*noky(i)-1)
        end do
    end if
    if(0<KOZ)then
        do i=1,KOZ
            read(11,*) nokz(i),rdis(3*nokz(i))
        end do
    end if
    !------------------------------------------------------------------------------------------
    !節点番号,x方向荷重,y方向荷重
    !------------------------------------------------------------------------------------------
    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
    close(11)

    !*****************************************************
    !入力確認出力
    !*****************************************************
    open(12, file='fnameW.csv', status='replace')
    write(12,'(a)') strcom
    write(12,'(a)') 'NODT,NELT,MATEL,KOX,KOY,KOZ,NF'
        write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',KOZ,',',NF
        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,rdis-x,rdis-y,rdis-z,deltaT'
    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,rdis(3*i-2),rdis(3*i-1),rdis(3*i),deltaT(i)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(12,'(a)') '*element characteristics'
    write(12,'(a)') 'element,node-1,node-2,E,A,I,gamma,kh,kv,alpha,matno'
    do ne=1,NELT
            write(linebuf,FMT=fmt2) ne,kakom(ne,1),kakom(ne,2),Em(ne),AA(ne),AI(ne),&
                gamma(ne),gkh(ne),gkv(ne),alpha(ne),matno(ne)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)
    !*****************************************************
    !メイン処理計算
    !*****************************************************
    !-----------------------------------------------------
    !全体剛性行列・全体荷重ベクトルクリア
    !-----------------------------------------------------
    do i=1,nt
        ftemp(i)=0.D00
        fmass(i)=0.0D0
        ftvec(i)=0.0D0
        reac(i)=0.0D0
        do j=1,nt
            tsm(i,j)=0.0D0
        end do
    end do
    !-----------------------------------------------------
    !要素剛性行列・荷重ベクトル組立
    !-----------------------------------------------------
    do ne=1,NELT
        !要素剛性行列
        call SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)

        !温度歪み相当荷重ベクトル，温度歪み相当応力(strsave)
        call TV_FRM(ne,kakom,deltaT,Em,AA,alpha,fevec,hen,NELT,NODT)
        call ASVEC(ne,nod,nhen,kakom,fevec,ftemp,NELT,NODT)

        !物体力ベクトル
        call MV_FRM(ne,kakom,x,y,AA,gamma,gkh,gkv,fevec,hen,NELT,NODT)
        call ASVEC(ne,nod,nhen,kakom,fevec,fmass,NELT,NODT)
    end do
    !-----------------------------------------------------
    !既知節点変位による荷重成分ベクトル作成
    !-----------------------------------------------------
    do i=1,nt
        fdis(i)=0.0D0
        do j=1,nt
            fdis(i)=fdis(i)+rdis(j)*tsm(i,j)
        end do
    end do
    !-----------------------------------------------------
    !荷重ベクトルセット
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=f(i)+fmass(i)+ftemp(i)-fdis(i)
    end do
    !-----------------------------------------------------
    !変位既知節点自由度処理
    !-----------------------------------------------------
    do i=1,nt
        index(i)=i
        reac(i)=ftvec(i)
    end do
    if(0<KOX)then
        do i=1,KOX
            n=3*nokx(i)-2
            index(n)=0
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            n=3*noky(i)-1
            index(n)=0
        end do
    end if
    if(0<KOZ)then
        do i=1,KOZ
            n=3*nokz(i)
            index(n)=0
        end do
    end if
    mm=0
    do i=1,nt
        if(index(i)>0)then
            mm=mm+1
            index(mm)=i
        end if
    end do
    do i=1,mm
        ia=index(i)
        reac(i)=reac(ia)
        do j=1,mm
            ja=index(j)
            tsm(i,j)=tsm(ia,ja)
        end do
    end do
    !-----------------------------------------------------
    !バンド幅計算とフルマトリックスのバンド化記憶
    !-----------------------------------------------------
    ib=IBAND(mm,tsm,fact0,nt)
    call BAND(mm,ib,tsm,nt)
    !-----------------------------------------------------
    !対角項確認と異常終了通知
    !-----------------------------------------------------
    do i=1,mm
        if(abs(tsm(i,1))<fact0)stop
    end do
    !-----------------------------------------------------
    !連立一次方程式（コレスキー法）
    !-----------------------------------------------------
    call CHBAND10(mm,ib,tsm,nt)!三角行列
    call CHBAND11(mm,ib,tsm,reac,nt)!前進消去・後退代入

    !*****************************************************
    !応力計算
    !*****************************************************
    do i=1,nt
        dis(i)=0.0D0
    end do
    do i=1,mm
        ia=index(i)
        dis(ia)=reac(i)
    end do
    do i=1,nt
        reac(i)=0.0D0
        dis(i)=dis(i)+rdis(i)
    end do
    do ne=1,NELT
        call SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen,NELT,NODT)
        call CALST_FRM(ne,kakom,dis,sm,hen,reac,tsm,Em,AA,alpha,deltaT,NELT,NODT,nt)
    end do
    !*****************************************************
    !結果出力
    !*****************************************************
    open(12, file='fnameW.csv', access='append')
    !節点変位
    write(12,'(a)') '*displacements and forces'
    write(12,'(a)') 'node,coord-x,coord-y,dir-x,dir-y,dir-z,reac-x,reac-y,reac-z,ftvec-x,ftvec-y,ftvec-z'
    do i=1,NODT
        write(linebuf,FMT=fmt3) i,x(i),y(i),dis(3*i-2),dis(3*i-1),dis(3*i),&
            reac(3*i-2),reac(3*i-1),reac(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,tsm(ne,1),tsm(ne,2),tsm(ne,3),tsm(ne,4),tsm(ne,5),tsm(ne,6)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do

    !時間計測完了
    call DATE_AND_TIME(date,time2,zone,iv)
    read(time2,*) t2
    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=',iv(1),'.',iv(2),'.',iv(3),'_',iv(5),':',iv(6),':',iv(7)
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    close(12)

stop

contains

    subroutine SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen,NELT,NODT)
        !要素剛性行列
        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)::sm(1:6,1:6),hen(1:6,1:6)
        integer :: i,j
        real(8) :: AL,EA,EI,cs,sn,w(1:6,1:6)

        !行列初期化
        do i=1,6
            do j=1,6
                hen(i,j)=0.0D0
                sm(i,j)=0.0D0
            end do
        end do

        !座標変換行列
        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

        !要素剛性行列
        EA=Em(ne)*AA(ne)
        EI=Em(ne)*AI(ne)
        sm(1,1)=EA/AL
        sm(2,1)=0.0D0
        sm(3,1)=0.0D0
        sm(4,1)=-EA/AL
        sm(5,1)=0.0D0
        sm(6,1)=0.0D0
        sm(2,2)=12.0D0*EI/AL/AL/AL
        sm(3,2)=6.0D0*EI/AL/AL
        sm(4,2)=0.0D0
        sm(5,2)=-12.0D0*EI/AL/AL/AL
        sm(6,2)=6.0D0*EI/AL/AL
        sm(3,3)=4.0D0*EI/AL
        sm(4,3)=0.0D0
        sm(5,3)=-6.0D0*EI/AL/AL
        sm(6,3)=2.0D0*EI/AL
        sm(4,4)=EA/AL
        sm(5,4)=0.0D0
        sm(6,4)=0.0D0
        sm(5,5)=12.0D0*EI/AL/AL/AL
        sm(6,5)=-6.0D0*EI/AL/AL
        sm(6,6)=4.0D0*EI/AL
        !対称項セット
        do i=1,5
            do j=i+1,6
                sm(i,j)=sm(j,i)
            end do
        end do

        !要素剛性行列の全体座標系への変換
        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
        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 SM_FRM


    subroutine MV_FRM(ne,kakom,x,y,AA,gamma,gkh,gkv,fevec,hen,NELT,NODT)
        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)::AA(1:NELT),gamma(1:NELT),gkh(1:NELT),gkv(1:NELT)
        real(8),intent(in)::hen(1:6,1:6)
        real(8),intent(out)::fevec(1:6)
        integer :: i,j
        real(8) :: AL,am(1:6,1:6),w(1:6,1:6),work(1:6)
        real(8) :: ggg=9.8D0

        !行列初期化
        do i=1,6
            do j=1,6
                am(i,j)=0.0D0
            end do
        end do

        !要素長さ
        i=kakom(ne,1)
        j=kakom(ne,2)
        AL=sqrt((x(j)-x(i))**2.0D0+(y(j)-y(i))**2.0D0)

        !質量行列
        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

        !全体座標系への変換
        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

        work(1)=ggg*gkh(ne)
        work(2)=ggg*gkv(ne)
        work(3)=0.0D0
        work(4)=ggg*gkh(ne)
        work(5)=ggg*gkv(ne)
        work(6)=0.0D0

        !全体座標系での慣性力ベクトル
        do i=1,6
            fevec(i)=0.0D0
            do j=1,6
                fevec(i)=fevec(i)+am(i,j)*work(j)
            end do
        end do
    end subroutine MV_FRM


    subroutine TV_FRM(ne,kakom,deltaT,Em,AA,alpha,fevec,hen,NELT,NODT)
        !温度荷重相当ベクトル
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:2)
        real(8),intent(in)::deltaT(1:NODT)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),alpha(1:NELT)
        real(8),intent(in)::hen(1:6,1:6)
        real(8),intent(out)::fevec(1:6)
        integer :: i,j
        real(8) :: tempe,work(1:6)

        i=kakom(ne,1)
        j=kakom(ne,2)
        tempe=0.5D0*(deltaT(i)+deltaT(j))
        !要素座標系での温度荷重相当ベクトル
        work(1)=-1.0D0*Em(ne)*AA(ne)*alpha(ne)*tempe
        work(2)=0.0D0
        work(3)=0.0D0
        work(4)=1.0D0*Em(ne)*AA(ne)*alpha(ne)*tempe
        work(5)=0.0D0
        work(6)=0.0D0
        !全体座標系での温度荷重相当ベクトル
        !（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる）
        do i=1,6
            fevec(i)=0.0D0
            do j=1,6
                fevec(i)=fevec(i)+hen(j,i)*work(j)
            end do
        end do
    end subroutine TV_FRM


    subroutine CALST_FRM(ne,kakom,dis,sm,hen,reac,tsm,Em,AA,alpha,deltaT,NELT,NODT,nt)
        !断面力算定
        integer,intent(in)::NELT,NODT,nt
        integer,intent(in)::ne,kakom(1:NELT,1:2)
        real(8),intent(in)::dis(1:nt),sm(1:6,1:6),hen(1:6,1:6)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),alpha(1:NELT),deltaT(1:NODT)
        real(8),intent(inout)::reac(1:nt)
        real(8),intent(out)::tsm(1:nt,1:nt)
        integer :: i,j,k,ia,ir,iw
        real(8) :: wd(1:6),ws(1:6),tempe

        !要素節点変位取得（全体座標系）
        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
        !全体座標系での要素節点力算定
        do i=1,6
            ws(i)=0.0D0
            do k=1,6
                ws(i)=ws(i)+sm(i,k)*wd(k)
            end do
        end do
        !要素断面力算定：座標変換行列×全体座標系要素節点力
        !今後使用しない全体剛性行列格納配列(tsm)に断面力を格納
        do i=1,6
            tsm(ne,i)=0.0D0
            do k=1,6
                tsm(ne,i)=tsm(ne,i)+hen(i,k)*ws(k)
            end do
        end do
        !軸力の温度応力補正
        i=kakom(ne,1)
        j=kakom(ne,2)
        tempe=0.5D0*(deltaT(i)+deltaT(j))
        tsm(ne,1)=tsm(ne,1)+Em(ne)*AA(ne)*alpha(ne)*tempe
        tsm(ne,4)=tsm(ne,4)-Em(ne)*AA(ne)*alpha(ne)*tempe
        !全体座標系での要素断面力
        !（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる）
        do i=1,6
            ws(i)=0.0D0
            do j=1,6
                ws(i)=ws(i)+hen(j,i)*tsm(ne,j)
            end do
        end do
        !要素節点力の全体系へのセット
        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 ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !荷重ベクトル組立(nod：1要素節点数，nhen：1節点の自由度)
        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)
        !剛性行列組立(nod：1要素節点数，nhen：1節点の自由度)
        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


end program f90_FRAME