      SUBROUTINE cooling_tower(time,delt,Xin,par,Out)
        
        
!-----------------------------------------------------------------------------------------------------------------------
!    DESCRIPTION: THIS COMPONENT MODELS A CLOSED CIRCUIT COOLING TOWER USING THE APPROACH OUTLINED IN THE
!    PAPER: "BUILDING ENERGY AND SYSTEM SIMULATION PROGRAMS: MODEL DEVELOPMENT, COUPLING AND INTEGRATION" 
!    BY G.ZWEIFEL, V.DORER, M.KOSCHENZ AND A.WEBER.  THE ALGORITHM ASSUMES A STEADY STATE CONDITION.
!
!    LAST MODIFIED:
!       SEPETMBER 2004 - JWT - ORIGINAL CODING
!       JULY 2005 - JWT - CHANGED THE WAY THAT THE NON-CONVERGED CASE IS HANDLED
!       OCTOBER 2005 - JWT - FIXED AN ERROR IN THE CALCULATION OF MDOT_EVAP
!       MAY 2008 - JWT - RENAMED THE DELT VARIABLE TO AVOID A PROBLEM WITH THE TIMESTEP VARIABLE DELT
!       JULY 2008 - JWT - Added the capability to to have the tower in "temperature-seeking" mode
!-----------------------------------------------------------------------------------------------------------------------
! 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
    
        integer :: MODE_AIR, N_COEFFS, status, stat
        real(DP) :: EWT_IN_DESIGN, EWT_OUT_DESIGN, FLOW_DESIGN, CP_FLUID, T_AIR_DESIGN
        real(DP) :: T_WB_DESIGN, FLOW_AIR_DESIGN, P_AIR_DESIGN, P_FAN_RATED, Y(4), psydat(9)
        real(DP) :: T_FLUID_IN, FLOW_FLUID, T_AIR_IN, W_AIR_IN, RH_AIR_IN, P_AIR_IN, GAMMA_FAN
        
        real(DP) :: T_FLUID_OUT, Q_TOWER, P_FAN, MDOT_EVAP,xin(7),par(4)
        !double precision, allocatable :: out(:)
        real(DP) :: out(5)
        real(DP) :: T_OUT_DESIRED
        real(DP) :: H_AIR_SAT_DESIGN, Q_DESIGN, ZETA, H_AIR_OUT_DESIGN
	    real(DP) :: WB_AIR_IN, H_AIR_IN, Y_NOW, Y_PREVIOUS, FY_PREVIOUS, FY_NOW
	    real(DP) :: T_FLUID_OUT_NEW, H_SAT_1, C, H_SAT_2, X, X_PREVIOUS, FX_PREVIOUS
	    real(DP) :: Y_POWER, T_AVE, HFG
	    real(DP) :: time,delt
	    integer :: psychmode, ITERS_TEMP, ITERS_AIRFLOW ,j
	    real(DP) :: converged,h_air_in_design,flow_air,x_now,fx_now
	    DATA CONVERGED/0.001/
	    logical :: ITERATE_AIRFLOW,first_run	    
	    save
	    namelist /ct/ MODE_AIR,EWT_IN_DESIGN,EWT_OUT_DESIGN,T_AIR_DESIGN,T_WB_DESIGN,P_air_design,N_COEFFS,Y  
 
	!allocate(out(5))   
    !if(.not.first_run)then
    if (time < delt) then
       ! first_run = .true.
        open(unit = 3, status = "old", file = "..\..\namelist_files\op_set.txt")
        read(3,nml = ct)
        close(3)
    endif
    
        FLOW_DESIGN = par(1)
        cp_fluid = par(2)
        flow_air_design = par(3)
        p_fan_rated = par(4)

        !       CHECK THE PARAMETERS FOR PROBLEMS AND RETURN FROM THE SUBROUTINE IF AN ERROR IS FOUND
        IF(EWT_OUT_DESIGN >= EWT_IN_DESIGN) write(*,10) "EWT_OUT_DESIGN greater than EWT_IN_DESIGN"
        10 format (a)
        IF(FLOW_DESIGN <= 0.) write(*,10) "Flow less than 0"
        IF(CP_FLUID <= 0.) write(*,10) "CP_fluid less than 0"
        IF(T_WB_DESIGN > T_AIR_DESIGN) write(*,10) "Wet bulb greater than dry bulb"
        IF(FLOW_AIR_DESIGN <= 0.) write(*,10) "Air design flow rate less than 0"
        !IF(P_AIR_DESIGN <= 0.) write(*,10) "Rated air pressure less than 0"
        !IF(P_FAN_RATED <= 0.) write(*,10) "Rated fan pressure less than 0"

        
        out = 0.
        out(1) = xin(1)
    !endif
    !else
    
    T_FLUID_IN = Xin(1)
    FLOW_FLUID = Xin(2) 
    T_AIR_IN = Xin(3) 
    W_AIR_IN = Xin(4)
    RH_AIR_IN = Xin(5)
    P_AIR_IN = xin(6) 
    GAMMA_FAN = Xin(7)   
    
    T_OUT_DESIRED = 999999.
	    
