module HyGSHP_calc

contains
subroutine HyGSHP_sub(in,LCC2)
use global
use opt_global, only: exit_early,exit_early_factor,first_val,logfile
use global_print
implicit none 

    integer :: i,j,its,k,ii,jj    
    real(DP) :: totallength,size_towerkw,tpc,dtpc
    integer,parameter :: max_its = 100
    integer,parameter :: n_items = 5 !number of items in the input file
    real(DP) :: etol,etol2,etol3,x(1000),fx(1000),dfx(1000),delta,delta2,delta3,last,yra,dthp(1000),Tdiff(1000), &
                        Tdiff2(1000),Tdiff3(1000),Tdiff4(1000),etol4,delta4,repeat,&
                        test1(1000), test2(1000)
    parameter tol = 0.001
    real(DP) :: time_gr
    character(len = 4) stop_run
    character(len = 500) :: makedirectory
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !General data variables
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: Maxcoolkjhr,Maxcoolkjhr_test,Maxheatkjhr,Maxheatkjhr_test,rho_fluid
    real(DP) :: cp_fluid,maxmassflow,maxflow,tboil,time_tot,time_step,time_step_day,timeold
    real(DP) :: temp
    real(DP) :: time2
    integer :: tempint
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Heat pump variables
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: QLOAD_C, QLOAD_H, onsig_fan
    real(DP) :: P_AIR_RETurn, DP_AIR_RETurn, T_AIR_RETurn, RH_AIR_RETurn, W_AIR_RETurn
    real(DP) :: P_AIR_OA, DP_AIR_OA, T_AIR_OA, RH_AIR_OA, W_AIR_OA
    integer :: Continuous, mode, status, FLOWMODE
    real(DP) :: T_AIR_IN, wb_air_in
    real(DP) :: T_FLUID_IN,  T_DHW_IN, Frac_OA,  P_AIR_IN
    real(DP) :: slopec, intc, COP_C, slopec_cop, intc_cop, Power_CF_c 
    real(DP) :: FRAC_CAP_C, FRAC_POW_C, LOAD_MODE, CONTPOW
    real(DP) :: slopeh, inth, COP_H, slopeh_cop, inth_cop, Power_CF_h 
    real(DP) :: FRAC_CAP_H, FRAC_POW_H,t_air_ret_c,t_air_ret_h,RetAirTemp 
    !!related to interpolation
    integer :: mc, nc, nh
    real(DP) ::  T_DHW_OUT, POWER_UNIT, Q_DHW, EER
    integer :: iso_loop
    real(DP) :: per_pg, psydat(9)
    real(DP) :: outc(6),outh(6),occup  
    real(DP) :: par_hp(6),xin_hpc(9),xin_hph(9),blowpow,lps_tot,bypassflow
    
    real(DP) :: P_AIR_RET, DP_AIR_RET, T_AIR_RET, RH_AIR_RET, W_AIR_RET
    real(DP) :: P_AIR_OAi, DP_AIR_OAi, T_AIR_OAi, RH_AIR_OAi, W_AIR_OAi
    real(DP) :: eff_furnace,Qsupp,Qheat_loop,Qf_cap !actually necessary for furnace, just easier to put it here
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Valve variables
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: out_mix(2),out_mix2(2),out_mix3(2),out_div(6),out_div2(4),out_div3(4)    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !GHX variables
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: par21,rlamlin,tfchs,qwchs,tachsh,tachs,circul,crin,out_ghx(10),par_ghx(15),&
                xin_ghx(5)       
    real(DP) :: FLOW_SIZED_GHX
    real(DP) :: maxdrilldepth,spacingbores,headerde,centerto,borehole,diff_g_m2day,k_g_w,innerrad,c_g_kJm3K_in
    real(DP) :: k_grout_w,outerrad,k_g_kjhr(1),k_grout_kjhr,c_g_kjm3k(1),numseries,maxnum,num_interm,num
    real(DP) :: storagevolume,boredepth,thicknesslayer,numvert,numrad,refflow
    real(DP) :: Tij(15,15),zdist(15),rdist(15),rdist1(1,15),FRRij(15,15)
    real(DP) :: qgc,qgh
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !CT variables
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    integer :: MODE_AIR, N_COEFfs
    real(DP) :: EWT_IN_DESIGN, EWT_OUT_DESIGN, FLOW_DESIGN, T_AIR_DESIGN
    real(DP) :: T_WB_DESIGN, FLOW_AIR_DESIGN, P_AIR_DESIGN, P_FAN_RATED, Y(4)
    real(DP) :: FLOW_FLUID, W_AIR_IN, RH_AIR_IN, GAMMA_FAN,t_air_in_Ct,w_air_in_ct
    real(DP) :: par_ct(4),xin_ct(7),out_ct(5),mULT_CTFLUIDFL,MULT_CTFANFL,MULT_CTFANP
    real(DP) :: flowdes_c,airdes_c,ccct_ratio_hp,PFANDES_C,spray_gpm,lo_spd,lo_spd_tset
    integer :: nvar,ntim
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !controller
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: xin_cont(32),out_cont(8),tct(1000),tghx(1000)
    real(DP) :: TH1,Lo_SPD_PC
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !pumps
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: par_p1(3),xin_p1(5),out_p1(5),pdroppctghx,pumpkwper,tons_rejection
    real(DP) :: pd_tot_max,ghx_off_fact,pumpflow,pumprating,pd_tot,eff_pump,eff_motor,eff_motor2
    real(DP) :: par_p2(11),xin_p2(4),out_p2(5)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !integrator stuff
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    logical :: first_run_455,first_run_dem,first_run_econ
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Misc
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: gpm_sized_ccct      
    real(DP) :: supplemental_device,qhcorplr,cnstr40_1min,ewt_penalty,tmin_ghp,tmax_ghp
    real(DP) :: GHX_load_total,Imbalance_pen,bldg_load_total,load_pen_fact,temp_pen_fact,ghx_Ratio
    real(DP) :: energy_pen,e1,e2
    real(DP) :: tmin_fluid,sizea,cost_ghxperm,cost_pgperm,cost_inside,tank_ratio_cost_ct
    real(DP) :: ccct_ratio_cost,first_closed10i,mult_fstclosed,cost_pg_bldg,cost_pg_ghx
    real(DP) :: first_gchx,first_sump,first_closed,cost_watperm3,evap_gph_c
    real(DP) :: cycle_conc,blow_gph_c,gallons_water,drift_gph_c,ft3_water,cost_water,cost_trtmnt_c
    
    real(DP) :: par_time(6),out_time(18),summer_first,summer_last,peak_first,peak_last,mode_lcc,nalternatives 
    
    real(DP) :: maint_pg,maint_closed10i,mult_mntclosed,ccct_ratio_costmnt,maint_gchx,maint_closed
    
    real(DP) :: continfan,cost_pkkwh,cost_offpkkwh,cost_kwh_onw,cost_kwh_offw,cost_kw,cost_kw_w
    real(DP) :: ccct_ratio_sp,pow_spray_kw,plr,fanonlyenergy,p_hp_int,pump_power,cost_kwh,cost_gas
    real(DP) :: cost_kw_,elec_kw,cost_elec,pow_spray_kjhr,peak_dur,out_455(100),xin_455(10),par_455(70)

    
    real(DP) :: time0,cost_kwr,cost_in_ann,q_abs,q_rej,q_cl_hp,q_ht_hp,years_analysis
    real(DP) :: cost_ratchet,cost_demand,cost_maint,cost_in_anntot,cost_yearly
    real(DP) :: reset_econ,xin_econ(19),out_econ(19),base_econ,store_econ(3),store_455(251)
    real(DP) :: par_dem(70),xin_dem(10),out_dem(100),store_dem(251),MAINT_CLOSED10,first_CLOSED10
    
    real(DP) :: max_flow_lps,pow_total,num_penalties,temp_annchange,optlimits,datalimits,ewt_hp
    real(DP) :: lwt_hp,lwt_ghx,flow_hp,power_tower,power_pumping,power_hps,xin_lcc(20),par_lcc(91)
    real(DP) :: out_lcc(70),years_loan,comm,fuel_infl,disc_rate,tax_rate,inflation,prop_tax,loan_int_rate
    real(DP) :: depr_life,down_paymt,rebate,salvage_frac,iincr_cool,iincr_heat,iincr_tons,lcc2,lastd,&
                        lastd2,lastd3
    real(DP) :: first_incr_hp_perton,in(10),LCC_pen,elec_total,gas_total,cost_utility
    integer :: fluid_type,nin,nin_2,ni_econ,REflag,EWT_penalty_count
    
    real(DP) :: timeprint,fluid_vol,tgk,CH_ratio,a1,a2,a3,a4,a5
    real(DP), save :: exit_early_factor2
    real(DP),allocatable,save :: store(:,:)    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Printing
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP),allocatable :: time_print(:),T_hp_h(:),F_hp_h(:),T_hp_C(:), &
                        F_hp_C(:),T_ghx(:),F_ghx(:),T_ct(:), &
                        F_ct(:),T_in(:),time_day(:),Ptower(:), &
                        Ppumping(:),Phps(:),T_ct_in(:),T_ghx_in(:)
    real(DP),allocatable :: Q_hp_h(:),Q_hp_c(:),Q_ghx(:),Q_ct(:), &
                        Q_cool(:),Q_heat(:),T_store(:),Y_closed(:), &
                        Y_fan_closed(:),Y_GHX(:),Y_night(:)                
    real(DP),allocatable :: ITS_cum(:),delta_end(:),etol_end(:), &
                        delta2_end(:),etol2_end(:),Tamb(:), &
                        delta3_end(:),etol3_end(:)  
    real(DP),allocatable :: Celec(:),Cwater(:),Cmaint(:),CGHX(:), &
                        Csump(:),CCT(:),Cfirst(:),intfuelcost(:), &
                        p1(:),p2(:),PCon(:),Qsupp_print(:),Qheat_loop_print(:), &
                        Egas(:),Eelect(:),Qratio(:),Twet(:)
    real(DP),allocatable :: Tg_rad1(:)
    real(DP),allocatable :: zprint(:)  
    real(DP),allocatable :: Temp_print(:,:),Flow_print(:,:),Power_print(:,:),Heat_print(:,:),Misc_print(:,:)
    character*33 :: temp_print_lbl(28)
    integer,allocatable :: pv_2(:)
    integer :: sumT,sumF,sumP,sumQ,sumM,sumTot
    real(DP) :: furnace_tot,ghx_tot,qgc_ts,qgh_ts    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Modified controls
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP) :: bypassflow_ct_ghx,fct_pc,pc_dur,pc_start,dt3,spd1,spd2,nspeeds,ct1,ct2,ghx1,m,b
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !trouble shooting
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real(DP),allocatable :: debal(:)
    real(DP),allocatable :: demandp(:),ratchetp(:),cost_utilityp(:),&
                        cost_yearlyp(:),cost_kwhp(:),cost_elecp(:),&
                        elec_kwp(:),y_closedh(:),y_fan_closedh(:),&
                        y_gchxh(:)
    integer :: its2
    real(DP) :: therm_store(2)                                  
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Interpolation for CT cost
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    integer :: nct
    real(DP), allocatable :: xct(:),yct(:),y2ct(:) 
    real(DP) :: yp1ct,ypnct
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Variables for run time
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    integer :: t1,t2,clock_rate,clock_max,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12
    real(DP) :: ti1,junk,tihp,ticont,tict,tighx,tiprint
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Namelist
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    namelist /economics/ mode_lcc,nalternatives,peak_last,peak_first,Summer_First,Summer_Last, &
	                     years_loan,COMM,FUEL_INFL,DISC_RATE,TAX_RATE,INFLATION,PROP_TAX,LOAN_INT_RATE, &
    	                 DEPR_LIFE,DOWN_PAYMT,REBATE,SALVAGE_FRAC,YEARS_ANALYSIS,FIRST_INCR_HP_PERTON, &
	                     COST_GHXperM,COST_PGperM,COST_WATperM3,FIRST_CLOSED10I,MULT_FSTCLOSED,MAINT_PG, &
    	                 MAINT_CLOSED10I,MULT_MNTCLOSED,COST_PKkwh,COST_OFFPKkwh,Cost_kWh_OnW,Cost_kWh_OffW, &
    	                 COST_kW,Cost_kW_W,COST_KWR,COST_IN_ANN,load_pen_fact,temp_pen_fact,cost_gas
    namelist /DST/ MaxDrillDepth,SpacingBores,K_g_W,c_g_kJm3K,tground
    namelist /HP/ Continuous,blowpow,LPS_TOT,t_air_ret_c,t_air_ret_h,TMIN_GHP,TMAX_GHP
    namelist /heat_pumps/ contpow,load_mode,flowmode,slopeh,inth,slopec,intc,slopeh_cop,inth_cop,slopec_cop, &
                          intc_cop,Frac_OA,P_air_OA,FRAC_POW_C,FRAC_CAP_C,FRAC_CAP_H,FRAC_POW_H,Per_PG,mode, &
                          DP_AIR_OAi,cop_c,cop_h,RH_AIR_RET,P_AIR_RET,DP_AIR_RET,w_AIR_RET,T_DHW_IN,P_AIR_IN, &
                          DP_AIR_OA,Flow_fluid
    namelist /CT/ MULT_CTFANP,MULT_CTFANFL,MULT_CTFLUIDFL,Lo_SPD_PC,spd1,spd2,nspeeds
    namelist /Fluid_props/ RHO_FLUID,cp_fluid,FLUID_TYPE,Iso_Loop
    namelist /circ_pump/ PumpKWPer,PdropPctGHX,PD_TOT_MAX,eff_pump,eff_motor,eff_motor2    
    namelist /Furnace_n/ eff_furnace 
    
    open(unit = 1,file = "..\..\namelist_files\component_inputs.txt")
    read(1,nml = CT)
    read(1,nml = HP)
    read(1,nml = DST) 
    read(1,nml = Fluid_props)
    read(1,nml = Furnace_n)
    read(1,nml = circ_pump)
    read(1,nml = economics)
    
    close(1)
    
    open(unit = 1,file = "..\..\namelist_files\op_set.txt")
    read(1,nml = heat_pumps)
    close(1)    
    
   
    ti1 = 0
    tihp = 0
    ticont = 0
    tict = 0
    tighx = 0
    tiprint = 0
    t1 = 0
    t2 = 0
    t3 = 0
    t4 = 0
    t5 = 0
    t6 = 0
    t7 = 0
    t8 = 0
    t9 = 0
    t10 = 0
    t11 = 0
    t12 = 0
    
    call system_clock(t1,clock_rate,clock_max)
    
    nct = 21
    allocate(xct(nct))
    allocate(yct(nct))
    allocate(y2ct(nct))
    xct = (/1,2,5,10,15,20,25,30,40,50,52,75,95,115,130,150,180,210, &
          240,270,300/);
    yct = (/1.1539,1.4356,2.2555,3.538,4.7155,5.788,6.7555,7.618, &
          9.028,10.018,10.1656,11.35,12,12.5,12.75,13,13.2,13.33,13.4,&
          13.45,13.5/);  
    yp1ct = 1e30
    ypnct = 1e30
     
    k = 1
    timeprint = 0
    first_run_455 = .false.
    first_run_dem = .false.
    first_run_econ = .false.
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Read in input data
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    !Note that this code assumes the data file provides hourly data; the hourly data
    !are broken into subhourly data (if step < 1) by simply using the same value for the
    !entire hour
    nin = 4
    nin_2 = 2
      
    time_tot = ilen*yr   !total simulation time
    time_step = ilen*yr/delt !total simulation time considering a time step  
    temp = delt
    allocate(Qratio(int(ilen*yr/delt)+1))

    allocate(time_print(int(ilen*yr/delt)+1))
    allocate(T_hp_h(int(ilen*yr/delt)+1))
    allocate(F_hp_h(int(ilen*yr/delt)+1))
    allocate(T_hp_C(int(ilen*yr/delt)+1))
    allocate(F_hp_C(int(ilen*yr/delt)+1))
    allocate(T_ghx(int(ilen*yr/delt)+1))
    allocate(F_ghx(int(ilen*yr/delt)+1))
    allocate(T_ct(int(ilen*yr/delt)+1))
    allocate(F_ct(int(ilen*yr/delt)+1))
    allocate(T_in(int(ilen*yr/delt)+1))
    allocate(time_day(int(ilen*yr/delt)+1))
    allocate(Ptower(int(ilen*yr/delt)+1))
    allocate(Ppumping(int(ilen*yr/delt)+1))
    allocate(Phps(int(ilen*yr/delt)+1))
    allocate(T_ct_in(int(ilen*yr/delt)+1))
    allocate(T_ghx_in(int(ilen*yr/delt)+1))
    allocate(Q_hp_h(int(ilen*yr/delt)+1))
    allocate(Q_hp_c(int(ilen*yr/delt)+1))
    allocate(Q_ghx(int(ilen*yr/delt)+1))
    allocate(Q_ct(int(ilen*yr/delt)+1))
    allocate(Q_cool(int(ilen*yr/delt)+1))
    allocate(Q_heat(int(ilen*yr/delt)+1))
    allocate(T_store(int(ilen*yr/delt)+1))
    allocate(Y_closed(int(ilen*yr/delt)+1))
    allocate(Y_fan_closed(int(ilen*yr/delt)+1))
    allocate(Y_GHX(int(ilen*yr/delt)+1))
    allocate(Y_night(int(ilen*yr/delt)+1)  )
    allocate(ITS_cum(int(ilen*yr/delt)+1))
    allocate(delta_end(int(ilen*yr/delt)+1))
    allocate(etol_end(int(ilen*yr/delt)+1))
    allocate(delta2_end(int(ilen*yr/delt)+1))
    allocate(etol2_end(int(ilen*yr/delt)+1))
    allocate(Tamb(int(ilen*yr/delt)+1))
    allocate(Twet(int(ilen*yr/delt)+1))
    allocate(delta3_end(int(ilen*yr/delt)+1))    
    allocate(etol3_end(int(ilen*yr/delt)+1))
    allocate(Celec(int(ilen*yr/delt)+1))
    allocate(Cwater(int(ilen*yr/delt)+1))
    allocate(Cmaint(int(ilen*yr/delt)+1))
    allocate(CGHX(int(ilen*yr/delt)+1))
    allocate(Csump(int(ilen*yr/delt)+1))
    allocate(CCT(int(ilen*yr/delt)+1))
    allocate(Cfirst(int(ilen*yr/delt)+1))   
    allocate(intfuelcost(int(ilen*yr/delt)+1)) 
    allocate(p1(int(ilen*yr/delt)+1))
    allocate(p2(int(ilen*yr/delt)+1))
    allocate(Tg_rad1(int(ilen*yr/delt)+1))
    allocate(zprint(15*(int(ilen*yr/delt)+1)) )
    allocate(debal(int(ilen*yr/delt)+1))
    allocate(demandp(int(ilen*yr/delt)+1))
    allocate(ratchetp(int(ilen*yr/delt)+1))
    allocate(cost_utilityp(int(ilen*yr/delt)+1))
    allocate(cost_yearlyp(int(ilen*yr/delt)+1))
    allocate(cost_kwhp(int(ilen*yr/delt)+1))
    allocate(cost_elecp(int(ilen*yr/delt)+1))
    allocate(elec_kwp(int(ilen*yr/delt)+1))
    allocate(y_closedh(int(ilen*yr/delt)+1))
    allocate(y_fan_closedh(int(ilen*yr/delt)+1))
    allocate(y_gchxh(int(ilen*yr/delt)+1))
    allocate(PCon(int(ilen*yr/delt)+1))
    allocate(Qsupp_print(int(ilen*yr/delt)+1))
    allocate(Qheat_loop_print(int(ilen*yr/delt)+1))
    allocate(Egas(int(ilen*yr/delt)+1))
    allocate(Eelect(int(ilen*yr/delt)+1))
    
    Qsupp_print = 0
    T_store = tground
    
    !Only allocate print files if the user has chosen to print certain items
    !temperature data
    sumT = sum(print_vector(1:11))
    sumF = sum(print_vector(12:16))
    sumP = sum(print_vector(17:19))
    sumQ = sum(print_vector(20:27))
    sumM = (print_vector(28))
    sumTot = sumT+sumF+sumP+sumQ+sumM
    if(sumTot > 0)then
        print_stuff = 1
        allocate(pv_2(sumTot))
        jj = 1
        do ii = 1,28
            if(print_vector(ii) > 0)then
                pv_2(jj) = ii
                jj = jj + 1
            endif
        enddo
        
        temp_print_lbl = (/"GHX_Inlet_Temperature_(C)", &
                           "GHX_Outlet_Temperature_(C)", &
                           "GHX_Mean_Storage_Temperature_(C)", &
                           "CT_Inlet_Temperature_(C)", &
                           "CT_Outlet_Temperature_(C)", &
                           "HP_Inlet_Temperature_(C)", &
                           "HP_C_Outlet_Temperature_(C)", &
                           "HP_H_Outlet_Temperature_(C)", &
                           "HP_Outlet_Temperature_(C)", &
                           "Wet_Bulb_Temperature_(C)", &
                           "Dry_Bulb_Temperature_(C)",&
                           "Pump_Flow_Rate_(kg/hr)", &
                           "HP_C_Flow_Rate_(kg/hr)", &
                           "HP_H_Flow_Rate_(kg/hr)", &
                           "CT_Flow_Rate_(kg/hr)", &
                           "GHX_Flow_Rate_(kg/hr)",&
                           "CT_Power_Consumption_(kW)", &
                           "HP_Power_Consumption_(kW)", &
                           "Pump_Power_Consumption_(kW)",&
                           "Heat_from_the_GHX_(kJ/hr)", &
                           "Heat_to_the_GHX_(kJ/hr)", &
                           "Heat_rejected_by_the_CT_(kJ/hr)", &
                           "Heating_by_the_HP_(kJ/hr)", &
                           "Cooling_by_the_HP_(kJ/hr)", &
                           "Heating_by_the_Furnace_(kJ/hr)", &
                           "Cooling_Load_(kJ/hr)", &
                           "Heating_Load_(kJ/hr)", &
                           "CT_Fan_Speed_(-)"/)
    endif
    allocate(Temp_print(int(ilen*yr/delt)+1,11))   
    allocate(Flow_print(int(ilen*yr/delt)+1,5))    
    allocate(Power_print(int(ilen*yr/delt)+1,3))          
    allocate(Heat_print(int(ilen*yr/delt)+1,8))     
    allocate(Misc_print(int(ilen*yr/delt)+1,1))         
                    
    if(sumT > 0)then              
        open(unit = 50,file = trim(file50),recl = 2140000)  
        write(50,*) temp_print_lbl(pv_2(1:sumT))
    endif
    !flow data
    if(sumF > 0)then    
        open(unit = 51,file = trim(file51),recl = 2140000)  
        write(51,*) temp_print_lbl(pv_2(sumT+1:sumF+sumT))
    endif
    !power data
    if(sumP > 0)then    
        open(unit = 52,file = trim(file52),recl = 2140000)  
        write(52,*) temp_print_lbl(pv_2(sumT+sumF+1:sumP+sumF+sumT))
    endif
    !heat data
    if(sumQ > 0)then    
        open(unit = 53,file = trim(file53),recl = 2140000)  
        write(53,*) temp_print_lbl(pv_2(sumT+sumF+sumP+1:sumQ+sumP+sumF+sumT))
    endif
    !miscellaneous
    if(sumM > 0)then     
        open(unit = 54,file = trim(file54),recl = 2140000) 
        write(54,*) temp_print_lbl(pv_2(sumT+sumF+sumP+sumQ+1:sumM+sumQ+sumP+sumF+sumT))
    endif
   
    if(.not.first_run_opt)then
          first_run_opt = .true.
    
    if((print_stuff == 1).and.(opt == 0))then         
        open(unit = 18, file = trim(file3),recl = 2140000)
        write(18,10) "time2","timeday","Thph","Thpc","Tg","Tct","EWT","Tct_in","Tg_in","Tamb","Fhph", &
                     "Fhpc","Fg","Fct","Tg_rad1","Twb"
        10 format(2x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x)
            
        open(unit = 19, file = trim(file4),recl = 2140000)
        write(19,11) "time2","timeday","Qhph","Qhpc","Qg","Qct","Qcool","Qheat","tstore","CT","FanSpeed", &
                     "GHX","Night","Qsupp","Qheat_loop","Heat_ratio","Qratio"
        11 format (2x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x)   
         
        open(unit = 20, file = trim(file5),recl = 2140000)
        write(20,12) "time_tot","time_ghx","time_hp_cool","time_hp_heat","time_ct","time_printer1"
        12 format (2x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x)  
        
        open(unit = 25, file = trim(file7),recl = 2140000)
    
        open(unit = 23, file = trim(file8),recl = 2140000)
        write(23,14) "time2","its","delta","etol","delta2","etol2","delta3","etol3"
        14 format (2x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x)
        
        open(unit = 26, file = trim(file9),recl = 2140000)
        write(26,15) "time2","timeday","Pct","Ppump","Php","EImbalance"
        15 format(2x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x)
        
        open(unit = 27, file = trim(file10),recl = 2140000)
        write(27,16) "time2","timeday","cost_elec","cost_water","cost_maint","first_gchx", &
                     "first_sump","first_closed","first_cost","intfuelcost","p1","p2"
        16 format(2x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x)
        
        open(unit = 28, file = trim(file11),recl = 2140000)
        write(28,17) "depth","temperature"
        17 format(2x,a,5x,a,5x)
        
        open(unit = 29, file = trim(file12),recl = 2140000)
        write(29,18) "time2","timeday","demand","ratchet","cost_utility","cost_yearly",&
                     "cost_kwh","cost_elec","elec_kw"
        18 format(2x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x,a,5x)
        
        open(unit = 30, file = trim(file13),recl = 21400000)
        write(30,17) "depth","flux"    
    endif
    
    open(unit = 16, file = trim(file6),status = "replace",recl = 2140000)
    write(16,13) "LCC2","FIRST_COST","COST_ELEC","ELEC_KW","Size_Tower_kW", &
                 "Qf_cap","CT1","CT2","TOTAL_LENGTH_m","GHX1", &
                 "TPC","DTPC","Temp_AnnChange","FUEL_COST","ANNUAL_FUEL", &
                 "MIN_GCHX","MAX_GCHX","PUMP_POWER","Pfan","Qct", &
                 "POW_SPRAY_KJHR","Qghx","p_hp_int","Q_REJ","Q_ABS", &    
                 "GALLONS_WATER","DRIFT_GPH_C","EVAP_GPH_C","BLOW_GPH_C","FIRST_GCHX", &                 
                 "FIRST_CLOSED","MAINT_GCHX","MAINT_CLOSED","EWT_PENALTY","COST_WATER", &                  
                 "COST_DEMAND","COST_MAINT","PD_TOT","TOTAL_ENERGY_kJ","storagevolume", &
                 "num","Total_PC","Imbalance_pen","Cost_Gas","GHX_temp_change", &
                 "qghx_ratio","Source_elec","Source_gas","Source_total","energy_pen", &                  
                 "furnace_tot(kJ/hr)","ghx_tot(kJ/hr)"    
    13 format (2x,a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x,a,5x,a,5x,a,5x, & 
                  a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x,a,5x,a,5x,a,5x, & 
                  a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x,a,5x,a,5x,a,5x, &
                  a,5x,a,5x)       
                  
