program f90_FEM4nodASTS
    !２次元no-tension軸対称応力解析
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    integer::i,j,k,l,kk,ne,n,ia,ja
    character::strcom*50                !書出用コメント
    integer::NODT                       !節点総数
    integer::NELT                       !要素総数
    integer::MATEL                      !材料種類数
    integer::KOZ                        !z方向変位既知節点数
    integer::KOR                        !r方向変位既知節点数
    integer::NF                         !載荷節点数
    integer,allocatable::kakom(:,:)     !要素構成接点番号
    integer,allocatable::matno(:)       !材料種別No
    real(8),allocatable::wm(:,:)        !材料物性作業用
    real(8),allocatable::Em(:)          !要素弾性係数
    real(8),allocatable::po(:)          !要素ポアソン比
    real(8),allocatable::alpha(:)       !要素線膨張係数
    real(8),allocatable::ts(:)          !要素引張強度
    real(8),allocatable::z(:)           !節点z座標
    real(8),allocatable::r(:)           !節点r座標
    real(8),allocatable::deltaT(:)      !節点温度変化量(＋：温度上昇)
    integer,allocatable::nokz(:)        !z方向変位既知節点番号
    integer,allocatable::nokr(:)        !r方向変位既知節点番号
    real(8),allocatable::f(:)           !全体外力ベクトル
    integer,allocatable::index(:)       !既知節点変位処理のためのインデックス
    real(8),allocatable::ftvec(:)       !全体荷重ベクトル
    real(8),allocatable::tsm(:,:)       !全体剛性マトリックス
    real(8),allocatable::dis(:)         !全体節点変位増分
    real(8),allocatable::rdis(:)        !既知節点変位変
    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):: sigz,sigr,sigt,tauzr,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*10,time1*20,time2*20,zone*10

    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,',',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,',',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
    !------------------------------------------------------------------------------------------
    !節点数，要素数，材料種類数,z方向変位既知節点数，r方向変位既知節点数，載荷点数
    !------------------------------------------------------------------------------------------
    read(11,*) NODT,NELT,MATEL,KOZ,KOR,NF,IPR
    !------------------------------------------------------------------------------------------
    !配列寸法宣言
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:4))
    allocate(Em(1:NELT))
    allocate(po(1:NELT))
    allocate(alpha(1:NELT))
    allocate(ts(1:NELT))
    allocate(z(1:NODT))
    allocate(r(1:NODT))
    allocate(deltaT(1:NODT))
    allocate(nokz(1:KOZ))
    allocate(nokr(1:KOR))
    allocate(f(1:nt))
    allocate(index(1:nt))
    allocate(ftvec(1:nt))
    allocate(tsm(1:nt,1:nt))
    allocate(dis(1:nt))
    allocate(rdis(1:nt))
    allocate(wdis(1:nt))
    allocate(dist(1:nt))
    allocate(reac(1:nt))
    allocate(strsave(1:NELT,1:4,1:4))
    allocate(noten(1:NELT,1:4))
    allocate(sm(1:nod*nhen,1:nod*nhen))
    allocate(fevec(1:nod*nhen))
    !------------------------------------------------------------------------------------------
    !材料特性（Em,po,alpha,ts）
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3),wm(i,4)
    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
                Em(ne)   =wm(i,1)
                po(ne)   =wm(i,2)
                alpha(ne)=wm(i,3)
                ts(ne)   =wm(i,4)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !節点座標（x,y），温度変化量
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) z(i),r(i),deltaT(i)
    end do
    !------------------------------------------------------------------------------------------
    !変位既知節点番号，既知節点変位量
    !------------------------------------------------------------------------------------------
    do i=1,nt
        rdis(i)=0.0D0
    end do
    if(0<KOZ)then
        do i=1,KOZ
            read(11,*) nokz(i),rdis(2*nokz(i)-1)
        end do
    end if
    if(0<KOR)then
        do i=1,KOR
            read(11,*) nokr(i),rdis(2*nokr(i))
        end do
    end if
    !------------------------------------------------------------------------------------------
    !節点番号,z方向荷重,r方向荷重
    !------------------------------------------------------------------------------------------
    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,*) 'NODT,NELT,MATEL,KOZ,KOR,NF,IPR'
        write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOZ,',',KOR,',',NF,',',IPR
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,*) '*node characteristics'
    write(12,*) 'node,z,r,fz,fr,fix-z,fix-r,rdis-z,rdis-r,deltaT'
    do i=1,NODT
        k=0
        l=0
        if(0<KOZ)then
            do j=1,KOZ
                if(i==nokz(j))k=1
            end do
        end if
        if(0<KOR)then
            do j=1,KOR
                if(i==nokr(j))l=1
            end do
        end if
            write(linebuf,FMT=fmt1) i,z(i),r(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,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),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*((z(k)-z(j))*r(i)+(z(i)-z(k))*r(j)+(z(j)-z(i))*r(k))
        if(elarea1<fact0)stop
        i=kakom(ne,1)
        j=kakom(ne,3)
        k=kakom(ne,4)
        elarea2=0.5*((z(k)-z(j))*r(i)+(z(i)-z(k))*r(j)+(z(j)-z(i))*r(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
        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_AX4(ne,kakom,z,r,Em,po,sm,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
    end do
    !-----------------------------------------------------
    !荷重ベクトルセット
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=f(i)
    end do
    !-----------------------------------------------------
    !変位既知節点自由度処理
    !-----------------------------------------------------
    do i=1,nt
        index(i)=i
    end do
    if(0<KOZ)then
        do i=1,KOZ
            n=2*nokz(i)-1
            index(n)=0
        end do
    end if
    if(0<KOR)then
        do i=1,KOR
            n=2*nokr(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_AX4(ne,kakom,z,r,deltaT,dist,Em,po,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)-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-z,coord-r,dist-z,dist-r,reac-z,reac-r,ftvec-z,ftvec-r'
    do i=1,NODT
        write(linebuf,FMT=fmt3) i,z(i),r(i),dist(2*i-1),dist(2*i),reac(2*i-1),reac(2*i),&
            f(2*i-1)-reac(2*i-1),f(2*i)-reac(2*i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !要素平均応力
    write(12,*) '*stresses'
    write(12,*) 'element,kk,sig-z,sig-r,sig-t,tau-zr,ps1,ps2,ang,noten,matno'
    select case (IPR)
        case(0)
            do ne=1,NELT
                do kk=1,4
                    sigz =strsave(ne,kk,1)
                    sigr =strsave(ne,kk,2)
                    sigt =strsave(ne,kk,3)
                    tauzr=strsave(ne,kk,4)
                    call PST_AX4(sigz,sigr,tauzr,ps1,ps2,ang)
                    write(linebuf,FMT=fmt4) ne,kk,sigz,sigr,sigt,tauzr,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
                sigz =0.0D0
                sigr =0.0D0
                sigt =0.0D0
                tauzr=0.0D0
                do kk=1,4
                    sigz =sigz +strsave(ne,kk,1)
                    sigr =sigr +strsave(ne,kk,2)
                    sigt =sigt +strsave(ne,kk,3)
                    tauzr=tauzr+strsave(ne,kk,4)
                    ntsum=ntsum+noten(ne,kk)
                end do
                kk=0
                sigz =0.25D0*sigz
                sigr =0.25D0*sigr
                sigt =0.25D0*sigt
                tauzr=0.25D0*tauzr
                call PST_AX4(sigz,sigr,tauzr,ps1,ps2,ang)
                write(linebuf,FMT=fmt4) ne,kk,sigz,sigr,sigt,tauzr,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_AX4(Eme,pzr,pot,d)
        !応力-ひずみ関係行列(De)
        real(8),intent(in):: Eme,pzr,pot
        real(8),intent(out)::d(1:4,1:4)
        integer::i,j

        do i=1,4
            do j=1,4
                d(i,j)=0.0D0
            end do
        end do
        d(1,1)=1.0D0-pzr
        d(1,2)=pzr
        d(1,3)=pot
        d(1,4)=0.0D0
        d(2,1)=pzr
        d(2,2)=1.0D0-pzr
        d(2,3)=pot
        d(2,4)=0.0D0
        d(3,1)=pot
        d(3,2)=pot
        d(3,3)=1.0D0-pot
        d(3,4)=0.0D0
        d(4,1)=0.0D0
        d(4,2)=0.0D0
        d(4,3)=0.0D0
        d(4,4)=0.5D0*(1.0D0-2.0D0*pzr)
        do i=1,4
            do j=1,4
                d(i,j)=d(i,j)*Eme/(1.0D0+pzr)/(1.0D0-2.0D0*pzr)
            end do
        end do

    end subroutine DMAT_AX4


    subroutine BMAT_AX4(ne,kakom,z,r,bm,detJ,a,b,NELT,NODT)
        !ひずみ-変位関係行列[B]
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::z(1:NODT),r(1:NODT),a,b
        real(8),intent(out)::bm(1:4,1:8),detJ

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

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

        N1=0.25D0*(1.0D0-a)*(1.0D0-b)
        N2=0.25D0*(1.0D0+a)*(1.0D0-b)
        N3=0.25D0*(1.0D0+a)*(1.0D0+b)
        N4=0.25D0*(1.0D0-a)*(1.0D0+b)
        rr=N1*r(i)+N2*r(j)+N3*r(k)+N4*r(l)

        ![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*z(i)+dn2a*z(j)+dn3a*z(k)+dn4a*z(l)
        J12=dn1a*r(i)+dn2a*r(j)+dn3a*r(k)+dn4a*r(l)
        J21=dn1b*z(i)+dn2b*z(j)+dn3b*z(k)+dn4b*z(l)
        J22=dn1b*r(i)+dn2b*r(j)+dn3b*r(k)+dn4b*r(l)
        detJ=J11*J22-J12*J21

        !(B)行列要素[dN/dx][dN/dy]
        do i=1,4
            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,2)=N1/rr*detJ
        bm(3,4)=N2/rr*detJ
        bm(3,6)=N3/rr*detJ
        bm(3,8)=N4/rr*detJ

        bm(4,1)=-J21*dn1a+J11*dn1b
        bm(4,2)= J22*dn1a-J12*dn1b
        bm(4,3)=-J21*dn2a+J11*dn2b
        bm(4,4)= J22*dn2a-J12*dn2b
        bm(4,5)=-J21*dn3a+J11*dn3b
        bm(4,6)= J22*dn3a-J12*dn3b
        bm(4,7)=-J21*dn4a+J11*dn4b
        bm(4,8)= J22*dn4a-J12*dn4b

        do i=1,4
            do j=1,8
                bm(i,j)=bm(i,j)/detJ
            end do
        end do

    end subroutine BMAT_AX4


    subroutine IPOINT4(kk,a,b)
        integer,intent(in)::kk
        real(8),intent(out)::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_AX4(ne,kakom,z,r,Em,po,sm,NELT,NODT)
        !要素剛性行列
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::z(1:NODT),r(1:NODT),Em(1:NELT),po(1:NELT)
        real(8),intent(out)::sm(1:8,1:8)

        integer:: kk,i,j,k
        integer i1,i2,i3,i4
        real(8):: d(1:4,1:4),bm(1:4,1:8),str(1:4,1:8),detJ,a,b,Eme,poe
        real(8):: N1,N2,N3,N4,rr

        !配列クリア
        do i=1,8
            do j=1,8
                sm(i,j)=0.0D0
            end do
        end do
        !剛性行列[sm]=[B]T[D][B]*r*detJ
        Eme = Em(ne)
        poe = po(ne)
        call DMAT_AX4(Eme,poe,poe,d)

        i1=kakom(ne,1)
        i2=kakom(ne,2)
        i3=kakom(ne,3)
        i4=kakom(ne,4)

        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_AX4(ne,kakom,z,r,bm,detJ,a,b,NELT,NODT)
            do i=1,4
                do j=1,8
                    str(i,j)=0.0D0
                    do k=1,4
                        str(i,j)=str(i,j)+d(i,k)*bm(k,j)
                    end do
                end do
            end do
            N1=0.25D0*(1.0D0-a)*(1.0D0-b)
            N2=0.25D0*(1.0D0+a)*(1.0D0-b)
            N3=0.25D0*(1.0D0+a)*(1.0D0+b)
            N4=0.25D0*(1.0D0-a)*(1.0D0+b)
            rr=N1*r(i1)+N2*r(i2)+N3*r(i3)+N4*r(i4)
            do i=1,8
                do j=1,8
                    do k=1,4
                        sm(i,j)=sm(i,j)+bm(k,i)*str(k,j)*rr*detJ
                    end do
                end do
            end do
        end do

    end subroutine SM_AX4


    subroutine CALST_AX4(ne,kakom,z,r,deltaT,dist,Em,po,alpha,ts,strsave,fevec,noten,NELT,NODT,nt)
        !応力計算
        integer,intent(in):: NELT,NODT,nt
        integer,intent(in):: ne,kakom(1:NELT,1:4)
        real(8),intent(in):: z(1:NODT),r(1:NODT),deltaT(1:NODT),dist(1:nt)
        real(8),intent(in):: Em(1:NELT),po(1:NELT),alpha(1:NELT),ts(1:NELT)
        real(8),intent(out)::fevec(1:8),strsave(1:NELT,1:4,1:4)
        integer,intent(inout)::noten(1:NELT,1:4)

        integer:: kk,i,j,i1,i2,i3,i4,ia
        real(8):: tem,d(1:4,1:4),bm(1:4,1:8),stress(1:4),epsilon(1:4),eps0(1:4),wd(1:8)
        real(8):: sigz,sigr,sigt,tauzr,ps1,ps2,ang,phi,detJ,a,b,Eme,pzr,pot
        real(8):: N1,N2,N3,N4,rr

        !要素座標節点変位増分
        do i=1,4
            ia=kakom(ne,i)
            wd(2*i-1)=dist(2*ia-1)
            wd(2*i)  =dist(2*ia)
        end do

        i1=kakom(ne,1)
        i2=kakom(ne,2)
        i3=kakom(ne,3)
        i4=kakom(ne,4)

        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_AX4(ne,kakom,z,r,bm,detJ,a,b,NELT,NODT)
            !ひずみ計算 {epsilon}=(B){u}
            do i=1,4
                epsilon(i)=0.0D0
                do j=1,8
                    epsilon(i)=epsilon(i)+bm(i,j)*wd(j)
                end do
            end do
            !積分点温度ひずみ
            N1=0.25D0*(1.0D0-a)*(1.0D0-b)
            N2=0.25D0*(1.0D0+a)*(1.0D0-b)
            N3=0.25D0*(1.0D0+a)*(1.0D0+b)
            N4=0.25D0*(1.0D0-a)*(1.0D0+b)
            tem=N1*deltaT(i1)+N2*deltaT(i2)+N3*deltaT(i3)+N4*deltaT(i4)
            eps0(1)=alpha(ne)*tem
            eps0(2)=alpha(ne)*tem
            eps0(3)=alpha(ne)*tem
            eps0(4)=0.0D0
            Eme=Em(ne)
            select case(noten(ne,kk))
                Case(0)
                    pzr=po(ne)
                    pot=po(ne)
                Case(1)
                    pzr=0.0D0
                    pot=po(ne)
                Case(2)
                    pzr=0.0D0
                    pot=po(ne)
                Case(3)
                    pzr=po(ne)
                    pot=0.0D0
                Case(4)
                    pzr=0.0D0
                    pot=0.0D0
                Case(5)
                    pzr=0.0D0
                    pot=0.0D0
            end select
            call DMAT_AX4(Eme,pzr,pot,d)
            !応力計算 {stress}=(D){epsilon}
            do i=1,4
                stress(i)=0.0D0
                do j=1,4
                    stress(i)=stress(i)+d(i,j)*(epsilon(j)-eps0(j))
                end do
             end do
            !現状応力保存
            sigz =stress(1)
            sigr =stress(2)
            sigt =stress(3)
            tauzr=stress(4)
            Call PST_AX4(sigz,sigr,tauzr,ps1,ps2,ang)
            !No-tension判定
            noten(ne,kk)=0
            if((ts(ne)<ps1   .and.ps2   <ts(ne)).and.sigt  <ts(ne))noten(ne,kk)=1
            if((ts(ne)<ps1   .and.ts(ne)<ps2   ).and.sigt  <ts(ne))noten(ne,kk)=2
            if((ps1   <ts(ne).and.ps2   <ts(ne)).and.ts(ne)<sigt  )noten(ne,kk)=3
            if((ts(ne)<ps1   .and.ps2   <ts(ne)).and.ts(ne)<sigt  )noten(ne,kk)=4
            if((ts(ne)<ps1   .and.ts(ne)<ps2   ).and.ts(ne)<sigt  )noten(ne,kk)=5
            !引張成分除去
            phi=-ang/180.0D0*pi
            if(noten(ne,kk)==1.or.noten(ne,kk)==2.or.noten(ne,kk)==4)then
                sigz =sigz -0.5D0*ps1*(1.0D0+cos(2.0D0*phi))
                sigr =sigr -0.5D0*ps1*(1.0D0-cos(2.0D0*phi))
                tauzr=tauzr+0.5D0*ps1*sin(2.0D0*phi)
            end if
            if(noten(ne,kk)==2)then
                sigz =sigz -0.5D0*ps2*(1.0D0-cos(2.0D0*phi))
                sigr =sigr -0.5D0*ps2*(1.0D0+cos(2.0D0*phi))
                tauzr=tauzr-0.5D0*ps2*sin(2.0D0*phi)
            end if
            if(noten(ne,kk)==3.or.noten(ne,kk)==4.or.noten(ne,kk)==5)then
                sigt=0.0D0
            end if
            !作業用応力保存
            strsave(ne,kk,1)=sigz
            strsave(ne,kk,2)=sigr
            strsave(ne,kk,3)=sigt
            strsave(ne,kk,4)=tauzr
        end do

        !内力ベクトル
        do i=1,8
            fevec(i)=0.0D0
        end do
        do kk=1,4
            call IPOINT4(kk,a,b)
            call BMAT_AX4(ne,kakom,z,r,bm,detJ,a,b,NELT,NODT)
            N1=0.25D0*(1.0D0-a)*(1.0D0-b)
            N2=0.25D0*(1.0D0+a)*(1.0D0-b)
            N3=0.25D0*(1.0D0+a)*(1.0D0+b)
            N4=0.25D0*(1.0D0-a)*(1.0D0+b)
            rr=N1*r(i1)+N2*r(i2)+N3*r(i3)+N4*r(i4)
            do i=1,8
                do j=1,4
                    fevec(i)=fevec(i)+bm(j,i)*strsave(ne,kk,j)*rr*detJ
                end do
            end do
        end do

    end subroutine CALST_AX4


    subroutine PST_AX4(sigz,sigr,tauzr,ps1,ps2,ang)
        !主応力計算
        real(8),intent(in):: sigz,sigr,tauzr
        real(8),intent(out)::ps1,ps2,ang

        ps1=0.5D0*(sigz+sigr)+sqrt(0.25D0*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr)
        ps2=0.5D0*(sigz+sigr)-sqrt(0.25D0*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr)
        if(sigz==sigr)then
            if(tauzr> 0.0D0)ang= 45.0D0
            if(tauzr< 0.0D0)ang=135.0D0
            if(tauzr==0.0D0)ang=  0.0D0
        else
            ang=0.5D0*atan(2.0D0*tauzr/(sigz-sigr))
            ang=180.0D0/pi*ang
            if((sigz>sigr).and.(tauzr>=0.0D0))ang=ang
            if((sigz>sigr).and.(tauzr< 0.0D0))ang=ang+180.0D0
            if(sigz<sigr)                     ang=ang+ 90.0D0
        end if

    end subroutine PST_AX4


    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_FEM4nodASTS