       SUBROUTINE flux_table_hur(u,v,sst,at,ref_ht,Hsig,u_orb,v_orb,qs,
     +                           q,press,tau_u,tau_v,shf,lhf)



C      DEVELOPED BY: Dr. Mark A. Bourassa
C                    Center for Ocean-Atmospheric Prediction Studies
C                    Florida State University
C                    Tallahassee, FL 32306-2840
C
C                    Email: bourassa@coaps.fsu.edu
C                    Phone: (850) 644-6923

                     
C                    Paul Hughes
C                    Center for Ocean-Atmospheric Prediction Studies
C                    Florida State University
C                    Tallahassee, FL 32306-2840 

C                    Email: hughes@coaps.fsu.edu

C      VERSION: 1.0

C      MODIFICATIONS:

C      Date       Programmer      Description of Change 
C      ====       ==========      =====================
C      12/19/05    P. Hughes       Original Code


C      PURPOSE: Subroutine computes wind stress (u and v components), 
C               sensible heat flux, and latent heat flux via the bulk flux 
C               formulas.  Transfer coefficients calculated using a 
C               surface flux model (adopted from Bourassa 2005)


C      INPUT VARIABLES: 'u' component of wind speed  [m/s]
C                       'v' component of wind speed  [m/s]
C                       sea surface temperature  [Celsius]
C                       air temperature  [Celsius]
C                       reference height  [m]
C                       significant wave height  [m]
C                       'u' orbital velocity  [m/s]
C                       'v' orbital velocity  [m/s]
C                       specific humidity (surface)  [kg/kg]
C                       specific humidity (reference height)  [kg/kg]
C                       surface pressure  [Pa]

C      OUTPUT VARIABLES: 'u' component of stress  [N/m^2]
C                        'v' component of stress  [N/m^2]
C                        sensible heat flux  [W/m^2]
C                        latent heat flux  [W/m^2] 


C      RANGE OF FLUX TABLE: 

C      wind speed (0.0 < wind speed < 90.0 m/s)
C      SST-Tair (-20.0 < SST-Tair < 20.0 degress C)
C      reference height: (10.0 <= adjusted referecne height < 300.0 m)

C     *** If additional ranges needed please contact ***

 
C      ADDITIONAL COMMENTS:

C     (1) Wind speed adjusted using orbital velocity
C         Uorb=pi*Hsig/Peak Period
C         adjusted wind speed = wind speed - (0.8*Uorb)

C     (1b) For more accurate fluxes the currents (if applicable) should
C          also be subtracted from the wind speed
C          adjusted wind speed = windspd -(0.8*Uorb) - current

C     (2) Reference height adjusted using significant wave height
C         adjusted reference height = ref_ht - (0.8*Hsig)

C     (3) If no wave information then set significant wave height,
C         'u' and 'v' orbital velocity equal to zero

C     (4) If using total wind speed (orbital velocity) instead of 
C         vector components then set 'v' ('v_orb') equal to 0.0 and 
C         'u' ('u_orb') equal to the total wind speed (orbital velocity).
C         The stress will be output in the 'tau_u' variable

C     (5) The surface flux model used to determine transfer coefficients
C         is a combination of Bourassa et al.(1999), Clayson et al.(1996),
C         and Bourassa (2005) 

C     (6) Postive sensible and latent heat flux => ocean to atmoshpere
C         Postive Stress => atmoshpere to ocean

C     (7) Validation of flux table in a paper in preparation


C      REFERENCES:

C      Bourassa, M.A., D.G. Vincent, and W.L. Wood, 1999: A flux 
C         parameterization including the effects of capillary waves
C         and sea state. J. Atmos. Sci., 56, 1123-1139.


C      Bourassa, M.A., 2005, Satellite-based observations of surface
C         turbulent stress during severe weather, Atmosphere-Ocean 
C         Interactions, Vol. 2., ed., W. Perrie, Wessex Institute of 
C         Technology, in press.