!    CHECK THE INPUTS FOR PROBLEMS
      IF(FLOW_FLUID < 0.) write(*,10) "Flow less than 0"
      IF(W_AIR_IN < 0.) write(*,10) "humidity ratio less than 0"
      IF(RH_AIR_IN < 0.) write(*,10) "relative humidity less than 0"
      IF(RH_AIR_IN > 1.) write(*,10) "relative humidity greater than 1"
      IF(P_AIR_IN <= 0.) write(*,10) "Pressure less than 0"
      IF(GAMMA_FAN < 0.) write(*,10) "Fan signal less than 0"

!    CHECK FOR FLOW
      IF(((GAMMA_FAN <= 0.).AND.(GAMMA_FAN > -0.5)) &
	  .OR.(FLOW_FLUID <= 0.)) THEN
	   T_FLUID_OUT = T_FLUID_IN
	   Q_TOWER = 0.
	   P_FAN = 0.
	   MDOT_EVAP = 0.
	   GAMMA_FAN = 0.
	    Out(1) = T_FLUID_OUT
        Out(2) = Flow_fluid
        Out(3) = Q_TOWER
        out(4) = P_FAN
        Out(5) = MDOT_EVAP
	   GOTO 100
	  ENDIF

!    GET THE INLET AIR PROPERTIES AT THE DESIGN CONDITION
      PSYDAT(1) = P_AIR_DESIGN
      PSYDAT(2) = T_AIR_DESIGN
      PSYDAT(3) = T_WB_DESIGN
      CALL MoistAirProperties(1,0,psydat,2,status)
	  H_AIR_IN_DESIGN = PSYDAT(7)
      
!    GET THE SATURATION ENTHALPY AT THE OUTLET DESIGN WATER TEMPERATURE
      PSYDAT(1) = P_AIR_DESIGN
      PSYDAT(2) = EWT_OUT_DESIGN
      PSYDAT(4) = 1.
      CALL MoistAirProperties(2,0,psydat,2,status)
  	  H_AIR_SAT_DESIGN = PSYDAT(7)
  
!   CALCULATE THE HEAT TRANSFER AT DESIGN CONDITIONS
      Q_DESIGN = FLOW_DESIGN*CP_FLUID*(EWT_IN_DESIGN-EWT_OUT_DESIGN)
		
!    CALCULATE THE OUTLET AIR PROPERTIES AT DESIGN CONDITIONS
      H_AIR_OUT_DESIGN = H_AIR_IN_DESIGN+Q_DESIGN/FLOW_AIR_DESIGN

!    MAKE SURE THE DENOMINATOR IS POSITIVE
      IF((H_AIR_OUT_DESIGN >= H_AIR_SAT_DESIGN)) THEN
         write(*,10) "h_air_out_design is greater than h_air_sat_design"
	   RETURN 
	 ENDIF
	        	
!    CALCULATE THE TOWER COEFFICIENT BASED ON THE DESIGN CONDITIONS
	 ZETA = DLOG((H_AIR_SAT_DESIGN-H_AIR_IN_DESIGN)/(H_AIR_SAT_DESIGN-H_AIR_OUT_DESIGN))

!    GET THE INLET AIR CONDITIONS
	  IF(MODE_AIR == 1) THEN
	   PSYCHMODE=4
	  ELSE
	   PSYCHMODE=2
      ENDIF

      PSYDAT(1) = P_AIR_IN
      PSYDAT(2) = T_AIR_IN
      PSYDAT(4) = RH_AIR_IN
      PSYDAT(6) = W_AIR_IN
      CALL MoistAirProperties(PSYCHMODE,1,psydat,2,status)
	  WB_AIR_IN = PSYDAT(3)
	  H_AIR_IN = PSYDAT(7)

