SUBROUTINE TYPE251 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
!************************************************************************
! Object: Detailed parabolic trough plant controller
! Simulation Studio Model: Type251
! 
! Author: Michael J. Wagner
! Editor: Michael J. Wagner
! Date:	 December 12, 2009
! Modified: March 19, 2010

! COPYRIGHT 2010 NATIONAL RENEWABLE ENERGY LABORATORY


!Parameters
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable            | Description                                                                | Input Units  | Internal Units
!-----------------------------------------------------------------------------------------------------------------------------------------
!  1|  LU_fl2              | Fluid property file logical unit                                           | none         | none
!  2|  TSHOURS             | [hr] Equivalent full-load thermal storage hours                            | hr           | hr
!  3|  NUMTOU              | [-] Number of time-of-use periods                                          | none         | none
!  4|  is_hx               | 1=yes, 0=no                                                                | none         | none
!  5|  dT_hot              | [K] Hot side heat exchanger approach temp                                  | C            | C
!  6|  dT_cold             | [K] Cold side heat exchanger approach temp                                 | C            | C
!  7|  hx_config           | [-] Heat exchanger configuration                                           | none         | none
!  8|  q_max_aux           | [MWt] Maximum heat rate of the auxiliary heater                            | MWt          | Wt
!  9|  LHV_eff             | [-] Fuel LHV efficiency (0..1)                                             | none         | none
! 10|  T_set_aux           | [C] Aux heater outlet temperature set point                                | C            | K
! 11|  store_fluid         | [-] material number for the storage fluid (see type229)                    | none         | none
! 12|  V_tank_hot0         | [m3] Initial hot tank fluid volume                                         | m3           | m3
! 13|  T_tank_hot0         | [K] Initial hot tank fluid temperature                                     | C            | K
! 14|  T_tank_cold0        | [K] Initial cold tank fluid temperature                                    | C            | K
! 15|  vol_tank            | [m3] Total tank volume, including the unusable volume of HTF at the bottom | m3           | m3
! 16|  h_tank              | [m] Total height of the tank (height of HTF when tank is full)             | m            | m
! 17|  h_tank_min          | [m] Minimum htf height in the storage tank                                 | m            | m
! 18|  u_tank              | [W/m2-K] Loss coefficient from the tank                                    | W/m2-K       | W/m2-K
! 19|  tank_pairs          | [-] Number of equivalent tank pairs                                        | none         | none
! 20|  tank_htr_set_point  | [C] Minimum allowable fluid temperature before auxiliary heater turns on   | C            | K
! 21|  tank_max_heat       | [W] Rated heater capacity for tank heating                                 | MW           | W
! 22|  eta_heater_tank     | [-] Conversion efficiency for electric heater                              | none         | none
! 23|  field_fluid         | [-] material number for the collector field (see type229)                  | none         | none
! 24|  T_field_in_des      | [K] Field design inlet temperature                                         | C            | K
! 25|  T_field_out_des     | [K] Field loop outlet design temperature                                   | C            | K
! 26|  q_pb_design         | [MWt] Design heat input to the power block                                 | MWt          | Wt
! 27|  W_pb_design         | [MWe]  Rated plant capacity                                                | MWe          | We
! 28|  cycle_max_fraction  | [-] Maximum turbine over design operation fraction                         | none         | none
! 29|  cycle_cutoff_frac   | [-] Minimum turbine operation fraction before shutdown                     | none         | none
! 30|  solarm              | [-] Solar multiple                                                         | none         | none
! 31|  PB_pump_coef        | [kJ/kg] Pumping power required to move 1kg of HTF through power block flow | kW/kg        | kW/kg
! 32|  TES_pump_coef       | [kJ/kg] Pumping power required to move 1kg of HTF through storage flow loop| kW/kg        | kW/kg
! 33|  PB_fixed_par        | [-] fraction of rated gross power consumed at all hours of the year        | none         | none
! 34|  BOP_parVal          | [MWe/MWcap] Balance of plant parasitic                                     | MWe/MWcap    | MWe/MWcap
! 35|  BOP_parPF           | [-] Balance of plant parasitic multiplier factor                           | none         | none
! 36|  BOP_par0            | Balance of plant parasitic polynomial constant                             | none         | none
! 37|  BOP_par1            | Balance of plant parasitic polynomial linear term                          | none         | none
! 38|  BOP_par2            | Balance of plant parasitic polynomial quadratic term                       | none         | none
! 39|  Aux_parVal          | [MWe/MWcap] Aux heater/boiler parasitic                                    | MWe/MWcap    | MWe/MWcap
! 40|  Aux_parPF           | [-] Aux heater/boiler parasitic multiplier factor                          | none         | none
! 41|  Aur_par0            | Aux heater/boiler parasitic polynomial constant                            | none         | none
! 42|  Aur_par1            | Aux heater/boiler parasitic polynomial linear term                         | none         | none
! 43|  Aur_par2            | Aux heater/boiler parasitic polynomial quadratic term                      | none         | none
! 44|  T_startup           | [C] Startup temperature (same as field startup)                            | C            | K
! 45|  TSLOGIC_0           | [-] Dispatch logic 1 - no solar                                            | none         | none
! 46|  TSLOGIC_1           | [-] Dispatch logic 1 - with solar                                          | none         | none
! 47|  TSLOGIC_2           | [-] Dispatch logic 1 - Turbine load fraction                               | none         | none
! 48|  TSLOGIC_3           | [-] Dispatch logic 2 - no solar                                            | none         | none
! 49|  TSLOGIC_4           | [-] Dispatch logic 2 - with solar                                          | none         | none
! 50|  TSLOGIC_5           | [-] Dispatch logic 2 - Turbine load fraction                               | none         | none
! 51|  TSLOGIC_6           | [-] Dispatch logic 3 - no solar                                            | none         | none
! 52|  TSLOGIC_7           | [-] Dispatch logic 3 - with solar                                          | none         | none
! 53|  TSLOGIC_8           | [-] Dispatch logic 3 - Turbine load fraction                               | none         | none
! 54|  TSLOGIC_9           | [-] Dispatch logic 4 - no solar                                            | none         | none
! 55|  TSLOGIC_10          | [-] Dispatch logic 4 - with solar                                          | none         | none
! 56|  TSLOGIC_11          | [-] Dispatch logic 4 - Turbine load fraction                               | none         | none
! 57|  TSLOGIC_12          | [-] Dispatch logic 5 - no solar                                            | none         | none
! 58|  TSLOGIC_13          | [-] Dispatch logic 5 - with solar                                          | none         | none
! 59|  TSLOGIC_14          | [-] Dispatch logic 5 - Turbine load fraction                               | none         | none
! 60|  TSLOGIC_15          | [-] Dispatch logic 6 - no solar                                            | none         | none
! 61|  TSLOGIC_16          | [-] Dispatch logic 6 - with solar                                          | none         | none
! 62|  TSLOGIC_17          | [-] Dispatch logic 6 - Turbine load fraction                               | none         | none
! 63|  FFRAC_0             | [-] Fossil dispatch logic 1                                                | none         | none
! 64|  FFRAC_1             | [-] Fossil dispatch logic 2                                                | none         | none
! 65|  FFRAC_2             | [-] Fossil dispatch logic 3                                                | none         | none
! 66|  FFRAC_3             | [-] Fossil dispatch logic 4                                                | none         | none
! 67|  FFRAC_4             | [-] Fossil dispatch logic 5                                                | none         | none
! 68|  FFRAC_5             | [-] Fossil dispatch logic 6                                                | none         | none

!Inputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable            | Description                                                                | Input Units  | Internal Units
!-----------------------------------------------------------------------------------------------------------------------------------------
!  1|  I_bn                | [W/m2] Beam normal radiation (input kJ/m2-hr)                              | kJ/m2-hr     | W/m2
!  2|  TOUPeriod           | Time of use period                                                         | none         | none
!  3|  m_dot_field         | [kg/s]  Flow rate from the field                                           | kg/hr        | kg/s
!  4|  m_dot_htf_ref       | [kg/s] Reference HTF flow rate at design conditions                        | kg/hr        | kg/s
!  5|  T_field_out         | [K] HTF Temperature from the field                                         | C            | K
!  6|  T_pb_out            | [K] Fluid temperature from the power block                                 | C            | K
!  7|  T_amb               | [K] Ambient temperature                                                    | C            | K
!  8|  m_pb_demand         | [kg/s] Demand htf flow from the power block                                | kg/hr        | kg/s
!  9|  q_startup           | [MWt] Startup energy reported by the collector field                       | MWt          | MWt

!Outputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable            | Description                                                                | Output Units | Internal Units
!-----------------------------------------------------------------------------------------------------------------------------------------
!  1|  defocus             | Defocus instigated by control decisions                                    | none         | none
!  2|  cycle_pl_control    | Part-load control flag - used by Type224                                   | none         | none
!  3|  standby_control     | Standby control flag - used by Type224                                     | none         | none
!  4|  m_dot_pb            | Mass flow rate to the power block                                          | kg/hr        | kg/s
!  5|  T_pb_in             | [C] HTF Temperature to the power block                                     | C            | K
!  6|  T_field_in          | [C] HTF temperature into the collector field header                        | C            | K
!  7|  no name             | [kg/hr] Mass flow rate of the collector side storage loop                  | kg/hr        | kg/hr
!  8|  no name             | [kg/hr] Mass flow rate of the tank side storage loop                       | kg/hr        | kg/hr
!  9|  Ts_hot              | [C] Fluid temperature of the hot HTF in the collector loop                 | C            | K
! 10|  Ts_cold             | [C] Fluid temperature of the cold HTF in the collector loop                | C            | K
! 11|  T_tank_hot_in       | [C] Hot tank inlet temperature                                             | C            | K
! 12|  T_tank_cold_in      | [C] Cold tank inlet temperature                                            | C            | K
! 13|  v_tank_hotf         | [m3] Available Hot tank HTF volume at the end of the timestep              | m3           | m3
! 14|  T_tank_hotf         | [C] Hot tank stored HTF temperature at the end of the timestep             | C            | K
! 15|  v_tank_coldf        | [m3] Available Cold tank HTF volume at the end of the timestep             | m3           | m3
! 16|  T_tank_coldf        | [C] Cold tank stored HTF temperature at the end of the timestep            | C            | K
! 17|  no name             | [MWe] Total parasitic power required for tank freeze protection            | MWe          | MWe
! 18|  m_dot_aux           | [kg/hr] Auxiliary heater mass flow rate                                    | kg/hr        | kg/s
! 19|  q_aux_fluid         | [MWt] Thermal energy provided to the fluid passing through the aux heater  | MWt          | MWt
! 20|  q_aux_fuel          | [MMBTU] Heat content of fuel required to provide aux firing                | MMBTU        | MMBTU
! 21|  no name             | [m3] Total accounted htf volume in storage                                 | m3           | m3
! 22|  hx_eff              | [-] Heat exchanger effectiveness                                           | none         | none
! 23|  m_tank_hotf         | [kg] mass of total fluid in the hot tank                                   | kg           | kg
! 24|  m_tank_coldf        | [kg] mass of total fluid in the cold tank                                  | kg           | kg
! 25|  no name             | [kg] total mass in the storage system                                      | kg           | kg
! 26|  htf_pump_power      | [MWe] pumping power for storage, power block loops                         | MWe          | MWe
! 27|  BOP_par             | [MWe] parasitic power as a function of power block load                    | MWe          | MWe
! 28|  no name             | [MWe] Fixed parasitic power losses.. for every hour of operation           | MWe          | MWe
! 29|  Aux_par             | [MWe] Parasitic power associated with operation of the aux boiler          | MWe          | MWe
! 30|  q_pb                | [MWt] Thermal energy to the power block                                    | MWt          | Wt


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

!-----------------------------------------------------------------------------------------------------------------------
! required by the multi-dll version of TRNSYS
!DEC$ATTRIBUTES DLLEXPORT :: TYPE251

implicit none 

!TRNSYS declarations
real(8):: time
integer*4:: info(15), iunit, itype, icntrl
integer*4,parameter::np=100,ni=9,nout=31,nd=0,ns=8

!Dimension the TRNSYS variables
real(8)::xin(ni),out(nout),par(np),stored(ns),T(nd),dtdt(nd) 

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


!---- User variable declarations----------------------------------------------------------------------------------------
!    PARAMETERS
real(8) :: TSLOGIC_A(6),TSLOGIC_B(6),TSLOGIC_C(6),FFRAC(6), TSLOGIC(18), TSHOURS, cycle_max_fraction, cycle_cutoff_frac,&
           V_tank_hot_ini, T_tank_hot_ini, V_tank_cold_ini, T_tank_cold_ini, T_field_in_des, T_field_out_des, &
           LHV_eff, T_set_aux, vol_tank, h_tank, u_tank, tank_pairs, tank_htr_set_point, eta_heater_tank,&
           q_max_aux, q_pb_design, solarm, dT_hot, dT_cold, store_fluid, field_fluid, h_tank_min, tank_max_heat,&
           W_pb_design, PB_pump_coef, TES_pump_coef
integer :: NUMTOU, hx_config, LU_fl
logical :: is_hx

! INPUTS
real(8) :: I_bn, m_dot_field, m_dot_htf_ref, T_field_out, T_pb_out, T_amb, m_pb_demand, q_startup
integer :: TOUPeriod


! OUTPUTS

! Locals
real(8) :: TSELECT(6), ms_charge_avail, ms_disch_avail, V_tank_hot0, T_tank_hot0, V_tank_cold0, T_tank_cold0, m_dot_aux_avail,&
           density, specheat, dt, c_htf_aux, c_htf_chg, c_htf_dis, c_htf_fld, T_tank_hot_in, T_tank_hot_out, T_tank_cold_in, T_tank_cold_out, &
           m_tank_disch_avail, m_tank_charge_avail, Ts_hot, Ts_cold, duty, hx_eff_des, hx_ua, m_avail_tot, hx_eff, q_hx,&
           V_tank_active, qs_disch_avail, qs_charge_avail, qaux_avail, qfield_avail, q_int, T_int, m_int, ms_disch, ms_charge, &
           m_dot_aux, q_aux, m_dot_pb, T_pb_in, defocus, cycle_pl_control, standby_control, T_field_in, T_field_inX, &
           err, tol, q_pb_demand, q_pb, m_tank_charge, m_tank_disch, T_tank_hot_ave, vol_hot_tank_ave,&
           q_loss_hot_tank, T_tank_hotf, V_tank_hotf, q_tank_hot_htr, T_tank_cold_ave, vol_cold_tank_ave, q_loss_cold_tank, &
           T_tank_coldf, v_tank_coldf, q_tank_cold_htr, T_aux_act, q_aux_fluid, q_aux_fuel, demand_var, &
           m_tank_hotf, m_tank_coldf, m_tank_hot0, m_tank_cold0, v_tank_hotfx, v_tank_coldfx, &
           errt, T_tank_hotX, T_tank_coldX, V_tank_hotA, V_tank_coldA, rho_tank_coldX, rho_tank_hotX, msghx1, msghx2, msghx3, msgiter,&
           derr, errx, defocusx, htf_pump_power, PB_fixed_par, BOP_parVal, BOP_parPF, BOP_par0, BOP_par1, BOP_par2, BOP_par,&
           Aux_parVal, Aux_parPF, Aux_par0, Aux_par1, Aux_par2, Aux_par, T_startup, T_field_out0
integer :: mode, up1, stored_mode, iter, itert
character::message(10)*300

integer :: i, dumi

!-----------------------------------------------------------------------------------------------------------------------
!set all the messages for this type
message(1) = "An non-physical solution was found for the heat exchanger size. Please correct temperature specifications."  !Fatal
message(2) = "The simulation has encountered a non-physical condition in the heat exchanger. Please adjust design values." !Warning
!message(3) reserved for iteration. Value is set at iteration calculation.  Warning. 

! --- Initial call to detect the TRNSYS version for which this Type is written -----------------------------------------
if (info(7) .eq. -2) then
    info(12) = 16   ! This component is a TRNSYS 16 Type
    return 1
endif


! --- Very last call in simulation -------------------------------------------------------------------------------------

if (info(8) == -1) then

    return 1    ! Exit 
    
endif


! --- Post-convergence call --------------------------------------------------------------------------------------------
if (info(13) == 1) then

    !******************************************************************
    ! Set the system state values for the next timestep 
    !******************************************************************
    stored(1) = v_tank_hotf     ![m3] Final hot tank HTF volume
    stored(2) = T_tank_hotf     ![K] Final hot tank HTF temperature
    stored(3) = v_tank_coldf    ![m3] Final cold tank HTF volume
    stored(4) = T_tank_coldf    ![K] Final cold tank HTF temperature
    stored(5) = dble(mode)      !operation mode
    stored(6) = m_tank_hotf     ![kg] hot tank available mass
    stored(7) = m_tank_coldf    ![kg] cold tank available mass
    stored(8) = xin(5)+273.15   ![K] T_field_out
            
    call setStorageVars(stored,nS,INFO)
    
    !Print out error messages
    if(msghx1 /= 0.) then
        call messages(-1,trim(message(msghx1)),"Fatal",info(1),info(2))
    endif
    if(msghx2 /= 0.) then
        call messages(-1,trim(message(msghx2)),"Warning",info(1),info(2))
    endif
    if(msghx3 /= 0.) then
        call messages(-1,trim(message(msghx3)),"Warning",info(1),info(2))
    endif
    if(msgiter /= 0.) then
        call messages(-1,trim(message(msgiter)),"Warning",info(1),info(2))
    endif


    return 1    ! Exit - End of the routine for post-convergence calls
    
