SUBROUTINE TYPE224 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) 
!************************************************************************
! Object: Solar Power Tower Rankine Component
! Simulation Studio Model: TYPE224
 
! Author: Michael J. Wagner
! Editor: -
! Date:	 July 03, 2008 
! Last modified: March 19th, 2010
! COPYRIGHT 2010 NATIONAL RENEWABLE ENERGY LABORATORY

! Parameters
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Input Units   |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  P_ref                         |  Reference output electric power at design condition            |  MW            |  kW
!  2|  eta_ref                       |  Reference conversion efficiency at design condition            |  none          |  none
!  3|  T_htf_hot_ref                 |  Reference HTF inlet temperature at design                      |  C             |  K
!  4|  T_htf_cold_ref                |  Reference HTF outlet temperature at design                     |  C             |  K
!  5|  dT_cw_ref                     |  Reference condenser cooling water inlet/outlet T diff          |  C             |  C
!  6|  HTF                           |  Integer flag identifying HTF in power block                    |  none          |  none
!  7|  q_sby_frac                    |  Fraction of thermal power required for standby mode            |  none          |  none
!  8|  P_boil                        |  Boiler operating pressure                                      |  bar           |  bar
!  9|  CT_array(1)                   |  Flag for using dry cooling or wet cooling system               |  none          |  none
! 10|  startup_time                  |  Time needed for power block startup                            |  hr            |  hr
! 11|  startup_frac                  |  Fraction of design thermal power needed for startup            |  none          |  none

! Inputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Input Units   |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  mode                          |  Cycle part load control, from plant controller                 |  none          |  none
!  2|  T_htf_hot                     |  Hot HTF inlet temperature, from storage tank                   |  C             |  K
!  3|  m_dot_htf                     |  HTF mass flow rate                                             |  kg/hr         |  kg/hr
!  4|  T_wb                          |  Ambient wet bulb temperature                                   |  C             |  C
!  5|  demand_var                    |  Control signal indicating operational mode                     |  none          |  none
!  6|  standby_control               |  Control signal indicating standby mode                         |  none          |  none
!  7|  T_db                          |  Ambient dry bulb temperature                                   |  C             |  C
!  8|  P_amb                         |  Ambient pressure                                               |  atm           |  Pa

! Outputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Output Units  |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  P_cycle                       |  Cycle power output                                             |  MWe           |  kWe
!  2|  eta                           |  Cycle thermal efficiency                                       |  none          |  none
!  3|  T_htf_cold                    |  Heat transfer fluid outlet temperature                         |  C             |  C
!  4|  m_dot_makeup                  |  Cooling water makeup flow rate                                 |  kg/hr         |  kg/s
!  5|  m_dot_demand                  |  HTF required flow rate to meet power load                      |  kg/hr         |  kg/hr
!  6|  m_dot_htf                     |  Actual HTF flow rate passing through the power cycle           |  kg/hr         |  kg/hr
!  7|  m_dot_htf_ref                 |  Calculated reference HTF flow rate at design                   |  kg/hr         |  kg/hr
!  8|  W_cool_par                    |  Cooling system parasitic load                                  |  MWe           |  MWe
!  9|  P_ref                         |  Reference power level output at design (mirror param)          |  MWe           |  kWe

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

!    TRNSYS acess functions (allow to acess TIME etc.) 
USE TrnsysConstants
USE TrnsysFunctions
use global_props  !user global module

!REQUIRED BY THE MULTI-DLL VERSION OF TRNSYS
!DEC$ATTRIBUTES DLLEXPORT :: TYPE224				!SET THE CORRECT TYPE NUMBER HERE

!-----------------------------------------------------------------------------------------------------------------------
implicit none !real(8) (a-z) 

!     Parameters
real(8):: P_ref, eta_ref,  T_htf_hot_ref, HTF, T_htf_cold_ref, dT_cw_ref, q_sby_frac, P_boil,Is_dry_bulb, &
          startup_time, startup_frac, T_amb_des, T_approach, T_ITD_des, P_cond_ratio

!    Inputs
real(8):: mode, T_htf_hot, T_wb, T_db, m_dot_cw, demand_var, m_dot_htf,standby_control, last_P_cycle, P_amb

!    Outputs / Local variables
real(8):: P_cycle, eta, T_htf_cold, m_dot_demand,m_dot_htf_ref, q_ND_tot, c_p_w, &
          q_sby_needed, m_dot_sby, c_htf, q_tot, specheat, T_cw_in, last_standby_control, tstep, startup_remain ,&
          tech_type, CT_array(3), W_cool_par, m_dot_makeup


real(8):: xin, out, time, par, stored, T, dTdt, prop(7), pb_bd_frac, h_st_hot, h_st_cold, dh_steam, m_dot_st_bd !,CT_array(3)
integer(4)::INFO(15), np, ni, nout, nd, npar, nin, nder, iunit, itype, icntrl, nstored, ierr
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
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!SET THE MAXIMUM NUMBER OF PARAMETERS (NP), INPUTS (NI),OUTPUTS (NOUT), AND DERIVATIVES (ND) 
PARAMETER (NP=17,NI=8,NOUT=9,ND=0,NSTORED=3)

!REQUIRED TRNSYS DIMENSIONS
DIMENSION XIN(NI),OUT(NOUT),PAR(NP),YCHECK(NI),OCHECK(NOUT),STORED(NSTORED),T(ND),DTDT(ND)
INTEGER NITEMS
logical called
!-----------------------------------------------------------------------------------------------------------------------

tstep=getSimulationTimeStep()

!-----------------------------------------------------------------------------------------------------------------------
!READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
P_ref = PAR(1)*1000. !Convert from MW to kW
eta_ref = PAR(2)
T_htf_hot_ref = PAR(3)
T_htf_cold_ref = PAR(4)
dT_cw_ref = PAR(5)
T_amb_des = PAR(6)
HTF=PAR(7)
q_sby_frac=PAR(8)
P_boil=PAR(9)
CT_array(1)=PAR(10) !Cooling type
startup_time=PAR(11)
startup_frac=PAR(12)
tech_type=PAR(13)
T_approach=PAR(14)
T_ITD_des=PAR(15)
P_cond_ratio=PAR(16)
pb_bd_frac=PAR(17)

!Set up the condenser info array 
if(CT_array(1)==1.) then
    CT_array(2)=dT_cw_ref
    CT_array(3)=T_approach
else
    CT_array(2)=T_ITD_des
    CT_array(3)=P_cond_ratio
endif


!-----------------------------------------------------------------------------------------------------------------------
!RETRIEVE THE CURRENT VALUES OF THE INPUTS TO THIS MODEL FROM THE XIN ARRAY IN SEQUENTIAL ORDER
mode = XIN(1)
T_htf_hot = XIN(2)
m_dot_htf = XIN(3)
T_wb = XIN(4) + 273.15      !Convert to K
demand_var = XIN(5)
standby_control = XIN(6)
T_db = XIN(7) + 273.15      !Convert to K
P_amb = xin(8)*101300. ![atm] -> [Pa]
    IUNIT=INFO(1)
    ITYPE=INFO(2)
    if(mode.eq.1.) demand_var = demand_var*1000.  !If the mode is to operate in power demand, convert from MW to kW

!-----------------------------------------------------------------------------------------------------------------------
!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(db)) deallocate(db)  !deallocate the power cycle coef. array
    RETURN 1
ENDIF
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!'AFTER-ITERATION' MANIPULATIONS 
IF (INFO(13).GT.0) THEN
    STORED(1)=standby_control
    STORED(2)=max(startup_remain-tstep,0.)
    STORED(3)=P_cycle

    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
    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)

    !RETURN TO THE CALLING PROGRAM
    RETURN 1

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

!DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE 
IF (TIME .LT. (getSimulationStartTime() + getSimulationTimeStep()/2.D0)) THEN
    !SET THE UNIT NUMBER FOR FUTURE CALLS
    IUNIT=INFO(1)
    ITYPE=INFO(2)

    !PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE
    !Cycle power output
    out(:)=0.d0
    !Heat transfer fluid outlet temp
    OUT(3)=T_htf_hot
    
    !Calculate the power block side steam enthalpy rise for blowdown calculations
    !Steam properties are as follows:
    !=======================================================================
    ! | T(C) | P(MPa) | h(kJ/kg) | s(kJ/kg-K) | x(-) | v(m3/kg) | U(kJ/kg) |
    !=======================================================================
    !hot steam
    prop = (/(T_htf_hot_ref - 25.d0),(P_boil*0.1d0),1.d0,1.d0,1.d0,1.d0,1.d0/)
    call steam_props("SI",prop,12,ierr)
    h_st_hot = prop(3) !Use the hot steam enthalpy
    prop = (/1.d0,(P_boil*0.1d0),1.d0,1.d0,0.d0,1.d0,1.d0/)
    call steam_props("SI",prop,25,ierr)
    !Use the cold steam enthalpy at x=0, subtract subcooled enthalpy with specific heat. Cp is based
    !on an integral of specific heat for water in the slightly subcooled range of ~100degC
    h_st_cold = prop(3) - 4.91*100. 
    dh_steam = (h_st_hot - h_st_cold)

    RETURN 1

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

!RETRIEVE THE VALUES IN THE STORAGE ARRAY FOR THIS ITERATION
CALL getStorageVars(STORED,NSTORED,INFO)
last_standby_control=STORED(1)
startup_remain=STORED(2)
last_P_cycle=STORED(3)

!-----------------------------------------------------------------------------------------------------------------------
!    *** PERFORM ALL THE CALCULATION HERE FOR THIS MODEL. ***
!-----------------------------------------------------------------------------------------------------------------------

!Specific heat of water
c_p_w=4.183  ![kJ/kg-C]

      
select case (int(standby_control))
case(1)  !The cycle is in normal operation

    call RankineCycle(time,P_ref, eta_ref, T_htf_hot_ref, T_htf_cold_ref, T_db, T_wb, P_amb, dT_cw_ref, HTF, c_p_w, &
                      T_htf_hot, m_dot_htf, mode, demand_var,P_boil, tech_type,T_amb_des, CT_array,&
                      P_cycle, eta, T_htf_cold, m_dot_demand, m_dot_htf_ref, m_dot_makeup, W_cool_par)

    !Check the output to make sure it's reasonable. If not, return zeros.
    if(((eta>1.).or.(eta<0.)).or.((T_htf_cold > T_htf_hot).or.(T_htf_cold < T_htf_cold_ref - 50.))) then
        P_cycle = 0.d0
        eta = 0.d0
        T_htf_cold = T_htf_cold_ref
    endif
    
    !-----Calculate the blowdown fraction-----
    m_dot_st_bd = P_cycle/(eta*dh_steam)*pb_bd_frac

case(2)  !The cycle is in standby operation
    c_htf=specheat(HTF,(T_htf_hot+T_htf_cold_ref)/2.+273.15)
    q_tot = P_ref*(1.-eta_ref)

    !Calculate the actual q_sby_needed from the reference flows
    q_sby_needed = q_tot*q_sby_frac

    !now calculate the mass flow rate knowing the inlet temperature of the salt,
    !..and holding the outlet temperature at the reference outlet temperature
    m_dot_sby = q_sby_needed/(c_htf*(T_htf_hot - T_htf_cold_ref))*3600.

    
    !Set other output values
    P_cycle = 0.
    eta = 0.
    T_htf_cold = T_htf_cold_ref
    m_dot_demand = m_dot_sby
    m_dot_makeup = 0.
    W_cool_par = 0.

case(3)  !The cycle has been completely shut down
    P_cycle = 0.
    eta = 0.
    T_htf_cold = T_htf_cold_ref  !Changed from T_htf_hot 12/18/2009 was causing problems with T250
    m_dot_demand = 0.0
    m_dot_makeup = 0.
    W_cool_par = 0.
    
end select
      
      
!If the cycle is going from completely shut down to starting up, set the remaining startup
!time to be equal to the designated startup time
if((last_standby_control.eq.3.).and.(standby_control.eq.1.)) then
    startup_remain=startup_time
endif

!If the cycle is starting up beginning in this time period, or it is continuing to start
!up from the last time period, then subtract the appropriate fraction of electric power
!from the output.  Note that during the startup time, not only is the cycle not producing power,
!but it is also consuming thermal energy
if(((last_standby_control.eq.3.).and.(standby_control.eq.1.)).or.(startup_remain.gt.0.)) then
    P_cycle = P_cycle*(1.-min(1.,startup_remain/tstep))          !The power cycle only operates after startup

    !Fraction of the timestep running at full capacity
    m_dot_demand = m_dot_demand*(1.-min(1.,startup_remain/tstep))+ m_dot_demand*min(startup_remain/tstep,1.)*startup_frac 

    eta = eta_ref
    T_htf_cold = T_htf_cold_ref
    m_dot_htf_ref = m_dot_htf_ref
    
    startup_remain=max(startup_remain-tstep,0.)
endif
      
      
!SET THE OUTPUTS FROM THIS MODEL

!Cycle power output
OUT(1)=P_cycle/1000.  !Convert from kW to MW
!Cycle thermal efficiency
OUT(2)=eta
!Heat transfer fluid outlet temp
OUT(3)=T_htf_cold
!Wet cooling makeup water flow rate
OUT(4)=(m_dot_makeup + m_dot_st_bd)*3600.
!Heat transfer fluid demand flow rate
OUT(5)=m_dot_demand
!Heat transfer fluid flow rate
OUT(6)=m_dot_htf			
!Calculated reference htf flow rate			
OUT(7)=m_dot_htf_ref
!Cooling tower parasitic load [MW]
OUT(8)=W_cool_par
!Reference power level output
OUT(9)=P_ref/1000.   !Convert from kW to MW

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


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

subroutine RankineCycle(time,P_ref, eta_ref, T_htf_hot_ref, T_htf_cold_ref, T_db, T_wb, P_amb, dT_cw_ref, HTF, c_p_w, &
                        T_htf_hot, m_dot_htf, mode, demand_var,P_boil, tech_type,T_amb_des,CT_array, &
                        P_cycle, eta, T_htf_cold, m_dot_demand, m_dot_htf_ref, m_dot_makeup,W_cool_par)

!implicit none
implicit none !real(8) (a-z)

real(8):: P_ref, eta_ref, T_htf_hot_ref, T_htf_cold_ref, T_cw_in, dT_cw_ref, c_htf, HTF, c_htf_ref, &
          T_htf_hot, m_dot_htf, m_dot_cw, mode, demand_var, P_boil, tech_type,T_amb_des,q_dot_ref, &
          c_p_w, T_htf_hot_ND, m_dot_htf_ND, m_dot_cw_ND,adj, err, P_dem_ND, P_ND_tot, Q_ND_tot
real(8):: P_ND(3), Q_ND(3), T_ref, specheat,Interpolate, P_AB, P_AC, P_BC, Q_AB, Q_AC, Q_BC, CT_array(3),&
          W_cool_par, P_cond, T_cond, q_reject_est, m_dot_makeup, q_reject, m_dot_air, T_ITD_des, P_cond_ratio,&
          T_approach, P_amb, T_db, T_wb, time

