SUBROUTINE TYPE229 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) 
!************************************************************************
! Object: Substance Properties
! Simulation Studio Model: Type229
! 
! Author: Michael J Wagner
! Editor: 
! Date:	 February 02, 2009 last modified: October 26, 2009
! COPYRIGHT 2009 NATIONAL RENEWABLE ENERGY LABORATORY

! Parameters
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Input Units   |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  Fluid_number                  |  Integer flag denoting HTF fluid                                |  none          |  none
!  2|  LU_FL                         |  Fluid data logical unit                                        |  None          |  none
!  3|  LuFlEr                        |  Fluid property error file unit                                 |  None          |  none


! Inputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Input Units   |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  Temperature                   |  Temperature of the HTF                                         |  C             |  K

! Outputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Output Units  |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  c                             |  Fluid specific heat                                            |  kJ/kg-K       |  kJ/kg-K
!  2|  rho                           |  Fluid density                                                  |  kg/m^3        |  kg/m^3
!  3|  mu                            |  Fluid viscosity                                                |  Pa-s          |  Pa-s
!  4|  k                             |  Fluid thermal conductivity                                     |  W/m-K         |  W/m-K


!************************************************************************

!    TRNSYS acess functions (allow to acess TIME etc.) 
USE TrnsysConstants
USE TrnsysFunctions
!    Special access functions for fluid properties
use global_props      

!-----------------------------------------------------------------------------------------------------------------------
!    REQUIRED BY THE MULTI-DLL VERSION OF TRNSYS
!DEC$ATTRIBUTES DLLEXPORT :: TYPE229				!SET THE CORRECT TYPE NUMBER HERE
!-----------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------
!    TRNSYS DECLARATIONS
IMPLICIT NONE			!REQUIRES THE USER TO DEFINE ALL VARIABLES BEFORE USING THEM

DOUBLE PRECISION XIN	!THE ARRAY FROM WHICH THE INPUTS TO THIS TYPE WILL BE RETRIEVED
DOUBLE PRECISION OUT	!THE ARRAY WHICH WILL BE USED TO STORE THE OUTPUTS FROM THIS TYPE
DOUBLE PRECISION TIME	!THE CURRENT SIMULATION TIME - YOU MAY USE THIS VARIABLE BUT DO NOT SET IT!
DOUBLE PRECISION PAR	!THE ARRAY FROM WHICH THE PARAMETERS FOR THIS TYPE WILL BE RETRIEVED
DOUBLE PRECISION STORED !THE STORAGE ARRAY FOR HOLDING VARIABLES FROM TIMESTEP TO TIMESTEP
DOUBLE PRECISION T		!AN ARRAY CONTAINING THE RESULTS FROM THE DIFFERENTIAL EQUATION SOLVER
DOUBLE PRECISION DTDT	!AN ARRAY CONTAINING THE DERIVATIVES TO BE PASSED TO THE DIFF.EQ. SOLVER
INTEGER*4 INFO(15)		!THE INFO ARRAY STORES AND PASSES VALUABLE INFORMATION TO AND FROM THIS TYPE
INTEGER*4 NP,NI,NOUT,ND	!VARIABLES FOR THE MAXIMUM NUMBER OF PARAMETERS,INPUTS,OUTPUTS AND DERIVATIVES
INTEGER*4 NPAR,NIN,NDER	!VARIABLES FOR THE CORRECT NUMBER OF PARAMETERS,INPUTS,OUTPUTS AND DERIVATIVES
INTEGER*4 IUNIT,ITYPE	!THE UNIT NUMBER AND TYPE NUMBER FOR THIS COMPONENT
INTEGER*4 ICNTRL		!AN ARRAY FOR HOLDING VALUES OF CONTROL FUNCTIONS WITH THE NEW SOLVER
INTEGER*4 NSTORED		!THE NUMBER OF VARIABLES THAT WILL BE PASSED INTO AND OUT OF STORAGE
CHARACTER*3 OCHECK		!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE OUTPUTS
CHARACTER*3 YCHECK		!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE INPUTS
! ----------------------------------------------------------------------------------------------------------------------

! ----------------------------------------------------------------------------------------------------------------------
!    USER DECLARATIONS - SET THE MAXIMUM NUMBER OF PARAMETERS (NP), INPUTS (NI),
!    OUTPUTS (NOUT), AND DERIVATIVES (ND) THAT MAY BE SUPPLIED FOR THIS TYPE
PARAMETER (NP=3,NI=1,NOUT=4,ND=0,NSTORED=1)
! ----------------------------------------------------------------------------------------------------------------------

! ----------------------------------------------------------------------------------------------------------------------
!    REQUIRED TRNSYS DIMENSIONS
DIMENSION XIN(NI),OUT(NOUT),PAR(NP),YCHECK(NI),OCHECK(NOUT),STORED(NSTORED),T(ND),DTDT(ND)
INTEGER NITEMS
! ----------------------------------------------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------------------------------------------
!    ADD DECLARATIONS AND DEFINITIONS FOR THE USER-VARIABLES HERE


!    PARAMETERS
DOUBLE PRECISION Fluid_number 
integer::LU_FL

!    INPUTS
DOUBLE PRECISION Temperature

!     Outputs
DOUBLE PRECISION c, rho, mu, k

!     Locals
DOUBLE PRECISION specheat, Density, Viscosity, Conductivity 
!      DOUBLE PRECISION,ALLOCATABLE:: fprop(:,:)
!      CHARACTER*40 Coolist(2)
CHARACTER checkname*500, tc*1
INTEGER ios,ntot, nufl, i, j, m, actb  !fl_bounds(25)
LOGICAL Is_there

      
! ----------------------------------------------------------------------------------------------------------------------
!       READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
Fluid_number=PAR(1)
LU_FL=int(PAR(2))
LuFlEr = int(par(3)) !Note LuFlEr is stored globally in the global_props module
ttime=time  !Store time as a global variable
! ----------------------------------------------------------------------------------------------------------------------
!    RETRIEVE THE CURRENT VALUES OF THE INPUTS TO THIS MODEL FROM THE XIN ARRAY IN SEQUENTIAL ORDER

Temperature=XIN(1)
IUNIT=INFO(1)
ITYPE=INFO(2)

! ----------------------------------------------------------------------------------------------------------------------
!    SET THE VERSION INFORMATION FOR TRNSYS
IF(INFO(7).EQ.-2) THEN
    INFO(12)=16
    RETURN 1
ENDIF
! ----------------------------------------------------------------------------------------------------------------------

! ----------------------------------------------------------------------------------------------------------------------
!    DO ALL THE VERY LAST CALL OF THE SIMULATION MANIPULATIONS HERE
IF (INFO(8).EQ.-1) THEN
    if(allocated(fprop)) deallocate(fprop)

        inquire(unit=LuFlEr,opened=is_there)
        if(is_there) close(LuFlEr)

    RETURN 1
ENDIF
! ----------------------------------------------------------------------------------------------------------------------

! ----------------------------------------------------------------------------------------------------------------------
!    PERFORM ANY 'AFTER-ITERATION' MANIPULATIONS THAT ARE REQUIRED HERE
!    e.g. save variables to storage array for the next timestep
IF (INFO(13).GT.0) THEN
    STORED(1)=Gjsav
    CALL setStorageVars(STORED,NSTORED,INFO)
    RETURN 1
ENDIF

! ----------------------------------------------------------------------------------------------------------------------

! ----------------------------------------------------------------------------------------------------------------------
!    DO ALL THE VERY FIRST CALL OF THE SIMULATION MANIPULATIONS HERE
IF (INFO(7).EQ.-1) THEN

    !       SET SOME INFO ARRAY VARIABLES TO TELL THE TRNSYS ENGINE HOW THIS TYPE IS TO WORK
    INFO(6)=NOUT				
    INFO(9)=1				
    INFO(10)=0	!STORAGE FOR VERSION 16 HAS BEEN CHANGED				

    !       SET THE REQUIRED NUMBER OF INPUTS, PARAMETERS AND DERIVATIVES THAT THE USER SHOULD SUPPLY IN THE INPUT FILE
    !       IN SOME CASES, THE NUMBER OF VARIABLES MAY DEPEND ON THE VALUE OF PARAMETERS TO THIS MODEL....
    NIN=NI
    NPAR=NP
    NDER=ND

    !       CALL THE TYPE CHECK SUBROUTINE TO COMPARE WHAT THIS COMPONENT REQUIRES TO WHAT IS SUPPLIED IN 
    !       THE TRNSYS INPUT FILE
    CALL TYPECK(1,INFO,NIN,NPAR,NDER)

    !       SET THE NUMBER OF STORAGE SPOTS NEEDED FOR THIS COMPONENT
    CALL setStorageSize(NSTORED,INFO)

    !MJW 7/09 Open the warnings file
    inquire(unit=LuFlEr,name=checkname,opened=is_there)
    !open the file for writing
    if(.not.is_there) then
        open(unit=LuFlEr,file=trim(checkname),status="REPLACE",position="REWIND")
    endif


    !       RETURN TO THE CALLING PROGRAM
    RETURN 1