endif


! --- Second call in simulation: initialization call (not a simulation call) -------------------------------------------
if (info(7) .eq. -1) then
    message(:) = ''

    dt = getSimulationTimeStep()*3600.  ![s]
    
    !initialize incremental counter function. This increments once each time
    !it is called with arg 1. handy for reading in parameters without having to manually
    !assign each one a parameter number
    dumi=up1(0)
    
    !Initialize TSLOGIC 
    TSLOGIC=0.; TSLOGIC_A=0.; TSLOGIC_B=0.; TSLOGIC_C=0.;
    
    !---- read static parameters---
    !Logic and control parameters
    LU_fl = int(par(up1(1)))                !Fluid property file logical unit
    TSHOURS=int(PAR(up1(1)))                ![hr] Equivalent full-load thermal storage hours
    NUMTOU=int(PAR(up1(1)))                 ![-] Number of time-of-use periods
    
    !heat exchanger parameters
    ![-] Does the system use a heat exchanger?
    if(par(up1(1)) == 1.) then; is_hx=.true.; else; is_hx = .false.; endif 
    dT_hot = par(up1(1))                    ![K] Hot side heat exchanger approach temp
    dT_cold = par(up1(1))                   ![K] Cold side heat exchanger approach temp
    hx_config = int(par(up1(1)))            ![-] Heat exchanger configuration
                                            !{1=Same fluid, 2=Counter-flow, 3=Parallel-flow, 4=Cross-flow}
    !Auxiliary fossil heater parameters
    q_max_aux = par(up1(1))*1.e6            ![Wt] Maximum heat rate of the auxiliary heater
    LHV_eff = par(up1(1))                   ![-] Fuel LHV efficiency (0..1)
    T_set_aux = par(up1(1))+273.15          ![K] Aux heater outlet temperature set point
    
    !Storage tank parameters
    store_fluid = par(up1(1))               ![-] material number for the storage fluid (see type229)
    V_tank_hot_ini = par(up1(1))            ![m3] Initial hot tank fluid volume
    T_tank_hot_ini = par(up1(1))+273.15     ![K] Initial hot tank fluid temperature   
    T_tank_cold_ini = par(up1(1))+273.15    ![K] Initial cold tank fluid temperature
    vol_tank = par(up1(1))                  ![m3] Total tank volume, including the unusable volume of HTF at the bottom
    h_tank = par(up1(1))                    ![m] Total height of the tank (height of HTF when tank is full)
    h_tank_min = par(up1(1))                ![m] Minimum htf height in the storage tank
    u_tank = par(up1(1))                    ![W/m2-K] Loss coefficient from the tank
    tank_pairs = par(up1(1))                ![-] Number of equivalent tank pairs
    tank_htr_set_point = par(up1(1))+273.15 ![K] Minimum allowable fluid temperature before auxiliary heater turns on
    tank_max_heat = par(up1(1))*1.e6        ![W] Rated heater capacity for tank heating
    eta_heater_tank = par(up1(1))           ![-] Conversion efficiency for electric heater
    
    V_tank_cold_ini = vol_tank - V_tank_hot_ini   ![m3] Initial cold tank fluid volume
    V_tank_active = vol_tank*(1.-2.*h_tank_min/h_tank) ![m3] Active tank volume.. that is, volume above the 
                                                       !minimum fluid level and below the maximum level
    
    !Field parameters
    field_fluid = par(up1(1))               ![-] material number for the collector field (see type229)
    T_field_in_des = par(up1(1))+273.15     ![K] Field design inlet temperature
    T_field_out_des = par(up1(1))+273.15    ![K] Field loop outlet design temperature
    !Check to make sure a fluid mismatch isn't present
    if(.not.is_hx) store_fluid = field_fluid

    !Power block parameters
    q_pb_design = par(up1(1))*1.e6          ![Wt] Design heat input to the power block
    W_pb_design = par(up1(1))*1.e6          ![We] Design electrical output of the power block
    cycle_max_fraction = par(up1(1))        ![-] Maximum turbine over design operation fraction
    cycle_cutoff_frac = par(up1(1))         ![-] Minimum turbine operation fraction
    solarm = par(up1(1))                    ![-] Solar multiple
    PB_pump_coef = par(up1(1))              ![kW/kg] Flow through power block pump power coeffcient
    TES_pump_coef = par(up1(1))             ![kW/kg] Flow through storage pump power coefficient\
    PB_fixed_par = par(up1(1))              ![-] Fraction.. fixed power block losses. Constant throughout year.
    BOP_parVal = par(up1(1))                ![MWe/MWcap] Balance of plant parasitic 
    BOP_parPF = par(up1(1))                 ![-] Balance of plant parasitic multiplier factor
    BOP_par0 = par(up1(1))                  !Balance of plant parasitic polynomial constant
    BOP_par1 = par(up1(1))                  !Balance of plant parasitic polynomial linear term
    BOP_par2 = par(up1(1))                  !Balance of plant parasitic polynomial quadratic term
    Aux_parVal = par(up1(1))                ![MWe/MWcap] Aux heater/boiler parasitic 
    Aux_parPF = par(up1(1))                 ![-] Aux heater/boiler parasitic multiplier factor
    Aux_par0 = par(up1(1))                  !Aux heater/boiler parasitic polynomial constant
    Aux_par1 = par(up1(1))                  !Aux heater/boiler parasitic polynomial linear term
    Aux_par2 = par(up1(1))                  !Aux heater/boiler parasitic polynomial quadratic term
    T_startup = par(up1(1))+273.15          ![K] Startup temperature
    
    !Logic and control parameters, cont...
    do i=1,18
        TSLOGIC(i)=PAR(up1(1))              ![-] Values of TSLOGIC array.. intermediate
    enddo
    do i=1,6
        FFRAC(i)=PAR(up1(1))                ![-] Fossil fill fraction control array
    enddo


    !Assign temporary logic array to appropriate structures
    do i=1,int(NUMTOU)
        TSLOGIC_A(i)=TSLOGIC(1+(i-1)*3)
        TSLOGIC_B(i)=TSLOGIC(2+(i-1)*3)
        !Limit turbine operation to max fraction
        TSLOGIC_C(i)=dmax1(dmin1(TSLOGIC(3+(i-1)*3),cycle_max_fraction),cycle_cutoff_frac)
    enddo

    !On the initial timestep call, read in the fluid property file
    call readFluidPropFile(LU_FL)    

    !Calculate heat exchanger size based on parameters
    !hx_ua will be used as the heat exchanger size for simulation calculations
    duty = (solarm - 1.0)*q_pb_design
    call hx_size(hx_config,duty,dT_hot,dT_cold,T_field_out_des,T_field_in_des,field_fluid,store_fluid,hx_eff_des,hx_ua,msghx1)
    

    !---- Initialize outputs
    out(1:nout) = 0.d0
    info(6)=nout


    !---- Set storage array size 
    CALL setStorageSize(nS,INFO)
    !Set initial storage values
    stored(1) = V_tank_hot_ini      ![m3] Initial hot tank HTF volume
    stored(2) = T_tank_hot_ini      ![K] Initial hot tank HTF temperature
    stored(3) = V_tank_cold_ini     ![m3] Initial cold tank HTF volume
    stored(4) = T_tank_cold_ini     ![K] Initial cold tank HTF temperature
    stored(5) = 1.
    stored(6) = V_tank_hot_ini*density(store_fluid,T_tank_hot_ini,1.d0) ![kg] Initial hot tank htf mass
    stored(7) = V_tank_cold_ini*density(store_fluid,T_tank_cold_ini,1.d0) ![kg] Initial cold tank htf mass
    stored(8) = T_field_in_des      ![K] last field outlet temperature
    call setStorageVars(stored,nS,info)

    return 1
endif



!******************************************************************************************************************************
!               Time-dependent conditions
!******************************************************************************************************************************
I_bn = xin(1)/3.6                   ![W/m2] Beam normal radiation (input kJ/m2-hr)
TOUPeriod = xin(2)                  !Time of use period
m_dot_field = xin(3)/3600.          ![kg/s]  Flow rate from the field
m_dot_htf_ref = xin(4)/3600.        ![kg/s] Reference HTF flow rate at design conditions
T_field_out = xin(5)+273.15         ![K] HTF Temperature from the field
T_pb_out = xin(6)+273.15            ![K] Fluid temperature from the power block
T_amb = xin(7)+273.15               ![K] Ambient temperature
m_pb_demand = xin(8)/3600.          ![kg/s] Demand htf flow from the power block
q_startup = xin(9)                  ![MWt] Startup energy reported by the collector field

!******* Read in stored variables every timestep*******
call getStorageVars(stored,nS,info)