C      Clayson, C.A., C.W. Fairall, and J.A. Curry, 1996: Evaluation of 
C         turbulent fluxes at the ocean surface using surface renewal
C         theory. J. Geophys. Res., 101, 28503-28513.        


       IMPLICIT NONE

       include "flux_hur.inc"


       INTEGER i, i_index, i_index2, i_u, i_v
       INTEGER j, j_index, j_index2, j_index3
       INTEGER k, k_index, k_index2     
       INTEGER u_index, u_index2
       INTEGER v_index, v_index2 
       REAL cd_uin(2,2,2), cd_vin(2,2,2)
       REAL ch_in(2,2,2), ce_in(2,2,2)
       REAL x_in(2), x_uin(2), x_vin(2), y_in(2), z_in(2) 
       REAL cdu, cdv, ch, ce 
       REAL u, v, sst, at, ref_ht, Hsig, u_orb, v_orb, qs, q, press
       REAL tau_u, tau_v, shf, lhf
       REAL deltaT,deltaq, adj_ref_ht, spdu, spdv, spd
       REAL step_spd, step_temp, step_ht
       REAL fast_trilin
       REAL rho, Cp, Lv, Rd, tvirt
       REAL temp_u, temp_v

       temp_u = 0.0
       temp_v = 0.0

       step_spd = 3.0
       step_temp = 1.0
       step_ht = 10.0


C      fast_trilin: trilinear interpolation function
C      cd_uin :  values of drag coefficient sent to fast_trilin
C      cd_vin :  values of drag coefficient sent to fast_trilin
C      ch_in :  values of heat transfer coefficient sent to fast_trilin
C      ce_in :  values of moisture transfer coefficient sent to fast_trilin
C      x_in: two values of wind speed sent to fast_trilin
C      x_uin: two values of 'u' component wind speed sent to fast_trilin
C      x_vin: two values of 'v' component wind speed sent to fast_trilin
C      y_in: two values of sea/air temp difference sent to fast_trilin
C      z_in: two values or adjusted referecne height sent to fast_trilin
C      u: 'u' component of wind speed (input variable)  [m/s]
C      v: 'v' component of wind speed (input variable)  [m/s]
C      at: air temperature at reference height (input variable)  [Celsius]
C      sst: sea surface temperature (input variable) [Celsius]
C      deltaT: sea surface temp - air temp  [Celsius]
C      ref_ht: reference height (input variable)  [m]
C      Hsig: significant wave height (input variable) [m]
C      u_orb: 'u' component of orbital velocity (input variable) [m/s]
C      v_orb: 'v' component of orbital velocity (input variable) [m/s]
C      qs: specific humidity at surface (input variable) [kg/kg]
C      q:  specific humidity at reference height (input variable) [kg/kg]
C      deltaq : qs-q  [kg/kg] 
C      press: surface pressure (input variable)  [Pa]
C      adj_ref_ht: adjusted reference height  [m]
C      spdu: adjusted 'u' component wind speed  [m/s]
C      spdv: adjusted 'v' component wind speed  [m/s]
C      spd: adjusted wind speed  [m/s]
C      cdu: 'u' component drag coefficient returned from fast_trilin
C      cdv: 'v' component drag coefficient returned from fast_trilin
C      ch: heat transfer coefficient returned from fast_trilin
C      ce: moisture transfer coefficient returned from fast_trilin
C      tau_u: 'u' component stress (output)  [N/m^2]
C      tau_v: 'v' component stress (output)  [N/m^2]
C      shf: sensible heat flux (output)  [W/m^2]
C      lhf: latent heat flux (output)  [W/m^2]
C      rho: density of air  [kg/m^3]
C      Cp: heat capcity of air  [J/kg K]
C      Lv: latent heat of vaporization [J/kg]
C      Rd: gas constant  [J/K kg]
C      tvirt: virtual temperature  [Kelvin]
C      step_spd: increment between windspeed values
C      step_temp: increment between tdif values 
C      step_ht: increment between adj_ht values 



       deltaT = sst - at
       deltaq = qs - q
       adj_ref_ht = ref_ht - (0.8*Hsig)
       spdu = u - (0.8*u_orb)
       spdv = v - (0.8*v_orb)
       spd = SQRT(spdu**2 + spdv**2)

       IF (spdu .lt. 0.0) temp_u=spdu
       IF (spdv .lt. 0.0) temp_v=spdv       

       i = INT(spd/step_spd)
       i_u = INT(spdu/step_spd)
       i_v = INT(spdv/step_spd)
       k = INT(adj_ref_ht/step_ht)
       j = INT(deltaT/step_temp)


       IF ((i .ge. 0) .AND. (i .lt. 30)) THEN

           i_index = i+1
           i_index2 = i_index+1
           x_in(1) = spd_table(i_index)
           x_in(2) = spd_table(i_index2)
           
       ELSE

           i_index = 31
           spd = spd_table(i_index)

       END IF

      

       IF ((i_u .ge. 0) .AND. (i_u .lt. 30)) THEN

          u_index = i_u+1
          u_index2 = u_index+1
          x_uin(1) = spd_table(u_index)
          x_uin(2) = spd_table(u_index2)

       ELSE IF ((i_u .lt. 0) .AND. (ABS(i_u) .lt. 30)) THEN

