      SUBROUTINE heat_pump_c(par,xin,out)

!-----------------------------------------------------------------------------------------------------------------------
!    DESCRIPTION:
!    THIS COMPONENT IS INTENDED TO MODEL A WATER SOURCE (GEOTHERMAL) HEAT PUMP FOR WHICH THE HEATING AND COOLING 
!    PERFORMANCE CAN BE READ FROM A FILE OF CATALOG PERFORMANCE DATA.  IT IS ASSUMED THAT THE CATALOG DATA ACCOUNTS 
!    FOR THE ADDITIONAL HEAT INPUT BY THE FAN.
!
!    THIS SUBROUTINE WAS WRITTEN BY JEFF THORNTON OF THERMAL ENERGY SYSTEM SPECIALISTS
!
!     Modified from Type199 (which was modified from TYPE504) 
!     by Amanda Pertzborn of the Solar Energy Laboratory, UW-MAdison, 5-4-10
!     Add stuff for fan and occupancy
!-----------------------------------------------------------------------------------------------------------------------
! Copyright © 2004 Thermal Energy System Specialists, LLC. All rights reserved.
!-----------------------------------------------------------------------------------------------------------------------
    implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision    
    !inputs
    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, Power_CF_c 
    real(DP) :: FRAC_CAP_C, FRAC_POW_C, LOAD_MODE, MAXFLOW, CONTPOW
    real(DP) :: cp_fluid, slopeh, inth, COP_HEATING, slopeh_cop, inth_cop, Power_CF_h 
    real(DP) :: FRAC_CAP_H, FRAC_POW_H, Per_PG,out(6),tboil 
    real(DP) :: cop_c,cop_h,par(6),xin(9)
    !!related to interpolation
    integer,parameter ::  mc = 6
    integer,parameter ::  nc = 5

    !outputs
    real(DP) ::  T_FLUID_OUT, T_DHW_OUT, POWER, POWER_UNIT, Q_DHW, COP, EER, Q_out    

    !internal
    real(DP) :: FanOnly, H_AIR_RETURN, H_AIR_OA
    real(DP) :: psydat(9)
    logical :: ONCOOL, onheat, onfan
    integer :: onsig_c, onsig_h, psychmode
    real(DP) :: flow_air, Q_TOT_C, Q_SENS_C, Q_LAT_C, Q_REJECTED, Q_TOT_H
    real(DP) :: Q_ABSORBED, Q_AUX, POWER_BLOWER, POWER_CONT, CP_AIR
    real(DP) :: H_AIR_IN, W_AIR_IN, RHO_AIR_DRY_IN
    real(DP) :: T_K,P_KPA,AIRPROPS(5)
    real(DP) :: lps_tot, LPS_FLUID, RUNTIME_FRAC, BLOWPOW, Q_evap	
    real(DP) :: rh_air_in, q_cond

    real(DP) :: CTdp(mc), CTwb(nc), Ccap(mc*nc), Csens(mc*nc), Cpow(mc*nc)
    real(DP) :: Ccap2(mc,nc), Csens2(mc,nc), Cpow2(mc,nc)
    integer, parameter :: NN = 100
    integer :: i,j,k,l,temp
    real(DP) :: y1c,y2c,y1h,y2h,yp1,ypn,y2a(mc,nc)
    logical :: first_run
    save
    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
    if(.not.first_run)then
        open(unit = 3, status = "old", file = "..\..\namelist_files\op_set.txt")
        read(3,nml = heat_pumps)
        close(3)
        Power_CF_c = -0.0006414*Per_PG + 1 
        Power_CF_h = -0.0012963*Per_PG + 1
        RHO_FLUID = par(1)
        CP_FLUID = par(2)
        BLOWPOW = par(3)
        LPS_TOT = par(4)
        MAXFLOW = par(5)
        

        open (unit = 15, status = "old", &
            file = "..\..\interp stuff\cool.txt")
        do j = 1,mc*nc+2
            if (j == 1) then
                READ(15,*) CTdp
            elseif (j == 2) then
                READ(15,*) CTwb
            else
                READ(15,*) Ccap(j-2),Csens(j-2),Cpow(j-2)  
            endif
        enddo           
        close(unit = 15) 
        
        !          GET THE CORRECTION FACTORS FOR AIR TEMPERATURES              
        l = 0
        do k = 1,nc
            i = 1
            do j = 1,mc
                Ccap2(j,k) = Ccap(i+l)
                Csens2(j,k) = Csens(i+l)
                Cpow2(j,k) = Cpow(i+l)
                i = i + nc
            enddo
            l = l+1
        enddo
        
        out(1) = xin(1)
        out(2) = 0.
        out(3) = T_DHW_IN
        out(4) = 0.
        out(5) = 0.
        
        first_run = .true.
    !else
    endif
        T_FLUID_IN=XIN(1)            !C
        FLOW_FLUID=XIN(2)    
        T_AIR_OAi=XIN(3)              !C
        W_AIR_OAi=XIN(4)              !%
        RH_AIR_OAi=XIN(5)      !FRACTION
        P_AIR_OAi=XIN(6)
        QLOAD_C=XIN(7)
        ONSIG_FAN = xin(8)
        T_air_ret = xin(9)
        
        !    PERFORM ALL THE CALCULATION HERE FOR THIS MODEL.
          FanOnly = 0.
          
        !    DETERMINE IF THE HEAT PUMP IS OPERATING IN COOLING MODE
          IF(QLOAD_C > 0.) THEN
             ONCOOL=.TRUE.
             ONSIG_C=1
          ELSE
             ONCOOL=.FALSE.
             ONSIG_C=0
          ENDIF

        !    DETERMINE IF THE FAN SHOULD BE ON REGARDLESS OF COMPRESSOR OPERATION
          IF(ONSIG_FAN >= 0.5) THEN
             ONSIG_FAN = 1. * Continuous
             ONFAN = .TRUE.
          ELSE
             ONSIG_FAN = 0.
             ONFAN = .FALSE.
          ENDIF

        !    SET THE CORRECT INLET STATE
          IF(MODE == 1) THEN
           PSYCHMODE = 4
          ELSE
           PSYCHMODE = 2
          ENDIF
          
        !    CALL THE PSYCH ROUTINE TO DETERMINE THE PROPERTIES OF THE RETURN AIR (AFTER THE DAMPER)
          PSYDAT(1) = P_AIR_RET-DP_AIR_RET
          PSYDAT(2) = T_AIR_RET
          PSYDAT(4) = RH_AIR_RET
          PSYDAT(6) = W_AIR_RET
          
          call MoistAirProperties(psychmode,0,psydat,2,status)
          
          P_AIR_RETURN = PSYDAT(1)
          T_AIR_RETURN = PSYDAT(2)
          RH_AIR_RETURN = PSYDAT(4)
          W_AIR_RETURN = PSYDAT(6)
          H_AIR_RETURN = PSYDAT(7)

        !    CALL THE PSYCH ROUTINE TO DETERMINE THE PROPERTIES OF THE OUTSDIE AIR (AFTER THE DAMPER)
          PSYDAT(1) = P_AIR_OAi-DP_AIR_OAi
          PSYDAT(2) = T_AIR_OAi
          PSYDAT(4) = RH_AIR_OAi
          PSYDAT(6) = W_AIR_OAi
          
          call MoistAirProperties(psychmode,0,psydat,2,status)
          
          P_AIR_OA = PSYDAT(1)
          T_AIR_OA = PSYDAT(2)
          RH_AIR_OA = PSYDAT(4)
          W_AIR_OA = PSYDAT(6)
          H_AIR_OA = PSYDAT(7)
          
        !    CHECK FOR NO-FLOW CONDITIONS
          IF((ONSIG_C+ONSIG_FAN) < 0.5) THEN

             T_FLUID_OUT = T_FLUID_IN
             FLOW_FLUID = 0.
             FLOW_AIR = 0.0

             T_DHW_OUT = T_DHW_IN
             
             Q_TOT_C = 0.
             Q_SENS_C = 0.
             Q_LAT_C = 0.
             Q_REJECTED = 0.

             Q_TOT_H = 0.
             Q_ABSORBED = 0.
             Q_AUX = 0.

             POWER = 0.
             POWER_UNIT = 0.
             POWER_BLOWER = 0.
             POWER_CONT = 0.
             Q_DHW = 0.
             COP = 0.
             EER = 0.
             CP_AIR = 1.007

             Q_OUT = 0.
           
          ELSE

        !       DETERMINE THE MIXED AIR ENTHALPY AND HUMIDITY RATIO
             H_AIR_IN = ((1.-FRAC_OA)*H_AIR_RETURN+FRAC_OA*H_AIR_OA)
             W_AIR_IN = ((1.-FRAC_OA)*W_AIR_RETURN+FRAC_OA*W_AIR_OA)
             P_AIR_IN = DMIN1(P_AIR_OA,P_AIR_RETURN)

        !       CALL THE PSYCH ROUTINE TO DETERMINE THE PROPERTIES OF THE MIXED AIR
             PSYDAT(1) = P_AIR_IN
             PSYDAT(6) = W_AIR_IN
             PSYDAT(7) = H_AIR_IN
             
             call MoistAirProperties(7,1,psydat,2,status)
          
             P_AIR_IN = PSYDAT(1)
             T_AIR_IN = PSYDAT(2)
             WB_AIR_IN = PSYDAT(3)
             RH_AIR_IN = PSYDAT(4)
             W_AIR_IN = PSYDAT(6)
             H_AIR_IN = PSYDAT(7)
             RHO_AIR_DRY_IN = PSYDAT(9)                         !KG DRY AIR/M^3

        !       CALCULATE THE SPECIFIC HEAT OF THE AIR STREAM
             T_K = T_AIR_IN+273.15
             P_KPA = P_AIR_IN*101.325
             CALL AIRPROP(T_K,P_KPA,AIRPROPS)
             CP_AIR=AIRPROPS(5)
            
        !       FIND PERFORMANCE FOR COOLING AND HEATING
             IF(ONCOOL) THEN		 
        !	     CALCULATE THE FLUID FLOW RATE BASED ON MAX CAPACITY
                LPS_FLUID = 0.000012397*999999999		!L/S

        !          CALCULATE THE MASS FLOW RATE OF THE LIQUID STREAM
                FLOW_FLUID = LPS_FLUID*3.6*RHO_FLUID   !KG/H

        !	     CALCULATE THE COOLING CAPACITY BASED ON ENTERING WATER TEMP,
               Q_TOT_C = MAX(1.,(slopec*T_FLUID_IN+intc)*999999999)	 !KJ/H
        	                  
        !		 CALCULATE THE POWER AND EFFICIENCY BASED ON INPUT AND ENTERING WATER TEMP
               COP_COOLING = MAX(0.001,COP_C*(slopec_cop*T_FLUID_IN+intc_cop))
               POWER = Q_TOT_C/COP_COOLING/Power_CF_c
              	      
                call linear(CTdp,CTwb,Ccap2,mc,nc,T_AIR_IN,WB_AIR_IN,y1c)