ENDIF
! ----------------------------------------------------------------------------------------------------------------------

! ----------------------------------------------------------------------------------------------------------------------
!    DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE - THERE ARE NO ITERATIONS AT THE INTIAL TIME
      IF (TIME .LT. (getSimulationStartTime() + getSimulationTimeStep()/2.D0)) THEN

!       SET THE UNIT NUMBER FOR FUTURE CALLS
         IUNIT=INFO(1)
         ITYPE=INFO(2)

        !On the initial timestep call, read in the fluid property file

        call readFluidPropFile(LU_FL)       
        fl_flag = 0. !Always reset the fluid warning flag to zero
        fl_ct = 0.  !and the warning counter too


!       PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE
!		 Specific heat capacity
			OUT(1)=1.53
!		 Density
			OUT(2)=1750.
!		 Viscosity
			OUT(3)=1e-5
!		 Conductivity
			OUT(4)=.5

!       PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL STORAGE VARIABLES HERE
	   STORED(1)=0. !Gjsav

!       PUT THE STORED ARRAY IN THE GLOBAL STORED ARRAY
         CALL setStorageVars(STORED,NSTORED,INFO)

!       RETURN TO THE CALLING PROGRAM
         RETURN 1

      ENDIF
! ----------------------------------------------------------------------------------------------------------------------

! ----------------------------------------------------------------------------------------------------------------------
!    *** ITS AN ITERATIVE CALL TO THIS COMPONENT ***
! ----------------------------------------------------------------------------------------------------------------------


! The following scheme is used to define the substance number.  Some fluids are
! blocked off for the trough model, and some are blocked off for future use. 
! Any fluid number specified greater than 35 will be a user-defined fluid.
!    1.) Air
!    2.) Stainless_AISI316
!    3.) Water (liquid)
!    4.) Steam
!    5.) CO2
!    6.) Salt (68% KCl, 32% MgCl2)
!    7.) Salt (8% NaF, 92% NaBF4)
!    8.) Salt (25% KF, 75% KBF4)
!    9.) Salt (31% RbF, 69% RbBF4)
!    10.) Salt (46.5% LiF, 11.5%NaF, 42%KF)
!    11.) Salt (49% LiF, 29% NaF, 29% ZrF4)
!    12.) Salt (58% KF, 42% ZrF4)
!    13.) Salt (58% LiCl, 42% RbCl)
!    14.) Salt (58% NaCl, 42% MgCl2)
!    15.) Salt (59.5% LiCl, 40.5% KCl)
!    16.) Salt (59.5% NaF, 40.5% ZrF4)
!    17.) Salt (60% NaNO3, 40% KNO3)
!    18.) Nitrate Salt**
!    19.) Caloria HT 43**
!    20.) Hitec XL**
!    21.) Therminol VP-1**
!    22.) Hitec**
!    23.) Dowtherm Q**
!    24.) Dowtherm RP**
!    26.) Argon (ideal gas properties)
!    27.) Hydrogen (ideal gas properties)
!    28.) -blank-
!    29.) -blank-
!    30.) -blank-
!    31.) -blank-
!    32.) -blank-
!    33.) -blank-
!    34.) -blank-
!    35.) -blank-
!    36+) User specified (lookup tables)

!     Note that the fluid properties are stored in the FPROP array in the following order:
!*****************************************************************
!|  #    |   1   |   2   |   3   |   4   |   5   |   6   |   7   |
!|-------|-------|-------|-------|-------|-------|-------|-------|
!| Name  |   T   |   Cp  |  rho  |   Mu  |   Nu  |   k   |   h   |
!| Units |   C   |kJ/kg-K| kg/m3 |  Pa-s |  m2-s | W/m-K |  J/kg |
!|-------|-------|-------|-------|-------|-------|-------|-------|
!|   1   |   :   |   :   |   :   |   :   |   :   |   :   |   :   |
!|   2   |   :   |   :   |   :   |   :   |   :   |   :   |   :   |
!*****************************************************************

     

!       GET THE STORED ARRAY FROM THE GLOBAL STORED ARRAY
         CALL getStorageVars(STORED,NSTORED,INFO)
	   Gjsav=STORED(1)  !Gjsav is a module variable, shared between this and the subfunctions


      c=specheat(Fluid_number,Temperature+273.15)
      rho=Density(Fluid_number,Temperature+273.15,dble(1.))
      mu=Viscosity(Fluid_number,Temperature+273.15)
      k=Conductivity(Fluid_number,Temperature+273.15)




!     Send the updated shared Gjsav back to storage
      STORED(1)=Gjsav
      call setStorageVars(STORED,NSTORED,INFO)



    !-----------------------------------------------------------------------------------------------------------------------

    !    SET THE OUTPUTS FROM THIS MODEL IN SEQUENTIAL ORDER AND GET OUT

    !Specific heat capacity
    OUT(1)=c
    !Density
    OUT(2)=rho
    !Viscosity
    OUT(3)=mu
    !Conductivity
    OUT(4)=k

!-----------------------------------------------------------------------------------------------------------------------
!    EVERYTHING IS DONE - RETURN FROM THIS SUBROUTINE AND MOVE ON
      RETURN 1
      END
!-----------------------------------------------------------------------------------------------------------------------



double precision function Density(fnumd,T,P)
use global_props

implicit none
!This function accepts as inputs temperature [K] and pressure [Pa]
!This function outputs in units of [kg/m^3]
double precision::v,R_air,xlo,xhi, Dens_fluid
double precision,intent(in)::T,P,fnumd
double precision,dimension(size(fprop(1,:)))::dxx,dyy !Create dummy arrays
integer::fnum,lb,ub,dum,t_warn
Density=1.
fnum=nint(fnumd)

select case(fnum)
case(1)   !    1.) Air
R_air = 287.                            !Gas constant [J/kg-K]
  v = R_air*T/P
Density = dmax1(1./v,1.e-6)
case(2)   !    2.) Stainless_AISI316
  Density=8349.38 - 0.341708*T - 0.0000865128*T**2  !EES
case(3)   !    3.) Water (liquid)
  continue 
case(4)   !    4.) Steam
  continue
case(5)   !    5.) CO2
  continue
case(6)   !    6.) Salt (68% KCl, 32% MgCl2)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (68% KCl, 32% MgCl2)")
Density = 1E-10*T**3 - 3E-07*T**2 - 0.4739*T + 2384.2
case(7)   !    7.) Salt (8% NaF, 92% NaBF4)
xlo=658.; xhi=969.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (8% NaF, 92% NaBF4)")
Density = 8E-09*T**3 - 2E-05*T**2 - 0.6867*T + 2438.5
case(8)   !    8.) Salt (25% KF, 75% KBF4)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (25% KF, 75% KBF4)")
Density = 2E-08*T**3 - 6E-05*T**2 - 0.7701*T + 2466.1
case(9)   !    9.) Salt (31% RbF, 69% RbBF4)
xlo=715.; xhi=1343.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (31% RbF, 69% RbBF4)")
Density = -1E-08*T**3 + 4E-05*T**2 - 1.0836*T + 3242.6
case(10)   !    10.) Salt (46.5% LiF, 11.5%NaF, 42%KF)
xlo=727.; xhi=1843.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (46.5% LiF, 11.5%NaF, 42%KF)")
Density =  -2E-09*T**3 + 1E-05*T**2 - 0.7427*T + 2734.7
case(11)   !    11.) Salt (49% LiF, 29% NaF, 29% ZrF4)
xlo=709.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (49% LiF, 29% NaF, 29% ZrF4)")
Density = -2E-11*T**3 + 1E-07*T**2 - 0.5172*T + 3674.3
case(12)   !    12.) Salt (58% KF, 42% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% KF, 42% ZrF4)")
Density =  -6E-10*T**3 + 4E-06*T**2 - 0.8931*T + 3661.3
case(13)   !    13.) Salt (58% LiCl, 42% RbCl)
xlo=586.; xhi=1323.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% LiCl, 42% RbCl)")
Density = -8E-10*T**3 + 1E-06*T**2 - 0.689*T + 2929.5
case(14)   !    14.) Salt (58% NaCl, 42% MgCl2)
xlo=718.; xhi=1738.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% NaCl, 42% MgCl2)")
Density = -5E-09*T**3 + 2E-05*T**2 - 0.5298*T + 2444.1
case(15)   !    15.) Salt (59.5% LiCl, 40.5% KCl)
xlo=628.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% LiCl, 40.5% KCl)")
Density = 1E-09*T**3 - 5E-06*T**2 - 0.864*T + 2112.6
case(16)   !    16.) Salt (59.5% NaF, 40.5% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% NaF, 40.5% ZrF4)")
Density =  -5E-09*T**3 + 2E-05*T**2 - 0.9144*T + 3837.
case(17)   !    17.) Salt (60% NaNO3, 40% KNO3)
xlo=533.; xhi=866.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (60% NaNO3, 40% KNO3)")
Density = dmax1(-1E-07*T**3 + 0.0002*T**2 - 0.7875*T + 2299.4,1000.d0)
case(18:25) !  18-25.) Call trough properties
Density = Dens_fluid((T-273.15),int(Fnumd)) !Trough calcs take fluid properties in degC, not K
case(26) !Argon
Density = dmax1(P/(208.13*T),1.e-10)
case(27) !Hydrogen
Density = dmax1(P/(4124.*T),1.e-10)
case(28:35) 
continue !no informaion
case(36:) !Any integer greater than 35
!Call the user-defined property table
lb=fl_bounds(fnum-35)
ub=fl_bounds(fnum-35+1)-1
if(ub.lt.lb) ub=size(fprop(1,:))
dxx(:)=fprop(1,lb:ub)
dyy(:)=fprop(3,lb:ub)
call interp(T,size(dxx),dxx,dyy,Gjsav,Density)
if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
end select

