[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [TRNSYS-users] Strange results from Type52



Hello, David.

Thanks for your reply. I have made a slight modification of the model code in order to trace the problem. Then I find that when the control signal drops below 0.4, the water flow inside the tube lies in the transition region. The resulting heat transfer coefficient on the tube inside surface decreases by more than 80% from the "just turbulent" situation to the "just laminar" condition. As a linear interpolation is used in the model to calculate the heat transfer coefficient within the transition region, the reduction in the coil capacity becomes much larger than the reduction in the water flow. Hence, the outlet water temperature starts to decrease when the control signal drops below 0.4. The simplest way is to rectify is to choose a smaller tube size. It helps increase the Reynolds number inside the tubes. Another problem is the adoption of a constant Nusselt number when the flow becomes laminar inside the tubes in the model code. I have made some modifications of the original model code. Here, the calculation of the Nusselt number in the laminar region is based on the formulations from "Nellis G, Klein S. Heat Transfer. New York: Cambridge University Press, 2009". For the Nusselt number in the turbulent region, the more accurate equations from "Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow. International Chemical Engineering 1976;16:359-68" is followed. These formulations are used in EES in calculating the heat transfer coefficient for internal pipe flow. Attached is my revised source code for your reference. The original equations are converted to comment statements so that you can easily revert to the original model. Try and see if it is better.

Best Regards
Chun Kwong LEE
City University of Hong Kong


David BRADLEY wrote:

Chun Kwong,
I suspect that when the water flowrate drops below a certain point, the heat transfer rate between the air and the water begins to decrease. I modified your simulation so that there is an equation that sets the pump control signal starting at a value of 1 and ending at 0 after 12 hours. That way, you can see the progression of the outlet water temperature (I also used a shorter timestep). As you say, the water outlet temperature starts out increasing but when the flow rate through your coils gets to about 4000 kg/h (1.1 L/s) the outlet temperature starts to go back down. At that point, the face velocity of water through a single tube in the coil is only about 0.2 m/s, which seems a bit slow. I don't think that you have reached a laminar/turbulent flow point because there is not a sharp break in the water outlet temperatures (which you usually see because the Reynolds number drops very suddenly). In any case, I am not sure that the model is wrong and I certainly don't see anything incorrect in your implementation.
Best,
 David
On 7/10/2013 02:30, a8304506@graduate.hku.hk wrote:

Hello.

I try to investigate the performance of a cooling coil at different chilled
water flow rates (see attached .tpf file).  When I adjust the control signal of
the water pump from 1 down to around 0.4, the chilled water leaving
temperatures increase as expected.  However, when I further reduce the control
signal, it appears that the leaving chilled water temperatures begin to
decrease. Is it the problem from the cooling coil model or any other issue? Please advise. I am using TRNSYS16.1. Thank you!

Best Regards
Chun Kwong LEE
City University of Hong Kong







_______________________________________________
TRNSYS-users mailing list
TRNSYS-users@cae.wisc.edu
https://mailman.cae.wisc.edu/listinfo/trnsys-users

--
***************************
David BRADLEY
Principal
Thermal Energy Systems Specialists, LLC
22 North Carroll Street - suite 370
Madison, WI  53703 USA

P:+1.608.274.2577
F:+1.608.278.1475
d.bradley@tess-inc.com

http://www.tess-inc.com
http://www.trnsys.com


      SUBROUTINE TYPE52(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)

C-----------------------------------------------------------------------------------------------------------------------
C  THIS ROUTINE DETERMINES THE THERMAL PERFORMANCE OF A COOLING COIL  *
C  USING AN EFFECTIVENESS MODEL AS DESCRIBED IN "METHODOLOGIES FOR    *
C  THE DESIGN AND CONTROL OF CHILLED WATER SYSTEMS", PHD THESIS,      *
C  SOLAR ENERGY LABORATORY, UNIVERSITY OF WISCONSIN, J.E. BRAUN, 1988 *
C****************************  PARAMETERS  ***************************
C
C    MODE   = SIMPLE OR DETAILED CALCULATIONS
C    IU     = UNITS: 1-SI, 2-ENGLISH
C    NROWS  = NUMBER OF ROWS OF TUBES
C    NTUBES = NUMBER TUBES PER ROW
C    LT     = LENGTH OF TUBES IN EACH ROW
C    WD     = WIDTH OF DUCT PERPENDICULAR TO THE TUBES
C    DO     = TUBE OUTER DIAMETER
C    DI     = TUBE INNER DIAMETER
C    KT     = TUBE THERMAL CONDUCTIVITY
C    FT     = FIN THICKNESS
C    FS     = FIN SPACING
C    NFINS  = NUMBER OF FINS ON ONE TUBE
C    KF     = FIN THERMAL CONDUCTIVITY
C    FMODE  = CONTINUOUS FLAT PLATE OR ANNULAR FINS
C    ADIST  = DISTANCE BETWEEN TUBE CENTERS ACROSS ROW
C    DE     = DIAMETER OF ANNULAR FIN
C    CDIST  = DISTANCE BETWEEN TUBE ROW CENTERLINES (FLOW DIRECTION)
C-----------------------------------------------------------------------------------------------------------------------
! Copyright © 2004 Solar Energy Laboratory, University of Wisconsin-Madison. All rights reserved.

C-----------------------------------------------------------------------------------------------------------------------
C    USE STATEMENTS
	USE TrnsysFunctions
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    REQUIRED BY THE MULTI-DLL VERSION OF TRNSYS
      !DEC$ATTRIBUTES DLLEXPORT :: TYPE52				!SET THE CORRECT TYPE NUMBER HERE
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    TRNSYS DECLARATIONS
      IMPLICIT NONE
	DOUBLE PRECISION XIN,OUT,TIME,PAR,T,DTDT,TIME0,TFINAL,DELT	
      INTEGER*4 INFO(15),NP,NI,NOUT,ND,NPAR,IUNIT,ITYPE,ICNTRL,
	1   NSTORED
      CHARACTER*3 YCHECK,OCHECK
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    USER DECLARATIONS - SET THE MAXIMUM NUMBER OF PARAMETERS (NP), INPUTS (NI),
C    OUTPUTS (NOUT), AND DERIVATIVES (ND) THAT MAY BE SUPPLIED FOR THIS TYPE
      PARAMETER (NP=16,NI=5,NOUT=12,ND=0,NSTORED=0)
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    REQUIRED TRNSYS DIMENSIONS
      DIMENSION XIN(NI),OUT(NOUT),PAR(NP),YCHECK(NI),OCHECK(NOUT),
	1   T(ND),DTDT(ND)
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    DECLARATIONS AND DEFINITIONS FOR THE USER-VARIABLES
      DOUBLE PRECISION MA,MW,MHUA,KW,MHUW,LT,KT,KF,JF,MCONV,NTUD,NTUW,
     & NTUO,NTUI,MSTAR,LCOIL,LIMITW,LIMITD,RNUM,TCONV,A1,A0,HCONV,H,B0,
     & B1,CTOF,C1,C0,TS,HS,TSAT,PRW,CPA,CPW,CPV,HFG,PRA,PATM,TOL,PI,HSAT
     & ,WD,DO,DI,FT,FS,ADIST,DE,CDIST,PSYDAT,RI,RO,RE,AC,FH,AX,AI,ATO,AF
     & ,AO,FRATIO,WA1,TW1,UM,TA1,TA2,WA2,TW2,QCOIL,QSENS,QLAT,RHA2,RHA1,
     & QDAIR,FDRY,GW,TPROP,REI,HI,HIT,HIL,DH,CJ1,CJ2,GA,REO,HDO,ALPHA,
     & BETA,EFF,FIN,EFFD,CPM,RA,TDP,HA1,HW1,CSTAR,TEMP,C,EPSD,TA2D,TS2D,
     & TW2W,HW2,CS,HWO,EFFW,ERAT,EPSW,G,G1,TOLD,G2,DGDT,TS1W,CK,TW2D,
     & TMAX,TMIN,GMAX,TWX,GMIN,SGNMIN,DMAX1,DMIN1,HA2,HSAVG,TSAVG,TAX,
     & HAX2,HAX,HSW,TSW,
     & GZ,NU_LAM,RE_TUR,FD,NU_TUR				!Added by Avatar Lee on 17/07/2013
	INTEGER MODE,FMODE,EMODE,ROUND,IU,IMAX,NROWS,NTUBES,NFINS,IWARN,
     & ISTAT,ITER
	LOGICAL DRY,WET,MAXDRY
      DIMENSION CPA(2),MHUA(2),CPW(2),KW(2),CPV(2),PSYDAT(9),HFG(2),
     &A0(2),A1(2),B0(2),B1(2),C0(2),C1(2),MCONV(2)
      DATA CPA/1.005, 0.240/, MHUA/0.0659, 0.0443/, CPW/4.186, 1.000/,   
     .      KW/2.043, 0.328/, CPV/1.884, 0.450/,  
     .     HFG/2452., 1054./, PRA/0.72/, PATM/1.0/
      DATA TOL/1.E-04/, IMAX/20/, PI/3.14159/, IUNIT/0/
      DATA A0/0., 32./, A1/1., 1.8/, B0/0., 7.687/, B1/1.,0.43/
     .     C0/32., 0./, C1/1.8, 1./, MCONV/1.488, 1.0/
C-----------------------------------------------------------------------------------------------------------------------
C  EMODE DETERMINES HOW ERRORS ARE HANDLED FROM CALLS TO THE PSYCH
C  SUBROUTINE. IF EMODE IS: 0 - NO ERROR MESSAGES WILL BE PRINTED,
C  1 - ERROR MESSAGES WILL BE PRINTED ONLY ONCE PER SIMULATION,
C  2 - ERROR MESSAGES WILL BE PRINTED EVERY TIMESTEMP THAT THEY OCCUR.
C
      DATA EMODE/1/
C
C*********************** STATEMENT FUNCTIONS **************************
C
C-- ROUND-OFF REAL NUMBERS TO NEAREST INTEGER
C
      ROUND(RNUM) = JFIX(RNUM + SIGN(0.5,RNUM))
C
C-- UNIT CONVERSIONS FOR TEMPERATURE AND ENTHALPY
C
C   TEMPERATURE:  SI TO ENGLISH (IU=2), SI TO SI (IU=1)
C   ENTHALPY:     ENGLISH TO SI (IU=2), SI TO SI (IU=1)
C   TEMPERATURE:  ENGLISH TO ENGLISH (IU=2), SI TO ENGLISH (IU = 1) 
C
      TCONV(T)  = T*A1(IU) + A0(IU)
      HCONV(H)  = (H - B0(IU))/B1(IU)
      CTOF(T) = T*C1(IU) + C0(IU)
C
C-- CORRELATION FOR SATURATION TEMPERATURE IN TERMS OF SATURATION
C   ENTHALPY: CORRELATION IS IN SI UNITS AND IS GOOD FROM 9.473 KJ/KG
C   TO 355.137 KJ/KG (0 C TO 55 C).
C   (CORRELATION IS FOR A TOTAL SYSTEM PRESSURE OF 1 ATMOSPHERE)
C
      TS(HS)     = -5.79013 + 6.64030E-01*HS - 5.07802E-03*HS**2 + 
     .              2.80381E-05*HS**3 - 9.47051E-08*HS**4 +
     .              1.72758E-10*HS**5 - 1.29547E-13*HS**6
C
      TSAT(HSAT) = TCONV(TS(HCONV(HSAT)))
C
C-- CORRELATION FOR KINEMATIC WATER VISCOSITY IN TERMS OF TEMPERATURE: 
C   CORRELATION IS IN ENGLISH UNITS (LBM/FT-HR AND DEG. F) AND
C   IS GOOD FROM 32 F TO 200 F.
C
      MHUW(T) = MCONV(IU)*(7.67201 - 1.39290E-01*T + 1.23322E-03*T**2
     .         - 5.32425E-06*T**3 + 8.87455E-09*T**4)
C
C-- CORRELATION FOR PRANTL NUMBER OF WATER IN TERMS OF TEMPERATURE: 
C   TEMPERATURE IS IN ENGLISH UNITS (DEG. F) AND THE CORRELATION IS
C   GOOD FROM 32 F TO 200 F.
C
      PRW(T)  = 2.73327E01 - 6.30480E-01*T + 7.56704E-03*T**2 -
     .          5.04074E-05*T**3 + 1.74546E-07*T**4 - 2.43923E-10*T**5

C-----------------------------------------------------------------------------------------------------------------------
C    TRNSYS FUNCTIONS
      TIME0=getSimulationStartTime()
      TFINAL=getSimulationStopTime()
      DELT=getSimulationTimeStep()
C    SET THE VERSION INFORMATION FOR TRNSYS
      IF(INFO(7).EQ.-2) THEN
	   INFO(12)=16
	   RETURN 1
	ENDIF
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    DO ALL THE VERY LAST CALL OF THE SIMULATION MANIPULATIONS HERE
      IF (INFO(8).EQ.-1) THEN
	   RETURN 1
	ENDIF
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    PERFORM ANY "AFTER-ITERATION" MANIPULATIONS THAT ARE REQUIRED
      IF(INFO(13).GT.0) THEN
	   RETURN 1
	ENDIF

C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    DO ALL THE VERY FIRST CALL OF THE SIMULATION MANIPULATIONS HERE
      IF (INFO(7).EQ.-1) THEN

C       RETRIEVE THE UNIT NUMBER AND TYPE NUMBER FOR THIS COMPONENT FROM THE INFO ARRAY
         IUNIT=INFO(1)
	   ITYPE=INFO(2)

C       SET SOME INFO ARRAY VARIABLES TO TELL THE TRNSYS ENGINE HOW THIS TYPE IS TO WORK
         INFO(6)=12			
         INFO(9)=1				
	   INFO(10)=0

C       SET THE REQUIRED NUMBER OF INPUTS, PARAMETERS AND DERIVATIVES THAT THE USER SHOULD SUPPLY IN THE INPUT FILE
	   NPAR=15
	       
C       CALL THE TYPE CHECK SUBROUTINE TO COMPARE WHAT THIS COMPONENT REQUIRES TO WHAT IS SUPPLIED 
	   CALL TYPECK(1,INFO,NI,NPAR,ND)

C       SET THE VARIABLE TYPES
         DATA YCHECK/'TE1','DM1','MF1','TE1','MF1'/
         DATA OCHECK/'TE1','DM1','MF1','TE1','MF1','PW1','PW1','PW1',
	1               'DM1','PC1','NAV','NAV'/

C       CALL THE INPUT-OUTPUT CHECK SUBROUTINE TO SET THE CORRECT INPUT AND OUTPUT UNITS
         CALL RCHECK(INFO,YCHECK,OCHECK)

C       RETURN TO THE CALLING PROGRAM
         RETURN 1

      ENDIF
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE - THERE ARE NO ITERATIONS AT THE INTIAL TIME
      IF (TIME.LT.(TIME0+DELT/2.D0)) THEN

C       SET THE UNIT NUMBER FOR FUTURE CALLS
         IUNIT=INFO(1)
         ITYPE=INFO(2)

C       READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
         MODE=ROUND(PAR(1))
         IF(MODE.NE.1 .AND. MODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
         IU=1
         NROWS  = ROUND(PAR(2))
         NTUBES = ROUND(PAR(3))
         LT     = PAR(4)
         WD     = PAR(5)
         DO     = PAR(6)
         DI     = PAR(7)
         KT     = PAR(8)
         FT     = PAR(9)
         FS     = PAR(10)
         NFINS  = ROUND(PAR(11))
         KF     = PAR(12)
         FMODE=ROUND(PAR(13))
         IF(FMODE.NE.1 .AND. FMODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
         ADIST  = PAR(14)
         IF (FMODE .EQ. 2) DE = ADIST
         CDIST  = PAR(15)
         PSYDAT(1) = PATM

C** COIL AREAS **
C
C   RI & RO - INSIDE AND OUTSIDE TUBE RADIUS
C        RE - ANNULAR FIN RADIUS OF EQUIVALENT RADIUS FOR 
C             CONTINUOUS FLAT FINS     
C        FH - FIN HEIGHT
C   AI & AO - INSIDE AND OUTSIDE COIL SURFACE AREAS
C    FRATIO - RATIO OF FIN AREA TO OUTSIDE COIL SURFACE AREA
C     LCOIL - LENGTH OF THE COIL SECTION IN THE AIR FLOW DIRECTION
C             (SEE KAYS & LONDON, "COMPACT HEAT EXCHANGERS")
C        AC - FLOW CROSS SECTIONAL AREA
C        AX - CROSS SECTIONAL AREA OF ONE TUBE
C
         RI  = DI/2.
         RO  = DO/2.
         IF (FMODE .EQ. 1) THEN
            RE = DSQRT(ADIST*CDIST/PI)
            AC = LT*WD-NFINS*FT*(NTUBES*ADIST+ADIST/2.)
     .           -(NTUBES*LT*DO-(NFINS*FT*DO))
         ELSE
            RE = DE/2.
            AC = LT*WD-NTUBES*LT*DO-NTUBES*NFINS*(DE-DO)*FT
         END IF

         LCOIL = NROWS * CDIST
         FH  = RE - RO
         AX  = PI*RI**2
         AI  = NTUBES*NROWS*PI*DI*LT
         ATO = NTUBES*NROWS*PI*DO*(LT-NFINS*FT)
         AF  = NFINS*NTUBES*NROWS*2.*PI*(RE**2-RO**2)
         AO  = ATO + AF
         FRATIO = AF/AO
C
C** TUBE HEAT TRANSFER COEFFICIENT ** 
C   BASED UPON INSIDE AREA
C
         UM = KT/RI/DLOG(RO/RI)

C       PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE
         OUT(1) = XIN(1)
         OUT(2) = XIN(2)
         OUT(3) = XIN(3)
         OUT(4) = XIN(4)
         OUT(5) = XIN(5)
         OUT(6) = 0.
         OUT(7) = 0.
         OUT(8) = 0.
         OUT(9) = 1.0
         OUT(10) = 0.0

C       RETURN TO THE CALLING PROGRAM
         RETURN 1

      ENDIF
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    *** ITS AN ITERATIVE CALL TO THIS COMPONENT ***
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    RE-READ THE PARAMETERS IF ANOTHER UNIT OF THIS TYPE HAS BEEN CALLED SINCE THE LAST TIME THE PARAMETERS
C    WERE READ IN
      IF(INFO(1).NE.IUNIT) THEN

C       RESET THE UNIT NUMBER
	   IUNIT=INFO(1)
	   ITYPE=INFO(2)
	    
C       READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
         MODE=ROUND(PAR(1))
         IF(MODE.NE.1 .AND. MODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
         IU=1
         NROWS  = ROUND(PAR(2))
         NTUBES = ROUND(PAR(3))
         LT     = PAR(4)
         WD     = PAR(5)
         DO     = PAR(6)
         DI     = PAR(7)
         KT     = PAR(8)
         FT     = PAR(9)
         FS     = PAR(10)
         NFINS  = ROUND(PAR(11))
         KF     = PAR(12)
         FMODE=ROUND(PAR(13))
         IF(FMODE.NE.1 .AND. FMODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
         ADIST  = PAR(14)
         IF (FMODE .EQ. 2) DE = ADIST
         CDIST  = PAR(15)
         PSYDAT(1) = PATM

C** COIL AREAS **
C
C   RI & RO - INSIDE AND OUTSIDE TUBE RADIUS
C        RE - ANNULAR FIN RADIUS OF EQUIVALENT RADIUS FOR 
C             CONTINUOUS FLAT FINS     
C        FH - FIN HEIGHT
C   AI & AO - INSIDE AND OUTSIDE COIL SURFACE AREAS
C    FRATIO - RATIO OF FIN AREA TO OUTSIDE COIL SURFACE AREA
C     LCOIL - LENGTH OF THE COIL SECTION IN THE AIR FLOW DIRECTION
C             (SEE KAYS & LONDON, "COMPACT HEAT EXCHANGERS")
C        AC - FLOW CROSS SECTIONAL AREA
C        AX - CROSS SECTIONAL AREA OF ONE TUBE
C
         RI  = DI/2.
         RO  = DO/2.
         IF (FMODE .EQ. 1) THEN
            RE = DSQRT(ADIST*CDIST/PI)
            AC = LT*WD-NFINS*FT*(NTUBES*ADIST+ADIST/2.)
     .           -(NTUBES*LT*DO-(NFINS*FT*DO))
         ELSE
            RE = DE/2.
            AC = LT*WD-NTUBES*LT*DO-NTUBES*NFINS*(DE-DO)*FT
         END IF

         LCOIL = NROWS * CDIST
         FH  = RE - RO
         AX  = PI*RI**2
         AI  = NTUBES*NROWS*PI*DI*LT
         ATO = NTUBES*NROWS*PI*DO*(LT-NFINS*FT)
         AF  = NFINS*NTUBES*NROWS*2.*PI*(RE**2-RO**2)
         AO  = ATO + AF
         FRATIO = AF/AO
C
C** TUBE HEAT TRANSFER COEFFICIENT ** 
C   BASED UPON INSIDE AREA
C
         UM = KT/RI/DLOG(RO/RI)

      ENDIF
C-----------------------------------------------------------------------------------------------------------------------


C******************************  INPUTS   ****************************
C
C  TA1     = ENTERING AIR DRY BULB TEMPERATURE (C, F)
C  WA1     = ENTERING AIR HUMIDITY RATIO
C  MA      = MASS FLOW RATE OF AIR (KG/HR, LBM/HR)
C  TW1     = ENTERING WATER TEMPERATURE (C, F)
C  MW      = MASS FLOW RATE OF WATER (KG/HR LBM/HR)
C
      TA1  = XIN(1)
      WA1  = XIN(2)
      MA   = XIN(3)
      TW1  = XIN(4)
      MW   = XIN(5)

C    CHECK THE INLET AIR STATE
      PSYDAT(1) = PATM
      PSYDAT(2) = TA1
      PSYDAT(6) = WA1
      CALL PSYCHROMETRICS(TIME,INFO,IU,4,0,PSYDAT,EMODE,ISTAT,*99)
      CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
99    TA1 = PSYDAT(2)
      WA1 = PSYDAT(6)
      RHA1 = PSYDAT(4)

C-----------------------------------------------------------------------------------------------------------------------
C************************* COIL ANALYSIS *************************
C
C  TA2     = DRY-BULB TEMPERATURE OF AIR LEAVING COIL
C  WA2     = ABSOLUTE HUMIDITY OF AIR LEAVING COIL
C  RHA2    = RELATIVE HUMIDITY OF AIR LEAVING COIL
C  TW2     = DRY-BULB TEMPERATURE OF WATER LEAVING COIL
C  QCOIL   = TOTAL HEAT TRANSFER TO COIL
C  QSENS   = SENSIBLE HEAT REMOVED FROM MOIST AIR
C  QLAT    = LATENT HEAT REMOVED FROM MOIST AIR
C  QDAIR   = SENSIBLE HEAT REMOVED FROM AIR ONLY (DRY AIR)
C  FDRY    = FRACTION OF COIL WHICH IS DRY
C  GW      = WATER MASS FLOW RATE DIVIDED BY CROSS-SECTIONAL AREA
C  REI     = INSIDE REYNOLDS NUMBER
C  HI      = INSIDE COEFFICIENT OF CONVECTION
C  GA      = AIR MASS FLOW RATE DIVIDED BY CROSS-SECTIONAL AREA 
C  REO     = OUTSIDE REYNOLDS NUMBER
C  CPM     = MOIST AIR HEAT CAPACITY
C  NTUO    = NTU'S OUTSIDE
C  NTUI    = NTU'S INSIDE
C  TDP     = AIR DEW-POINT
C  HA1     = ENTHALPY OF ENTERING AIR
C  HA2     = ENTHALPY OF EXITING AIR
    
C ***  NO FLOW  ***
      IF (MA .LE. 1.0E-06 .OR. MW .LE. 1.0E-06) THEN
         TA2   = TA1
         IF (MA .LE. 1.0E-02 .AND. MW .GT. 1.0E-02) TA2 = TW1
         WA2   = WA1
         RHA2  = RHA1
         TW2   = TW1
         QCOIL = 0.0
         QSENS = 0.0
         QLAT  = 0.0
         QDAIR = 0.0
         FDRY  = 1.0
         GOTO 500      
      END IF
C ***   FLOW    ***
C ** INSIDE FLUID COEFFICIENT ** 
C  FIND THE REYNOLDS NUMBER
      GW    = MW/NTUBES/AX
      TPROP = CTOF(TW1)
      IF ((32.0.GT.TPROP .OR. TPROP.GT.200.) .AND. 
     .    (TIME-OUT(12)) .GT. DELT/2.) THEN
	   CALL MESSAGES(-1,'THE CORRELATIONS FOR THE VISCOSITY AND PRANDT
     &L # OF WATER USED IN THE COOLING COIL MODEL WERE USED WITH A VALUE
     & OUTSIDE OF THEIR INTENDED RANGE','WARNING',INFO(1),INFO(2))
         IWARN=IWARN+1
         OUT(12) = TIME
      END IF
	REI  = GW*DI/MHUW(TPROP)
C  Added by Avatar Lee on 17/07/2013
C  Calculation of laminar Nusselt number based on Nellis G, Klein S. Heat Transfer. New York: Cambridge University Press, 2009.
	GZ = DI*MIN(REI,2300.0D0)*PRW(TPROP)/LT
	NU_LAM = 3.66+(0.049+0.02/PRW(TPROP))*GZ**1.12/(1.0+0.065*GZ**0.7)
C  Calculation of turbulent Nusselt number based on Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow.
C  International Chemical Engineering 1976;16:359-68.
	RE_TUR = MAX(REI,3000.0D0)
	FD = 4.0*(-0.0015702/LOG(RE_TUR)+0.3942031/(LOG(RE_TUR))**2.0
     &+2.5341533/(Log(RE_TUR))**3.0)
      NU_TUR = (FD/8.0)*(RE_TUR-1000.0)*PRW(TPROP)*(1.0+(DI/LT)**0.7)
     &/(1.0+12.7*SQRT(FD/8.0)*(PRW(TPROP)**(2.0/3.0)-1.0))
C  Revision ended...
C  IF THE REYNOLDS NUMBER IS GREATER THEN 3000, USE THE CORRELATION
C  FOR TURBULENT WATER FLOW IN A TUBE.  IF THE REYNOLDS NUMBER IS
C  LESS THAN 2000, USE THE CORRELATION FOR LAMINAR WATER FLOW IN A
C  TUBE.  FOR A REYNOLDS NUMBER BETWEEN 2000 AND 3000, A LINEAR
C  RELATIONSHIP BETWEEN THE CONVECTION COEFFICIENTS FOR THE TURBULENT 
C  AND LAMINAR CASES IS USED.  THE INLET TEMPERATURE IS USED FOR 
C  PROPERTY EVALUATIONS.
      IF (REI .GE. 3000.) THEN
C        HI = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4
         HI = NU_TUR*KW(IU)/DI		!Added by Avatar Lee on 17/07/2013
      ELSE IF (REI .LE. 2300.) THEN
C        HI = 4.36*KW(IU)/DI
	   HI = NU_LAM*KW(IU)/DI		!Added by Avatar Lee on 17/07/2013
      ELSE
C        HIT = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4
	   HIT = NU_TUR*KW(IU)/DI		!Added by Avatar Lee on 17/07/2013
C        HIL = 4.36*KW(IU)/DI
	   HIL = NU_LAM*KW(IU)/DI		!Added by Avatar Lee on 17/07/2013
         HI = HIL + (REI - 2000.0) * (HIT - HIL) / 1000.0
      END IF
C ** OUTSIDE COEFFICIENT FOR DRY SURFACES **
C   CORRELATION FOR OUTSIDE HEAT TRANSFER COEFFICIENT IS FROM
C   "FINNED TUBE HEAT EXCHANGER: CORRELATION OF DRY SURFACE HEAT 
C   TRANSFER DATA," A.H. ELMAHDY, ASHRAE TRANSACTIONS.  THE
C   HYDRAULIC DIAMETER IS: 2 * THE FLOW LENGTH OF THE HEAT EXCHANGER *
C   THE FLOW CROSS SECTIONAL AREA / TOTAL HEAT TRANSFER AREA.
C   THE FIN EFFICIENCY EQUATION IS FOR ANNULAR FINS.  STRAIGHT FINS
C   ARE TREATED AS ANNULAR FINS BY FINDING AN EQUIVALENT ANNULAR FIN
C   WITH THE SAME SURFACE AREA.
      DH  = 4.*LCOIL*AC/AO
      CJ1 = 0.159*(FT/FH)**0.141*(DH/FT)**0.065
      CJ2 = -0.323*(FT/FH)**0.049*(FS/FT)**0.077
      GA  = MA/AC
      REO = GA*DH/MHUA(IU)
      JF  = CJ1*REO**CJ2
      HDO = GA*CPA(IU)*JF/PRA**0.6667
C FIN EFFICIENCY CALCULATIONS
      ALPHA  = RO/RE
      BETA   = RE*DSQRT(2.*HDO/KF/FT)
      EFF = FIN(ALPHA,BETA)
      EFFD = 1. - FRATIO*(1.-EFF)
C ** DETERMINE NTU'S **
C   THE AIR SPECIFIC HEAT IS FOR MOIST AIR
      CPM  = CPA(IU) + WA1*CPV(IU)
      NTUO = EFFD*HDO*AO/(MA*CPM)
      NTUI = AI/(1./HI + 1./UM)/(MW*CPW(IU))
      RA   = MW/MA
C ** ENTHALPIES FOR ENTERING CONDITIONS **
      PSYDAT(2) = TA1
      PSYDAT(6) = WA1
      CALL PSYCHROMETRICS(TIME,INFO,IU,4,0,PSYDAT,EMODE,ISTAT,*101)
      CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
101   TDP = PSYDAT(5)
      HA1 = PSYDAT(7)
      PSYDAT(2) = TW1
      PSYDAT(3) = TW1
      CALL PSYCHROMETRICS(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*102)
      CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
102   HW1 = PSYDAT(7)
C ** DRY ANALYSIS **
C   DETERMINES THE AIR EFFECTIVENESS (EPSD) AND OUTLET WATER TEMPERATURE
C   (TW2D) ASSUMING THE COIL IS DRY
C   IF THE SURFACE TEMPERATURE AT THE OUTLET IS GREATER THAN THE AIR
C   ENTERING DEWPOINT THEN COIL IS COMPLETELY DRY, OTHERWISE DO THE WET
C   ANALYSIS
      CSTAR = CPM/RA/CPW(IU)
      NTUD   = NTUO/(1. + NTUO/NTUI*CSTAR)
      TEMP = -NTUD*(1.-CSTAR)
      IF (TEMP .LT. 50.0) THEN
         C     = DEXP(TEMP)
         EPSD  = (1. - C)/(1. - CSTAR*C)
      ELSE
         EPSD = 1/CSTAR
      ENDIF
      TW2D  = TW1 + EPSD*CSTAR*(TA1 - TW1)
      TA2D  = TA1 - EPSD*(TA1 - TW1)
      TS2D  = TW1 + CSTAR*NTUD/NTUI*(TA2D - TW1)
C
      DRY   = TS2D .GT. TDP
      FDRY = 1.0
      IF(.NOT. DRY) THEN
         FDRY = 0.0
C ** WET COIL **
C   DETERMINES THE AIR EFFECTIVENESS (EPSW) AND OUTLET WATER TEMPERATURE
C   (TW2W) ASSUMING THE COIL IS WET.
C   THE AVERAGE SATURATION SPECIFIC HEAT (CS) DEPENDS UPON THE OUTLET 
C   CONDITIONS; THEREFORE AN ITERATIVE SOLUTION.
C   THE WET SURFACE HEAT TRANSFER COEFFICIENT (HWO) AND THEREFORE THE FIN 
C   EFFICIENCY (EFFW) DEPEND UPON THE SATURATION SPECIFIC HEAT.
C   IF THE SURFACE TEMPERATURE AT WATER OUTLET IS LESS THAN THE DEWPOINT 
C   OF THE ENTERING AIR THEN THE COIL IS COMPLETELY WET.
         TW2W = TW2D
         ITER = 0
10       ITER = ITER + 1
         PSYDAT(2) = TW2W
         PSYDAT(3) = TW2W
         CALL PSYCHROMETRICS(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*104)
         CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
104      HW2 = PSYDAT(7)
         CS    = (HW2 - HW1)/(TW2W - TW1)
         MSTAR = CS/RA/CPW(IU)
         HWO   = CS*HDO/CPM
         ALPHA = RO/RE
         BETA  = RE*DSQRT(2.*HWO/KF/FT)
         EFF   = FIN(ALPHA,BETA)
         EFFW  = 1. - FRATIO*(1.-EFF)
         ERAT  = EFFW/EFFD
         NTUW  = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR)
         TEMP  = -NTUW*(1.-MSTAR)
         IF (TEMP .LT. 50.0) THEN
            C     = DEXP(TEMP)
            EPSW  = (1. - C)/(1. - MSTAR*C)
         ELSE
            EPSW = 1/MSTAR
         ENDIF
C  USE A SECANT METHOD TO CONVERGE ON TW2
         G=TW2W-(TW1 + EPSW*(HA1 - HW1)/RA/CPW(IU))
         IF (DABS(G) .GT. TOL .AND. ITER .LT. IMAX) THEN
            IF (ITER .EQ. 1) THEN
               G1    = G
               TOLD  = TW2W
               TW2W  = TOLD - G1
            ELSE
               G2    = G
               DGDT  = (G2 - G1)/(TW2W - TOLD)
               DGDT  = SIGN((DMAX1(1.0E-06,DABS(DGDT))),DGDT)
               TOLD  = TW2W
               TW2W  = TOLD - G2/DGDT
               IF (DABS(TW2W-TOLD) .LT. TOL) TW2W = TOLD - G2
               G1    = G2
            END IF
            GOTO 10
          END IF
          TS1W = TW2W + CSTAR/CPM*NTUW/NTUI*(HA1 - HW2)
      ENDIF
C  CHECK TO SEE IF COIL COMPLETELY WET
      WET = ((.NOT. DRY) .AND. TS1W .LT. TDP)
      IF(.NOT.(DRY) .AND. .NOT.(WET)) THEN
         IF (MODE .EQ. 1) THEN
            DRY = TW2D .GT. TW2W
            FDRY = -1.0
         ELSE IF (MODE .EQ. 2) THEN
C ** COMBINED WET AND DRY COIL PERFORMANCE **
C
C  AN ITERATIVE BY-SECTION METHOD IS USEDD TO FIND TW2. A FUNCTION
C  G, WHICH EQUALS ZERO AT THE SOLUTION, IS DEFINED FROM AN ENERGY
C  BALANCE.  THE VALUE OF G IS FOUND FOR TW2 CORRESPONDING TO FDRY=0
C  AND FDRY=1.  IF BOTH VALUES OF G ARE POSITIVE OR BOTH ARE NEGATIVE,
C  THE LOGICAL VARIABLES "WET" AND "DRY" ARE APPROPRIATELY SET.  (IE.
C  THE COIL WILL BE CONSIDERED WET IF THE WET LIMIT GIVES THE G VALUE
C  CLOSEST TO 0, OR THE COIL WILL BE CONSIDERED DRY IF THE DRY LIMIT 
C  GIVES THE G VALUE CLOSEST TO 0.)
            CK = NTUD*(1. - CSTAR)
            LIMITW=(TA1*(CSTAR-1.+CK/NTUO)+TDP*(1-CSTAR))/(CK/NTUO)
            LIMITD=(TA1*(CSTAR-(1.-CK/NTUO)*DEXP(-CK))+TDP*(1.-CSTAR))/
     &             (1.-(1.-CK/NTUO)*DEXP(-CK))
            IF(LIMITD .GT. LIMITW) THEN
               TMAX = 0.999999 * LIMITD
               TMIN = 1.000001 * LIMITW
               MAXDRY = .TRUE.
            ELSE
               TMAX = 0.999999 * LIMITW
               TMIN = 1.000001 * LIMITD
               MAXDRY = .FALSE.
            ENDIF
            CALL CALCG(TMAX,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
     &                 NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,INFO,
     &                 IU,GMAX,FDRY,EPSD,TWX,ERAT,EMODE)
            CALL CALCG(TMIN,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
     &                 NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,INFO,
     &                 IU,GMIN,FDRY,EPSD,TWX,ERAT,EMODE)
            IF ((GMAX .GT. 0.0 .AND. GMIN .GT. 0.0) .OR.
     .          (GMAX .LT. 0.0 .AND. GMIN .LT. 0.0)) THEN
               IF (DABS(GMAX) .LT. DABS(GMIN) ) THEN
                  DRY = MAXDRY
                  WET = .NOT. MAXDRY
               ELSE
                  WET = MAXDRY
                  DRY = .NOT. MAXDRY
               ENDIF
            ELSE
               TW2 = TMIN
               ITER = -1
C  ITERATION LOOP TO DETERMINE TW2.
15             ITER = ITER + 1
               CALL  CALCG(TW2,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
     &                     NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,
     &                     INFO,IU,G,FDRY,EPSD,TWX,ERAT,EMODE)
               IF ((DABS(G) .GT. TOL .AND. ITER .LT. IMAX) .AND.
     .             (DABS(TMIN - TMAX) .GT. TOL)) THEN
                  IF (ITER .EQ. 0) THEN
                     SGNMIN = SIGN(1.,G)
C  USE THE MAXIMUM TW2 FROM THE WET AND DRY ANALYSIS AS AN INITIAL GUESS
                     TW2 = DMAX1(TMIN,DMIN1(DMAX1(TW2W,TW2D),TMAX))
                     GOTO 15
                  ELSE IF (ITER .EQ. 1) THEN
                     G1    = G
                     TOLD  = TW2
                     TW2   = TOLD - G1
                  ELSE
                     G2    = G
                     DGDT  = (G2 - G1)/(TW2 - TOLD)
                     DGDT  = SIGN((DMAX1(1.0E-06,DABS(DGDT))),DGDT)
                     TOLD  = TW2
                     TW2   = TOLD - G2/DGDT
                     G1    = G2
                  END IF
C  NARROW THE BOUNDS TO THE SOLUTION
                  IF (G1 .LT. 0.) THEN
                     IF (SGNMIN .LT. 0.) THEN
                        TMIN = TOLD
                     ELSE
                        TMAX = TOLD
                     END IF
                  ELSE
                     IF (SGNMIN .GT. 0.) THEN
                        TMIN = TOLD
                     ELSE
                        TMAX = TOLD
                     END IF
                  END IF
C  USE A BI-SECTION IF THE NEW VALUE OF TW2 IS OUTSIDE THE BOUNDS
                  IF (TW2.GE.TMAX.OR.TW2.LE.TMIN) TW2 = (TMAX + TMIN)/2.
                  GOTO 15
               END IF
            END IF
C * END OF COMBINED COIL ANALYSIS
         END IF
      END IF
C ** END OF COIL ANALYSIS
C
C *** NET COIL PERFORMANCE ***
C   FOR MODE 1, UTILIZE THE RESULTS OF THE ANALYSIS THAT GIVES THE
C   GREATEST HEAT TRANSFER (I.E. HIGHEST WATER RETURN TEMPERATURE) 
C   IF COIL IS NOT COMPLETELY WET OR DRY
C ** DRY **
      IF(DRY) THEN
         TW2  = TW2D
         TA2  = TA2D
         HA2  = HA1 - CPM*(TA1 - TA2)
         WA2  = WA1
         RHA2 = RHA1
      ELSE IF ((WET) .OR. MODE .EQ. 1) THEN
C ** WET **
         TW2   = TW2W
         HA2   = HA1 - EPSW*(HA1 - HW1)
         HSAVG = HA1 - (HA1 - HA2)/(1. - DEXP(-ERAT*NTUO))
         IF ((9.473.GT.HCONV(HSAVG) .OR. HCONV(HSAVG).GT.355.137)
     .       .AND. (TIME-OUT(11)) .GT. DELT/2.) THEN
	      CALL MESSAGES(-1,'THE CORRELATION FOR THE SATURATION TEMPERA
     &TURE USED IN THE COOLING COIL MODEL WAS USED WITH AN ENTHALPY OUTS
     &IDE OF ITS INTENDED RANGE','WARNING',INFO(1),INFO(2))
            IWARN=IWARN+1
            OUT(11) = TIME
         END IF
         TSAVG = TSAT(HSAVG)
         TA2  = TSAVG + (TA1 - TSAVG)*DEXP(-ERAT*NTUO)
         PSYDAT(2) = TA2
         PSYDAT(7) = HA2
         CALL PSYCHROMETRICS(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*105)
         CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
105      WA2  = PSYDAT(6)
         RHA2 = PSYDAT(4)
      ELSE
C ** COMBINED **
         HA2 = HA1 - RA*CPW(IU)*(TW2 - TW1)
         TAX = TA1 - EPSD*(TA1 - TWX)
         HAX2 = HA1 - RA*CPW(IU)*(TW2 - TWX)
         HAX = HA1 - CPM*(TA1 - TAX)
         HSW = HAX + (HA2 - HAX)/(1. - DEXP(-(1. - FDRY)*ERAT*NTUO))
         IF ((9.473.GT.HCONV(HSW) .OR. HCONV(HSW).GT.355.137)
     .          .AND. (TIME-OUT(11)) .GT. DELT/2.) THEN
	      CALL MESSAGES(-1,'THE CORRELATION FOR THE SATURATION TEMPERA
     &TURE USED IN THE COOLING COIL MODEL WAS USED WITH AN ENTHALPY OUTS
     &IDE OF ITS INTENDED RANGE','WARNING',INFO(1),INFO(2))
            IWARN=IWARN+1
            OUT(11) = TIME
         END IF
         TSW = TSAT(HSW)
         TA2 = TSW + (TAX - TSW)*DEXP(-(1.-FDRY)*ERAT*NTUO)
         PSYDAT(2) = TA2
         PSYDAT(7) = HA2
         CALL PSYCHROMETRICS(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*106)
         CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
106      WA2  = PSYDAT(6)
         RHA2 = PSYDAT(4)
      ENDIF

C-----------------------------------------------------------------------------------------------------------------------
C    SET THE OUTPUTS FROM THIS MODEL IN SEQUENTIAL ORDER AND GET OUT
C****************************** OUTPUTS *****************************
C
      QCOIL = MW*CPW(IU)*(TW2 - TW1)
      QLAT  = MA*(WA1 - WA2)*HFG(IU)
      QSENS = QCOIL - QLAT
C      
500   OUT(1) = TA2
      OUT(2) = WA2
      OUT(3) = MA
      OUT(4) = TW2
      OUT(5) = MW
      OUT(6) = QCOIL
      OUT(7) = QSENS
      OUT(8) = QLAT
      OUT(9) = FDRY
	OUT(10) = 100.*DMIN1(1.,(DMAX1(0.,RHA2)))

C-----------------------------------------------------------------------------------------------------------------------
C    EVERYTHING IS DONE - RETURN FROM THIS SUBROUTINE AND MOVE ON
      RETURN 1
      END
C-----------------------------------------------------------------------------------------------------------------------
C  THIS SUBROUTINE CALCULATES G (WHICH IS BASED ON AN ENERGY BALANCE ON THE
C    WET SIDE OF THE COIL), FDRY AND EPSD.
C    SEE VARIABLE DEFINITIONS IN THE MAIN PROGRAM.  CPW IS ASSUMED TO BE
C    PASSED IN THE CORRECT UNITS.
C
      SUBROUTINE CALCG(TW2,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
     .                 NTUI,NTUD,HDO,CPM,CSTAR,CPW,EFFD,FRATIO,
     .                 INFO,IU,G,FDRY,EPSD,TWX,ERAT,EMODE)
	IMPLICIT NONE
      INTEGER*4 INFO
      DIMENSION PSYDAT(9),INFO(15)
      DOUBLE PRECISION KF,NTUD,NTUW,NTUO,NTUI,MSTAR,PATM,EXPK,TDP,TW2,
     & CSTAR,TA1,FDRY,TEMP,C,EPSD,TWX,TW1,PSYDAT,HWX,CS,HW1,RA,HWO,HDO,
     & CPM,ALPHA,RO,RE,BETA,FT,EFF,FIN,EFFW,FRATIO,ERAT,EFFD,CPW,EPSW,G,
     & HA1
      INTEGER EMODE,IU,ISTAT
      DATA PATM/1.0/
      EXPK = ((TDP-TW2)+CSTAR*(TA1-TDP))/
     &        (1.-NTUD*(1.-CSTAR)/NTUO)/(TA1-TW2)
      EXPK = DMAX1(1.0E-6,EXPK)
      FDRY = -1./NTUD/(1.-CSTAR)*DLOG(EXPK)
      FDRY = DMAX1(0.D0,DMIN1(1.D0,FDRY))
      TEMP = -FDRY*NTUD*(1. - CSTAR)
      IF (TEMP .LT. 50.0) THEN
         C     = DEXP(TEMP)
         EPSD  = (1. - C)/(1. - CSTAR*C)
      ELSE
         EPSD = 1/CSTAR
      ENDIF
      TWX  = (TW2 - CSTAR*EPSD*TA1)/(1 - CSTAR*EPSD)
      TWX  = DMAX1(TWX,TW1+1.E-04)
      PSYDAT(1) = PATM
      PSYDAT(2) = TWX
      PSYDAT(3) = TWX
      CALL PSYCHROMETRICS(1.0D0,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*107)
      CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
107   HWX  = PSYDAT(7)
      CS   = (HWX - HW1)/(TWX - TW1)
      MSTAR= CS/RA/CPW
      HWO  = CS*HDO/CPM
      ALPHA  = RO/RE
      BETA   = RE*DSQRT(2.*HWO/KF/FT)
      EFF = FIN(ALPHA,BETA)
      EFFW = 1. - FRATIO*(1.-EFF)
      ERAT = EFFW/EFFD
      NTUW = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR)
      TEMP = -(1. - FDRY)*NTUW*(1. - MSTAR)
      IF (TEMP .LT. 50.0) THEN
         C    = DEXP(TEMP)
         EPSW = (1. - C)/(1. - MSTAR*C)
      ELSE
         EPSW = 1/MSTAR
      ENDIF
      G=TWX*(1.-CSTAR*EPSW*EPSD)-TW1-CSTAR/CPM*EPSW*(HA1-EPSD*CPM*TA1
     &  -HW1)
      RETURN
      END
C******************************************************************************
C  THIS FUNCTION CALCULATES THE FIN EFFICIENCY (EFFECTIVENESS)
C  OF AN ANNULAR FIN OF CONSTANT THICKNESS.
C
C   ALPHA = RADIUS AT FIN BASE / RADIUS AT FIN TIP
C   BETA  = RADIUS AT FIN TIP *
C           (SQRT (2 * CONVECTION COEFFICIENT / 
C                  FIN CONDUCTIVITY * FIN THICKNESS))
      FUNCTION FIN(ALPHA,BETA)
	   IMPLICIT NONE
         DOUBLE PRECISION FIN,I0,I1,K0,K1,ALPBET,ALPHA,BETA,XI0,XI1,XK0,
     & XK1,YI0,YI1,YK0,YK1
         ALPBET = ALPHA * BETA
         CALL BESSEL(ALPBET,I0,I1,K0,K1)
         XI0 = I0
         XI1 = I1
         XK0 = K0
         XK1 = K1
         CALL BESSEL(BETA,I0,I1,K0,K1)
         YI0 = I0
         YI1 = I1
         YK0 = K0
         YK1 = K1
         FIN=2.*ALPHA/BETA/(1.- ALPHA**2)*(XK1*YI1-XI1*YK1)/(XK0*YI1
     &       +XI0*YK1)
         RETURN
      END
C****************************************************************************
C  THIS SUBROUTINE USES POLYNOMIAL APPROXIMATIONS TO EVALUATE
C  THE BESSEL FUNCTIONS.  THE APPROXIMATIONS ARE FROM ABRAMOWITZ
C  AND STEGUN, HANDBOOD OF MATHEMATICAL FUNCTIONS, DOVER 
C  PUBLICATIONS, INC., NEW YORK, NY.
C
      SUBROUTINE BESSEL(X,I0,I1,K0,K1)
	IMPLICIT NONE
      DOUBLE PRECISION X,I0,I1,K0,K1,IT,A0,A1,A2,A3,A4,A5,A6,B0,B1,B2,
     & B3,B4,B5,B6,B7,B8,C0,C1,C2,C3,C4,C5,C6,D0,D1,D2,D3,D4,D5,D6,D7,
     & D8,E0,E1,E2,E3,E4,E5,E6,F0,F1,F2,F3,F4,F5,F6,G0,G1,G2,G3,G4,G5,
     & G6,H0,H1,H2,H3,H4,H5,H6,T,TT,X1,X2
C  THE FOLLOWING DATA STATEMENTS CONTAIN THE COEFFICIENTS TO 
C  THE POLYNOMIALS.
C  I0 
      DATA A0/1.0/,A1/3.5156229/,A2/3.0899424/,A3/1.2067492/
      DATA A4/0.2659732/,A5/0.0360768/,A6/0.0045813/
C  I0
      DATA B0/0.39894228/,B1/0.01328592/,B2/0.00225319/
      DATA B3/-0.00157565/,B4/0.00916281/,B5/-0.02057706/
      DATA B6/0.02635537/,B7/-0.01647633/,B8/0.00392377/
C  I1
      DATA C0/0.5/,C1/0.87890594/,C2/0.51498869/,C3/0.15084934/
      DATA C4/0.02658733/,C5/0.00301532/,C6/0.00032411/
C  I1
      DATA D0/0.39894228/,D1/-0.03988024/,D2/-0.00362018/
      DATA D3/0.00163801/,D4/-0.01031555/,D5/0.02282967/
      DATA D6/-0.02895312/,D7/0.01787654/,D8/-0.00420059/
C  K0
      DATA E0/-0.57721566/,E1/0.4227842/,E2/0.23069756/
      DATA E3/0.0348859/,E4/0.00262698/,E5/0.0001075/,E6/0.0000074/
C  K0
      DATA F0/1.25331414/,F1/-0.07832358/,F2/0.02189568/
      DATA F3/-0.01062446/,F4/0.00587872/,F5/-0.0025154/
      DATA F6/0.00053208/
C  K1
      DATA G0/1.0/,G1/0.15443144/,G2/-0.67278579/,G3/-0.18156897/
      DATA G4/-0.01919402/,G5/-0.00110404/,G6/-0.00004686/
C  K1
      DATA H0/1.25331414/,H1/0.23498619/,H2/-0.0365562/
      DATA H3/0.01504268/,H4/-0.00780353/,H5/0.00325614/
      DATA H6/-0.00068245/
C
      IF (X .LT. -3.75) THEN
	   CALL MESSAGES(-1,'THE BESSEL FUNCTION CALLED FROM THE COOLING C
     &OIL SUBROUTINE COULD NOT BE EVALUATED AT THE GIVEN VALUE','FATAL',
     &-1,-1)
         RETURN
      END IF
      T=X/3.75
      TT=T*T
C I0
      IF (X .LE. 3.75) THEN
         I0=A0+TT*(A1+TT*(A2+TT*(A3+TT*(A4+TT*(A5+TT*A6)))))
      ELSE
         IT=1/T
         I0=(B0+IT*(B1+IT*(B2+IT*(B3+IT*(B4+IT*(B5+IT*(B6+IT*
     .      (B7+IT*B8))))))))/(DSQRT(X)*DEXP(-X))
      END IF
C  I1
      IF (X .LE. 3.75) THEN
         I1=(C0+TT*(C1+TT*(C2+TT*(C3+TT*(C4+TT*(C5+TT*C6))))))*X
      ELSE
         IT=1/T
         I1=(D0+IT*(D1+IT*(D2+IT*(D3+IT*(D4+IT*(D5+IT*(D6+IT*
     .      (D7+IT*D8))))))))/(DSQRT(X)*DEXP(-X))
      END IF
C  K0
      IF (X .LE. 0.0) THEN
	   CALL MESSAGES(-1,'THE BESSEL FUNCTION CALLED FROM THE COOLING C
     &OIL SUBROUTINE COULD NOT BE EVALUATED AT THE GIVEN VALUE','FATAL',
     &-1,-1)
         RETURN
      END IF
      X1 = (X/2.)**2
      X2 = 2./X
      IF (X .LE. 2.0) THEN
         K0=-DLOG(X/2)*I0+E0+X1*(E1+X1*(E2+X1*(E3+X1*(E4+X1*
     .      (E5+X1*E6)))))
      ELSE
         K0=(F0+X2*(F1+X2*(F2+X2*(F3+X2*(F4+X2*(F5+X2*F6))))))
     .      /(DSQRT(X)*DEXP(X))
      END IF
C  K1
      IF (X .LE. 2.0) THEN
         K1=(X*DLOG(X/2.)*I1+G0+X1*(G1+X1*(G2+X1*(G3+X1*(G4+X1*
     .      (G5+X1*G6))))))/X
      ELSE
         K1=(H0+X2*(H1+X2*(H2+X2*(H3+X2*(H4+X2*(H5+X2*H6))))))
     .      /(DSQRT(X)*DEXP(X))
      END IF
      RETURN
      END