!    GUESS AN OUTLET FLUID TEMPERATURE FOR THE FIRST ITERATION
      T_FLUID_OUT_NEW = T_FLUID_IN

!    INITIALIZE A FEW VARIABLES
      ITERS_TEMP = 1

!    SET SOME INTIAL CONDITIONS
      IF(GAMMA_FAN > 0.) THEN
	   ITERATE_AIRFLOW = .FALSE.
      ELSE IF(T_FLUID_IN <= T_OUT_DESIRED) THEN
	   T_FLUID_OUT = T_FLUID_IN
	   Q_TOWER = 0.
	   P_FAN = 0.
	   MDOT_EVAP = 0.
	   GAMMA_FAN = 0.
	   GOTO 100
      ELSE IF(DABS(T_FLUID_IN-T_OUT_DESIRED) <= CONVERGED) THEN
	   T_FLUID_OUT = T_FLUID_IN
	   Q_TOWER = 0.
	   P_FAN = 0.
	   MDOT_EVAP = 0.
	   GAMMA_FAN = 0.
	   GOTO 100
      ELSE IF(T_FLUID_IN <= WB_AIR_IN) THEN
	   T_FLUID_OUT = T_FLUID_IN
	   Q_TOWER = 0.
	   P_FAN = 0.
	   MDOT_EVAP = 0.
	   GAMMA_FAN = 0.
	   GOTO 100
	  ELSE
	   ITERATE_AIRFLOW = .TRUE.
	   Y_NOW = 1.
	   GAMMA_FAN = Y_NOW
	   Y_PREVIOUS = 0.
	   FY_PREVIOUS = T_FLUID_OUT_NEW-T_OUT_DESIRED
	  ENDIF
	  ITERS_AIRFLOW = 1

!    CALCULATE THE HEAT TRANSFER RATE
50    Q_TOWER = FLOW_FLUID*CP_FLUID*(T_FLUID_IN-T_FLUID_OUT_NEW)

!    SET THE AIR FLOW RATE
      FLOW_AIR = FLOW_AIR_DESIGN*GAMMA_FAN

!    GET THE SATURATION ENTHALPY OF THE AIR AT THE OUTLET FLUID TEMPERATURE
      PSYDAT(1) = P_AIR_IN
      PSYDAT(2) = T_FLUID_OUT_NEW
      PSYDAT(4) = 1.
      CALL MoistAirProperties(2,1,psydat,2,status)
	  H_SAT_1=PSYDAT(7)

!    GET ANOTHER GUESS OF THE SATURATION ENTHALPY OF THE AIR
      C = 0.6                                                  !RECOMMENDED FOR TUBE BUNDLES
      X = 1.-DEXP(-ZETA*(GAMMA_FAN**(C-1.)))
      H_SAT_2 = H_AIR_IN+Q_TOWER/FLOW_AIR/X

!    COMPARE THE TWO RESULTS AND INCREMENT THE GUESS
      IF(ITERS_TEMP == 1) THEN
	   X_PREVIOUS = T_FLUID_OUT_NEW
	   FX_PREVIOUS = H_SAT_1-H_SAT_2
    
       IF(DABS(FX_PREVIOUS)<= CONVERGED) THEN
            GOTO 60
	   ELSE
	        T_FLUID_OUT_NEW = WB_AIR_IN
            ITERS_TEMP = ITERS_TEMP+1
            GOTO 50
	   ENDIF
	ELSE
	   X_NOW = T_FLUID_OUT_NEW
	   FX_NOW = H_SAT_1-H_SAT_2

       IF(DABS(FX_PREVIOUS) <= CONVERGED) THEN
            GOTO 60
	   ELSE IF(ITERS_TEMP > 100) THEN
            T_FLUID_OUT_NEW = T_FLUID_OUT_NEW/2.+X_NOW/2.
            GOTO 60
	   ELSE
            T_FLUID_OUT_NEW = X_NOW-FX_NOW*(X_NOW-X_PREVIOUS)/(FX_NOW-FX_PREVIOUS)
            X_PREVIOUS = X_NOW
	        FX_PREVIOUS = FX_Now
            ITERS_TEMP = ITERS_TEMP+1
            GOTO 50
	   ENDIF
	ENDIF