V_tank_hot0 = stored(1)             ![m3] Previous total hot tank HTF volume
T_tank_hot0 = stored(2)             ![K] Previous hot tank HTF temperature
V_tank_cold0 = stored(3)            ![m3] Previous total cold tank HTF volume
T_tank_cold0 = stored(4)            ![K] Previous cold tank HTF temperature
stored_mode = int(stored(5))        ![] operation mode
m_tank_hot0 = stored(6)             ![kg] total mass in the hot tank
m_tank_cold0 = stored(7)            ![kg] total mass in the cold tank
T_field_out0 = stored(8)            ![K] Field outlet temp at last timestep


!******************************************************************************************************************************
!               Timestep calculations
!******************************************************************************************************************************
if(I_bn.gt.1.e-6) then   !The solar resource is available above 0. W/m^2
    TSELECT=TSLOGIC_B
else                  !Use the dispatch for no solar resource
    TSELECT=TSLOGIC_A
endif


!-Account for any available thermal storage
itert = 1
T_tank_hotX = T_tank_hot0
T_tank_coldX = T_tank_cold0
30 continue !iteration loop for storage
itert = itert+1

if(TSHOURS > 0.01) then
    !We know the total mass in the tanks, but need to calculate available mass
    !which depends on fluid temperature and minimum volume requirements
    rho_tank_hotX = density(store_fluid,T_tank_hotX,1.d0)
    rho_tank_coldX = density(store_fluid,T_tank_coldX,1.d0)
    V_tank_hot0 = m_tank_hot0/rho_tank_hotX                                 ![m3] The total fluid volume in the hot tank, including unusable volume
    V_tank_cold0 = m_tank_cold0/rho_tank_coldX                              ![m3] The total fluid volume in the cold tank, including unusable volume
    V_tank_hotA = dmax1(V_tank_hot0 - vol_tank*h_tank_min/h_tank,0.d0)      ![m3] The available fluid volume in the hot tank
    V_tank_coldA = dmax1(V_tank_cold0 - vol_tank*h_tank_min/h_tank,0.d0)    ![m3] The available fluid volume in the cold tank
    
    !Calculate the available mass flow rate from the hot and cold tanks
    m_tank_disch_avail = dmax1((V_tank_hotA - V_tank_active*TSELECT(TOUperiod))*rho_tank_hotX/dt, 0.d0)         ![kg/s] Available mass flow rate during discharge cycle
    m_tank_charge_avail = dmax1(V_tank_coldA*rho_tank_coldX/dt,0.d0)      ![kg/s] Available mass flow rate during charge cycle
    
    !-Adjust the available storage to correct for the intermediate heat exchanger
    if(is_hx) then
        T_tank_hot_out = T_tank_hotX
        T_tank_cold_out = T_tank_coldX
        !Calculate the field-equivalent mass flow rate produced if all the storage is dispatched.
        !Discharge cycle.. hot flow comes from storage. Cold flow comes from PB outlet
        if(m_tank_disch_avail > 0.) then
            call hx_perf(hx_ua, hx_config, store_fluid, field_fluid, m_tank_disch_avail, T_tank_hot_out, &
                         T_pb_out, dT_hot, dT_cold,hx_eff, T_tank_cold_in, Ts_hot, q_hx, ms_disch_avail, msghx2)
        else
            ms_disch_avail = 0.
        endif
        !Calculate the field-equivalent mass flow rate produced if all the cold storage is charged. 
        !Charging cycle.. hot flow comes from field. Cold flow comes from storage.
        if(m_tank_charge_avail > 0.) then
            call hx_reverse(hx_ua, hx_config, field_fluid, store_fluid, m_tank_charge_avail, T_field_out, T_tank_cold_out, dT_hot, dT_cold,&
                            hx_eff, Ts_cold, T_tank_hot_in, q_hx, ms_charge_avail, msghx3)
            !Don't allow message triggering for this hypothetical calculation. We really only care about the performance call.
            msghx3=0.d0
        else
            ms_charge_avail = 0.
        endif

    else
        T_tank_hot_out = T_tank_hotX
        T_tank_cold_out = T_tank_coldX
        ms_disch_avail = m_tank_disch_avail         ![kg/s] Available mass flow rate during discharge cycle
        ms_charge_avail = m_tank_charge_avail      ![kg/s] Available mass flow rate during charge cycle
        Ts_hot = T_tank_hot_out
        Ts_cold = T_tank_cold_out
    endif
else
    ms_disch_avail = 0.
    ms_charge_avail = 0.
endif

!-Account for any available auxiliary heater flow
c_htf_aux = specheat(field_fluid,(T_field_out + T_pb_out)/2.,1.d0)*1000. ![J/kg-K] average specific heat
if(ffrac(TOUPeriod) > 0.) then
    m_dot_aux_avail = dmin1(q_max_aux,ffrac(TOUPeriod)*q_pb_design)/(c_htf_aux*(T_field_out - T_pb_out))
else
    m_dot_aux_avail = 0.
endif

!Calculate the maximum potentially available mass flow
m_avail_tot = ms_disch_avail + m_dot_aux_avail + m_dot_field

!Calculate the demanded values from the power block
q_pb_demand = q_pb_design*TSLOGIC_C(TOUperiod) !Adjust the power block demand energy for the TOU period 

!Calculate the maximum potentially available heat flows.  Flows are relative to the
!power block outlet temperature, except in the case of the storage charge level and
!the field inlet temperature. This is because the field inlet temp is a function of
!the calculated values of the other temperatures and flow rates, and because mass 
!flow in this case is calculated relative to the field-side cold outlet temperature.

!A small degree of iteration is needed to get the field inlet temperature to converge
!in the case where storage is charging and an intermediate heat exchanger is used.
!For cases with iteration, set the initial guess values
if(m_avail_tot > 1.e-3) then
    T_field_inX = T_pb_out
else
    T_field_inX = T_field_out
endif
err = 999.; tol = 0.0001; iter = 0; derr=999.; errx = 9999.; defocusX=1.
50 continue  !iterative loop header
iter=iter+1
T_field_in = T_field_inX

!Calculate available heat flows
c_htf_dis = specheat(field_fluid,((Ts_hot + T_field_in_des)/2.d0),1.d0)*1000.  !Specific heat for discharge flow  T_pb_out
c_htf_chg = specheat(field_fluid,((T_field_out + Ts_cold)/2.d0),1.d0)*1000.  !Specific heat for charge flow
c_htf_fld = specheat(field_fluid,((T_field_out + T_field_in)/2.d0),1.d0)*1000.  !Specific heat for field flow


qs_disch_avail = dmax1(ms_disch_avail*c_htf_dis*(Ts_hot - T_pb_out),0.d0)  
qs_charge_avail = dmax1(ms_charge_avail*c_htf_chg*(T_field_out - Ts_cold),0.d0) 
qaux_avail = dmax1(m_dot_aux_avail*c_htf_aux*(T_set_aux - T_pb_out),0.d0) 
qfield_avail = dmax1(m_dot_field*c_htf_fld*(T_field_out - T_field_in),0.d0) 

!Initial values
q_int = 0.d0    !Intermediate equivalent power output from the various heat sources
T_int = T_pb_out !Set intermediate temperature to outlet temperature for now
m_int = 0.d0    !Intermediate mass flow from the various heat sources

!---Allocate the flows from the various sources

!Start for the case where the energy available from the field is less than  
!or equal to the design requirement
if(qfield_avail <= q_pb_demand) then
    !First add the available energy from the field
    q_int = qfield_avail  
    m_int = m_dot_field
    !If energy is coming from the field, adjust the intermediate temperature
    if(q_int > 0.) T_int = T_field_out
    !Next, try to add energy from storage
    if((qfield_avail + qs_disch_avail) < q_pb_demand) then
        !If the total amount of energy is less than the design condition, 
        !add all of the remaining energy from storage
        q_int = q_int + qs_disch_avail
        ms_disch = ms_disch_avail
        ms_charge = 0.
        m_int = m_dot_field + ms_disch
        !Adjust the intermediate temperature, if needed
        if(q_int > 0.) T_int = (T_int*m_dot_field + Ts_hot*ms_disch)/(m_int)

        !Lastly, add any available auxiliary heat to get to the design point
        if((q_int + qaux_avail) < q_pb_demand) then
            q_int = q_int + qaux_avail
            m_dot_aux = m_dot_aux_avail
            m_int = m_int + m_dot_aux
            !Adjust the intermediate temperature, if needed
            if(q_int > 0.) T_int = (T_int*(m_dot_field + ms_disch) + T_set_aux*m_dot_aux)/m_int
        else
            !The aux heater was able to produce more energy than what was needed
            q_aux = q_pb_demand - q_int
            q_int = q_pb_demand
            !Set the required mass flow rate
            m_dot_aux = q_aux/(c_htf_aux*(T_set_aux - T_pb_out))  
            m_int = m_int + m_dot_aux
            !Adjust the intermediate temperature
            T_int = (m_dot_field*T_field_out + ms_disch*Ts_hot + m_dot_aux*T_set_aux)/m_int
        endif

    else
        !Between the energy available from the field and from storage, the dispatch requirement
        !has been met.  No aux load will be added in this condition
        q_int = q_pb_demand
        ms_disch = (q_pb_demand - qfield_avail)/(c_htf_dis * (Ts_hot - T_pb_out))  
        m_int = m_int + ms_disch
        ms_charge = 0.
        m_dot_aux = 0.
        T_int = (m_dot_field*T_field_out + ms_disch*Ts_hot)/m_int
    endif

    !Run a check to see if the energy produced is above the cycle cutout fraction
    if(q_int < q_pb_design*cycle_cutoff_frac) then
        mode = 1
    else
        mode = 2
    endif
    