C      Accounts for situations when the orbital velocity is greater 
C      than the u-component of the wind speed.  This results in a 
C      negative stress or a flux of momentum from the ocean to the 
C      atmosphere        
            
          u_index = ABS(i_u)+1
          u_index2 = u_index+1
          x_uin(1) = spd_table(u_index)
          x_uin(2) = spd_table(u_index2)

       ELSE

          u_index = 31
          spdu = spd_table(u_index)

       END IF 



       IF ((i_v .ge. 0) .AND. (i_v .lt. 30)) THEN

          v_index = i_v+1
          v_index2 = v_index+1
          x_vin(1) = spd_table(v_index)
          x_vin(2) = spd_table(v_index2)

       ELSE IF ((i_v .lt. 0) .AND. (ABS(i_v) .lt. 30)) THEN

C      Accounts for situations when the orbital velocity is greater 
C      than the v-component of the wind speed.  This results in a 
C      negative stress or a flux of momentum from the ocean to the 
C      atmosphere
            
          v_index = ABS(i_v)+1
          v_index2 = v_index+1
          x_vin(1) = spd_table(v_index)
          x_vin(2) = spd_table(v_index2)

       ELSE

          v_index = 31
          spdv = spd_table(v_index)

       END IF 



       IF ((k .ge. 1) .AND. (k .lt. 30)) THEN 

          k_index = k
          k_index2 = k_index+1
          z_in(1) = ref_ht_table(k_index)
          z_in(2) = ref_ht_table(k_index2)
          
       ELSE IF (k .ge. 30) THEN

          k_index = 30

       ELSE

          k_index = 1

       END IF



       IF ((j .gt. -20) .AND. (j .lt. 20)) THEN

          j_index = j+21
          j_index2 = j_index-1
          j_index3 = j_index+1

       ELSE IF (j .ge. 20) THEN

          j_index = 41
          deltaT = tdif_table(j_index)

       ELSE

          j_index = 1
          deltaT = tdif_table(j_index)

       END IF 


       IF ((i .ge. 0) .AND. (i .lt. 30) .AND. (ABS(i_u) .lt. 30)
     +    .AND. (ABS(i_v) .lt. 30) .AND. (k .ge. 1) 
     +    .AND. (k .lt. 30) .AND. (j .gt. -20) .AND. 
     +    (j .lt. 20)) THEN

          IF (deltaT .lt. 0.0) THEN

             y_in(1) = tdif_table(j_index)
             y_in(2) = tdif_table(j_index2)
             cd_uin(1,1,1) = cd_table(u_index,j_index,k_index)
             cd_uin(2,1,1) = cd_table(u_index2,j_index,k_index)
             cd_uin(1,2,1) = cd_table(u_index,j_index2,k_index)
             cd_uin(2,2,1) = cd_table(u_index2,j_index2,k_index)
             cd_uin(1,1,2) = cd_table(u_index,j_index,k_index2)
             cd_uin(2,1,2) = cd_table(u_index2,j_index,k_index2)
             cd_uin(1,2,2) = cd_table(u_index,j_index2,k_index2)
             cd_uin(2,2,2) = cd_table(u_index2,j_index2,k_index2)

             cd_vin(1,1,1) = cd_table(v_index,j_index,k_index)
             cd_vin(2,1,1) = cd_table(v_index2,j_index,k_index)
             cd_vin(1,2,1) = cd_table(v_index,j_index2,k_index)
             cd_vin(2,2,1) = cd_table(v_index2,j_index2,k_index)
             cd_vin(1,1,2) = cd_table(v_index,j_index,k_index2)
             cd_vin(2,1,2) = cd_table(v_index2,j_index,k_index2)
             cd_vin(1,2,2) = cd_table(v_index,j_index2,k_index2)
             cd_vin(2,2,2) = cd_table(v_index2,j_index2,k_index2)

             ch_in(1,1,1) = ch_table(i_index,j_index,k_index)
             ch_in(2,1,1) = ch_table(i_index2,j_index,k_index)
             ch_in(1,2,1) = ch_table(i_index,j_index2,k_index)
             ch_in(2,2,1) = ch_table(i_index2,j_index2,k_index)
             ch_in(1,1,2) = ch_table(i_index,j_index,k_index2)
             ch_in(2,1,2) = ch_table(i_index2,j_index,k_index2)
             ch_in(1,2,2) = ch_table(i_index,j_index2,k_index2)
             ch_in(2,2,2) = ch_table(i_index2,j_index2,k_index2)

             ce_in(1,1,1) = ce_table(i_index,j_index,k_index)
             ce_in(2,1,1) = ce_table(i_index2,j_index,k_index)
             ce_in(1,2,1) = ce_table(i_index,j_index2,k_index)
             ce_in(2,2,1) = ce_table(i_index2,j_index2,k_index)
             ce_in(1,1,2) = ce_table(i_index,j_index,k_index2)
             ce_in(2,1,2) = ce_table(i_index2,j_index,k_index2)
             ce_in(1,2,2) = ce_table(i_index,j_index2,k_index2)
             ce_in(2,2,2) = ce_table(i_index2,j_index2,k_index2)

          ELSE

             y_in(1) = tdif_table(j_index)
             y_in(2) = tdif_table(j_index3)
             cd_uin(1,1,1) = cd_table(u_index,j_index,k_index)
             cd_uin(2,1,1) = cd_table(u_index2,j_index,k_index)
             cd_uin(1,2,1) = cd_table(u_index,j_index3,k_index)
             cd_uin(2,2,1) = cd_table(u_index2,j_index3,k_index)
             cd_uin(1,1,2) = cd_table(u_index,j_index,k_index2)
             cd_uin(2,1,2) = cd_table(u_index2,j_index,k_index2)
             cd_uin(1,2,2) = cd_table(u_index,j_index3,k_index2)
             cd_uin(2,2,2) = cd_table(u_index2,j_index3,k_index2)

             cd_vin(1,1,1) = cd_table(v_index,j_index,k_index)
             cd_vin(2,1,1) = cd_table(v_index2,j_index,k_index)
             cd_vin(1,2,1) = cd_table(v_index,j_index3,k_index)
             cd_vin(2,2,1) = cd_table(v_index2,j_index3,k_index)
             cd_vin(1,1,2) = cd_table(v_index,j_index,k_index2)
             cd_vin(2,1,2) = cd_table(v_index2,j_index,k_index2)
             cd_vin(1,2,2) = cd_table(v_index,j_index3,k_index2)
             cd_vin(2,2,2) = cd_table(v_index2,j_index3,k_index2)

             ch_in(1,1,1) = ch_table(i_index,j_index,k_index)
             ch_in(2,1,1) = ch_table(i_index2,j_index,k_index)
             ch_in(1,2,1) = ch_table(i_index,j_index3,k_index)
             ch_in(2,2,1) = ch_table(i_index2,j_index3,k_index)
             ch_in(1,1,2) = ch_table(i_index,j_index,k_index2)
             ch_in(2,1,2) = ch_table(i_index2,j_index,k_index2)
             ch_in(1,2,2) = ch_table(i_index,j_index3,k_index2)
             ch_in(2,2,2) = ch_table(i_index2,j_index3,k_index2)

             ce_in(1,1,1) = ce_table(i_index,j_index,k_index)
             ce_in(2,1,1) = ce_table(i_index2,j_index,k_index)
             ce_in(1,2,1) = ce_table(i_index,j_index3,k_index)
             ce_in(2,2,1) = ce_table(i_index2,j_index3,k_index)
             ce_in(1,1,2) = ce_table(i_index,j_index,k_index2)
             ce_in(2,1,2) = ce_table(i_index2,j_index,k_index2)
             ce_in(1,2,2) = ce_table(i_index,j_index3,k_index2)
             ce_in(2,2,2) = ce_table(i_index2,j_index3,k_index2)

          END IF 

