! ============================================================
! ThermoFluidCalc — Manning's Equation Open Channel Flow
! Solve for Q (discharge) or y_n (normal depth)
! Supports: rectangular, trapezoidal, triangular, circular
! ============================================================
program manning_equation
    implicit none

    real(8), parameter :: g = 9.81d0
    real(8), parameter :: PI = 3.141592653589793d0

    integer :: shape, solve_mode
    real(8) :: b, z, D, mann_n, S0, known_val
    real(8) :: A, P_wet, R_h, T, y_h, V, Q, Fr, E, y_c, y_n, Re
    real(8) :: nu_water
    integer :: i, n_points
    real(8) :: y_plot, Q_plot, y_max
    character(len=20) :: shape_name, solve_name, regime_name

    ! Read inputs
    read(*,*) shape        ! 1=rect, 2=trap, 3=tri, 4=circular
    read(*,*) b            ! bottom width [m] (rect/trap)
    read(*,*) z            ! side slope H:V (trap/tri)
    read(*,*) D            ! diameter [m] (circular)
    read(*,*) mann_n       ! Manning's n
    read(*,*) S0           ! bed slope
    read(*,*) solve_mode   ! 1=find Q from y, 2=find y from Q
    read(*,*) known_val    ! depth y [m] or discharge Q [m³/s]

    nu_water = 1.0d-6  ! kinematic viscosity of water [m²/s]

    ! Map shape names
    select case (shape)
    case (1)
        shape_name = 'Rectangular'
    case (2)
        shape_name = 'Trapezoidal'
    case (3)
        shape_name = 'Triangular'
    case (4)
        shape_name = 'Circular'
    case default
        shape_name = 'Unknown'
    end select

    ! Map solve mode names
    if (solve_mode == 1) then
        solve_name = 'Q from Depth y_n'
    else
        solve_name = 'Depth y_n from Q'
    end if

    if (solve_mode == 1) then
        ! === MODE 1: Given depth y, find Q ===
        y_n = known_val
        call geometry(shape, b, z, D, y_n, A, P_wet, R_h, T)
        if (A <= 0.0d0 .or. P_wet <= 0.0d0) then
            write(*,'(A)') "ERROR: Invalid geometry or zero depth"
            stop
        end if
        Q   = (1.0d0/mann_n) * A * R_h**(2.0d0/3.0d0) * sqrt(S0)
        V   = Q / A
        y_h = A / T
        Fr  = V / sqrt(g * y_h)
        E   = y_n + V**2 / (2.0d0*g)
        Re  = V * R_h / nu_water
        ! Critical depth by bisection
        call critical_depth(shape, b, z, D, Q, y_c)
    else
        ! === MODE 2: Given Q, find normal depth y_n by bisection ===
        Q = known_val
        call normal_depth(shape, b, z, D, mann_n, S0, Q, y_n)
        call geometry(shape, b, z, D, y_n, A, P_wet, R_h, T)
        V   = Q / A
        y_h = A / T
        Fr  = V / sqrt(g * y_h)
        E   = y_n + V**2 / (2.0d0*g)
        Re  = V * R_h / nu_water
        call critical_depth(shape, b, z, D, Q, y_c)
    end if

    ! Determine flow regime
    if (Fr < 0.95d0) then
        regime_name = 'Subcritical'
    else if (Fr > 1.05d0) then
        regime_name = 'Supercritical'
    else
        regime_name = 'Critical'
    end if

    ! === OUTPUT ASCII REPORT ===
    write(*,'(A)') '============================================================'
    write(*,'(A)') '      MANNINGS EQUATION OPEN CHANNEL FLOW CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT CONFIGURATION ---'
    write(*,'(A,A)')       '  Channel Shape           = ', trim(shape_name)
    select case (shape)
    case (1)
        write(*,'(A,F12.4,A)') '  Bottom Width (b)        = ', b, ' m'
    case (2)
        write(*,'(A,F12.4,A)') '  Bottom Width (b)        = ', b, ' m'
        write(*,'(A,F12.4)')   '  Side Slope (z)          = ', z
    case (3)
        write(*,'(A,F12.4)')   '  Side Slope (z)          = ', z
    case (4)
        write(*,'(A,F12.4,A)') '  Pipe Diameter (D)       = ', D, ' m'
    end select
    write(*,'(A,F12.4)')     '  Mannings Roughness (n)  = ', mann_n
    write(*,'(A,F12.6,A)')   '  Bed Slope (S0)          = ', S0, ' m/m'
    write(*,'(A,A)')       '  Calculation Mode        = ', trim(solve_name)
    write(*,*)
    write(*,'(A)') '--- HYDRAULIC RESULTS ---'
    write(*,'(A,F12.4,A)')   '  Normal Depth (y_n)      = ', y_n, ' m'
    write(*,'(A,F12.4,A)')   '  Flow Area (A)           = ', A, ' m2'
    write(*,'(A,F12.4,A)')   '  Wetted Perimeter (P)    = ', P_wet, ' m'
    write(*,'(A,F12.4,A)')   '  Hydraulic Radius (R_h)  = ', R_h, ' m'
    write(*,'(A,F12.4,A)')   '  Top Width (T)           = ', T, ' m'
    write(*,'(A,F12.4,A)')   '  Hydraulic Depth (y_h)   = ', y_h, ' m'
    write(*,'(A,F12.4,A)')   '  Mean Velocity (V)       = ', V, ' m/s'
    write(*,'(A,F12.4,A)')   '  Discharge (Q)           = ', Q, ' m3/s'
    write(*,'(A,F12.4)')     '  Froude Number (Fr)      = ', Fr
    write(*,'(A,F12.4,A)')   '  Specific Energy (E)     = ', E, ' m'
    write(*,'(A,ES12.4)')    '  Reynolds Number (Re)    = ', Re
    write(*,'(A,F12.4,A)')   '  Critical Depth (y_c)    = ', y_c, ' m'
    write(*,'(A,A)')         '  Flow Regime             = ', trim(regime_name)
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ---'
    write(*,'(A)') '  Manning Equation: Q = (1/n) * A * R_h^(2/3) * sqrt(S0)'
    write(*,'(A)') '  Hydraulic Radius: R_h = A / P'
    write(*,'(A)') '  Froude Number:    Fr = V / sqrt(g * A / T)'
    write(*,'(A)') '  Specific Energy:  E = y + V^2 / (2g)'
    write(*,'(A)') '  Critical Depth:   Solved by A^3 / T = Q^2 / g'
    write(*,'(A)') '============================================================'
    write(*,*)

    ! === CHART DATA (Rating curve Q vs y) ===
    write(*,'(A)') "--- CHART_DATA ---"
    if (shape == 4) then
        y_max = D * 0.95d0
    else
        y_max = max(y_n * 2.0d0, y_c * 2.0d0)
        if (y_max < 0.1d0) y_max = 0.5d0
    end if
    n_points = 30
    do i = 1, n_points
        y_plot = y_max * real(i, 8) / real(n_points, 8)
        call geometry(shape, b, z, D, y_plot, A, P_wet, R_h, T)
        if (A > 0.0d0 .and. P_wet > 0.0d0) then
            Q_plot = (1.0d0/mann_n) * A * R_h**(2.0d0/3.0d0) * sqrt(S0)
        else
            Q_plot = 0.0d0
        end if
        write(*,'(F10.4,A,F10.4)') y_plot, ",", Q_plot
    end do
    write(*,'(A)') "--- END_CHART_DATA ---"