else 
    !The flow from the field is sufficient to meet the power block demand. If the system has storage
    !the excess flow will be used to charge it.
    
    !Check if the excess flow can be dropped into storage
    if((qfield_avail - q_pb_demand) > qs_charge_avail) then
        mode = 3
    else
        mode = 4
    endif
    
endif

!Check to make sure mode iteration is not happening
if(info(7)>15) mode = stored_mode

!Do mode-specific calculations here
select case(mode)
case(1)
    !If below the cutout fraction, check to see if extra energy can be dumped into storage
    ms_charge = dmax1(dmin1(ms_charge_avail, m_dot_field),0.d0)
    ms_disch = 0.d0
    m_dot_pb = 0.d0
    T_pb_in = T_pb_out
    if(tshours > 0.) then
        if(m_dot_field > 0.) then
            defocus = dmax1(0.d0,ms_charge/dmax1(m_dot_field,1.e-6)*0.95)  !Defocus according the the amount of charge accepted as compared to what is produced
        else
            defocus = 1.d0
        endif
    else
        !Here we need to allow the field to run normally and artificially force the power block to dump
        !energy.  
        defocus = 1.d0
        m_dot_field = 0.d0
    endif
    cycle_pl_control = 2. !1=demand, 2=part load
    standby_control = 3. !!1=normal operation, 2=standby, 3=shut down
case(2)
    !Operation below design point. Continue with calculated operation values
    defocus = 1.            !No defocusing in this case
    m_dot_pb = m_int        !The mass flow rate to the power cycle
    T_pb_in = T_int         ![K] The HTF temperature to the power block
    cycle_pl_control = 2.   !1=demand, 2=part load
    standby_control = 1.    !1=normal operation, 2=standby, 3=shut down
case(3)
    !More thermal energy is available than what can be put into storage
    ms_charge = qs_charge_avail/(c_htf_fld*(T_field_out - Ts_cold)) 
    !Defocus the remaining energy
    select case(iter)
    case(0:5)
        defocus = .92*(qs_charge_avail + q_pb_demand)/qfield_avail
    case(6:10)
        defocus = .80*(qs_charge_avail + q_pb_demand)/qfield_avail
    case(11:)
        defocus = .65*(qs_charge_avail + q_pb_demand)/qfield_avail
    end select
    
    m_dot_pb = dmax1(q_pb_demand/(c_htf_fld*(T_field_out - T_pb_out)),0.d0)
    T_pb_in = T_field_out
    ms_disch = 0.
    cycle_pl_control = 2. !1=demand, 2=part load
    standby_control = 1. !1=normal operation, 2=standby, 3=shut down
    m_dot_aux = 0.
case(4)
    !More energy is available than can be used in the PB, but all of the extra thermal energy can be dropped into storage
    ms_charge = dmax1((qfield_avail - q_pb_demand)/(c_htf_fld*dmax1(T_field_out - Ts_cold,1.e-6)),0.d0)  
    defocus = 1.0
    m_dot_pb = dmax1(q_pb_demand/(c_htf_fld*dmax1(T_field_out - T_pb_out,1.e-6)),0.d0)
    T_pb_in = T_field_out
    ms_disch = 0.
    cycle_pl_control = 2. !1=demand, 2=part load
    standby_control = 1.  !1=normal operation, 2=standby, 3=shut down
    m_dot_aux = 0.
end select
!Set mode here
stored(5) = mode
call setStorageVars(stored,nS,info)


!Calculate the new field inlet temperature 
if((m_dot_pb + ms_charge) > 1.e-3) then
    T_field_inX = (m_dot_pb*T_pb_out + ms_charge*Ts_cold)/(m_dot_pb + ms_charge)
else
    T_field_inX = dmin1(T_field_out,T_startup)
endif

!convergence error is based on actual field mass flow rate compared to balance of mass flow rates
err = sqrt(abs(((ms_charge + m_dot_pb - ms_disch) - m_dot_field)/dmax1(m_dot_field,1.e-6))**2 + ((T_field_in - T_field_inX)/T_field_in)**2)
derr = abs((err - errX)/errX)
errX = err

!Limit number of iterations
if(iter > 20) then
    !Write the error message in preparation for failure, but don't call the printer yet.. 
    write(message(3),fmt=200) iter,dmin1(err,999.d0),tol
    200 format("After ",I2," iterations, the controller convergence error of ",E8.3," could not resolve below the tolerance of ",F6.4)
    msgiter=3.
    !continue on, don't iterate further
else
    !If the defocus control has changed during the calculations, bypass iteration and recalculate field output.
    !Otherwise, the iteration will continue until the max # has been reached.
    if((derr < 0.001) .and. (defocus /= defocusX)) then
        defocusX = defocus
        continue
    else
        if(err > tol) goto 50  !Iterate if error is too high
        !no errors occured, convergence acheived. Move on
        msgiter=0.
    endif
endif

if(is_hx) then
    !-----Calculate flows through the heat exchanger
    !Charging cycle: Collector is hot side - Tanks are cold side
    if(ms_charge > 0.001) then
        call hx_perf(hx_ua, hx_config, field_fluid, store_fluid, ms_charge, T_field_out, T_tank_cold_out, &
                     dT_hot, dT_cold, hx_eff, Ts_cold, T_tank_hot_in, q_hx, m_tank_charge, msghx2)
        m_tank_disch = 0.
    !Discharging cycle: Tanks are hot side - Collector loop is cold side
    elseif(ms_disch > 0.001) then
        call hx_reverse(hx_ua, hx_config,  store_fluid, field_fluid, ms_disch, T_tank_hot_out, T_pb_out, &
                        dT_hot, dT_cold, hx_eff, T_tank_cold_in, Ts_hot, q_hx, m_tank_disch, msghx3)
        m_tank_charge = 0.
    else
        !No flow across the storage tank
        hx_eff = 1.
        Ts_cold = T_pb_out
        Ts_hot = T_field_out
        T_tank_cold_in = T_field_in
        T_tank_hot_in = T_field_out
        m_tank_charge = 0. 
        m_tank_disch = 0.   
    endif
else
    !No heat exchanger
    hx_eff = 0.
    T_tank_cold_in = T_field_in
    T_tank_hot_in = T_field_out
    Ts_cold = T_tank_cold_out
    Ts_hot = T_tank_hot_out
    m_tank_charge = ms_charge
    m_tank_disch = ms_disch
endif

!-----Update the tanks
if(tshours > 0.) then
    !hot tank
    call mixed_tank(& !Inputs and parameters
                          vol_tank,h_tank,h_tank_min,u_tank,tank_pairs,tank_htr_set_point,tank_max_heat,eta_heater_tank,&
                          m_tank_hot0,T_tank_hot0, m_tank_charge,m_tank_disch,T_tank_hot_in,T_amb,store_fluid,&
                            !outputs
                          T_tank_hot_ave,vol_hot_tank_ave,q_loss_hot_tank,T_tank_hotf,&
                          v_tank_hotf,m_tank_hotf,q_tank_hot_htr)
    !Cold tank
    call mixed_tank(& !Inputs and parameters
                          vol_tank,h_tank,h_tank_min,u_tank,tank_pairs,tank_htr_set_point,tank_max_heat,eta_heater_tank,&
                          m_tank_cold0,T_tank_cold0, m_tank_disch,m_tank_charge,T_tank_cold_in,T_amb,store_fluid,&
                            !outputs
                          T_tank_cold_ave,vol_cold_tank_ave,q_loss_cold_tank,T_tank_coldf,&
                          v_tank_coldf,m_tank_coldf,q_tank_cold_htr)

    errt = sqrt(((T_tank_hotX - T_tank_hot_ave)/dmax1(T_tank_hot_ave,1.e-6))**2 + ((T_tank_coldX - T_tank_cold_ave)/dmax1(T_tank_cold_ave,1.e-6))**2)
    T_tank_hotX = T_tank_hot_ave
    T_tank_coldX = T_tank_cold_ave
    if((errt>1.e-5).and.(itert<10)) goto 30 !Iterate to get tank temperature right
    