60    IF(.NOT.ITERATE_AIRFLOW) GOTO 75

!    COMPARE THE TWO AIRFLOW RESULTS AND INCREMENT THE GUESS
	Y_NOW = GAMMA_FAN
	FY_NOW = T_FLUID_OUT_NEW-T_OUT_DESIRED

    IF(DABS(FY_PREVIOUS) <= CONVERGED) THEN
         GOTO 75
	ELSE IF(ITERS_AIRFLOW > 100) THEN
         GOTO 75
	ELSE IF((T_FLUID_OUT_NEW >= T_OUT_DESIRED).AND.(ITERS_AIRFLOW == 1)) THEN
         GOTO 75
	ELSE
	     IF(FY_NOW == FY_PREVIOUS) THEN
	        GAMMA_FAN = Y_NOW
	     ELSE
            GAMMA_FAN = Y_NOW-FY_NOW*(Y_NOW-Y_PREVIOUS)/(FY_NOW-FY_PREVIOUS)
         ENDIF
	     GAMMA_FAN = DMAX1(1.D-05,DMIN1(1.,GAMMA_FAN))
         Y_PREVIOUS = Y_NOW
	     FY_PREVIOUS = FY_NOW
         ITERS_AIRFLOW = ITERS_AIRFLOW+1
         ITERS_TEMP = 1
	     T_FLUID_OUT_NEW = T_FLUID_IN
         GOTO 50
      ENDIF
	   
!    CALCULATE THE HEAT TRANSFER
75    T_FLUID_OUT=T_FLUID_OUT_NEW
      Q_TOWER=FLOW_FLUID*CP_FLUID*(T_FLUID_IN-T_FLUID_OUT)

!    CALCULATE THE POWER MULTIPLIER BASED ON CONTROL SIGNAL
      Y_POWER = 0.
	  DO J = 1,N_COEFFS
	   Y_POWER = Y_POWER+Y(J)*(GAMMA_FAN**DBLE(J-1))
      ENDDO

!    SET THE POWER CONSUMPTION
	  P_FAN=DMAX1(0.,(P_FAN_RATED*Y_POWER))

!    CALCULATE THE WATER CONSUMPTION - THIS ASSUMES THAT THE AIR AND SPRAY WATER TEMPERATURES ARE RELATIVELY
!   CLOSE AND THE HEAT TRANSFER IS COMPLETELY DUE TO EVAPORATION
      T_AVE = (T_FLUID_OUT+T_FLUID_IN)/2.
      call HFGWATER(T_AVE,STAT,HFG)
      MDOT_EVAP = Q_TOWER/HFG   !estimate water mass as the heat transfer divided by the change in enthalpy of the water
                                !in other words, Q = mdot*Dh

!    WARN THE USER IF THE MODEL RAN INTO ITERATION PROBLEMS
      IF((ITERS_AIRFLOW >= 100).OR.(ITERS_TEMP >= 100)) THEN
         write(*,10) "More than 100 iterations"
      ENDIF      
              
        Out(1) = T_FLUID_OUT
        Out(2) = Flow_fluid
        Out(3) = Q_TOWER
        out(4) = P_FAN
        Out(5) = max(MDOT_EVAP,0.)
        
!    write(25,*) t_fluid_in,flow_fluid,t_air_in,w_air_in,rh_air_in,gamma_fan, &
!             flow_design,EWT_IN_DESIGN,EWT_OUT_DESIGN
    !endif
    
100   RETURN
      END subroutine cooling_tower


      subroutine HFGWATER(T1,STAT,HFG)
	  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
        
	  real(DP) :: T1, HFG  !C
	  INTEGER :: STAT

!      CHECK FOR OUT OF CORRELATION RANGE
	  IF(T1 < -164.) THEN
	     STAT=1
	     T1=-164.
	  ELSE IF(T1 > 300.) THEN
	     STAT=1
	     T1=300.
	  ELSE
	     STAT=0
	  ENDIF

!      USE A CORRELATION FROM EES TO GET THE VALUE
      IF(T1 < 0.) THEN
	     HFG = 2834.27-0.240108*T1-0.00354147*T1*T1
	  ELSE
	     HFG = 2505.03 - 2.63922*T1 + 0.00394873*T1*T1  - 0.0000242865*T1*T1*T1
	  ENDIf
	  
	  
	    
	  return
      END subroutine HFGWATER
      
      
