program lid_driven_cavity_re
    implicit none
    integer :: iostat_val
    double precision :: L_cavity_mm, Ulid_ms, nu_cSt, rho_kgm3
    double precision :: L_m, nu_m2s, Re_cavity
    double precision :: xv_norm, yv_norm, psi_min, delta_BR_norm, delta_BL_norm
    double precision :: delta_bl_mm, max_vorticity, drag_force_N
    character(len=32) :: eddy_regime

    ! Read inputs
    read(*,*,iostat=iostat_val) L_cavity_mm   ! Cavity Side Length [mm] (e.g. 100.0)
    read(*,*,iostat=iostat_val) Ulid_ms       ! Top Lid Velocity [m/s] (e.g. 1.0)
    read(*,*,iostat=iostat_val) nu_cSt        ! Kinematic Viscosity [cSt] (e.g. 1.0)
    read(*,*,iostat=iostat_val) rho_kgm3      ! Fluid Density [kg/m3] (e.g. 998.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for lid-driven cavity calculation.'
        stop
    end if

    if (L_cavity_mm <= 0.0d0 .or. Ulid_ms <= 0.0d0 .or. nu_cSt <= 0.0d0) then
        write(*,*) 'ERROR: Cavity length, lid velocity, and viscosity must be positive.'
        stop
    end if

    L_m = L_cavity_mm * 1.0d-3
    nu_m2s = nu_cSt * 1.0d-6

    ! Reynolds number Re = U_lid * L / nu
    Re_cavity = (Ulid_ms * L_m) / nu_m2s

    ! Ghia et al. Benchmark empirical fits for Primary Vortex Center (xv/L, yv/L)
    if (Re_cavity <= 100.0d0) then
        xv_norm = 0.617d0
        yv_norm = 0.734d0
        psi_min = -0.103d0
        delta_BR_norm = 0.12d0
        delta_BL_norm = 0.08d0
        eddy_regime = 'Weak Corner Eddies'
    else if (Re_cavity <= 400.0d0) then
        xv_norm = 0.617d0 - 0.060d0 * (log10(Re_cavity) - 2.0d0) / 0.602d0
        yv_norm = 0.734d0 - 0.128d0 * (log10(Re_cavity) - 2.0d0) / 0.602d0
        psi_min = -0.113d0
        delta_BR_norm = 0.25d0
        delta_BL_norm = 0.16d0
        eddy_regime = 'Distinct Bottom Corner Eddies'
    else if (Re_cavity <= 1000.0d0) then
        xv_norm = 0.557d0 - 0.026d0 * (log10(Re_cavity) - 2.602d0) / 0.398d0
        yv_norm = 0.606d0 - 0.041d0 * (log10(Re_cavity) - 2.602d0) / 0.398d0
        psi_min = -0.118d0
        delta_BR_norm = 0.33d0
        delta_BL_norm = 0.22d0
        eddy_regime = 'Strong Bottom-Corner Eddies'
    else if (Re_cavity <= 5000.0d0) then
        xv_norm = 0.531d0 - 0.016d0 * (log10(Re_cavity) - 3.0d0) / 0.699d0
        yv_norm = 0.565d0 - 0.030d0 * (log10(Re_cavity) - 3.0d0) / 0.699d0
        psi_min = -0.120d0
        delta_BR_norm = 0.37d0
        delta_BL_norm = 0.27d0
        eddy_regime = 'Top-Left Tertiary Eddy Forms'
    else
        xv_norm = 0.515d0 - 0.003d0 * min(1.0d0, (log10(Re_cavity) - 3.699d0) / 0.301d0)
        yv_norm = 0.535d0 - 0.005d0 * min(1.0d0, (log10(Re_cavity) - 3.699d0) / 0.301d0)
        psi_min = -0.121d0
        delta_BR_norm = 0.39d0
        delta_BL_norm = 0.29d0
        eddy_regime = 'Complex Multi-Eddy Recirculation'
    end if

    delta_bl_mm = (L_m / sqrt(max(1.0d0, Re_cavity))) * 1000.0d0
    max_vorticity = (Ulid_ms / max(0.001d0, delta_bl_mm * 1.0d-3))
    drag_force_N = (rho_kgm3 * nu_m2s * Ulid_ms * L_m) / max(0.0001d0, delta_bl_mm * 1.0d-3)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — 2D LID-DRIVEN CAVITY BENCHMARK ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Cavity Side / Lid Velocity= ', L_cavity_mm, ' mm / ', Ulid_ms, ' m/s'
    write(*,'(A,F10.2,A,F10.2,A)') 'Kinematic Viscosity / rho = ', nu_cSt, ' cSt / ', rho_kgm3, ' kg/m3'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F12.1)')   'REYNOLDS NUMBER (Re)      = ', Re_cavity
    write(*,'(A,A)')       'Vortex Structure Regime   = ', trim(eddy_regime)
    write(*,'(A,F8.4,A,F8.4,A)') 'Primary Vortex Core (x,y) = (', xv_norm, ' L, ', yv_norm, ' L)'
    write(*,'(A,F10.4)')   'Streamfunction Min (psi)  = ', psi_min
    write(*,'(A,F8.3,A,F8.3,A)') 'Corner Eddy Lengths BR/BL = ', delta_BR_norm*L_cavity_mm, ' mm / ', delta_BL_norm*L_cavity_mm, ' mm'
    write(*,'(A,F10.3,A)') 'Wall Boundary Layer Thick = ', delta_bl_mm, ' mm'
    write(*,'(A,F10.4,A)') 'Top Lid Viscous Drag Force= ', drag_force_N, ' N per m depth'
    write(*,'(A)') '============================================================'

end program lid_driven_cavity_re