program impinging_jet_cooling_array
    implicit none
    integer :: iostat_val
    double precision :: d_nozzle_mm, S_pitch_mm, H_dist_mm, U_jet_ms, Ts_target_C, Tj_fluid_C
    double precision :: rho_dens, mu_cP, k_cond, cp_spec
    double precision :: d_m, S_m, H_m, mu_Pas, Re_d, Pr_num, f_area, H_d_ratio
    double precision :: Nu_avg, Nu_peak, h_avg, h_peak, q_flux_kWm2, DeltaP_jet_kPa
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) d_nozzle_mm   ! Orifice Diameter [mm] (e.g. 3.0)
    read(*,*,iostat=iostat_val) S_pitch_mm    ! Nozzle Pitch [mm] (e.g. 15.0)
    read(*,*,iostat=iostat_val) H_dist_mm     ! Nozzle-to-Plate Distance [mm] (e.g. 12.0)
    read(*,*,iostat=iostat_val) U_jet_ms      ! Jet Discharge Velocity [m/s] (e.g. 25.0)
    read(*,*,iostat=iostat_val) Ts_target_C   ! Target Plate Temperature [deg C] (e.g. 85.0)
    read(*,*,iostat=iostat_val) Tj_fluid_C    ! Coolant Jet Temperature [deg C] (e.g. 25.0)
    read(*,*,iostat=iostat_val) rho_dens      ! Fluid Density [kg/m3] (e.g. 1.184 for air or 997 for water)
    read(*,*,iostat=iostat_val) mu_cP         ! Fluid Viscosity [cP] (e.g. 0.0185 for air or 0.89 for water)
    read(*,*,iostat=iostat_val) k_cond        ! Thermal Conductivity [W/(m.K)] (e.g. 0.026 for air or 0.60 for water)
    read(*,*,iostat=iostat_val) cp_spec       ! Specific Heat [J/(kg.K)] (e.g. 1005.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for impinging jet array calculation.'
        stop
    end if

    if (d_nozzle_mm <= 0.0d0 .or. S_pitch_mm <= d_nozzle_mm .or. H_dist_mm <= 0.0d0 .or. U_jet_ms <= 0.0d0) then
        write(*,*) 'ERROR: Nozzle geometry (S > d > 0, H > 0) and velocity must be positive.'
        stop
    end if

    d_m = d_nozzle_mm * 1.0d-3
    S_m = S_pitch_mm * 1.0d-3
    H_m = H_dist_mm * 1.0d-3
    mu_Pas = mu_cP * 1.0d-3

    ! Dimensionless Parameters
    Re_d = (rho_dens * U_jet_ms * d_m) / mu_Pas
    Pr_num = (cp_spec * mu_Pas) / k_cond
    f_area = (PI / 4.0d0) / ((S_pitch_mm / d_nozzle_mm)**2)
    H_d_ratio = H_dist_mm / d_nozzle_mm

    ! Martin (1977) Correlation for Array of Round Impinging Jets
    Nu_avg = 2.0d0 * sqrt(f_area) * ((1.0d0 - 2.2d0 * sqrt(f_area)) / &
             (1.0d0 + 0.2d0 * (H_d_ratio - 6.0d0) * sqrt(f_area))) * &
             (2.0d0 * sqrt(Re_d) * sqrt(1.0d0 + 0.005d0 * (Re_d**0.55d0))) * (Pr_num**0.42d0)

    Nu_peak = 0.8d0 * sqrt(Re_d) * (Pr_num**0.40d0) * (max(1.0d0, H_d_ratio)**(-0.05d0))

    h_avg = (Nu_avg * k_cond) / d_m
    h_peak = (Nu_peak * k_cond) / d_m
    q_flux_kWm2 = (h_avg * (Ts_target_C - Tj_fluid_C)) / 1000.0d0

    ! Orifice Pressure Drop DeltaP [kPa]
    DeltaP_jet_kPa = (0.5d0 * rho_dens * (U_jet_ms**2) * 1.5d0) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — IMPINGING JET ARRAY HEAT TRANSFER ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Orifice d / Pitch S (H/d) = ', d_nozzle_mm, ' mm / ', S_pitch_mm, ' mm (H/d = ', H_d_ratio, ')'
    write(*,'(A,F10.2,A,F10.2,A)') 'Jet Velocity / Reynolds Re= ', U_jet_ms, ' m/s / ', Re_d, ' (f_area = ', f_area*100.0d0, ' %)'
    write(*,'(A,F10.1,A,F10.1,A)') 'Plate Temp / Jet Temp     = ', Ts_target_C, ' deg C / ', Tj_fluid_C, ' deg C'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2)')   'AVERAGE NUSSELT (Nu_avg)  = ', Nu_avg, ' (Peak Nu_0 = ', Nu_peak, ')'
    write(*,'(A,F10.1,A,F10.1,A)') 'Heat Transfer Coeff h_avg = ', h_avg, ' W/m2.K (Peak = ', h_peak, ' W/m2.K)'
    write(*,'(A,F10.2,A)')  'AVERAGE HEAT FLUX q"      = ', q_flux_kWm2, ' kW/m2'
    write(*,'(A,F10.2,A)')  'Orifice Jet Pressure Drop = ', DeltaP_jet_kPa, ' kPa'
    write(*,'(A)') '============================================================'

end program impinging_jet_cooling_array