!Declare inputs
!real(8):: P_ND(3), Q_ND(3), CT_array(3)
integer::i, qq, TT, CT
logical:: float_pres
!Declare output variables
real(8),intent(out):: P_cycle, eta, T_htf_cold, m_dot_demand, m_dot_htf_ref

!______________________________________________________________________________

!Select the technology type
!	tech_type=1. !TEMP
if(tech_type.eq.1.) TT=1  !high temp
if(tech_type.eq.2.) TT=2  !low temp

!Adjust for the temperature design condition.  We need to ensure that the design cycle efficiency is
!acheieved at the specified temperature. The cycle adjustment for condenser pressure is equal to 1.0 
!at roughly 8600 [Pa], which corresponds to a condenser temperature of about 43degC. 
if(CT_array(1)==1) then
    !Wet cooled case
    T_wb = T_wb + (43. - (CT_array(2) + 3. + CT_array(3))) - T_amb_des
    T_db = T_db + (43. - (CT_array(2) + 3. + CT_array(3))) - T_amb_des
else
    !Dry cooled case
    T_db = T_db + (43. - CT_array(2)) - T_amb_des
    !Wet bulb isn't used in this situation
endif

!Calculate the specific heat before converting to Kelvin
c_htf_ref = specheat(HTF,(T_htf_hot_ref+T_htf_cold_ref)/2.+273.15)
c_htf = specheat(HTF,(T_htf_hot+T_htf_cold_ref)/2.+273.15)

!Convert units
!**Temperatures from Celcius to Kelvin
T_htf_hot = T_htf_hot + 273.15  ![K]
T_htf_hot_ref = T_htf_hot_ref + 273.15 ![K]
T_htf_cold_ref = T_htf_cold_ref + 273.15 ![K]
!Mass flow rates from kg/hr to kg/s
m_dot_htf = m_dot_htf/3600. ![kg/s]

!****Calculate the reference values
q_dot_ref = P_ref/eta_ref   !The reference heat flow
m_dot_htf_ref = q_dot_ref/(c_htf_ref*(T_htf_hot_ref - T_htf_cold_ref))  !The HTF mass flow rate [kg/s]

T_ref = T_sat(P_boil) !The saturation temp at the boiler

!Calculate the htf hot temperature, in non-dimensional form
T_htf_hot_ND = (T_htf_hot - T_ref)/(T_htf_hot_ref - T_ref)

!float_pres = .false.
!if(float_pres) then
!    T_ref = T_htf_hot - (T_htf_hot_ref - T_ref)  !Make the new reference temperature track the inlet HTF temperature
!    
!    P_boil= dmax1(P_sat(T_ref),60.)  !The saturation pressure at the new reference temperature, limit to 60 bar
!    T_ref_ll = T_sat(P_boil)  !Recalculate to make sure the reference temp hasn't fallen below the lower limit (ll)
!    
!    !Recalculate the hot reference temperature
!    T_htf_hot_ref = T_htf_hot + (T_ref - T_ref_ll)
!    
!    T_htf_hot_ND = (T_htf_hot - dmax1(T_ref_ll,T_ref))/(T_htf_hot_ref - dmax1(T_ref_ll,T_ref))
!    continue
!endif
!****Determine the non-dimensional parameters needed for the correlations

!Calculate the htf mass flow rate in non-dimensional form
m_dot_htf_ND = m_dot_htf/m_dot_htf_ref

!Do an initial cooling tower call to estimate the turbine back pressure. 
q_reject_est = q_dot_ref*1000.*(1.-eta_ref)*m_dot_htf_ND*T_htf_hot_ND
!Sort out the cooling tower info array
CT = int(CT_array(1))  !Cooling technology type {1=evaporative cooling, 2=air cooling}
if(CT==1) then
    !For a wet-cooled system
    dT_cw_ref = CT_array(2)       ![C] Cooling water temp rise at design
    T_approach = CT_array(3)      ![C] Approach temperature
    call evap_tower(dT_cw_ref, T_approach, (P_ref*1000.), eta_ref, T_db, T_wb, P_amb, q_reject_est, m_dot_makeup, W_cool_par, P_cond, T_cond)
else
    !For a dry-cooled system
    T_ITD_des = CT_array(2)       ![C] ACC initial temperature difference at design
    P_cond_ratio = CT_array(3)    ![-] Condenser inlet/outlet pressure ratio at design
    call ACC(T_ITD_des, P_cond_ratio, (P_ref*1000.), eta_ref, T_db, P_amb, q_reject_est, m_dot_air, W_cool_par, P_cond, T_cond)
    m_dot_makeup = 0.d0
endif



!  Set initial values
ADJ = 1. ; err = 1. ; qq=0

!Do a quick check to see if there is actually a mass flow being supplied 
!  to the cycle. If not, go to the end.
if(abs(m_dot_htf_ND).lt.1.e-3) then
    P_cycle = 0.
    eta = 0.
    T_htf_cold = T_htf_hot_ref
    m_dot_demand = m_dot_htf_ref
    W_cool_par = 0.
    m_dot_makeup = 0.
    !Set the error to zero, since we don't want to iterate
    err=0.
endif