end function

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

double precision function specheat(fnumd,T,P)
use global_props
implicit none
!This function accepts as inputs temperature [K] and pressure [Pa]
!This function outputs in units of [kJ/kg-K]
double precision::xlo,xhi, Cp_fluid
double precision,intent(in)::T,P,fnumd
double precision,dimension(size(fprop(1,:)))::dxx,dyy !Create dummy arrays
integer::fnum,lb,ub,dum,t_warn
specheat=1.
fnum=nint(fnumd)

select case(fnum)
case(1)   !    1.) Air
specheat = 1.03749 - 0.000305497*T + 7.49335E-07*T**2 - 3.39363E-10*T**3
case(2)   !    2.) Stainless_AISI316
specheat = 0.368455 + 0.000399548*T - 1.70558E-07*T**2 !EES
case(3)   !    3.) Water (liquid)
  continue 
case(4)   !    4.) Steam
  continue
case(5)   !    5.) CO2
  continue
case(6)   !    6.) Salt (68% KCl, 32% MgCl2)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (68% KCl, 32% MgCl2)")
specheat = 1.156
case(7)   !    7.) Salt (8% NaF, 92% NaBF4)
xlo=658.; xhi=969.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (8% NaF, 92% NaBF4)")
specheat = 1.507
case(8)   !    8.) Salt (25% KF, 75% KBF4)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (25% KF, 75% KBF4)")
specheat = 1.306
case(9)   !    9.) Salt (31% RbF, 69% RbBF4)
xlo=715.; xhi=1343.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (31% RbF, 69% RbBF4)")
specheat = 9.127
case(10)   !    10.) Salt (46.5% LiF, 11.5%NaF, 42%KF)
xlo=727.; xhi=1843.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (46.5% LiF, 11.5%NaF, 42%KF)")
specheat = 2.010
case(11)   !    11.) Salt (49% LiF, 29% NaF, 29% ZrF4)
xlo=709.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (49% LiF, 29% NaF, 29% ZrF4)")
specheat = 1.239
case(12)   !    12.) Salt (58% KF, 42% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% KF, 42% ZrF4)")
specheat = 1.051
case(13)   !    13.) Salt (58% LiCl, 42% RbCl)
xlo=586.; xhi=1323.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% LiCl, 42% RbCl)")
specheat = 8.918
case(14)   !    14.) Salt (58% NaCl, 42% MgCl2)
xlo=718.; xhi=1738.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% NaCl, 42% MgCl2)")
specheat = 1.080
case(15)   !    15.) Salt (59.5% LiCl, 40.5% KCl)
xlo=628.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% LiCl, 40.5% KCl)")
specheat = 1.202
case(16)   !    16.) Salt (59.5% NaF, 40.5% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% NaF, 40.5% ZrF4)")
specheat = 1.172
case(17)   !    17.) Salt (60% NaNO3, 40% KNO3)
xlo=533.; xhi=866.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (60% NaNO3, 40% KNO3)")
specheat = -1E-10*T**3 + 2E-07*T**2 + 5E-06*T + 1.4387
case(18:25) !  18-25.) trough properties
specheat = Cp_fluid((T-273.15),int(fnumd))/1000.d0
case(26)    ! Argon
specheat = 0.5203 !Cp only, Cv is different
case(27)    ! Hydrogen
specheat = dmin1(dmax1(-45.4022 + 0.690156*T - 0.00327354*T**2 + 0.00000817326*T**3 - 1.13234E-08*T**4 + 8.24995E-12*T**5 - 2.46804E-15*T**6,11.3d0),14.7d0)
case(28:35)
continue
case(36:) !Any integer greater than 35
!Call the user-defined property table
lb=fl_bounds(fnum-35)
ub=fl_bounds(fnum-35+1)-1
if(ub.lt.lb) ub=size(fprop(1,:))
dxx(:)=fprop(1,lb:ub)
dyy(:)=fprop(2,lb:ub)
call interp(T,size(dxx),dxx,dyy,Gjsav,specheat)
        if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
end select


end function

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

double precision function Cv(fnumd,T)
use global_props
implicit none
!This function accepts as inputs temperature [K] and pressure [Pa]
!This function outputs in units of kJ/kg-K
double precision::xlo,xhi
double precision,intent(in)::T,fnumd
double precision,dimension(size(fprop(1,:)))::dxx,dyy !Create dummy arrays
integer::fnum,lb,ub,dum,t_warn
Cv=1.
fnum=nint(fnumd)

select case(fnum)
case(1)     !Air
Cv = 0.750466 - 0.000305497*T + 7.49335E-07*T**2 - 3.39363E-10*T**3
case(2:25) 
continue !no information
case(26)
Cv = 0.3122  !Argon
case(27)     !Hydrogen
Cv = dmin1(dmax1(-49.5264 + 0.690156*T - 0.00327354*T**2 + 0.00000817326*T**3 - 1.13234E-08*T**4 + 8.24995E-12*T**5 - 2.46804E-15*T**6,7.2d0),10.6d0)
case(28:35)
continue !no information
case(36:) !Any integer greater than 35
!Call the user-defined property table
!lb=fl_bounds(fnum-35)
!ub=fl_bounds(fnum-35+1)-1
!if(ub.lt.lb) ub=size(fprop(1,:))
!dxx(:)=fprop(1,lb:ub)
!dyy(:)=fprop(2,lb:ub)
!call interp(T,size(dxx),dxx,dyy,Gjsav,Cv)
!        if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")

!NOTE: no column selected yet for Cv

end select


end function



!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


double precision function Viscosity(fnumd,T,P)
use global_props
implicit none
!This function accepts as inputs temperature [K] and pressure [Pa]
!This function outputs in units of [Pa-s]
double precision::Tx,xlo,xhi, visc_fluid
double precision,intent(in)::T,P,fnumd
double precision,dimension(size(fprop(1,:)))::dxx,dyy !Create dummy arrays
integer::fnum,lb,ub,dum,t_warn
viscosity=1.
fnum=nint(fnumd)

select case(fnum)
case(1)   !    1.) Air
Viscosity = dmax1(0.0000010765 + 7.15173E-08*T - 5.03525E-11*T**2 + 2.02799E-14*T**3,1.e-6)
case(2)   !    2.) Stainless_AISI316
  continue
case(3)   !    3.) Water (liquid)
  continue 
case(4)   !    4.) Steam
  continue
case(5)   !    5.) CO2
  continue
