💻 Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
Manning Open Channel Uniform Flow
Core Numerical Engine in Fortran 90 • 0 total downloads
manning_equation.f90
! =========================================================================
! Source File: manning_equation.f90
! =========================================================================
! ============================================================
! 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
Solver Description
Solves for normal depth ($y_n$) in open channels. Formulates Manning's equation ($Q = \frac{1}{n} A R_h^{2/3} S_0^{1/2}$) and iteratively solves the non-linear equation for flow depth ($y$) using the Newton-Raphson method.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -ffree-line-length-none manning_equation.f90 -o manning_equation_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
manning_equation_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Channel Geometry Code (1=Rectangular, 2=Trapezoidal, 3=Triangular, 4=Circular)
2
! Width ($b$) / Diameter ($D$) [m]
6.0
! Side Slope ($z$) [H:V] (ignored for rect/circular)
2.0
! Manning n Roughness
0.025
! Bed Slope ($S_0$) [m/m]
0.001
! Flow Discharge ($Q$) [m³/s]
30.0
2
! Width ($b$) / Diameter ($D$) [m]
6.0
! Side Slope ($z$) [H:V] (ignored for rect/circular)
2.0
! Manning n Roughness
0.025
! Bed Slope ($S_0$) [m/m]
0.001
! Flow Discharge ($Q$) [m³/s]
30.0