C A FORTRAN PROGRAM TO COMPUTE SURFACE HEAT EXCHANGE AND EQUILIBRIUM
C TEMPERATURE FROM DAILY OBSERVATIONS  FEB 13, 2002

      DATA PI /3.1416/, SOLAR_CONST /1395/

      REAL LATITUDE, KE 

      OPEN(10,FILE='heat_exchange.npt',STATUS='OLD')
      OPEN(11,FILE='heat_exchange.opt',STATUS='UNKNOWN')

C CONVERT LATITUDE IN DEGREES TO RADIANS

      LATITUDE = 2. * PI * (38.6)/360.

C INITIAL ESTIMATES OF PARAMETERS NEEDED TO ITERATE KE AND TE

      TENEW = 8.3
      
 1    READ (10,*,END=999) YEAR, MONTH, DAY, D_OF_YR, 
     $  DRY_BULB, WET_BULB, DEW_POINT, WIND, SKY_COVER

C MULTIPLY WIND BY 1.5 TO ACCOUNT FOR LARGER WIND OVER WATER
	WIND = WIND * 1.5

C DAILY AVERAGE SOLAR RADIATION FROM COMPUTED HOURLY VALUES
C COMPUTATIONS BASED ON TVA LAB REPORT 14, PAGES 2.2 - 2.7

      DECL = 23.45*PI/180. * COS(2.*PI*(172. - D_OF_YR)/365.)
      RAD_VEC = 1. + 0.017*COS(2.*PI*(186. - D_OF_YR)/365.)
      SUMRAD = 0

      DO J=1,24

        HOUR = FLOAT(J)-0.5
        TAU = 2.*PI*(HOUR-12.)/24.
        SIN_ALPH = SIN(LATITUDE)*SIN(DECL) + 
     $    COS(LATITUDE)*COS(DECL)*COS(TAU)
        HS = SOLAR_CONST*SIN_ALPH/RAD_VEC/RAD_VEC
        HS = MAX(HS,0.)
        SUMRAD = SUMRAD + HS

      END DO
      HS = SUMRAD/24.

C CORRECT FOR REFLECTION AND EFFECT OF CLOUDS AS PER R&H PAGE 31
C SKY COVER IS REPORTED IN EIGHTHS.  CONVERT TO FRACTION.

      SKY_COVER = MIN(SKY_COVER,8.)/8.
      HS = 0.94 * HS * (1.-0.65*SKY_COVER*SKY_COVER)

C CONVERT WATTS TO EINSTEINS USING MY OWN EMPIRICAL FORMULA
      PAR = 0.143 * HS
 
C FRACTIONAL DAYLENGTH AS PER TVA LAB REPORT 14, PAGE 2.10

      FD = ACOS(-SIN(LATITUDE)*SIN(DECL)/COS(LATITUDE)/COS(DECL))
     $  / PI

C ATMOSPHERIC RADIATION
C USE SWINBANK'S FORMULA AS PER TVA LAB REPORT 14, PAGE 3.18
C THEIR UNITS ARE KJ/M**2/HR.  0.278 CONVERTS TO WATTS/M**2
C ASSUME 3% REFLECTION  

      TABS = DRY_BULB +273.15
      HA = 0.937E-5 * 2.0411E-7 * TABS**6 *
     $  (1.+0.17*SKY_COVER*SKY_COVER)
      HA = 0.97 * 0.278 * HA

C TO USE DEWPOINT AS INITIAL GUESS AT TE, UN-COMMENT THE FOLLOWING LINE
C      TENEW = DEW_POINT

C ITERATE UNTIL TE CHANGES BY LESS THAN 0.1 C
C USE FORMULAE IN APPENDIX K OF EDINGER, BRADY, & GEYER

      FW = 9.2 + 0.46 * WIND * WIND

      DO I=1,100

        TE = TENEW
        TM = (TE + DEW_POINT)/2.
        BETA = 0.35 + 0.015 * TM + 0.0012 * TM * TM
        HN = HS + HA

        KE = 4.48 + (BETA + 0.47) * FW + 0.05 * TE
        TDSTAR = DEW_POINT + 0.47 * (DRY_BULB - DEW_POINT) / 
     $    (BETA + 0.47)
        TENEW = (HN - 306. + (KE - 4.48) * TDSTAR) / 
     $    (KE + 0.05 * TDSTAR - 0.025 * TE)
        IF (ABS(TENEW - TE) .LT. 0.1) GO TO 2
        
      END DO

 2    WRITE(11,3) D_OF_YR -1, KE, TENEW, PAR, FD, WIND
 3    FORMAT(6F8.2)

      GO TO 1

 999  END
  
            
      
