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

program f90_fem_plnt
    !*****************************************************************
    !2D Stress Analysis with the function of no-tension analysis
    !*****************************************************************
    use defpi
    implicit none
    integer::i,j,k,l,kk,ne,n,ia,ja
    character::strcom*50                !Comment
    integer::NSTRES                     !Stress state (0: plane strain, 1: plane stress)
    integer::NODT                       !Number of nodes
    integer::NELT                       !Number of elements
    integer::MATEL                      !Number of material sets
    integer::KOX                        !Number of restricted nodes in x-direction
    integer::KOY                        !Number of restricted nodes in y-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::t(:)           !Thickness of element
    real(8),allocatable::gamma(:)       !Unit weight of element
    real(8),allocatable::gkh(:)         !Horizontal acceleration of element in the ratio to 'g'
    real(8),allocatable::gkv(:)         !Vertical acceleration of element in the ratio of 'g'
    real(8),allocatable::alpha(:)       !Coefficient of thermal expantion of element
    real(8),allocatable::ts(:)          !Tensile strength of element
    real(8),allocatable::x(:)           !x-coordinate of node
    real(8),allocatable::y(:)           !y-coordinate of node
    real(8),allocatable::deltaT(:)      !Temperature change of node (+: temperature rize)
    integer,allocatable::nokx(:)        !Restricted node number in x-direction
    integer,allocatable::noky(:)        !Restricted node number in y-direction
    real(8),allocatable::f(:)           !External load vector
    integer,allocatable::mindex(:)      !Index for treatment of restricted nodes
    real(8),allocatable::ftvec(:)       !Total load vector
    real(8),allocatable::fmass(:)       !Body force vector
    real(8),allocatable::rdis(:)        !Displacement of restricted nodes
    real(8),allocatable::tsm(:,:)       !Stiffness matrix
    real(8),allocatable::dis(:)         !Displacement increment of nodes
    real(8),allocatable::wdis(:)        !work for displacement
    real(8),allocatable::dist(:)        !Cumulated displacement of nodes
    real(8),allocatable::reac(:)        !Internal force vector & work
    real(8),allocatable::strsave(:,:,:) !Memory for stresses
    integer,allocatable::noten(:,:)     !Memory for stress state (no-tension or not)
    real(8),allocatable::sm(:,:)        !Element stiffness matrix
    real(8),allocatable::fevec(:)       !Element load vector
    integer::nt                         !Number of degrees of freedom of all nodes (NODT x nhen)
    integer::nod                        !Number of nodes per element
    integer::nhen=2                     !Number of degree of freedom per node
    integer::mm                         !Number of degrees of freedom of reduced FE equation
    integer::ib                         !Band width of reduced matrix
    integer::kkmax                      !Number of integral points
    real(8)::elarea                     !Area of element
    real(8)::fact0=1.0D-30              !0 judgement
    real(8):: sigx,sigy,tauxy,ps1,ps2,ang !Stresses
    real(8)::xg,yg                      !Typical coordinate for 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 displacement, sum total of unbalanced force
    integer:: IPR                       !Flag for print of stresses (0: all Gauss points, 1: element mean stresses)
    character::fnameR*50,fnameW*50      !Input file name & output file name

    character :: linebuf*1000
    character :: fmt1*200,fmt2_3*200,fmt2_4*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_3="(i5,3(',',i5),8(',',e15.7),',',i5)"                            !Element(3 nodes)
    fmt2_4="(i5,4(',',i5),8(',',e15.7),',',i5)"                            !Element(4 nodes)
    fmt3="(i5,2(',',e15.7),2(',',e15.7),2(',',e15.7),2(',',e15.7))"        !Displacement & reaction
    fmt4="(i5,',',i5,2(',',e15.7),3(',',e15.7),3(',',e15.7),2(',',i5))"    !Stresses

    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,*) nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR
    if(nod==3)IPR=0
    !------------------------------------------------------------------------------------------
    !To declare array size
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    if(nod==3)kkmax=1
    if(nod==4)kkmax=4
    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(mindex(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:kkmax,1:3))
    allocate(noten(1:NELT,1:kkmax))
    allocate(sm(1:nod*nhen,1:nod*nhen))
    allocate(fevec(1:nod*nhen))
    !------------------------------------------------------------------------------------------
    !Material properties (t,Em,po,gamma,gkh,gkv,alpha,ts)
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) (wm(i,j),j=1,8)
    end do
    !------------------------------------------------------------------------------------------
    !Element-nodes relationship, material set number
    !------------------------------------------------------------------------------------------
    do ne=1,NELT
        read(11,*) (kakom(ne,j),j=1,nod),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)
                if(NSTRES==0)t(ne)=1.0D0
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !(x,y) coordinates, temperature change
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),y(i),deltaT(i)
    end do
    !------------------------------------------------------------------------------------------
    !Restricted node number, displacement
    !------------------------------------------------------------------------------------------
    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
    !------------------------------------------------------------------------------------------
    !Loaded node number, load in x-direction, load in y-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)') 'nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR'
        write(linebuf,*) nod,',',NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',NF,',',NSTRES,',',IPR
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,'(a)') '*node characteristics'
    write(12,'(a)') '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,'(a)') '*element characteristics'

    select case(nod)
        case(3)
            write(12,'(a)') 'element,node-1,node-2,node-3,E,po,t,gamma,kh,kv,alpha,ts,matno'
            do ne=1,NELT
                write(linebuf,FMT=fmt2_3) ne,kakom(ne,1),kakom(ne,2),kakom(ne,3), &
                    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
        case(4)
            write(12,'(a)') '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_4) 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
    end select
    !To verify element area & To stop if negative area
    do ne=1,NELT
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        elarea=0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        if(elarea<fact0)then
            write(12,'(a)') '***program stop due to area<0'
            close(12)
            stop
        end if
        if(nod==4)then
            i=kakom(ne,1)
            j=kakom(ne,3)
            k=kakom(ne,4)
            elarea=0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
            if(elarea<fact0)then
                write(12,'(a)') '***program stop due to area<0'
                close(12)
                stop
            end if
        end if
    end do
    close(12)
    !*****************************************************
    !Main calculation
    !*****************************************************
    !-----------------------------------------------------
    !Initialization of array for stress saving
    !-----------------------------------------------------
    do ne=1,NELT
        do kk=1,kkmax
            do i=1,3
                strsave(ne,kk,i)=0.0D0
            end do
            noten(ne,kk)=0
        end do
    end do
    !-----------------------------------------------------
    !Initialization of array and vector
    !-----------------------------------------------------
    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
    !-----------------------------------------------------
    !Assembly of stiffness matrix & load vector
    !-----------------------------------------------------
    select case(nod)
        case(3)
            do ne=1,NELT
                !Stiffness matrix
                call SM_PL3(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
                call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
                !Body force vector
                call MV_PL3(ne,kakom,x,y,t,gamma,gkh,gkv,fevec,NELT,NODT)
                call ASVEC(ne,nod,nhen,kakom,fevec,fmass,NELT,NODT)
            end do
        case(4)
            do ne=1,NELT
                !Stiffness matrix
                call SM_PL4(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
                call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
                !Body force vector
                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
    end select
    !-----------------------------------------------------
    !To set load vector
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=f(i)+fmass(i)
    end do
    !-----------------------------------------------------
    !Treatment of restricted nodes
    !-----------------------------------------------------
    do i=1,nt
        mindex(i)=i
    end do
    if(0<KOX)then
        do i=1,KOX
            n=2*nokx(i)-1
            mindex(n)=0
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            n=2*noky(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)
    !-----------------------------------------------------
    !To verify diagonal element & to stop if zero value
    !-----------------------------------------------------
    do i=1,mm
        if(abs(tsm(i,1))<fact0)stop
    end do
    !-----------------------------------------------------
    !Cholesky method
    !-----------------------------------------------------
    call CHBAND10(mm,ib,tsm,nt)!Triangle matrix


    !****************************************************
    !Iterative process 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 displacement increment
        !-----------------------------------------------------
        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 & assembly of internal force vector
        !-----------------------------------------------------
        select case(nod)
        case(3)
            do ne=1,NELT
                call CALST_PL3TS(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)!Assembly of internal force vector
            end do
        case(4)
            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)!Assembly of internal force vector
            end do
        end select
        !-----------------------------------------------------
        !Calculation of unbalanced force vector
        !-----------------------------------------------------
        do i=1,nt
            ftvec(i)=f(i)+fmass(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 Gauss points with no-tension state
        ntsum = 0
        do ne=1,NELT
            do kk=1,kkmax
                ntsum=ntsum+abs(noten(ne,kk))
            end do
        end do
        !Sum total of displacement increment and unbalanced force
        icount=0
        dtest=0.0
        ftest=0.0
        do i=1,mm
            dtest=dtest+abs(dis(i))   !Sum total of displacement increment
            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 cumulated displacements
    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 the results
    !*****************************************************
    open(12, file=fnameW, access='append')
    !Displacement of nodes
    write(12,'(a)') '*displacements and forces'
    write(12,'(a)') '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
    !Stresses of elements
    write(12,'(a)') '*stresses'
    write(12,'(a)') 'element,kk,coord-x,coord-y,sig-x,sig-y,tau-xy,ps1,ps2,ang,noten,matno'
    select case (IPR)
        case(0)
            do ne=1,NELT
                do kk=1,kkmax
                    sigx=strsave(ne,kk,1)
                    sigy=strsave(ne,kk,2)
                    tauxy=strsave(ne,kk,3)
                    call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
                    if(nod==3)call XGYG3M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
                    if(nod==4)call XGYG4I(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
                    write(linebuf,FMT=fmt4) ne,kk,xg,yg,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,kkmax
                    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=1
                sigx=0.25D0*sigx
                sigy=0.25D0*sigy
                tauxy=0.25D0*tauxy
                call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
                call XGYG4M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
                write(linebuf,FMT=fmt4) ne,kk,xg,yg,sigx,sigy,tauxy,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()
    !Output of Summary
    write(12,'("#,Summary")')
    write(12,'("#,NELT=",i0," NODT=",i0," nt=",i0," mm=",i0," ib=",i0)'),NELT,NODT,nt,mm,ib
    write(12,'("#,Convergence nnn=",i0," icount=",i0," dtest=",E10.4," ftest=",E10.4)'),nnn,icount,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 XGYG3M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
        !*****************************************************************
        !Coordinates of centroid of 3 node element
        !*****************************************************************
        integer,intent(in)::ne,kk,NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::xg,yg

        integer::i,j,k

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        xg=(x(i)+x(j)+x(k))/3.0D0
        yg=(y(i)+y(j)+y(k))/3.0D0
    end subroutine XGYG3M


    subroutine XGYG4M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
        !*****************************************************************
        !Coordinates of centroid of 4 node element
        !*****************************************************************
        integer,intent(in)::ne,kk,NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::xg,yg

        integer::i,j,k,l

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        l=kakom(ne,4)
        xg=0.25D0*(x(i)+x(j)+x(k)+x(l))
        yg=0.25D0*(y(i)+y(j)+y(k)+y(l))
    end subroutine XGYG4M


    subroutine XGYG4I(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
        !*****************************************************************
        !Coordinates of Gauss points of 4 node element
        !*****************************************************************
        integer,intent(in)::ne,kk,NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::xg,yg

        integer::i,j,k,l
        real(8)::a,b,N1,N2,N3,N4

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        l=kakom(ne,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)
        xg=N1*x(i)+N2*x(j)+N3*x(k)+N4*x(l)
        yg=N1*y(i)+N2*y(j)+N3*y(k)+N4*y(l)
    end subroutine XGYG4I


    subroutine DMAT_PL(NSTRES,Eme,poe,d)
        !*****************************************************************
        !Stress-strain relationship [De]
        !*****************************************************************
        integer,intent(in)::NSTRES
        real(8),intent(in)::Eme,poe
        real(8),intent(out)::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) !Plane strain
                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) !Plane stress
                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


    real(8) function CALA_PL3(ne,kakom,x,y,NELT,NODT)
        !*****************************************************************
        !Area of 3 node element
        !*****************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        integer::i,j,k

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        CALA_PL3=0.5D0*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))

    end function CALA_PL3


    subroutine BMAT_PL3(ne,kakom,x,y,b,NELT,NODT)
        !*****************************************************************
        !Displacement-strain relationship [B] for 3 node element
        !*****************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::b(1:3,1:6)

        integer::i,j,k
        real(8)::xi,xj,xk,yi,yj,yk,area

        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        xi=x(i)
        xj=x(j)
        xk=x(k)
        yi=y(i)
        yj=y(j)
        yk=y(k)

        do i=1,3
            do j=1,6
                b(i,j)=0.0D0
            end do
        end do
        b(1,1)=yj-yk
        b(1,3)=yk-yi
        b(1,5)=yi-yj
        b(2,2)=xk-xj
        b(2,4)=xi-xk
        b(2,6)=xj-xi
        b(3,1)=b(2,2)
        b(3,2)=b(1,1)
        b(3,3)=b(2,4)
        b(3,4)=b(1,3)
        b(3,5)=b(2,6)
        b(3,6)=b(1,5)
        do i=1,3
            do j=1,6
                b(i,j)=b(i,j)*0.5D0/area
            end do
        end do

    end subroutine BMAT_PL3


    subroutine SM_PL3(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
        !*****************************************************************
        !Element stiffness matrix for 3 node element
        !*****************************************************************
        integer,intent(in)::ne,NELT,NODT,NSTRES
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT)
        real(8),intent(out)::sm(1:6,1:6)

        integer::i,j,k
        real(8)::b(1:3,1:6),d(1:3,1:3),str(1:3,1:6)
        real(8)::area,Eme,poe

        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)
        call BMAT_PL3(ne,kakom,x,y,b,NELT,NODT)
        Eme=Em(ne)
        poe=po(ne)
        call DMAT_PL(NSTRES,Eme,poe,d)
        !Stiffness matrix [sm]=A*t*[B]T[D][B]
        do i=1,3
            do j=1,6
                str(i,j)=0.0D0
                do k=1,3
                    str(i,j)=str(i,j)+d(i,k)*b(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                sm(i,j)=0.0D0
                do k=1,3
                    sm(i,j)=sm(i,j)+b(k,i)*str(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                sm(i,j)=sm(i,j)*area*t(ne)
            end do
        end do

    end subroutine SM_PL3


    subroutine MV_PL3(ne,kakom,x,y,t,gamma,gkh,gkv,fevec,NELT,NODT)
        !*****************************************************************
        !Body force vector for 3 node element
        !*****************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::t(1:NELT),gamma(1:NELT),gkh(1:NELT),gkv(1:NELT)
        real(8),intent(out)::fevec(1:6)

        integer::i,j,k
        real(8)::mmat(1:6,1:6),w(1:6),area

        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)

        mmat(1,1)=0.50D0
        mmat(1,2)=0.00D0
        mmat(1,3)=0.25D0
        mmat(1,4)=0.00D0
        mmat(1,5)=0.25D0
        mmat(1,6)=0.00D0
        mmat(2,1)=0.00D0
        mmat(2,2)=0.50D0
        mmat(2,3)=0.00D0
        mmat(2,4)=0.25D0
        mmat(2,5)=0.00D0
        mmat(2,6)=0.25D0
        mmat(3,1)=0.25D0
        mmat(3,2)=0.00D0
        mmat(3,3)=0.50D0
        mmat(3,4)=0.00D0
        mmat(3,5)=0.25D0
        mmat(3,6)=0.00D0
        mmat(4,1)=0.00D0
        mmat(4,2)=0.25D0
        mmat(4,3)=0.00D0
        mmat(4,4)=0.50D0
        mmat(4,5)=0.00D0
        mmat(4,6)=0.25D0
        mmat(5,1)=0.25D0
        mmat(5,2)=0.00D0
        mmat(5,3)=0.25D0
        mmat(5,4)=0.00D0
        mmat(5,5)=0.50D0
        mmat(5,6)=0.00D0
        mmat(6,1)=0.00D0
        mmat(6,2)=0.25D0
        mmat(6,3)=0.00D0
        mmat(6,4)=0.25D0
        mmat(6,5)=0.00D0
        mmat(6,6)=0.50D0

        do i=1,6
            do j=1,6
                mmat(i,j)=mmat(i,j)*gamma(ne)*area*t(ne)/3.0D0
            end do
        end do
        do i=1,3
            w(i*2-1)=gkh(ne)
            w(i*2)=gkv(ne)
        end do
        do i=1,6
            fevec(i)=0.0D0
            do j=1,6
                fevec(i)=fevec(i)+mmat(i,j)*w(j)
            end do
        end do

    end subroutine MV_PL3


    subroutine CALST_PL3TS(ne,kakom,x,y,deltaT,dist,NSTRES,Em,po,t,alpha,ts,strsave,fevec,noten,NELT,NODT,nt)
        !*****************************************************************
        !Calculation of stresses & Internal force vector for 3 node element
        !*****************************************************************
        integer,intent(in)::ne,NELT,NODT,nt,NSTRES
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),y(1:NODT),deltaT(1:NODT),dist(1:nt)
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT),alpha(1:NELT),ts(1:NELT)
        real(8),intent(out)::strsave(1:NELT,1:3),fevec(1:6)
        integer,intent(inout)::noten(1:NELT)

        integer::i,j,ia,k1,k2,k3
        real(8)::b(1:3,1:6),d(1:3,1:3),epsilon(1:3),eps0(1:3),area
        real(8)::sigx,sigy,tauxy,Te,ps1,ps2,ang,phi
        real(8)::wd(1:6),ws(1:6),Eme,poe

        !Displacement increments of element
        do i=1,3
            ia=kakom(ne,i)
            wd(2*i-1)=dist(2*ia-1)
            wd(2*i)=dist(2*ia)
        end do
        !Mean temperature change of element
        k1=kakom(ne,1)
        k2=kakom(ne,2)
        k3=kakom(ne,3)
        Te=(deltaT(k1)+deltaT(k2)+deltaT(k3))/3.0D0

        !Stresses due to nodal displacement
        area=CALA_PL3(ne,kakom,x,y,NELT,NODT)
        call BMAT_PL3(ne,kakom,x,y,b,NELT,NODT)
        Eme=Em(ne)
        poe=po(ne)
        if(noten(ne)/=0)poe=0.0D0
        call DMAT_PL(NSTRES,Eme,poe,d)
        do i=1,3
            epsilon(i)=0.0D0
            do j=1,6
                epsilon(i)=epsilon(i)+b(i,j)*wd(j)
            end do
        end do
        !Thermal strain
        select case(NSTRES)
            case(0)!Plane strain
                eps0(1)=(1.0+poe)*alpha(ne)*Te
            case(1)!Plane stress
                eps0(1)=alpha(ne)*Te
        end select
        eps0(2)=eps0(1)
        eps0(3)=0.0D0
        !Stresses {sigma}=[D]{epsilon-eps0}
        do i=1,3
            ws(i)=0.0D0
            do j=1,3
                ws(i)=ws(i)+d(i,j)*(epsilon(j)-eps0(j))
            end do
        end do

        !Judgement of no-tension state
        sigx=ws(1)
        sigy=ws(2)
        tauxy=ws(3)
        call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
        !Judgement of no-tension state
        if(ts(ne)<ps1)noten(ne)=1
        if(ts(ne)<ps2)noten(ne)=2
        !Elimination of tensile stress component
        phi=-ang/180.0*pi
        if(noten(ne)/=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)==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
        !Memory of stresses
        strsave(ne,1)=sigx
        strsave(ne,2)=sigy
        strsave(ne,3)=tauxy

        !Internal force vector
        do i=1,6
            fevec(i)=0.0D0
            do j=1,3
                fevec(i)=fevec(i)+b(j,i)*strsave(ne,j)*area*t(ne)
            end do
        end do

    end subroutine CALST_PL3TS


    subroutine BMAT_PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
        !*****************************************************************
        !Displacement-strain relationship [B] for 4 node element
        !*****************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(NELT,4)
        real(8),intent(in)::x(NODT),y(NODT),a,b
        real(8),intent(out)::bm(3,8),detJ

        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 matrix and 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/dx][dN/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)
        !*****************************************************************
        !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_PL4(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
        !*****************************************************************
        !Element stiffness matrix
        !*****************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        integer,intent(in)::NSTRES
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT)
        real(8),intent(out)::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
        !Stiffness matrix [sm]=[B]T[D][B]*t*det(J)
        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)
        !*****************************************************************
        !Body force vector for 4 node element
        !*****************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::t(1:NELT),gamma(1:NELT),gkh(1:NELT),gkv(1:NELT)
        real(8),intent(out)::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 matrix and 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
        !Body force
        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)
        !*****************************************************************
        !Calculation of stresses & Internal force vectorfor 4 node element
        !*****************************************************************
        integer,intent(in)::NELT,NODT,nt
        integer,intent(in)::ne,kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT),deltaT(1:NODT),dist(1:nt)
        integer,intent(in)::NSTRES
        real(8),intent(in)::Em(1:NELT),po(1:NELT),t(1:NELT),alpha(1:NELT),ts(1:NELT)
        real(8),intent(out)::fevec(1:8),strsave(1:NELT,1:4,1:3)
        integer,intent(inout)::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

        !Nodal displacement of element
        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)

        !Stresses of element
        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)
            !Strain {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
            !Temperature change of Gauss point
            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)
            !Thermal strain
            select case(NSTRES)
                case(0) !Plane strain
                    eps0(1)=tem*(1.0D0+poe)*alpha(ne)
                    eps0(2)=eps0(1)
                    eps0(3)=0.0D0
                case(1) !Plane stress
                    eps0(1)=tem*alpha(ne)
                    eps0(2)=eps0(1)
                    eps0(3)=0.0D0
            end select
            !Stresses {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
            !Memory of current stresses
            sigx =stress(1)
            sigy =stress(2)
            tauxy=stress(3)
            call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
            !Judgement of no-tension state
            noten(ne,kk)=0
            if(ts(ne)<ps1)noten(ne,kk)=1
            if(ts(ne)<ps2)noten(ne,kk)=2
            !Elimination of tensile stress component
            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
            !Memory of stresses
            strsave(ne,kk,1)=sigx
            strsave(ne,kk,2)=sigy
            strsave(ne,kk,3)=tauxy
        end do

        !Internal force vector
        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)
        !*****************************************************************
        !Calculation of principal stress
        !*****************************************************************
        real(8),intent(in)::sigx,sigy,tauxy
        real(8),intent(out)::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 ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !*****************************************************************
        !Assembly of vector
        !*****************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,nod,nhen
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::fevec(1:nod*nhen)
        real(8),intent(inout)::ftvec(1:NODT*nhen)

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

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


    subroutine ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
        !*****************************************************************
        !Assembly of matrix
        !*****************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,nod,nhen
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::sm(1:nod*nhen,1:nod*nhen)
        real(8),intent(inout)::tsm(1:NODT*nhen,1:NODT*nhen)

        integer :: i,j,k,l,ki,kj,ks,js,kik,kjl,ksk,jsl

        do i=1,nod
            do j=1,nod
                ki=(kakom(ne,i)-1)*nhen
                kj=(kakom(ne,j)-1)*nhen
                ks=(i-1)*nhen
                js=(j-1)*nhen
                do k=1,nhen
                    do l=1,nhen
                        kik=ki+k
                        kjl=kj+l
                        ksk=ks+k
                        jsl=js+l
                        tsm(kik,kjl)=tsm(kik,kjl)+sm(ksk,jsl)
                    end do
                end do
            end do
        end do
    end subroutine ASMAT


    subroutine 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
        !-----------------------------------------------------
        !Calculation of band width
        !-----------------------------------------------------
        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,k,mn
        real(8) :: work(1:nt)
        !-----------------------------------------------------
        !Memody of band matrix
        !-----------------------------------------------------
        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_plnt
