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

program f90_fem_ini4
    use defpi
    implicit none
    integer:: i,j,k,l,kk,kkk,ne,n,ia,ja,mn
    integer:: NSTRES
    integer:: NODT
    integer:: NELT
    integer:: MATEL
    integer:: KOX
    integer:: KOY
    integer:: NF
    integer:: IPR
    integer:: nt
    integer:: nod=4
    integer:: nhen=2
    integer:: mm
    integer:: ib
    real(8):: fact0=1.0E-30
    real(8):: elarea1
    real(8):: elarea2
    real(8):: xg,yg,sigx,sigy,tauxy,ps1,ps2,ang,cc,fa,st,SF
    real(8):: gamma,sy0,r
    real(8):: ps1max=-1.0D010,ps1min=1.0D010
    real(8):: ps2max=-1.0D010,ps2min=1.0D010
    real(8):: sfmax=-1.0D010,sfmin=1.0D010
    real(8):: dismax=0.0D0

    integer,allocatable::kakom(:,:)
    real(8),allocatable:: Em(:)
    real(8),allocatable:: po(:)
    real(8),allocatable:: t(:)
    real(8),allocatable:: chs(:)
    real(8),allocatable:: phi(:)
    real(8),allocatable:: sigt(:)
    integer,allocatable:: matno(:)
    real(8),allocatable:: wm(:,:)
    real(8),allocatable:: strsave(:,:,:)
    real(8),allocatable:: x(:)
    real(8),allocatable:: y(:)
    integer,allocatable:: nokx(:)
    integer,allocatable:: noky(:)
    integer,allocatable:: nonf(:)
    integer,allocatable:: mindex(:)
    real(8),allocatable:: tsm(:,:)
    real(8),allocatable:: f(:)
    real(8),allocatable:: ftvec(:)
    real(8),allocatable:: finit(:)
    real(8),allocatable:: fdis(:)
    real(8),allocatable:: rdis(:)
    real(8),allocatable:: dis(:)
    real(8),allocatable:: reac(:)
    real(8):: sm(1:8,1:8)
    real(8):: fevec(1:8)
    character(len=50)::fnameR,fnameW,strcom
    character(len=20)::dummy

    call getarg(1,dummy); read(dummy,*) gamma ;dummy=''! Unit weight of the rock
    call getarg(2,dummy); read(dummy,*) sy0   ;dummy=''! Initial vertical stress at y=0 on the model
    call getarg(3,dummy); read(dummy,*) r     ;dummy=''! Lateral pressure ratio
    call getarg(4,fnameR) !Input file name
    call getarg(5,fnameW) !Output file name

    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR
    nt=NODT*2
    allocate(kakom(1:NELT,1:4),matno(1:NELT))
    allocate(Em(1:NELT),po(1:NELT),t(1:NELT),chs(1:NELT),phi(1:NELT),sigt(1:NELT))
    allocate(wm(1:MATEL,1:6))
    allocate(strsave(1:NELT,1:4,1:3))
    allocate(x(1:NODT),y(1:NODT))
    allocate(nokx(1:KOX),noky(1:KOY),nonf(1:NF))
    allocate(mindex(1:nt),tsm(1:nt,1:nt),f(1:nt),ftvec(1:nt),finit(1:nt))
    allocate(fdis(1:nt),rdis(1:nt),dis(1:nt),reac(1:nt))
    do i=1,MATEL
        read(11,*) (wm(i,j),j=1,6)
    end do
    do ne=1,NELT
        read(11,*) (kakom(ne,j),j=1,4),matno(ne)
        do i=1,MATEL
            if(i==matno(ne))then
                t(ne)    =wm(i,1)
                Em(ne)   =wm(i,2)
                po(ne)   =wm(i,3)
                chs(ne)  =wm(i,4)
                phi(ne)  =wm(i,5)
                sigt(ne) =wm(i,6)
            end if
        end do
    end do
    do i=1,NODT
        read(11,*) x(i),y(i)
    end do
    do i=1,nt
        rdis(i)=0.0
    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
    if(0<NF)then
        do i=1,NF
            read(11,*) nonf(i)
        end do
    end if
    close(11)

    do ne=1,NELT
        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        elarea1=0.5*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        if(elarea1<fact0)then
            write(6,'(a)') 'area of element is less or equal to 0'
            stop
        end if
        i=kakom(ne,1)
        j=kakom(ne,3)
        k=kakom(ne,4)
        elarea2=0.5*((x(k)-x(j))*y(i)+(x(i)-x(k))*y(j)+(x(j)-x(i))*y(k))
        if(elarea2<fact0)then
            write(6,'(a)') 'area of element is less or equal to 0'
            stop
        end if
    end do

    do i=1,nt
        finit(i)=0.0
        ftvec(i)=0.0
        reac(i)=0.0
        do j=1,nt
            tsm(i,j)=0.0
        end do
    end do
    call MAKE_INISTR(gamma,sy0,r,kakom,x,y,strsave,NELT,NODT,nod)
    do ne=1,NELT
        call PLANESM(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)

        call INITSTR(ne,kakom,x,y,NSTRES,t,strsave,fevec,NELT,NODT)
        call ASVEC(ne,nod,nhen,kakom,fevec,finit,NELT,NODT)
    end do

    do i=1,nt
        fdis(i)=0.0
        do j=1,nt
            fdis(i)=fdis(i)+rdis(j)*tsm(i,j)
        end do
    end do
    do i=1,nt
        f(i)=0.0
    end do
    if(0<NF)then
        do i=1,NF
            n=nonf(i)
            f(2*n-1)=-finit(2*n-1)
            f(2*n)=-finit(2*n)
        end do
    end if
    do i=1,nt
        reac(i)=f(i)-fdis(i)
    end do
    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)
        reac(i)=reac(ia)
        do j=1,mm
            ja=mindex(j)
            tsm(i,j)=tsm(ia,ja)
        end do
    end do
    ib=IBAND(mm,tsm,fact0,nt)
    call BAND(mm,ib,tsm,nt)
    do i=1,mm
        if(abs(tsm(i,1))<fact0)then
            write(6,'(a)') 'diagonal value is less or equal to 0'
            stop
        end if
    end do
    call CHBAND10(mm,ib,tsm,nt)
    call CHBAND11(mm,ib,tsm,reac,nt)

    do i=1,nt
        dis(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
        dis(i)=dis(i)+rdis(i)
    end do
    do i=1,nt
        if(dismax<sqrt(dis(2*i-1)**2+dis(2*i)**2))dismax=sqrt(dis(2*i-1)**2+dis(2*i)**2)
        ftvec(i)=0.0
    end do
    do ne=1,NELT
        call PRST(ne,kakom,x,y,dis,NSTRES,Em,po,t,strsave,fevec,NELT,NODT,nt)
        call ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
    end do

    open(12, file=fnameW, status='replace')
    write(12,'(a)') strcom
    write(12,'(9a8)') 'nod','NODT','NELT','MATEL','KOX','KOY','NF','NSTRES','IPR'
    write(12,'(9i8)') nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR
    write(12,'(a)') '*node characteristics'
    write(12,'(a8,4a15,2a8,2a15)') 'node','x','y','fx','fy','fix-x','fix-y','rdis-x','rdis-y'
    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(12,'(i8,4e15.7,2i8,2e15.7)') i,x(i),y(i),f(2*i-1),f(2*i),k,l,rdis(2*i-1),rdis(2*i)
    end do
    write(12,'(a)') '*element characteristics'
    write(12,'(5a8,6a15,a8)') 'element','node-1','node-2','node-3','node-4','E','po','t','chs','phi','sigt','matno'
    do ne=1,NELT
        write(12,'(5i8,6e15.7,i8)') ne,kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4), &
              Em(ne),po(ne),t(ne),chs(ne),phi(ne),sigt(ne),matno(ne)
    end do
    write(12,'(a)') '*displacements and forces'
    write(12,'(a8,8a15)') 'node','coord-x','coord-y','dis-x','dis-y','reac-x','reac-y','ftvec-x','ftvec-y'
    do i=1,NODT
        write(12,'(i8,8e15.7)') i,x(i),y(i),dis(2*i-1),dis(2*i),reac(2*i-1),reac(2*i),ftvec(2*i-1),ftvec(2*i)
    end do
    write(12,'(a)') '*stresses'
    write(12,'(2a8,9a15,a8)') 'element','kk','coord-x','coord-y','sig-x','sig-y','tau-xy','ps1','ps2','ang','SF','matno'
    do ne=1,NELT
        cc=chs(ne)
        fa=phi(ne)
        st=sigt(ne)
        select case(IPR)
        case(0)
            do kk=1,4
                sigx=strsave(ne,kk,1)
                sigy=strsave(ne,kk,2)
                tauxy=strsave(ne,kk,3)
                call PSTRESS(sigx,sigy,tauxy,ps1,ps2,ang)
                call XGYG4I(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
                SF=CALSF(ps1,ps2,cc,fa,st)
                write(12,'(2i8,9e15.7,i8)') ne,kk,xg,yg,sigx,sigy,tauxy,ps1,ps2,ang,SF,matno(ne)
                if(ps1max<ps1)ps1max=ps1
                if(ps1<ps1min)ps1min=ps1
                if(ps2max<ps2)ps2max=ps2
                if(ps2<ps2min)ps2min=ps2
                if(sfmax<SF)sfmax=SF
                if(SF<sfmin)sfmin=SF
            end do
        case(1)
            kk=1
            sigx =0.25D0*(strsave(ne,1,1)+strsave(ne,2,1)+strsave(ne,3,1)+strsave(ne,4,1))
            sigy =0.25D0*(strsave(ne,1,2)+strsave(ne,2,2)+strsave(ne,3,2)+strsave(ne,4,2))
            tauxy=0.25D0*(strsave(ne,1,3)+strsave(ne,2,3)+strsave(ne,3,3)+strsave(ne,4,3))
            call PSTRESS(sigx,sigy,tauxy,ps1,ps2,ang)
            call XGYG4M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
            SF=CALSF(ps1,ps2,cc,fa,st)
            write(12,'(2i8,9e15.7,i8)') ne,kk,xg,yg,sigx,sigy,tauxy,ps1,ps2,ang,SF,matno(ne)
                if(ps1max<ps1)ps1max=ps1
                if(ps1<ps1min)ps1min=ps1
                if(ps2max<ps2)ps2max=ps2
                if(ps2<ps2min)ps2min=ps2
                if(sfmax<SF)sfmax=SF
                if(SF<sfmin)sfmin=SF
        end select
    end do
    write(12,'("#,Summary")')
    write(12,'("#,NELT=",i0," NODT=",i0," nt=",i0," mm=",i0," ib=",i0)'),NELT,NODT,nt,mm,ib
    write(12,'("#,ps1_range,",e15.7,",",e15.7)') ps1min,ps1max
    write(12,'("#,ps2_range,",e15.7,",",e15.7)') ps2min,ps2max
    write(12,'("#,sf_range,",e15.7,",",e15.7)') sfmin,sfmax
    write(12,'("#,dis_max,",e15.7)') dismax
    close(12)

stop

contains

    subroutine MAKE_INISTR(gamma,sy0,r,kakom,x,y,strsave,NELT,NODT,nod)
        real(8),intent(in)::gamma,sy0,r
        integer,intent(in)::NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(out)::strsave(1:NELT,1:4,1:3)

        integer::ne,kk
        real(8)::xg,yg
        real(8)::sx,sy,tu

        do ne=1,NELT
            do kk=1,4
                call XGYG4I(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod)
                sy=-(sy0-gamma*yg)
                sx=r*sy
                tu=0.0D0
                strsave(ne,kk,1)=sx
                strsave(ne,kk,2)=sy
                strsave(ne,kk,3)=tu
            end do
        end do
    end subroutine MAKE_INISTR


    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 IPOINT(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(ne,NSTRES,Em,po,d,NELT)
        integer,intent(in)::ne,NSTRES,NELT
        real(8),intent(in)::Em(1:NELT),po(1:NELT)
        real(8),intent(out)::d(1:3,1:3)
        integer:: i,j

        do i=1,3
            do j=1,3
                d(i,j)=0.0
            end do
        end do
        if(NSTRES==0)then
            d(1,1)=Em(ne)*(1.0-po(ne))/((1.0+po(ne))*(1.0-2.0*po(ne)))
            d(1,2)=Em(ne)*po(ne)/((1.0+po(ne))*(1.0-2.0*po(ne)))
            d(2,1)=d(1,2)
            d(2,2)=d(1,1)
            d(3,3)=Em(ne)*0.5/(1.0+po(ne))
        else
            d(1,1)=Em(ne)/(1.0-po(ne)*po(ne))
            d(1,2)=po(ne)*d(1,1)
            d(2,1)=d(1,2)
            d(2,2)=d(1,1)
            d(3,3)=Em(ne)*0.5/(1.0+po(ne))
        end if
    end subroutine DMAT


    subroutine BMAT(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::a,b
        real(8),intent(out)::bm(1:3,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

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

        dn1a=-0.25*(1.0-b)
        dn2a= 0.25*(1.0-b)
        dn3a= 0.25*(1.0+b)
        dn4a=-0.25*(1.0+b)
        dn1b=-0.25*(1.0-a)
        dn2b=-0.25*(1.0+a)
        dn3b= 0.25*(1.0+a)
        dn4b= 0.25*(1.0-a)

        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 i=1,3
            do j=1,8
                bm(i,j)=0.0
            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


    subroutine IPOINT(kk,a,b)
        integer,intent(in)::kk
        real(8),intent(out)::a,b

        a=0.0
        b=0.0
        select case(kk)
            case(1)
                a=-0.5773502692
                b=-0.5773502692
            case(2)
                a= 0.5773502692
                b=-0.5773502692
            case(3)
                a= 0.5773502692
                b= 0.5773502692
            case(4)
                a=-0.5773502692
                b= 0.5773502692
        end select
    end subroutine IPOINT


    subroutine PLANESM(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
        integer,intent(in)::ne,NSTRES,NODT,NELT
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::Em(1:NELT),po(1:NELT)
        real(8),intent(inout)::t(1:NELT)
        real(8),intent(out)::sm(1:8,1:8)

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

        do i=1,8
            do j=1,8
                sm(i,j)=0.0
            end do
        end do

        if(NSTRES==0)t(ne)=1.0
        call DMAT(ne,NSTRES,Em,po,d,NELT)

        do kk=1,4
            call IPOINT(kk,a,b)
            call BMAT(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 PLANESM


    subroutine INITSTR(ne,kakom,x,y,NSTRES,t,strsave,fevec,NELT,NODT)
        integer,intent(in)::ne,NSTRES,NODT,NELT
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::strsave(1:NELT,4,3)
        real(8),intent(inout)::t(1:NELT)
        real(8),intent(out)::fevec(1:8)

        integer:: kk,i,j
        real(8):: detJ,a,b
        real(8):: bm(1:3,1:8)

        if(NSTRES==0)t(ne)=1.0
        do i=1,8
            fevec(i)=0.0
        end do
        do kk=1,4
            call IPOINT(kk,a,b)
            call BMAT(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 INITSTR


    subroutine PRST(ne,kakom,x,y,dis,NSTRES,Em,po,t,strsave,fevec,NELT,NODT,nt)
        integer,intent(in)::ne,NSTRES,NODT,NELT,nt
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),y(1:NODT),dis(1:nt)
        real(8),intent(in)::Em(1:NELT),po(1:NELT)
        real(8),intent(inout)::t(1:NELT)
        real(8),intent(out)::strsave(1:NELT,1:4,1:3),fevec(1:8)

        integer:: kk,i,j,ia
        real(8):: detJ,a,b
        real(8):: d(1:3,1:3)
        real(8):: bm(1:3,1:8)
        real(8):: epsilon(1:3)
        real(8):: stress(1:3)
        real(8):: wd(1:8)

        do i=1,4
            ia=kakom(ne,i)
            wd(2*i-1)=dis(2*ia-1)
            wd(2*i)=dis(2*ia)
        end do
        if(NSTRES==0)t(ne)=1.0
        call DMAT(ne,NSTRES,Em,po,d,NELT)
        do kk=1,4
            call IPOINT(kk,a,b)
            call BMAT(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
            do i=1,3
                epsilon(i)=0.0
                do j=1,8
                    epsilon(i)=epsilon(i)+bm(i,j)*wd(j)
                end do
            end do
            do i=1,3
                stress(i)=0.0
                do j=1,3
                    stress(i)=stress(i)+d(i,j)*epsilon(j)
                end do
            end do
            strsave(ne,kk,1)=stress(1)+strsave(ne,kk,1)
            strsave(ne,kk,2)=stress(2)+strsave(ne,kk,2)
            strsave(ne,kk,3)=stress(3)+strsave(ne,kk,3)
        end do

        do i=1,8
            fevec(i)=0.0
        end do
        do kk=1,4
            call IPOINT(kk,a,b)
            call BMAT(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 PRST


    subroutine PSTRESS(sigx,sigy,tauxy,ps1,ps2,ang)
        real(8),intent(in)::sigx,sigy,tauxy
        real(8),intent(out)::ps1,ps2,ang

        ps1=0.5D0*(sigx+sigy)+sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        ps2=0.5D0*(sigx+sigy)-sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        ang=0.0D0
        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=180.0D0/pi*(0.5D0*atan(2.0*tauxy/(sigx-sigy)))
            if(sigx>sigy.and.tauxy>=0.0D0) ang=ang+0.0D0
            if(sigx>sigy.and.tauxy<0.0D0)  ang=ang+180.0D0
            if(sigx<sigy)                  ang=ang+90.0D0
        end if
    end subroutine PSTRESS


    real(8) function CALSF(ps1,ps2,cc,fa,st)
        real(8),intent(in)::ps1,ps2,cc,fa,st

        real(8):: D1,D2,dd1,dd2,SF1,SF2,SF
        real(8):: factL=1.0e10

        D1=(cc-0.5D0*(ps1+ps2)*tan(fa/180.0D0*pi))*cos(fa/180.0D0*pi)
        dd1=D1-0.5D0*(ps1-ps2)
        D2=st-0.5D0*(ps1+ps2)
        dd2=st-ps1
        if(ps1==ps2)then
            if(0.0D0<cc-ps1*tan(fa/180.0D0*pi))then
                SF1=factL
            else
                SF1=0.0D0
            end if
            if(ps1<=0.0D0)then
                SF2=factL
            else
                SF2=st/ps1
            end if
        else
            SF1=D1/(D1-dd1)
            SF2=D2/(D2-dd2)
        end if
        CALSF=min(SF1,SF2)
    end function CALSF


    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

end program f90_fem_ini4