C         Trilinear interpolate single value of transfer coefficient
    
          cdu = fast_trilin(x_uin,y_in,z_in,cd_uin,ABS(spdu),
     +                      deltaT,adj_ref_ht)

          cdv = fast_trilin(x_vin,y_in,z_in,cd_vin,ABS(spdv),
     +                      deltaT,adj_ref_ht)

          ch = fast_trilin(x_in,y_in,z_in,ch_in,spd,
     +                     deltaT,adj_ref_ht)

          ce = fast_trilin(x_in,y_in,z_in,ce_in,spd,
     +                     deltaT,adj_ref_ht)

       ELSE

C      WARNING: Wind speed (u-component, v-component, or overall) and/or
C               adjusted reference height and/or Sea-Air temperature 
C               difference outside range of table


          cdu = cd_table(u_index,j_index,k_index)
          cdv = cd_table(v_index,j_index,k_index)
          ch = ch_table(i_index,j_index,k_index)
          ce = ce_table(i_index,j_index,k_index)

       END IF 



C      Calculate fluxes 



       Rd = 287.0       
       Cp = 1004.0*(1.0 + 0.9 * qs)
       Lv = 4186.8*(597.31 - 0.56525 * sst) 
       tvirt = (at + 273.15)*(1.0 + 0.61*q)
       rho = press/(Rd*tvirt)


      tau_u = rho*cdu*(spdu)**2

      tau_v = rho*cdv*(spdv)**2 

      shf = rho*Cp*ch*(deltaT)*spd

      lhf = rho*Lv*ce*(deltaq)*spd

      IF (temp_u .lt. 0.0) tau_u=-1.0*tau_u
      IF (temp_v .lt. 0.0) tau_v=-1.0*tau_v

      
     
      RETURN
      END


