program pcm_latent_thermal_energy_storage
    implicit none
    integer :: iostat_val, pcm_type
    double precision :: M_pcm_kg, Tm_melt_C, T_htf_in_C, T_initial_C, m_dot_htf_kgs
    double precision :: Lf_latent_kJkg, cp_s, cp_l, k_pcm, rho_pcm, cp_htf
    double precision :: UA_overall_WK, NTU_val, epsilon_eff, Q_power_kW, E_latent_kWh, E_total_kWh
    double precision :: t_charge_hr, t_charge_min

    ! Read inputs
    read(*,*,iostat=iostat_val) pcm_type        ! 1=Paraffin RT58, 2=Inorganic Salt Hydrate, 3=Erythritol, 4=Nitrate Solar Salt
    read(*,*,iostat=iostat_val) M_pcm_kg        ! Total PCM Mass [kg] (e.g. 500.0)
    read(*,*,iostat=iostat_val) T_htf_in_C      ! HTF Charging Inlet Temp [deg C] (e.g. 75.0)
    read(*,*,iostat=iostat_val) T_initial_C     ! Initial Storage Temp [deg C] (e.g. 20.0)
    read(*,*,iostat=iostat_val) m_dot_htf_kgs   ! HTF Mass Flow Rate [kg/s] (e.g. 1.50)
    read(*,*,iostat=iostat_val) UA_overall_WK   ! Overall Heat Transfer UA [W/K] (e.g. 850.0)

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

    ! Thermophysical properties of PCMs
    if (pcm_type == 1) then ! Paraffin RT58
        Tm_melt_C = 58.0d0; Lf_latent_kJkg = 180.0d0; cp_s = 2.2d0; cp_l = 2.4d0; k_pcm = 0.22d0; rho_pcm = 880.0d0
    else if (pcm_type == 2) then ! Salt Hydrate (CaCl2.6H2O)
        Tm_melt_C = 29.5d0; Lf_latent_kJkg = 190.0d0; cp_s = 1.4d0; cp_l = 2.1d0; k_pcm = 0.54d0; rho_pcm = 1500.0d0
    else if (pcm_type == 3) then ! Erythritol (Medium Temp)
        Tm_melt_C = 118.0d0; Lf_latent_kJkg = 340.0d0; cp_s = 1.8d0; cp_l = 2.8d0; k_pcm = 0.73d0; rho_pcm = 1450.0d0
    else ! Solar Nitrate Salt (KNO3-NaNO3)
        Tm_melt_C = 220.0d0; Lf_latent_kJkg = 160.0d0; cp_s = 1.3d0; cp_l = 1.6d0; k_pcm = 0.55d0; rho_pcm = 1900.0d0
    end if

    cp_htf = 4.184d0 ! Water HTF [kJ/(kg.K)]
    if (T_htf_in_C > 95.0d0) cp_htf = 2.10d0 ! Thermal oil

    ! NTU & Effectiveness of Latent Storage Module
    NTU_val = UA_overall_WK / (m_dot_htf_kgs * cp_htf * 1000.0d0)
    epsilon_eff = 1.0d0 - exp(-NTU_val)

    ! Storage Energy (kWh)
    E_latent_kWh = (M_pcm_kg * Lf_latent_kJkg) / 3600.0d0
    E_total_kWh = (M_pcm_kg * (cp_s * max(0.0d0, Tm_melt_C - T_initial_C) + Lf_latent_kJkg + &
                  cp_l * max(0.0d0, T_htf_in_C - Tm_melt_C))) / 3600.0d0

    ! Average Charging Thermal Power [kW]
    Q_power_kW = m_dot_htf_kgs * cp_htf * (T_htf_in_C - Tm_melt_C) * epsilon_eff
    if (Q_power_kW < 0.1d0) Q_power_kW = 0.1d0

    ! Total Charging Time
    t_charge_hr = E_total_kWh / Q_power_kW
    t_charge_min = t_charge_hr * 60.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — PCM LATENT HEAT THERMAL ENERGY STORAGE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'PCM Mass / Melting Temp Tm = ', M_pcm_kg, ' kg / ', Tm_melt_C, ' deg C'
    write(*,'(A,F10.1,A,F10.2,A)') 'HTF Inlet Temp / Mass Flow = ', T_htf_in_C, ' deg C / ', m_dot_htf_kgs, ' kg/s'
    write(*,'(A,F10.3,A,F10.3)')   'Storage NTU / Effectiveness= ', NTU_val, ' / ', epsilon_eff
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'TOTAL STORED ENERGY       = ', E_total_kWh, ' kWh'
    write(*,'(A,F10.2,A)')  'Latent Portion (Lf)       = ', E_latent_kWh, ' kWh'
    write(*,'(A,F10.2,A)')  'AVERAGE CHARGING POWER    = ', Q_power_kW, ' kW'
    write(*,'(A,F10.2,A,F10.1,A)') 'FULL CHARGING DURATION    = ', t_charge_hr, ' hours (', t_charge_min, ' min)'
    write(*,'(A)') '============================================================'

end program pcm_latent_thermal_energy_storage