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