!                      13 format (2x,a,5x,a,5x,a,5x,a,5x,a,5x, &
!                  a,5x,a,5x,a,5x,a,5x,a,5x, &
!                  a,5x,a,5x,a,5x,           &
!                  a,5x,a,5x,a,5x,a,5x,a,5x, &
!                  a,5x,a,5x,a,5x,a,5x,a,5x, &
!                  a,5x,a,5x,a,5x,a,5x,a,5x, &
!                  a,5x,a,5x,a,5x,a,5x,      &
!                  a,5x,                     &
!                  a,5x,a,5x,a,5x,a,5x,a,5x, &
!                  a,5x,a,5x,                &
!                  a,5x,a,5x,a,5x,           &
!                  a,5x,a,5x,a,5x)              
  
    endif     
   
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    totallength = in(1)
    size_towerkw = in(2)
    ct1 = in(3)
    ct2 = in(4)
    ghx1 = in(5)
    tpc = in(6)
    dtpc = in(7)
    Qset_ratio = in(8)
    m = 0.838    !in(9)
    b = -0.552   !in(10)
    
    if((totallength > 8550).and.(size_towerkw > 450).and.(ct1 > 40).and.(ct2 > 28).and.(ghx1 > 13).and.(ghx1 < 14))then
        temp = 1
    endif
    
    xin_econ=0
    out_econ=0
    base_econ=0
    store_econ=0
    xin_455 = 0
    out_455 = 0
    store_455=0
    par_dem=0
    xin_dem=0
    out_dem=0
    store_dem=0
    out_lcc=0
    XIN_lcc = 0
    PCon = 0
    
    time_step = 8760*yr/delt    !total simulation time considering time step
    time_tot = 8760*yr
       
    !!!Data needed by several components!!!
    MaxcoolkJhr_test = maxval(Qcool)
    MaxheatkJhr_test = maxval(Qheat)
    
    !write(*,*) MaxcoolkJhr_test
    !write(*,*) MaxheatkJhr_test
    
    qgc = 0
    qgh = 0
    
    MaxCoolkJhr = -MaxcoolkJhr_test  !-5132.6 * 1055.1
    MaxheatkJhr = MaxheatkJhr_test  !7465.1*1055.1
    Qset = Qset_ratio*MaxheatkJhr
    maxflow = max(-0.04463*MaxCoolkJhr,.03513*MaxheatkJhr) 
    TONS_REJECTION = -MaxCoolkJhr*0.08333*0.0009478

    par_hp(1) = RHO_FLUID
    par_hp(2) = CP_FLUID
    par_hp(3) = BLOWPOW
    par_hp(4) = LPS_TOT
    par_hp(5) = MAXFLOW 
                  
    xin_cont(1) = 5 !dead band
    xin_cont(2) = size_towerkW
    xin_cont(3) = lo_spd_pc              
    xin_cont(10) = ct1   
    xin_cont(11) = ct2  
    xin_cont(12) = ghx1   
    xin_cont(13) = -15
    xin_cont(20:22) = 0 !initialize whether or not the device was on in the prior time step
    xin_cont(23) = TPC
    xin_cont(24) = DTPC    
    xin_cont(26) = spd1
    xin_cont(27) = spd2
    xin_cont(28) = nspeeds
    xin_cont(30) = m
    xin_cont(31) = b
    
    PEAK_DUR = Peak_Last - Peak_first
    
    par_455(1) = 1
    par_455(2) = Peak_First
    par_455(3) = Peak_dur
    par_455(4) = 24
    par_455(5) = -1
    par_455(6) = 0
    par_455(7) = 876000    !time_tot
    
    par_455(8) = 1
    par_455(9) = 0
    par_455(10) = 8760
    par_455(11) = 8760
    par_455(12) = 8760
    par_455(13) = 0
    par_455(14) = 876000  !time_tot
    
    par_455(15) = -3
    par_455(16) = -1
    par_455(17) = time_tot
    par_455(18) = time_tot
    par_455(19) = time_tot
    par_455(20) = 0
    par_455(21) = time_tot
    
    par_455(22) = -4
    par_455(23) = -1
    par_455(24) = time_tot
    par_455(25) = time_tot
    par_455(26) = time_tot
    par_455(27) = 0
    par_455(28) = time_tot
    
    par_dem(1) = -1
    par_dem(2) = 0
    par_dem(3) = time_step
    par_dem(4) = time_step
    par_dem(5) = time_step
    par_dem(6) = 0
    par_dem(7) = time_tot
    
    par_dem(8) = -2
    par_dem(9) = 0
    par_dem(10) = time_tot
    par_dem(11) = time_tot
    par_dem(12) = time_tot
    par_dem(13) = 0
    par_dem(14) = time_tot
       
    par_lcc(1) = 2
    par_lcc(2) = 1
    par_lcc(3) = YEARS_ANALYSIS
    par_lcc(4) = YEARS_LOAN
    par_lcc(5) = COMM
    
    par_lcc(6) = FUEL_INFL
    par_lcc(7) = DISC_RATE
    par_lcc(8) = TAX_RATE
    par_lcc(9) = INFLATION
    par_lcc(10) = PROP_TAX
    
    par_lcc(11) = LOAN_INT_RATE
    par_lcc(12) = 1    
    par_lcc(14) = DEPR_LIFE
    par_lcc(15) = DOWN_PAYMT
    
    par_lcc(16) = 0
    par_lcc(17) = REBATE
    par_lcc(18) = SALVAGE_FRAC
    par_lcc(19) = 0  
    
    PumpRating = max(PumpKWPer,0.)*TONS_REJECTION*0.01*3600
    
    time0 = 0.
    ni_econ = 19
    reset_econ = time_tot
    base_econ = 1
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
    !!!first cost
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
    Sizea = (0.000001+Size_TowerkW * 0.028436) !(/(10*12*1055/3600))

    COST_INSIDE = 0 !FIRST_INCR_HP_PERTON*TONS_REJECTION
    TANK_RATIO_COST_CT = .6909*sizea+.3162

    call spline(xct,yct,nct,yp1ct,ypnct,y2ct)
    call splint(xct,yct,y2ct,nct,sizea,CCCT_RATIO_COST)
    
    COST_PG_BLDG = abs(FLUID_TYPE == 2)*TONS_REJECTION*1.2*abs(Iso_Loop == 0)
    COST_PG_GHX = abs(FLUID_TYPE == 2)*TotalLength*COST_PGperM
    CCCT_RATIO_COSTMNT = -0.00001*sizea*sizea + 0.055*sizea + 1.059
    MAINT_CLOSED10 = MAINT_CLOSED10I*MULT_MNTCLOSED*abs(size_towerkW > 0.)
    FIRST_CLOSED10 = FIRST_CLOSED10I*MULT_FSTCLOSED*abs(size_towerkW > 0.)
    
    FIRST_GCHX = (TotalLength*COST_GHXperM+COST_PG_GHX)
    FIRST_SUMP = TANK_RATIO_COST_CT*960*abs(size_towerkW > 0)
    FIRST_CLOSED = ((CCCT_RATIO_COST*FIRST_CLOSED10+FIRST_SUMP)+COST_PG_BLDG)           
    FIRST_COST = COST_INSIDE+FIRST_GCHX+FIRST_CLOSED
    
    par_lcc(13) = FIRST_COST
    
    ContinFan= continuous
    CCCT_RATIO_SP = 0.9884*sizea

    xin_cont(14) = Tground + 1     !TH1
    k_g_kJhr(1) = k_g_W*3600/1000
    diff_g_m2day = k_g_kJhr(1)/(c_g_kJm3K(1)/24)
    NumSeries = abs(MaxDrillDepth >= 45)+2*abs(MaxDrillDepth < 45)
    MaxNum = ceiling((MaxFlow*8/5675)*NumSeries)
    Num_interm = int(TotalLength/MaxDrillDepth)+1
    Num = min(Num_interm,MaxNum)      !Num_interm
    REflag = abs(Num > MaxNum)
    FLOW_SIZED_GHX = 1600*num !this assumes 7.0 gpm max through each borehole - seems taht most projects target 5-6gpm
    BoreDepth = TotalLength/max(Num,1.)
    StorageVolume = (SpacingBores*.525)*(SpacingBores*.525)*3.14159*BoreDepth*Num
    ThicknessLayer = 2*BoreDepth
    NumVert = min(50,int(120/NumSeries))
    NumRad = NumSeries
    RefFlow = FLOW_SIZED_GHX*0.57/max((Num/max(NumSeries,1.)),1.) 
       
    par_ghx(1) = StorageVolume
    par_ghx(2) = BoreDepth
    par_ghx(3) = num
    par_ghx(4) = RefFlow
    par_ghx(5) = cp_fluid
    par_ghx(6) = rho_fluid
    par_ghx(7) = NumSeries
    par_ghx(8) = NumRad
    par_ghx(9) = NumVert
    par_ghx(10) = thicknesslayer
    par_ghx(11) = Tground
    par_ghx(12) = k_g_kJhr(1)
    par_ghx(13) = c_g_kJm3K(1)
    par_ghx(14) = k_g_kJhr(1)
    par_ghx(15) = c_g_kJm3K(1)
    par21 = -1
    
    FLOWDES_C = MULT_CTFLUIDFL*25*60*rho_FLUID/264.172*sizea
    AIRDES_C = MULT_CTFANFL*4970*60/35.2*1.1614*1.0808*sizea
    CCCT_RATIO_HP = max(abs(sizea <= 18.64)*(0.0183*sizea*sizea+1.4886*sizea),abs(sizea > 18.64)*(-0.014*sizea*sizea+2.0907*sizea))  
    PFANDES_C = MULT_CTFANP*2*745.7*3.6*CCCT_RATIO_HP
    GPM_SIZED_CCCT = 25*(0.000001+size_towerkw/(10*12*1055/3600))
    par_ct(1) = flowdes_c
    par_ct(2) = cp_fluid
    par_ct(3) = airdes_c
    par_ct(4) = PFANDES_C

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
    !!!Initial guess temperature - ground temperature
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
    EWT_penalty_count = 0
    its = 1
    x(its) = 0   !Tground
    its2 = 2
    therm_store = 0
    Tdiff(its) = 0
    Tdiff2(its) = 0
    Tdiff3(its) = 0
    Tdiff4(its) = 0
    repeat = 1

    outc(4) = 0
    outh(4) = 0
    out_Ct = 0
          
    time2 = 0    !delt
    out_ghx(3) = 0  !Tground 
     
    out_div = 0
    out_div(1) = 0  !Tground
    out_div(3) = 0  !Tground
    out_div(5) = 0  !Tground
    out_div2 = 0
    out_div2(1) = 0  !Tground
    out_div2(3) = 0  !Tground  
    out_div3 = 0
    out_div3(1) = 0  !Tground
    out_div3(3) = 0  !Tground  
    
    out_mix(1) = 0  !Tground  
    out_mix(2) = 0   
    out_mix2(1) = 0  !Tground  
    out_mix2(2) = 0 
    out_mix3(1) = 0  !Tground  
    out_mix3(2) = 0 

    if(opt == 0)then
        write(*,*) "*****Running main calculation loop*****"  
        write(*,*)  
        call print_bar(1,time_step+1)
    endif
    do j = 1,time_step+1
        if(j == 8760)then
            temp = 1
        endif
        i = j       
        if(opt == 0)then
            call delete_bar(i,time_step+1)
            call print_bar(i,time_step+1)
        endif

        psydat(2) = Tdb(i)
        psydat(3) = Twb(i)
        psydat(1) = 1
        call MoistAirProperties(1,0,psydat,2,status)
        T_AIR_IN_ct = psydat(2)
	    W_AIR_IN_ct = psydat(6)
	    RH_AIR_OA = psydat(4)
	    
	    T_AIR_IN = Tdb(i)        
        T_air_oa = Tdb(i) 
	     
        delta = 5   !Set the difference between the guess value and the calculated value to 1
        repeat = 1
        
        Tdiff(1) = Tdiff(its)
        Tdiff2(1) = Tdiff2(its)
        Tdiff3(1) = Tdiff3(its)
        Tdiff4(1) = Tdiff4(its)
        x(1) = x(its)
        
        its = 1     !initialize the iteration    
        its2 = 2
        etol = tol  !initialize etol to tolerance
        etol2 = tol
        etol3 = tol

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !Time calculations
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                 
        par_time(1) = 0
        par_time(2) = delt
        par_time(3) = peak_last
        par_time(4) = peak_first
        par_time(5) = summer_first
        par_time(6) = summer_last
        call timing_calcs(time2,par_time,out_time)
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !Furnace load
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if(furnace)then
            Qsupp = max(qheat(i) - Qset,0.)    !If there is a furnace and the heating load
                                              !is greater than Qset, then this is the load met by the furnace
            if(Qsupp > 0)then
                temp  = 1
            endif
            Qheat_loop = min(Qset,qheat(i))   !If there is a furnace, then the loop 
                                              !load is Qset or qheat if qheat is less than Qset            
        else
            Qheat_loop = qheat(i)
            Qsupp = 0
        endif
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !Controller
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
        call system_clock(t5,clock_rate,clock_max)
        xin_cont(4) = tdb(i)
        xin_cont(5) = twb(i) 
        xin_cont(6) = qcool(i)
        xin_cont(7) = Qheat_loop
        xin_cont(8) = out_mix(2) !flow
        xin_cont(9) = out_div3(1)   !tghx_in
        xin_cont(15) = out_mix(1)   !tfl 
        xin_cont(16) = strat   
        xin_cont(17) = out_time(14)
        xin_cont(18) = qcool(i)
        xin_cont(19) = Qheat_loop
        xin_cont(25) = out_ghx(1)
        xin_cont(29) = size_towerkw*3600/qcool(i)*abs(qcool(i) > 0)
        xin_cont(32) = out_ct(1)
        call control_main(time2,xin_cont,out_cont) 
        
        call system_clock(t6,clock_rate,clock_max)
        ticont = ticont + real(t6 - t5)/real(clock_rate)
        
        RetAirTemp = abs(t_air_ret_h*(out_time(2)<=4)+ &
                     t_air_ret_c*(out_time(2)<=9)*(out_time(2)>4)+ &
                     t_air_ret_h*(out_time(2)>9))     
                       
    do while(repeat == 1)
