!Subroutine for estimating the size of the GHX and CT
subroutine ASHRAE_Sizing
!use consts
use global
implicit none
real(DP) :: Par21,q_a,h_c,h_h
integer,parameter :: headerlines = 2
integer :: i
real(DP) :: QLOAD_C, QLOAD_H, onsig_fan
real(DP) :: P_AIR_RET, DP_AIR_RET, T_AIR_RET, RH_AIR_RET, W_AIR_RET
real(DP) :: P_AIR_RETURN, DP_AIR_RETurn, T_AIR_RETURN, RH_AIR_RETURN, W_AIR_RETURN
real(DP) :: P_AIR_OAi, DP_AIR_OAi, T_AIR_OAi, RH_AIR_OAi, W_AIR_OAi
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, FLOW_FLUID, T_DHW_IN, Frac_OA,  P_AIR_IN
real(DP) :: RHO_FLUID , slopec, intc, COP_COOLING, slopec_cop, intc_cop 
real(DP) :: FRAC_CAP_C, FRAC_POW_C, LOAD_MODE, MAXFLOW, CONTPOW
real(DP) :: cp_fluid, slopeh, inth, COP_HEATING, slopeh_cop, inth_cop 
real(DP) :: FRAC_CAP_H, FRAC_POW_H, Per_PG,out(6),tboil 
real(DP) :: cop_c,cop_h
real(DP) :: maxdrilldepth,spacingbores,K_g_W,c_g_kJm3K(1),k_g,c_g 
real(DP) :: alpha_g,d_b,k_grout,tau_1,tau_2,tau_f,Fo_1,Fo_2,Fo_f,G_1,G_2,G_f
real(DP) :: R_ga,R_gm,R_gd,F_sc,PLF_c,PLF_h,W_c,W_h,t_g,t_wi_c,t_wi_h,Q_w
real(DP) :: t_wo_c,t_wo_h,t_p,R_b,L_c,L_c_m,L_h,L_h_m,q_cooler,L_hyb
real(DP) :: qhmean,qset_guess,sum_qheata,sum_qcoola
character(len = 500) :: loadfilename

namelist /dst_params/ DEPTH,R0,NPIPES,Rpo,Rpi,DPIPES,Lb,Lpipe,Lgap,Gap,Tref,Par21,ISO,FRISO, &
                      THISO,RISLAM,TIMO3,TSTMAX,TGRAD,IPRE,TCMAX,TCMIN,TCPH,TAIRM,ATAIR,TAPH, &
                      ILAY,LFPRT,LPRT
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 /DST/ MaxDrillDepth,SpacingBores,K_g_W,c_g_kJm3K,tground       
                         
                      
open(unit = 3, status = "old", file = "..\..\namelist_files\op_set.txt")
read(3,nml = dst_params)
close(3)

open(unit = 3, status = "old", file = "..\..\namelist_files\op_set.txt")
read(3,nml = heat_pumps)
close(3)

open(unit = 3,file = "..\..\namelist_files\component_inputs.txt")
read(3,nml = DST) 
close(3)



h_c = 0
h_h = 0
sum_qheata = 0
sum_qcoola = 0
do i = 1,ilen
    if(qcoola(i) > 0)then   !hours of cooling
        h_c = h_c + 1
        sum_qcoola = sum_qcoola + qcoola(i)
    else
        h_c = h_c  
    endif
    if(qheata(i) > 0)then   !hours of heating
        h_h = h_h + 1
        sum_qheata = sum_qheata + qheata(i)
    else
        h_h = h_h
    endif
enddo
if(furnace == 1)then
    qhmean = sum_qheata/h_h
    Qset_guess = 0.75*((q_lh) - qhmean)
    Qset_ash = max(qhmean,Qset_guess)/(q_lh)
    q_lh = (q_lh)*(Qset_ash)    !load on the ground
endif
!q_a = -sum_qcoola + sum_qheata
q_a = q_lc + q_lh

k_g = k_g_w*0.57778924  !W/m-K converted to Btu/h-ft-F
c_g = c_g_kJm3K(1)*0.01491066   !kJ/m3-K converted to Btu/ft3-F
alpha_g = k_g/c_g*24    !convert ft2/h to ft2/day
d_b = 2*R0*3.2808399      !borehole diameter convert to ft
k_grout = Lb*0.57778924  !W/m-K converted to Btu/h-ft-F
t_g = tground*9/5+32    !C to F