else
    T_tank_hot_ave = 0. ; vol_hot_tank_ave = 0. ; q_loss_hot_tank = 0. ; T_tank_hotf = 0. ; v_tank_hotf = 0. ; q_tank_hot_htr = 0.
    T_tank_cold_ave = 0. ; vol_cold_tank_ave = 0. ; q_loss_cold_tank = 0. ; T_tank_coldf = 0. ; v_tank_coldf = 0. ; q_tank_cold_htr = 0.
endif 

!-----Update the aux heater
if(qaux_avail>0.) then
                    !Inputs                                                         | Outputs
    call aux_heater(field_fluid, q_max_aux, LHV_eff, T_set_aux, T_pb_out, m_dot_aux, T_aux_act, q_aux_fluid, q_aux_fuel)
else
    T_aux_act = 0. ; q_aux_fluid = 0. ; q_aux_fuel = 0.
endif

!Store the active mode
stored(5) = mode
call setStorageVars(stored,nS,info)

100 continue

!******************************************************************
! Calculate final output values
!******************************************************************
!Pumping power
htf_pump_power = (TES_pump_coef*abs(m_tank_disch-m_tank_charge) + PB_pump_coef*(abs(ms_disch-ms_charge) + m_dot_pb))/1000.  !MW
!Variable parasitic power
q_pb = m_dot_pb*(c_htf_fld*(T_pb_in - T_pb_out)) 
if(q_pb>0.) then
    BOP_par = W_pb_design*BOP_parVal*BOP_parPF*dmin1(BOP_par0 + BOP_par1*(q_pb/q_pb_design)+BOP_par2*(q_pb/q_pb_design)**2,cycle_max_fraction)/1.e6 !MW
else
    BOP_par = 0.d0
endif
!Aux load parasitic power
if(m_dot_aux > 0.) then
    Aux_par = W_pb_design*Aux_parVal*Aux_parPF*dmin1(Aux_par0 + Aux_par1*(q_aux_fluid/q_pb_design)+Aux_par2*(q_aux_fluid/q_pb_design)**2,cycle_max_fraction)/1.e6 !MW
else
    Aux_par = 0.d0
endif

!outputs
out(1) = defocus 
out(2) = cycle_pl_control 
out(3) = standby_control 
out(4) = m_dot_pb*3600.             !demand var - [kg/hr]  Mass flow rate to the power block
out(5) = T_pb_in - 273.15           ![C] HTF Temperature to the power block
out(6) = T_field_in  - 273.15       ![C] HTF temperature into the collector field header
out(7) = (ms_disch-ms_charge)*3600. ![kg/hr] Mass flow rate of the collector side storage loop
out(8) = (m_tank_disch-m_tank_charge)*3600. ![kg/hr] Mass flow rate of the tank side storage loop
out(9) = Ts_hot - 273.15            ![C] Fluid temperature of the hot HTF in the collector loop
out(10) = Ts_cold - 273.15          ![C] Fluid temperature of the cold HTF in the collector loop
out(11) = T_tank_hot_in - 273.15    ![C] Hot tank inlet temperature
out(12) = T_tank_cold_in - 273.15   ![C] Cold tank inlet temperature
out(13) = v_tank_hotf               ![m3] Available Hot tank HTF volume at the end of the timestep
out(14) = T_tank_hotf - 273.15      ![C] Hot tank stored HTF temperature at the end of the timestep
out(15) = v_tank_coldf              ![m3] Available Cold tank HTF volume at the end of the timestep
out(16) = T_tank_coldf - 273.15     ![C] Cold tank stored HTF temperature at the end of the timestep
out(17) = (q_tank_hot_htr + q_tank_cold_htr)  ![MWe] Total parasitic power required for tank freeze protection
out(18) = m_dot_aux*3600.           ![kg/hr] Auxiliary heater mass flow rate
out(19) = q_aux_fluid               ![MWt] Thermal energy provided to the fluid passing through the aux heater
out(20) = q_aux_fuel                ![MMBTU] Heat content of fuel required to provide aux firing
out(21) = v_tank_hotf + v_tank_coldf  ![m3] Total accounted htf volume in storage
out(22) = hx_eff                    ![-] Heat exchanger effectiveness
out(23) = m_tank_hotf               ![kg] mass of total fluid in the hot tank
out(24) = m_tank_coldf              ![kg] mass of total fluid in the cold tank
out(25) = sum(out(23:24))           ![kg] total mass in the storage system
out(26) = htf_pump_power            ![MWe] pumping power for storage, power block loops
out(27) = BOP_par                   ![MWe] parasitic power as a function of power block load
out(28) = W_pb_design*PB_fixed_par/1.e6  ![MWe] Fixed parasitic power losses.. for every hour of operation
out(29) = Aux_par                   ![MWe] Parasitic power associated with operation of the aux boiler
out(30) = dmax1(q_pb/1.e6,0.d0)     ![MWt] Thermal energy to the power block
out(31) = (q_loss_cold_tank+q_loss_hot_tank) ![MWt] Tank thermal losses

return 1

end subroutine


!Include all this stuff because we need to know the condition of the storage tanks and the 
!heat exchanger during the control decision process.  This will keep the control strategy from
!calling and iterating needlessly. Store the ending tank volume levels for the next timestep, 
!and the available storage level will always be known. The aux heater might also need to be included
!in this component.
!The rankine cycle component might need to be run in supply mode, with the controller trying to get
!as close to design as possible based on the calculations.

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

!incrementing function: when called, this function increments one integer value.
!if called with an argument, the counter resets to 0. 
integer function up1(i)
    implicit none
    integer::j,i ; save j
    j=j+1
    if(i==0) j = 0
    up1=j
end function


!Subroutine for storage tanks based on type230
!**************************************************************************************************
!**************************************************************************************************
subroutine mixed_tank(& !Inputs and parameters
                      vol,h,h_min,u,tank_pairs,htr_set_point,max_heat_rate,eta_heater,&
                      m0,T0,m_dot_in,m_dot_out,T_in,T_amb,fluid,&
                        !outputs
                      T_ave,vol_ave,q_loss,T_fin,vol_fin,m_fin,q_heater)
use TrnsysFunctions

!----------Inputs and parameters-------------------------------------------------------------------
!   * vol       - [m3] Total tank volume, including the unusable volume of HTF at the bottom
!   * h         - [m] Total height of the tank (height of HTF when tank is full)
!   * u         - [W/m2-K] Loss coefficient from the tank
!   * tank_pairs- [-] Number of equivalent tank pairs
!   * htr_set_point - [K] Minimum allowable fluid temperature before auxiliary heater turns on
!   * eta_heater- [-] Conversion efficiency for electric heater
!   * m0        - [m3] total HTF volume level at the end of the last time step
!   * T0        - [K] Temperature of the HTF at the end of the last time step
!   * m_dot_in  - [kg/s] mass flow rate of HTF into the tank
!   * m_dot_out - [kg/s] mass flow rate leaving the tank
!   * T_in      - [K] Temperature of the HTF entering the tank
!   * T_amb     - [K] Temperature of the ambient air
!----------Outputs---------------------------------------------------------------------------------
!   * T_ave     - [K] Average HTF temperature throughout the timestep
!   * vol_ave   - [m3] Average HTF volume level during the timestep
!   * q_loss    - [MWt] Total thermal loss rate from the tank
!   * T_fin     - [K] Temperature of the HTF at the end of the timestep
!   * vol_fin   - [m3] Volume of the HTF at the end of the timestep (total volume)
!   * m_fin     - [kg] total mass at the end of the timestep
!   * q_heater  - [MWt] Total energy consumed by the freeze protection heater
!--------------------------------------------------------------------------------------------------
implicit none

!Inputs
real(8),intent(in):: vol,h,h_min,u,tank_pairs,htr_set_point,max_heat_rate,eta_heater,&
                      m0,T0,m_dot_in,m_dot_out,T_in,T_amb,fluid
!Outputs
real(8),intent(out):: T_ave,vol_ave,q_loss,T_fin,vol_fin,m_fin,q_heater

!Locals
real(8):: pi, dt, acs, vol_min, vol_max, dia, ua, density, rho, m_ave, vol0, &
          h_ave, cp, specheat, B, C, D, CC, DD, AA, BB, G, H1, A1, E, &
          Q_vol, Q_flow
!--------------------------------------------------------------------------------------------------
pi=3.1415926
dt= getSimulationTimeStep()*3600.  ![s]

!Geometrical calculations
vol_min = vol*h_min/h       ![m3] Minimum HTF volume among all tanks
vol_max = vol*(1.-h_min/h)  ![m3] Maximum HTF volume among all tanks
acs = vol/(h*tank_pairs)    ![m2] Cross-sectional area of a single tank
dia = (acs/pi)**0.5*2.      ![m] The circumference of a single tank

!Calculate heat loss coefficient
ua = u*(acs + pi*dia*h)*tank_pairs  !u [W/m2-K]