case(6)   !    6.) Salt (68% KCl, 32% MgCl2)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (68% KCl, 32% MgCl2)")
Viscosity = .0146*exp(2230./T)*0.001 !convert cP to kg/m-s
case(7)   !    7.) Salt (8% NaF, 92% NaBF4)
xlo=658.; xhi=969.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (8% NaF, 92% NaBF4)")
Viscosity = .0877*exp(2240./T)*0.001 !convert cP to kg/m-s
case(8)   !    8.) Salt (25% KF, 75% KBF4)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (25% KF, 75% KBF4)")
Viscosity = .0431*exp(3060./T)*0.001 !convert cP to kg/m-s
case(9)   !    9.) Salt (31% RbF, 69% RbBF4)
xlo=715.; xhi=1343.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (31% RbF, 69% RbBF4)")
Viscosity = .0009
case(10)   !    10.) Salt (46.5% LiF, 11.5%NaF, 42%KF)
xlo=727.; xhi=1843.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (46.5% LiF, 11.5%NaF, 42%KF)")
Viscosity = .0400*exp(4170./T)*0.001 !convert cP to kg/m-s
case(11)   !    11.) Salt (49% LiF, 29% NaF, 29% ZrF4)
xlo=709.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (49% LiF, 29% NaF, 29% ZrF4)")
Viscosity = .0069
case(12)   !    12.) Salt (58% KF, 42% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% KF, 42% ZrF4)")
Viscosity = .0159*exp(3179./T)*0.001 !convert cP to kg/m-s
case(13)   !    13.) Salt (58% LiCl, 42% RbCl)
xlo=586.; xhi=1323.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% LiCl, 42% RbCl)")
Viscosity = .0861*exp(2517./T)*0.001 !convert cP to kg/m-s          !
case(14)   !    14.) Salt (58% NaCl, 42% MgCl2)
xlo=718.; xhi=1738.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% NaCl, 42% MgCl2)")
Viscosity = .0286*exp(1441./T)*0.001 !convert cP to kg/m-s
case(15)   !    15.) Salt (59.5% LiCl, 40.5% KCl)
xlo=628.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% LiCl, 40.5% KCl)")
Viscosity = .0861*exp(2517./T)*0.001 !convert cP to kg/m-s          !
case(16)   !    16.) Salt (59.5% NaF, 40.5% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% NaF, 40.5% ZrF4)")
Viscosity = .0767*exp(3977./T)*0.001 !convert cP to kg/m-s
case(17)   !    17.) Salt (60% NaNO3, 40% KNO3)
xlo=533.; xhi=866.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (60% NaNO3, 40% KNO3)")
Tx=T-273.15  !This particular equation is in terms of degrees celsius
Viscosity = dmax1(-1.473302E-10*Tx**3 + 2.279989E-07*Tx**2 - 1.199514E-04*Tx + 2.270616E-02,.0001d0)
case(18:25) !  18-35.) trough properties
Viscosity = visc_fluid((T-273.15),int(fnumd))
case(26)   !Argon 
Viscosity = 4.4997e-6 + 6.38920E-08*T - 1.24550E-11*T**2
case(27)  !Hydrogen
Viscosity=0.00000231 + 2.37842E-08*T - 5.73624E-12*T**2
case(28:35)
continue  !no information
case(36:) !Any integer greater than 35
!Call the user-defined property table
lb=fl_bounds(fnum-35)
ub=fl_bounds(fnum-35+1)-1
if(ub.lt.lb) ub=size(fprop(1,:))
dxx(:)=fprop(1,lb:ub)
dyy(:)=fprop(4,lb:ub)
call interp(T,size(dxx),dxx,dyy,Gjsav,Viscosity)
if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
end select


continue
end function

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

double precision function Conductivity(fnumd,T,P)
use global_props	
implicit none
!This function accepts as inputs temperature [K] and pressure [Pa]
!This function outputs in units of [W/m-K]
double precision::xlo,xhi, k_fluid
double precision,intent(in)::T,P,fnumd
double precision,dimension(size(fprop(1,:)))::dxx,dyy !Create dummy arrays
integer::fnum,lb,ub,dum,t_warn
conductivity=1.
fnum=nint(fnumd)

select case(fnum)
case(1)   !    1.) Air
Conductivity = dmax1(0.00145453 + 0.0000872152*T - 2.20614E-08*T**2,1.e-4)
case(2)   !    2.) Stainless_AISI316
Conductivity = 3E-09*T**3 - 8E-06*T**2 + 0.0177*T + 7.7765
case(3)   !    3.) Water (liquid)
  continue 
case(4)   !    4.) Steam
  continue
case(5)   !    5.) CO2
  continue
case(6)   !    6.) Salt (68% KCl, 32% MgCl2)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (68% KCl, 32% MgCl2)")
Conductivity = 0.39
case(7)   !    7.) Salt (8% NaF, 92% NaBF4)
xlo=658.; xhi=969.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (8% NaF, 92% NaBF4)")
Conductivity = 0.5
case(8)   !    8.) Salt (25% KF, 75% KBF4)
xlo=699.; xhi=1691.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (25% KF, 75% KBF4)")
Conductivity = 0.4
case(9)   !    9.) Salt (31% RbF, 69% RbBF4)
xlo=715.; xhi=1343.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (31% RbF, 69% RbBF4)")
Conductivity = 0.28
case(10)   !    10.) Salt (46.5% LiF, 11.5%NaF, 42%KF)
xlo=727.; xhi=1843.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (46.5% LiF, 11.5%NaF, 42%KF)")
Conductivity = 0.92
case(11)   !    11.) Salt (49% LiF, 29% NaF, 29% ZrF4)
xlo=709.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (49% LiF, 29% NaF, 29% ZrF4)")
Conductivity = 0.53
case(12)   !    12.) Salt (58% KF, 42% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% KF, 42% ZrF4)")
Conductivity = 0.45
case(13)   !    13.) Salt (58% LiCl, 42% RbCl)
xlo=586.; xhi=1323.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% LiCl, 42% RbCl)")
Conductivity = 0.39
case(14)   !    14.) Salt (58% NaCl, 42% MgCl2)
xlo=718.; xhi=1738.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (58% NaCl, 42% MgCl2)")
Conductivity = 0.43
case(15)   !    15.) Salt (59.5% LiCl, 40.5% KCl)
xlo=628.; xhi=1673.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% LiCl, 40.5% KCl)")
Conductivity = 0.43
case(16)   !    16.) Salt (59.5% NaF, 40.5% ZrF4)
xlo=773.; xhi=1623.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (59.5% NaF, 40.5% ZrF4)")
Conductivity = 0.49
case(17)   !    17.) Salt (60% NaNO3, 40% KNO3)
xlo=533.; xhi=866.;
if(T.lt.xlo.or.T.gt.xhi) dum=t_warn(T,xlo,xhi,"Salt (60% NaNO3, 40% KNO3)")
Conductivity = -1E-11*T**3 + 3E-08*T**2 + 0.0002*T + 0.3922
case(18:25) !  18-25.) Trough properties
Conductivity = k_fluid((T-273.15),int(fnumd))
case(26)    ! Argon
Conductivity = 0.00548 + 0.0000438969*T - 6.81410E-09*T**2
case(27)   !Hydrogen
Conductivity = dmax1(0.0302888 + 0.00053634*T - 1.59604E-07*T**2,.01d0)
case(28:35)
continue  !no information
case(36:) !Any integer greater than 35
!Call the user-defined property table
lb=fl_bounds(fnum-35)
ub=fl_bounds(fnum-35+1)-1
if(ub.lt.lb) ub=size(fprop(1,:))
dxx(:)=fprop(1,lb:ub)
dyy(:)=fprop(6,lb:ub)
call interp(T,size(dxx),dxx,dyy,Gjsav,Conductivity)
if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
end select


end function


!Density of Fluid [kg/m3]
Double PRecision Function dens_fluid(T, fluidnum) 
implicit none
Double Precision T, dens_salt, dens_caloria, dens_salt_xl, dens_therminol,&
                 dens_salt_hitec, dens_Dowtherm_Q, dens_Dowtherm_RP, dens_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); dens_fluid = 1; !no props for this model
    case(18); dens_fluid = dens_salt(T) ;
    case(19); dens_fluid = dens_caloria(T) ;
    case(20); dens_fluid = dens_salt_xl(T) ;
    case(21); dens_fluid = dens_therminol(T) ;
    case(22); dens_fluid = dens_salt_hitec(T) ;
    case(23); dens_fluid = dens_Dowtherm_Q(T) ;
    case(24); dens_fluid = dens_Dowtherm_RP(T) ;
    case(25); dens_fluid = dens_salt_xl(T) ;
    case(26:35); dens_fluid = 1; !no props for this model
    case(36:); dens_fluid = dens_user(T,Fluidnum-35) ;
    end select
    
End Function

!Conductivity of Fluid [W/mK]
Double PRecision Function k_fluid(T, fluidnum)
implicit none
Double Precision T, k_salt, k_caloria, k_salt_xl, k_therminol,&
                 k_salt_hitec, k_Dowtherm_Q, k_Dowtherm_RP, k_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); k_fluid = 1; !no props for this model
    case(18); k_fluid = k_salt(T); 
    case(19); k_fluid = k_caloria(T); 
    case(20); k_fluid = k_salt_xl(T); 
    case(21); k_fluid = k_therminol(T);
    case(22); k_fluid = k_salt_hitec(T);
    case(23); k_fluid = k_Dowtherm_Q(T);
    case(24); k_fluid = k_Dowtherm_RP(T);
    case(25); k_fluid = k_salt_xl(T);
    case(26:35); k_fluid = 1; !no props for this model
    case(36:); k_fluid = k_user(T,Fluidnum-35) ;
    end select
    
