module defg
    implicit none
    real(8),parameter::g=9.80D0
end module defg

program f90_NUF
    ! Non uniform flow analysis (subcritical flow, supercritical flow)
    use defg
    implicit none
    integer,parameter::nn=200
    integer,parameter::mm=100
    integer::nas(1:nn)
    real(8)::ninp(1:nn)
    real(8)::xx(1:nn)        ! Distance
    real(8)::zz(1:nn)      ! Elevation of river bed
    real(8)::hh(1:nn,1:mm)   ! Water depth
    real(8)::area(1:nn,1:mm) ! Section area of Flow
    real(8)::wper(1:nn,1:mm) ! Wetted perimeter
    real(8)::wids(1:nn,1:mm) ! Width of water surface
    character(len=50)::fnameR
    character(len=10)::dummy1,dummy2,dummy3,dummy4
    integer::nd,i,j,k,nask,it
    integer::icase           ! case
    real(8)::Q0              ! Discharge at Weir
    real(8)::h0              ! Water depth at Weir
    real(8)::DQ(1:nn)        ! Changeof Discharge due to Intake
    real(8)::Q1,n1,x1,a1,r1,t1,h1,z1,v1,hcrit,cs,ds
    real(8)::Q2,n2,x2,a2,r2,t2,h2,z2
    integer,parameter::itlim=100
    real(8)::ha,hb,hc,fa,fb,fc
    real(8),parameter::eps=1.0D-6

    call getarg(1,dummy1)
    call getarg(2,dummy2)
    call getarg(3,dummy3)
    call getarg(4,fnameR)

    read(dummy1,*) icase
    read(dummy2,*) Q0
    read(dummy3,*) h0

    open(11,file=fnameR,status='old')
    read(11,*) nd
    do i=1,nd
        read(11,'()')
        read(11,*) k,xx(i),zz(i),nas(i),ninp(i)
        do j=1,nas(i)
            read(11,*) hh(i,j),area(i,j),wper(i,j),wids(i,j)
        end do
    end do
    close(11)

    if(icase==11.or.icase==21)then
        fnameR='inp_DQ.csv'
        open(11,file=fnameR,status='old')
        do i=1,nd
            read(11,*) xx(i),DQ(i)
        end do
        close(11)
    end if

    k=1
    nask=nas(k)
    select case (icase)
        case (10,20); Q1=Q0
        case (11,21); Q1=Q0+DQ(1)
    end select
    ds=sqrt((xx(k+1)-xx(k))**2.0D0+(zz(k+1)-zz(k))**2.0D0)
    cs=(xx(k+1)-xx(k))/ds
    hcrit=CALHC(icase,k,nn,mm,nask,hh,area,wper,wids,Q1,cs)
    h1=h0
    select case (icase)
    case(10,11)
        if(h1<hcrit)h1=hcrit+0.001D0
    case(20,21)
        if(hcrit<h1)h1=hcrit-0.001D0
    end select
    n1=ninp(k)
    x1=xx(k)
    z1=zz(k)
    call CALART(k,nn,mm,nask,hh,area,wper,wids,h1,a1,r1,t1)
    v1=Q1/a1
    write(6,'(11(a9),a5)') 'x1','Q1','a1','t1','h1','z1','h1+z1','v1','E1','Fr','hc','it'
    write(6,'(11(f9.3))') x1,Q1,a1,t1,h1,z1,h1+z1,v1,h1+z1+v1*v1/2.0D0/g,v1/sqrt(g*a1/t1),hcrit
    do k=2,nd
        nask=nas(k)
        n2=ninp(k)
        select case (icase)
            case (10,20); Q1=Q0; Q2=Q0
            case (11,12); Q1=Q0+DQ(k-1); Q2=Q0+DQ(k)
        end select
        x2=xx(k)
        z2=zz(k)
        ds=sqrt((x2-x1)**2.0D0+(z2-z1)**2.0D0)
        cs=(x2-x1)/ds
        hcrit=CALHC(icase,k,nn,mm,nask,hh,area,wper,wids,Q2,cs)
        select case (icase)
        case(10,11) ! Subcritical flow
            ha=1.0D0*hcrit
            hb=ha+0.5D0*hcrit
            call CALART(k,nn,mm,nask,hh,area,wper,wids,ha,a2,r2,t2)
            fa=CALFF(icase,Q1,Q2,n1,n2,x1,a1,r1,h1,z1,x2,a2,r2,ha,z2)
            do it=1,itlim
                call CALART(k,nn,mm,nask,hh,area,wper,wids,hb,a2,r2,t2)
                fb=CALFF(icase,Q1,Q2,n1,n2,x1,a1,r1,h1,z1,x2,a2,r2,hb,z2)
                if(fa*fb<0.0D0)exit
                hb=hb+0.5D0*hcrit
            end do
        case(20,21) ! Supercritical flow
            ha=0.001D0*hcrit
            hb=1.0D0*hcrit
        end select
        do it=1,itlim
            hc=0.5D0*(ha+hb)
            call CALART(k,nn,mm,nask,hh,area,wper,wids,ha,a2,r2,t2)
            fa=CALFF(icase,Q1,Q2,n1,n2,x1,a1,r1,h1,z1,x2,a2,r2,ha,z2)
            call CALART(k,nn,mm,nask,hh,area,wper,wids,hb,a2,r2,t2)
            fb=CALFF(icase,Q1,Q2,n1,n2,x1,a1,r1,h1,z1,x2,a2,r2,hb,z2)
            call CALART(k,nn,mm,nask,hh,area,wper,wids,hc,a2,r2,t2)
            fc=CALFF(icase,Q1,Q2,n1,n2,x1,a1,r1,h1,z1,x2,a2,r2,hc,z2)
            if(fa*fc<0.0D0)hb=hc
            if(fc*fb<0.0D0)ha=hc
            if(abs(fc)<eps)exit
        end do

        if(itlim<it)cycle

        x1=x2
        a1=a2
        r1=r2
        t1=t2
        h1=hc
        z1=z2
        Q1=Q2
        n1=n2
        v1=Q1/a1
        write(6,'(11(f9.3),i5)') x1,Q1,a1,t1,h1,z1,h1+z1,v1,h1+z1+v1*v1/2.0D0/g,v1/sqrt(g*a1/t1),hcrit,it
    end do