!Begin iterations
do while ((err.gt.1.e-6).and.(qq.lt.100)) 
    qq=qq+1
    !Now use the constrained variable to calculate the demand mass flow rate
    if(mode.eq.1.) then
        P_dem_ND = demand_var/P_ref
        if(qq.eq.1) m_dot_htf_ND = P_dem_ND   !An initial guess (function of power)
        if(qq.gt.1) m_dot_htf_ND = m_dot_htf_ND*ADJ
    elseif(mode.eq.2.) then
        continue     ! do nothing
    endif

    !++++++++++++++Correlations++++++++++++++++++
    !Calculate the correlations
    !++++++++++++++++++++++++++++++++++++++++++++
    !POWER  
    !Main effects
    P_ND(1) = Interpolate(TT,11,1,T_htf_hot_ND)-1.
    P_ND(2) = Interpolate(TT,12,2,P_cond)-1.
    P_ND(3) = Interpolate(TT,13,3,m_dot_htf_ND)-1.

    !Interactions
    P_AB = Interpolate(TT,112,12,P_cond)-1.
    P_AC = Interpolate(TT,113,13,m_dot_htf_ND)-1.  !*(adtx*T_htf_hot_ND + (1-adtx))-1
    P_BC = Interpolate(TT,123,23,m_dot_htf_ND)-1.  !*(acmx*m_dot_htf_ND+(1-acmx))*(actx*T_htf_hot_ND+(1-actx))-1.

    P_ND(1) = P_ND(1)*(1.+P_AB+P_AC)
    P_ND(2) = P_ND(2)*(1.+P_BC)
    P_ND(3) = P_ND(3)   !*(1.+Interpolate(TT,134,34,m_dot_htf_ND)-1.)

    !HEAT               
    !Main effects
    Q_ND(1) = Interpolate(TT,21,1,T_htf_hot_ND)-1.
    Q_ND(2) = Interpolate(TT,22,2,P_cond)-1.
    Q_ND(3) = Interpolate(TT,23,3,m_dot_htf_ND)-1.

    !Interactions
    Q_AB = Interpolate(TT,212,12,P_cond)-1.
    Q_AC = Interpolate(TT,213,13,m_dot_htf_ND)-1.  
    Q_BC = Interpolate(TT,223,23,m_dot_htf_ND)-1.  

    Q_ND(1) = Q_ND(1)*(1.+Q_AB+Q_AC)
    Q_ND(2) = Q_ND(2)*(1.+Q_BC)
    Q_ND(3) = Q_ND(3)   

    !Calculate the cumulative values
    P_ND_tot = 1.
    Q_ND_tot = 1.

    !Add up main effects, which are normalized about 1.. so subtract 1 from each effect
    do i=1,size(P_ND)
        P_ND_tot = P_ND_tot + P_ND(i)
        Q_ND_tot = Q_ND_tot + Q_ND(i)
    enddo
    !Add up interactions, which have not been normalized about 1, and add directly

    !Calculate the output values:
    P_cycle = P_ND_tot*P_ref
    T_htf_cold = T_htf_hot-Q_ND_tot*q_dot_ref/(m_dot_htf*c_htf)
    eta = P_cycle/(Q_ND_tot*q_dot_ref)
    m_dot_demand = dmax1(m_dot_htf_ND*m_dot_htf_ref,0.00001d0)   ![kg/s]
    
    !Call the cooling tower model to update the condenser pressure
    q_reject = (1.-eta)*q_dot_ref*Q_ND_tot*1000.       !P_cycle*1000.*(1./eta - 1.)
    if(CT==1) then
        call evap_tower(dT_cw_ref, T_approach, (P_ref*1000.), eta_ref, T_db, T_wb, P_amb, q_reject, m_dot_makeup, W_cool_par, P_cond, T_cond)
    else
        call ACC(T_ITD_des, P_cond_ratio, (P_ref*1000.), eta_ref, T_db, P_amb, q_reject, m_dot_air, W_cool_par, P_cond, T_cond)
    endif
    
    
    !Check to see if the calculated and demand values match
    !If they don't match, calculate the "ADJ" factor
    if(mode.eq.1.) then
        err = (P_cycle - demand_var)/demand_var
        ADJ = 1.+(demand_var-P_cycle)/(3.*demand_var)
    elseif(mode.eq.2.) then
        err = 0.
    endif

    if(qq.eq.99) then
        print *,'Did not converge after 100 iterations'
        P_cycle = 0.
        eta = 0.
        T_htf_cold = T_htf_hot_ref
        m_dot_demand = m_dot_htf_ref
        exit
    endif
    !If this is not true, the cycle has not yet converged, and we should return
    ! to continue in the iterations
enddo

!Finally, convert the values back to their original units
T_htf_cold = T_htf_cold - 273.15 ![K]-->[C]
T_htf_cold_ref = T_htf_cold_ref - 273.15 ![K]->[C]
T_htf_hot_ref = T_htf_hot_ref - 273.15 ![K]->[C]
m_dot_demand = m_dot_demand*3600.  ![kg/s]->[kg/hr]
m_dot_htf = m_dot_htf*3600.        ![kg/s]->[kg/hr]
m_dot_htf_ref = m_dot_htf_ref*3600. ![kg/s]->[kg/hr]

contains

    real(8) function T_sat(P)
        implicit none !real(8) (a-z)
        real(8)::P
        !P is in [bar], T is in [K]
        !T_sat = 450.214 + 2.2691*P - 0.0123064*P**2 + 0.0000302829*P**3  !The saturation temp at the boiler
        T_sat = 439.486188 + 2.88644991*P - 0.0243308856*P**2 + 0.000125910226*P**3 - 2.66233987E-07*P**4
    end function

    real(8) function P_sat(T)
        implicit none !real(8) (a-z)
        real(8)::T
        !T is in [K], P is in [bar]
        !P_sat = -1551.35 + 9.94705*T - 0.0218711*T**2 + 0.0000165747*T**3
        P_sat = 1965.19859 - 15.1689701*T + 0.0452997046*T**2 - 0.000063150801*T**3 + 3.54340123E-08*T**4
    end function

end subroutine
	
	
!Interpolation subroutine
double precision function Interpolate(TT,YT,XT,X)
use global_props
implicit none

!This function interpolates the data of one data list based on the value provided to a 
!..corresponding list of the same length. 
!TT: The flag indicating whether the High temperature or Low temperature tables should be used
!YT: The name of the dependent Y variable being interpolated
!XT: The name of the independent X variable which the Y variable is a function of
!X: The value of the X variable
!Interpolate: is returned as a double precision value of Y

integer,intent(in)::YT,XT,TT
real(8),intent(in)::X
real(8)::Y,ind
!double precision,allocatable::db(:,:)  !The database of interpolation values. Generally of the form A, PA, QA, B, PB, QB...
real(8),dimension(30)::datx,daty
integer::XI,YI,i,j,lbi,ubi,jsav,n

if(.not.allocated(db)) then
    ! On first call, allocate the array and define coefficients
    allocate(db(18,30))
