! ============================================================================
! ThermoFluidCalc — IAPWS-IF97 Water and Steam Properties Calculator
! Self-contained Fortran 90 implementation generated programmatically
! ============================================================================
module iapws_if97_mod
    implicit none

    double precision, parameter :: R = 0.461526D0 ! kJ/(kg K)
    double precision, parameter :: Tc = 647.096D0 ! K
    double precision, parameter :: Pc = 22.064D0  ! MPa
    double precision, parameter :: rhoc = 322.0D0 ! kg/m3
    double precision, parameter :: Pmin = 0.000611212677444D0 ! MPa
    double precision, parameter :: Ps_623 = 16.5291642526D0 ! MPa

    ! Region 1 coefficients
    double precision, dimension(34) :: r1_n = (/ &
          1.463297121316700D-01,  -8.454818716911400D-01,  -3.756360367204000D+00, &
          3.385516916838500D+00,  -9.579196338787200D-01,   1.577203851322800D-01, &
         -1.661641719950100D-02,   8.121462998356800D-04,   2.831908012380400D-04, &
         -6.070630156587400D-04,  -1.899006821841900D-02,  -3.252974877050500D-02, &
         -2.184171717541400D-02,  -5.283835796993000D-05,  -4.718432107326700D-04, &
         -3.000178079302600D-04,   4.766139390698700D-05,  -4.414184533084600D-06, &
         -7.269499629759400D-16,  -3.167964484505400D-05,  -2.827079798531200D-06, &
         -8.520512812010300D-10,  -2.242528190800000D-06,  -6.517122289560100D-07, &
         -1.434172993792400D-13,  -4.051699686011700D-07,  -1.273430174164100D-09, &
         -1.742487123063400D-10,  -6.876213129553100D-19,   1.447830782852100D-20, &
          2.633578166279500D-23,  -1.194762264007100D-23,   1.822809458140400D-24, &
         -9.353708729245800D-26 /)
    integer, dimension(34) :: r1_li = (/ &
        0, 0, 0, 0, 0, 0, 0, 0, &
        1, 1, 1, 1, 1, 1, 2, 2, &
        2, 2, 2, 3, 3, 3, 4, 4, &
        4, 5, 8, 8, 21, 23, 29, 30, &
        31, 32 /)
    integer, dimension(34) :: r1_lj = (/ &
        -2, -1, 0, 1, 2, 3, 4, 5, &
        -9, -7, -1, 0, 1, 3, -3, 0, &
        1, 3, 17, -4, 0, 6, -5, -2, &
        10, -8, -11, -6, -29, -31, -38, -39, &
        -40, -41 /)

    ! Region 2 coefficients
    double precision, dimension(43) :: r2_n = (/ &
         -1.773174247321300D-03,  -1.783486229235800D-02,  -4.599601369636500D-02, &
         -5.758125908343200D-02,  -5.032527872793000D-02,  -3.303264167020300D-05, &
         -1.894898751631500D-04,  -3.939277724335500D-03,  -4.379729565057300D-02, &
         -2.667454791408700D-05,   2.048173769230900D-08,   4.387066728443500D-07, &
         -3.227767723857000D-05,  -1.503392454214800D-03,  -4.066825356264900D-02, &
         -7.884730955936700D-10,   1.279071785228500D-08,   4.822537271850700D-07, &
          2.292207633766100D-06,  -1.671476645106100D-11,  -2.117147232135500D-03, &
         -2.389574193410400D+01,  -5.905956432427000D-18,  -1.262180889910100D-06, &
         -3.894684243573900D-02,   1.125621136045900D-11,  -8.231134089799800D+00, &
          1.980971280208800D-08,   1.040696521017400D-19,  -1.023474709592900D-13, &
         -1.001817937951100D-09,  -8.088290864698500D-11,   1.069303187940900D-01, &
         -3.366225057417100D-01,   8.918584535542100D-25,   3.062931687623200D-13, &
         -4.200246769820800D-06,  -5.905602968563900D-26,   3.782694761345700D-06, &
         -1.276860893468100D-15,   7.308761059506100D-29,   5.541471535077800D-17, &
         -9.436970724121000D-07 /)
    integer, dimension(43) :: r2_li = (/ &
        1, 1, 1, 1, 1, 2, 2, 2, &
        2, 2, 3, 3, 3, 3, 3, 4, &
        4, 4, 5, 6, 6, 6, 7, 7, &
        7, 8, 8, 9, 10, 10, 10, 16, &
        16, 18, 20, 20, 20, 21, 22, 23, &
        24, 24, 24 /)
    integer, dimension(43) :: r2_lj = (/ &
        0, 1, 2, 3, 6, 1, 2, 4, &
        7, 36, 0, 1, 3, 6, 35, 1, &
        2, 3, 7, 3, 16, 35, 0, 11, &
        25, 8, 36, 13, 4, 10, 14, 29, &
        50, 57, 20, 35, 48, 21, 53, 39, &
        26, 40, 58 /)

    ! Region 2 Ideal gas cp0 coefficients
    double precision, dimension(9) :: r2_cp0_no = (/ &
         -9.692768650021700D+00,   1.008665596801800D+01,  -5.608791128302000D-03, &
          7.145273808145500D-02,  -4.071049822392800D-01,   1.424081917144400D+00, &
         -4.383951131945000D+00,  -2.840863246077200D-01,   2.126846375330700D-02 /)
    integer, dimension(9) :: r2_cp0_jo = (/ &
        0, 1, -5, -4, -3, -2, -1, 2, &
        3 /)

    ! Region 3 coefficients
    double precision, dimension(39) :: r3_n = (/ &
         -1.573284529023900D+01,   2.094439697430700D+01,  -7.686770787871600D+00, &
          2.618594778795400D+00,  -2.808078114862000D+00,   1.205336969651700D+00, &
         -8.456681281250201D-03,  -1.265431547771400D+00,  -1.152440780668100D+00, &
          8.852104398431800D-01,  -6.420776518160700D-01,   3.849346018667100D-01, &
         -8.521470882420600D-01,   4.897228154187700D+00,  -3.050261725696500D+00, &
          3.942053687915400D-02,   1.255840842430800D-01,  -2.799932969871000D-01, &
          1.389979956946000D+00,  -2.018991502357000D+00,  -8.214763717396300D-03, &
         -4.759603573492300D-01,   4.398407447350000D-02,  -4.447643542873900D-01, &
          9.057207071973300D-01,   7.052245008796700D-01,   1.077051262633200D-01, &
         -3.291362325895400D-01,  -5.087106204115800D-01,  -2.217540087309600D-02, &
          9.426075166509200D-02,   1.643627844796100D-01,  -1.350337224134800D-02, &
         -1.483434535247200D-02,   5.792295362808399D-04,   3.230890470371100D-03, &
          8.096480299621501D-05,  -1.655767979503700D-04,  -4.492389906181500D-05 /)
    integer, dimension(39) :: r3_li = (/ &
        0, 0, 0, 0, 0, 0, 0, 1, &
        1, 1, 1, 2, 2, 2, 2, 2, &
        2, 3, 3, 3, 3, 3, 4, 4, &
        4, 4, 5, 5, 5, 6, 6, 6, &
        7, 8, 9, 9, 10, 10, 11 /)
    integer, dimension(39) :: r3_lj = (/ &
        0, 1, 2, 7, 10, 12, 23, 2, &
        6, 15, 17, 0, 2, 6, 7, 22, &
        26, 0, 2, 4, 16, 26, 0, 2, &
        4, 26, 1, 3, 26, 0, 2, 26, &
        2, 26, 2, 26, 0, 1, 26 /)