End Function

!Kinematic Viscosity of Fluid [m2/s]
Double PRecision Function kin_visc_fluid(T, fluidnum)
implicit none
Double Precision T, kin_visc_salt, kin_visc_caloria, kin_visc_salt_xl, kin_visc_therminol,&
                 kin_visc_salt_hitec, kin_visc_Dowtherm_Q, kin_visc_Dowtherm_RP, kin_visc_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); kin_visc_fluid = 1; !no props for this model
    case(18); kin_visc_fluid = kin_visc_salt(T); 
    case(19); kin_visc_fluid = kin_visc_caloria(T); 
    case(20); kin_visc_fluid = kin_visc_salt_xl(T); 
    case(21); kin_visc_fluid = kin_visc_therminol(T);
    case(22); kin_visc_fluid = kin_visc_salt_hitec(T);
    case(23); kin_visc_fluid = kin_visc_Dowtherm_Q(T);
    case(24); kin_visc_fluid = kin_visc_Dowtherm_RP(T);
    case(25); kin_visc_fluid = kin_visc_salt_xl(T);
    case(26:35); kin_visc_fluid = 1; !no props for this model
    case(36:); kin_visc_fluid = kin_visc_user(T,Fluidnum-35) ;
    end select
    
End Function

!Absolute Viscosity of Fluid [kg-s/m]
Double PRecision Function visc_fluid(T, fluidnum) 
implicit none
Double Precision T, visc_salt, visc_caloria, visc_salt_xl, visc_therminol,&
                 visc_salt_hitec, visc_Dowtherm_Q, visc_Dowtherm_RP, visc_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); visc_fluid = 1; !no props for this model
    case(18); visc_fluid = visc_salt(T); 
    case(19); visc_fluid = visc_caloria(T); 
    case(20); visc_fluid = visc_salt_xl(T); 
    case(21); visc_fluid = visc_therminol(T);
    case(22); visc_fluid = visc_salt_hitec(T);
    case(23); visc_fluid = visc_Dowtherm_Q(T);
    case(24); visc_fluid = visc_Dowtherm_RP(T);
    case(25); visc_fluid = visc_salt_xl(T);
    case(26:35); visc_fluid = 1; !no props for this model
    case(36:); visc_fluid = visc_user(T,Fluidnum-35) ;
    end select

End Function

!Specific Heat of Fluid [J/kgC]
Double PRecision Function Cp_fluid(T, fluidnum)
implicit none
Double Precision T, Cp_salt, Cp_caloria, Cp_salt_xl, Cp_therminol,&
                 Cp_salt_hitec, Cp_Dowtherm_Q, Cp_Dowtherm_RP, Cp_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); Cp_fluid = 1; !no props for this model
    case(18); Cp_fluid = Cp_salt(T); 
    case(19); Cp_fluid = Cp_caloria(T); 
    case(20); Cp_fluid = Cp_salt_xl(T); 
    case(21); Cp_fluid = Cp_therminol(T);
    case(22); Cp_fluid = Cp_salt_hitec(T);
    case(23); Cp_fluid = Cp_Dowtherm_Q(T);
    case(24); Cp_fluid = Cp_Dowtherm_RP(T);
    case(25); Cp_fluid = Cp_salt_xl(T);
    case(26:35); Cp_fluid = 1; !no props for this model
    case(36:); Cp_fluid = Cp_user(T,Fluidnum-35) ;
    end select

End Function

!Enthalpy of Fluid [J/kg]
Double PRecision Function H_fluid(T, fluidnum)
implicit none
Double Precision T, H_salt, H_caloria, H_salt_xl, H_therminol,&
                 H_salt_hitec, H_Dowtherm_Q, H_Dowtherm_RP, H_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); H_fluid = 1; !no props for this model
    case(18); H_fluid = H_salt(T); 
    case(19); H_fluid = H_caloria(T); 
    case(20); H_fluid = H_salt_xl(T); 
    case(21); H_fluid = H_therminol(T);
    case(22); H_fluid = H_salt_hitec(T);
    case(23); H_fluid = H_Dowtherm_Q(T);
    case(24); H_fluid = H_Dowtherm_RP(T);
    case(25); H_fluid = H_salt_xl(T);
    case(26:35); H_fluid = 1; !no props for this model
    case(36:); H_fluid = H_user(T,Fluidnum-35) ;
    end select

End Function

!Temperature of Fluid [C]
Double PRecision Function T_fluid(H, fluidnum) 
implicit none
Double Precision H, T_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); T_fluid = 1; !no props for this model
    case(18)
        T_fluid = -0.0000000000262 * H**2 + 0.0006923 * H + 0.03058
    case(19)
        T_fluid = 6.4394E-17 * H**3 - 0.00000000023383 * H**2 + 0.0005821 * H + 1.2744
    case(20)
        T_fluid = 0.00000000005111 * H**2 + 0.0006466 * H + 0.2151
    case(21)
        T_fluid = 7.4333E-17 * H**3 - 0.00000000024625 * H**2 + 0.00063282 * H + 12.403
    case(22)
        T_fluid = -3.309E-24 * H**2 + 0.000641 * H + 0.000000000001364
    case(23)
        T_fluid = 6.186E-17 * H**3 - 0.00000000022211 * H**2 + 0.00059998 * H + 0.77742
    case(24)
        T_fluid = 6.6607E-17 * H**3 - 0.00000000023347 * H**2 + 0.00061419 * H + 0.77419
    case(25)
        T_fluid = 0.00000000005111 * H**2 + 0.0006466 * H + 0.2151
    case(26:35); T_fluid = 1; !no props for this model
    case(36:); T_fluid = T_user(H,Fluidnum-35) ;
    End select
    
End Function

!Prandtl Number of Fluid []
Double PRecision Function Pr_fluid(T, fluidnum) 
implicit none
Double Precision T, Pr_salt, Pr_caloria, Pr_salt_xl, Pr_therminol,&
                 Pr_salt_hitec, Pr_Dowtherm_Q, Pr_Dowtherm_RP, Pr_user
Integer Fluidnum

    select case(Fluidnum)
    case(1:17); Pr_fluid = 1; !no props for this model
    case(18); Pr_fluid = Pr_salt(T); 
    case(19); Pr_fluid = Pr_caloria(T); 
    case(20); Pr_fluid = Pr_salt_xl(T); 
    case(21); Pr_fluid = Pr_therminol(T);
    case(22); Pr_fluid = Pr_salt_hitec(T);
    case(23); Pr_fluid = Pr_Dowtherm_Q(T);
    case(24); Pr_fluid = Pr_Dowtherm_RP(T);
    case(25); Pr_fluid = Pr_salt_xl(T);
    case(26:35); Pr_fluid = 1; !no props for this model
    case(36:); Pr_fluid = Pr_user(T,Fluidnum-35) ;
    end select

End Function

!
! **************** Caloria HT 43 *****************************
!
!Density of Caloria HT 43 [kg/m3]
Double Precision Function dens_caloria(T) !T[C]
implicit none
Double Precision T
dens_caloria = dmax1(885 - 0.6617 * T - 0.0001265 * T**2,100.d0)
End Function

!Conductivity of Caloria HT 43 [W/mK]
Double Precision Function k_caloria(T) !T[C]
implicit none
Double Precision T
k_caloria = dmax1(-0.00014 * T + 0.1245,.01d0)
End Function

!Kinematic Viscosity of Caloria HT 43 [m2/s]
Double Precision Function kin_visc_caloria(T) !T[C]
implicit none
Double Precision T
kin_visc_caloria = 0.040439268 * T**-1.946401872
End Function

!Absolute Viscosity of Caloria HT 43 [m2/s]
Double Precision Function visc_caloria(T) !T[C]
implicit none
Double Precision T, dens_caloria
visc_caloria = (0.040439268 * T**-1.946401872) * dens_caloria(T)
End Function

!Specific Heat of Caloria HT 43 [J/kgC]
Double Precision Function Cp_caloria(T) !T[C]
implicit none
Double Precision T
Cp_caloria = 3.88 * T + 1606.0
End Function

!Enthalpy of Caloria HT 43 [J/kg]
Double Precision Function H_caloria(T) !T[C]
implicit none
Double Precision T
H_caloria = 1.94 * T**2 + 1606.0 * T
End Function