if(TT==2)then !Low temperature parabolic trough applications
	db(1,:)=(/ &
		0.10000, 0.14830, 0.19660, 0.24480, 0.29310, 0.34140,  &
		0.38970, 0.43790, 0.48620, 0.53450, 0.58280, 0.63100,  &
		0.67930, 0.72760, 0.77590, 0.82410, 0.87240, 0.92070,  &
		0.96900, 1.01720, 1.06550, 1.11380, 1.16210, 1.21030,  &
		1.25860, 1.30690, 1.35520, 1.40340, 1.45170, 1.50000  &
		/)
	db(2,:)=(/ &
		0.08427, 0.12766, 0.17257, 0.21883, 0.26630, 0.31525,  &
		0.36420, 0.41430, 0.46494, 0.51597, 0.56721, 0.61853,  &
		0.66966, 0.72078, 0.77145, 0.82169, 0.87141, 0.92056,  &
		0.96901, 1.01714, 1.06579, 1.11448, 1.16321, 1.21199,  &
		1.26084, 1.30977, 1.35878, 1.40788, 1.45709, 1.50641  &
		/)
	db(3,:)=(/ &
		0.09945, 0.14752, 0.19569, 0.24395, 0.29234, 0.34088,  &
		0.38942, 0.43806, 0.48674, 0.53541, 0.58405, 0.63261,  &
		0.68102, 0.72942, 0.77760, 0.82567, 0.87361, 0.92143,  &
		0.96915, 1.01677, 1.06431, 1.11178, 1.15919, 1.20655,  &
		1.25389, 1.30120, 1.34850, 1.39580, 1.44310, 1.49043  &
		/)
	db(4,:)=(/ &
		2800.0000, 3566.0000, 4331.0000, 5097.0000, 5862.0000, 6628.0000,  &
		7393.0000, 8159.0000, 8924.0000, 9690.0000, 10455.0000, 11221.0000,  &
		11986.0000, 12752.0000, 13517.0000, 14283.0000, 15048.0000, 15814.0000,  &
		16579.0000, 17345.0000, 18110.0000, 18876.0000, 19641.0000, 20407.0000,  &
		21172.0000, 21938.0000, 22703.0000, 23469.0000, 24234.0000, 25000.0000  &
		/)
	db(5,:)=(/ &
		1.09549, 1.07656, 1.06098, 1.04772, 1.03613, 1.02583,  &
		1.01654, 1.00823, 1.00045, 0.99324, 0.98653, 0.98024,  &
		0.97433, 0.96874, 0.96344, 0.95841, 0.95360, 0.94901,  &
		0.94462, 0.94041, 0.93635, 0.93244, 0.92868, 0.92504,  &
		0.92151, 0.91810, 0.91480, 0.91159, 0.90847, 0.90544  &
		/)
	db(6,:)=(/ &
		1.01320, 1.01059, 1.00843, 1.00658, 1.00495, 1.00349,  &
		1.00217, 1.00098, 0.99987, 0.99883, 0.99786, 0.99695,  &
		0.99609, 0.99528, 0.99451, 0.99377, 0.99306, 0.99239,  &
		0.99174, 0.99112, 0.99052, 0.98994, 0.98938, 0.98884,  &
		0.98831, 0.98780, 0.98731, 0.98683, 0.98636, 0.98590  &
		/)
	db(7,:)=(/ &
		0.10000, 0.14830, 0.19660, 0.24480, 0.29310, 0.34140,  &
		0.38970, 0.43790, 0.48620, 0.53450, 0.58280, 0.63100,  &
		0.67930, 0.72760, 0.77590, 0.82410, 0.87240, 0.92070,  &
		0.96900, 1.01720, 1.06550, 1.11380, 1.16210, 1.21030,  &
		1.25860, 1.30690, 1.35520, 1.40340, 1.45170, 1.50000  &
		/)
	db(8,:)=(/ &
		0.09100, 0.13673, 0.18349, 0.23106, 0.27952, 0.32862,  &
		0.37824, 0.42817, 0.47846, 0.52892, 0.57940, 0.62970,  &
		0.67989, 0.72977, 0.77922, 0.82803, 0.87628, 0.92376,  &
		0.97037, 1.01636, 1.06286, 1.10926, 1.15557, 1.20168,  &
		1.24780, 1.29382, 1.33976, 1.38552, 1.43129, 1.47697  &
		/)
	db(9,:)=(/ &
		0.10522, 0.15521, 0.20482, 0.25402, 0.30304, 0.35181,  &
		0.40036, 0.44860, 0.49675, 0.54472, 0.59253, 0.64009,  &
		0.68759, 0.73496, 0.78219, 0.82920, 0.87619, 0.92306,  &
		0.96982, 1.01637, 1.06292, 1.10936, 1.15571, 1.20188,  &
		1.24803, 1.29411, 1.34010, 1.38591, 1.43174, 1.47748  &
		/)
	db(10,:)=(/ &
		2500.0000, 3684.0000, 4868.0000, 6053.0000, 7237.0000, 8421.0000,  &
		9605.0000, 10789.5000, 11974.0000, 13158.0000, 14342.0000, 15526.0000,  &
		16711.0000, 17895.0000, 19079.0000, 20263.0000, 21447.0000, 22632.0000,  &
		23816.0000, 25000.0000, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(11,:)=(/ &
		1.09243, 1.06413, 1.04308, 1.02623, 1.01208, 1.00000,  &
		0.98924, 0.97955, 0.97084, 0.96283, 0.95543, 0.94856,  &
		0.94214, 0.93609, 0.93039, 0.92501, 0.91990, 0.91502,  &
		0.91037, 0.90592, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(12,:)=(/ &
		1.01085, 1.00757, 1.00511, 1.00312, 1.00144, 1.00000,  &
		0.99871, 0.99754, 0.99649, 0.99552, 0.99461, 0.99378,  &
		0.99299, 0.99225, 0.99155, 0.99088, 0.99025, 0.98964,  &
		0.98907, 0.98851, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(13,:)=(/ &
		0.20000, 0.26840, 0.33680, 0.40530, 0.47370, 0.54210,  &
		0.61050, 0.67890, 0.74740, 0.81580, 0.88420, 0.95260,  &
		1.02110, 1.08950, 1.15790, 1.22630, 1.29475, 1.36320,  &
		1.43160, 1.50000, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(14,:)=(/ &
		0.19165, 0.26272, 0.33519, 0.40833, 0.48106, 0.55267,  &
		0.62238, 0.68939, 0.75301, 0.81228, 0.87086, 0.93417,  &
		1.00000, 1.06387, 1.12749, 1.19089, 1.25405, 1.31712,  &
		1.37988, 1.44244, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(15,:)=(/ &
		0.20521, 0.27354, 0.34131, 0.40870, 0.47559, 0.54211,  &
		0.60830, 0.67419, 0.73989, 0.80524, 0.87035, 0.93524,  &
		1.00000, 1.06447, 1.12874, 1.19283, 1.25675, 1.32059,  &
		1.38418, 1.44761, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(16,:)=(/ &
		0.20000, 0.26840, 0.33680, 0.40530, 0.47370, 0.54210,  &
		0.61050, 0.67890, 0.74740, 0.81580, 0.88420, 0.95260,  &
		1.02110, 1.08950, 1.15790, 1.22630, 1.29470, 1.36320,  &
		1.43160, 1.50000, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(17,:)=(/ &
		0.18667, 0.25349, 0.32139, 0.39017, 0.45935, 0.52876,  &
		0.59814, 0.66723, 0.73590, 0.80369, 0.87047, 0.93600,  &
		1.00000, 1.06432, 1.12843, 1.19234, 1.25607, 1.31969,  &
		1.38304, 1.44623, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(18,:)=(/ &
		0.21169, 0.28118, 0.34959, 0.41730, 0.48419, 0.55040,  &
		0.61606, 0.68110, 0.74579, 0.80995, 0.87364, 0.93697,  &
		1.00000, 1.06258, 1.12480, 1.18674, 1.24836, 1.30983,  &
		1.37083, 1.43167, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)

else  !Power tower applications
    ! On first call, allocate the array and define coefficients
	db(1,:)=(/ &
		0.10000, 0.14830, 0.19660, 0.24480, 0.29310, 0.34140,  &
		0.38970, 0.43790, 0.48620, 0.53450, 0.58280, 0.63100,  &
		0.67930, 0.72760, 0.77590, 0.82410, 0.87240, 0.92070,  &
		0.96900, 1.01720, 1.06550, 1.11380, 1.16210, 1.21030,  &
		1.25860, 1.30690, 1.35520, 1.40340, 1.45170, 1.50000  &
		/)
	db(2,:)=(/ &
		0.07890, 0.12063, 0.16439, 0.20998, 0.25720, 0.30584,  &
		0.35562, 0.40628, 0.45752, 0.50910, 0.56081, 0.61249,  &
		0.66417, 0.71567, 0.76693, 0.81763, 0.86795, 0.91785,  &
		0.96751, 1.01781, 1.06827, 1.11892, 1.16975, 1.22077,  &
		1.27199, 1.32341, 1.37502, 1.42677, 1.47900, 1.53140  &
		/)
	db(3,:)=(/ &
		0.09714, 0.14487, 0.19310, 0.24178, 0.29085, 0.34018,  &
		0.38966, 0.43915, 0.48852, 0.53771, 0.58666, 0.63533,  &
		0.68386, 0.73218, 0.78028, 0.82794, 0.87538, 0.92262,  &
		0.96967, 1.01656, 1.06330, 1.10993, 1.15645, 1.20288,  &
		1.24923, 1.29551, 1.34173, 1.38792, 1.43430, 1.48066  &
		/)
	db(4,:)=(/ &
		2800.0000, 3565.5172, 4331.0345, 5096.5517, 5862.0690, 6627.5862,  &
		7393.1034, 8158.6207, 8924.1379, 9689.6552, 10455.1724, 11220.6897,  &
		11986.2069, 12751.7241, 13517.2414, 14282.7586, 15048.2759, 15813.7931,  &
		16579.3103, 17344.8276, 18110.3448, 18875.8621, 19641.3793, 20406.8966,  &
		21172.4138, 21937.9310, 22703.4483, 23468.9655, 24234.4828, 25000.0000  &
		/)
	db(5,:)=(/ &
		1.07879, 1.06191, 1.04803, 1.03621, 1.02588, 1.01670,  &
		1.00842, 1.00100, 0.99406, 0.98764, 0.98165, 0.97605,  &
		0.97077, 0.96578, 0.96106, 0.95656, 0.95228, 0.94818,  &
		0.94426, 0.94050, 0.93688, 0.93339, 0.93003, 0.92678,  &
		0.92364, 0.92059, 0.91764, 0.91478, 0.91199, 0.90929  &
		/)
	db(6,:)=(/ &
		1.00936, 1.00736, 1.00569, 1.00427, 1.00302, 1.00190,  &
		1.00088, 0.99997, 0.99912, 0.99832, 0.99758, 0.99688,  &
		0.99622, 0.99560, 0.99501, 0.99445, 0.99391, 0.99339,  &
		0.99289, 0.99242, 0.99196, 0.99151, 0.99109, 0.99067,  &
		0.99027, 0.98988, 0.98950, 0.98913, 0.98878, 0.98843  &
		/)
	db(7,:)=(/ &
		0.10000, 0.14830, 0.19660, 0.24480, 0.29310, 0.34140,  &
		0.38970, 0.43790, 0.48620, 0.53450, 0.58280, 0.63100,  &
		0.67930, 0.72760, 0.77590, 0.82410, 0.87240, 0.92070,  &
		0.96900, 1.01720, 1.06550, 1.11380, 1.16210, 1.21030,  &
		1.25860, 1.30690, 1.35520, 1.40340, 1.45170, 1.50000  &
		/)
	db(8,:)=(/ &
		0.09403, 0.14058, 0.18795, 0.23606, 0.28481, 0.33408,  &
		0.38376, 0.43374, 0.48389, 0.53410, 0.58426, 0.63425,  &
		0.68395, 0.73327, 0.78207, 0.83025, 0.87771, 0.92432,  &
		0.97009, 1.01632, 1.06244, 1.10844, 1.15434, 1.20013,  &
		1.24581, 1.29140, 1.33689, 1.38229, 1.42760, 1.47282  &
		/)
	db(9,:)=(/ &
		0.10659, 0.15680, 0.20654, 0.25589, 0.30491, 0.35365,  &
		0.40213, 0.45039, 0.49843, 0.54628, 0.59394, 0.64145,  &
		0.68879, 0.73598, 0.78303, 0.82995, 0.87674, 0.92341,  &
		0.96996, 1.01640, 1.06273, 1.10896, 1.15508, 1.20111,  &
		1.24705, 1.29290, 1.33865, 1.38433, 1.42992, 1.47542  &
		/)
	db(10,:)=(/ &
		1750.0000, 3078.9474, 4407.8947, 5736.8421, 7065.7895, 8394.7368,  &
		9723.6842, 11052.6316, 12381.5789, 13710.5263, 15039.4737, 16368.4211,  &
		17697.3684, 19026.3158, 20355.2632, 21684.2105, 23013.1579, 24342.1053,  &
		25671.0526, 27000.0000, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(11,:)=(/ &
		1.09728, 1.06400, 1.04185, 1.02504, 1.01140, 1.00000,  &
		0.99000, 0.98116, 0.97321, 0.96598, 0.95935, 0.95322,  &
		0.94751, 0.94194, 0.93689, 0.93211, 0.92757, 0.92323,  &
		0.91909, 0.91511, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(12,:)=(/ &
		1.00857, 1.00568, 1.00373, 1.00224, 1.00102, 1.00000,  &
		0.99910, 0.99830, 0.99758, 0.99692, 0.99631, 0.99575,  &
		0.99522, 0.99473, 0.99427, 0.99383, 0.99341, 0.99302,  &
		0.99264, 0.99228, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(13,:)=(/ &
		0.20000, 0.26840, 0.33680, 0.40530, 0.47370, 0.54210,  &
		0.61050, 0.67890, 0.74740, 0.81580, 0.88420, 0.95260,  &
		1.02110, 1.08950, 1.15790, 1.22630, 1.29470, 1.36320,  &
		1.43160, 1.50000, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(14,:)=(/ &
		0.19414, 0.26485, 0.33672, 0.40891, 0.48073, 0.55147,  &
		0.62040, 0.68683, 0.75003, 0.80931, 0.87036, 0.93545,  &
		1.00000, 1.06378, 1.12731, 1.19061, 1.25370, 1.31657,  &
		1.37925, 1.44173, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(15,:)=(/ &
		0.20308, 0.27115, 0.33879, 0.40606, 0.47301, 0.53968,  &
		0.60608, 0.67225, 0.73819, 0.80392, 0.86946, 0.93482,  &
		1.00000, 1.06502, 1.12987, 1.19458, 1.25914, 1.32357,  &
		1.38786, 1.45202, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(16,:)=(/ &
		0.10000, 0.17370, 0.24740, 0.32110, 0.39470, 0.46840,  &
		0.54210, 0.61580, 0.68950, 0.76320, 0.83680, 0.91050,  &
		0.98420, 1.05790, 1.13160, 1.20530, 1.27890, 1.35260,  &
		1.42630, 1.50000, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(17,:)=(/ &
		0.09619, 0.16874, 0.24290, 0.31831, 0.39481, 0.47196,  &
		0.54941, 0.62680, 0.70378, 0.77998, 0.85508, 0.92870,  &
		1.00000, 1.07170, 1.14315, 1.21436, 1.28532, 1.35606,  &
		1.42660, 1.49693, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
	db(18,:)=(/ &
		0.11843, 0.19976, 0.27851, 0.35537, 0.43084, 0.50495,  &
		0.57807, 0.65018, 0.72157, 0.79217, 0.86204, 0.93136,  &
		1.00000, 1.06827, 1.13590, 1.20309, 1.26983, 1.33620,  &
		1.40215, 1.46768, 0.00000, 0.00000, 0.00000, 0.00000,  &
		0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000  &
		/)
endif
endif 

!Now select which to interpolate
n=30  !Initially set the array length to 30. Resets to 20 in special case of interaction
select case(XT)
case(1) ; XI=1    !A
case(2) ; XI=4    !B
case(3) ; XI=7    !C
case(12) ; XI=10 ; n=20   !AB
case(13) ; XI=13 ; n=20   !AC
case(23) ; XI=16 ; n=20   !BC
end select

select case(YT)
case(11) ;  YI=2      !PA
case(12) ;  YI=5      !PB
case(13) ;  YI=8      !PC
case(112) ;  YI=11    !PAB
case(113) ;  YI=14    !PAC
case(123) ;  YI=17    !PBC
case(21) ;  YI=3      !QA
case(22) ;  YI=6      !QB
case(23) ;  YI=9      !QC
case(212) ;  YI=12    !QAB
case(213) ;  YI=15    !QAC
case(223) ;  YI=18    !QBC
end select

!Set the data to be interpolated
datx(:)=db(XI,1:n)
daty(:)=db(YI,1:n)
	
!Use brute force interpolation.. it is faster in this case than
!..bisection or hunting methods used in the user-specified HTF case		
do i=1,n
    !if the index is equal to n, then set bounds and end loop
    if(i.eq.n) then
        lbi=n
        ubi=n
        exit
    endif

    !if the x variable is outside the table range, set the bounds and get out
    if(i.eq.1) then
        if(datx(2).gt.datx(1)) then !The table is in ascending order
            if(X.le.datx(1)) then
                lbi=1; ubi=1; exit;
            endif
            if(X.ge.datx(n)) then
                lbi=n; ubi=n; exit;
            endif
        else  !the table is in descending order
            if(x.ge.datx(1)) then
                lbi=1; ubi=1; exit;
            endif
            if(X.le.datx(n)) then
                lbi=n; ubi=n; exit;
            endif
        endif
    endif


    if(((X.ge.datx(i)).and.(X.lt.datx(i+1))).or.((X.le.datx(i)).and.(X.gt.datx(i+1)))) then
        lbi=i
        ubi=i+1
        exit
    endif
enddo

if(datx(ubi).eq.datx(lbi)) then
    ind=0.
else
    ind=(X-datx(lbi))/(datx(ubi)-datx(lbi))
endif

Y = daty(lbi)+ind*(daty(ubi)-daty(lbi))

Interpolate=Y


end function Interpolate 


!************************************************************************************************************
!*                                      Cooling subroutines                                                 *
!************************************************************************************************************

subroutine evap_tower(DeltaT_cw_des, T_approach, P_cycle, eta_ref, T_db, T_wb, P_amb, q_reject, m_dot_water, &
                      W_dot_tot, P_cond, T_cond)

!use wet_cooling_prop_funcs
implicit none

real(8):: C_AIR, C_CW, DELTAH_EVAP, DELTAT_CW, DeltaT_cw_des, DP_EVAP, DRIFT_LOSS_FRAC, DT_OUT, eta_ref, ETA_FAN, ETA_FAN_S, &
          ETA_PCW_S, ETA_PUMP, H_FAN_IN, H_FAN_IN_CHK, H_FAN_OUT, H_FAN_OUT_S, H_PCW_IN, H_PCW_OUT, &
          H_PCW_OUT_S, M_DOT_AIR, M_DOT_BLOWDOWN, M_DOT_CW, m_dot_cw_des, M_DOT_DRIFT, blowdown_frac, &
          M_DOT_EVAP, M_DOT_WATER, MASS_RATIO_FAN, P_AMB, P_COND, P_CYCLE, P_RATIO_FAN, Q_REJECT, Q_reject_des, R, RHO_CW, &
          S_FAN_IN, S_PCW_IN, T_APPROACH, T_COND, T_DB, T_FAN_IN, T_FAN_IN_K, T_FAN_OUT, T_FAN_OUT_K, &
          T_WB, W_DOT_CW_PUMP, W_DOT_FAN, W_DOT_TOT

!------------------------------------------------------------------------------------------------------------
!--Inputs
!   * DeltaT_cw_des [K]     Cooling water temperature rise across condenser
!   * T_approach    [K]     Cooling tower approach temperature, difference between cw out and wet bulb temp
!   * P_cycle       [W]     Rated power block capacity
!   * eta_ref       [-]     Rated gross conversion efficiency
!   * T_db          [K]     Dry bulb temperature (converted to C)
!   * P_amb         [Pa]    Atmospheric pressure
!------------------------------------------------------------------------------------------------------------
!--Output
!   * m_dot_water   [kg/s]  Total cooling tower water usage
!   * W_dot_tot     [MW]    Total parasitic power for cooling tower model
!   * P_cond        [Pa]    Condenser steam pressure
!   * T_cond        [K]     Condenser steam temperature
!------------------------------------------------------------------------------------------------------------
!Unit conversions
T_db = T_db - 273.15    ![C] Converted dry bulb temp
T_wb = T_wb - 273.15    ![C] Converted wet bulb temp

! Values that can be estimated
dt_out = 3.                 !Temperature difference at hot side of the condenser
drift_loss_frac = 0.001     !Drift loss fraction
blowdown_frac = 0.003       !Blowdown fraction
dP_evap = 0.37*1.e5         ![Pa] Pressure drop across the condenser and cooling tower
eta_pump = 0.75             !Total pump efficiency
eta_pcw_s = 0.8             !Isentropic cooling water pump efficiency
eta_fan = 0.75              !Fan mechanical efficiency
eta_fan_s = 0.8             !Fan isentropic efficiency
P_ratio_fan = 1.0025        !Fan pressure ratio
mass_ratio_fan = 1.01       !Ratio of air flow to water flow in the cooling tower

!Cooling water specific heat
c_cw = f_c_psat(P_amb)

!**** Calculations for design conditions
Q_reject_des = P_cycle*(1./eta_ref-1.)    	    !Heat rejection from the cycle
m_dot_cw_des = Q_reject_des/(c_cw*DeltaT_cw_des)	!Mass flow rate of cooling water required to absorb the rejected heat
!**** 

!Based on the required heat rejection load, use either full flow rate or 50% flow rate
if(Q_reject >= Q_reject_des/2.) then
    m_dot_cw = m_dot_cw_des
else
    m_dot_cw = m_dot_cw_des/2.
endif

!**** Calculations for performance
DeltaT_cw = Q_reject/(m_dot_cw*c_cw)

!Condenser saturation temperature
T_cond = T_wb + DeltaT_cw + dt_out + T_approach

!Condenser back pressure
P_cond = f_psat_T(T_cond)

!Circulating water pump power
h_pcw_in = f_hw_psat(P_amb)     ![J/kg] cw pump inlet enthalpy 
s_pcw_in = f_s_hw_psat(P_amb)     ![J/kg-K] cw pump inlet entropy
rho_cw = f_rho_P(P_amb)         ![kg/m3] cooling water density in the pump
h_pcw_out_s = dP_evap/rho_cw + h_pcw_in                         ![J/kg] isentropic outlet enthalpy.. incompressible fluid
h_pcw_out = h_pcw_in + (h_pcw_out_s - h_pcw_in)/eta_pcw_s       ![J/kg] Outlet enthalpy accounting for irreversibility
W_dot_cw_pump = (h_pcw_out - h_pcw_in)*m_dot_cw/eta_pump*1.e-6  ![MW] Cooling water circulating pump power

!Fan power
m_dot_air = m_dot_cw*mass_ratio_fan
T_fan_in = (T_db + T_wb + T_approach)/2
h_fan_in = f_h_air_T(T_fan_in)

c_air = 1003.               ![J/kg-K] specific heat of air (This is relatively constant)
R = 8314./28.97             ![J/kmol-K]/[kg/kmol] Gas constant over the molar mass of air

T_fan_in_K = T_fan_in + 273.15  !Fan inlet temp, in K
T_fan_out_K = T_fan_in_K*P_ratio_fan**(R/C_air)    ![K] isentropic temperature rise
T_fan_out = T_fan_out_K - 273.15    ![C] Convert isentropic temperature rise to deg C
h_fan_out_s = f_h_air_T(T_fan_out)  ![J/kg] Calculate isentropic enthalpy at fan outlet
h_fan_out = h_fan_in + (h_fan_out_s - h_fan_in)/eta_fan_s   ![J/kg] Actual enthalpy, accounting for irreversibility

W_dot_fan = (h_fan_out - h_fan_in)*m_dot_air/eta_fan*1.e-6  ![MW] Fan parasitic power

!Total cooling tower parasitic power
W_dot_tot = W_dot_cw_pump + W_dot_fan   ![MW] 

!Enthalpy of evaporation
DELTAh_evap = f_dh_evap(P_amb)

!Evaporative water loss
m_dot_evap = Q_reject/Deltah_evap

!Other water losses
m_dot_drift = drift_loss_frac*m_dot_cw	!Drift loss fraction, based on cooling water mass flow rate
m_dot_blowdown = blowdown_frac*m_dot_cw	!Blow down fraction

!Total power block water usage
m_dot_water = m_dot_evap + m_dot_drift + m_dot_blowdown

!Unit conversions
T_db = T_db + 273.15    ![C] Converted dry bulb temp
T_wb = T_wb + 273.15    ![C] Converted wet bulb temp
T_cond = T_cond + 273.15    ![K] Convert to K for output


!--------- Declare the property functions used in this subroutine --------------
contains

    real(8) function f_dh_evap(P)
        real(8),intent(in)::P
	    !Calculates enthalpy of evaporation given an atmospheric pressure
	    !Input: P -> Pressure [Pa]
	    !Output: Enthalpy [J/kg]
	    f_dh_evap = 2.36230E+06 - 1.35459*P + 0.00000308492*P**2
    end function

    real(8) function f_c_psat(P)
        real(8),intent(in)::P
	    !calculates the specific heat of water at the saturated liquid state as a real(8) function of pressure
	    !Input: P-> Pressure [Pa]
	    !Output: Specific heat [J/kg-K]
	    f_c_psat= 4170.46 + 0.000538088*P - 7.73437E-10*P**2
    end function

    real(8) function f_psat_T(T)
	    real(8),intent(in)::T
	    !Calculates the saturation pressure of steam given a certain temperature
	    !Input: T -> Temperature [K]
	    !Output: Pressure [Pa]
	    f_psat_T = 1125.09 - 19.6444*T + 4.42596*T**2 - 0.0391851*T**3 + 0.000965517*T**4
    end function

    real(8) function f_hw_psat(P)
	    real(8),intent(in)::P
	    !Calculates the enthalpy of water at liquid saturation near ambient pressures
	    !Input: P -> Pressure [Pa]
	    !Output: enthalpy [J/kg]
	    f_hw_psat = 229628.719 + 2.78471579*P - 0.0000111907252*P**2 + 2.12030100E-11*P**3
    end function

    real(8) function f_s_hw_psat(P)
	    real(8),intent(in)::P
	    !Calculates the entropy of water at liquid saturation near ambient pressures
	    !Input: P -> Pressure [Pa]
	    !Output: entropy [J/kg-K]
	    f_s_hw_psat = 779.989872 + 0.00791597131*P - 3.33033640E-08*P**2 + 6.38602593E-14*P**3
    end function

    real(8) function f_rho_P(P)
	    real(8),intent(in)::P
	    !Calculates density of water at liquid saturation near ambient pressures
	    !Input: P -> Pressure [Pa]
	    !Output: density [kg/m3]
	    f_rho_P = 984.079732 - 0.000307058016*P + 5.32272340E-10*P**2
    end function

    real(8) function f_h_air_T(T)
	    real(8),intent(in)::T
	    !Calculates enthalpy of air as a real(8) function of temperature
	    !Input: T -> Temperature [C]
	    !Output: enthalpy [J/kg]
	    f_h_air_T=273474.659 + 1002.9404*T + 0.0326819988*T**2
    end function


end subroutine

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


subroutine ACC(T_ITD_des, P_cond_ratio, P_cycle, eta_ref, T_db, P_amb, q_reject, m_dot_air, W_dot_fan, P_cond, T_cond)

implicit none

real(8)::dt_fan, eta_ref, eta_fan, eta_fan_s, h_fan_in, h_fan_out, h_fan_out_s, m_dot_air, &
         mm, p_amb, p_cond, p_cond_ratio, p_cycle, q_reject, r, rh, t_cond, t_db, t_fan_in_k, t_fan_out, t_fan_out_k,&
         t_hot_diff, t_itd, w_dot_fan, t_itd_des, c_air, q_reject_des, m_dot_air_des

!------------------------------------------------------------------------------------------------------------
!--Inputs
!   * T_ITD         [K]     ACC initial temperature difference, difference between dry bulb and steam inlet temp
!   * P_cond_ratio  [-]     Condenser air inlet/outlet pressure ratio
!   * P_cycle       [W]     Rated power block capacity
!   * eta_ref       [-]     Rated gross conversion efficiency
!   * T_db          [K]     Dry bulb temperature (converted to C)
!   * P_amb         [Pa]    Atmospheric pressure
!------------------------------------------------------------------------------------------------------------
!--Output
!   * m_dot_air     [kg/s]  Total ACC air mass flow rate
!   * W_dot_fan     [MW]    Total parasitic power for ACC model
!   * P_cond        [Pa]    Condenser steam pressure
!   * T_cond        [K]     Condenser steam temperature
!------------------------------------------------------------------------------------------------------------
!Unit conversions
T_db = T_db - 273.15        ![C] Converted dry bulb temp

! Values that can be estimated
T_hot_diff = 3.             ![C] Temperature difference between saturation steam and condenser outlet air temp
eta_fan_s = 0.8             ![-] Fan isentropic efficiency
eta_fan = .98**3            ![-] Fan mechanical efficiency
c_air = 1005.               ![J/kg-K] Specific heat of air, relatively constant over dry bulb range

!**** Calculations for design conditions
Q_reject_des = P_cycle*(1./eta_ref-1.)    	    !Heat rejection from the cycle
m_dot_air_des = Q_reject_des/(c_air*(T_ITD_des - T_hot_diff))	
!**** 

!Based on the required heat rejection load, use either full flow rate or 50% flow rate
if(Q_reject >= Q_reject_des/2.) then
    m_dot_air = m_dot_air_des
else
    m_dot_air = m_dot_air_des/2.
endif

!Fan power
T_ITD = T_hot_diff + q_reject/(m_dot_air*c_air)  ![C] Calculate the actual ITD during off-design operation

! Calculated output
T_cond = T_db + T_ITD	    !Condensation temperature
P_cond =  f_psat_T(T_cond)	!Turbine back pressure

h_fan_in = f_h_air_T(T_db)  ![J/kg] Fan inlet enthalpy

mm = 28.97                  ![kg/kmol] molar mass of air
R = 8314./mm                ![J/kg-K] Gas constant for air

!These temperature calculations are for the isentropic expansion across the fan, not accounting for heat gain in the ACC
T_fan_in_K = T_db + 273.15  ![K] Fan inlet temperature
T_fan_out_K = T_fan_in_K*P_cond_ratio**(R/C_air)
T_fan_out = T_fan_out_K - 273.15    ![C] Fan outlet temperature
dT_fan = T_fan_out - T_db   ![C] Difference in temperature caused by irreversibilities in fan

h_fan_out_s = f_h_air_T(T_fan_out)	![J/kg] Isentropic fan outlet temperature
h_fan_out = h_fan_in + (h_fan_out_s - h_fan_in)/eta_fan_s   ![J/kg] Actual fan outlet temperature
!Total ACC parasitic power
W_dot_fan = (h_fan_out - h_fan_in)*m_dot_air/eta_fan*1.e-6  ![MW] Fan power

!Unit conversions
T_db = T_db + 273.15    ![C] Converted dry bulb temp
T_cond = T_cond + 273.15    ![K] Convert to K for output

!--------- Declare the property functions used in this subroutine --------------
contains

    real(8) function f_psat_T(T)
        real(8),intent(in)::T
	    !Calculates the saturation pressure of steam given a certain temperature
	    !Input: T -> Temperature [K]
	    !Output: Pressure [Pa]
	    f_psat_T = 1125.09 - 19.6444*T + 4.42596*T**2 - 0.0391851*T**3 + 0.000965517*T**4
    end function

    real(8) function f_h_air_T(T)
	    real(8),intent(in)::T
	    !Calculates enthalpy of air as a function of temperature
	    !Input: T -> Temperature [C]
	    !Output: enthalpy [J/kg]
	    f_h_air_T=273474.659 + 1002.9404*T + 0.0326819988*T**2
    end function

end subroutine