101           T_FLUID_IN = x(its)
              out_div(1) = T_fluid_in
              out_div(3) = T_fluid_in
              out_div(5) = T_fluid_in
              
              call system_clock(t3,clock_rate,clock_max)
                                                  
                XIN_hpc(1) = out_div(1)           !C
                XIN_hpc(2) = 1000.  
                XIN_hpc(3) = T_AIR_OA            !C
                XIN_hpc(4) = 0.008       !W_AIR_OA            !%
                XIN_hpc(5) = RH_AIR_OA     !FRACTION
                XIN_hpc(6) = 1.   !P_AIR_OA
                XIN_hpc(7) = Qcool(i) 
                XIN_hpc(8) = continuous
                XIN_hpc(9) = RetAirTemp
                !!!Cooling hp!!!
                call heat_pump_c(par_hp,xin_hpc,outc)    
                
                xin_hph = Xin_hpc
                xin_hph(7) = Qheat_loop
                !!!Heating hp!!!      
                call heat_pump_h(par_hp,xin_hph,outh)                 
                
                call system_clock(t4,clock_rate,clock_max)
                tihp = tihp + real(t4 - t3)/real(clock_rate)                              
                              
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !Mixing valve, downstream of HPs
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                bypassflow = max(0.,0.3*maxflow - outc(2) - outh(2))* &
                             abs(outc(2) + outh(2) + out_ghx(2) + out_ct(2) > 0.)
                out_mix(2) = outc(2) + outh(2) + bypassflow
                if(out_mix(2) > 0.) then
                    out_mix(1) = (outc(1)*outc(2) + outh(1)*outh(2) + bypassflow*T_FLUID_IN)/out_mix(2)
                else
                    out_mix(1) = T_FLUID_IN
                endif
                
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !Diverting valve at inlet to cooling tower
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
                out_div2(1) = out_mix(1)
                out_div2(3) = out_mix(1)
                out_div2(2) = min(GPM_SIZED_CCCT*60*rho_fluid/264.172,out_mix(2))*out_cont(1)              
                out_div2(4) = max(0.,out_mix(2) - out_div2(2))
                
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !Cooling Tower
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                      
                GAMMA_FAN = out_cont(2) 
        	    
        	    call system_clock(t7,clock_rate,clock_max)                
                
                Xin_ct(1) = out_div2(1)
                Xin_ct(2) = out_div2(2)
                Xin_ct(3) = T_AIR_IN_ct
                Xin_ct(4) = W_AIR_IN_ct
                Xin_ct(5) = RH_AIR_oa
                xin_ct(6) = 1.
                Xin_ct(7) = GAMMA_FAN
 
                call cooling_tower(time2,delt,Xin_ct,par_ct,Out_ct)
                
                call system_clock(t8,clock_rate,clock_max)
                tict = tict + real(t8 - t7)/real(clock_rate)
                
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !Mixing valve, downstream of CT
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
                out_mix2(2) = out_mix(2)
                if(out_mix2(2) > 0.) then
                    out_mix2(1) = (out_div2(3)*out_div2(4) + out_ct(1)*out_ct(2))/out_mix2(2)
                else
                    out_mix2(1) = out_mix(1)
                endif
                
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !Diverting valve, inlet to GHX, determine how much to bypass GHX
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
                out_div3(1) = out_mix2(1)
                out_div3(3) = out_mix2(1)
                out_div3(2) = min(flow_sized_ghx,out_mix2(2))*out_cont(3)
                out_div3(4) = out_mix2(2) - out_div3(2)  
                
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !ghx
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                call system_clock(t9,clock_rate,clock_max)
                                
                TFCHS = out_div3(1)
                Tgk = TFCHS
                QWCHS = out_div3(2)
                CIRCUL = 0. 
                TACHSH = Tdb(i)
                TACHS = Tdb(i)            

                call GHX( Par21,par_ghx, TFCHS, QWCHS, TACHSH, TACHS, CIRCUL,time2,out_ghx)
                !,&
                          !Tij,zdist,rdist,FRRij  )   
                !timeprint(k,1) = time2    
                call system_clock(t10,clock_rate,clock_max)
                tighx = tighx + real(t10 - t9)/real(clock_rate)         
                                
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !Mixing valve downstream of GHX
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
                out_mix3(2) = out_mix2(2)  !out_ghx(2) + out_div3(4)
                if (out_mix3(2) > 0.) then
                    out_mix3(1) = (out_ghx(1)*out_ghx(2) + out_div3(3)*out_div3(4))/out_mix3(2)
                else
                    out_mix3(1) = out_mix2(1)
                endif              
                                                               
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            !Pump
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
            GHX_Off_Fact = (PdropPctGHX)*abs(out_ghx(2) <= 0)
            PumpFlow = max(out_mix3(2),0.3*MaxFlow)*abs(out_mix3(2) > 0.)
            PD_TOT = PD_TOT_MAX*PumpFlow/MaxFlow*(1-GHX_Off_Fact) 
            
            if (pumpflow > 0) then
                temp = 1
            endif
            
            par_p1(1) = cp_fluid
            par_p1(2) = rho_fluid
            par_p1(3) = 0.
            xin_p1(1) = out_mix3(1)
            xin_p1(2) = pumpflow
            xin_p1(3) = eff_pump
            xin_p1(4) = eff_motor
            xin_p1(5) = pd_tot
            call vspump742(par_p1,xin_p1,out_p1) 
            
            par_p2(1) = maxflow
            par_p2(2) = cp_fluid
            par_p2(3) = pumprating
            par_p2(4) = 0.
            par_p2(5) = 3
            par_p2(6) = 0.1
            par_p2(7) = 0.05
            par_p2(8) = 0.88
            xin_p2(1) = out_mix3(1)
            xin_p2(2) = pumpflow
            xin_p2(3) = eff_pump
            xin_p2(4) = eff_motor2
            call vspump743(par_p2,xin_p2,out_p2)
            
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            !calculations
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                x(its+1) = out_mix3(1)
    
                Tdiff(its+1) = out_mix(1) - out_mix3(1)
                etol = max(tol,abs(tol*(out_mix(1) - out_mix3(1))))
                delta = abs(Tdiff(its+1) - Tdiff(its))

                Tdiff2(its+1) = out_ghx(1) - Tgk    
                etol2 = max(tol,abs(tol*(out_ghx(1) - Tgk )))   !calculate relative tolerance
                delta2 = abs(Tdiff2(its+1) - Tdiff2(its))      !calculate delta   
                                                   
                Tdiff3(its+1) = out_ct(1) - out_div2(1)   
                etol3 = max(tol,abs(tol*(out_ct(1) - out_div2(1))))   !calculate relative tolerance
                delta3 = abs(Tdiff3(its+1) - Tdiff3(its))      !calculate delta   
                           
            if((delta > etol).or.(delta2 > etol2).or.(delta3 > etol3))then
                repeat = 1
            else
                repeat = 0
            endif
                
                 
