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