! ============================================================ ! 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