program porous_media_convection_darcy
    implicit none
    integer :: iostat_val
    double precision :: H_bed_m, W_bed_m, L_bed_m, Perm_m2, Porosity_phi, U_darcy_ms
    double precision :: Th_hot_C, Tc_cold_C, rho_dens, mu_cP, kf_fluid, ks_solid, cp_spec, beta_exp
    double precision :: mu_Pas, ke_eff, alpha_e, DeltaT, Ra_K, Pe_K, CF_forch
    double precision :: DeltaP_Pa, Nu_porous, h_porous, Q_dot_W, P_pump_W
    double precision, parameter :: G_GRAV = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) H_bed_m        ! Bed Height / Length [m] (e.g. 0.50)
    read(*,*,iostat=iostat_val) W_bed_m        ! Bed Width [m] (e.g. 0.30)
    read(*,*,iostat=iostat_val) L_bed_m        ! Bed Depth [m] (e.g. 0.30)
    read(*,*,iostat=iostat_val) Perm_m2        ! Permeability K [m2] (e.g. 1.0e-9)
    read(*,*,iostat=iostat_val) Porosity_phi   ! Porosity (0.0 to 1.0) (e.g. 0.40)
    read(*,*,iostat=iostat_val) U_darcy_ms     ! Superficial Velocity [m/s] (e.g. 0.02)
    read(*,*,iostat=iostat_val) Th_hot_C       ! Hot Temperature [deg C] (e.g. 80.0)
    read(*,*,iostat=iostat_val) Tc_cold_C      ! Cold Temperature [deg C] (e.g. 20.0)
    read(*,*,iostat=iostat_val) rho_dens       ! Fluid Density [kg/m3] (e.g. 998.0 for water)
    read(*,*,iostat=iostat_val) mu_cP          ! Fluid Viscosity [cP] (e.g. 1.0)
    read(*,*,iostat=iostat_val) kf_fluid       ! Fluid Conductivity [W/(m.K)] (e.g. 0.60)
    read(*,*,iostat=iostat_val) ks_solid       ! Solid Matrix Conductivity [W/(m.K)] (e.g. 1.5)
    read(*,*,iostat=iostat_val) cp_spec        ! Specific Heat [J/(kg.K)] (e.g. 4180.0)
    read(*,*,iostat=iostat_val) beta_exp       ! Fluid Expansion Coeff [1/K] (e.g. 0.00021)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for porous media convection calculation.'
        stop
    end if

    if (H_bed_m <= 0.0d0 .or. W_bed_m <= 0.0d0 .or. Perm_m2 <= 0.0d0 .or. Porosity_phi <= 0.0d0) then
        write(*,*) 'ERROR: Bed geometry, permeability, and porosity must be positive.'
        stop
    end if

    DeltaT = abs(Th_hot_C - Tc_cold_C)
    if (DeltaT <= 0.0d0) DeltaT = 1.0d-3

    mu_Pas = mu_cP * 1.0d-3
    ke_eff = Porosity_phi * kf_fluid + (1.0d0 - Porosity_phi) * ks_solid
    alpha_e = ke_eff / (rho_dens * cp_spec)

    ! Forchheimer Inertial Form Drag Coefficient
    CF_forch = 1.75d0 / (sqrt(150.0d0) * (Porosity_phi**1.5d0))

    ! Forchheimer-Darcy Pressure Drop across porous bed
    DeltaP_Pa = ((mu_Pas / Perm_m2) * U_darcy_ms + rho_dens * (CF_forch / sqrt(Perm_m2)) * (U_darcy_ms**2)) * H_bed_m

    ! Dimensionless Numbers
    Ra_K = (rho_dens * G_GRAV * beta_exp * DeltaT * Perm_m2 * H_bed_m) / (mu_Pas * alpha_e)
    Pe_K = (U_darcy_ms * H_bed_m) / alpha_e

    ! Mixed Porous Convection Nusselt
    if (U_darcy_ms > 1.0d-5) then
        Nu_porous = sqrt(1.0d0 + 0.318d0 * Pe_K)
    else
        Nu_porous = max(1.0d0, sqrt(max(1.0d0, Ra_K)))
    end if

    h_porous = (Nu_porous * ke_eff) / H_bed_m
    Q_dot_W = h_porous * (W_bed_m * L_bed_m) * DeltaT
    P_pump_W = (U_darcy_ms * W_bed_m * L_bed_m) * DeltaP_Pa

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — POROUS MEDIA CONVECTION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,E12.4,A,F8.3)')    'Permeability K / Porosity phi = ', Perm_m2, ' m2 / ', Porosity_phi
    write(*,'(A,F10.3,A,F10.2,A)') 'Effective Conductivity ke_eff = ', ke_eff, ' W/(m.K) (Th-Tc = ', DeltaT, ' C)'
    write(*,'(A,E12.4,A,F10.2)')   'Porous Rayleigh Ra_K / Peclet = ', Ra_K, ' / ', Pe_K
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2)')           'POROUS NUSSELT NUMBER Nu      = ', Nu_porous
    write(*,'(A,F10.2,A)')  'Heat Transfer Coeff h_eff     = ', h_porous, ' W/(m2.K)'
    write(*,'(A,F10.2,A,F10.2,A)') 'TOTAL HEAT FLOW Q_dot         = ', Q_dot_W, ' W (', Q_dot_W/1000.0d0, ' kW)'
    write(*,'(A,F10.2,A,F10.2,A)') 'Forchheimer Pressure Drop     = ', DeltaP_Pa, ' Pa (', DeltaP_Pa/1000.0d0, ' kPa)'
    write(*,'(A,F10.3,A)')  'Filtration Pumping Power      = ', P_pump_W, ' W'
    write(*,'(A)') '============================================================'

end program porous_media_convection_darcy