module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
end module defpi

program f90_FEM_ASNT
    !****************************************************************
    !Axisymmetric Analysis with function of no-tension analysis
    !****************************************************************
    use defpi
    implicit none
    integer::i,j,k,l,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::KOZ                        !Number of restricted nodes in z-direction (axial direction)
    integer::KOR                        !Number of restricted nodes in r-direction (radius direction)
    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::po(:)          !Poisson's ratio of element
    real(8),allocatable::gamma(:)       !Unit weight of element
    real(8),allocatable::gkz(:)         !Coefficient of inertia force (ratio to gravity acceleration)
    real(8),allocatable::alpha(:)       !Coefficient of thermal expansion of element
    real(8),allocatable::ts(:)          !Tensile strength of element
    real(8),allocatable::z(:)           !z-coordinate of node
    real(8),allocatable::r(:)           !r-coordinate of node
    real(8),allocatable::deltaT(:)      !Temperature change of node (+: temperature rise)
    integer,allocatable::nokz(:)        !Restricted node number in z-direction
    integer,allocatable::nokr(:)        !Restricted node number in r-direction
    real(8),allocatable::f(:)           !External force vector
    real(8),allocatable::fmass(:)       !Inertia force vector
    integer,allocatable::mindex(:)      !Index for treatment of restricted nodes
    real(8),allocatable::ftvec(:)       !Nodal force vector
    real(8),allocatable::tsm(:,:)       !Stiffness matrix
    real(8),allocatable::dis(:)         !Increment of nodal displacement
    real(8),allocatable::rdis(:)        !Displacement of restricted node
    real(8),allocatable::wdis(:)        !work for displacement
    real(8),allocatable::dist(:)        !Cumulative displacement of node
    real(8),allocatable::reac(:)        !Internal force vector
    real(8),allocatable::strsave(:,:,:) !Memory of stresses
    integer,allocatable::noten(:,:)     !Memory of stress state
    real(8),allocatable::sm(:,:)        !Stiffness matrix of element
    real(8),allocatable::fevec(:)       !Nodal force vector of element
    integer::nt                         !Number of degrees of freedom of all nodes (NODT x nhen)
    integer::nod=4                      !Number of nodes per element
    integer::nhen=2                     !Number of degrees of freedom per node
    integer::mm                         !Number of degrees of freedom of reduced FE equation
    integer::ib                         !Band width of redused matrix
    real(8)::elarea1                    !Half area of element
    real(8)::elarea2                    !Half area of element
    real(8)::fact0=1.0D-30              !zero judgement
    real(8):: sigz,sigr,sigt,tauzr,ps1,ps2,ang !Stresses
    integer:: nnn                       !Increment for iterative calculation
    integer:: maxita=2000               !Limit value of iterative calculation
    integer:: ntsum                     !Number of no-tension state gauss points
    integer:: icount                    !Number of converged degrees of freedom
    real(8):: dtest,ftest               !Sum total of increment of nodal displacement, um total of unbalanced force
    integer:: IPR                       !Output format of stresses (0: all Gauss points, 1: average of element)
    character::fnameR*50,fnameW*50      !Input file name, output file name

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

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

    !To input file names
    call getarg(1,fnameR) !Input file name
    call getarg(2,fnameW) !Output file name

    !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,KOZ,KOR,NF,IPR
    !------------------------------------------------------------------------------------------
    !To Declare array size
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:6))
    allocate(Em(1:NELT))
    allocate(po(1:NELT))
    allocate(gamma(1:NELT))
    allocate(gkz(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(fmass(1:nt))
    allocate(mindex(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))
    !------------------------------------------------------------------------------------------
    !Material properties (Em,po,gamma,gkz,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)
    end do
    !------------------------------------------------------------------------------------------
    !Element-nodes relationship, material set No.
    !------------------------------------------------------------------------------------------
    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)
                gamma(ne)=wm(i,3)
                gkz(ne)  =wm(i,4)
                alpha(ne)=wm(i,5)
                ts(ne)   =wm(i,6)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !Coordinate (z,r), temperature change of node
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) z(i),r(i),deltaT(i)
    end do
    !------------------------------------------------------------------------------------------
    !Restricted node No., displacement of restricted node
    !------------------------------------------------------------------------------------------
    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
    !------------------------------------------------------------------------------------------
    !Node No., load in z-direction, load in r-direction
    !------------------------------------------------------------------------------------------
    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)

    !*****************************************************
    !To print input data
    !*****************************************************
    open(12, file=fnameW, status='replace')
    write(12,'(a)') strcom
    write(12,'(a)') '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,'(a)') '*node characteristics'
    write(12,'(a)') '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,'(a)') '*element characteristics'
    write(12,'(a)') 'element,node-1,node-2,node-3,node-4,E,po,gamma,gkz,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),gamma(ne),gkz(ne),alpha(ne),ts(ne),matno(ne)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)

    !to verify the area of element & to stop if finding the element with negative area
    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

    !*****************************************************
    !Main calculation
    !*****************************************************
    !-----------------------------------------------------
    !Initialization of array for stresses
    !-----------------------------------------------------
    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
    !-----------------------------------------------------
    !Initialization of array & vector
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=0.0D0
        fmass(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
    !-----------------------------------------------------
    !To Assemble stiffness matrix
    !-----------------------------------------------------
    do ne=1,NELT
        !Element stiffness matrix
        call SM_AX4(ne,kakom,z,r,Em,po,sm,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
        call MV_AX4(ne,kakom,z,r,gamma,gkz,fevec,NELT,NODT)
        call ASVEC(ne,nod,nhen,kakom,fevec,fmass,NELT,NODT)
    end do
    !-----------------------------------------------------
    !To set load vector
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=fmass(i)+f(i)
    end do
    !-----------------------------------------------------
    !Treatment of restricted nodes
    !-----------------------------------------------------
    do i=1,nt
        mindex(i)=i
    end do
    if(0<KOZ)then
        do i=1,KOZ
            n=2*nokz(i)-1
            mindex(n)=0
        end do
    end if
    if(0<KOR)then
        do i=1,KOR
            n=2*nokr(i)
            mindex(n)=0
        end do
    end if
    mm=0
    do i=1,nt
        if(mindex(i)>0)then
            mm=mm+1
            mindex(mm)=i
        end if
    end do
    do i=1,mm
        ia=mindex(i)
        ftvec(i)=ftvec(ia)
        do j=1,mm
            ja=mindex(j)
            tsm(i,j)=tsm(ia,ja)
        end do
    end do
    !-----------------------------------------------------
    !Calculation of band width & memory of band matrix
    !-----------------------------------------------------
    ib=IBAND(mm,tsm,fact0,nt)
    call BAND(mm,ib,tsm,nt)
    !-----------------------------------------------------
    !Checking of diagonal element of band matrix
    !-----------------------------------------------------
    do i=1,mm
        if(abs(tsm(i,1))<fact0)stop
    end do
    !-----------------------------------------------------
    !Cholesky method
    !-----------------------------------------------------
    call CHBAND10(mm,ib,tsm,nt) !Triangle matrix


    !****************************************************
    !Iterative calculation for No-tension analysis
    !****************************************************
    nnn=0
    do while(nnn<=maxita)
        nnn=nnn+1
        !-----------------------------------------------------
        !Cholesky method
        !-----------------------------------------------------
        do i=1,mm
            reac(i)=ftvec(i)
        end do
        call CHBAND11(mm,ib,tsm,reac,nt) !forward elimination & backward substitution
        !-----------------------------------------------------
        !Calculation of increment of displacement
        !-----------------------------------------------------
        do i=1,nt
            dis(i)=0.0
            dist(i)=0.0
        end do
        do i=1,mm
            ia=mindex(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
        !-----------------------------------------------------
        !Calculation of stresses & assembling internal force
        !-----------------------------------------------------
        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) !To assemble Internal force vector
        end do
        !-----------------------------------------------------
        !calculation of unbalanced force vector
        !-----------------------------------------------------
        do i=1,nt
            ftvec(i)=fmass(i)+f(i)-reac(i)
        end do
        !-----------------------------------------------------
        !Treatment of restricted nodes for unbalanced force vector
        !-----------------------------------------------------
        do i=1,mm
            ia=mindex(i)
            ftvec(i)=ftvec(ia)
            dis(i)=dis(ia)
            dist(i)=dist(ia)
        end do
        !-----------------------------------------------------
        !Judgement of convergence
        !-----------------------------------------------------
        !Number of No-tension state Gauss points
        ntsum = 0
        do ne=1,NELT
            do kk=1,4
                ntsum=ntsum+abs(noten(ne,kk))
            end do
        end do
        !Sum total of increment of displacement & unbalanced force
        icount=0
        dtest=0.0
        ftest=0.0
        do i=1,mm
            dtest=dtest+abs(dis(i))   !Sum total of increment of displacement
            ftest=ftest+abs(ftvec(i)) !Sum total of unbalanced force
            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

    ! Arrangement of cumulative displacement 
    do i=1,nt
        dis(i)=0.0
    end do
    do i=1,mm
        ia=mindex(i)
        dis(ia)=dist(i)
    end do
    do i=1,nt
        dist(i)=dis(i)+rdis(i)
    end do

    !*****************************************************
    !To print out of results
    !*****************************************************
    open(12, file=fnameW, access='append')
    !Displacements of nodes
    write(12,'(a)') '*displacements and forces'
    write(12,'(a)') '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
    !Average stresses of elements
    write(12,'(a)') '*stresses'
    write(12,'(a)') '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

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

stop

contains

    subroutine DMAT_AX4(Eme,pzr,pot,d)
        !***************************************************
        !Stress-Strain relationship [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)
        !***************************************************
        !Strain-displacement relationship [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 matrix & 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

        !Element of (B)matrix [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)
        !***************************************************
        !Gauss points for integration of element area
        !***************************************************
        integer,intent(in)::kk
        real(8),intent(out)::a,b

        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)
        !***************************************************
        !Element stiffness matrix
        !***************************************************
        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

        !Initialization of array
        do i=1,8
            do j=1,8
                sm(i,j)=0.0D0
            end do
        end do
        !Stiffness matrix [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)
        !***************************************************
        !Calculation of stresses & Internal force vector
        !***************************************************
        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)
            !Calculation of strains {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
            !Strains due to temperature change at Gauss points
            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)
            !Calculation of stresses {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
            !Memory of stresses (temporary)
            sigz =stress(1)
            sigr =stress(2)
            sigt =stress(3)
            tauzr=stress(4)
            Call PST_AX4(sigz,sigr,tauzr,ps1,ps2,ang)
            !Judgement of No-tension state
            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
            !Elimination of tensile component
            phi=-ang/180.0D0*pi
            if(noten(ne,kk)==1.or.noten(ne,kk)==2.or.noten(ne,kk)==4.or.noten(ne,kk)==5)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.or.noten(ne,kk)==5)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
            !Memory of stresses
            strsave(ne,kk,1)=sigz
            strsave(ne,kk,2)=sigr
            strsave(ne,kk,3)=sigt
            strsave(ne,kk,4)=tauzr
        end do

        !Calculation of internal force vector
        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)
        !***************************************************
        !Calculation of principal stress
        !***************************************************
        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 MV_AX4(ne,kakom,z,r,gamma,gkz,fevec,NELT,NODT)
        !***************************************************
        !Inertia force in z-direction
        !***************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::z(1:NODT),r(1:NODT),gamma(1:NELT),gkz(1:NELT)
        real(8),intent(out)::fevec(1:8)

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

        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)
            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 matrix & 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
            do ii=1,2
                do jj=1,8
                    nm(ii,jj)=0.0D0
                end do
            end do
            nm(1,1)=N1
            nm(1,3)=N2
            nm(1,5)=N3
            nm(1,7)=N4
            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)*rr*detJ
                    end do
                end do
            end do
        end do
        !Body force
        w(1)=gkz(ne)
        w(2)=0.0D0
        w(3)=gkz(ne)
        w(4)=0.0D0
        w(5)=gkz(ne)
        w(6)=0.0D0
        w(7)=gkz(ne)
        w(8)=0.0D0
        do i=1,8
            fevec(i)=0.0D0
            do j=1,8
                fevec(i)=fevec(i)+ntn(i,j)*w(j)
            end do
        end do
    end subroutine MV_AX4


    subroutine ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !***************************************************
        !Assembly of vector 
        !***************************************************
        !(nod: number of nodes for 1 element, nhen: degree of freedom of 1 node)
        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
        !***************************************************
        !(nod: number of nodes for 1 element, nhen: degree of freedom of 1 node)
        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)
        !***************************************************
        !Cholesky method (1)
        !Triangle matric
        !***************************************************
        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)
        !***************************************************
        !Cholesky method (2)
        !Forward elimination & backward substitution
        !***************************************************
        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)
        !***************************************************
        !Calculation of band width
        !***************************************************
        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)
        !***************************************************
        !Memory of band matrix
        !***************************************************
        integer,intent(in)::mm,ib,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        integer :: i,j,k,mn
        real(8) :: work(1:nt)

        do i=2,mm
            if(i+ib-1<mm)then
                mn=ib
            else
                mn=mm-i+1
            end if
            do j=1,ib
                work(j)=0.0D0
            end do
            do k=1,mn
                work(k)=array(i,i+k-1)
            end do
            do j=1,ib
                array(i,j)=work(j)
            end do
        end do
    end subroutine BAND


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


    real(4) function DIFT()
        integer::count,count_rate,count_max
        call system_clock(count,count_rate,count_max)
        DIFT=real(count)/real(count_rate)
    end function DIFT


    character(len=50) function CALLDATE()
        character(len=8)::date
        character(len=10)::time
        character(len=5)::zone
        integer::values(1:8)
        integer::i
        character(len=4)::ye
        character(len=4)::da
        character(len=10)::mo
        character(len=5)::ti
        character(len=10)::str
        character(len=50)::datime

        call DATE_AND_TIME(date,time,zone,values)
        ye=date(1:4)
        mo=date(5:6)
        da=date(7:8)
        ti=time(1:2)//':'//time(3:4)

        if(trim(adjustl(mo))=='01')mo='January'
        if(trim(adjustl(mo))=='02')mo='February'
        if(trim(adjustl(mo))=='03')mo='March'
        if(trim(adjustl(mo))=='04')mo='April'
        if(trim(adjustl(mo))=='05')mo='May'
        if(trim(adjustl(mo))=='06')mo='June'
        if(trim(adjustl(mo))=='07')mo='July'
        if(trim(adjustl(mo))=='08')mo='August'
        if(trim(adjustl(mo))=='09')mo='September'
        if(trim(adjustl(mo))=='10')mo='October'
        if(trim(adjustl(mo))=='11')mo='November'
        if(trim(adjustl(mo))=='12')mo='December'

        CALLDATE=trim(adjustl(da))//' '//trim(adjustl(mo))//' '//trim(adjustl(ye))//' at '//trim(adjustl(ti))
    end function CALLDATE


end program f90_FEM_ASNT