!"Fourier numbers for different time periods: 1 - 10 years (3650 days),
!2 - 1 month (30 days), f - 6 hour pulse (0.25 days), "
tau_1 = 3650 ![days]	"10 years"
tau_2 = 3680 ![days]	"month"
tau_f = 3680.25 ![days]	"0.25 day"

Fo_1 = 4*alpha_g*(tau_f - tau_1)/d_b**2
Fo_2 = 4*alpha_g*(tau_f - tau_2)/d_b**2
Fo_f = 4*alpha_g*tau_f/d_b**2

!"Calculate a G factor (not to be confused with g-function from Eskilson). 
!The equations here are based on a curve fit of a chart of Fo vs. G factor, Fig. 15 in Ch. 34"
G_1 =  (log(Fo_1)+1.3165)/13.166
G_2 = (log(Fo_2)+1.3165)/13.166
G_f = (log(Fo_f)+1.3165)/13.166

!"Estimate effective thermal resistances using G factors: 
!ga - annual pulse, gm - month, gd - day"
R_ga = (G_f - G_1)/k_g
R_gm = (G_1 - G_2)/k_g
R_gd = G_2/k_g

!"Short-circuit heat loss factor accounts for the degradation in heat transfer 
!with the ground that occurs because of the heat transferred between the legs 
!of the u-tube. The value is 1.06 for 2 gpm/ton and 1.04 for 3 gpm/ton"
F_sc = 1.04

!"Part load factors are just estimates at this point"
PLF_c = 0.85
PLF_h = 0.85

!PLF_c = 0.5
!PLF_h = 0.5

!"The power input at the design load is based on the peak load and the baseline COP"
W_c = PLF_c*abs(q_lc/cop_c)
W_h = PLF_h*q_lh/cop_h

!"Based on recommendations in the text of Ch. 34"
t_wi_c = t_g + 25   ![F]
t_wi_h = t_g - 15   ![F]
!"Based on Eq. 3.9 in Kavanaugh and Rafferty"
Q_w = 3*abs(q_lh)*0.0000833334   ![gpm/ton]*convert(Btu/hr,ton)	!"flow rate assuming designing to 3 gpm/ton"
t_wo_c =t_wi_c -(q_lc-W_c)/(500*Q_w)  !units on 500 are [Btu-min/gal-hr-F]
t_wo_h = t_wi_h -(q_lh-W_h)/(500*Q_w) 

!"Based on Table 7, fourth row, the temperature penalty is estimated; this is a fuzzy number"
t_p = 3.8   ![F]

!"Borehole thermal resistance from the DST model"
R_b = 0.0939*1.731   !convert(m-K/W,h-ft-F/Btu)

!"Equation for the estimate of the length required for cooling; 
!t_p is added because it is negative for cooling"
L_c = ( q_a*R_ga+(q_lc - W_c)*(R_b + PLF_c*R_gm + R_gd*F_sc) )/ &
	     ( t_g - (t_wi_c + t_wo_c)*0.5 + t_p )
L_c_m = abs(L_c)*0.3048    !convert(ft,m)

!"Equation for the estimate of the length required for heating"
L_h = ( q_a*R_ga+(q_lh - W_h)*(R_b + PLF_h*R_gm + R_gd*F_sc) )/ &
	     ( t_g - (t_wi_h + t_wo_h)*0.5 - t_p )
L_h_m = abs(L_h)*0.3048    !convert(ft,m)

!"If it's a hybrid"
q_cooler = (L_c - L_h)/L_c*abs(q_lc)*0.000293071   !convert(Btu/hr,kW)
L_hyb = 254*q_lh*0.001/(t_g - t_wo_h)  !the units on 254 are [ft-h-F/kBtu], convert(Btu/h,kBtu/h)  
L_hyb_m = L_hyb*0.3048*abs(strat /= 0) + max(L_c_m,L_h_m)*abs(strat == 0)    !convert(ft,m)
C_CCCT = max(1.3*( abs(q_lc)*0.0000833334 - L_hyb/( 4.72*(t_g - 0) ) )*3.516852,30.)*abs(strat /= 0)
    !convert(Btu/hr,ton), units on 4.72 are [ft/ton-F], convert(ton,kW)
end subroutine ASHRAE_Sizing