500         if(its > max_its)then
                repeat = 0
            endif 
            
            if(its > 50) then
                temp = 1
            endif
            
            its = its + 1   !increment the iteration   
            its2 = its2 + 1
    enddo

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    xin_cont(20) = out_cont(1)  !set the value for whether or not the CT was on 
    xin_cont(21) = out_cont(4)  !set the value for whether or not the CT was on hispeed
    xin_cont(22) = out_cont(3)  !set the value for whther or not the GHX was on 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Solution is converged at the current time step; print out whatever you want
    !and go to the next time step
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    POW_SPRAY_KW = 0.333*CCCT_RATIO_SP*out_cont(1)*.7457*abs(size_towerkw > 0.)
    plr = max((qcool(i) + Qheat_loop)/abs(MaxCoolkJhr),(qcool(i) + Qheat_loop)/MaxHeatkJhr)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Temperature violation
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    Supplemental_Device = (out_ct(3))/1055
    QHCorPLR=abs(Qheat_loop/MaxHeatkJhr)*abs(abs(Qheat_loop) >= abs(qcool(i)))
    CNSTR40_1min=max(abs(Time2 < 96),abs(QHCorPLR < 0.1)*abs(out_mix3(1) < Tground))*Tground &
                     +((1-abs(QHCorPLR < 0.1)*abs(out_mix3(1) < Tground))*abs(Time2 >= 96))*out_mix3(1)
    
    EWT_PENALTY = (9000*TONS_REJECTION)*MAX(abs(CNSTR40_1min < TMIN_GHP),abs(CNSTR40_1min > TMAX_GHP))

    if(.not.first_thru)then
        if((EWT_PENALTY > 1).and.(opt > 0))then
            LCC2 = 2*exit_early_factor2*time_tot/time2  
            Qf_cap = maxval(Qsupp_print)/eff_furnace*0.0009478
            exit_early = 1
            T_in(j:(int(ilen*yr/delt)+1)) = 0
            goto 501
        else
            exit_early = 0
        endif
    endif
    
    TMIN_FLUID = 40*abs(FLUID_TYPE == 1)+20*abs(FLUID_TYPE == 2)   
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !water cost
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    EVAP_GPH_C = 0.87*(out_ct(5)*0.26516)       !2.20462/8.3144)  
    CYCLE_CONC = 2
    BLOW_GPH_C = EVAP_GPH_C/(CYCLE_CONC-1)  
    SPRAY_GPM = 35*SIZEa
    DRIFT_GPH_C = .0001*SPRAY_GPM*60*abs(out_ct(2)>0.) 
    GALLONS_WATER = DRIFT_GPH_C+EVAP_GPH_C+BLOW_GPH_C
    FT3_WATER = GALLONS_WATER*0.13368       !/7.48052
    COST_WATER = FT3_WATER*0.02831685*COST_WATperM3
    COST_TRTMNT_C = out_time(13)*COST_WATER*0.556*abs(size_towerkw > 0.)     !(10/18)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !maintenance cost
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    MAINT_GCHX = abs(FLUID_TYPE == 2)*((COST_PG_GHX+COST_PG_BLDG)/10+MAINT_PG)
    MAINT_CLOSED = (COST_TRTMNT_C+MAINT_CLOSED10*CCCT_RATIO_COSTMNT)*abs(size_towerkw > 0.)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !energy cost
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
    if(qcool(i) < Qheat_loop)then
        FanOnlyEnergy = (1-PLR)*0.15*ContinFan*out_time(18)*abs(MaxHeatkJhr)/COP_h
    else
        FanOnlyEnergy = (1-PLR)*0.15*ContinFan*out_time(18)*abs(MaxCoolkJhr)/COP_c
    endif
    temp = abs(MaxCoolkJhr)/COP_c
    temp = abs(MaxHeatkJhr)/COP_h
    P_HP_int = (outc(3)+outh(3))+FanOnlyEnergy
    PUMP_POWER=abs(PumpKWPer < 0)*out_p1(3)+abs(PumpKWPer >= 0)*out_p2(3)*(1-GHX_Off_Fact)
    COST_KWH = COST_PKkwh*out_time(17)*out_time(16)+COST_OFFPKkwh*(1-out_time(17))*out_time(16) &
               + Cost_kWh_OnW*(out_time(17))*(1-out_time(16)) &
               + Cost_kWh_OffW*(1-out_time(17))*(1-out_time(16))  
    COST_KW_ = COST_kW*(out_time(16))+Cost_kW_W*(1-out_time(16))
    ELEC_KW = POW_SPRAY_KW+(P_HP_int+out_ct(4)+PUMP_POWER)/3600
    COST_ELEC = ELEC_KW*COST_KWH
    POW_SPRAY_KJHR = POW_SPRAY_KW*3600
        
    xin_455(1) = elec_kw
    xin_455(2) = elec_kw
    xin_455(3) = cnstr40_1min
    xin_455(4) = out_div3(2)    !flow rate into the GHX  
    call integrator455(time2,delt,time0,time_tot,nin,par_455,xin_455,out_455,store_455,&
                       first_run_455)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !heat pump outputs
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    Q_ABS = -MIN(0.,outc(4))-MIN(0.,outh(4))
    Q_REJ = MAX(0.,outc(4))+MAX(0.,outh(4))
    Q_CL_HP = -MIN(0.,outc(6))-MIN(0.,outh(6))
    Q_HT_HP = MAX(0.,outc(6))+MAX(0.,outh(6))
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !econ integ inputs
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    COST_UTILITY = COST_ELEC+COST_WATER+EWT_PENALTY+cost_gas*qsupp/eff_furnace
    
    COST_RATCHET = out_455(13)*COST_KWR*out_time(13)*12
    COST_DEMAND = out_455(3)*COST_KW_*out_time(12)
    COST_MAINT = (MAINT_GCHX+MAINT_CLOSED)*out_time(13)
    COST_IN_ANNTOT = (COST_IN_ANN)*out_time(13)
    COST_YEARLY = COST_DEMAND+COST_MAINT+COST_IN_ANNTOT+COST_RATCHET
    
    XIN_econ(1) = PUMP_POWER
    XIN_econ(2) = out_ct(4)
    XIN_econ(3) = out_ct(3)
    XIN_econ(4) = POW_SPRAY_KJHR
    XIN_econ(5) = out_ghx(4)
    
    XIN_econ(6) = p_hp_int
    XIN_econ(7) = Q_REJ
    XIN_econ(8) = Q_ABS
    XIN_econ(9) = Q_HT_HP
    XIN_econ(10) = Q_CL_HP
    
    XIN_econ(11) = GALLONS_WATER
    XIN_econ(12) = DRIFT_GPH_C
    XIN_econ(13) = EVAP_GPH_C
    XIN_econ(14) = BLOW_GPH_C
    XIN_econ(15) = ELEC_KW
    
    XIN_econ(16) = COST_ELEC
    XIN_econ(17) = EWT_PENALTY
    XIN_econ(18) = COST_WATER
    
    Xin_econ(19) = cost_gas*qsupp/eff_furnace
    
    call integrator24(time2,delt,time_tot,ni_econ,XIN_econ,reset_econ,base_econ,out_econ,store_econ,&
                      first_run_econ)
    
    xin_dem(1) = cost_demand
    xin_dem(2) = cost_maint
    call integrator455_2(time2,delt,time0,time_tot,nin_2,par_dem,xin_dem,out_dem,store_dem,first_run_dem)
    
    MAX_FLOW_LPS = out_455(33)*0.0002782278
    POW_TOTAL = out_econ(1)+out_econ(2)+out_econ(4)+out_econ(6)
    NUM_PENALTIES = out_econ(17)/max(0.0001,(200*TONS_REJECTION))
    TEMP_ANNCHANGE = (out_ghx(3)-TGROUND)*8760/(time_tot-0.)
    OptLimits = abs(size_towerkw > 0.)*(abs(Sizea > 67)+abs(Sizea < 0.4))
    DataLimits = abs(size_towerkw > 0.)*abs(Sizea > 52)
    EWT_HP=out_mix3(1)*9/5+32
    LWT_HP=out_mix(1)*9/5+32
    LWT_GHX=out_ghx(1)*9/5+32
    Flow_HP=out_mix(2)/60/rho_FLUID*264.172 
    Power_Tower=POW_SPRAY_KW+out_ct(4)/3600
    Power_Pumping=PUMP_POWER/3600
    Power_HPs=P_HP_int/3600
    
    XIN_lcc(1) = COST_UTILITY 
    XIN_lcc(2) = COST_YEARLY
    
    if(time2 >= time_tot - 2*delt)then
        temp = 1
    endif
    call econ_analysis(TIME2,delt,time_tot,time0,Par_lcc,XIN_lcc,OUT_lcc)

        
    time_print(j) = time2
    time_day(j) = out_time(14)
    T_hp_h(j) = outh(1)
    F_hp_h(j) = outh(2)
    T_hp_C(j) = outc(1)
    F_hp_C(j) = outc(2)
    T_ghx(j) = out_ghx(1)
    F_ghx(j) = out_ghx(2)
    T_ct(j) = out_ct(1)
    F_ct(j) = out_ct(2)
    T_in(j) = out_mix3(1)
    Tamb(j) = Tdb(i)
    T_ct_in(j) = out_div2(1)
    T_ghx_in(j) = out_div3(1)
    Tg_rad1(j) = out_ghx(10)
    Twet(j) = Twb(i)

    Q_hp_h(j) = outh(4)
    Q_hp_c(j) = outc(4)
    Q_ghx(j) = out_ghx(4)
    if(Q_ghx(j) < 0) qgh = qgh + Q_ghx(j)    
    if(Q_ghx(j) > 0) qgc = qgc + Q_ghx(j)
    
    if(Q_ghx(j) > 0)then 
        qgc_ts = abs(Q_ghx(j))
        qgh_ts = 0
    elseif(Q_ghx(j) < 0)then
        qgh_ts = abs(Q_ghx(j))
        qgc_ts = 0
    else
        qgc_ts = 0
        qgh_ts = 0
    endif
    
    Q_ct(j) = out_ct(3)
    Q_cool(j) = qcool(i)
    Q_heat(j) = qheat(i)
    T_store(j) = out_ghx(3)
    Y_closed(j) = out_cont(1)
    Y_fan_closed(j) = out_cont(2)
    Y_GHX(j) = out_cont(3)
    PCon(j) = out_cont(8)
    Qratio(j) = size_towerkw*3600/qcool(i)*abs(qcool(i) > 0)
    
    Ptower(j) = Power_Tower
    Ppumping(j) = Power_Pumping
    Phps(j) = Power_HPs   
    debal(j) = outc(4) - out_ghx(4) - out_ct(3)

    Celec(j) = cost_elec
    Cwater(j) = cost_water
    Cmaint(j) = cost_maint
    CGHX(j) = first_gchx
    Csump(j) = first_sump
    CCT(j) = first_closed
    Cfirst(j) = first_cost
    intfuelcost(j) = OUT_lcc(1)
    p1(j) = OUT_lcc(7)
    p2(j) = OUT_lcc(8)

    ITS_cum(j) = its
    delta_end(j) = delta
    etol_end(j) = etol
    delta2_end(j) = delta2
    etol2_end(j) = etol2
    delta3_end(j) = delta3
    etol3_end(j) = etol3
    
    demandp(j) = cost_demand   !out_455(3)
    ratchetp(j) = cost_ratchet  !out_455(13)
    cost_utilityp(j) = cost_utility
    cost_yearlyp(j) = cost_yearly
    cost_kwhp(j) = cost_kwh
    cost_elecp(j) = cost_elec
    elec_kwp(j) = elec_kw
    
    Qsupp_print(j) = Qsupp
    Qheat_loop_print(j) = Qheat_loop
       
    zprint(k:k+14) = zdist
    !tprint(k:k+14,:) = Tij
   ! frrprint(k:k+14,1:15)  = FRRij
    k = k+15   
    
    !Energy in kJ/hr
    Egas(j) = qsupp/eff_furnace
    Eelect(j) = elec_kw*3600
    
    Temp_Print(j,:) = (/T_ghx_in(j),T_ghx(j),T_store(j),T_ct_in(j),T_ct(j),T_in(j),T_hp_C(j),T_hp_h(j), &
                       T_ct_in(j),Tdb(j),Twb(j)/)
    Flow_Print(j,:) = (/out_mix3(2),F_hp_C(j),F_hp_h(j),F_ct(j),F_ghx(j)/)        
    Power_Print(j,:) = (/Power_Tower,Power_HPs,Power_Pumping/)   
    Heat_print(j,:) = (/qgh_ts,qgc_ts,Q_ct(j),abs(Q_hp_h(j)),Q_hp_c(j),Qsupp,Q_cool(j),Q_heat(j)/)
    Misc_print(j,:) = (/Y_fan_closed(j)/)   
    !Misc_print(j,:) = (/Y_fan_closed(j),PLR/)          
      
    time2 = time2 + delt
    end do
    
    !Total source energy
    e_source = 3.34
    g_source = 1.047
    elec_total = sum(Eelect)*e_source
    gas_total = sum(Egas)*g_source  
    
    furnace_tot = sum(Qsupp_print)
    ghx_tot = sum(Qheat_loop_print)
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Done at all time steps. Print out final data, such as LCC, close any open files
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
    
    if(opt == 0)then
        write(*,*)
        write(*,*)
        write(*,'(a)') '*****Main calculation loop complete*****'
        write(*,*)
