program f90_fem_seep
    !******************************************************************************************
    !2D Unsteady Seepage Flow Analysis with the function of saturation-unsaturation analysis
    !******************************************************************************************
    implicit none
    integer::i,j,k,l,k1,k2,k3,n,ia,ja,ne,nnn
    character::strcom*256           !Comment
    integer::NODT                   !Number of nodes
    integer::NELT                   !Number of elements
    integer::MATEL                  !Number of material sets
    integer::KOH                    !Number of nodes with given total head
    integer::KOQ                    !Number of nodes with given discharge
    integer::KOU                    !Number of nodes with the condition of seepage point
    integer::KOP                    !Number of nodes with zero pressure head in seepage points
    integer::nt                     !Number of degrees of freedom of all nodes (NODT x nhen)
    integer::mm                     !Number of degrees of freedom of reduced FE equation
    integer::ib                     !Band width of matrix
    integer::nod                    !Number of nodes per element
    integer::nhen=1                 !Number of degree of freedom per node
    integer::idan                   !Section (0: vertical, 1: horizontal)
    integer,allocatable::kakom(:,:) !Element connectivity
    real(8),allocatable::Ak0(:)     !Saturated permeability coefficient of element
    real(8),allocatable::Akr(:)     !Relative hydraulic conductivity (ratio of permeability: unsat/sat)
    real(8),allocatable::alpha(:)   !Unsaturated property parameter 'alpha' of element
    real(8),allocatable::em(:)      !Unsaturated property parameter 'm' of element
    integer,allocatable::matno(:)   !Material set number
    real(8),allocatable::wm(:,:)    !work for material properties
    real(8),allocatable::x(:)       !x-coordinate of node
    real(8),allocatable::z(:)       !z-coordinate of node
    integer,allocatable::nokh(:)    !Node number with given total head
    integer,allocatable::nokq(:)    !Node number with given discharge
    integer,allocatable::noku(:)    !Node number with the condition of seepage point
    integer,allocatable::nokp(:)    !Node number with zero pressure head in seepage points
    real(8),allocatable::Hinp(:)    !Given total head
    real(8),allocatable::Qinp(:)    !Given discharge
    real(8),allocatable::tak(:,:)   !Permeability matrix
    real(8),allocatable::wak(:,:)   !work for nodes with giventotal head
    real(8),allocatable::hvec(:)    !Total head vector
    real(8),allocatable::h1vec(:)   !Total head vector
    real(8),allocatable::h2vec(:)   !Total head vector
    real(8),allocatable::qvec(:)    !Discharge vector
    real(8),allocatable::q0vec(:)   !Discharge vector
    real(8),allocatable::reac(:)    !work vector
    real(8),allocatable::vx(:)      !Element velocity vector in x-direction
    real(8),allocatable::vz(:)      !Element velocity vector in z-direction
    real(8),allocatable::vm(:)      !Element mean velocity vector
    real(8),allocatable::mindex(:)  !Index for treatment of nodes with given total head
    real(8)::eak(1:4, 1:4)          !Element permeability matrix
    real(8)::fact0=1.0E-30          !0 judgement
    real(8)::elarea                 !Half area of element
    integer::itmax=100              !Limit value of iterative calculation
    integer::icount                 !Number of converged degrees of freedom
    real(8)::QTi,QTo                !Discharge of in-flow, discharge of out-flow
    real(8)::vmmax,xgm,zgm          !Max.velocity, coordinates of centroid of element with max.velocity
    real(8)::vimax,xgi,zgi          !Max.velocity in in-flow area, coordinates of centroid of element with max.velocity
    real(8)::vomax,xgo,zgo          !Max.velocity in out-flow area, coordinates of centroid of element with max.velocity
    real(8)::xg,zg                  !Coordinates of centroid of element
    integer::kmmax,kimax,komax      !Element number with max.velocity
    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),3(',',i2),2(',',e15.7))"             !Node
    fmt2_3="(i5,3(',',i5),3(',',e15.7),',',i5)"                 !3 node element
    fmt2_4="(i5,4(',',i5),3(',',e15.7),',',i5)"                 !4 node element
    fmt3="(i5,2(',',e15.7),3(',',e15.7))"                       !head and discharge
    fmt4="(i5,2(',',e15.7),3(',',e15.7),2(',',e15.7),',',i5)"   !Velocity and permeability

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

    !To start the measurement of calculation taime
    t1=DIFT()
    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    !!--------------------------------------------------------------------------------------
    !!Comment
    !!--------------------------------------------------------------------------------------
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
    !Basic value for analysis
    !------------------------------------------------------------------------------------------
    read(11,*) nod,NODT,NELT,MATEL,KOH,KOQ,KOU,idan
    !------------------------------------------------------------------------------------------
    !To declare array size
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:3))
    allocate(Ak0(1:NELT))
    allocate(Akr(1:NELT))
    allocate(alpha(1:NELT))
    allocate(em(1:NELT))
    allocate(x(1:NODT))
    allocate(z(1:NODT))
    if(0<KOH)allocate(nokh(1:KOH))
    if(0<KOQ)allocate(nokq(1:KOQ))
    if(0<KOU)allocate(noku(1:KOU))
    if(0<KOU)allocate(nokp(1:KOU))
    allocate(Hinp(1:KOH))
    allocate(Qinp(1:KOQ))
    allocate(vx(1:NELT))
    allocate(vz(1:NELT))
    allocate(vm(1:NELT))
    allocate(mindex(1:nt))
    allocate(tak(1:nt,1:nt))
    allocate(h1vec(1:nt))
    allocate(h2vec(1:nt))
    allocate(hvec(1:nt))
    allocate(qvec(1:nt))
    allocate(q0vec(1:nt))
    allocate(reac(1:nt))
    if(0<KOH.and.KOU<=KOH)allocate(wak(1:nt,1:KOH))
    if(0<KOU.and.KOH<KOU)allocate(wak(1:nt,1:KOU))
    !------------------------------------------------------------------------------------------
    !Material properties (k0,alpha,m)
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3)
    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
                Ak0(ne)  =wm(i,1)
                alpha(ne)=wm(i,2)
                em(ne)   =wm(i,3)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !(x,z) coordinates, initial value of total head
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),z(i),hvec(i)
    end do
    !------------------------------------------------------------------------------------------
    !Node number with given total head, value of given total head
    !Node number with given discharge, value of given discharge
    !Node number with the condition of seepage point
    !------------------------------------------------------------------------------------------
    if(0<KOH)then
        do i=1,KOH
            read(11,*) nokh(i),Hinp(i)
        end do
    end if
    if(0<KOQ)then
        do i=1,KOQ
            read(11,*) nokq(i),Qinp(i)
        end do
    end if
    if(0<KOU)then
        do i=1,KOU
            read(11,*) noku(i)
        end do
    end if
    close(11)

    !*****************************************************
    !To print input data
    !*****************************************************
    do i=1,nt
        qvec(i)=0.0D0
    end do
    if(0<KOH)then
        do k=1,KOH
            hvec(nokh(k))=Hinp(k) !Amendment of initial total head to given value
        end do
    end if
    if(0<KOQ)then
        do k=1,KOQ
            qvec(nokq(k))=Qinp(k) !To set given discharge
        end do
    end if

    open(12, file=fnameW, status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(12,'(a)') 'nod,NODT,NELT,MATEL,KOH,KOQ,KOU,idan'
        write(linebuf,*) nod,',',NODT,',',NELT,',',MATEL,',',KOH,',',KOQ,',',KOU,',',idan
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,'(a)') '*node characteristics'
    write(12,'(a)') 'node,x,z,KOH,KOQ,KOU,Hinp,Qinp'
    do i=1,NODT
        k1=0
        k2=0
        k3=0
        if(0<KOH)then
            do j=1,KOH
                if(i==nokh(j))k1=1
            end do
        end if
        if(0<KOQ)then
            do j=1,KOQ
                if(i==nokq(j))k2=1
            end do
        end if
        if(0<KOU)then
            do j=1,KOU
                if(i==noku(j))k3=1
            end do
        end if
            write(linebuf,FMT=fmt1) i,x(i),z(i),k1,k2,k3,hvec(i),qvec(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,Ak0,alpha,em,matno'
            do ne=1,NELT
                write(linebuf,FMT=fmt2_3) ne,kakom(ne,1),kakom(ne,2),kakom(ne,3), &
                    Ak0(ne),alpha(ne),em(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,Ak0,alpha,em,matno'
            do ne=1,NELT
                write(linebuf,FMT=fmt2_4) ne,kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4),&
                    Ak0(ne),alpha(ne),em(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))*z(i)+(x(i)-x(k))*z(j)+(x(j)-x(i))*z(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))*z(i)+(x(i)-x(k))*z(j)+(x(j)-x(i))*z(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
    !*****************************************************
    do i=1,nt
        h1vec(i)=hvec(i)
    end do


    do nnn=1,itmax
        !-----------------------------------------------------
        !Initialization of matrix and vector
        !-----------------------------------------------------
        do i=1,nt
            reac(i)=0.0D0
            do j=1,nt
                tak(i,j)=0.0D0
            end do
        end do
        !-----------------------------------------------------
        !Assembly of permeability matrix, [k] is calculated using {h1}
        !-----------------------------------------------------
        select case(nod)
            case(3)
                do ne=1,NELT
                    call STIFF_S3(ne,kakom,x,z,Ak0,alpha,em,Akr,h1vec,eak,NELT,NODT,nt,idan)
                    call ASMAT(ne,nod,nhen,kakom,eak,tak,NELT,NODT)
                end do
            case(4)
                do ne=1,NELT
                    call STIFF_S4(ne,kakom,x,z,Ak0,alpha,em,Akr,h1vec,eak,NELT,NODT,nt,idan)
                    call ASMAT(ne,nod,nhen,kakom,eak,tak,NELT,NODT)
                end do
        end select
        !-----------------------------------------------------
        !calculation of nodal discharge {q}=[k]{h1}
        !-----------------------------------------------------
        do i=1,nt
            qvec(i)=0.0D0
            do j=1,nt
                qvec(i)=qvec(i)+tak(i,j)*h1vec(j)
            end do
        end do
        !-----------------------------------------------------
        !Uodate of nodal discharge {q}->{q0}
        !(In-flow is not allowed at seepage point. In-flow means positive discharge.)
        !(If discharge has positive value at seepage point, the discharge is set to zero)
        !-----------------------------------------------------
        do i=1,nt
            q0vec(i)=0.0D0
        end do
        if(0<KOH)then
            do k=1,KOH
                q0vec(nokh(k))=qvec(nokh(k))
            end do
        end if
        if(0<KOQ)then
            do k=1,KOQ
                q0vec(nokq(k))=Qinp(k)
            end do
        end if
        if(0<KOU)then
            do k=1,KOU
                q0vec(noku(k))=qvec(noku(k))
                if(0.0D0<q0vec(noku(k)))q0vec(noku(k))=0.0D0
            end do
        end if
        !-----------------------------------------------------
        !Memory of node numbers with positive pressure head in seepage points
        !-----------------------------------------------------
        KOP=0
        if(0<KOU)then
            do i=1,KOU
                k=noku(i)
                if(z(k)<h1vec(k))then
                    KOP=KOP+1
                    nokp(KOP)=k
                end if
            end do
        end if
        !-----------------------------------------------------
        !Amendment of nodal values {q}={q0}-[k]{h(known)}
        !-----------------------------------------------------
        do i=1,nt
            hvec(i)=0.0D0
            qvec(i)=q0vec(i)
            reac(i)=0.0D0
        end do
        !To set total head at the node with given total head
        if(0<KOH)then
            do k=1,KOH
                do i=1,nt
                    wak(i,k)=tak(i,nokh(k))
                end do
            end do
            do k=1,KOH
                do i=1,nt
                    hvec(i)=hvec(i)+Hinp(k)*wak(i,k)
                end do
            end do
        end if
        !To set zero pressure head at the node with positive presure head in seepage points
        !(Total head = position head)
        if(0<KOP)then
            do k=1,KOP
                do i=1,nt
                    wak(i,k)=tak(i,nokp(k))
                end do
            end do
            do k=1,KOP
                do i=1,nt
                    reac(i)=reac(i)+z(nokp(k))*wak(i,k)
                end do
            end do
        end if
        do i=1,nt
            qvec(i)=qvec(i)-hvec(i)-reac(i)
        end do
        !-----------------------------------------------------
        !Treatment of nodes with given total head
        !-----------------------------------------------------
        do i=1,nt
            mindex(i)=i
            reac(i)=0.0D0
        end do
        if(0<KOH)then
            do i=1,KOH
                n=nokh(i)
                mindex(n)=0
            end do
        end if
        if(0<KOP)then
            do i=1,KOP
                n=nokp(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)=qvec(ia)
            do j=1,mm
                ja=mindex(j)
                tak(i,j)=tak(ia,ja)
            end do
        end do
        !-----------------------------------------------------
        !Calculation of band width & memory of band matrix
        !-----------------------------------------------------
        ib=IBAND(mm,tak,fact0,nt)
        call BAND(mm,ib,tak,nt)
        !-----------------------------------------------------
        !To verify diagonal elements & to stop if zero value
        !-----------------------------------------------------
        do i=1,mm
            if(abs(tak(i,1))<fact0)then
                open(12, file='fnameW.csv', access='append')
                write(12,'(a)') '***program stop due to diagonal member=0'
                close(12)
                stop
            end if
        end do
        !-----------------------------------------------------
        !Cholesky method
        !-----------------------------------------------------
        call CHBAND10(mm,ib,tak,nt)!Triangle matrix
        call CHBAND11(mm,ib,tak,reac,nt)!Forward elimination & backward substitution
        !-----------------------------------------------------
        !Calculation of total head {h}=[k(-1)]{q}
        !-----------------------------------------------------
        do i=1,nt
            hvec(i)=0.0D0
        end do
        do i=1,mm
            ia=mindex(i)
            hvec(ia)=reac(i)
        end do
        do k=1,KOH
            hvec(nokh(k))=Hinp(k)
        end do
        do k=1,KOP
            hvec(nokp(k))=z(nokp(k))
        end do
        do i=1,nt
            h2vec(i)=hvec(i)
        end do
        !-----------------------------------------------------
        !Judgement of convergence |h2-h1|<1e-6
        !-----------------------------------------------------
        icount=0
        do i=1,nt
            if(abs(h2vec(i)-h1vec(i))<1.0D-6)icount=icount+1
        end do
            print "('nnn=',i0,' icount=',i5,' KOP=',i5)",nnn,icount,KOP
        if(icount==nt)exit
        !-----------------------------------------------------
        ! Update of total head {h1}={h2}
        !-----------------------------------------------------
        do i=1,nt
            h1vec(i)=h2vec(i)
        end do

    end do


    !-----------------------------------------------------
    !Calculation of total in-flow and out-flow
    !-----------------------------------------------------
    QTi=0.0D0  !Total in-flow
    QTo=0.0D0  !Total out-flow
    do i=1,nt
        if(0<=q0vec(i))QTi=QTi+q0vec(i)
        if(q0vec(i)<0)QTo=QTo+q0vec(i)
    end do
    !-----------------------------------------------------
    !Calculation of element velocity & to search max.velocity
    !-----------------------------------------------------
    vmmax=0.0D0  !Max.velocity
    vimax=0.0D0  !Max.velocity on in-flow area
    vomax=0.0D0  !max.velocity on out-flow area
    kmmax=0      !Element number with max.velocity
    kimax=0      !Element number with max.velocity on in-flow area
    komax=0      !Element number with max.velocity on out-flow area
    select case(nod)
        case(3)
            do ne=1,NELT
                call CALV_S3(ne,kakom,x,z,hvec,Ak0,alpha,em,Akr,vx,vz,NELT,NODT,nt,idan)
                vm(ne)=sqrt(vx(ne)**2.0D0+vz(ne)**2.0D0)
                if(vmmax<vm(ne))then
                    kmmax=ne
                    vmmax=vm(ne)
                end if
                i=kakom(ne,1)
                j=kakom(ne,2)
                k=kakom(ne,3)
                if(0.0D0<q0vec(i).or.0.0D0<q0vec(j).or.0.0D0<q0vec(k))then
                    if(vimax<vm(ne))then
                        kimax=ne
                        vimax=vm(ne)
                    end if
                end if
                if(0.0D0>q0vec(i).or.0.0D0>q0vec(j).or.0.0D0>q0vec(k))then
                    if(vomax<vm(ne))then
                        komax=ne
                        vomax=vm(ne)
                    end if
                end if
            end do
            !Coordinates of centroid of element with max.velocity
            call XGZG3M(kmmax,0,kakom,x,z,xgm,zgm,NELT,NODT,nod)
            call XGZG3M(kimax,0,kakom,x,z,xgi,zgi,NELT,NODT,nod)
            call XGZG3M(komax,0,kakom,x,z,xgo,zgo,NELT,NODT,nod)
        case(4)
            do ne=1,NELT
                call CALV_S4(ne,kakom,x,z,hvec,Ak0,alpha,em,Akr,vx,vz,NELT,NODT,nt,idan)
                vm(ne)=sqrt(vx(ne)**2.0D0+vz(ne)**2.0D0)
                if(vmmax<vm(ne))then
                    kmmax=ne
                    vmmax=vm(ne)
                end if
                i=kakom(ne,1)
                j=kakom(ne,2)
                k=kakom(ne,3)
                l=kakom(ne,4)
                if(0.0D0<q0vec(i).or.0.0D0<q0vec(j).or.0.0D0<q0vec(k).or.0.0D0<q0vec(l))then
                    if(vimax<vm(ne))then
                        kimax=ne
                        vimax=vm(ne)
                    end if
                end if
                if(0.0D0>q0vec(i).or.0.0D0>q0vec(j).or.0.0D0>q0vec(k).or.0.0D0>q0vec(l))then
                    if(vomax<vm(ne))then
                        komax=ne
                        vomax=vm(ne)
                    end if
                end if
            end do
            !Coordinates of centroid of element with max.velocity
            call XGZG4M(kmmax,0,kakom,x,z,xgm,zgm,NELT,NODT,nod)
            call XGZG4M(kimax,0,kakom,x,z,xgi,zgi,NELT,NODT,nod)
            call XGZG4M(komax,0,kakom,x,z,xgo,zgo,NELT,NODT,nod)
    end select

    !*****************************************************
    !To print the result
    !*****************************************************
    open(12, file=fnameW, access='append')
    !Node number, x-coordinate, y-coordinate, total head, discharge, pressure head
    !(positive: in-flow, negative:out-flow in discharge)
    write(12,'(a)') '*node value'
    write(12,'(a)') 'node,coord-x,coord-z,hvec,qvec,pvec'
    do i=1,NODT
        if(idan==0)write(linebuf,FMT=fmt3) i,x(i),z(i),hvec(i),q0vec(i),hvec(i)-z(i)
        if(idan==1)write(linebuf,FMT=fmt3) i,x(i),z(i),hvec(i),q0vec(i),hvec(i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !Mean velocity of element
    write(12,'(a)') '*element value'
    write(12,'(a)') 'element,xg,zg,vx,vz,vm,Ak0,Akr,matno'
    do ne=1,NELT
        if(nod==3)call XGZG3M(ne,0,kakom,x,z,xg,zg,NELT,NODT,nod)
        if(nod==4)call XGZG4M(ne,0,kakom,x,z,xg,zg,NELT,NODT,nod)
        write(linebuf,FMT=fmt4) ne,xg,zg,vx(ne),vz(ne),vm(ne),Ak0(ne),Akr(ne),matno(ne)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !To finish themeasurement of calculation time
    t2=DIFT()
    !Summary出力
    write(12,'("#,Summary")')
    write(12,'("#,item,QTi,QTo")')
    write(12,'("#,discharge,",e15.7,",",e15.7)'),QTi,QTo
    write(12,'("#,item,ne,vmax,xg,zg")')
    write(12,'("#,vmmax,",i5,",",e15.7,",",e15.7,",",e15.7)'),kmmax,vmmax,xgm,zgm
    write(12,'("#,vimax,",i5,",",e15.7,",",e15.7,",",e15.7)'),kimax,vimax,xgi,zgi
    write(12,'("#,vomax,",i5,",",e15.7,",",e15.7,",",e15.7)'),komax,vomax,xgo,zgo
    write(12,'("#,NELT=",i0," NODT=",i0," nt=",i0," mm=",i0," ib=",i0)'),NELT,NODT,nt,mm,ib
    write(12,'("#,Convergence nnn=",i0," icount=",i0)'),nnn,icount
    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 XGZG3M(ne,kk,kakom,x,z,xg,zg,NELT,NODT,nod)
        !*************************************************************
        !Calculation 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),z(1:NODT)
        real(8),intent(out)::xg,zg

        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
        zg=(z(i)+z(j)+z(k))/3.0D0
    end subroutine XGZG3M


    subroutine XGZG4M(ne,kk,kakom,x,z,xg,zg,NELT,NODT,nod)
        !*************************************************************
        !Calculation 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),z(1:NODT)
        real(8),intent(out)::xg,zg

        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))
        zg=0.25D0*(z(i)+z(j)+z(k)+z(l))
    end subroutine XGZG4M


    real(8) function CAL_AKR(alpha_el,em_el,sucp_el)
        !*************************************************************
        !Calculation of relative hydraulic conductivity
        !*************************************************************
        real(8),intent(in)::alpha_el,em_el,sucp_el
        real(8)::Se,en

        CAL_AKR=1.0D0
        if(0.0<sucp_el)then
            en=1.0D0/(1.0D0-em_el)
            Se=(1.0D0/(1.0D0+(alpha_el*sucp_el)**en))**em_el
            CAL_AKR=sqrt(Se)*(1.0D0-(1.0D0-Se**(1.0D0/em_el))**em_el)**2.0D0
        end if

    end function CAL_AKR


    subroutine STIFF_S3(ne,kakom,x,z,Ak0,alpha,em,Akr,h1vec,eak,NELT,NODT,nt,idan)
        !*************************************************************
        !Element permeability matrix for 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,nt,idan
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),z(1:NODT)
        real(8),intent(in)::Ak0(1:NELT),alpha(1:NELT),em(1:NELT)
        real(8),intent(in)::h1vec(1:nt)
        real(8),intent(out)::Akr(1:NELT),eak(1:3,1:3)

        integer::i,j,k
        real(8)::sucp,Akre1,Akre2,Akre3,Akk_el
        real(8)::b1,b2,b3,c1,c2,c3,Delta
        real(8)::wx(1:3,1:3),wz(1:3,1:3)

        do i=1,3
            do j=1,3
                eak(i,j)= 0.0D0
                wx(i,j) = 0.0D0
                wz(i,j) = 0.0D0
            end do
        end do

        i=kakom(ne,1)
        j=kakom(ne,2)
        k=kakom(ne,3)
        select case(idan)
            case(0)
                sucp=-(h1vec(i)-z(i))
                Akre1=CAL_AKR(alpha(ne),em(ne),sucp)
                sucp=-(h1vec(j)-z(j))
                Akre2=CAL_AKR(alpha(ne),em(ne),sucp)
                sucp=-(h1vec(k)-z(k))
                Akre3=CAL_AKR(alpha(ne),em(ne),sucp)
                Akr(ne)=(Akre1+Akre2+Akre3)/3.0D0
            case(1)
                Akr(ne)=1.0D0
        end select
        Akk_el=Ak0(ne)*Akr(ne)
        Delta=0.5D0*((x(k)-x(j))*z(i)+(x(i)-x(k))*z(j)+(x(j)-x(i))*z(k))
        b1=z(j)-z(k)
        b2=z(k)-z(i)
        b3=z(i)-z(j)
        c1=x(k)-x(j)
        c2=x(i)-x(k)
        c3=x(j)-x(i)
        wx(1,1)=b1*b1
        wx(1,2)=b1*b2
        wx(1,3)=b1*b3
        wx(2,1)=b2*b1
        wx(2,2)=b2*b2
        wx(2,3)=b2*b3
        wx(3,1)=b3*b1
        wx(3,2)=b3*b2
        wx(3,3)=b3*b3
        wz(1,1)=c1*c1
        wz(1,2)=c1*c2
        wz(1,3)=c1*c3
        wz(2,1)=c2*c1
        wz(2,2)=c2*c2
        wz(2,3)=c2*c3
        wz(3,1)=c3*c1
        wz(3,2)=c3*c2
        wz(3,3)=c3*c3
        do i=1,3
            do j=1,3
                eak(i,j)=0.25*Akk_el/Delta*(wx(i,j)+wz(i,j))
            end do
        end do

    end subroutine STIFF_S3


    subroutine CALV_S3(ne,kakom,x,z,hvec,Ak0,alpha,em,Akr,vx,vz,NELT,NODT,nt,idan)
        !*************************************************************
        !Calculation of velocity of 3 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,nt,idan
        integer,intent(in)::kakom(1:NELT,1:3)
        real(8),intent(in)::x(1:NODT),z(1:NODT),hvec(1:nt)
        real(8),intent(in)::Ak0(1:NELT),alpha(1:NELT),em(1:NELT)
        real(8),intent(out)::Akr(1:NELT),vx(1:NELT),vz(1:NELT)

        integer::i,j,k
        real(8)::b1,b2,b3,c1,c2,c3,Delta,h1,h2,h3
        real(8)::sucp,Akre1,Akre2,Akre3

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

        select case(idan)
            case(0)
                sucp=-(hvec(i)-z(i))
                Akre1=CAL_AKR(alpha(ne),em(ne),sucp)
                sucp=-(hvec(j)-z(j))
                Akre2=CAL_AKR(alpha(ne),em(ne),sucp)
                sucp=-(hvec(k)-z(k))
                Akre3=CAL_AKR(alpha(ne),em(ne),sucp)
                Akr(ne)=(Akre1+Akre2+Akre3)/3.0D0
            case(1)
                Akr(ne)=1.0D0
        end select

        Delta=0.5D0*((x(k)-x(j))*z(i)+(x(i)-x(k))*z(j)+(x(j)-x(i))*z(k))
        b1=z(j)-z(k)
        b2=z(k)-z(i)
        b3=z(i)-z(j)
        c1=x(k)-x(j)
        c2=x(i)-x(k)
        c3=x(j)-x(i)
        h1=hvec(i)
        h2=hvec(j)
        h3=hvec(k)
        vx(ne)=-0.5D0*Ak0(ne)*Akr(ne)/Delta*(b1*h1+b2*h2+b3*h3)
        vz(ne)=-0.5D0*Ak0(ne)*Akr(ne)/Delta*(c1*h1+c2*h2+c3*h3)

    end subroutine CALV_S3


    subroutine DNMAT_S4(ne,kakom,x,z,dnx,dnz,detJ,a,b,NELT,NODT)
        !*************************************************************
        !Differential of shape function [N]
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),z(1:NODT)
        real(8),intent(in)::a,b
        real(8),intent(out)::dnx(1:4),dnz(1:4)
        real(8),intent(out)::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)

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

        ![dN/dx][dN/dy]
        dnx(1)=( J22*dn1a-J12*dn1b)/detJ
        dnx(2)=( J22*dn2a-J12*dn2b)/detJ
        dnx(3)=( J22*dn3a-J12*dn3b)/detJ
        dnx(4)=( J22*dn4a-J12*dn4b)/detJ
        dnz(1)=(-J21*dn1a+J11*dn1b)/detJ
        dnz(2)=(-J21*dn2a+J11*dn2b)/detJ
        dnz(3)=(-J21*dn3a+J11*dn3b)/detJ
        dnz(4)=(-J21*dn4a+J11*dn4b)/detJ
    end subroutine DNMAT_S4


    subroutine IPOINT4(kk,a,b)
        !*************************************************************
        !Gauss points for integration of element area
        !*************************************************************
        integer,intent(in)::kk
        real(8),intent(out)::a,b

        selectcase(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 STIFF_S4(ne,kakom,x,z,Ak0,alpha,em,Akr,h1vec,eak,NELT,NODT,nt,idan)
        !*************************************************************
        !Element permeability matrix for 4 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,nt,idan
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),z(1:NODT)
        real(8),intent(in)::Ak0(1:NELT),alpha(1:NELT),em(1:NELT)
        real(8),intent(in)::h1vec(1:nt)
        real(8),intent(out)::Akr(1:NELT)
        real(8),intent(out)::eak(1:4,1:4)

        integer::kk,i,j
        integer::k1,k2,k3,k4
        real(8)::sucp,Akre,Akk_el
        real(8)::a,b,detJ
        real(8)::wx(1:4,1:4),wz(1:4,1:4)
        real(8)::dnx(1:4),dnz(1:4)
        real(8)::nm1,nm2,nm3,nm4

        k1=kakom(ne,1)
        k2=kakom(ne,2)
        k3=kakom(ne,3)
        k4=kakom(ne,4)

        Akr(ne)=0.0D0
        do i=1,4
            do j=1,4
                eak(i,j)=0.0D0
            end do
        end do
        do kk=1,4
            call IPOINT4(kk,a,b)
            call DNMAT_S4(ne,kakom,x,z,dnx,dnz,detJ,a,b,NELT,NODT)
            select case(idan)
                case(0)
                    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)
                    sucp=nm1*(h1vec(k1)-z(k1))+nm2*(h1vec(k2)-z(k2))+nm3*(h1vec(k3)-z(k3))+nm4*(h1vec(k4)-z(k4))
                    Akre=CAL_AKR(alpha(ne),em(ne),-sucp)
                case(1)
                    Akre=1.0D0
            end select
            Akr(ne)=Akr(ne)+Akre
            Akk_el=Ak0(ne)*Akre
            do i=1,4
                do j=1,4
                    wx(i,j)=0.0D0
                    wz(i,j)=0.0D0
                end do
            end do
            do i=1,4
                do j=1,4
                    wx(i,j)=wx(i,j)+dnx(i)*dnx(j)
                    wz(i,j)=wz(i,j)+dnz(i)*dnz(j)
                end do
            end do
            do i=1,4
                do j=1,4
                    eak(i,j)=eak(i,j)+Akk_el*(wx(i,j)+wz(i,j))*detJ
                end do
            end do
        end do
        Akr(ne)=0.25D0*Akr(ne)

    end subroutine STIFF_S4


    subroutine CALV_S4(ne,kakom,x,z,hvec,Ak0,alpha,em,Akr,vx,vz,NELT,NODT,nt,idan)
        !*************************************************************
        !Calculation of velocity of 4 node element
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT,nt,idan
        integer,intent(in)::kakom(1:NELT,1:4)
        real(8),intent(in)::x(1:NODT),z(1:NODT)
        real(8),intent(in)::Ak0(1:NELT),alpha(1:NELT),em(1:NELT)
        real(8),intent(in)::hvec(1:nt)
        real(8),intent(out)::Akr(1:NELT),vx(1:NELT),vz(1:NELT)

        integer::j,kk,inod(1:4)
        real(8)::dnx(1:4),dnz(1:4)
        real(8)::a,b,detJ,wvx,wvz
        real(8)::nm1,nm2,nm3,nm4
        real(8)::sucp,Akre

        vx(ne)=0.0D0
        vz(ne)=0.0D0
        Akr(ne)=0.0D0
        do j=1,4
            inod(j)=kakom(ne,j)
        end do
        do kk=1,4
            call IPOINT4(kk,a,b)
            call DNMAT_S4(ne,kakom,x,z,dnx,dnz,detJ,a,b,NELT,NODT)
            select case(idan)
                case(0)
                    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)
                    sucp=nm1*(hvec(inod(1))-z(inod(1)))&
                        +nm2*(hvec(inod(2))-z(inod(2)))&
                        +nm3*(hvec(inod(3))-z(inod(3)))&
                        +nm4*(hvec(inod(4))-z(inod(4)))
                    Akre=CAL_AKR(alpha(ne),em(ne),-sucp)
                case(1)
                    Akre=1.0D0
            end select
            Akr(ne)=Akr(ne)+Akre
            wvx=0.0
            wvz=0.0
            do j=1,4
                wvx=wvx+dnx(j)*hvec(inod(j))
                wvz=wvz+dnz(j)*hvec(inod(j))
            end do
            vx(ne)=vx(ne)+wvx
            vz(ne)=vz(ne)+wvz
        end do
        Akr(ne)=0.25D0*Akr(ne)
        vx(ne)=-0.25*Ak0(ne)*Akr(ne)*vx(ne)
        vz(ne)=-0.25*Ak0(ne)*Akr(ne)*vz(ne)
    end subroutine CALV_S4


    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)
        !-----------------------------------------------------
        !Memory 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_seep