!Calculate conditions at last time step
rho = density(fluid,T0,1.d0)  !density
cp = specheat(fluid,T0,1.d0)*1000.  ![J/kg-K] specific heat
!m0 = ma0 + vol_min*rho              ![kg] mass of total HTF, Add back on the portion of discarded volume 
vol0 = m0/rho                       ![m3] volume of total htf

!Calculate ending volume levels
m_fin = m0 + dt*(m_dot_in - m_dot_out)  ![kg] Available mass at the end of the timestep
m_ave = (m0 + m_fin)/2.d0               ![kg] Average mass 
vol_fin = m_fin/rho                    ![m3] Available volume at the end of the timestep
vol_ave = m_ave/rho                    ![m3] Average volume
h_ave = vol_ave/acs                     ![m] Average HTF height

!Check for no flow
if((dabs(m_dot_in-m_dot_out) < ((m_dot_in + ua/cp)/1.e5)).or.&
   ((m_dot_in < 0.001).and.(m_dot_out < 0.001))) then
    
    !Equations for no flow or zero net flow
    B = m_dot_in + ua/cp
    D = m_dot_in*T_in + (UA/cp)*T_amb
    G = -B/m0
    H1 = 1.0/(dt*(-B))
    A1 = D - B*T0
    E = A1 * exp(dt*G)
    T_fin = (E-D)/(-B)
    T_ave = H1*((E-A1)/G)+D/B
else
    !Equations for unbalanced flow
    B = m_dot_in + ua/cp
    C = m_dot_in - m_dot_out
    D = m_dot_in*T_in + ua/cp*T_amb
    CC = T0 - D/B
    DD = (1.+(C*dt)/m0)**(-B/C)
    T_fin = CC*DD+D/B
    AA = (T0 - D/B)/(C-B)
    BB = (1.+(C*dt)/m0)**(1.-B/C)
    T_ave = AA*(m0/dt)*(BB-1.0) + D/B
endif

!If the temperature of the fluid in the tank is less than the 
!..aux heater set point, heat the fluid in - and passing through - 
!..the tank. MJW
if(T_fin.lt.htr_set_point) then
    Q_vol=cp*vol_fin*rho/(dt*3600.)*(htr_set_point-T_fin)/1000.  !MW
    Q_flow=cp*m_dot_out/3600.*dt/(dt*3600.)*(htr_set_point-T_fin)/1000.  !MW
    Q_heater = dmin1(Q_flow+Q_vol,max_heat_rate) !MW
    T_fin=T0+dt*dmin1(Q_vol*1.e6,max_heat_rate)/(cp*rho*vol_fin)
    T_ave=(T_fin+T0)/2.
else
    Q_heater = 0.
endif

!Calculate losses
q_loss = UA*(T_ave - T_amb)/1.e6    ![MW]


end subroutine

!**************************************************************************************************
!**************************************************************************************************
!Subroutine for sizing a heat exchanger
subroutine hx_size(config,duty,dT_hot,dT_cold,T_h_in,T_h_out,fluid_h,fluid_c,eta,UA,mesg)

use global_props
implicit none

!This subroutine provides the design size for the collector field to storage heat 
!exchanger.
!Inputs required:
!--------------------------------------------------------------------------------------------------
!   * config    - Heat exchanger configuration  {1=Same fluid, 2=Counter-flow, 3=Parallel-flow, 4=Cross-flow, 5=shell-and-tube}
!   * duty      - The heat exchanger duty (MWt)
!   * dT_hot    - Hot side approach temperature [K]
!   * dT_cold   - Cold side approach temperature [K]
!   * T_h_in    - Design hot side inlet temperature [K]
!   * T_h_out   - Design hot side outlet temperature [K]
!   * fluid_h   - Hot side fluid material number [-]
!   * fluid_c   - Cold side fluid material number [-]
!Outputs:
!--------------------------------------------------------------------------------------------------
!   * eta       - Heat exchanger effectiveness at design [-]
!   * UA        - Heat exchanger size [W/K]
!--------------------------------------------------------------------------------------------------

!inputs
real(8),intent(in)::duty,dT_hot,dT_cold,T_h_in,T_h_out,fluid_h,fluid_c
integer,intent(in)::config
!outputs
real(8),intent(out)::eta, UA, mesg
!locals
real(8):: T_ave, specheat, c_c, c_h, T_c_in, T_c_out, m_dot_h, m_dot_c, q_trans,&
          c_dot_max, c_dot_min, c_dot_h, c_dot_c, CR, q_max, NTU, EE, EE1, FF, xx
mesg=0.d0

q_trans = duty     ![W] heat exchanger duty

T_ave = (T_h_in + T_h_out)/2.

!Fluid properties
c_h = specheat(fluid_h,T_ave,1.d0)*1000.
c_c = specheat(fluid_c,T_ave,1.d0)*1000.
!Required temperatures
T_c_out = T_h_in - dT_hot
T_c_in = T_h_out - dT_cold
!Mass flow rates
m_dot_h = q_trans/(c_h*(T_h_in - T_h_out))
m_dot_c = q_trans/(c_c*(T_c_out - T_c_in))
!Capacitance rates
c_dot_h = m_dot_h*c_h
c_dot_c = m_dot_c*c_c
c_dot_max = dmax1(c_dot_h,c_dot_c)
c_dot_min = dmin1(c_dot_h,c_dot_c)
!Capacitance ratio
CR = c_dot_min/c_dot_max
!Maximum possible energy flow
q_max = c_dot_min*(T_h_in - T_c_in)
!Effectiveness
eta = q_trans/q_max

!Heat exchanger configuration correlations
select case(config)
case(1) !Same-fluid
    NTU = -log(1-eta)
case(2) !counter-flow
    !check for no solution
    xx=(1-eta*CR)/(1-eta)
    if(xx < 1.) then
        mesg=1.
    endif
    if (CR /= 1.) then
        NTU = log((1-eta*CR)/(1-eta))/(1-CR)
    else
        NTU = eta/(1-eta)
    endif
case(3) !Parallel-flow
    xx=1.-eta*(1.+CR)
    if(xx < 1.) then
        mesg=1.
    endif

    NTU = log(1.-eta*(1.+CR))/(1.+CR)
case(4) !Cross-flow, unmixed
    xx=1.+log(1.-eta*CR)/CR
    if(xx <= 0.) then
        mesg=1.
    endif

    NTU = -log(1.+log(1.-eta*CR)/CR)
case(5) !Shell-and-tube (assume 1 pass)
    xx=(EE+1.)/(EE-1.)
    if(xx < 1.) then
        mesg=1.
    endif
    
    FF = (eta*CR-1.)/(eta - 1.)
    EE1 = (FF-1.)/(FF-CR)
    EE = (2.-EE1*(1.+CR))/(EE1*sqrt(1.+CR**2))
    NTU = log((EE+1.)/(EE-1.))/sqrt(1.+CR**2)
end select

!Calculate heat exchanger size
UA = NTU*c_dot_min

end subroutine

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


!Subroutine for heat exchanger performance. Provide the hot side mass flow rate, inlet temps, delta temps.
subroutine hx_perf(UA, config, fluid_h, fluid_c, m_dot_h, T_h_in, T_c_in, dT_hot, dT_cold,&
                 !outputs
                 eta, T_h_out, T_c_out, q_trans, m_dot_c,msg)

implicit none

!This subroutine models the heat exchanger between the field HTF and the storage HTF
!Inputs required:
!--------------------------------------------------------------------------------------------------
!   * UA        -   Heat exchanger conductance.  Calculated in the SAM interface [W/K]
!   * config    -   Configuration.. {1=Same fluid, 2=Counter-flow, 3=Parallel-flow, 4=Cross-flow, 5=shell-and-tube}
!   * fluid_h   -   Hot side fluid number (see type 229)
!   * fluid_c   -   Cold side fluid number (see type 229)
!   * m_dot_h   -   Hot side mass flow rate [kg/s]
!   * T_h_in    -   Hot side inlet temperature [K]
!   * T_c_in    -   Cold side inlet temperature [K]
!   * dT_hot    -   Hot side approach temperature
!   * dT_cold   -   Cold side approach temperature
!--------------------------------------------------------------------------------------------------
!Outputs provided:
!   * eta       -   Heat exchanger effectiveness [-]
!   * T_h_out   -   Hot side outlet temperature [K]
!   * T_c_out   -   Cold side outlet temperature [K]   
!   * q_trans   -   Total heat transferred between fluids [MWt]
!--------------------------------------------------------------------------------------------------

!Inputs
real(8),intent(in):: UA, fluid_h, fluid_c, m_dot_h, T_h_in, T_c_in, dT_hot, dT_cold
integer,intent(in):: config
!Outputs
real(8),intent(out):: eta, T_h_out, T_c_out, q_trans, m_dot_c, msg
!Locals
real(8):: NTU, CR, c_dot_min, c_dot_max, c_h, c_c, c_dot_h, c_dot_c, specheat, T_ave, q_max, eta1
msg=0.d0

