๐ป Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Gradually Varied Flow Profile
Core Numerical Engine in Fortran 90 โข 23 total downloads
gvf_profile.f90
! =========================================================================
! Source File: gvf_profile.f90
! =========================================================================
program gvf_profile
implicit none
integer :: i, n_steps, channelType, iostat_val, profileClass
double precision :: Q, b, S0, n_mann, z_side, y_start, x_length, g
double precision :: yc, yn, A, P_wet, Rh, Sf, Fr2, dydx
double precision :: x, y, dx, y_new, A2, P2, Rh2, Sf2, Fr2b, dydx2
double precision :: Ac, Pc, Rhc, An, Pn, Rhn
character(len=40) :: className, channelName
read(*,*,iostat=iostat_val) channelType
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid channel type input.'
stop
end if
read(*,*,iostat=iostat_val) Q
read(*,*,iostat=iostat_val) b
read(*,*,iostat=iostat_val) S0
read(*,*,iostat=iostat_val) n_mann
read(*,*,iostat=iostat_val) z_side
read(*,*,iostat=iostat_val) y_start
read(*,*,iostat=iostat_val) x_length
read(*,*,iostat=iostat_val) g
read(*,*,iostat=iostat_val) n_steps
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all GVF inputs.'
stop
end if
if (Q <= 0.0d0 .or. b <= 0.0d0 .or. n_mann <= 0.0d0) then
write(*,*) 'ERROR: Q, b, and Manning n must be positive.'
stop
end if
if (y_start <= 0.0d0) then
write(*,*) 'ERROR: Starting depth must be positive.'
stop
end if
if (g <= 0.0d0) g = 9.81d0
if (n_steps < 10) n_steps = 500
if (x_length <= 0.0d0) x_length = 1000.0d0
select case(channelType)
case(1)
channelName = 'Rectangular'
case(2)
channelName = 'Trapezoidal'
case default
channelName = 'Rectangular'
channelType = 1
end select
! Critical depth by bisection
yc = find_critical_depth(Q, b, z_side, g, channelType)
! Normal depth by bisection
yn = find_normal_depth(Q, b, z_side, S0, n_mann, channelType)
! Classify profile
call classify_profile(S0, y_start, yn, yc, profileClass, className)
write(*,'(A)') '============================================================'
write(*,'(A)') ' GRADUALLY VARIED FLOW (GVF) PROFILE ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Channel Shape = ', trim(channelName)
write(*,'(A,ES12.4,A)') ' Discharge Q = ', Q, ' m3/s'
write(*,'(A,ES12.4,A)') ' Bottom Width b = ', b, ' m'
write(*,'(A,ES12.4)') ' Bed Slope S0 = ', S0
write(*,'(A,ES12.4)') ' Manning n = ', n_mann
write(*,'(A,ES12.4)') ' Side Slope z = ', z_side
write(*,'(A,ES12.4,A)') ' Starting Depth y_start = ', y_start, ' m'
write(*,'(A,ES12.4,A)') ' Reach Length = ', x_length, ' m'
write(*,'(A,I8)') ' Number of Steps = ', n_steps
write(*,*)
write(*,'(A)') '--- CRITICAL / NORMAL DEPTHS --------------------------------'
write(*,'(A,ES12.4,A)') ' Critical Depth yc = ', yc, ' m'
write(*,'(A,ES12.4,A)') ' Normal Depth yn = ', yn, ' m'
write(*,'(A,A)') ' Profile Classification = ', trim(className)
write(*,*)
! Compute specific energy and Froude at critical depth
call channel_geometry(channelType, yc, b, z_side, Ac, Pc, Rhc)
call channel_geometry(channelType, yn, b, z_side, An, Pn, Rhn)
write(*,'(A)') '--- HYDRAULIC DATA ------------------------------------------'
write(*,'(A,ES12.4,A)') ' Critical Area Ac = ', Ac, ' m2'
write(*,'(A,ES12.4,A)') ' Normal Area An = ', An, ' m2'
write(*,'(A,ES12.4,A)') ' Critical Rh = ', Rhc, ' m'
write(*,'(A,ES12.4,A)') ' Normal Rh = ', Rhn, ' m'
write(*,*)
! Euler integration of dy/dx = (S0 - Sf)/(1 - Fr^2)
dx = x_length / dble(n_steps)
x = 0.0d0
y = y_start
write(*,'(A)') '--- WATER SURFACE PROFILE -----------------------------------'
write(*,'(A)') ' x[m] y[m] Sf Fr E[m]'
write(*,'(A)') ' ----------------------------------------------------------------'
do i = 0, n_steps
call channel_geometry(channelType, y, b, z_side, A, P_wet, Rh)
if (A > 0.0d0) then
Sf = (Q * n_mann / (A * Rh**(2.0d0/3.0d0)))**2
Fr2 = Q**2 * (b + 2.0d0*z_side*y) / (g * A**3)
else
Sf = S0
Fr2 = 0.0d0
end if
write(*,'(F12.3,2X,F12.6,2X,ES12.4,2X,F10.5,2X,F12.6)') &
x, y, Sf, sqrt(max(Fr2,0.0d0)), &
y + Q**2/(2.0d0*g*max(A,1.0d-30)**2)
if (i == n_steps) exit
! Euler predictor
if (abs(1.0d0 - Fr2) > 1.0d-8) then
dydx = (S0 - Sf) / (1.0d0 - Fr2)
else
dydx = 0.0d0
end if
! Heun corrector
y_new = y + dydx * dx
if (y_new <= 0.001d0) y_new = 0.001d0
call channel_geometry(channelType, y_new, b, z_side, A2, P2, Rh2)
if (A2 > 0.0d0) then
Sf2 = (Q * n_mann / (A2 * Rh2**(2.0d0/3.0d0)))**2
Fr2b = Q**2 * (b + 2.0d0*z_side*y_new) / (g * A2**3)
else
Sf2 = S0
Fr2b = 0.0d0
end if
if (abs(1.0d0 - Fr2b) > 1.0d-8) then
dydx2 = (S0 - Sf2) / (1.0d0 - Fr2b)
else
dydx2 = 0.0d0
end if
y = y + 0.5d0*(dydx + dydx2)*dx
if (y <= 0.001d0) y = 0.001d0
x = x + dx
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' dy/dx = (S0 - Sf) / (1 - Fr^2).'
write(*,'(A)') ' Sf = [Q n / (A Rh^(2/3))]^2 (Manning).'
write(*,'(A)') ' Fr^2 = Q^2 T / (g A^3); T = b + 2 z y.'
write(*,'(A)') ' Heun (improved Euler) integration scheme.'
contains
subroutine channel_geometry(ctype, depth, bw, zs, Aout, Pout, Rhout)
implicit none
integer, intent(in) :: ctype
double precision, intent(in) :: depth, bw, zs
double precision, intent(out) :: Aout, Pout, Rhout
if (ctype == 2) then
Aout = (bw + zs*depth)*depth
Pout = bw + 2.0d0*depth*sqrt(1.0d0 + zs**2)
else
Aout = bw * depth
Pout = bw + 2.0d0*depth
end if
if (Pout > 0.0d0) then
Rhout = Aout / Pout
else
Rhout = 0.0d0
end if
end subroutine channel_geometry
double precision function find_critical_depth(Qin, bw, zs, gin, ctype)
implicit none
double precision, intent(in) :: Qin, bw, zs, gin
integer, intent(in) :: ctype
double precision :: lo, hi, mid, Amid, Tmid, fval
integer :: it
lo = 0.0001d0
hi = 50.0d0
do it = 1, 120
mid = 0.5d0*(lo+hi)
if (ctype == 2) then
Amid = (bw + zs*mid)*mid
Tmid = bw + 2.0d0*zs*mid
else
Amid = bw*mid
Tmid = bw
end if
fval = Qin**2*Tmid - gin*Amid**3
if (fval > 0.0d0) then
lo = mid
else
hi = mid
end if
end do
find_critical_depth = 0.5d0*(lo+hi)
end function find_critical_depth
double precision function find_normal_depth(Qin, bw, zs, S0in, nin, ctype)
implicit none
double precision, intent(in) :: Qin, bw, zs, S0in, nin
integer, intent(in) :: ctype
double precision :: lo, hi, mid, Amid, Pmid, Rhmid, Qmann, fval
integer :: it
if (S0in <= 0.0d0) then
find_normal_depth = 1.0d30
return
end if
lo = 0.0001d0
hi = 50.0d0
do it = 1, 120
mid = 0.5d0*(lo+hi)
if (ctype == 2) then
Amid = (bw + zs*mid)*mid
Pmid = bw + 2.0d0*mid*sqrt(1.0d0 + zs**2)
else
Amid = bw*mid
Pmid = bw + 2.0d0*mid
end if
Rhmid = Amid / max(Pmid, 1.0d-30)
Qmann = (1.0d0/nin)*Amid*Rhmid**(2.0d0/3.0d0)*sqrt(S0in)
fval = Qmann - Qin
if (fval < 0.0d0) then
lo = mid
else
hi = mid
end if
end do
find_normal_depth = 0.5d0*(lo+hi)
end function find_normal_depth
subroutine classify_profile(S0in, ystart, ynorm, ycrit, pclass, pname)
implicit none
double precision, intent(in) :: S0in, ystart, ynorm, ycrit
integer, intent(out) :: pclass
character(len=40), intent(out) :: pname
if (S0in <= 0.0d0) then
! Adverse or horizontal
if (ystart > ycrit) then
pclass = 10; pname = 'A2 or H2 - Adverse/Horizontal'
else
pclass = 11; pname = 'A3 or H3 - Adverse/Horizontal'
end if
else if (ynorm > ycrit) then
! Mild slope
if (ystart > ynorm) then
pclass = 1; pname = 'M1 - Mild backwater'
else if (ystart > ycrit) then
pclass = 2; pname = 'M2 - Mild drawdown'
else
pclass = 3; pname = 'M3 - Mild supercritical'
end if
else if (abs(ynorm-ycrit) < 1.0d-6) then
pclass = 7; pname = 'C1/C3 - Critical slope'
else
! Steep slope
if (ystart > ycrit) then
pclass = 4; pname = 'S1 - Steep backwater'
else if (ystart > ynorm) then
pclass = 5; pname = 'S2 - Steep drawdown'
else
pclass = 6; pname = 'S3 - Steep supercritical'
end if
end if
end subroutine classify_profile
end program gvf_profile
Solver Description
Simulate gradually varied open-channel water surface profiles (M1, M2, S1, S2 etc.) using Heun integration of the GVF equation.
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 -O3 gvf_profile.f90 -o gvf_profile
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
gvf_profile < input.txt
๐ฅ Downloads & Local Files
Preview of the required input file (input.txt):
! ct_init\nDischarge Q [mรยณ/s]\nBottom width b [m]\nBed slope Sรขโโฌ\nManning n\nSide slope z (H:V, 0 = rectangular)\nys_init\nxl_init\nGravity g [m/sรยฒ]\nManning n
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0
! Parameter 10
0.0
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0
! Parameter 10
0.0