!Prandtl Number of Caloria HT 43 []
Double Precision Function Pr_caloria(T) !T[C]
implicit none
Double Precision T, kin_visc_caloria, dens_caloria, Cp_caloria, k_caloria
Pr_caloria = kin_visc_caloria(T) * dens_caloria(T) * Cp_caloria(T) / k_caloria(T)
End Function

!
! ********************* Nitrate Salt ***********************
!
!Density of Nitrate Salt, [kg/m3]
Double Precision Function dens_salt(T) !T [C]
implicit none
Double Precision T
dens_salt = dmax1(2090 - 0.636 * T,1000.d0)
End Function

!Heat Capacity of Nitrate Salt, [J/kg/K]
Double Precision Function Cp_salt(T) !T [C]
implicit none
Double Precision T
Cp_salt = 1443 + 0.172 * T
End Function

!Enthalpy of Salt (English Units), [Btu/lbm]
Double Precision Function H_engl_salt(T) !T [F]
implicit none
Double Precision T
H_engl_salt = 0.345 * T + 0.0000114 * T**2
End Function

!Enthalpy of Salt (Metric Units), [J/kg]
Double Precision Function H_salt(T)  !T[C]
implicit none
Double Precision T
H_salt = 1443 * T + 0.086 * T**2
End Function

!Nondimensional Penetration Distance for Nitrate Salt, z/D
Double Precision Function pen_salt(V, Tsalt, D, Ttube) !V=salt velocity [m/s], Tsalt=temperature of salt [C], D = tube inside diameter [m], Ttube = temperature of tube [C]
implicit none
Double Precision V,Tsalt,D,Ttube, Pr_salt, Re_salt, diff_salt, Cp_salt
pen_salt = 0.23 * Pr_salt(Tsalt)**0.5 * Re_salt(V, D, Tsalt)**0.75 * (diff_salt(Tsalt) / 0.000000151)**(1 / 9) * (161000 / (1481 * (221 - Ttube)))**(1 / 3) * (1 + 0.7 * Cp_salt(Tsalt) * (Tsalt - 221) / 161000)
End Function

!Thermal Conductivity of Nitrate Salt, W/m/K
Double Precision Function k_salt(T) !T [C]
implicit none
Double Precision T
k_salt = 0.443 + 0.00019 * T
End Function

!Absolute Viscosity of Nitrate Salt, [Pa s]
Double Precision Function visc_salt(T) !T [C]
implicit none
Double Precision T
visc_salt = dmax1((22.714 - 0.12 * T + 0.0002281 * T * T - 0.0000001474 * T * T * T) / 1000,1.e-6)
End Function

!Kinematic Viscosity of Nitrate Salt, [m2/s]
Double Precision Function kin_visc_salt(T) !T [C]
implicit none
Double Precision T, visc_salt, dens_salt
kin_visc_salt = visc_salt(T) / dens_salt(T)
End Function

!Thermal Diffusivity of Nitrate Salt, [m2/s]
Double Precision Function diff_salt(T) !T [C]
implicit none
Double Precision T, dens_salt, Cp_salt, k_salt
diff_salt = k_salt(T) / (dens_salt(T) * Cp_salt(T))
End Function

!Prandtl Number for Nitrate Salt
Double Precision Function Pr_salt(T) !T [C]
implicit none
Double Precision T, visc_salt, dens_salt, diff_salt
Pr_salt = visc_salt(T) / (dens_salt(T) * diff_salt(T))
End Function

!Reynolds number for Nitrate Salt
Double Precision Function Re_salt(V, D, T) ! V=salt velocity [m/s], D=inside diameter [m], T=salt temperature [C]
implicit none
Double Precision V,D,T, visc_salt, dens_salt
Re_salt = V * D / (visc_salt(T) / dens_salt(T))
End Function

!Volumetric Thermal Expansion Coefficient for Nitrate Salt, [1/C]
Double Precision Function beta_salt(T) !T [C]
implicit none
Double Precision T, dens_salt
beta_salt = 0.636 / dens_salt(T)
End Function
!
! ********************** HITEC XL ***********************
!
!Density of HITEC XL Nitrate Salt, [kg/m3]
Double Precision Function dens_salt_xl(T) !T [C]
implicit none
Double Precision T
dens_salt_xl = dmax1(2240 - 0.8266 * T,800.d0)
End Function

!Heat Capacity of HITEC XL Nitrate Salt, [J/kg/K]
Double Precision Function Cp_salt_xl(T) !T [C]
implicit none
Double Precision T
Cp_salt_xl = dmax1(1536 - 0.2624 * T - 0.0001139 * T * T,1000.d0)
End Function

!Enthalpy of HITEC XL Nitrate Salt, [J/kg]
Double Precision Function H_salt_xl(T) !T [C]
implicit none
Double Precision T
H_salt_xl = 1536 * T - 0.1312 * T**2 - 0.0000379667 * T**3
End Function

!Thermal Conductivity of HITEC XL Nitrate Salt, W/m/K
Double Precision Function k_salt_xl(T) !T [C]
implicit none
Double Precision T
k_salt_xl = 0.519
End Function

!Absolute Viscosity of HITEC XL Nitrate Salt, [Pa s]
Double Precision Function visc_salt_xl(T) !T [C]
implicit none
Double Precision T
visc_salt_xl = 1372000 * T**-3.364
End Function

!Kinematic Viscosity of HITEC XL Nitrate Salt, [m2/s]
Double Precision Function kin_visc_salt_xl(T) !T [C]
implicit none
Double Precision T, visc_salt_xl, dens_salt_xl
kin_visc_salt_xl = visc_salt_xl(T) / dens_salt_xl(T)
End Function

!Thermal Diffusivity of HITEC XL Nitrate Salt, [m2/s]
Double Precision Function diff_salt_xl(T) !T [C]
implicit none
Double Precision T, k_salt_xl, dens_salt_xl, Cp_salt_xl
diff_salt_xl = k_salt_xl(T) / (dens_salt_xl(T) * Cp_salt_xl(T))
End Function

!Prandtl Number for HITEC XL Nitrate Salt
Double Precision Function Pr_salt_xl(T) !T [C]
implicit none
Double Precision T, visc_salt_xl, dens_salt_xl, diff_salt_xl
Pr_salt_xl = visc_salt_xl(T) / (dens_salt_xl(T) * diff_salt_xl(T))
End Function

!Reynolds number for HITEC XL Nitrate Salt
Double Precision Function Re_salt_xl(V, D, T) ! V=salt velocity [m/s], D=inside diameter [m], T=salt temperature [C]
implicit none
Double Precision T, V, D, visc_salt_xl, dens_salt_xl
Re_salt_xl = V * D / (visc_salt_xl(T) / dens_salt_xl(T))
End Function

!
! ********************************* HITEC Salt **************************
!
!Density of HITEC Salt, [kg/m3]
Double Precision Function dens_salt_hitec(T) !T [C]
implicit none
Double Precision T
dens_salt_hitec = dmax1(2080 - 0.733 * T,1000.d0)
End Function

!Heat Capacity of HITEC Salt, [J/kg/K]
Double Precision Function Cp_salt_hitec(T) !T [C]
implicit none
Double Precision T
Cp_salt_hitec = 1560 - 0.0 * T
End Function

!Enthalpy of HITEC Salt, [J/kg]
Double Precision Function H_salt_hitec(T) !T [C]
implicit none
Double Precision T
H_salt_hitec = 1560 * T
End Function

!Thermal Conductivity of HITEC Salt, W/m/K
Double Precision Function k_salt_hitec(T) !T [C]
implicit none
Double Precision T
k_salt_hitec = 0.588 - 0.000647 * T
End Function

!Absolute Viscosity of HITEC Salt, [Pa s]
Double Precision Function visc_salt_hitec(T) !T [C] !Poor fit
implicit none
Double Precision T
visc_salt_hitec = dmax1(0.00622 - 0.0000102 * T,1.e-6)
End Function

!Kinematic Viscosity of HITEC Salt, [m2/s]
Double Precision Function kin_visc_salt_hitec(T) !T [C]
implicit none
Double Precision T, visc_salt_hitec, dens_salt_hitec
kin_visc_salt_hitec = visc_salt_hitec(T) / dens_salt_hitec(T)
End Function

!Thermal Diffusivity of HITEC Salt, [m2/s]
Double Precision Function diff_salt_hitec(T) !T [C]
implicit none
Double Precision T, k_salt_hitec, dens_salt_hitec, Cp_salt_hitec
diff_salt_hitec = k_salt_hitec(T) / (dens_salt_hitec(T) * Cp_salt_hitec(T))
End Function

!Prandtl Number for HITEC Salt
Double Precision Function Pr_salt_hitec(T) !T [C]
implicit none
Double Precision T, visc_salt_hitec, dens_salt_hitec, diff_salt_hitec
Pr_salt_hitec = visc_salt_hitec(T) / (dens_salt_hitec(T) * diff_salt_hitec(T))
End Function