contains

    ! -------------------------------------------------------------------------
    ! Saturation Pressure: P = f(T) in MPa
    ! -------------------------------------------------------------------------
    double precision function psat(T)
        double precision, intent(in) :: T
        double precision :: n(0:10)
        double precision :: tita, A, B, C

        n(0) = 0.0D0
        n(1) = 0.11670521452767D+04
        n(2) = -0.72421316703206D+06
        n(3) = -0.17073846940092D+02
        n(4) = 0.12020824702470D+05
        n(5) = -0.32325550322333D+07
        n(6) = 0.14915108613530D+02
        n(7) = -0.48232657361591D+04
        n(8) = 0.40511340542057D+06
        n(9) = -0.23855557567849D+00
        n(10) = 0.65017534844798D+03

        if (T < 273.15D0 .or. T > Tc) then
            psat = -1.0D0
            return
        end if

        tita = T + n(9) / (T - n(10))
        A = tita**2 + n(1)*tita + n(2)
        B = n(3)*tita**2 + n(4)*tita + n(5)
        C = n(6)*tita**2 + n(7)*tita + n(8)
        psat = (2.0D0 * C / (-B + sqrt(B**2 - 4.0D0*A*C)))**4
    end function psat

    ! -------------------------------------------------------------------------
    ! Saturation Temperature: T = f(P) in K
    ! -------------------------------------------------------------------------
    double precision function tsat(P)
        double precision, intent(in) :: P
        double precision :: n(0:10)
        double precision :: beta, E, F, G, D

        n(0) = 0.0D0
        n(1) = 0.11670521452767D+04
        n(2) = -0.72421316703206D+06
        n(3) = -0.17073846940092D+02
        n(4) = 0.12020824702470D+05
        n(5) = -0.32325550322333D+07
        n(6) = 0.14915108613530D+02
        n(7) = -0.48232657361591D+04
        n(8) = 0.40511340542057D+06
        n(9) = -0.23855557567849D+00
        n(10) = 0.65017534844798D+03

        if (P < Pmin .or. P > Pc) then
            tsat = -1.0D0
            return
        end if

        beta = P**0.25D0
        E = beta**2 + n(3)*beta + n(6)
        F = n(1)*beta**2 + n(4)*beta + n(7)
        G = n(2)*beta**2 + n(5)*beta + n(8)
        D = 2.0D0 * G / (-F - sqrt(F**2 - 4.0D0*E*G))
        tsat = (n(10) + D - sqrt((n(10) + D)**2 - 4.0D0 * (n(9) + n(10)*D))) / 2.0D0
    end function tsat

    ! -------------------------------------------------------------------------
    ! Boundary Region 2 / 3: P = f(T) in MPa
    ! -------------------------------------------------------------------------
    double precision function p23(T)
        double precision, intent(in) :: T
        p23 = 0.34805185628969D+03 - 0.11671859879975D+01 * T + 0.10192970039326D-02 * T**2
    end function p23

    ! -------------------------------------------------------------------------
    ! Region 1 Properties Solver
    ! -------------------------------------------------------------------------
    subroutine solve_region1(T, P, v, h, s, cp, cv, w)
        double precision, intent(in) :: T, P
        double precision, intent(out) :: v, h, s, cp, cv, w
        double precision :: Tr, Pr, g, gp, gpp, gt, gtt, gpt
        integer :: i

        Tr = 1386.0D0 / T
        Pr = P / 16.53D0

        g = 0.0D0
        gp = 0.0D0
        gpp = 0.0D0
        gt = 0.0D0
        gtt = 0.0D0
        gpt = 0.0D0

        do i = 1, 34
            g = g + r1_n(i) * (7.1D0 - Pr)**r1_li(i) * (Tr - 1.222D0)**r1_lj(i)
            
            if (r1_li(i) > 0) then
                gp = gp - r1_n(i) * r1_li(i) * (7.1D0 - Pr)**(r1_li(i)-1) * (Tr - 1.222D0)**r1_lj(i)
            end if
            
            if (r1_li(i) > 1) then
                gpp = gpp + r1_n(i) * r1_li(i) * (r1_li(i)-1) * (7.1D0 - Pr)**(r1_li(i)-2) * (Tr - 1.222D0)**r1_lj(i)
            end if
            
            if (r1_lj(i) /= 0) then
                gt = gt + r1_n(i) * r1_lj(i) * (7.1D0 - Pr)**r1_li(i) * (Tr - 1.222D0)**(r1_lj(i)-1)
            end if
            
            if (r1_lj(i) /= 0 .and. r1_lj(i) /= 1) then
                gtt = gtt + r1_n(i) * r1_lj(i) * (r1_lj(i)-1) * (7.1D0 - Pr)**r1_li(i) * (Tr - 1.222D0)**(r1_lj(i)-2)
            end if
            
            if (r1_li(i) > 0 .and. r1_lj(i) /= 0) then
                gpt = gpt - r1_n(i) * r1_li(i) * r1_lj(i) * (7.1D0 - Pr)**(r1_li(i)-1) * (Tr - 1.222D0)**(r1_lj(i)-1)
            end if
        end do

        v = Pr * gp * R * T / P / 1000.0D0
        h = Tr * gt * R * T
        s = R * (Tr * gt - g)
        cp = -R * Tr**2 * gtt
        cv = R * (-Tr**2 * gtt + (gp - Tr * gpt)**2 / gpp)
        w = sqrt(R * T * 1000.0D0 * gp**2 / ((gp - Tr * gpt)**2 / (Tr**2 * gtt) - gpp))
    end subroutine solve_region1

    ! -------------------------------------------------------------------------
    ! Region 2cp0 Helper
    ! -------------------------------------------------------------------------
    subroutine region2_cp0(Tr, Pr, go, gop, gopp, got, gott, gopt)
        double precision, intent(in) :: Tr, Pr
        double precision, intent(out) :: go, gop, gopp, got, gott, gopt
        integer :: i

        go = log(Pr)
        gop = 1.0D0 / Pr
        gopp = -1.0D0 / Pr**2
        gopt = 0.0D0

        got = 0.0D0
        gott = 0.0D0

        do i = 1, 9
            go = go + r2_cp0_no(i) * Tr**r2_cp0_jo(i)
            got = got + r2_cp0_no(i) * r2_cp0_jo(i) * Tr**(r2_cp0_jo(i)-1)
            gott = gott + r2_cp0_no(i) * r2_cp0_jo(i) * (r2_cp0_jo(i)-1) * Tr**(r2_cp0_jo(i)-2)
        end do
    end subroutine region2_cp0

    ! -------------------------------------------------------------------------
    ! Region 2 Properties Solver
    ! -------------------------------------------------------------------------
    subroutine solve_region2(T, P, v, h, s, cp, cv, w)
        double precision, intent(in) :: T, P
        double precision, intent(out) :: v, h, s, cp, cv, w
        double precision :: Tr, Pr, go, gop, gopp, got, gott, gopt
        double precision :: gr, grp, grpp, grt, grtt, grpt
        integer :: i

        Tr = 540.0D0 / T
        Pr = P / 1.0D0

        call region2_cp0(Tr, Pr, go, gop, gopp, got, gott, gopt)

        gr = 0.0D0
        grp = 0.0D0
        grpp = 0.0D0
        grt = 0.0D0
        grtt = 0.0D0
        grpt = 0.0D0

        do i = 1, 43
            gr = gr + r2_n(i) * Pr**r2_li(i) * (Tr - 0.5D0)**r2_lj(i)
            
            if (r2_li(i) > 0) then
                grp = grp + r2_n(i) * r2_li(i) * Pr**(r2_li(i)-1) * (Tr - 0.5D0)**r2_lj(i)
            end if
            
            if (r2_li(i) > 1) then
                grpp = grpp + r2_n(i) * r2_li(i) * (r2_li(i)-1) * Pr**(r2_li(i)-2) * (Tr - 0.5D0)**r2_lj(i)
            end if
            
            if (r2_lj(i) > 0) then
                grt = grt + r2_n(i) * r2_lj(i) * Pr**r2_li(i) * (Tr - 0.5D0)**(r2_lj(i)-1)
            end if
            
            if (r2_lj(i) > 1) then
                grtt = grtt + r2_n(i) * r2_lj(i) * (r2_lj(i)-1) * Pr**r2_li(i) * (Tr - 0.5D0)**(r2_lj(i)-2)
            end if
            
            if (r2_li(i) > 0 .and. r2_lj(i) > 0) then
                grpt = grpt + r2_n(i) * r2_li(i) * r2_lj(i) * Pr**(r2_li(i)-1) * (Tr - 0.5D0)**(r2_lj(i)-1)
            end if
        end do

        v = Pr * (gop + grp) * R * T / P / 1000.0D0
        h = Tr * (got + grt) * R * T
        s = R * (Tr * (got + grt) - (go + gr))
        cp = -R * Tr**2 * (gott + grtt)
        cv = R * (-Tr**2 * (gott + grtt) - (1.0D0 + Pr*grp - Tr*Pr*grpt)**2 / (1.0D0 - Pr**2 * grpp))
        w = (R*T*1000.0D0*(1.0D0 + 2.0D0*Pr*grp + Pr**2 * grp**2) / &
            (1.0D0 - Pr**2 * grpp + (1.0D0 + Pr*grp - Tr*Pr*grpt)**2 / (Tr**2 * (gott + grtt))))**0.5D0
    end subroutine solve_region2

    ! -------------------------------------------------------------------------
    ! Region 3 Properties Solver (Inputs: density rho in kg/m3 and T in K)
    ! -------------------------------------------------------------------------
    subroutine solve_region3(rho, T, P, v, h, s, cp, cv, w)
        double precision, intent(in) :: rho, T
        double precision, intent(out) :: P, v, h, s, cp, cv, w
        double precision :: d, Tr, g, gd, gdd, gt, gtt, gdt
        integer :: i

        d = rho / rhoc
        Tr = Tc / T

        ! 1.0658070028513D0 is the coefficient of log(d) in standard Region 3 Helmholtz
        g = 1.0658070028513D0 * log(d)
        gd = 1.0658070028513D0 / d
        gdd = -1.0658070028513D0 / d**2
        gt = 0.0D0
        gtt = 0.0D0
        gdt = 0.0D0

        do i = 1, 39
            g = g + r3_n(i) * d**r3_li(i) * Tr**r3_lj(i)
            
            if (r3_li(i) > 0) then
                gd = gd + r3_n(i) * r3_li(i) * d**(r3_li(i)-1) * Tr**r3_lj(i)
            end if
            
            if (r3_li(i) > 1) then
                gdd = gdd + r3_n(i) * r3_li(i) * (r3_li(i)-1) * d**(r3_li(i)-2) * Tr**r3_lj(i)
            end if
            
            if (r3_lj(i) > 0) then
                gt = gt + r3_n(i) * r3_lj(i) * d**r3_li(i) * Tr**(r3_lj(i)-1)
            end if
            
            if (r3_lj(i) > 1) then
                gtt = gtt + r3_n(i) * r3_lj(i) * (r3_lj(i)-1) * d**r3_li(i) * Tr**(r3_lj(i)-2)
            end if
            
            if (r3_li(i) > 0 .and. r3_lj(i) > 0) then
                gdt = gdt + r3_n(i) * r3_li(i) * r3_lj(i) * d**(r3_li(i)-1) * Tr**(r3_lj(i)-1)
            end if
        end do

        P = d * gd * R * T * rho / 1000.0D0
        v = 1.0D0 / rho
        h = R * T * (Tr * gt + d * gd)
        s = R * (Tr * gt - g)
        cp = R * (-Tr**2 * gtt + (d*gd - d*Tr*gdt)**2 / (2.0D0*d*gd + d**2 * gdd))
        cv = -R * Tr**2 * gtt
        w = sqrt(R * T * 1000.0D0 * (2.0D0*d*gd + d**2 * gdd - (d*gd - d*Tr*gdt)**2 / (Tr**2 * gtt)))
    end subroutine solve_region3

    ! -------------------------------------------------------------------------
    ! Region 3 Solver from P and T (finds density first using bisection)
    ! -------------------------------------------------------------------------
    subroutine solve_region3_PT(T, P, v, h, s, cp, cv, w)
        double precision, intent(in) :: T, P
        double precision, intent(out) :: v, h, s, cp, cv, w
        double precision :: rho_min, rho_max, rho_mid, P_calc
        double precision :: d_dummy, h_dummy, s_dummy, cp_dummy, cv_dummy, w_dummy
        integer :: iter

        ! Bisection search for density rho
        rho_min = 1.0D0
        rho_max = 1200.0D0

        do iter = 1, 50
            rho_mid = 0.5D0 * (rho_min + rho_max)
            call solve_region3(rho_mid, T, P_calc, d_dummy, h_dummy, s_dummy, cp_dummy, cv_dummy, w_dummy)
            if (P_calc < P) then
                if (T < Tc) then
                    ! Saturated liquid region has higher density, so higher density yields higher pressure
                    rho_min = rho_mid
                else
                    rho_min = rho_mid
                end if
            else
                rho_max = rho_mid
            end if
        end do

        rho_mid = 0.5D0 * (rho_min + rho_max)
        call solve_region3(rho_mid, T, P_calc, v, h, s, cp, cv, w)
        v = 1.0D0 / rho_mid
    end subroutine solve_region3_PT

    ! -------------------------------------------------------------------------
    ! Determine Region ID based on T (K) and P (MPa)
    ! -------------------------------------------------------------------------
    integer function get_region(T, P)
        double precision, intent(in) :: T, P
        double precision :: ps

        if (T < 623.15D0) then
            ps = psat(T)
            if (P >= ps) then
                get_region = 1
            else
                get_region = 2
            end if
        else if (T <= 863.15D0) then
            if (P >= p23(T)) then
                get_region = 3
            else
                get_region = 2
            end if
        else
            get_region = 2
        end if
    end function get_region

    ! -------------------------------------------------------------------------
    ! Master Properties Solver from Pressure and Temperature
    ! -------------------------------------------------------------------------
    subroutine solve_PT(P, T, v, h, s, cp, cv, w, region_id)
        double precision, intent(in) :: P, T
        double precision, intent(out) :: v, h, s, cp, cv, w
        integer, intent(out) :: region_id

        region_id = get_region(T, P)
        
        select case (region_id)
        case (1)
            call solve_region1(T, P, v, h, s, cp, cv, w)
        case (2)
            call solve_region2(T, P, v, h, s, cp, cv, w)
        case (3)
            call solve_region3_PT(T, P, v, h, s, cp, cv, w)
        case default
            call solve_region2(T, P, v, h, s, cp, cv, w)
        end select
    end subroutine solve_PT

