SUBROUTINE lookup_kara(tau_over_rho, t_sst, t_air, uten_out)

! PURPOSE: Subroutine determines 10 meter wind speed from a specified
!          stress/density, air temperature, and sea surface temperature
!
! INPUT VARIABLES:  (1) Stress divided by density of air (tau_over_rho)
!
!                               tau/rho = ustar^2
!
!                   (2) Sea surface temperature (t_sst) [degrees C]
!
!                   (3) Air temperature @ 10 meters (t_air) [degrees C] 
!
!
! OUTPUT VARIABLES: (1) Corresponding wind speed at 10 meters (uten_out)
!                       (m/s)
!
!
! RANGE OF LOOKUP TABLE:  0.0 <= ustar < 3.0 [increment 0.1]
!                        -15.0  < t_sst - t_air < 15.0 [increment 1.0]
!
!            *** If additional ranges needed please contact ***


USE lookup_kara_mod

IMPLICIT NONE

INTEGER :: i, i_index1, i_index2
INTEGER :: j, j_index1, j_index2, j_index3
REAL :: step1, step2
REAL, DIMENSION(2) :: x_in, y_in
REAL, DIMENSION(2,2) :: uten_in
REAL :: fast_blin, theta_sst, theta_air, theta_diff, ref_ht, g_over_cp
REAL :: u_star
REAL, INTENT(IN) :: tau_over_rho, t_sst, t_air
REAL, INTENT(OUT) :: uten_out

ref_ht = 10.0                    
g_over_cp = 0.0098 

u_star = SQRT(tau_over_rho)

theta_sst = t_sst + 273.15
theta_air = (t_air + 273.15) + g_over_cp * ref_ht
theta_diff = theta_sst - theta_air

!WRITE(*,'(A,1X,F4.1)') 'theta_diff = ', theta_diff

step1 = 0.1
step2 = 1.0

i = INT(u_star/step1)
j = INT(theta_diff/step2)


IF((i .GE. 0) .AND. (i .LT. 30)) THEN

    i_index1 = i + 1
    i_index2 = i_index1 + 1
    
    x_in(1) = ustar_lookup(i_index1)
    x_in(2) = ustar_lookup(i_index2)

!    WRITE(*,'(A,1X,F4.2)') 'x_in(1) = ', x_in(1)
!    WRITE(*,'(A,1X,F4.2)') 'x_in(2) = ', x_in(2)

ELSE

    i_index1 = 31
    i_index2 = i_index1
 
!    WRITE(*,*) 'OUTSIDE BOUNDS'

ENDIF




IF((j .GT. -15.0) .AND. (j .LT. 15.0)) THEN

    j_index1 = j + 16
    j_index2 = j_index1 + 1 
    j_index3 = j_index1 - 1

!    WRITE(*,'(A,1X,I2,1X,F4.1)') 'j_index1 = ', j_index1, tdif_lookup(j_index1)
!    WRITE(*,'(A,1X,I2,1X,F4.1)') 'j_index2 = ', j_index2, tdif_lookup(j_index2)
!    WRITE(*,'(A,1X,I2,1X,F4.1)') 'j_index3 = ', j_index3, tdif_lookup(j_index3)


ELSE IF (j .GE. 15.0) THEN 

    j_index1 = 31
    j_index2 = j_index1
    j_index3 = j_index1

!    WRITE(*,*) 'OUTSIDE BOUNDS'

ELSE

    j_index1 = 1
    j_index2 = j_index1
    j_index3 = j_index1

!    WRITE(*,*) 'OUTSIDE BOUNDS'

ENDIF


IF((i .GE. 0) .AND. (i .LT. 30) .AND. (j .GT. -15) .AND. &
   (j .LT. 15)) THEN

    IF(theta_diff .LT. 0.0) THEN

       y_in(1) = tdif_lookup(j_index3)
       y_in(2) = tdif_lookup(j_index1)

!       WRITE(*,'(A,1X,F4.1)') 'y_in(1) = ', y_in(1)
!       WRITE(*,'(A,1X,F4.1)') 'y_in(2) = ', y_in(2)

       uten_in(1,1) = uten_lookup(i_index1,j_index3)
       uten_in(2,1) = uten_lookup(i_index2,j_index3)
       uten_in(1,2) = uten_lookup(i_index1,j_index1)
       uten_in(2,2) = uten_lookup(i_index2,j_index1)

!       WRITE(*,*) 'uten_in = ', uten_in

    ELSE

       y_in(1) = tdif_lookup(j_index1)
       y_in(2) = tdif_lookup(j_index2)

!       WRITE(*,'(A,1X,F4.1)') 'y_in(1) = ', y_in(1)
!       WRITE(*,'(A,1X,F4.1)') 'y_in(2) = ', y_in(2)

       uten_in(1,1) = uten_lookup(i_index1,j_index1)
       uten_in(2,1) = uten_lookup(i_index2,j_index1)
       uten_in(1,2) = uten_lookup(i_index1,j_index2)
       uten_in(2,2) = uten_lookup(i_index2,j_index2)

!       WRITE(*,*) 'uten_in = ', uten_in

    ENDIF

!   Bilinear interploate single value of wind speed


    uten_out = fast_blin(x_in,y_in,uten_in,u_star,&
                         theta_diff)


ELSE

!   WARNING: SQRT(tau_over_rho) and/or Tdif outside range of table

    uten_out = uten_lookup(i_index1,j_index1)

END IF



RETURN
END 


!----------------------- FUNCTION fast_bliin --------------------

REAL FUNCTION fast_blin(x,y,val,xp,yp)

!Bilinearly interpolates from one regular grid to another regular grid.
!The regular grids allow the neighboring points to be known, thus 
!eliminating the search fro these points.

!Given a location xp, yp and an array val with 2 locations in the vector 
!x and 2 locations in y, return val(xp,yp) using bilinear interpolation.


IMPLICIT NONE 

REAL, INTENT(IN), DIMENSION(2) :: x,y
REAL, INTENT(IN), DIMENSION(2,2) :: val	
REAL, INTENT(IN) :: xp, yp
REAL :: t,u
REAL :: y1,y2,y3,y4


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

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

y1=val(1,1)
y2=val(2,1)
y3=val(2,2)
y4=val(1,2)


fast_blin=(1.0-t)*(1.0-u)*y1+t*(1.0-u)*y2+t*u*y3+(1.0-t)*u*y4  

return
END FUNCTION fast_blin