program f90_SEEP4nodSU
    !２次元飽和不飽和定常浸透流解析
    implicit none
    integer::i,j,k,l,k1,k2,k3,n,ia,ja,ne,nnn
    character::strcom*256           !書出用コメント
    integer::NODT                   !節点総数
    integer::NELT                   !要素総数
    integer::MATEL                  !材料種類数
    integer::KOH                    !全水頭指定節点数
    integer::KOQ                    !流量指定節点数
    integer::KOU                    !浸潤境界指定節点数
    integer::KOP                    !浸潤面圧力水頭0点数
    integer::nt                     !全自由度（総節点数×１[１節点自由度]）
    integer::mm                     !温度固定点処理後自由度
    integer::ib                     !バンド幅
    integer::nod=4                  !1要素節点数
    integer::nhen=1                 !1節点自由度
    integer,allocatable::kakom(:,:) !要素構成節点番号
    real(8),allocatable::Ak0(:)     !要素透水係数
    real(8),allocatable::Akr(:)     !要素透水係数
    real(8),allocatable::alpha(:)   !パラメータalpha
    real(8),allocatable::em(:)      !パラメータm
    integer,allocatable::matno(:)   !材料種別No
    real(8),allocatable::wm(:,:)    !材料物性作業用
    real(8),allocatable::x(:)       !節点x座標
    real(8),allocatable::z(:)       !節点y座標
    integer,allocatable::nokh(:)    !全水頭指定節点番号
    integer,allocatable::nokq(:)    !流量指定節点番号
    integer,allocatable::noku(:)    !浸潤境界節点番号
    integer,allocatable::nokp(:)    !浸潤面圧力水頭0点節点番号
    real(8),allocatable::Hinp(:)    !指定全水頭入力
    real(8),allocatable::Qinp(:)    !指定流量入力
    real(8),allocatable::tak(:,:)   !全体剛性行列
    real(8),allocatable::wak(:,:)   !水頭指定節点対応列記憶
    real(8),allocatable::hvec(:)    !全体全水頭ベクトル
    real(8),allocatable::h1vec(:)   !全水頭ベクトル
    real(8),allocatable::h2vec(:)   !全水頭ベクトル
    real(8),allocatable::qvec(:)    !全体節点流量ベクトル
    real(8),allocatable::q0vec(:)   !全体節点流量ベクトル
    real(8),allocatable::reac(:)    !作業用ベクトル
    real(8),allocatable::vx(:)      !x方向流速ベクトル
    real(8),allocatable::vz(:)      !y方向流速ベクトル
    real(8),allocatable::vm(:)      !流速ベクトル
    real(8),allocatable::index(:)   !温度固定点処理用
    real(8)::eak(1:4, 1:4)          !要素透水行列
    real(8)::fact0=1.0E-30          !0判定
    real(8)::elarea
    integer::itmax=100
    integer::icount
    real(8)::QTi,QTo,vmmax,xgm,zgm,vimax,xgi,zgi,vomax,xgo,zgo
    integer::kmmax,kimax,komax

    character :: linebuf*1000
    character :: fmt1*200,fmt2*200,fmt3*200,fmt4*200
    real(4) :: t1,t2
    integer :: iv(1:8)
    character date*8,time1*10,time2*10,zone*5

    !時間計測開始
    call DATE_AND_TIME(date,time1,zone,iv)
    read(time1,*) t1
    !***************************
    !data input
    !***************************
    open(11,file='fnameR.csv',status='old')
    !!--------------------------------------------------------------------------------------
    !!書出用コメント入力
    !!--------------------------------------------------------------------------------------
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
    !節点数，要素数，材料種類数,全水頭指定節点数，流量指定節点数，浸潤境界節点数
    !------------------------------------------------------------------------------------------
    read(11,*) NODT,NELT,MATEL,KOH,KOQ,KOU
    !------------------------------------------------------------------------------------------
    !配列寸法宣言
    !------------------------------------------------------------------------------------------
    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(index(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))
    !------------------------------------------------------------------------------------------
    !材料特性（k0,alpha,m）
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3)
    end do
    !------------------------------------------------------------------------------------------
    !要素構成節点と材料種別（node-1,node-2,node-3,node-4,matno）
    !------------------------------------------------------------------------------------------
    do ne=1,NELT
        read(11,*) kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4),matno(ne)
        do i=1,MATEL
            if(i==matno(ne))then
                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），全水頭初期値
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),z(i),hvec(i)
    end do
    !------------------------------------------------------------------------------------------
    !全水頭指定節点番号，指定全水頭
    !流量指定節点番号，指定流量
    !浸潤境界節点番号
    !------------------------------------------------------------------------------------------
    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)

    !*****************************************************
    !入力確認出力
    !*****************************************************
    do i=1,nt
        qvec(i)=0.0D0
    end do
    if(0<KOH)then
        do k=1,KOH
            hvec(nokh(k))=Hinp(k) !全水頭初期値修正設定
        end do
    end if
    if(0<KOQ)then
        do k=1,KOQ
            qvec(nokq(k))=Qinp(k) !流量初期値設定
        end do
    end if

    fmt1="(i5,',',e15.7,',',e15.7,',',i2,',',i2,',',i2,',',e15.7,',',e15.7)"      !節点入力値
    fmt2="(i5,',',i5,',',i5,',',i5,',',i5,',',e15.7,',',e15.7,',',e15.7,',',i5)"  !要素特性

    open(12, file='fnameW.csv', status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(12,'(a)') 'NODT,NELT,MATEL,KOH,KOQ,KOU'
        write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOH,',',KOQ,',',KOU
        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'
    write(12,'(a)') 'element,node-1,node-2,node-3,node-4,Ak0,alpha,em,matno'
    do ne=1,NELT
            write(linebuf,FMT=fmt2) 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
    !要素面積確認と異常終了通知
    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
        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 do
    close(12)

    !*****************************************************
    !メイン処理計算(収束計算)
    !*****************************************************
    do i=1,nt
        h1vec(i)=hvec(i)
    end do


    do nnn=1,itmax
        !-----------------------------------------------------
        !全体行列クリア
        !-----------------------------------------------------
        do i=1,nt
            reac(i)=0.0D0
            do j=1,nt
                tak(i,j)=0.0D0
            end do
        end do
        !-----------------------------------------------------
        !全体行列組立 [k]は{h1}より計算
        !-----------------------------------------------------
        do ne=1,NELT
            call STIFF_S4(ne,kakom,x,z,Ak0,alpha,em,Akr,h1vec,eak,NELT,NODT,nt)
            call ASMAT(ne,nod,nhen,kakom,eak,tak,NELT,NODT)
        end do
        !-----------------------------------------------------
        !節点流量計算  {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
        !-----------------------------------------------------
        !節点流量更新 {q}->{q0}
        !浸潤境界節点では流入は許容しない
        !(浸潤境界節点での流入量を0セット)
        !-----------------------------------------------------
        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
        !-----------------------------------------------------
        !浸潤境界正圧力保持節点記憶(圧力水頭0セット準備)
        !-----------------------------------------------------
        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
        !-----------------------------------------------------
        !全水頭指定節点対応列記憶，前処理 {q}={q0}-[k]{h(known)}
        !-----------------------------------------------------
        do i=1,nt
            hvec(i)=0.0D0
            qvec(i)=q0vec(i)
            reac(i)=0.0D0
        end do
        !全水頭指定節点での全水頭セット
        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
        !正圧保持浸潤境界節点での圧力水頭0セット
        !(全水頭に位置水頭をセット)
        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
        !-----------------------------------------------------
        !全水頭指定節点処理
        !-----------------------------------------------------
        do i=1,nt
            index(i)=i
            reac(i)=0.0D0
        end do
        if(0<KOH)then
            do i=1,KOH
                n=nokh(i)
                index(n)=0
            end do
        end if
        if(0<KOP)then
            do i=1,KOP
                n=nokp(i)
                index(n)=0
            end do
        end if
        mm=0
        do i=1,nt
            if(index(i)>0)then
                mm=mm+1
                index(mm)=i
            end if
        end do
        do i=1,mm
            ia=index(i)
            reac(i)=qvec(ia)
            do j=1,mm
                ja=index(j)
                tak(i,j)=tak(ia,ja)
            end do
        end do
        !-----------------------------------------------------
        !バンド幅計算とフルマトリックスのバンド化記憶
        !-----------------------------------------------------
        ib=IBAND(mm,tak,fact0,nt)
        call BAND(mm,ib,tak,nt)
        !-----------------------------------------------------
        !対角項確認と異常終了通知
        !-----------------------------------------------------
        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
        !-----------------------------------------------------
        !連立一次方程式（コレスキー法）
        !-----------------------------------------------------
        call CHBAND10(mm,ib,tak,nt)!三角行列
        call CHBAND11(mm,ib,tak,reac,nt)!前進消去・後退代入
        !-----------------------------------------------------
        !全水頭算定 {h}=[k(-1)]{q}
        !-----------------------------------------------------
        do i=1,nt
            hvec(i)=0.0D0
        end do
        do i=1,mm
            ia=index(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
        !-----------------------------------------------------
        !収束判定 |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
        !-----------------------------------------------------
        ! 全水頭更新 {h1}={h2}
        !-----------------------------------------------------
        do i=1,nt
            h1vec(i)=h2vec(i)
        end do

    end do


    !-----------------------------------------------------
    !全流入・流出量算定
    !-----------------------------------------------------
    QTi=0.0D0  !全流入量
    QTo=0.0D0  !全流出量
    do i=1,nt
        if(0<=q0vec(i))QTi=QTi+q0vec(i)
        if(q0vec(i)<0)QTo=QTo+q0vec(i)
    end do
    !-----------------------------------------------------
    !要素平均流速算定・最大流速検索
    !-----------------------------------------------------
    vmmax=0.0D0  !最大流速
    vimax=0.0D0  !流入部最大流速
    vomax=0.0D0  !流出部最大流速
    kmmax=0      !最大流速要素
    kimax=0      !流入部最大流速要素
    komax=0      !流出部最大流速要素
    do ne=1,NELT
        call CALV_S4(ne,kakom,x,z,hvec,Ak0,alpha,em,Akr,vx,vz,NELT,NODT,nt)
        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
    !最大流速節点重心座標
    i=kakom(kmmax,1)
    j=kakom(kmmax,2)
    k=kakom(kmmax,3)
    l=kakom(kmmax,4)
    xgm=0.25D0*(x(i)+x(j)+x(k)+x(l))
    zgm=0.25D0*(z(i)+z(j)+z(k)+z(l))
    i=kakom(kimax,1)
    j=kakom(kimax,2)
    k=kakom(kimax,3)
    l=kakom(kimax,4)
    xgi=0.25D0*(x(i)+x(j)+x(k)+x(l))
    zgi=0.25D0*(z(i)+z(j)+z(k)+z(l))
    i=kakom(komax,1)
    j=kakom(komax,2)
    k=kakom(komax,3)
    l=kakom(komax,4)
    xgo=0.25D0*(x(i)+x(j)+x(k)+x(l))
    zgo=0.25D0*(z(i)+z(j)+z(k)+z(l))

    !*****************************************************
    !結果出力
    !*****************************************************
    fmt3="(i5,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7)"     !節点圧力・流量出力
    fmt4="(i5,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7,',',i5)"   !流速・透水係数出力

    open(12, file='fnameW.csv', access='append')
    !節点番号,x座標,y座標,圧力水頭,節点流量
    !節点流量は，流入が正・流出が負の値
    write(12,'(a)') '*node value'
    write(12,'(a)') 'node,coord-x,coord-z,hvec,qvec,pvec'
    do i=1,NODT
        write(linebuf,FMT=fmt3) i,x(i),z(i),hvec(i),q0vec(i),hvec(i)-z(i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !要素平均流速
    write(12,'(a)') '*element value'
    write(12,'(a)') 'element,vx,vz,vm,Ak0,Akr,matno'
    do ne=1,NELT
        write(linebuf,FMT=fmt4) ne,vx(ne),vz(ne),vm(ne),Ak0(ne),Akr(ne),matno(ne)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    !時間計測完了
    call DATE_AND_TIME(date,time2,zone,iv)
    read(time2,*) t2
    !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=',iv(1),'.',iv(2),'.',iv(3),'_',iv(5),':',iv(6),':',iv(7)
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    close(12)

stop


contains


    real(8) function CAL_AKR(alpha_el,em_el,sucp_el)
        real(8)::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 DNMAT_S4(ne,kakom,x,z,dnx,dnz,detJ,a,b,NELT,NODT)
        integer::ne,NELT,NODT
        integer::kakom(1:NELT,1:4)
        real(8)::x(1:NODT),z(1:NODT)
        real(8)::dnx(1:4),dnz(1:4)
        real(8)::detJ,a,b

        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マトリックスと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/dxdN/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)
        integer::kk
        real(8)::a,b

        !Gauss積分点
        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)
        !要素透水行列
        integer::ne,NELT,NODT,nt
        integer::kakom(1:NELT,1:4)
        real(8)::x(1:NODT),z(1:NODT)
        real(8)::Ak0(1:NELT),alpha(1:NELT),em(1:NELT),Akr(1:NELT)
        real(8)::h1vec(1:nt)
        real(8)::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)
            !積分点圧力水頭
            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)
            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)
    !流速ベクトル算定
        integer::ne,NELT,NODT,nt
        integer::kakom(1:NELT,1:4)
        real(8)::x(1:NODT),z(1:NODT)
        real(8)::Ak0(1:NELT),alpha(1:NELT),em(1:NELT),Akr(1:NELT)
        real(8)::hvec(1:nt)
        real(8)::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)
            !積分点圧力水頭
            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)
            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)
        !剛性行列組立(nod：1要素節点数，nhen：1節点の自由度)
        integer :: NELT,NODT
        integer :: ne,nod,nhen
        integer :: kakom(1:NELT,1:nod)
        real(8) :: sm(1:nod*nhen,1:nod*nhen)
        real(8) :: 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 :: n1,m1,nt
        real(8) :: 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 :: n1,m1,nt
        real(8) :: array(1:nt,1:nt),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 :: mm,nt
        real(8) :: array(1:nt,1:nt),fact0
        integer :: i,j
        !-----------------------------------------------------
        !バンド幅計算
        !-----------------------------------------------------
        IBAND=1
        do i=1,mm-1
            do j=mm,i+1,-1
                if(fact0<=abs(array(i,j)))exit
            end do
            if(IBAND<=j-i+1)IBAND=j-i+1
        end do
        if(mm<IBAND)IBAND=mm
    end function IBAND


    subroutine BAND(mm,ib,array,nt)
        integer :: mm,ib,nt
        real(8) :: array(1:nt,1:nt)
        integer :: i,j,mn
        real(8) :: work(1:nt)
        !-----------------------------------------------------
        !フルマトリックスのバンド化記憶
        !-----------------------------------------------------
        do i=2,mm
            if(i+ib-1<mm)then
                mn=ib
            else
                mn=mm-i+1
            end if
            do j=1,ib
                work(j)=0.0D0
            end do
            do k=1,mn
                work(k)=array(i,i+k-1)
            end do
            do j=1,ib
                array(i,j)=work(j)
            end do
        end do
    end subroutine BAND


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


end program f90_SEEP4nodSU