!Reynolds number for HITEC Salt
Double Precision Function Re_salt_hitec(V, D, T) ! V=salt velocity [m/s], D=inside diameter [m], T=salt temperature [C]
implicit none
Double Precision V,D,T, visc_salt_hitec, dens_salt_hitec
Re_salt_hitec = V * D / (visc_salt_hitec(T) / dens_salt_hitec(T))
End Function

!
! ********************* Therminol VP-1 **************************
!
!Density of Therminol Oil [kg/m3]
Double Precision Function dens_therminol(T) !T [C]
implicit none
Double Precision T
dens_therminol = dmax1(1074.0 - 0.6367 * T - 0.0007762 * T**2,400.d0)
End Function

!Specific Heat of Therminol Oil, J/kg/K
Double Precision Function Cp_therminol(T) !T [C]
implicit none
Double Precision T
Cp_therminol = 1000 * (1.509 + 0.002496 * T + 0.0000007888 * T**2)
End Function

!Absoute Viscosity of Therminol Oil [Pa s]
Double Precision Function visc_therminol(T) !T [C]
implicit none
Double Precision T
T = dmax1(T,25.) !limit range of T
visc_therminol = 0.001 * (10**0.8703 * T**(0.2877 + Log10(T**-0.3638)))
End Function

!Thermal conductivity of Therminol Oil [W/mK]
Double Precision Function k_therminol(T) !T [C]
implicit none
Double Precision T
k_therminol = dmax1(0.1381 - 0.00008708 * T - 0.0000001729 * T**2,.001d0)
End Function

!Enthalpy of Therminol Oil [J/kg]
Double Precision Function H_therminol(T) !T [C]
implicit none
Double Precision T
H_therminol = 1000 * (-18.34 + 1.498 * T + 0.001377 * T**2)
End Function

!Temperature of Therminol Oil [C]
Double Precision Function T_therminol(H) !H [J/kg]
implicit none
Double Precision H
T_therminol = -0.000000000158 * H**2 + 0.0006072 * H + 13.37
End Function

!Vapor Pressure of Therminol Oil [kPa]
Double Precision Function Pvap_therminol(T) !T [C]
implicit none
Double Precision T
Pvap_therminol = 10**-11.13 * T**5.448
End Function

!Kinematic Viscosity of Therminol Oil, [m2/s]
Double Precision Function kin_visc_therminol(T) !T [C]
implicit none
Double Precision T, visc_therminol, dens_therminol
kin_visc_therminol = visc_therminol(T) / dens_therminol(T)
End Function

!Thermal Diffusivity of Therminol Oil, [m2/s]
Double Precision Function diff_therminol(T) !T [C]
implicit none
Double Precision T, k_therminol, dens_therminol, Cp_therminol
diff_therminol = k_therminol(T) / (dens_therminol(T) * Cp_therminol(T))
End Function

!Prandtl Number for Therminol Oil, []
Double Precision Function Pr_therminol(T) !T [C]
implicit none
Double Precision T, visc_Therminol, dens_therminol, diff_therminol
Pr_therminol = visc_therminol(T) / (dens_therminol(T) * diff_therminol(T))
End Function

!Reynolds number for Therminol Oil, []
Double Precision Function Re_therminol(V, D, T) ! V=salt velocity [m/s], D=inside diameter [m], T=salt temperature [C]
implicit none
Double Precision V,D,T, visc_therminol, dens_therminol
Re_therminol = V * D / (visc_therminol(T) / dens_therminol(T))
End Function
!
! ******************** Dowtherm Q ****************************
!
!Density of Dowtherm Q [kg/m3]
Double Precision Function dens_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T
dens_Dowtherm_Q = dmax1(-0.757332 * T + 980.787,100.d0)                               ! Russ 10-2-03
End Function

!Specific Heat of Dowtherm Q, J/kg/K
Double Precision Function Cp_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T
Cp_Dowtherm_Q = -0.00053943 * T**2 + 3.2028 * T + 1589.2               ! Russ 10-2-03
End Function

!Absoute Viscosity of Dowtherm Q [Pa s]
Double Precision Function visc_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T
visc_Dowtherm_Q = 1 / (132.40658 + 4.36107 * T + 0.0781417 * T**2 - 0.00011035416 * T**3) !Hank 10-2-03
End Function

!Thermal conductivity of Dowtherm Q [W/mK]
Double Precision Function k_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T
k_Dowtherm_Q = dmax1(-0.0000000626555 * T**2 - 0.000124864 * T + 0.124379,1.e-5)    !Russ 10-2-03
End Function

!Enthalpy of Dowtherm Q [J/kg]
Double Precision Function H_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T
H_Dowtherm_Q = (0.00151461 * T**2 + 1.59867 * T - 0.0250596) * 1000    !Hank 10-2-03
End Function

!Temperature of Dowtherm Q [C]
Double Precision Function T_Dowtherm_Q(H) !H [J/kg]
implicit none
Double Precision H
T_Dowtherm_Q = -0.0000000001136 * H**2 + 0.000552 * H + 4.318    !Hank 10-8-03
End Function

!Vapor Pressure of Dowtherm Q [W/mK]
Double Precision Function Pvap_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T
Pvap_Dowtherm_Q = -6238.9915 + 0.000000037388014 * T**5.5254072     !Hank 10-2-03
End Function

!Kinematic Viscosity of Dowtherm Q, [m2/s]
Double Precision Function kin_visc_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T, visc_Dowtherm_Q, dens_Dowtherm_Q
kin_visc_Dowtherm_Q = visc_Dowtherm_Q(T) / dens_Dowtherm_Q(T)
End Function

!Thermal Diffusivity of Dowtherm Q, [m2/s]
Double Precision Function diff_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T, k_dowtherm_Q, dens_Dowtherm_Q, Cp_Dowtherm_Q
diff_Dowtherm_Q = k_Dowtherm_Q(T) / (dens_Dowtherm_Q(T) * Cp_Dowtherm_Q(T))
End Function

!Prandtl Number for Dowtherm Q, []
Double Precision Function Pr_Dowtherm_Q(T) !T [C]
implicit none
Double Precision T, visc_Dowtherm_Q, dens_Dowtherm_Q, diff_Dowtherm_Q
Pr_Dowtherm_Q = visc_Dowtherm_Q(T) / (dens_Dowtherm_Q(T) * diff_Dowtherm_Q(T))
End Function

!Reynolds number for Dowtherm Q, []
Double Precision Function Re_Dowtherm_Q(V, D, T) ! V=salt velocity [m/s], D=inside diameter [m], T=salt temperature [C]
implicit none
Double Precision V,D,T, visc_Dowtherm_Q, dens_Dowtherm_Q
Re_Dowtherm_Q = V * D / (visc_Dowtherm_Q(T) / dens_Dowtherm_Q(T))
End Function
!
!*************** Dowtherm RP **************************
!
!Density of Dowtherm RP [kg/m3]
Double Precision Function dens_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T
dens_Dowtherm_RP = dmax1(-0.000186495 * T**2 - 0.668337 * T + 1042.11,200.d0)  !Russ 10-2-03
End Function

!Specific Heat of Dowtherm RP, J/kg/K
Double Precision Function Cp_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T
Cp_Dowtherm_RP = -0.0000031915 * T**2 + 2.977 * T + 1560.8       !Russ 10-2-03
End Function

!Absoute Viscosity of Dowtherm RP [Pa s]
Double Precision Function visc_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T
visc_Dowtherm_RP = 1 / (4.523003 + 0.39156855 * T + 0.028604206 * T**2)  !Hank 10-2-03
End Function

!Thermal conductivity of Dowtherm RP [W/mK]
Double Precision Function k_Dowtherm_RP(T) !T [C]
Double Precision T
k_Dowtherm_RP = -0.00012963 * T + 0.13397                         !Russ 10-2-03
End Function

!Enthalpy of Dowtherm RP [J/kg]
Double Precision Function H_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T
H_Dowtherm_RP = (0.0014879 * T**2 + 1.5609 * T - 0.0024798) * 1000      !Hank 10-2-03
End Function

!Temperature of Dowtherm RP [C]
Double Precision Function T_Dowtherm_RP(H) !H [J/kg]
implicit none
Double Precision H
T_Dowtherm_RP = -0.0000000001192 * H**2 + 0.0005647 * H + 4.343 !Hank 10-8-03
End Function

!Vapor Pressure of Dowtherm RP [Pa]
Double Precision Function Pvap_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T
Pvap_Dowtherm_RP = -2863.8678 + 2.7661827E-12 * T**6.894446     !Hank 10-2-03                       !Hank 10-2-03
End Function