contains

    ! ---- Cross-section geometry ----
    subroutine geometry(shape, b, z, D, y, A, P, R, T)
        integer, intent(in)  :: shape
        real(8), intent(in)  :: b, z, D, y
        real(8), intent(out) :: A, P, R, T
        real(8) :: theta

        select case (shape)
        case (1)  ! Rectangular
            A = b * y
            P = b + 2.0d0*y
            T = b
        case (2)  ! Trapezoidal
            A = (b + z*y) * y
            P = b + 2.0d0*y*sqrt(1.0d0 + z**2)
            T = b + 2.0d0*z*y
        case (3)  ! Triangular
            A = z * y**2
            P = 2.0d0*y*sqrt(1.0d0 + z**2)
            T = 2.0d0*z*y
        case (4)  ! Circular
            if (y >= D) then
                theta = 2.0d0*PI
            else
                theta = 2.0d0 * acos(1.0d0 - 2.0d0*y/D)
            end if
            A = (D**2/8.0d0) * (theta - sin(theta))
            P = D * theta / 2.0d0
            T = D * sin(theta/2.0d0)
        case default
            A = b * y; P = b + 2.0d0*y; T = b
        end select

        if (P > 0.0d0) then
            R = A / P
        else
            R = 0.0d0
        end if
    end subroutine geometry

    ! ---- Normal depth by bisection ----
    subroutine normal_depth(shape, b, z, D, n_mann, S0, Q_target, yn)
        integer, intent(in)  :: shape
        real(8), intent(in)  :: b, z, D, n_mann, S0, Q_target
        real(8), intent(out) :: yn
        real(8) :: y_lo, y_hi, y_mid, Q_mid, A_, P_, R_, T_
        integer :: iter

        y_lo = 1.0d-4
        if (shape == 4) then
            y_hi = D * 0.99d0
        else
            y_hi = 50.0d0
        end if

        do iter = 1, 100
            y_mid = (y_lo + y_hi) / 2.0d0
            call geometry(shape, b, z, D, y_mid, A_, P_, R_, T_)
            if (A_ > 0.0d0 .and. P_ > 0.0d0) then
                Q_mid = (1.0d0/n_mann) * A_ * R_**(2.0d0/3.0d0) * sqrt(S0)
            else
                Q_mid = 0.0d0
            end if
            if (Q_mid < Q_target) then
                y_lo = y_mid
            else
                y_hi = y_mid
            end if
            if ((y_hi - y_lo) < 1.0d-6) exit
        end do
        yn = (y_lo + y_hi) / 2.0d0
    end subroutine normal_depth

    ! ---- Critical depth by bisection: A³/T = Q²/g ----
    subroutine critical_depth(shape, b, z, D, Q, yc)
        integer, intent(in)  :: shape
        real(8), intent(in)  :: b, z, D, Q
        real(8), intent(out) :: yc
        real(8) :: y_lo, y_hi, y_mid, A_, P_, R_, T_, f_mid
        integer :: iter

        y_lo = 1.0d-4
        if (shape == 4) then
            y_hi = D * 0.99d0
        else
            y_hi = 50.0d0
        end if

        do iter = 1, 100
            y_mid = (y_lo + y_hi) / 2.0d0
            call geometry(shape, b, z, D, y_mid, A_, P_, R_, T_)
            if (T_ > 0.0d0) then
                f_mid = A_**3 / T_ - Q**2 / g
            else
                f_mid = -Q**2 / g
            end if
            if (f_mid < 0.0d0) then
                y_lo = y_mid
            else
                y_hi = y_mid
            end if
            if ((y_hi - y_lo) < 1.0d-6) exit
        end do
        yc = (y_lo + y_hi) / 2.0d0
    end subroutine critical_depth

end program manning_equation
