program f90_FEM4nodTS
    !２次元四角形要素応力解析
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    integer::i,j,k,l,kk,ne,n,ia,ja
    character::strcom*50                !書出用コメント
    integer::NSTRES                     !応力状態（0:平面歪み,1:平面応力）
    integer::NODT                       !節点総数
    integer::NELT                       !要素総数
    integer::MATEL                      !材料種類数
    integer::KOX                        !x方向変位既知節点数
    integer::KOY                        !y方向変位既知節点数
    integer::NF                         !載荷節点数
    integer,allocatable::kakom(:,:)     !要素構成接点番号
    integer,allocatable::matno(:)       !材料種別No
    real(8),allocatable::wm(:,:)        !材料物性作業用
    real(8),allocatable::Em(:)          !要素弾性係数
    real(8),allocatable::po(:)          !要素ポアソン比
    real(8),allocatable::t(:)           !要素板厚
    real(8),allocatable::gamma(:)       !要素単位体積重量(質量ではない)
    real(8),allocatable::gkh(:)         !要素の水平方向物体力係数(重力加速度に対する比率)
    real(8),allocatable::gkv(:)         !要素の鉛直方向物体力係数(重力加速度に対する比率)
    real(8),allocatable::alpha(:)       !要素線膨張係数
    real(8),allocatable::ts(:)          !要素引張強度
    real(8),allocatable::x(:)           !節点x座標
    real(8),allocatable::y(:)           !節点y座標
    real(8),allocatable::deltaT(:)      !節点温度変化量(＋：温度上昇)
    integer,allocatable::nokx(:)        !x方向変位既知節点番号
    integer,allocatable::noky(:)        !y方向変位既知節点番号
    real(8),allocatable::f(:)           !全体外力ベクトル
    integer,allocatable::index(:)       !既知節点変位処理のためのインデックス
    real(8),allocatable::ftvec(:)       !全体荷重ベクトル
    real(8),allocatable::fmass(:)       !全体物体力ベクトル
    real(8),allocatable::rdis(:)        !既知節点変位変
    real(8),allocatable::tsm(:,:)       !全体剛性マトリックス
    real(8),allocatable::dis(:)         !全体節点変位増分
    real(8),allocatable::wdis(:)        !全体節点変位作業用
    real(8),allocatable::dist(:)        !全体節点変位
    real(8),allocatable::reac(:)        !内力ベクトル
    real(8),allocatable::strsave(:,:,:) !応力記憶用
    integer,allocatable::noten(:,:)     !応力状態記憶用
    real(8),allocatable::sm(:,:)        !要素剛性マトリックス
    real(8),allocatable::fevec(:)       !要素荷重ベクトル
    integer::nt                         !全自由度（総節点数×２(１節点自由度)）
    integer::nod=4                      !1要素節点数
    integer::nhen=2                     !1節点自由度
    integer::mm                         !既知節点変位処理後自由度（逆行列の寸法）
    integer::ib                         !バンド幅
    real(8)::elarea1                    !要素半面積
    real(8)::elarea2                    !要素半面積
    real(8)::fact0=1.0D-30              !0判定
    real(8):: sigx,sigy,tauxy,ps1,ps2,ang
    integer:: nnn
    integer:: maxita=2000
    integer:: ntsum
    integer:: icount
    real(8):: dtest,ftest

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

    fmt1="(i5,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',i2,',',i2,',',e15.7,',',e15.7,',',e15.7)"    !節点入力値
    fmt2="(i5,',',i5,',',i5,',',i5,',',i5,',',e15.7,',',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)"     !節点変位・反力出力
    fmt4="(i5,',',i5,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',i5,',',i5)"   !応力出力


    !時間計測開始
    call DATE_AND_TIME(date,time1,zone,iv)
    read(time1,*) t1
    !***************************
    !data input
    !***************************
    open(11,file='fnameR.csv',status='old')
    !!--------------------------------------------------------------------------------------
    !!書出用コメント入力
    !!--------------------------------------------------------------------------------------
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
    !応力状態(0:平面ひずみ)，節点数，要素数，材料種類数,X方向変位既知節点数，Y方向変位既知節点数，載荷点数
    !------------------------------------------------------------------------------------------
    read(11,*) NSTRES,NODT,NELT,MATEL,KOX,KOY,NF,IPR
    !------------------------------------------------------------------------------------------
    !配列寸法宣言
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:8))
    allocate(Em(1:NELT))
    allocate(po(1:NELT))
    allocate(t(1:NELT))
    allocate(gamma(1:NELT))
    allocate(gkh(1:NELT))
    allocate(gkv(1:NELT))
    allocate(alpha(1:NELT))
    allocate(ts(1:NELT))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(deltaT(1:NODT))
    allocate(nokx(1:KOX))
    allocate(noky(1:KOY))
    allocate(f(1:nt))
    allocate(index(1:nt))
    allocate(ftvec(1:nt))
    allocate(fmass(1:nt))
    allocate(rdis(1:nt))
    allocate(tsm(1:nt,1:nt))
    allocate(dis(1:nt))
    allocate(wdis(1:nt))
    allocate(dist(1:nt))
    allocate(reac(1:nt))
    allocate(strsave(1:NELT,1:4,1:3))
    allocate(noten(1:NELT,1:4))
    allocate(sm(1:nod*nhen,1:nod*nhen))
    allocate(fevec(1:nod*nhen))
    !------------------------------------------------------------------------------------------
    !材料特性（t,Em,po,gamma,gkh,gkv,alpha,ts）
    !------------------------------------------------------------------------------------------
    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),wm(i,8)
    end do
    !------------------------------------------------------------------------------------------
    !要素構成節点と材料種別（node-1,node-2,node-3,node-4,matno）
    !------------------------------------------------------------------------------------------
    do ne=1,NELT
        read(11,*) kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4),matno(ne)
        do i=1,MATEL
            if(i==matno(ne))then
                t(ne)    =wm(i,1)
                Em(ne)   =wm(i,2)
                po(ne)   =wm(i,3)
                gamma(ne)=wm(i,4)
                gkh(ne)  =wm(i,5)
                gkv(ne)  =wm(i,6)
                alpha(ne)=wm(i,7)
                ts(ne)   =wm(i,8)
            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(2*nokx(i)-1)
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            read(11,*) noky(i),rdis(2*noky(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(2*n-1),f(2*n)
        end do
    end if
    close(11)

    !*****************************************************
    !入力確認出力
    !*****************************************************
    open(12, file='fnameW.csv', status='replace')
    write(12,*) strcom
    write(12,*) 'NSTRES,NODT,NELT,MATEL,KOX,KOY,NF,IPR'
        write(linebuf,*) NSTRES,',',NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',NF,',',IPR
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,*) '*node characteristics'
    write(12,*) 'node,x,y,fx,fy,fix-x,fix-y,rdis-x,rdis-y,deltaT'
    do i=1,NODT
        k=0
        l=0
        if(0<KOX)then
            do j=1,KOX
                if(i==nokx(j))k=1
            end do
        end if
        if(0<KOY)then
            do j=1,KOY
                if(i==noky(j))l=1
            end do
        end if
            write(linebuf,FMT=fmt1) i,x(i),y(i),f(2*i-1),f(2*i),k,l,rdis(2*i-1),rdis(2*i),deltaT(i)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(12,*) '*element characteristics'
    write(12,*) 'element,node-1,node-2,node-3,node-4,E,po,t,gamma,kh,kv,alpha,ts,matno'
    do ne=1,NELT
            write(linebuf,FMT=fmt2) ne,kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4), &
                Em(ne),po(ne),t(ne),gamma(ne),gkh(ne),gkv(ne),alpha(ne),ts(ne),matno(ne)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)
    !要素面積確認と異常終了通知
    do ne=1,NELT
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        elarea1=0.5*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        if(elarea1<fact0)stop
        i=kakom(ne,1)
        j=kakom(ne,3)
        k=kakom(ne,4)
        elarea2=0.5*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        if(elarea2<fact0)stop
    end do
    !*****************************************************
    !メイン処理計算
    !*****************************************************
    !-----------------------------------------------------
    !保存応力クリア
    !-----------------------------------------------------
    do ne=1,NELT
        do kk=1,4
            do i=1,3 
                strsave(ne,kk,i)=0.0D0
            end do
            noten(ne,kk)=0
        end do
    end do
    !-----------------------------------------------------
    !全体剛性行列・全体荷重ベクトルクリア
    !-----------------------------------------------------
    do i=1,nt
        fmass(i)=0.0D0
        ftvec(i)=0.0D0
        wdis(i)=0.0D0
        dist(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_PL4(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)

        !物体力ベクトル
        call MV_PL4(ne,kakom,x,y,t,gamma,gkh,gkv,fevec,NELT,NODT)
        call ASVEC(ne,nod,nhen,kakom,fevec,fmass,NELT,NODT)
    end do
    !-----------------------------------------------------
    !荷重ベクトルセット
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=f(i)+fmass(i)
    end do
    !-----------------------------------------------------
    !変位既知節点自由度処理
    !-----------------------------------------------------
    do i=1,nt
        index(i)=i
    end do
    if(0<KOX)then
        do i=1,KOX
            n=2*nokx(i)-1
            index(n)=0
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            n=2*noky(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)
        ftvec(i)=ftvec(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)!三角行列


    !****************************************************
    !No-tension収束計算
    !****************************************************
    nnn=0
    do while(nnn<=maxita)
        nnn=nnn+1
        !-----------------------------------------------------
        !連立一次方程式（コレスキー法）
        !-----------------------------------------------------
        do i=1,mm
            reac(i)=ftvec(i)
        end do
        call CHBAND11(mm,ib,tsm,reac,nt) !前進消去・後退代入
        !-----------------------------------------------------
        !変位増分計算
        !-----------------------------------------------------
        do i=1,nt
            dis(i)=0.0
            dist(i)=0.0
        end do
        do i=1,mm
            ia=index(i)
            dis(ia)=reac(i)
        end do
        do i=1,nt
            reac(i)=0.0
            wdis(i)=wdis(i)+dis(i)
            dist(i)=wdis(i)+rdis(i)
        end do
        !-----------------------------------------------------
        !応力計算・内力ベクトル組立
        !-----------------------------------------------------
        do ne=1,NELT
            call CALST_PL4TS(ne,kakom,x,y,deltaT,dist,NSTRES,Em,po,t,alpha,ts,strsave,fevec,noten,NELT,NODT,nt)
            call ASVEC(ne,nod,nhen,kakom,fevec,reac,NELT,NODT)!内力ベクトル組立
        end do
        !-----------------------------------------------------
        !不平衡力ベクトル算定
        !-----------------------------------------------------
        do i=1,nt
            ftvec(i)=f(i)+fmass(i)-reac(i)
        end do
        !-----------------------------------------------------
        !不平衡力ベクトル拘束節点処理
        !-----------------------------------------------------
        do i=1,mm
            ia=index(i)
            ftvec(i)=ftvec(ia)
            dis(i)=dis(ia)
            dist(i)=dist(ia)
        end do
        !-----------------------------------------------------
        !収束判定
        !-----------------------------------------------------
        !No-tension積分点の数
        ntsum = 0
        do ne=1,NELT
            do kk=1,4
                ntsum=ntsum+abs(noten(ne,kk))
            end do
        end do
        !変位増分絶対値・不平衡力絶対値総和
        icount=0
        dtest=0.0
        ftest=0.0
        do i=1,mm
            dtest=dtest+abs(dis(i))   !変位増分絶対値総和
            ftest=ftest+abs(ftvec(i)) !不平衡力絶対値総和
            if(abs(dis(i)/dist(i))<1e-6)icount=icount+1
        end do

        print "('nnn=',i0,' dtest=',D10.4,' ftest=',D10.4,' icount=',i0)",nnn,dtest,ftest,icount

        !****************************************************
        if(2<=nnn)then
            if(ntsum==0)exit
            if(icount==mm)exit
        end if
        !****************************************************
    end do

    ! 出力用全変位調整 
    do i=1,nt
        dis(i)=0.0
    end do
    do i=1,mm
        ia=index(i)
        dis(ia)=dist(i)
    end do
    do i=1,nt
        dist(i)=dis(i)+rdis(i)
    end do

    !*****************************************************
    !結果出力
    !*****************************************************
    open(12, file='fnameW.csv', access='append')
    !節点変位
    write(12,*) '*displacements and forces'
    write(12,*) 'node,coord-x,coord-y,dist-x,dist-y,reac-X,reac-y,ftvec-x,ftvec-y'
    do i=1,NODT
        write(linebuf,FMT=fmt3) i,x(i),y(i),dist(2*i-1),dist(2*i),reac(2*i-1),reac(2*i),&
            f(2*i-1)+fmass(2*i-1)-reac(2*i-1),f(2*i)+fmass(2*i)-reac(2*i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !要素平均応力
    write(12,*) '*stresses'
    write(12,*) 'element,kk,sig-x,sig-y,tau-xy,ps1,ps2,ang,noten,matno'
    select case (IPR)
        case(0)
            do ne=1,NELT
                do kk=1,4
                    sigx=strsave(ne,kk,1)
                    sigy=strsave(ne,kk,2)
                    tauxy=strsave(ne,kk,3)
                    call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
                    write(linebuf,FMT=fmt4) ne,kk,sigx,sigy,tauxy,ps1,ps2,ang,noten(ne,kk),matno(ne)
                    call del_spaces(linebuf)
                    write (12,'(a)') trim(linebuf)
                end do
            end do
        case(1)
            do ne=1,NELT
                ntsum=0
                sigx=0.0D0
                sigy=0.0D0
                tauxy=0.0D0
                do kk=1,4
                    sigx=sigx+strsave(ne,kk,1)
                    sigy=sigy+strsave(ne,kk,2)
                    tauxy=tauxy+strsave(ne,kk,3)
                    ntsum=ntsum+noten(ne,kk)
                end do
                kk=0
                sigx=0.25D0*sigx
                sigy=0.25D0*sigy
                tauxy=0.25D0*tauxy
                call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
                write(linebuf,FMT=fmt4) ne,kk,sigx,sigy,tauxy,ps1,ps2,ang,ntsum,matno(ne)
                call del_spaces(linebuf)
                write (12,'(a)') trim(linebuf)
            end do
    end select

    !時間計測完了
    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,'("nnn=",i0," dtest=",D10.4," ftest=",D10.4)') nnn,dtest,ftest
    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 DMAT_PL(NSTRES,Eme,poe,d)
        !応力-ひずみ関係行列(De)
        integer :: NSTRES
        real(8) :: Eme,poe,d(1:3,1:3)

        integer i,j

        do i=1,3
            do j=1,3
                d(i,j)=0.0D0
            end do
        end do
        select case (NSTRES)
            case (0) !平面ひずみ
                d(1,1)=Eme*(1.0D0-poe)/((1.0D0+poe)*(1.0D0-2.0D0*poe))
                d(1,2)=Eme*poe/((1.0D0+poe)*(1.0D0-2.0D0*poe))
                d(2,1)=d(1,2)
                d(2,2)=d(1,1)
                d(3,3)=Eme*0.5D0/(1.0D0+poe)
            case (1) !平面応力
                d(1,1)=Eme/(1.0D0-poe*poe)
                d(1,2)=poe*d(1,1)
                d(2,1)=d(1,2)
                d(2,2)=d(1,1)
                d(3,3)=Eme*0.5D0/(1.0D0+poe)
        end Select
    end subroutine DMAT_PL


    subroutine BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
        !ひずみ-変位関係行列(B)
        integer :: NELT,NODT
        integer :: ne,kakom(NELT,4)
        real(8) :: x(NODT),y(NODT)
        real(8) :: bm(3,8),detJ,a,b

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

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

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

        !Jacobiマトリックスとdet(J)
        J11=dn1a*x(i)+dn2a*x(j)+dn3a*x(k)+dn4a*x(l)
        J12=dn1a*y(i)+dn2a*y(j)+dn3a*y(k)+dn4a*y(l)
        J21=dn1b*x(i)+dn2b*x(j)+dn3b*x(k)+dn4b*x(l)
        J22=dn1b*y(i)+dn2b*y(j)+dn3b*y(k)+dn4b*y(l)
        detJ=J11*J22-J12*J21

        !(B)行列要素(dN/dxdN/dy)
        do i=1,3
            do j=1,8
                bm(i,j)=0.0D0
            end do
        end do
        bm(1,1)= J22*dn1a-J12*dn1b
        bm(1,3)= J22*dn2a-J12*dn2b
        bm(1,5)= J22*dn3a-J12*dn3b
        bm(1,7)= J22*dn4a-J12*dn4b
        bm(2,2)=-J21*dn1a+J11*dn1b
        bm(2,4)=-J21*dn2a+J11*dn2b
        bm(2,6)=-J21*dn3a+J11*dn3b
        bm(2,8)=-J21*dn4a+J11*dn4b
        bm(3,1)=-J21*dn1a+J11*dn1b
        bm(3,2)= J22*dn1a-J12*dn1b
        bm(3,3)=-J21*dn2a+J11*dn2b
        bm(3,4)= J22*dn2a-J12*dn2b
        bm(3,5)=-J21*dn3a+J11*dn3b
        bm(3,6)= J22*dn3a-J12*dn3b
        bm(3,7)=-J21*dn4a+J11*dn4b
        bm(3,8)= J22*dn4a-J12*dn4b
        do i=1,3
            do j=1,8
                bm(i,j)=bm(i,j)/detJ
            end do
        end do
    end subroutine BMAT_PL4


    subroutine IPOINT4(kk,a,b)
        integer :: kk
        real(8) a,b

        !Gauss積分点
        select case (kk)
            case (1)
                a=-0.5773502692D0
                b=-0.5773502692D0
            case (2)
                a= 0.5773502692D0
                b=-0.5773502692D0
            case (3)
                a=0.5773502692D0
                b=0.5773502692D0
            case (4)
                a=-0.5773502692D0
                b= 0.5773502692D0
        end select
    end subroutine IPOINT4


    subroutine SM_PL4(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
        !要素剛性行列
        integer :: NELT,NODT
        integer :: ne,kakom(1:NELT,1:4)
        real(8) :: x(1:NODT),y(1:NODT)
        integer :: NSTRES
        real(8) :: Em(1:NELT),po(1:NELT),t(1:NELT),sm(1:8,1:8)

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

        !配列クリア
        do i=1,8
            do j=1,8
                sm(i,j)=0.0D0
            end do
        end do
        !剛性行列(sm)=(B)T(D,B)*t*detJ
        if(NSTRES==0)t(ne)=1.0D0
        call DMAT_PL(NSTRES,Em(ne),po(ne),d)
        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
            do i=1,3
                do j=1,8
                    str(i,j)=0.0
                    do k=1,3
                        str(i,j)=str(i,j)+d(i,k)*bm(k,j)
                    end do
                end do
            end do
            do i=1,8
                do j=1,8
                    do k=1,3
                        sm(i,j)=sm(i,j)+bm(k,i)*str(k,j)*t(ne)*detJ
                    end do
                end do
            end do
        end do
    end subroutine SM_PL4


    subroutine MV_PL4(ne,kakom,x,y,t,gamma,gkh,gkv,fevec,NELT,NODT)
        !物体力ベクトル
        integer :: NELT,NODT
        integer :: ne,kakom(1:NELT,1:4)
        real(8) :: x(1:NODT),y(1:NODT)
        real(8) :: t(1:NELT),gamma(1:NELT),gkh(1:NELT),gkv(1:NELT),fevec(1:8)

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

        do ii=1,8
            do jj=1,8
                ntn(ii,jj)=0.0D0
            end do
        end do
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        l=kakom(ne,4)
        do kk=1,4
            call IPOINT4(kk,a,b)
            ![dN/da][dN/db]
            dn1a=-0.25D0*(1.0D0-b)
            dn2a= 0.25D0*(1.0D0-b)
            dn3a= 0.25D0*(1.0D0+b)
            dn4a=-0.25D0*(1.0D0+b)
            dn1b=-0.25D0*(1.0D0-a)
            dn2b=-0.25D0*(1.0D0+a)
            dn3b= 0.25D0*(1.0D0+a)
            dn4b= 0.25D0*(1.0D0-a)
            !Jacobiマトリックスとdet(J)
            J11=dn1a*x(i)+dn2a*x(j)+dn3a*x(k)+dn4a*x(l)
            J12=dn1a*y(i)+dn2a*y(j)+dn3a*y(k)+dn4a*y(l)
            J21=dn1b*x(i)+dn2b*x(j)+dn3b*x(k)+dn4b*x(l)
            J22=dn1b*y(i)+dn2b*y(j)+dn3b*y(k)+dn4b*y(l)
            detJ=J11*J22-J12*J21
            do ii=1,2
                do jj=1,8
                    nm(ii,jj)=0.0D0
                end do
            end do
            nm(1,1)=0.25D0*(1.0D0-a)*(1.0D0-b)
            nm(1,3)=0.25D0*(1.0D0+a)*(1.0D0-b)
            nm(1,5)=0.25D0*(1.0D0+a)*(1.0D0+b)
            nm(1,7)=0.25D0*(1.0D0-a)*(1.0D0+b)
            nm(2,2)=nm(1,1)
            nm(2,4)=nm(1,3)
            nm(2,6)=nm(1,5)
            nm(2,8)=nm(1,7)
            do ii=1,8
                do jj=1,8
                    do ll=1,2
                        ntn(ii,jj)=ntn(ii,jj)+nm(ll,ii)*nm(ll,jj)*gamma(ne)*t(ne)*detJ
                    end do
                end do
            end do
        end do
        !物体力{ρg}
        w(1)=gkh(ne)
        w(2)=gkv(ne)
        w(3)=gkh(ne)
        w(4)=gkv(ne)
        w(5)=gkh(ne)
        w(6)=gkv(ne)
        w(7)=gkh(ne)
        w(8)=gkv(ne)
        do i=1,8
            fevec(i)=0.0
            do j=1,8
                fevec(i)=fevec(i)+ntn(i,j)*w(j)
            end do
        end do
    end subroutine MV_PL4


    subroutine CALST_PL4TS(ne,kakom,x,y,deltaT,dist,NSTRES,Em,po,t,alpha,ts,strsave,fevec,noten,NELT,NODT,nt)

        integer :: NELT,NODT,nt
        integer :: ne,kakom(1:NELT,1:4)
        real(8) :: x(1:NODT),y(1:NODT),deltaT(1:NODT)
        real(8) :: dist(1:nt)
        integer :: NSTRES
        real(8) :: Em(1:NELT),po(1:NELT),t(1:NELT),alpha(1:NELT),ts(1:NELT),strsave(1:NELT,1:4,1:3)
        real(8) :: fevec(1:8)
        integer :: noten(1:NELT,1:4)

        !応力および内力ベクトル計算 
        integer :: kk,i,j,ia
        integer :: k1,k2,k3,k4
        real(8) :: detJ,a,b
        real(8) :: d(1:3,1:3),bm(1:3,1:8),epsilon(1:3),eps0(1:3),stress(1:3),wd(1:8)
        real(8) :: Eme,poe
        real(8) :: tem,nm1,nm2,nm3,nm4
        real(8) :: sigx,sigy,tauxy,ps1,ps2,ang,phi

        !要素座標節点変位 
        do i=1,4
            ia=kakom(ne,i)
            wd(2*i-1)=dist(2*ia-1)
            wd(2*i)=dist(2*ia)
        end do
        k1=kakom(ne,1)
        k2=kakom(ne,2)
        k3=kakom(ne,3)
        k4=kakom(ne,4)

        !要素応力
        if(NSTRES==0)t(ne)=1.0D0
        Eme=Em(ne)
        poe=po(ne)
        do kk=1,4
            if(noten(ne,kk)/=0)poe=0.0D0
        end do
        call DMAT_PL(NSTRES,Eme,poe,d)

        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
            !ひずみ計算 {epsilon}=(B){u} 
            do i=1,3
                epsilon(i)=0.0D0
                do j=1,8
                    epsilon(i)=epsilon(i)+bm(i,j)*wd(j)
                end do
            end do
            !積分点温度変化量
            nm1=0.25D0*(1.0D0-a)*(1.0D0-b)
            nm2=0.25D0*(1.0D0+a)*(1.0D0-b)
            nm3=0.25D0*(1.0D0+a)*(1.0D0+b)
            nm4=0.25D0*(1.0D0-a)*(1.0D0+b)
            tem=nm1*deltaT(k1)+nm2*deltaT(k2)+nm3*deltaT(k3)+nm4*deltaT(k4)
            !温度ひずみ計算
            select case(NSTRES)
                case(0) !平面歪み 
                    eps0(1)=tem*(1.0D0+poe)*alpha(ne)
                    eps0(2)=eps0(1)
                    eps0(3)=0.0D0
                case(1) !平面応力 
                    eps0(1)=tem*alpha(ne)
                    eps0(2)=eps0(1)
                    eps0(3)=0.0D0
            end select
            !応力計算 {stress}=(D){epsilon-eps0} 
            do i=1,3
                stress(i)=0.0D0
                do j=1,3
                    stress(i)=stress(i)+d(i,j)*(epsilon(j)-eps0(j))
                end do
            end do
            !現状応力保存
            sigx =stress(1)
            sigy =stress(2)
            tauxy=stress(3)
            call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
            !No-tension判定
            noten(ne,kk)=0
            if(ts(ne)<ps1)noten(ne,kk)=1
            if(ts(ne)<ps2)noten(ne,kk)=2
            !引張成分除去
            phi=-ang/180.0D0*pi
            if(noten(ne,kk)/=0)then
                sigx =sigx -0.5D0*ps1*(1.0D0+cos(2.0D0*phi))
                sigy =sigy -0.5D0*ps1*(1.0D0-cos(2.0D0*phi))
                tauxy=tauxy+0.5D0*ps1*sin(2.0D0*phi)
            end if
            if(noten(ne,kk)==2)then
                    sigx =sigx -0.5D0*ps2*(1.0D0-cos(2.0D0*phi))
                    sigy =sigy -0.5D0*ps2*(1.0D0+cos(2.0D0*phi))
                    tauxy=tauxy-0.5D0*ps2*sin(2.0D0*phi)
            end if
            !作業用応力保存
            strsave(ne,kk,1)=sigx
            strsave(ne,kk,2)=sigy
            strsave(ne,kk,3)=tauxy
        end do

        !内力ベクトル 
        do i=1,8
            fevec(i)=0.0D0
        end do
        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
            do i=1,8
                do j=1,3
                    fevec(i)=fevec(i)+bm(j,i)*strsave(ne,kk,j)*t(ne)*detJ
                end do
            end do
        end do

    end subroutine CALST_PL4TS


    subroutine PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
        !主応力計算
        real(8) :: sigx,sigy,tauxy,ps1,ps2,ang

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


    subroutine ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
        !剛性行列組立(nod：1要素節点数，nhen：1節点の自由度)
        integer :: NELT,NODT
        integer :: ne,nod,nhen
        integer :: kakom(1:NELT,1:nod)
        real(8) :: sm(1:nod*nhen,1:nod*nhen)
        real(8) :: 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 ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !荷重ベクトル組立(nod：1要素節点数，nhen：1節点の自由度)
        integer :: NELT,NODT
        integer :: ne,nod,nhen
        integer :: kakom(1:NELT,1:nod)
        real(8) :: fevec(1:nod*nhen)
        real(8) :: ftvec(1:NODT*nhen)

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

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


    subroutine CHBAND10(n1,m1,array,nt)
        integer :: n1,m1,nt
        real(8) :: 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 :: n1,m1,nt
        real(8) :: array(1:nt,1:nt),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 :: mm,nt
        real(8) :: 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 :: mm,ib,nt
        real(8) :: 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_FEM4nodTS