program atmospheric_boundary_layer_wind
    implicit none
    integer :: iostat_val
    double precision :: Uref_ms, zref_m, z0_rough_m, alpha_exp, Hbld_m, bld_width_m, rho_air
    double precision :: u_star, U_top_log, U_top_pow, q_top_Pa, q_top_kPa, Iz_top_pct
    double precision :: total_shear_kN, overturn_moment_kNm
    double precision, parameter :: KAPPA = 0.40d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Uref_ms        ! Reference Wind Speed at 10m [m/s] (e.g. 25.0)
    read(*,*,iostat=iostat_val) zref_m         ! Reference Height [m] (e.g. 10.0)
    read(*,*,iostat=iostat_val) z0_rough_m     ! Terrain Roughness Length z0 [m] (e.g. 0.30)
    read(*,*,iostat=iostat_val) alpha_exp      ! Power-Law Exponent alpha (e.g. 0.22)
    read(*,*,iostat=iostat_val) Hbld_m         ! Building Height [m] (e.g. 150.0)
    read(*,*,iostat=iostat_val) bld_width_m    ! Building Facade Width [m] (e.g. 35.0)
    read(*,*,iostat=iostat_val) rho_air        ! Air Density [kg/m3] (e.g. 1.225)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for atmospheric boundary layer calculation.'
        stop
    end if

    if (Uref_ms <= 0.0d0 .or. z0_rough_m <= 0.0d0 .or. Hbld_m <= 0.0d0) then
        write(*,*) 'ERROR: Wind speed, roughness, and building height must be positive.'
        stop
    end if

    ! Friction velocity u* [m/s]
    u_star = (KAPPA * Uref_ms) / log(max(1.1d0, zref_m / z0_rough_m))

    ! Wind Velocity at Top of Building
    U_top_log = (u_star / KAPPA) * log(max(1.05d0, Hbld_m / z0_rough_m))
    U_top_pow = Uref_ms * ((Hbld_m / zref_m)**alpha_exp)

    ! Dynamic Stagnation Pressure at Top [Pa]
    q_top_Pa = 0.5d0 * rho_air * (U_top_log**2)
    q_top_kPa = q_top_Pa / 1000.0d0

    ! Turbulence Intensity at Building Top [%]
    Iz_top_pct = (1.0d0 / log(max(1.05d0, Hbld_m / z0_rough_m))) * 100.0d0

    ! Total Integrated Base Shear Force [kN] and Overturning Moment [kN.m] (approximated)
    total_shear_kN = (0.5d0 * rho_air * (Uref_ms**2) * bld_width_m * Hbld_m * (1.0d0 / (2.0d0 * alpha_exp + 1.0d0)) * &
                     ((Hbld_m / zref_m)**(2.0d0 * alpha_exp))) / 1000.0d0
    overturn_moment_kNm = total_shear_kN * (Hbld_m * ((2.0d0 * alpha_exp + 1.0d0) / (2.0d0 * alpha_exp + 2.0d0)))

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — ATMOSPHERIC BOUNDARY LAYER WIND ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.1,A)') 'Reference Wind Uref / zref = ', Uref_ms, ' m/s / ', zref_m, ' m'
    write(*,'(A,F10.4,A,F10.3)')   'Roughness z0 / Power exp a= ', z0_rough_m, ' m / ', alpha_exp
    write(*,'(A,F10.1,A,F10.1,A)') 'Building Height / Width   = ', Hbld_m, ' m / ', bld_width_m, ' m'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Friction Velocity (u*)    = ', u_star, ' m/s'
    write(*,'(A,F10.2,A)')  'TOP WIND SPEED (Log-Law)  = ', U_top_log, ' m/s'
    write(*,'(A,F10.2,A)')  'TOP WIND SPEED (Power-Law)= ', U_top_pow, ' m/s'
    write(*,'(A,F10.2,A,F8.3,A)') 'Top Dynamic Wind Pressure = ', q_top_Pa, ' Pa (', q_top_kPa, ' kPa)'
    write(*,'(A,F10.2,A)')  'Top Turbulence Intensity  = ', Iz_top_pct, ' %'
    write(*,'(A,F10.2,A)')  'Base Overturning Moment   = ', overturn_moment_kNm, ' kN.m'
    write(*,'(A,F10.2,A)')  'Total Wind Base Shear     = ', total_shear_kN, ' kN'
    write(*,'(A)') '============================================================'

end program atmospheric_boundary_layer_wind