!Kinematic Viscosity of Dowtherm RP, [m2/s]
Double Precision Function kin_visc_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T, visc_Dowtherm_RP, dens_Dowtherm_RP
kin_visc_Dowtherm_RP = visc_Dowtherm_RP(T) / dens_Dowtherm_RP(T)
End Function

!Thermal Diffusivity of Dowtherm RP, [m2/s]
Double Precision Function diff_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T, k_Dowtherm_RP, dens_Dowtherm_RP, Cp_Dowtherm_RP
diff_Dowtherm_RP = k_Dowtherm_RP(T) / (dens_Dowtherm_RP(T) * Cp_Dowtherm_RP(T))
End Function

!Prandtl Number for Dowtherm RP, []
Double Precision Function Pr_Dowtherm_RP(T) !T [C]
implicit none
Double Precision T, visc_Dowtherm_RP, dens_Dowtherm_RP, diff_Dowtherm_RP
Pr_Dowtherm_RP = visc_Dowtherm_RP(T) / (dens_Dowtherm_RP(T) * diff_Dowtherm_RP(T))
End Function

!Reynolds number for Dowtherm RP, []
Double Precision Function Re_Dowtherm_RP(V, D, T) ! V=salt velocity [m/s], D=inside diameter [m], T=salt temperature [C]
implicit none
Double Precision V,D,T, visc_Dowtherm_RP, dens_Dowtherm_RP
Re_Dowtherm_RP = V * D / (visc_Dowtherm_RP(T) / dens_Dowtherm_RP(T))
End Function



!************************* User Fluid ****************************
!| Note that the fluid properties are stored in the FPROP array  *
!| in the following order:                                       *
!*****************************************************************
!|  #    |   1   |   2   |   3   |   4   |   5   |   6   |   7   |
!|-------|-------|-------|-------|-------|-------|-------|-------|
!| Name  |   T   |   Cp  |  rho  |   Mu  |   Nu  |   k   |   h   |
!| Units |   C   |kJ/kg-K| kg/m3 |  Pa-s |  m2-s | W/m-K |  J/kg |
!|-------|-------|-------|-------|-------|-------|-------|-------|
!|   1   |   :   |   :   |   :   |   :   |   :   |   :   |   :   |
!|   2   |   :   |   :   |   :   |   :   |   :   |   :   |   :   |
!*****************************************************************
!Density of User Fluid, [kg/m3]
Double Precision Function dens_user(T,fn) !T [C]
use global_props	
implicit none
Double Precision:: T
integer,intent(in)::fn
double precision,dimension(size(fprop(1,fl_bounds(fn):(fl_bounds(fn+1)-1))))::dxx,dyy !Create dummy arrays
integer::lb,ub,t_warn,dum
    !Call the user-defined property table
    lb=fl_bounds(fn)
    ub=fl_bounds(fn+1)-1
    if(ub.lt.lb) ub=size(fprop(1,:))
    dxx(:)=fprop(1,lb:ub)
    dyy(:)=fprop(3,lb:ub)
    call interp(T,size(dxx),dxx,dyy,Gjsav,dens_user)
    if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
End Function

!Heat Capacity of User Fluid, [J/kg/K]
Double Precision Function Cp_user(T,fn) !T [C]
use global_props	
implicit none
Double Precision T
integer,intent(in)::fn
double precision,dimension(size(fprop(1,fl_bounds(fn):(fl_bounds(fn+1)-1))))::dxx,dyy !Create dummy arrays
integer::lb,ub,t_warn,dum
    !Call the user-defined property table
    lb=fl_bounds(fn)
    ub=fl_bounds(fn+1)-1
    if(ub.lt.lb) ub=size(fprop(1,:))
    dxx(:)=fprop(1,lb:ub)
    dyy(:)=fprop(2,lb:ub)
    call interp(T,size(dxx),dxx,dyy,Gjsav,Cp_user)
    !Convert from kJ/kg-K to J/kg-K
    Cp_user = Cp_user*1000.
    if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
End Function

!Enthalpy of User Fluid, [J/kg]
Double Precision Function H_user(T,fn) !T [C]
use global_props	
implicit none
Double Precision T
integer,intent(in)::fn
double precision,dimension(size(fprop(1,fl_bounds(fn):(fl_bounds(fn+1)-1))))::dxx,dyy !Create dummy arrays
integer::lb,ub,t_warn,dum
    !Call the user-defined property table
    lb=fl_bounds(fn)
    ub=fl_bounds(fn+1)-1
    if(ub.lt.lb) ub=size(fprop(1,:))
    dxx(:)=fprop(1,lb:ub)
    dyy(:)=fprop(7,lb:ub)
    call interp(T,size(dxx),dxx,dyy,Gjsav,H_user)
    if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
End Function

!Temperature of User fluid, [C]
Double Precision Function T_user(H,fn) !H [J/kg-K]
use global_props	
implicit none
Double Precision H
integer,intent(in)::fn
double precision,dimension(size(fprop(1,fl_bounds(fn):(fl_bounds(fn+1)-1))))::dxx,dyy !Create dummy arrays
integer::lb,ub,t_warn,dum
    !Call the user-defined property table
    lb=fl_bounds(fn)
    ub=fl_bounds(fn+1)-1
    if(ub.lt.lb) ub=size(fprop(1,:))
    dxx(:)=fprop(7,lb:ub)
    dyy(:)=fprop(1,lb:ub)
    call interp(H,size(dxx),dxx,dyy,Gjsav,T_user)
    if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(H,dxx(lb),dxx(ub),"User-specified fluid")
end function

!Thermal Conductivity of User Fluid, W/m/K
Double Precision Function k_user(T,fn) !T [C]
use global_props	
implicit none
Double Precision T
integer,intent(in)::fn
double precision,dimension(size(fprop(1,fl_bounds(fn):(fl_bounds(fn+1)-1))))::dxx,dyy !Create dummy arrays
integer::lb,ub,t_warn,dum
    !Call the user-defined property table
    lb=fl_bounds(fn)
    ub=fl_bounds(fn+1)-1
    if(ub.lt.lb) ub=size(fprop(1,:))
    dxx(:)=fprop(1,lb:ub)
    dyy(:)=fprop(6,lb:ub)
    call interp(T,size(dxx),dxx,dyy,Gjsav,k_user)
    if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
End Function

!Absolute Viscosity of User Fluid, [Pa s]
Double Precision Function visc_user(T,fn) !T [C] !Poor fit
use global_props	
implicit none
Double Precision T
integer,intent(in)::fn
double precision,dimension(size(fprop(1,fl_bounds(fn):(fl_bounds(fn+1)-1))))::dxx,dyy !Create dummy arrays
integer::lb,ub,t_warn,dum
    !Call the user-defined property table
    lb=fl_bounds(fn)
    ub=fl_bounds(fn+1)-1
    if(ub.lt.lb) ub=size(fprop(1,:))
    dxx(:)=fprop(1,lb:ub)
    dyy(:)=fprop(4,lb:ub)
    call interp(T,size(dxx),dxx,dyy,Gjsav,visc_user)
    if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
End Function

!Kinematic Viscosity of User Fluid, [m2/s]
Double Precision Function kin_visc_user(T,fn) !T [C]
use global_props	
implicit none
Double Precision T
integer,intent(in)::fn
double precision,dimension(size(fprop(1,fl_bounds(fn):(fl_bounds(fn+1)-1))))::dxx,dyy !Create dummy arrays
integer::lb,ub,t_warn,dum
    !Call the user-defined property table
    lb=fl_bounds(fn)
    ub=fl_bounds(fn+1)-1
    if(ub.lt.lb) ub=size(fprop(1,:))
    dxx(:)=fprop(1,lb:ub)
    dyy(:)=fprop(5,lb:ub)
    call interp(T,size(dxx),dxx,dyy,Gjsav,kin_visc_user)
    if((Gjsav.eq.ub).or.(Gjsav.eq.lb)) dum=t_warn(T,dxx(lb),dxx(ub),"User-specified fluid")
End Function

!Thermal Diffusivity of User Fluid, [m2/s]
Double Precision Function diff_user(T,fn) !T [C]
implicit none
Double Precision T, k_user, dens_user, Cp_user
integer::fn
    diff_user = k_user(T,fn) / (dens_user(T,fn) * Cp_user(T,fn))
End Function

!Prandtl Number for User Fluid
Double Precision Function Pr_user(T,fn) !T [C]
implicit none
Double Precision T, visc_user, dens_user, diff_user
integer::fn
    Pr_user = visc_user(T,fn) / (dens_user(T,fn) * diff_user(T,fn))
End Function

!*******************************************************************************