!        IF ((time2-0.)/3600. < 8760.d0) THEN
!          write(*,'(a)') "Simulation < 1 year; annual fuel cost computed based on fuel cost during the simulation"
!        ELSEIF ((time2-0.)/3600. > 8760.d0) THEN
!          write(*,'(a)') "Simulation > 1 year; annual fuel cost computed based on average fuel cost during the simulation"
!        ENDIF
        write(*,*)
    endif
    
    Qf_cap = maxval(Qsupp_print)/eff_furnace*0.0009478  !kBTU

    IINCR_COOL = 1/(1-0.00551*(out_455(23)*9/5+32-95))-1
    IINCR_HEAT = 1/(1+0.01390*(out_455(22)*9/5+32-35))-1
    IINCR_TONS = ((abs(MaxCoolkJhr) > abs(MaxHeatkJhr))*(abs(MaxCoolkJhr/1055.1)*IINCR_COOL) &
                + (abs(MaxHeatkJhr) > abs(MaxCoolkJhr))*(abs(MaxHeatkJhr/1055.1)*IINCR_HEAT))/12

    GHX_load_total = T_store(time_step+1) - T_store(1)
    !(qgc + qgh)*delt/(c_g_kjm3k(1)*StorageVolume)  
    ghx_Ratio = abs(qgc/qgh)
    if(abs(GHX_load_total) > load_pen_fact)then
        Imbalance_pen = (abs(GHX_load_total) - load_pen_fact)*out_lcc(3)*10        !(GHX_load_total-(1+load_pen_fact/100))*100*Out_lcc(3) 
    else
        Imbalance_pen = 0
    endif
    
    if((energy_reduction > 0).and.(energy_reduction < 1000))then
        e1 = elec_total + gas_total
        e2 = comparison_results(1,2)/e1
        if(e2 < energy_reduction)then
            energy_pen = abs(energy_reduction - e2)*out_lcc(3)*10
        else
            energy_pen = 0
        endif
    else
        energy_pen = 0
    endif

    LCC2 = OUT_lcc(3)+OUT_lcc(8)*FIRST_INCR_HP_PERTON*IINCR_TONS+Imbalance_pen+energy_pen
    total_fuel_cost = out_lcc(1)
    
    if(first_thru) exit_early_factor2 = min(LCC2,1e8)
    first_thru = .false.
    
