program steam_calc
    implicit none

    integer, parameter :: rk8 = kind(1.0D+00)
    
    ! Input variables
    integer :: input_type
    double precision :: val1, val2
    
    ! Conversions
    double precision :: t_c, t_k, p_kpa, p_mpa
    double precision :: x_qual, spec_vol, spec_u, spec_h, spec_s, density
    double precision :: cp_val, cv_val, speed_sound, visc_eta, therm_lambda
    
    ! Saturated values
    double precision :: t_sat_k, p_sat_mpa, rhol, rhov
    double precision :: h_liq, h_vap, s_liq, s_vap, u_liq, u_vap
    double precision :: cp_liq, cp_vap, cv_liq, cv_vap, sound_liq, sound_vap
    double precision :: visc_liq, visc_vap, therm_liq, therm_vap
    
    ! Calculation values
    character(len=30) :: state_name
    double precision :: rho, rho_start, dpdr, dpdt
    double precision :: ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub
    double precision :: ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur
    double precision :: a_term, cjth, cjtt, cp_term, cv_term, g_term, h_term, p_term, s_term, u_term
    
    ! Solver variables
    double precision :: t_min, t_max, t_mid, h_mid, s_mid, err
    integer :: iter
    logical :: is_mix
    
    ! Functions from steam.f90
    double precision :: gascon
    
    ! Read from stdin
    read(*,*) input_type
    read(*,*) val1
    read(*,*) val2
    
    ! Default quality set to -1 (meaning not in mixture region, single phase)
    x_qual = -1.0d0
    is_mix = .false.
    
    select case(input_type)
    
    ! -------------------------------------------------------------
    ! 1. Temperature and Pressure (val1 = T [C], val2 = P [kPa])
    ! -------------------------------------------------------------
    case(1)
        t_c = val1
        t_k = t_c + 273.15d0
        p_kpa = val2
        p_mpa = p_kpa / 1000.0d0
        
        ! Check boundaries
        if (t_k < 273.16d0 .or. t_k > 1273.15d0) then
            write(*,*) 'ERROR: Temperature must be between 0.01 C and 1000 C'
            stop
        end if
        if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
            write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
            stop
        end if
        
        if (p_mpa < 22.055d0) then
            call tsat(p_mpa, t_sat_k, rhol, rhov)
            if (abs(t_k - t_sat_k) < 1.0d-3) then
                ! On the saturation boundary, assume superheated or mixture
                ! But T & P is not enough to define mixture state, so default to saturated liquid
                state_name = 'Saturated boundary'
                rho_start = rhol
            else if (t_k < t_sat_k) then
                state_name = 'Subcooled Liquid'
                rho_start = 1.11d0 - 0.0004d0 * t_k
            else
                state_name = 'Superheated Vapor'
                rho_start = p_mpa / (gascon() * t_k)
            end if
        else
            state_name = 'Supercritical Fluid'
            if (t_k < 647.13d0) then
                rho_start = 1.11d0 - 0.0004d0 * t_k
            else
                rho_start = p_mpa / (gascon() * t_k)
            end if
        end if
        
        call dense(p_mpa, t_k, rho_start, rho, dpdr)
        call therm(t_k, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
                   g_term, spec_h, p_term, spec_s, spec_u)
        call viscosity(t_k, rho, visc_eta)
        call thercon(t_k, rho, therm_lambda)
        call sound(t_k, p_mpa, speed_sound)
        
    ! -------------------------------------------------------------
    ! 2. Temperature and Quality (val1 = T [C], val2 = x)
    ! -------------------------------------------------------------
    case(2)
        t_c = val1
        t_k = t_c + 273.15d0
        x_qual = val2
        
        if (t_k < 273.16d0 .or. t_k >= 647.13d0) then
            write(*,*) 'ERROR: Saturated state requires T between 0.01 C and 373.9 C'
            stop
        end if
        if (x_qual < 0.0d0 .or. x_qual > 1.0d0) then
            write(*,*) 'ERROR: Quality x must be between 0 and 1'
            stop
        end if
        
        call psat(t_k, p_sat_mpa, rhol, rhov)
        p_mpa = p_sat_mpa
        p_kpa = p_mpa * 1000.0d0
        is_mix = .true.
        
        if (x_qual == 0.0d0) then
            state_name = 'Saturated Liquid'
        else if (x_qual == 1.0d0) then
            state_name = 'Saturated Vapor'
        else
            state_name = 'Saturated Mixture'
        end if
        
        ! Calculate properties of saturated liquid
        call therm(t_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
                   g_term, h_liq, p_term, s_liq, u_liq)
        call viscosity(t_k, rhol, visc_liq)
        call thercon(t_k, rhol, therm_liq)
        call sound(t_k, p_mpa, sound_liq)
        
        ! Calculate properties of saturated vapor
        call therm(t_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
                   g_term, h_vap, p_term, s_vap, u_vap)
        call viscosity(t_k, rhov, visc_vap)
        call thercon(t_k, rhov, therm_vap)
        call sound(t_k, p_mpa, sound_vap)
        
        ! Blend properties
        spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
        spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
        spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
        spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
        density = 1.0d0 / (spec_vol * 1000.0d0)
        cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
        cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
        speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
        visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
        therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
        rho = density
        
    ! -------------------------------------------------------------
    ! 3. Pressure and Quality (val1 = P [kPa], val2 = x)
    ! -------------------------------------------------------------
    case(3)
        p_kpa = val1
        p_mpa = p_kpa / 1000.0d0
        x_qual = val2
        
        if (p_mpa < 0.000611d0 .or. p_mpa >= 22.055d0) then
            write(*,*) 'ERROR: Saturated state requires P between 0.61 kPa and 22055 kPa'
            stop
        end if
        if (x_qual < 0.0d0 .or. x_qual > 1.0d0) then
            write(*,*) 'ERROR: Quality x must be between 0 and 1'
            stop
        end if
        
        call tsat(p_mpa, t_sat_k, rhol, rhov)
        t_k = t_sat_k
        t_c = t_k - 273.15d0
        is_mix = .true.
        
        if (x_qual == 0.0d0) then
            state_name = 'Saturated Liquid'
        else if (x_qual == 1.0d0) then
            state_name = 'Saturated Vapor'
        else
            state_name = 'Saturated Mixture'
        end if
        
        ! Calculate properties of saturated liquid
        call therm(t_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
                   g_term, h_liq, p_term, s_liq, u_liq)
        call viscosity(t_k, rhol, visc_liq)
        call thercon(t_k, rhol, therm_liq)
        call sound(t_k, p_mpa, sound_liq)
        
        ! Calculate properties of saturated vapor
        call therm(t_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
                   g_term, h_vap, p_term, s_vap, u_vap)
        call viscosity(t_k, rhov, visc_vap)
        call thercon(t_k, rhov, therm_vap)
        call sound(t_k, p_mpa, sound_vap)
        
        ! Blend properties
        spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
        spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
        spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
        spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
        density = 1.0d0 / (spec_vol * 1000.0d0)
        cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
        cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
        speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
        visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
        therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
        rho = density
        
    ! -------------------------------------------------------------
    ! 4. Pressure and Enthalpy (val1 = P [kPa], val2 = h [kJ/kg])
    ! -------------------------------------------------------------
    case(4)
        p_kpa = val1
        p_mpa = p_kpa / 1000.0d0
        spec_h = val2
        
        if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
            write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
            stop
        end if
        
        if (p_mpa < 22.055d0) then
            ! Check if it falls in saturation boundary
            call tsat(p_mpa, t_sat_k, rhol, rhov)
            call therm(t_sat_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
                       g_term, h_liq, p_term, s_liq, u_liq)
            call therm(t_sat_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
                       g_term, h_vap, p_term, s_vap, u_vap)
            
            if (spec_h >= h_liq .and. spec_h <= h_vap) then
                ! Saturated mixture!
                x_qual = (spec_h - h_liq) / (h_vap - h_liq)
                t_k = t_sat_k
                t_c = t_k - 273.15d0
                is_mix = .true.
                state_name = 'Saturated Mixture'
                
                call viscosity(t_k, rhol, visc_liq)
                call thercon(t_k, rhol, therm_liq)
                call sound(t_k, p_mpa, sound_liq)
                
                call viscosity(t_k, rhov, visc_vap)
                call thercon(t_k, rhov, therm_vap)
                call sound(t_k, p_mpa, sound_vap)
                
                spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
                spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
                spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
                density = 1.0d0 / (spec_vol * 1000.0d0)
                cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
                cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
                speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
                visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
                therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
                rho = density
            else if (spec_h < h_liq) then
                ! Subcooled liquid
                state_name = 'Subcooled Liquid'
                t_min = 273.16d0
                t_max = t_sat_k
                rho_start = 1.0d0
                
                ! Bisection solver for T
                do iter = 1, 40
                    t_mid = 0.5d0 * (t_min + t_max)
                    call dense(p_mpa, t_mid, rho_start, rho, dpdr)
                    call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
                               g_term, h_mid, p_term, s_mid, u_term)
                    if (h_mid < spec_h) then
                        t_min = t_mid
                    else
                        t_max = t_mid
                    end if
                end do
                t_k = 0.5d0 * (t_min + t_max)
                t_c = t_k - 273.15d0
                spec_h = h_mid
                spec_s = s_mid
                spec_u = u_term
                call viscosity(t_k, rho, visc_eta)
                call thercon(t_k, rho, therm_lambda)
                call sound(t_k, p_mpa, speed_sound)
            else
                ! Superheated vapor
                state_name = 'Superheated Vapor'
                t_min = t_sat_k
                t_max = 1273.15d0
                
                ! Bisection solver for T
                do iter = 1, 40
                    t_mid = 0.5d0 * (t_min + t_max)
                    rho_start = p_mpa / (gascon() * t_mid)
                    call dense(p_mpa, t_mid, rho_start, rho, dpdr)
                    call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
                               g_term, h_mid, p_term, s_mid, u_term)
                    if (h_mid < spec_h) then
                        t_min = t_mid
                    else
                        t_max = t_mid
                    end if
                end do
                t_k = 0.5d0 * (t_min + t_max)
                t_c = t_k - 273.15d0
                spec_h = h_mid
                spec_s = s_mid
                spec_u = u_term
                call viscosity(t_k, rho, visc_eta)
                call thercon(t_k, rho, therm_lambda)
                call sound(t_k, p_mpa, speed_sound)
            end if
        else
            ! Supercritical region
            state_name = 'Supercritical Fluid'
            t_min = 273.16d0
            t_max = 1273.15d0
            
            do iter = 1, 40
                t_mid = 0.5d0 * (t_min + t_max)
                if (t_mid < 647.13d0) then
                    rho_start = 1.0d0
                else
                    rho_start = p_mpa / (gascon() * t_mid)
                end if
                call dense(p_mpa, t_mid, rho_start, rho, dpdr)
                call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
                           g_term, h_mid, p_term, s_mid, u_term)
                if (h_mid < spec_h) then
                    t_min = t_mid
                else
                    t_max = t_mid
                end if
            end do
            t_k = 0.5d0 * (t_min + t_max)
            t_c = t_k - 273.15d0
            spec_h = h_mid
            spec_s = s_mid
            spec_u = u_term
            call viscosity(t_k, rho, visc_eta)
            call thercon(t_k, rho, therm_lambda)
            call sound(t_k, p_mpa, speed_sound)
        end if
        
    ! -------------------------------------------------------------
    ! 5. Pressure and Entropy (val1 = P [kPa], val2 = s [kJ/(kg K)])
    ! -------------------------------------------------------------
    case(5)
        p_kpa = val1
        p_mpa = p_kpa / 1000.0d0
        spec_s = val2
        
        if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
            write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
            stop
        end if
        
        if (p_mpa < 22.055d0) then
            ! Check if it falls in saturation boundary
            call tsat(p_mpa, t_sat_k, rhol, rhov)
            call therm(t_sat_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
                       g_term, h_liq, p_term, s_liq, u_liq)
            call therm(t_sat_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
                       g_term, h_vap, p_term, s_vap, u_vap)
            
            if (spec_s >= s_liq .and. spec_s <= s_vap) then
                ! Saturated mixture!
                x_qual = (spec_s - s_liq) / (s_vap - s_liq)
                t_k = t_sat_k
                t_c = t_k - 273.15d0
                is_mix = .true.
                state_name = 'Saturated Mixture'
                
                call viscosity(t_k, rhol, visc_liq)
                call thercon(t_k, rhol, therm_liq)
                call sound(t_k, p_mpa, sound_liq)
                
                call viscosity(t_k, rhov, visc_vap)
                call thercon(t_k, rhov, therm_vap)
                call sound(t_k, p_mpa, sound_vap)
                
                spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
                spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
                spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
                density = 1.0d0 / (spec_vol * 1000.0d0)
                cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
                cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
                speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
                visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
                therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
                rho = density
            else if (spec_s < s_liq) then
                ! Subcooled liquid
                state_name = 'Subcooled Liquid'
                t_min = 273.16d0
                t_max = t_sat_k
                rho_start = 1.0d0
                
                ! Bisection solver for T
                do iter = 1, 40
                    t_mid = 0.5d0 * (t_min + t_max)
                    call dense(p_mpa, t_mid, rho_start, rho, dpdr)
                    call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
                               g_term, h_mid, p_term, s_mid, u_term)
                    if (s_mid < spec_s) then
                        t_min = t_mid
                    else
                        t_max = t_mid
                    end if
                end do
                t_k = 0.5d0 * (t_min + t_max)
                t_c = t_k - 273.15d0
                spec_h = h_mid
                spec_s = s_mid
                spec_u = u_term
                call viscosity(t_k, rho, visc_eta)
                call thercon(t_k, rho, therm_lambda)
                call sound(t_k, p_mpa, speed_sound)
            else
                ! Superheated vapor
                state_name = 'Superheated Vapor'
                t_min = t_sat_k
                t_max = 1273.15d0
                
                ! Bisection solver for T
                do iter = 1, 40
                    t_mid = 0.5d0 * (t_min + t_max)
                    rho_start = p_mpa / (gascon() * t_mid)
                    call dense(p_mpa, t_mid, rho_start, rho, dpdr)
                    call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
                               g_term, h_mid, p_term, s_mid, u_term)
                    if (s_mid < spec_s) then
                        t_min = t_mid
                    else
                        t_max = t_mid
                    end if
                end do
                t_k = 0.5d0 * (t_min + t_max)
                t_c = t_k - 273.15d0
                spec_h = h_mid
                spec_s = s_mid
                spec_u = u_term
                call viscosity(t_k, rho, visc_eta)
                call thercon(t_k, rho, therm_lambda)
                call sound(t_k, p_mpa, speed_sound)
            end if
        else
            ! Supercritical region
            state_name = 'Supercritical Fluid'
            t_min = 273.16d0
            t_max = 1273.15d0
            
            do iter = 1, 40
                t_mid = 0.5d0 * (t_min + t_max)
                if (t_mid < 647.13d0) then
                    rho_start = 1.0d0
                else
                    rho_start = p_mpa / (gascon() * t_mid)
                end if
                call dense(p_mpa, t_mid, rho_start, rho, dpdr)
                call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
                           g_term, h_mid, p_term, s_mid, u_term)
                if (s_mid < spec_s) then
                    t_min = t_mid
                else
                    t_max = t_mid
                end if
            end do
            t_k = 0.5d0 * (t_min + t_max)
            t_c = t_k - 273.15d0
            spec_h = h_mid
            spec_s = s_mid
            spec_u = u_term
            call viscosity(t_k, rho, visc_eta)
            call thercon(t_k, rho, therm_lambda)
            call sound(t_k, p_mpa, speed_sound)
        end if
    
    end select
    
    ! Density conversion: rho is in g/cm3. Output specific volume v in m3/kg.
    if (.not. is_mix) then
        density = rho * 1000.0d0     ! kg/m3
        spec_vol = 1.0d0 / density   ! m3/kg
    end if
    
    ! If it's single phase, compute dynamic viscosity in Pa s
    ! visc_eta is in microPascal-seconds. Let's convert to Pa s: * 1.0E-6
    visc_eta = visc_eta * 1.0d-6
    
    ! Convert thermal conductivity to W/(m K): therm_lambda is in mW/(m K) => / 1000.d0
    therm_lambda = therm_lambda / 1000.d0
    
    ! Convert quality
    if (.not. is_mix) then
        if (state_name == 'Subcooled Liquid') then
            x_qual = 0.0d0
        else if (state_name == 'Superheated Vapor') then
            x_qual = 1.0d0
        else
            x_qual = -1.0d0  ! Supercritical has no quality defined
        end if
    end if
    
    ! Print output in parsable format
    write(*,'(A,A)') 'STATE=', trim(state_name)
    write(*,'(A,F14.4)') 'T=', t_c
    write(*,'(A,F14.4)') 'P=', p_kpa
    write(*,'(A,F14.6)') 'V=', spec_vol
    write(*,'(A,F14.4)') 'H=', spec_h
    write(*,'(A,F14.6)') 'S=', spec_s
    write(*,'(A,F14.4)') 'U=', spec_u
    write(*,'(A,F14.4)') 'CP=', cp_val
    write(*,'(A,F14.4)') 'CV=', cv_val
    write(*,'(A,F14.2)') 'SOUND=', speed_sound
    write(*,'(A,E14.6)') 'VISC=', visc_eta
    write(*,'(A,F14.6)') 'THERM=', therm_lambda
    write(*,'(A,F14.4)') 'QUALITY=', x_qual

end program steam_calc