stop

contains

    real(8) function CALFF(icase,Q1,Q2,n1,n2,x1,a1,r1,h1,z1,x2,a2,r2,h2,z2)
        integer,intent(in)::icase
        real(8),intent(in)::Q1,Q2,n1,n2,x1,a1,r1,h1,z1,x2,a2,r2,h2,z2
        real(8)::c1,c2,c3,cs,ds
        CALFF=0.0D0
        ds=sqrt((x2-x1)**2.0D0+(z2-z1)**2.0D0)
        cs=(x2-x1)/ds
        c2=h2*cs + z2 + 1.0D0/2.0D0/g * Q2*Q2/a2/a2
        c1=h1*cs + z1 + 1.0D0/2.0D0/g * Q1*Q1/a1/a1
        c3=0.5D0 * (n1*n1*Q1*Q1/r1**(4.0D0/3.0D0)/a1/a1 + n2*n2*Q2*Q2/r2**(4.0D0/3.0D0)/a2/a2) * ds
        if(icase==10.or.icase==11)CALFF=c2-c1-c3 ! Subcritical flow
        if(icase==20.or.icase==21)CALFF=c1-c2-c3 ! Supercritical flow
    end function CALFF

    real(8) function CALY(x1,y1,x2,y2,xx)
        real(8),intent(in)::x1,y1,x2,y2,xx
        real(8)::aa,bb
        aa=(y2-y1)/(x2-x1)
        bb=(x2*y1-x1*y2)/(x2-x1)
        CALY=aa*xx+bb
    end function CALY

    subroutine CALART(k,nn,mm,nask,hh,area,wper,wids,wd,aa,rr,tt)
        integer,intent(in)::k,nn,mm,nask
        real(8),intent(in)::hh(1:nn,1:mm),area(1:nn,1:mm),wper(1:nn,1:mm),wids(1:nn,1:mm),wd
        real(8),intent(out)::aa,rr,tt
        integer::i
        do i=1,nask-1
            if(hh(k,i)<=wd.and.wd<=hh(k,i+1))then
                aa=CALY(hh(k,i),area(k,i),hh(k,i+1),area(k,i+1),wd)
                rr=aa/CALY(hh(k,i),wper(k,i),hh(k,i+1),wper(k,i+1),wd)
                tt=CALY(hh(k,i),wids(k,i),hh(k,i+1),wids(k,i+1),wd)
            end if
        end do
    end subroutine CALART

    real(8) function CALHC(icase,k,nn,mm,nask,hh,area,wper,wids,Q,cs)
        integer,intent(in)::icase,k,nn,mm,nask
        real(8),intent(in)::hh(1:nn,1:mm),area(1:nn,1:mm),wper(1:nn,1:mm),wids(1:nn,1:mm),Q,cs
        real(8)::aa,rr,tt
        real(8)::ha,hb,hc,fa,fb,fc
        real(8)::dh=0.01D0
        integer::limit=100
        real(8)::eps=1.0D-6

        ha=0.01D0
        hb=ha+dh
        call CALART(k,nn,mm,nask,hh,area,wper,wids,ha,aa,rr,tt)
        fa=CALFU(Q,aa,tt,cs)
        do i=1,limit*100
            call CALART(k,nn,mm,nask,hh,area,wper,wids,hb,aa,rr,tt)
            fb=CALFU(Q,aa,tt,cs)
            if(fa*fb<0.0D0)exit
            hb=hb+dh
        end do
        do i=1,limit
            hc=0.5D0*(ha+hb)
            call CALART(k,nn,mm,nask,hh,area,wper,wids,ha,aa,rr,tt)
            fa=CALFU(Q,aa,tt,cs)
            call CALART(k,nn,mm,nask,hh,area,wper,wids,hb,aa,rr,tt)
            fb=CALFU(Q,aa,tt,cs)
            call CALART(k,nn,mm,nask,hh,area,wper,wids,hc,aa,rr,tt)
            fc=CALFU(Q,aa,tt,cs)
            if(fa*fc<0.0D0)hb=hc
            if(fc*fb<0.0D0)ha=hc
            if(abs(fc)<eps)exit
        end do
        CALHC=hc
    end function CALHC

    real(8) function CALFU(Q,a,t,cs)
        real(8),intent(in)::Q,a,t,cs
        CALFU=a*a*a/t-Q*Q/g/cs
    end function CALFU


end program f90_NUF