! deallocate(qcool)
!    deallocate(qheat)
!    deallocate(tdb)
!    deallocate(twb) 
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Check for exit code
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
501    open(unit = 40, file = 'stop_code.txt')
    read(40,*)
    read(40,*)
    read(40,*) stop_run
    if(stop_run == 'stop')  stop 5555
    close(unit = 40)

    if(PC_index == 0) PC_index = 1
    PCon_sum(PC_index,1) = sum(PCon)
    PC_index = PC_index + 1
    
    write(21,31) totallength,size_towerkw,Qset_ratio,ct1,ct2, &
                 ghx1,TPC,DTPC,LCC2,Imbalance_pen, &
                 energy_pen
    31 format(f14.5,',',f14.5,',',f4.3,',',f14.5,',',f14.5,',', &
              f14.5,',',f14.5,',',f14.5,',',f20.3,',',f20.3,',', &
              f20.3)
              
!    write(21,31) totallength,size_towerkw,ct1,ct2,ghx1, &
!                 TPC,DTPC,LCC2,exit_early,PCon_sum(PC_index-1,1), &
!                 Imbalance_pen,Qset_ratio,Qset,GHX_load_total,ghx_ratio, &
!                 energy_pen,m,b
!    31 format(f14.5,',',f14.5,',',f14.5,',',f14.5,',',f14.5,',', &
!              f14.5,',',f14.5,',',f20.3,',',i2,',',f20.3,',', &
!              f20.3,',',f4.3,',',f20.3,',',f20.3,',',f20.3,',', &
!              f20.3,',',f8.4,',',f8.4)               
     
    !write(*,*) "totallength   ","size_towerkw   ","tc1   ","dt1   ","dt2   ","dt3   ","tc2   ","TPC   ","DTPC   ","LCC2   ","Qset_ratio   " 
