#ifdef SPHUM subroutine turbhf(tair, shum, sst, wspd, press, qse, qla, icec, tfrezc) #else subroutine turbhf(tair, tdew, sst, wspd, press, qse, qla, icec, tfrezc) #endif ! ! Introduction: ! ============= ! ! FORTRAN90 Code to calculate turbulent, i.e. sensible and latent heat fluxes. ! The code has kindly been provided by: ! ! A. Birol Kara ---> U.S. Naval Research Laboratory (NRL) ! Alan J. Wallcraft ---> U.S. Naval Research Laboratory (NRL) ! Harley E. Hurlburt ---> U.S. Naval Research Laboratory (NRL) ! ! The algorithm of Kara et al. (in press) yields air-sea fluxes ! that are highly accurate relative to those obtained with the ! standard algorithm used by the Tropical Ocean Global Atmosphere ! Coupled Ocean-Atmosphere Response Experiment (TOGA COARE). ! ! The transfer coefficients were calculated using latest algorithm ! by running TOGA COARE 3.0 (Fairall et al., 2003) for various ! temperature differences and winds. ! The humidity effect on the transfer coefficients were also taken ! into account via the COARE 3.0. Because stability depends on ! temperature and humidity differences, as well, the algorithm ! includes full stability. ! ! The code was expanded by using sea ice concentration (Roske, in press). ! By using the preprocessor keyword SPHUM the option to use 2 m specific ! humidity can be activated. Default is to use 2 m dew point temperature. ! ! The Kara-Algorithm needs input parameters at standard height (10m). ! It was applied to the second version of the "OMIP-Forcing" (365 days) ! which is based on the ECMWF Re-Analyses (ERA). ! (http://www.omip.zmaw.de/omip.php) ! ! ERA do not provide air temperature at 10m but at 2m. ! An attempt to derive T(10m) from T(2m) by using the original TOGA COARE ! algorithm was unfortunately not successful. The problem was not pursued. ! 2m air temperature was used instead. ! ! ! Further details or earlier version of algorithms can be obtained ! directly from Frank Roske (roeske@dkrz.de). ! ! Frank Roske, Max-Planck-Institute for Meteorology, April 2005. ! ! ! ! References: ! =========== ! ! Buck, A.L. (1981) New Equations for Computing Vapor Pressure and ! Enhancement Factor. Journal of Applied Meteorology, Vol.20, pp.1527-1532. ! ! Fairall, C. W., Bradley, E. F., Hare, J. E., Grachev, A. A., ! Edson, J. B., 2003: Bulk Parameterization of Air-Sea Fluxes: ! Updates and Verification for the COARE Algorithm. ! J. Climate 16, 571-591. ! ! Kara, A. B., Hurlburt, H. E., Wallcraft, A. J., in press. ! Stability-Dependent Exchange Coefficients of Air-Sea Fluxes. ! Journal of Atmospheric Oceanic Technology. ! ! Roske, F., in press. A global heat and freshwater forcing ! dataset for ocean models. Ocean Modelling. ! ! ! Input parameters: ! ================= ! ! tair: 2m Air Temperature [Kelvin] #ifdef SPHUM ! shum: 2m Specific Humidity [kg/kg] #else ! tdew: 2m Dewpoint Temperature [Kelvin] #endif ! sst: Sea Surface Temperature [Kelvin] ! wspd: Scalar Wind Speed [m/s] ! press: Mean Sea Level Pressure [Pa} ! icec: Ice concentration [0-1] ! tfrezc: Freezing point temperature [degree Celsius] ! ! ! ! Output parameters: ! ================== ! ! qla: Latent Heat Flux [W/m**2] ! qse: Sensible Heat Flux [W/m**2] ! ! !********** ! implicit none ! physical constants. those having no dependence are assigned values. real , parameter :: & cpa = 1004.67, & ! specific heat of dry air (j/kg/K, businger 1982) epsilon = 0.62197, & ! parameter for saturated specific humidity rgas = 287.1, & ! gas const. dry air (j/kg/K) tzerok = 273.16, & ! tok - 0 deg. C in units of Kelvin xlatent_heat_sub = 2.834e6 real :: tfrezc real , parameter :: rmiss=-999.0,eps=1.e-7 real :: cl, cs, qla, qse, xlatent_heat_vap, rhoa, & tdiff, tdiffw, clw, csw, qdiff, qdiffw #ifdef SPHUM real :: tair , sst , shum , press , wspd, icec real :: tairc, sstc, pressmb, qa, qsat, esat, esatw, qsatw #else real :: tair , sst , tdew , press , wspd, icec real :: tairc, sstc, tdewc, pressmb, qa, qsat, ee, esat, esatw, qsatw #endif sstc = sst - tzerok ! [deg. Celsius] tairc = tair - tzerok ! [deg. Celsius] #ifndef SPHUM tdewc = tdew - tzerok ! [deg. Celsius] #endif pressmb = press/100. ! [mb or hPa] !--------------------------------------------------------------------- ! Vapor pressure formulae according to Buck (1981): ! Ref. see above. !--------------------------------------------------------------------- ! #ifdef SPHUM qa = shum #else ee = 6.1121*exp((18.729-tdewc/227.3)*tdewc/(tdewc+257.87)) ! vapor pressure in mb qa = epsilon*(ee / (pressmb - 0.378*ee) ) ! spec. humidity kg/kg #endif if(sstc.le.tfrezc)then esat = 6.1115*exp((23.036-min(sstc ,0.)/333.7)*min(sstc ,0.)/(max(sstc ,-80.)+279.82)) esatw = 6.1115*exp((23.036-min(tfrezc ,0.)/333.7)*min(tfrezc ,0.)/(max(tfrezc ,-80.)+279.82)) else esat = 0.9815*6.1121*exp((18.729-sstc/227.3)*sstc/(sstc+257.87)) endif qsat = epsilon*( esat /(pressmb - 0.378*esat ) ) ! spec. humidity kg/kg qsatw = epsilon*( esatw /(pressmb - 0.378*esatw ) ) ! spec. humidity kg/kg tdiff = tairc - sstc tdiffw = tairc - tfrezc qdiff = qa - qsat qdiffw = qa - qsatw call kara(tairc, sstc, wspd, qa, cl, pressmb) cs = 0.9554*cl call kara(tairc, tfrezc, wspd, qa, clw, pressmb) csw = 0.9554*clw rhoa = press/(rgas*tair*(1.0 + 0.61*qa) ) ! Sensible heat flux (upward fluxes) if(sstc.le.tfrezc) then if(abs(icec-rmiss).gt.eps) then qse = (cpa*rhoa*wspd)*(cs*tdiff*icec + csw*tdiffw*(1.-icec)) else qse = (cpa*rhoa*wspd)*(cs*tdiff) endif else qse = cpa*rhoa*wspd*cs*tdiff endif ! Latent heat flux if(sstc.le.tfrezc) then if(abs(icec-rmiss).gt.eps) then xlatent_heat_vap = (2.501 - 0.00237*tfrezc)*1e+6 qla = rhoa*wspd*(xlatent_heat_sub*cl*qdiff*icec + xlatent_heat_vap*clw*qdiffw*(1.-icec)) else qla = xlatent_heat_sub*rhoa*wspd*(cl*qdiff) endif else xlatent_heat_vap = (2.501 - 0.00237*sstc)*1e+6 qla = xlatent_heat_vap*rhoa*wspd*cl*qdiff endif end subroutine turbhf subroutine kara(tairc, sstc, wspd, qa, cl, pressmb ) real cl,qq,va,qva,tamts,tairc,sstc,wspd,qa,ee,pressmb real, parameter :: tzerok = 273.16 ! ! --- Ta-Ts ! --- ============== ! --- STABLE: 7 to 0.75 degC ! --- NEUTRAL: 0.75 to -0.75 degC ! --- UNSTABLE: -0.75 to -8 degC ! ! --- Va ! --- ============== ! --- Low: 1 to 5 m/s ! --- High: 5 to 40 m/s ! real, parameter :: vamin= 1.2, vamax=40.0 real, parameter :: tdmin=-8.0, tdmax= 7.0 real, parameter :: & as0_00=-0.2925, as0_10= 0.07272, as0_20=-0.006948, & as0_01= 0.5498, as0_11=-0.1740, as0_21= 0.01637, & as0_02=-0.05544, as0_12= 0.02489, as0_22=-0.002618 real, parameter :: & as5_00= 1.023, as5_10=-0.002672, as5_20= 0.001546, & as5_01= 0.009657, as5_11= 0.2103, as5_21=-0.06228, & as5_02=-2.281e-5, as5_12=-5.329, as5_22= 0.5094 real, parameter :: & au0_00= 2.077, au0_10=-0.2899, au0_20=-0.01954, & au0_01=-0.3933, au0_11= 0.07350, au0_21= 0.005483, & au0_02= 0.03971, au0_12=-0.006267, au0_22=-0.0004867 real, parameter :: & au5_00= 1.074, au5_10= 0.006912, au5_20= 0.0001849, & au5_01= 0.005579, au5_11=-0.2244, au5_21=-0.002167, & au5_02= 5.263e-5, au5_12=-1.027, au5_22=-0.1010 real, parameter :: & an0_00= 1.14086, an5_00= 1.073, & an0_01=-0.003120, an5_01= 0.005531, & an0_02=-0.0009300, an5_02= 5.433e-05 real, parameter :: & ap0_00= as0_00 + as0_10*0.75 + as0_20*0.75**2, & ap0_01= as0_01 + as0_11*0.75 + as0_21*0.75**2, & ap0_02= as0_02 + as0_12*0.75 + as0_22*0.75**2 real, parameter :: & ap5_00= as5_00 + as5_10*0.75 + as5_20*0.75**2, & ap5_01= as5_01, & ap5_02= as5_02, & ap5_11= as5_11*0.75 + as5_21*0.75**2, & ap5_12= as5_12*0.75 + as5_22*0.75**2 real, parameter :: & am0_00= au0_00 - au0_10*0.75 + au0_20*0.75**2, & am0_01= au0_01 - au0_11*0.75 + au0_21*0.75**2, & am0_02= au0_02 - au0_12*0.75 + au0_22*0.75**2 real, parameter :: & am5_00= au5_00 - au5_10*0.75 + au5_20*0.75**2, & am5_01= au5_01, & am5_02= au5_02, & am5_11= - au5_11*0.75 + au5_21*0.75**2, & am5_12= - au5_12*0.75 + au5_22*0.75**2 real :: qsatur real, parameter :: epsilon = 0.62197 ! parameter for saturated specific humidity ee = 6.1121*exp((18.729-tairc/227.3)*tairc/(tairc+257.87)) ! vapor pressure in mb qsatur = epsilon*(ee / (pressmb - 0.378*ee) ) ! spec. humidity kg/kg tamts = tairc - sstc - 0.61*(tairc+tzerok) * (qsatur - qa) tamts = min( tdmax, max( tdmin, tamts) ) va = min( vamax, max( vamin, wspd ) ) if (va.le.5.0) then if (tamts.ge. 0.75) then !stable cl = & (as0_00 + as0_01* va + as0_02* va**2) & + (as0_10 + as0_11* va + as0_12* va**2)*tamts & + (as0_20 + as0_21* va + as0_22* va**2)*tamts**2 elseif (tamts.le.-0.75) then !unstable cl = & (au0_00 + au0_01* va + au0_02* va**2) & + (au0_10 + au0_11* va + au0_12* va**2)*tamts & + (au0_20 + au0_21* va + au0_22* va**2)*tamts**2 elseif (tamts.ge.-0.098) then qq = (tamts-0.75)/0.848 !linear between 0.75 and -0.098 qq = qq**2 !favor 0.75 cl = (1.0-qq)*(ap0_00 + ap0_01* va + ap0_02* va**2) & + qq *(an0_00 + an0_01* va + an0_02* va**2) else qq = (tamts+0.75)/0.652 !linear between -0.75 and -0.098 qq = qq**2 !favor -0.75 cl = (1.0-qq)*(am0_00 + am0_01* va + am0_02* va**2) & + qq *(an0_00 + an0_01* va + an0_02* va**2) endif !tamts else !va>5 qva = 1.0/va if (tamts.ge. 0.75) then !stable cl = & (as5_00 + as5_01* va + as5_02* va**2) & + (as5_10 + as5_11*qva + as5_12*qva**2)*tamts & + (as5_20 + as5_21*qva + as5_22*qva**2)*tamts**2 elseif (tamts.le.-0.75) then !unstable cl = & (au5_00 + au5_01* va + au5_02* va**2) & + (au5_10 + au5_11*qva + au5_12*qva**2)*tamts & + (au5_20 + au5_21*qva + au5_22*qva**2)*tamts**2 elseif (tamts.ge.-0.098) then qq = (tamts-0.75)/0.848 !linear between 0.75 and -0.098 qq = qq**2 !favor 0.75 cl = (1.0-qq)*(ap5_00 + ap5_01* va + ap5_02* va**2 & + ap5_11*qva + ap5_12*qva**2) & + qq *(an5_00 + an5_01* va + an5_02* va**2) else qq = (tamts+0.75)/0.652 !linear between -0.75 and -0.098 qq = qq**2 !favor -0.75 cl = (1.0-qq)*(am5_00 + am5_01* va + am5_02* va**2 & + am5_11*qva + am5_12*qva**2) & + qq *(an5_00 + an5_01* va + an5_02* va**2) endif !tamts endif !va cl = cl * 1.e-3 end subroutine kara