program submerged_orifice_weir
    implicit none
    integer :: weirType, i, iostat_val
    double precision :: H, P, b, theta_deg, theta_rad, Cd_user, g, tailwater
    double precision :: Cd, Q, V_approach, E_specific, A_flow, Fr
    double precision :: H_i, Q_i, Cd_i
    double precision, parameter :: PI = 3.141592653589793d0

    read(*,*,iostat=iostat_val) weirType
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid weir type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) H
    read(*,*,iostat=iostat_val) P
    read(*,*,iostat=iostat_val) b
    read(*,*,iostat=iostat_val) theta_deg
    read(*,*,iostat=iostat_val) Cd_user
    read(*,*,iostat=iostat_val) g
    read(*,*,iostat=iostat_val) tailwater
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all weir/orifice inputs.'
        stop
    end if
    if (H < 0.0d0 .or. P < 0.0d0 .or. b < 0.0d0) then
        write(*,*) 'ERROR: Head, crest height and width must be non-negative.'
        stop
    end if
    if (g <= 0.0d0) g = 9.81d0
    theta_rad = theta_deg * PI / 180.0d0

    select case(weirType)
    case(1)
        ! Rectangular sharp-crested weir (Rehbock)
        if (Cd_user > 0.01d0) then
            Cd = Cd_user
        else
            Cd = rehbock_cd(H, P)
        end if
        Q = Cd * (2.0d0/3.0d0) * b * sqrt(2.0d0*g) * H**1.5d0
    case(2)
        ! V-notch (triangular) weir
        if (Cd_user > 0.01d0) then
            Cd = Cd_user
        else
            Cd = 0.58d0
        end if
        Q = Cd * (8.0d0/15.0d0) * sqrt(2.0d0*g) * tan(theta_rad/2.0d0) * H**2.5d0
    case(3)
        ! Trapezoidal (Cipolletti) weir
        if (Cd_user > 0.01d0) then
            Cd = Cd_user
        else
            Cd = 0.63d0
        end if
        Q = Cd * (2.0d0/3.0d0) * sqrt(2.0d0*g) * H**1.5d0 * &
            (b + H*tan(theta_rad/2.0d0))
    case(4)
        ! Submerged sluice gate orifice
        if (Cd_user > 0.01d0) then
            Cd = Cd_user
        else
            Cd = 0.61d0
        end if
        ! H = upstream head above gate bottom, P = gate opening height
        ! tailwater = downstream depth
        if (tailwater > 0.0d0 .and. tailwater < H) then
            Q = Cd * b * P * sqrt(2.0d0*g*(H - tailwater))
        else
            Q = Cd * b * P * sqrt(2.0d0*g*H)
        end if
    case default
        write(*,*) 'ERROR: Invalid weir type. Use 1-4.'
        stop
    end select

    ! Approach velocity and specific energy
    if (b > 0.0d0 .and. H > 0.0d0) then
        A_flow = b * (H + P)
        V_approach = Q / max(A_flow, 1.0d-30)
    else
        V_approach = 0.0d0
    end if
    E_specific = H + V_approach**2 / (2.0d0*g)
    if (H > 0.0d0 .and. b > 0.0d0) then
        Fr = V_approach / sqrt(g * H)
    else
        Fr = 0.0d0
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   SUBMERGED ORIFICE & WEIR FLOW ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,I8)')       '  Weir Type                = ', weirType
    write(*,'(A,ES12.4,A)') '  Head H                   = ', H, ' m'
    write(*,'(A,ES12.4,A)') '  Crest Height P           = ', P, ' m'
    write(*,'(A,ES12.4,A)') '  Width b                  = ', b, ' m'
    write(*,'(A,ES12.4,A)') '  Notch Angle              = ', theta_deg, ' deg'
    write(*,'(A,ES12.4,A)') '  Tailwater Depth          = ', tailwater, ' m'
    write(*,'(A,ES12.4,A)') '  Gravity                  = ', g, ' m/s2'
    write(*,*)
    write(*,'(A)') '--- FLOW RESULTS --------------------------------------------'
    write(*,'(A,ES12.4)')   '  Discharge Coefficient Cd = ', Cd
    write(*,'(A,ES12.4,A)') '  Flow Rate Q              = ', Q, ' m3/s'
    write(*,'(A,ES12.4,A)') '  Approach Velocity        = ', V_approach, ' m/s'
    write(*,'(A,ES12.4,A)') '  Specific Energy          = ', E_specific, ' m'
    write(*,'(A,ES12.4)')   '  Froude Number            = ', Fr
    write(*,*)

    write(*,'(A)') '--- Q VS HEAD SWEEP -----------------------------------------'
    write(*,'(A)') '  H[m]          Q[m3/s]       Cd'
    write(*,'(A)') '  -------------------------------------------'
    do i = 1, 50
        H_i = 0.01d0 + 2.0d0 * dble(i-1) / 49.0d0
        select case(weirType)
        case(1)
            Cd_i = rehbock_cd(H_i, P)
            Q_i = Cd_i*(2.0d0/3.0d0)*b*sqrt(2.0d0*g)*H_i**1.5d0
        case(2)
            Cd_i = 0.58d0
            Q_i = Cd_i*(8.0d0/15.0d0)*sqrt(2.0d0*g)*tan(theta_rad/2.0d0)*H_i**2.5d0
        case(3)
            Cd_i = 0.63d0
            Q_i = Cd_i*(2.0d0/3.0d0)*sqrt(2.0d0*g)*H_i**1.5d0*(b+H_i*tan(theta_rad/2.0d0))
        case(4)
            Cd_i = 0.61d0
            Q_i = Cd_i*b*P*sqrt(2.0d0*g*H_i)
        end select
        write(*,'(F10.4,2X,ES12.4,2X,ES12.4)') H_i, Q_i, Cd_i
    end do
    write(*,*)

    write(*,'(A)') '--- SPECIFIC ENERGY SWEEP -----------------------------------'
    write(*,'(A)') '  H[m]          E[m]          Fr'
    write(*,'(A)') '  -------------------------------------------'
    do i = 1, 40
        H_i = 0.02d0 + 2.5d0 * dble(i-1) / 39.0d0
        select case(weirType)
        case(1)
            Q_i = rehbock_cd(H_i,P)*(2.0d0/3.0d0)*b*sqrt(2.0d0*g)*H_i**1.5d0
        case(2)
            Q_i = 0.58d0*(8.0d0/15.0d0)*sqrt(2.0d0*g)*tan(theta_rad/2.0d0)*H_i**2.5d0
        case(3)
            Q_i = 0.63d0*(2.0d0/3.0d0)*sqrt(2.0d0*g)*H_i**1.5d0*(b+H_i*tan(theta_rad/2.0d0))
        case(4)
            Q_i = 0.61d0*b*P*sqrt(2.0d0*g*H_i)
        end select
        if (b > 0.0d0 .and. (H_i+P) > 0.0d0) then
            V_approach = Q_i/(b*(H_i+P))
        else
            V_approach = 0.0d0
        end if
        E_specific = H_i + V_approach**2/(2.0d0*g)
        Fr = V_approach / sqrt(g*max(H_i,1.0d-30))
        write(*,'(F10.4,2X,ES12.4,2X,ES12.4)') H_i, E_specific, Fr
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Rectangular weir: Q = Cd (2/3) b sqrt(2g) H^1.5.'
    write(*,'(A)') '  V-notch: Q = Cd (8/15) sqrt(2g) tan(theta/2) H^2.5.'
    write(*,'(A)') '  Trapezoidal: Q = Cd (2/3) sqrt(2g) H^1.5 (b + H tan(theta/2)).'
    write(*,'(A)') '  Sluice gate: Q = Cd b a sqrt(2g (H - Htail)).'
    write(*,'(A)') '  Rehbock Cd: Cd = 0.602 + 0.083 H/P.'

contains
    double precision function rehbock_cd(Hin, Pin)
        implicit none
        double precision, intent(in) :: Hin, Pin
        if (Pin > 0.0d0) then
            rehbock_cd = 0.602d0 + 0.083d0 * Hin / Pin
        else
            rehbock_cd = 0.62d0
        end if
        if (rehbock_cd > 0.85d0) rehbock_cd = 0.85d0
        if (rehbock_cd < 0.55d0) rehbock_cd = 0.55d0
    end function rehbock_cd
end program submerged_orifice_weir