SUBROUTINE flux_table(u,v,sst,at,ref_ht,Hsig,u_orb,v_orb,qs,q, &
press,tau_u,tau_v,shf,lhf)
! DEVELOPED BY: Dr. Mark A. Bourassa
! Center for Ocean-Atmospheric Prediction Studies
! Florida State University
! Tallahassee, FL 32306-2840
!
! Email: bourassa@coaps.fsu.edu
! Phone: (850) 644-6923
! Paul Hughes
! Center for Ocean-Atmospheric Prediction Studies
! Florida State University
! Tallahassee, FL 32306-2840
! Email: hughes@coaps.fsu.edu
! VERSION: 1.0
! MODIFICATIONS:
! Date Programmer Description of Change
! ==== ========== =====================
! 12/19/05 P. Hughes Original Code
! PURPOSE: Subroutine computes wind stress (u and v components),
! sensible heat flux, and latent heat flux via the bulk flux
! formulas. Transfer coefficients calculated using a
! surface flux model (adopted from Bourassa 2005)
! INPUT VARIABLES: 'u' component of wind speed [m/s]
! 'v' component of wind speed [m/s]
! sea surface temperature [Celsius]
! air temperature [Celsius]
! reference height [m]
! significant wave height [m]
! 'u' orbital velocity [m/s]
! 'v' orbital velocity [m/s]
! specific humidity (surface) [kg/kg]
! specific humidity (reference height) [kg/kg]
! surface pressure [Pa]
! OUTPUT VARIABLES: 'u' component of stress [N/m^2]
! 'v' component of stress [N/m^2]
! sensible heat flux [W/m^2]
! latent heat flux [W/m^2]
! RANGE OF FLUX TABLE:
! wind speed (0.0 < wind speed < 60.0 m/s)
! SST-Tair (-20.0 < SST-Tair < 20.0 degress C)
! reference height: (2.0 <= adjusted referecne height < 40.0 m)
! *** If additional ranges needed please contact ***
! ADDITIONAL COMMENTS:
! (1) Wind speed adjusted using orbital velocity
! Uorb=pi*Hsig/Peak Period
! adjusted wind speed = wind speed - (0.8*Uorb)
! (1b) For more accurate fluxes the currents (if applicable) should
! also be subtracted from the wind speed
! adjusted wind speed = windspd -(0.8*Uorb) - current
! (2) Reference height adjusted using significant wave height
! adjusted reference height = ref_ht - (0.8*Hsig)
! (3) If no wave information then set significant wave height,
! 'u' and 'v' orbital velocity equal to zero
! (4) If using total wind speed (orbital velocity) instead of
! vector components then set 'v' ('v_orb') equal to 0.0 and
! 'u' ('u_orb') equal to the total wind speed (orbital velocity).
! The stress will be output in the 'tau_u' variable
! (5) The surface flux model used to determine transfer coefficients
! is a combination of Bourassa et al.(1999), Clayson et al.(1996),
! and Bourassa (2005)
! (6) Postive sensible and latent heat flux => ocean to atmoshpere
! Postive Stress => atmoshpere to ocean
! (7) Validation of flux table in a paper in preparation
! REFERENCES:
! Bourassa, M.A., D.G. Vincent, and W.L. Wood, 1999: A flux
! parameterization including the effects of capillary waves
! and sea state. J. Atmos. Sci., 56, 1123-1139.
! Bourassa, M.A., 2005, Satellite-based observations of surface
! turbulent stress during severe weather, Atmosphere-Ocean
! Interactions, Vol. 2., ed., W. Perrie, Wessex Institute of
! Technology, in press.
! Clayson, C.A., C.W. Fairall, and J.A. Curry, 1996: Evaluation of
! turbulent fluxes at the ocean surface using surface renewal
! theory. J. Geophys. Res., 101, 28503-28513.
USE flux
IMPLICIT NONE
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_spd, u_index, u_index2
INTEGER :: v_spd, v_index, v_index2
REAL, DIMENSION(2,2,2) :: cd_uin, cd_vin
REAL, DIMENSION(2,2,2) :: ch_in, ce_in
REAL, DIMENSION(2) :: x_in, x_uin, x_vin, y_in, z_in
REAL, INTENT(IN) :: u, v, sst, at, ref_ht, Hsig, u_orb, v_orb, qs, q
REAL, INTENT(IN) :: press
REAL, INTENT(OUT) :: tau_u, tau_v, shf, lhf
REAL :: cdu, cdv, ch, ce
REAL :: deltaT,deltaq, adj_ref_ht, spdu, spdv, spd
REAL :: step, step2
REAL :: fast_trilin
REAL :: rho, Cp, Lv, Rd, tvirt
REAL :: temp_u, temp_v
temp_u = 0.0
temp_v = 0.0
step = 1.0
step2 = 2.0
! fast_trilin: trilinear interpolation function
! cd_uin : values of drag coefficient sent to fast_trilin
! cd_vin : values of drag coefficient sent to fast_trilin
! ch_in : values of heat transfer coefficient sent to fast_trilin
! ce_in : values of moisture transfer coefficient sent to fast_trilin
! x_in: two values of wind speed sent to fast_trilin
! x_uin: two values of 'u' component wind speed sent to fast_trilin
! x_vin: two values of 'v' component wind speed sent to fast_trilin
! y_in: two values of sea/air temp difference sent to fast_trilin
! z_in: two values or adjusted referecne height sent to fast_trilin
! u: 'u' component of wind speed (input variable) [m/s]
! v: 'v' component of wind speed (input variable) [m/s]
! at: air temperature at reference height (input variable) [Celsius]
! sst: sea surface temperature (input variable) [Celsius]
! deltaT: sea surface temp - air temp [Celsius]
! ref_ht: reference height (input variable) [m]
! Hsig: significant wave height (input variable) [m]
! u_orb: 'u' component of orbital velocity (input variable) [m/s]
! v_orb: 'v' component of orbital velocity (input variable) [m/s]
! qs: specific humidity at surface (input variable) [kg/kg]
! q: specific humidity at reference height (input variable) [kg/kg]
! deltaq: qs-q [kg/kg]
! press: surface pressure (input variable) [Pa]
! adj_ref_ht: adjusted reference height [m]
! spdu: adjusted 'u' component wind speed [m/s]
! spdv: adjusted 'v' component wind speed [m/s]
! spd: adjusted wind speed [m/s]
! cdu: 'u' component drag coefficient returned from fast_trilin
! cdv: 'v' component drag coefficient returned from fast_trilin
! ch: heat transfer coefficient returned from fast_trilin
! ce: moisture transfer coefficient returned from fast_trilin
! tau_u: 'u' component stress (output) [N/m^2]
! tau_v: 'v' component stress (output) [N/m^2]
! shf: sensible heat flux (output) [W/m^2]
! lhf: latent heat flux (output) [W/m^2]
! rho: density of air [kg/m^3]
! Cp: heat capcity of air [J/kg K]
! Lv: latent heat of vaporization [J/kg]
! Rd: gas constant [J/K kg]
! tvirt: virtual temperature [Kelvin]
! step: increment between windspeed and tdif values
! step2: increment between reference height values
adj_ref_ht = ref_ht - (0.8*Hsig)
deltaq = qs - q
deltaT = sst - at
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)
i_u = INT(spdu/step)
i_v = INT(spdv/step)
k = INT(adj_ref_ht/step2)
j = INT(deltaT/step)
IF ((i .GE. 0) .AND. (i .LT. 60)) 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 = 61
spd = spd_table(i_index)
END IF
IF ((i_u .GE. 0) .AND. (i_u .LT. 60)) 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. 60)) THEN
! Accounts for situations when the orbital velocity is greater
! than the u-component of the wind speed. This results in a
! negative stress or a flux of momentum from the ocean to the
! 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 = 61
spdu = spd_table(u_index)
END IF
IF ((i_v .GE. 0) .AND. (i_v .LT. 60)) 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. 60)) THEN
! Accounts for situations when the orbital velocity is greater
! than the v-component of the wind speed. This results in a
! negative stress or a flux of momentum from the ocean to the
! 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 = 61
spdv = spd_table(v_index)
END IF
IF ((k .GE. 1) .AND. (k .LT. 20)) 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. 20) THEN
k_index = 20
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. 60) .AND. (ABS(i_u) .LT. 60) &
.AND. (ABS(i_v) .LT. 60) .AND. (k .GE. 1) &
.AND. (k .LT. 20) .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
! 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
! WARNING: Wind speed (u-component, v-component, or overall) and/or
! adjusted reference height and/or sea-air temperature
! 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
! 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
! ------------------------ FUNCTION fast_trilin ------------------
REAL FUNCTION fast_trilin(x,y,z,val,xp,yp,zp)
! Trilinearly interpolates from one regular grid to another regular
! grid. The regular grids allow the neighbooring points to be
! known, thus elimenating the search for these points.
! Given a location xp,yp,zp and an array val with 2 locations in
! the vector x, 2 locations in y, and 2 locations in the vector z
! return val(xp,yp,zp) using trilinear interpolation
IMPLICIT NONE
REAL, INTENT(IN), DIMENSION(2) :: x,y,z
REAL, INTENT(IN), DIMENSION(2,2,2) :: val
REAL, INTENT(IN) :: xp, yp, zp
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 FUNCTION fast_trilin