end module iapws_if97_mod

! ============================================================================
! Main Program for iapws_if97_calc
! ============================================================================
program iapws_if97_calc
    use iapws_if97_mod
    implicit none

    integer :: input_type
    double precision :: val1, val2
    
    ! Conversions
    double precision :: T_c, T_k, P_bar, P_mpa
    double precision :: v, h, s, u, cp, cv, w, density, quality
    character(len=30) :: phase_name
    integer :: region_id

    ! Saturation checks
    double precision :: t_sat, h_f, h_g, s_f, s_g, v_f, v_g, cp_f, cp_g, cv_f, cv_g, w_f, w_g
    double precision :: t_min, t_max, t_mid, h_mid, s_mid
    integer :: iter, reg_dummy

    ! Read inputs from stdin
    read(*,*) input_type
    read(*,*) val1
    read(*,*) val2

    quality = -1.0D0
    phase_name = "Single Phase"

    select case (input_type)

    ! -------------------------------------------------------------------------
    ! 1. Pressure & Temperature (val1 = P [bar], val2 = T [C])
    ! -------------------------------------------------------------------------
    case (1)
        P_bar = val1
        P_mpa = P_bar / 10.0D0
        T_c = val2
        T_k = T_c + 273.15D0

        if (T_k < 273.15D0 .or. T_k > 1073.15D0) then
            write(*,*) 'ERROR: Temperature must be between 0 C and 800 C'
            stop
        end if
        if (P_mpa < Pmin .or. P_mpa > 100.0D0) then
            write(*,*) 'ERROR: Pressure must be between 0.006 bar and 1000 bar'
            stop
        end if

        call solve_PT(P_mpa, T_k, v, h, s, cp, cv, w, region_id)
        
        if (region_id == 1) then
            phase_name = "Subcooled Liquid"
            quality = 0.0D0
        else if (region_id == 2) then
            phase_name = "Superheated Vapor"
            quality = 1.0D0
        else if (region_id == 3) then
            phase_name = "Supercritical Fluid"
            quality = -1.0D0
        end if

    ! -------------------------------------------------------------------------
    ! 2. Pressure & Quality (val1 = P [bar], val2 = x)
    ! -------------------------------------------------------------------------
    case (2)
        P_bar = val1
        P_mpa = P_bar / 10.0D0
        quality = val2

        if (P_mpa < Pmin .or. P_mpa > Pc) then
            write(*,*) 'ERROR: Saturated state requires P between 0.006 bar and 220.64 bar'
            stop
        end if
        if (quality < 0.0D0 .or. quality > 1.0D0) then
            write(*,*) 'ERROR: Quality x must be between 0 and 1'
            stop
        end if

        T_k = tsat(P_mpa)
        T_c = T_k - 273.15D0
        region_id = 4

        call solve_PT(P_mpa, T_k, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
        ! For vapor line, we evaluate Region 2 at Saturation Temperature and Pressure
        call solve_region2(T_k, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)

        ! Blend properties
        v = (1.0D0 - quality) * v_f + quality * v_g
        h = (1.0D0 - quality) * h_f + quality * h_g
        s = (1.0D0 - quality) * s_f + quality * s_g
        cp = (1.0D0 - quality) * cp_f + quality * cp_g
        cv = (1.0D0 - quality) * cv_f + quality * cv_g
        w = (1.0D0 - quality) * w_f + quality * w_g

        if (quality == 0.0D0) then
            phase_name = "Saturated Liquid"
        else if (quality == 1.0D0) then
            phase_name = "Saturated Vapor"
        else
            phase_name = "Saturated Mixture"
        end if

    ! -------------------------------------------------------------------------
    ! 3. Temperature & Quality (val1 = T [C], val2 = x)
    ! -------------------------------------------------------------------------
    case (3)
        T_c = val1
        T_k = T_c + 273.15D0
        quality = val2

        if (T_k < 273.15D0 .or. T_k >= Tc) then
            write(*,*) 'ERROR: Saturated state requires T between 0 C and 373.94 C'
            stop
        end if
        if (quality < 0.0D0 .or. quality > 1.0D0) then
            write(*,*) 'ERROR: Quality x must be between 0 and 1'
            stop
        end if

        P_mpa = psat(T_k)
        P_bar = P_mpa * 10.0D0
        region_id = 4

        call solve_PT(P_mpa, T_k, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
        call solve_region2(T_k, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)

        ! Blend properties
        v = (1.0D0 - quality) * v_f + quality * v_g
        h = (1.0D0 - quality) * h_f + quality * h_g
        s = (1.0D0 - quality) * s_f + quality * s_g
        cp = (1.0D0 - quality) * cp_f + quality * cp_g
        cv = (1.0D0 - quality) * cv_f + quality * cv_g
        w = (1.0D0 - quality) * w_f + quality * w_g

        if (quality == 0.0D0) then
            phase_name = "Saturated Liquid"
        else if (quality == 1.0D0) then
            phase_name = "Saturated Vapor"
        else
            phase_name = "Saturated Mixture"
        end if

    ! -------------------------------------------------------------------------
    ! 4. Pressure & Enthalpy (val1 = P [bar], val2 = h [kJ/kg])
    ! -------------------------------------------------------------------------
    case (4)
        P_bar = val1
        P_mpa = P_bar / 10.0D0
        h = val2

        if (P_mpa < Pmin .or. P_mpa > 100.0D0) then
            write(*,*) 'ERROR: Pressure must be between 0.006 bar and 1000 bar'
            stop
        end if

        ! 1. Check if we are in two-phase region (Region 4)
        if (P_mpa <= Pc) then
            T_sat = tsat(P_mpa)
            call solve_PT(P_mpa, T_sat, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
            call solve_region2(T_sat, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)

            if (h >= h_f .and. h <= h_g) then
                ! Saturation mixture
                quality = (h - h_f) / (h_g - h_f)
                T_k = T_sat
                T_c = T_k - 273.15D0
                region_id = 4
                
                v = (1.0D0 - quality) * v_f + quality * v_g
                s = (1.0D0 - quality) * s_f + quality * s_g
                cp = (1.0D0 - quality) * cp_f + quality * cp_g
                cv = (1.0D0 - quality) * cv_f + quality * cv_g
                w = (1.0D0 - quality) * w_f + quality * w_g
                
                phase_name = "Saturated Mixture"
            end if
        end if

        ! 2. Single-phase region (either P > Pc or h is outside saturation limits)
        if (quality < 0.0D0) then
            ! Check if it is liquid (h < h_f) or vapor (h > h_g)
            if (P_mpa <= Pc .and. h < h_f) then
                phase_name = "Subcooled Liquid"
                t_min = 273.15D0
                t_max = T_sat
            else if (P_mpa <= Pc .and. h > h_g) then
                phase_name = "Superheated Vapor"
                t_min = T_sat
                t_max = 1073.15D0
            else
                phase_name = "Supercritical Fluid"
                t_min = 273.15D0
                t_max = 1073.15D0
            end if

            ! Bisection search for Temperature
            do iter = 1, 40
                t_mid = 0.5D0 * (t_min + t_max)
                call solve_PT(P_mpa, t_mid, v, h_mid, s, cp, cv, w, region_id)
                if (h_mid < 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
            call solve_PT(P_mpa, T_k, v, h, s, cp, cv, w, region_id)
            if (region_id == 1) quality = 0.0D0
            if (region_id == 2) quality = 1.0D0
        end if

    ! -------------------------------------------------------------------------
    ! 5. Pressure & Entropy (val1 = P [bar], val2 = s [kJ/kg-K])
    ! -------------------------------------------------------------------------
    case (5)
        P_bar = val1
        P_mpa = P_bar / 10.0D0
        s = val2

        if (P_mpa < Pmin .or. P_mpa > 100.0D0) then
            write(*,*) 'ERROR: Pressure must be between 0.006 bar and 1000 bar'
            stop
        end if

        ! 1. Check if we are in two-phase region (Region 4)
        if (P_mpa <= Pc) then
            T_sat = tsat(P_mpa)
            call solve_PT(P_mpa, T_sat, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
            call solve_region2(T_sat, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)

            if (s >= s_f .and. s <= s_g) then
                ! Saturation mixture
                quality = (s - s_f) / (s_g - s_f)
                T_k = T_sat
                T_c = T_k - 273.15D0
                region_id = 4
                
                v = (1.0D0 - quality) * v_f + quality * v_g
                h = (1.0D0 - quality) * h_f + quality * h_g
                cp = (1.0D0 - quality) * cp_f + quality * cp_g
                cv = (1.0D0 - quality) * cv_f + quality * cv_g
                w = (1.0D0 - quality) * w_f + quality * w_g
                
                phase_name = "Saturated Mixture"
            end if
        end if

        ! 2. Single-phase region (either P > Pc or s is outside saturation limits)
        if (quality < 0.0D0) then
            if (P_mpa <= Pc .and. s < s_f) then
                phase_name = "Subcooled Liquid"
                t_min = 273.15D0
                t_max = T_sat
            else if (P_mpa <= Pc .and. s > s_g) then
                phase_name = "Superheated Vapor"
                t_min = T_sat
                t_max = 1073.15D0
            else
                phase_name = "Supercritical Fluid"
                t_min = 273.15D0
                t_max = 1073.15D0
            end if

            ! Bisection search for Temperature
            do iter = 1, 40
                t_mid = 0.5D0 * (t_min + t_max)
                call solve_PT(P_mpa, t_mid, v, h, s_mid, cp, cv, w, region_id)
                if (s_mid < 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
            call solve_PT(P_mpa, T_k, v, h, s, cp, cv, w, region_id)
            if (region_id == 1) quality = 0.0D0
            if (region_id == 2) quality = 1.0D0
        end if

    case default
        write(*,*) 'ERROR: Unknown input combination'
        stop
    end select

    ! Density conversion: v in m3/kg => rho = 1/v kg/m3
    density = 1.0D0 / v

    ! Specific Internal Energy u = h - P * v
    ! P is in bar. 1 bar = 100 kPa.
    ! P * v is in (100 kPa) * (m3/kg) = 100 kJ/kg.
    u = h - P_bar * 100.0D0 * v

    ! Print output in standard key-value format for PHP wrapper
    write(*,'(A,A)') 'STATE=', trim(phase_name)
    write(*,'(A,I1)') 'REGION=', region_id
    write(*,'(A,F14.4)') 'T=', T_c
    write(*,'(A,F14.4)') 'P=', P_bar
    write(*,'(A,F14.6)') 'V=', v
    write(*,'(A,F14.4)') 'RHO=', density
    write(*,'(A,F14.4)') 'H=', h
    write(*,'(A,F14.6)') 'S=', s
    write(*,'(A,F14.4)') 'U=', u
    write(*,'(A,F14.4)') 'CP=', cp
    write(*,'(A,F14.4)') 'CV=', cv
    write(*,'(A,F14.2)') 'SOUND=', w
    write(*,'(A,F14.4)') 'QUALITY=', quality

end program iapws_if97_calc