!                call splie2(CTdp,CTwb,Ccap2,mc,nc,y2a)
!                call splin2(CTdp,CTwb,Ccap2,y2a,mc,nc,T_AIR_IN,WB_AIR_IN,y1c)
                Q_TOT_C=Q_TOT_C*y1c
                call linear(CTdp,CTwb,Cpow2,mc,nc,T_AIR_IN,WB_AIR_IN,y2c)
 !               call splie2(CTdp,CTwb,Cpow2,mc,nc,y2a)
 !               call splin2(CTdp,CTwb,Cpow2,y2a,mc,nc,T_AIR_IN,WB_AIR_IN,y2c)
                POWER=POWER*y2c
                
        !          MODIFY THE CATALOG VALUES ACCORDING TO THE USER PARAMETERS
                Q_TOT_C = Q_TOT_C*FRAC_CAP_C
                POWER = POWER*FRAC_POW_C
        	 
        !		 SET HEATING VALUES TO ZERO
                Q_AUX=0.
                Q_TOT_H=0.
                Q_ABSORBED=0.
                Q_DHW=0.

        !		 CALCULATE THE FRACTION OF TIME THE DEVICE IS OPERATING
                IF(LOAD_MODE > 0.) THEN
                  RUNTIME_FRAC = QLOAD_C/Q_TOT_C
                ELSE
                  RUNTIME_FRAC = DMIN1(1.,(QLOAD_C/Q_TOT_C))
                ENDIF

        !          RESET THE ENERGY FLOWS BASED ON THE ACTUAL RUN-TIME
                Q_TOT_C = Q_TOT_C*RUNTIME_FRAC
                POWER = POWER*RUNTIME_FRAC
        !		  ALSO RESET THE FLUID FLOW RATE, UNLESS ISO LOOP EXISTS
                IF(FLOWMODE == 1) THEN
	                FLOW_FLUID = MAXFLOW
                ELSE
	                FLOW_FLUID = FLOW_FLUID*RUNTIME_FRAC
                ENDIF

        !          CALCULATE THE COMPRESSOR (INCL. BLOWER) POWER FOR FULL-LOAD OPERATION
                POWER_UNIT = POWER-CONTPOW
                IF(POWER_UNIT < 0) POWER_UNIT=0.
        !          CALCULATE THE CONDENSER HEAT REJECTION
                Q_REJECTED = Q_TOT_C+POWER_UNIT

        !          CALCULATE THE OUTLET FLUID TEMPERATURE
                T_FLUID_OUT = DMIN1(70.,T_FLUID_IN + Q_REJECTED/FLOW_FLUID/CP_FLUID)
         
        !          CALCULATE THE COP AND EER
                COP = Q_TOT_C/MAX(0.0001,(POWER+Q_AUX))
                EER = COP*3.414

                Q_OUT = -Q_TOT_C
        !		  FanOnly = (1-PLR)*(-0.15)*continuous*onsig_fan*	
        !    1          -maxcool/COP_cool/2		  
           	
           ELSE

        !          JUST THE FAN IS RUNNING AND NOT THE REFRIGERANT CYCLE
                FLOW_FLUID=0.
                BLOWPOW= 0   
                Q_SENS_C=0.
                Q_LAT_C=0.
              
                POWER_UNIT = BLOWPOW
                POWER = CONTPOW

                Q_ABSORBED = 0.
                Q_REJECTED = 0.
                Q_AUX = 0.
                Q_COND = 0.
                Q_EVAP = 0.
                Q_DHW = 0.

                Q_OUT = POWER

        !          CALCULATE THE OUTLET FLUID TEMPERATURES
                T_FLUID_OUT = T_FLUID_IN
                T_DHW_OUT = T_DHW_IN

        !          CALCULATE THE COP AND EER
                COP = (Q_OUT)/MAX(0.0001,(POWER+Q_AUX))
                EER = COP*3.414

             ENDIF

        ENDIF
        if(flow_fluid>0)then
            temp = 1
        endif
        
        out(1) = T_FLUID_OUT
        out(2) = flow_fluid
        out(3) = POWER + Q_aux
        out(4) = FLOW_FLUID*CP_FLUID*(T_FLUID_OUT-T_FLUID_IN)
        out(5) = COP
        out(6) = Q_out