!    write(*,30) totallength,size_towerkw,Qset_ratio,ct1,ct2, &
!                ghx1,m,b,TPC,DTPC, &
!                LCC2  
    write(*,*) "GHX Length (m): ", totallength
    write(*,*) "CT size (kW): ",size_towerkw
    write(*,*) "Furnace Capacity (kBtu/hr):",Qf_cap
    write(*,*) "CT1 Set Point (C):",ct1
    write(*,*) "CT2 Set Point (C):",ct2
    write(*,*) "GHX1 Set Point (C):",ghx1
    if((strat == 5).or.(strat == 7).or.(strat == 8).or.(strat == 13))then
        write(*,*) "TPC Set Point (C):",TPC
        write(*,*) "DTPC Set Point (C):",DTPC
        write(*,*) "Life Cycle Cost ($):",LCC2  
    else
        write(*,*) "Life Cycle Cost ($):",LCC2  
    endif     
    write(*,*)
    write(*,*)      
    30 format(f14.0,',',f14.1,',',f5.3,',',f4.1,',',f4.1,',', &
              f4.1,',',f8.4,',',f8.4,',',f4.1,','f4.1,',', &
              f14.1) 
    
    call system_clock(t11,clock_rate,clock_max)
    
    comparison_results(2,1) = lcc2
    comparison_results(2,2) = elec_total + gas_total
    comparison_results(2,3:10) = (/totallength,size_towerkw,ct1,ct2,ghx1,tpc,dtpc,Qf_cap/)
    
    write(16,1000) lcc2,first_closed,out_econ(16),out_econ(15),size_towerkw, &
                Qf_cap,ct1,ct2,totallength,ghx1, &
                tpc,dtpc,temp_annchange,out_lcc(1),out_lcc(2), &   
                out_455(22),out_455(23),out_econ(1),out_econ(2),out_econ(3), &
                out_econ(4),out_econ(5),out_econ(6),out_econ(7),out_econ(8), &                
                out_econ(11),out_econ(12),out_econ(13),out_econ(14),first_gchx, &
                first_cost,maint_gchx,maint_closed,num_penalties,out_econ(18), &
                out_dem(1),out_dem(11),pd_tot,pow_total,storagevolume, &                
                num,PCon_sum(1,1)*delt,Imbalance_pen,out_econ(19),GHX_load_total, &
                ghx_ratio,elec_total, gas_total,elec_total+gas_total,energy_pen, & 
                furnace_tot,ghx_tot
                
