he_prop Subroutine

subroutine he_prop(ro, T, cv, cp, mu, lambda, dPdT) ! TODO not elemental due to set_error 
    !! Viscosity and thermal conductivity correlations for helium. The critical-region
    !! enhancement term in lambda needs dP/dRo and dP/dT;
    real(dp), intent(in)  :: ro, T
    real(dp), intent(out) :: cv, cp, mu, lambda, dPdT
    real(dp) :: dPdRo, dTdRo

    real(dp), parameter :: CC(4) = [3.739232544_dp, -2.620316969e1_dp, 5.982252246e1_dp, -4.926397634e1_dp]
    real(dp), parameter :: Coef(11) = [1.862970530e-4_dp, -7.275964435e-7_dp, -1.427549651e-4_dp, &
                    3.290833592e-5_dp, -5.213335363e-8_dp, 4.492659933e-8_dp, -5.924416513e-9_dp, &
                    7.087321137e-6_dp, -6.013335678e-6_dp, 8.067145814e-7_dp,  3.995125013e-7_dp]
    real(dp) :: dRodP
    real(dp) :: ro_unit
    real(dp) :: x, EtaPrim0, BBB, CCC, DDD, EtaPrimE, Eta0, EtaE
    real(dp) :: Lambd0, LambdC, LambdE, SumCjSurTj
    real(dp) :: AA, BB, x0, E1, E2, Beta, Gama, Delta, roc, Tc, Pc, RR, m, k
    real(dp) :: DeltaT, DeltaRo, Eta, KT, rostar, LongTerm, KTprim, CoefCrit
    real(dp) :: W, hhh, dhhhdx
    integer  :: j

    call jacobian_roT(ro, T, cv, dPdT, dTdRo=dTdRo)
    dPdRo = -dTdRo*dPdT
    cp = cv + T*(dPdT**2)/(dPdRo*(ro**2))

    if (cp < 0.0_dp .or. cv < 0.0_dp) then
        call set_error('negative cp or cv in he_prop')
        return
    endif

    dRodP = 1.0_dp/dPdRo
    ro_unit = ro/1000.0_dp
    if (T <= 300.0_dp) then
        x = log(T)
    else
        x = 5.7037825_dp  ! log(300.0_dp)
    endif
    EtaPrim0 = -0.135311743_dp/x+1.00347841_dp+1.20654649_dp*x-0.149564551_dp*(x**2)+0.0125208416_dp*(x**3)
    BBB = -47.5295259_dp/x+87.6799309_dp-42.0741589_dp*x+8.33128289_dp*(x**2)-0.589252385_dp*(x**3)
    CCC = 547.309267_dp/x-904.870586_dp+431.404928_dp*x-81.4504854_dp*(x**2)+5.37008433_dp*(x**3)
    DDD = -1684.39324_dp/x+3331.08630_dp-1632.19172_dp*x+308.804413_dp*(x**2)-20.2936367_dp*(x**3)
    EtaPrimE = ro_unit*BBB+(ro_unit**2)*CCC+(ro_unit**3)*DDD
    if (EtaPrimE > 100.0_dp) then
        call set_error('EtaPrimE > 100 during viscosity calculation in he_prop')
        return
    endif
    if (T <= 100.0_dp) then
        mu = exp(EtaPrim0+EtaPrimE)
    else
        Eta0 = 196.0_dp*(T**0.71938_dp)*exp(12.451_dp/T-295.67_dp/(T**2)-4.1249_dp)
        if (T < 110.0_dp) then
            Eta0 = exp(EtaPrim0)+(Eta0-exp(EtaPrim0))*(T-100.0_dp)/10.0_dp
        endif
        EtaE = exp(EtaPrim0+EtaPrimE)-exp(EtaPrim0)
        mu = Eta0+EtaE
    endif
    mu = 1.0e-7_dp*mu

    AA = 2.7870034e-3_dp
    BB = 7.034007057e-1_dp
    SumCjSurTj = 0.0_dp
    do j = 1, 4
        SumCjSurTj = SumCjSurTj + CC(j)/(T**j)
    enddo
    Lambd0 = AA*(T**BB)*exp(SumCjSurTj)

    LambdC = 0.0_dp
    if (T >= 3.5_dp .and. T <= 12.0_dp) then
        x0 = 0.392_dp
        E1 = 2.8461_dp
        E2 = 0.27156_dp
        Beta = 0.3554_dp
        Gama = 1.1743_dp
        Delta = 4.304_dp
        roc = 69.158_dp
        Tc = 5.18992_dp
        Pc = 2.2746e5_dp
        RR = 4.633e-10_dp
        m = 6.6455255e-27_dp
        k = 1.38066e-23_dp
        DeltaT = abs(1.0_dp-T/Tc)
        DeltaRo = abs(1.0_dp-ro/roc)
        Eta = mu
        rostar = ro/roc

        ! sqrt(m*k)/(6*Pi*RR), as in REFPROP:
        CoefCrit = 3.4685233e-17_dp

        KT = dRodP/ro
        W = (DeltaT/0.2_dp)**2+(DeltaRo/0.25_dp)**2

        if (W < 1.0_dp) then
            x = DeltaT/(DeltaRo**(1.0_dp/Beta))
            hhh = E1*(1.0_dp+x/x0)*(1.0_dp+E2*(1.0_dp+x/x0)**(2.0_dp*Beta))**((Gama-1.0_dp)/(2.0_dp*Beta))

            dhhhdx = E1*(Gama-1.0_dp)*E2*(1.0_dp+x/x0)**(2.0_dp*Beta)*(1.0_dp+ &
                E2*(1.0_dp+x/x0)**(2.0_dp*Beta))**((Gama-1.0_dp)/(2.0_dp*Beta)-1.0_dp)+E1*(1.0_dp+ &
                E2*(1.0_dp+x/x0)**(2.0_dp*Beta))**((Gama-1.0_dp)/(2.0_dp*Beta))
            dhhhdx = dhhhdx/x0

            LongTerm = (DeltaRo**(Delta-1.0_dp))*(Delta*hhh-x*dhhhdx/Beta)
            KTprim = 1.0_dp/(rostar**2)/Pc/LongTerm
            KT = W*KT+(1.0_dp-W)*KTprim
        endif

        LambdC = CoefCrit*sqrt(KT*(T**3)/ro)*(dPdT**2)*exp(-18.66_dp*(DeltaT**2)-4.25_dp*(DeltaRo**4))/Eta
    endif

    rostar = ro/68.0_dp
    LambdE=(Coef(1)+Coef(2)*T+Coef(3)*T**(1.0_dp/3.0_dp)+Coef(4)*T**(2.0_dp/3.0_dp))*ro+&
           (Coef(5)+Coef(6)*T**(1.0_dp/3.0_dp)+Coef(7)*T**(2.0_dp/3.0_dp))*(ro**3)+&
           (Coef(8)+Coef(9)*T**(1.0_dp/3.0_dp)+Coef(10)*T**(2.0_dp/3.0_dp)+Coef(11)/T)*(ro**2)*log(rostar)

    lambda = Lambd0 + LambdC + LambdE
end subroutine he_prop