C ------------------------ FUNCTION fast_trilin ------------------


       function fast_trilin(x,y,z,val,xp,yp,zp)

C      Trilinearly interpolates from one regular grid to another regular
C      grid. The regular grids allow the neighbooring points to be 
C      known, thus elimenating the search for these points. 
C
C
C      Given a location xp,yp,zp and an array val with 2 locations in 
C      the vector x, 2 locations in y, and 2 locations in the vector z
C      return val(xp,yp,zp) using trilinear interpolation



       IMPLICIT NONE 
   
       real x(2), y(2),z(2), xp, yp, zp 
       real val(2,2,2), fast_trilin      
       real t,u,v
       real y1,y2,y3,y4,y5,y6,y7,y8 
  

       t = (xp-x(1))/(x(2)-x(1))

       u = (yp-y(1))/(y(2)-y(1))

       v = (zp-z(1))/(z(2)-z(1))




       y1 = val(1,1,1)
       y2 = val(2,1,1)
       y3 = val(1,2,1)
       y4 = val(2,2,1)
       y5 = val(1,1,2)
       y6 = val(2,1,2)
       y7 = val(1,2,2)
       y8 = val(2,2,2)


       fast_trilin = y1*(1.0-t)*(1.0-u)*(1.0-v) + y2*t*(1.0-u)*(1.0-v) +
     +               y3*(1.0-t)*u*(1.0-v) + y4*t*u*(1.0-v) +
     +               y5*(1.0-t)*(1.0-u)*v + y6*t*(1.0-u)*v +
     +               y7*(1.0-t)*u*v + y8*t*u*v
  

       return
       end