1000    format(f15.2,',',f15.2,',',f15.2,',',f15.2,',',f10.2,',', &
               f6.2,',',f6.2,',',f6.2,',',f15.2,',',f6.2,',', &
               f6.2,',',f6.2,',',f6.3,',',f15.2,',',f15.2,',', &
               f6.2,',',f6.2,',',f20.2,',',f20.2,',',f20.2,',', &
               f20.2,',',f20.2,',',f20.2,',',f20.2,',',f20.2,',', &
               f15.2,',',f15.2,',',f20.2,',',f20.2,',',f15.2,',', &
               f15.2,',',f10.2,',',f10.2,',',f20.2,',',f15.2,',', &
               f10.2,',',f10.2,',',f10.4,',',f20.2,',',f15.2,',', &
               f6.0,',',f10.2,',',f20.2,',',f15.2,',',f6.2,',', &
               f6.3,',',f20.2,',',f20.2,',',f20.2,',',f20.2,',', &
               f20.2,',',f20.2)              
    
    if(opt == 0)then        
        makedirectory = 'copy' // ' ' // trim(file6) // " ..\..\results\results_day.txt"
        call system(trim(makedirectory))     
    endif
            
    if(opt == 0)then
        write(*,*)
        if(sumT > 0)then            
            write(*,'(a)') 'Writing temperature file'             
            call print_bar(1,time_step+1)  
            do i = 1,time_step+1
                call delete_bar(i,time_step+1)
                call print_bar(i,time_step+1)  
                write(50,*) Temp_print(i,pv_2(1:sumT))
            enddo
            write(*,*)
        endif    
        close(unit = 50)
        if(sumF > 0)then
            write(*,'(a)') 'Writing flow data file'  
            call print_bar(1,time_step+1)  
            do i = 1,time_step+1
                call delete_bar(i,time_step+1)
                call print_bar(i,time_step+1) 
                write(51,*) Flow_print(i,pv_2(sumT+1:sumT+sumF)-11)
            enddo
            write(*,*)
        endif 
        close(unit = 51)
        if(sumP > 0)then
            write(*,'(a)') 'Writing power consumption file' 
            do i = 1,time_step+1
                call delete_bar(i,time_step+1)
                call print_bar(i,time_step+1) 
                write(52,*) Power_print(i,pv_2(sumT+sumF+1:sumT+sumF+sumP)-16)
            enddo
            write(*,*)
        endif 
        close(unit = 52)
        if(sumQ > 0)then
            write(*,'(a)') 'Writing heat transfer file'
            do i = 1,time_step+1
                call delete_bar(i,time_step+1)
                call print_bar(i,time_step+1) 
                write(53,*) Heat_print(i,pv_2(sumT+sumF+sumP+1:sumT+sumF+sumP+sumQ)-19)
            enddo
            write(*,*)
        endif 
        close(unit = 53)
        if(sumM > 0)then
            write(*,'(a)') 'Writing miscellaneous data file'
            do i = 1,time_step+1
                call delete_bar(i,time_step+1)
                call print_bar(i,time_step+1) 
                write(54,*) Misc_print(i,pv_2(sumT+sumF+sumP+sumQ+1:sumT+sumF+sumP+sumQ+sumM)-27)
            enddo
            write(*,*)
        endif 
        close(unit = 54)
    endif        
            

!    if((opt == 0).and.(print_stuff == 1)) then 
!            
!        write(*,'(a)') 'Writing detailed temperature file'  
!        call print_bar(1,time_step+1)  
!        do i = 1,time_step+1         
!            call delete_bar(i,time_step+1)
!            call print_bar(i,time_step+1)  
!            write(18,*) time_print(i),time_day(i),T_hp_h(i),T_hp_C(i),T_ghx(i), &
!                        T_ct(i),T_in(i),T_ct_in(i),T_ghx_in(i),Tamb(i),F_hp_h(i), &
!                        F_hp_C(i),F_ghx(i),F_ct(i),Tg_rad1(i),Twet(i)
!        enddo
!        write(*,*)
!        close(unit = 18)
!        
!        write(*,'(a)') 'Writing detailed load and control file'  
!        call print_bar(1,time_step+1)  
!        do i = 1,time_step+1     
!            call delete_bar(i,time_step+1)
!            call print_bar(i,time_step+1)       
!            write(19,*) time_print(i),time_day(i),Q_hp_h(i),Q_hp_c(i),Q_ghx(i),Q_ct(i),Q_cool(i),Q_heat(i), &
!                        T_store(i),Y_closed(i),Y_fan_closed(i),Y_GHX(i),PCon(i),Qsupp_print(i),Qheat_loop_print(i), &
!                        Qsupp_print(i)/Q_heat(i),Qratio(i)
!        enddo
!        write(*,*)
!        close(unit = 19)
!        
!        write(*,'(a)') 'Writing detailed convergence file'  
!        call print_bar(1,time_step+1)  
!        do i = 1,time_step+1           
!            call delete_bar(i,time_step+1)
!            call print_bar(i,time_step+1) 
!            write(23,*) time_print(i),ITS_cum(i),delta_end(i),etol_end(i),delta2_end(i),etol2_end(i), &
!                        delta3_end(i),etol3_end(i)
!        enddo    
!        write(*,*)
!        close(unit = 23)
!        
!        write(*,'(a)') 'Writing detailed power consumption file'  
!        call print_bar(1,time_step+1)  
!        do i = 1,time_step+1        
!            call delete_bar(i,time_step+1)
!            call print_bar(i,time_step+1)    
!            write(26,*) time_print(i),time_day(i),Ptower(i),Ppumping(i),Phps(i),debal(i)
!        enddo
!        write(*,*)
!        close(unit = 26)
!        
!        write(*,'(a)') 'Writing detailed economics file'  
!        call print_bar(1,time_step+1)  
!        do i = 1,time_step+1       
!            call delete_bar(i,time_step+1)
!            call print_bar(i,time_step+1)     
!            write(27,*) time_print(i),time_day(i),Celec(i),Cwater(i),Cmaint(i),CGHX(i), &
!                        Csump(i),CCT(i),Cfirst(i),intfuelcost(i),p1(i),p2(i)
!        enddo
!        write(*,*)
!        close(unit = 27)        
!        
!        close(unit = 28)
!        
!        write(*,'(a)') 'Writing another detailed economics file' 
!        call print_bar(1,time_step+1)   
!        do i = 1,time_step+1        
!            call delete_bar(i,time_step+1)
!            call print_bar(i,time_step+1)    
!            write(29,*) time_print(i),time_day(i),demandp(i),ratchetp(i),cost_utilityp(i),&
!                        cost_yearlyp(i),cost_kwhp(i),cost_elecp(i),elec_kwp(i)
!        enddo
!        write(*,*)
!        close(unit = 29)    
!        close(unit = 30)
!        
!    endif    
    
    if(sumTot > 0) deallocate(pv_2)
    deallocate(temp_print)
    deallocate(flow_print)
    deallocate(power_print)
    deallocate(heat_print)
    deallocate(misc_print)
    
    call system_clock(t12,clock_rate,clock_max)
    tiprint = tiprint + real(t12 - t11)/real(clock_rate)
    
    close(unit = 25)
    close(unit = 20)
    close(unit = 22)
    close(unit = 16)
    close(unit = 15) 
    
    call system_clock(t2,clock_rate,clock_max)
    ti1 = ti1 + real(t2 - t1)/real(clock_rate)  
    
!    open(unit = 31, file = trim(file14),recl = 2140000)
!    write(31,24) "total_time","Printer","Controller","HP","CT","GHX"
!    24 format(2x,a,5x,a,5x,a,5x,a,5x,a,5x,a)
!    write(31,*) ti1,tiprint ,ticont,tihp,tict,tighx   
!    close(unit = 31)        
!    write(*,*) ti1,tiprint,ticont,tihp,tict,tighx
!    read(*,*) junk  
return    
end subroutine HyGSHP_sub
end module HyGSHP_calc