!endif   
return
end subroutine heat_pump_c


!
!       !inputs
!        double precision :: QLOAD_C, QLOAD_H, onsig_fan
!        double precision :: P_AIR_RETURN, DP_AIR_RET, T_AIR_RETURN, RH_AIR_RETURN, W_AIR_RETURN
!        double precision :: P_AIR_OA, DP_AIR_OA, T_AIR_OA, RH_AIR_OA, W_AIR_OA
!        integer :: Continuous, mode, status, FLOWMODE
!        real :: T_AIR_IN, wb_air_in
!        double precision :: T_FLUID_IN, FLOW_FLUID, T_DHW_IN, Frac_OA,  P_AIR_IN
!        double precision :: RHO_FLUID , slopec, intc, COP_COOLING, slopec_cop, intc_cop, Power_CF_c 
!        double precision :: FRAC_CAP_C, FRAC_POW_C, LOAD_MODE, MAXFLOW, CONTPOW
!        double precision :: cp_fluid, slopeh, inth, COP_HEATING, slopeh_cop, inth_cop, Power_CF_h 
!        double precision :: FRAC_CAP_H, FRAC_POW_H 
!        !!related to interpolation
!        integer :: LU_ATC, mc, nc, LU_ATH, nh, i 
!        double precision ::  T_DHW_OUT, POWER_UNIT, Q_DHW, EER
!        double precision ::  time(35), Qcool(35), Qheat(35), Tdb(35), Twb(35), Thp_in(35), FhpC(35), FhpH(35)   
!        integer :: iso_loop
!        double precision :: maxcoolkjhr, per_pg, psydat(9), dummy
!        double precision :: T_FLUID_OUT(35),POWER(35),COP(35),Q_out(35), occup(35)
!        mc = 6
!        nc = 5
!        nh = 6 
!        LU_ATC = 15
!        LU_ATH = 16
!        !outputs
!        
!        open (unit = 17, status = "old", &
!        file = "C:\Documents and Settings\apertzbo\My Documents\Research\standalone\Main_testing\Main_testing\hp_test_data.txt")
!        do i = 1,35
!            read(17,*) time(i), Qcool(i), Qheat(i), Tdb(i), Twb(i), Thp_in(i), FhpC(i), FhpH(i),occup(35)
!        enddo
!
!        Qcool = abs(Qcool) * 1055.06
!        Qheat = abs(Qheat) * 1055.06
!        Tdb = (Tdb - 32) * 5/9
!        Twb = (Twb - 32) * 5/9
!        RHO_FLUID = 1000.
!        cp_fluid = 4.19
!        
!        CONTPOW = 36.
!        LOAD_MODE = 1       
!        Iso_Loop = 0
!        MaxCoolkJhr = -5132.6 * 1055.1
!        MAXFLOW = max(-0.04463*MaxCoolkJhr,.03513*7465.1*1055.1)
!        Continuous = 0
!        slopeh = 0.022507
!        inth = 1.
!        slopec = -0.00856  
!        intc = 1.213
!        slopeh_cop = 0.014366 
!        inth_cop = 1.
!        slopec_cop = -0.02521
!        intc_cop = 1.6298       
!        
!        Frac_OA = 0.
!        P_AIR_OA = 1.
!        FRAC_POW_C = 1.
!        FRAC_CAP_C = 1.
!        FRAC_CAP_H = 1.
!        FRAC_POW_H = 1.
!        Per_PG = 0.
!        dummy = 0.
!        
!        mode = 2
!        
!        do i = 1,35
!            T_FLUID_IN = Thp_in(i)
!            T_AIR_IN = Tdb(i)
!            Flow_fluid = 1000. 
!            T_air_oa = Tdb(i)
!            
!            cop_cooling = 16.22/3.413
!            cop_heating = 4.35
!            T_AIR_RETURN = 21.
!            RH_AIR_RETURN = .42
!            P_AIR_RETURN = 1.
!            DP_AIR_RET = 0.
!            T_DHW_IN = 40
!        
!            P_AIR_IN = 1.
!            
!            psydat(2) = Tdb(i)
!            psydat(3) = Twb(i)
!            psydat(1) = 1
!            call MoistAirProperties(1,0,psydat,2,status)
!            RH_AIR_OA = psydat(4)
!            W_AIR_OA = psydat(6)
!        
!            open (unit = LU_ATC, status = "old", &
!            file = "C:\Documents and Settings\apertzbo\My Documents\Research\standalone\Main_testing\Main_testing\interp stuff\cool.txt")
!
!            open (unit = LU_ATH, status = "old", &
!            file = "C:\Documents and Settings\apertzbo\My Documents\Research\standalone\Main_testing\Main_testing\interp stuff\heat.txt")
!            
!            !Cooling hp
!            call heat_pumps(Qcool(i), dummy, occup(i), P_AIR_RETURN, DP_AIR_RET, T_AIR_RETURN, & 
!                  RH_AIR_RETURN, W_AIR_RETURN, P_AIR_OA, DP_AIR_OA, T_AIR_OA, RH_AIR_OA, W_AIR_OA, Continuous, &
!                  mode, T_FLUID_IN, FLOW_FLUID, T_DHW_IN, Frac_OA, T_AIR_IN, P_AIR_IN, RHO_FLUID, slopec, intc, &
!                  COP_COOLING, slopec_cop, intc_cop, Per_PG, FRAC_CAP_C, FRAC_POW_C, LOAD_MODE, FLOWMODE, &
!                  MAXFLOW, CONTPOW, cp_fluid, slopeh, inth, COP_HEATING, slopeh_cop, inth_cop,  &
!                  FRAC_CAP_H, FRAC_POW_H, LU_ATC, mc, nc, LU_ATH, nh,T_FLUID_OUT(i),POWER(i),&
!                  COP(i),Q_out(i))
!            
!            open(unit = 18, file = "cool.txt")
!            write(18,10) Qcool(i),T_FLUID_OUT(i),POWER(i),COP(i),Q_out(i)
!            
!            !Heating hp      
!            call heat_pumps(dummy, Qheat(i), occup(i), P_AIR_RETURN, DP_AIR_RET, T_AIR_RETURN, & 
!                  RH_AIR_RETURN, W_AIR_RETURN, P_AIR_OA, DP_AIR_OA, T_AIR_OA, RH_AIR_OA, W_AIR_OA, Continuous, &
!                  mode, T_FLUID_IN, FLOW_FLUID, T_DHW_IN, Frac_OA, T_AIR_IN, P_AIR_IN, RHO_FLUID, slopec, intc, &
!                  COP_COOLING, slopec_cop, intc_cop, Per_PG, FRAC_CAP_C, FRAC_POW_C, LOAD_MODE, FLOWMODE, &
!                  MAXFLOW, CONTPOW, cp_fluid, slopeh, inth, COP_HEATING, slopeh_cop, inth_cop,  &
!                  FRAC_CAP_H, FRAC_POW_H, LU_ATC, mc, nc, LU_ATH, nh,T_FLUID_OUT(i),POWER(i),&
!                  COP(i),Q_out(i))
!                  
!            open(unit = 19, file = "heat.txt")
!            write(19,10) Qheat(i),T_FLUID_OUT(i),POWER(i),COP(i),Q_out(i)
!            10 format (f20.4,2x,f20.4,2x,f20.4,2x,f20.4,2x,f20.4)
!
!            
!         enddo
!    close(unit = 18)
!    close(unit = 19)