!-------------All configurations
!Calculate outlet temps 
T_h_out = T_c_in + dT_cold
T_c_out = T_h_in - dT_hot

!look up fluid properties. Assume average specific heat based on inlet temperatures of alternate streams
T_ave = (T_h_in + T_c_in)/2.
c_h = specheat(fluid_h,T_ave,1.d0)*1000. ![J/kg-K]
c_c = specheat(fluid_c,T_ave,1.d0)*1000. ![J/kg-K]

!Calculate heat flow
q_trans = dmax1(m_dot_h * c_h * (T_h_in - T_h_out),0.d0)

!Calculate cold side flow
m_dot_c = q_trans/(c_c * (T_c_out - T_c_in))

!Calculate capacity rates
c_dot_c = m_dot_c*c_c
c_dot_h = m_dot_h*c_h

!Calculate min/max rates, ratio
c_dot_min = dmin1(c_dot_c,c_dot_h)
c_dot_max = dmax1(c_dot_c,c_dot_h)
CR = c_dot_min/c_dot_max

!Calculate non-dimensional HX size
NTU = UA/C_dot_min

!Calculate the effectiveness of the heat exchanger 
select case(config)
case(1) !Same fluid
    eta = 1. - exp(-NTU)
case(2) !Counter-flow
    if(CR < 1.) then
        eta = (1.-exp(-NTU*(1.-CR)))/(1.-CR*exp(-NTU*(1.-CR)))
    else
        eta = NTU/(1.+NTU)
    endif
case(3) !Parallel-flow
    eta = (1.-exp(-NTU*(1.+CR)))/(1.+CR)
case(4) !Cross-flow (unmixed)
    eta = 1.-exp(NTU**0.22/CR*(exp(-CR*NTU**0.78)-1.))
case(5) !Shell-and-tube (assume 1 pass)
    eta1 = 2.*(1.+CR+sqrt(1.+CR**2)*((1.+exp(-NTU*sqrt(1.+CR**2)))/(1.-exp(-NTU*sqrt(1.+CR**2)))))**(-1)
    eta = ((1.-eta1*CR)/(1.-eta1)-1.)/((1.-eta1*CR)/(1.-eta1)-CR)
end select

if ((eta <= 0.) .or. (eta > 1.)) then
    msg=2.
endif

q_trans = q_trans*1.e-6     ![MWt]

end subroutine      
              
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------

!Subroutine for heat exchanger performance. Provide the COLD side mass flow rate, inlet temps, delta temps.
!The subroutine will calculate the hot side mass flow rate required to match cold conditions.  Given an 
!amount of cold side fluid, this will calculate the corresponding hot side flow.
subroutine hx_reverse(UA, config, fluid_h, fluid_c, m_dot_c, T_h_in, T_c_in, dT_hot, dT_cold,&
                      !outputs
                      eta, T_h_out, T_c_out, q_trans, m_dot_h,msg)

implicit none

!This subroutine models the heat exchanger between the field HTF and the storage HTF
!Inputs required:
!--------------------------------------------------------------------------------------------------
!   * UA        -   Heat exchanger conductance.  Calculated in the SAM interface [W/K]
!   * config    -   Configuration.. {1=Same fluid, 2=Counter-flow, 3=Parallel-flow, 4=Cross-flow, 5=shell-and-tube}
!   * fluid_h   -   Hot side fluid number (see type 229)
!   * fluid_c   -   Cold side fluid number (see type 229)
!   * m_dot_c   -   Cold side mass flow rate [kg/s]
!   * T_h_in    -   Hot side inlet temperature [K]
!   * T_c_in    -   Cold side inlet temperature [K]
!   * dT_hot    -   Hot side approach temperature
!   * dT_cold   -   Cold side approach temperature
!--------------------------------------------------------------------------------------------------
!Outputs provided:
!   * eta       -   Heat exchanger effectiveness [-]
!   * T_h_out   -   Hot side outlet temperature [K]
!   * T_c_out   -   Cold side outlet temperature [K]   
!   * q_trans   -   Total heat transferred between fluids [MWt]
!   * m_dot_h   -   Hot side mass flow rate [kg/s]
!--------------------------------------------------------------------------------------------------

!Inputs
real(8),intent(in):: UA, fluid_h, fluid_c, m_dot_c, T_h_in, T_c_in, dT_hot, dT_cold
integer,intent(in):: config
!Outputs
real(8),intent(out):: eta, T_h_out, T_c_out, q_trans, m_dot_h, msg
!Locals
real(8):: NTU, CR, c_dot_min, c_dot_max, c_h, c_c, c_dot_h, c_dot_c, specheat, T_ave, q_max, eta1
msg=0.d0

!-------------All configurations
!Calculate outlet temps 
T_h_out = T_c_in + dT_cold
T_c_out = T_h_in - dT_hot

!look up fluid properties. Assume average specific heat based on inlet temperatures of alternate streams
T_ave = (T_h_in + T_c_in)/2.
c_h = specheat(fluid_h,T_ave,1.d0)*1000. ![J/kg-K]
c_c = specheat(fluid_c,T_ave,1.d0)*1000. ![J/kg-K]

!Calculate heat flow
q_trans = dmax1(m_dot_c * c_c * (T_c_out - T_c_in),1.e-6)

!Calculate hot side flow
m_dot_h = q_trans/(c_h * (T_h_in - T_h_out))

!Calculate capacity rates
c_dot_c = m_dot_c*c_c
c_dot_h = m_dot_h*c_h

!Calculate min/max rates, ratio
c_dot_min = dmin1(c_dot_c,c_dot_h)
c_dot_max = dmax1(c_dot_c,c_dot_h)
CR = c_dot_min/dmax1(c_dot_max,1.e-6)

!Calculate non-dimensional HX size
NTU = UA/C_dot_min

!Calculate the effectiveness of the heat exchanger 
select case(config)
case(1) !Same fluid
    eta = 1. - exp(-NTU)
case(2) !Counter-flow
    if(CR < 1.) then
        eta = (1.-exp(-NTU*(1.-CR)))/(1.-CR*exp(-NTU*(1.-CR)))
    else
        eta = NTU/(1.+NTU)
    endif
case(3) !Parallel-flow
    eta = (1.-exp(-NTU*(1.+CR)))/(1.+CR)
case(4) !Cross-flow (unmixed)
    eta = 1.-exp(NTU**0.22/CR*(exp(-CR*NTU**0.78)-1.))
case(5) !Shell-and-tube (assume 1 pass)
    eta1 = 2.*(1.+CR+sqrt(1.+CR**2)*((1.+exp(-NTU*sqrt(1.+CR**2)))/(1.-exp(-NTU*sqrt(1.+CR**2)))))**(-1)
    eta = ((1.-eta1*CR)/(1.-eta1)-1.)/((1.-eta1*CR)/(1.-eta1)-CR)
end select

if ((eta <= 0.) .or. (eta > 1.)) then
    msg=2.
endif

q_trans = q_trans*1.e-6     ![MWt]

end subroutine      

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

!Subroutine for auxiliary heater - simplified/modified version of TRNSYS Type6 aux heater
subroutine aux_heater(& !Inputs                                | Outputs
                      fluid, qmax, lhv_eff, T_set, T_in, m_dot, T_out, q_fluid, q_fuel)

!--------------------------------------------------------------------------------------------------
! Model inputs
!   * fluid     - [-] Fluid type
!   * qmax      - [MWt] Maximum heating rate 
!   * lhv_eff   - [-] Lower heating value efficiency of the fuel
!   * T_set     - [K] Fluid outlet temperature set point
!   * T_in      - [K] Inlet fluid temperature
!   * m_dot     - [kg/s] Fluid mass flow rate
!--------------------------------------------------------------------------------------------------
! Model outputs
!   * T_out     - [K] Achieved fluid outlet temperature
!   * q_fluid   - [MW] Thermal energy delivered to the fluid
!   * q_fuel    - [MMBTU] Thermal energy content required in fuel
!--------------------------------------------------------------------------------------------------
use TrnsysFunctions

implicit none

!Inputs
real(8),intent(in):: fluid, qmax, lhv_eff, T_in, m_dot, T_set
!Outputs
real(8),intent(out)::T_out, q_fluid, q_fuel
!Locals
real(8):: specheat, cp, dthr
!--------------------------------------------------------------------------------------------------
dthr = getSimulationTimeStep()

cp = specheat(fluid,T_in,1.d0)*1000. ![J/kg-K]

T_out = dmin1(T_set, (T_in + (qmax*1.e6)/(m_dot*cp)))  ![K]
q_fluid = m_dot*cp*(T_out - T_in)/1.e6          ![MWt]
q_fuel = q_fluid/lhv_eff*dthr*3.41214116        ![MMBTU]

end subroutine
