; $Id: HeatFluxptt.pro,v 5.4 ; ; Copyright (c) 2001 Javier G Corripio. ; jgc@geo.ed.ac.uk ;+ ; NAME: ; HEATFLUX ; ; PURPOSE: ; Computes heat fluxes (sensible and turbulent) to the atmosphere ; ; CATEGORY: ; Atmospheric Sciences ; ; CALLING SEQUENCE: ; Result = HeatFlux() ; ; INPUTS: ; Ta air temperature (K) ; Ts snow temperature (K) ; RH relative humidity (%) ; u_bar wind speed (m/s) ; Z_u windspeed measurement height (m) ; Z_T Temperature measurement height (m) ; Z_RH RH measurement height (m) ; z altitude a.s.l. (m) ; z_0 roughness length (m) ; ; ; KEYWORDS: ; ZZZ makes all roughness lengths equal ; z_0 = z_q = z_T ; RiNum uses abulk Richardson Number instead of Monin Obukhov Length ; ; ; OUTPUTS: ; ; ; SUBROUTINES: ; ; REQUIREMENTS: ; FUNCTIONS wvapsat, z2p, rh2sh,rh2mxrt,potential_t ; ; NOTES: ; ; ; COMMENTS: ; ; ; REFERENCES: ; Andreas, E. L.: 1987, A theory for the scalar roughness and the scalar ; transfer coefficients over snow and sea ice, ; Boundary layer Meteorology 38, 159-184. ; ; Brutsaert, W.: 1982, Evaporation into the atmosphere : theory, history, ; and applications, Reidel, Dordrecht. ; ; Jacobson, M. Z.: 1999, Fundamentals of Atmospheric Modeling, ; Cambridge University Press, Cambridge. ; ; Marks, D. and Dozier, J.: 1992b, Climate and energy exchange at the ; snow surface in the alpine region of the Sierra Nevada. 2. Snow cover ; energy balance, Water Resources Research 28(11), 3043 3054. ; ; Martin, E. and Lejeune, Y.: 1998, Turbulent fluxes above the snow ; surface, Annals of Glaciology 26, 179 -183. ; ; ; EXAMPLE: ; Q_HL=heatflux(280.0,273.15,65.,5.,2.,1.8,1.8,2900.,.0026) ; ; ; MODIFICATION HISTORY: ; JGC, 13 November, 2001 ; jgc July, 2002 z_0 /= z_q /= z_t (Andreas 1987) ; jgc August,2002 modular call to atmospheric variables ; ;- @wvapsat @z2p @rh2sh @rh2mxrt @p2rho @potential_t FUNCTION heatflux, Ta,Ts,RH,u_bar,Z_u,Z_T,Z_RH,z,z_0,FUNCS=funcs,RiNum=rinum,$ VERBOSE=verbose,ZZZ=zzz, WATER=water ON_ERROR, 2 IF (n_params() EQ 0) THEN BEGIN message, 'USAGE: result = heatflux(Ta(K),Ts(K),RH(%),u_bar'+$ ',Z_u,Z_T,Z_RH,z,z_0[FUNCS=funcs,/RiNum])',/INFO message, 'RETURNS: [H,Q_EV,u_star,L]',/INFO message, 'RETURNS: [H,Q_EV,Ri,Cn] if /RiNum set',/INFO RETURN,!values.F_NaN ENDIF ;-------- NOMENCLATURE ----------------------------------------------- ; theta potential temperature. Brutsaert 3.35 (Jacobson 2.86) ; q_a rho_v/rho specific humidity air. Brutsaert 3.2 ; e_v partial pressure water vapour ; e_v_star vapour pressure at saturation ; rho_d density dry air ; rho_v density water vapour ; L_s latent heat of sublimation L_e + L_m ; L_e latent heat of evaporation f(T) ; L_m latent heat of melting f(T) ; Psism ; Psish ; Psisv ; rho_d density dry air ; rho_v density water vapor ; rho density moist air ; q_s specific humidity air over snow ; q_a specific humidity air at screen level ; L Obukhov stability length ; u_star friction velocity ; Z_u windspeed measurement height (m) ; Z_RH RH measurement height (m) ; Z_T Temperature measurement height (m) ; u_bar average wind speed (m/s) ; RH relative humidity (%) ; Ta air temperature (K) ; Ts snow temperature (K) ; z altitude a.s.l. (m) ; z_0 roughness length (m) ; d_0 zero-plane displacement height (m) ; Ri Richardson bulk number ; Cn ; Ch ; p_v_a ; p_v_s ; mu coefficient of dynamic viscosity for air ; nu kinamatic viscosity of air ; b0_Andreas1,b1_,b2_ cofficients relating z_0, z_h and z_q ; Re roughness Reynolds number = u_star z_0 nu^-1 ;-------- CONSTANTS ------------------------------------------------ ;** most values from Jacobson99 R_star = 8.3145 ;Universal gas constant J/molK R_d = 287.04 ;Gas constant for dry air J/KgK R_v = 461.40 ;Gas constant for water vapor J/KgK R_Rdv = R_d/R_v ;Ratio Rd/Rv = 0.622 karm = 0.41 ;von Karman's constant C_p = 1004.67 ;Specific heat of dry air @ ct. P J/KgK ah = 1.0 ;ratio of eddy diffusivity and viscosity for heat ae = 1.0 ;ratio of eddy diffusivity and viscosity for water vapor sigma = 5.67051E-08 ;Stephan-Boltzman constant W/m^2K^4 P_0 = 1000.0 ;standard pressure level (hPa) Mv = 18.016 ;molecular weight water vapor Md = 28.966d0 ;Molecular weight of dry air REarth = 6.3756766D6 ;Average earth's radius (m) atmos_P0 = 1.013250D5 ;Standard sea-level atmospheric pressure Pa atmos_T0 = 0288.15d0 ;Standard sea-level Temperature atmos_R = 8.31432D0 ;Universal gas constant (J deg-1 kmol-1) earth_G = 9.80665D0 ;Acceleration due to gravity (m s-2) stlapse = -0.0065d0 ;standard lapse rate K/m mu = 1.78e-5 ;kg m^(-1)s^(-1) coefficient of dynamic viscosity for air ;-------- VARIABLES --------------------------------------------------- d_0 = (2/3.0)*7.35*z_0 ; zero displacement height Brutsaert, p.116 Tsc = Ts-273.15 L_s = 2.83458 * 1e6 - Tsc*(340+10.46*Tsc) ; latent heat of sublimation IF Keyword_set(water) THEN BEGIN L_e = 2.501 * 1e6 + Tsc * (2030 - 10.46*Tsc) ; latent heat of evaporation L_s = L_e ENDIF ;-------- Initial values for iterations u_star = 1. u_star0 = 0. ;friction velocity_0 EV = 1. EV0 = 0. ;Evaporation_0 H = 1. H0 = 0. ;sensible heat_0 L = 1e37 ;Obukhov stability length initial v. large value L0 = 0. Psism = 0. ;stability function for momentum Psish = 0. ;stability function for heat Psisv = 0. ;stability function for water vapor Zetasm = 0. Zetash = 0. Zetasv = 0. chi = 1. maxchg = 1. ; ------------------- CALCULATIONS ----------------------------------------- Pz = z2p(z) ;pressure at current altitude rho = P2rho(Pz,Ta,RH) ;air density q_a = rh2sh(RH,Ta,Pz) ;specific humidity of air (g/kg) IF Keyword_set(water) THEN BEGIN q_s = rh2sh(100.,Ts,Pz) ;specific humidity of air over saturated water ENDIF ELSE BEGIN q_s = rh2sh(100.,Ts,Pz,/ice) ;specific humidity of air over saturated snow/ice ENDELSE q_a = q_a/1000. ;specific humidity in kg/kg q_s = q_s/1000. ;specific humidity in kg/kg theta_a = potential_t(Ta,Pz) ; potential temperature theta_s = potential_t(Ts,Pz) ; potential temperatur air over snow nu = mu/rho ;kinematic viscosity of air mxratio = rh2mxrt(RH,Ta,Pz) ;mass mixing ratio of water vapour mxratio = mxratio/1000. ;mass mixing in kg/kg C_pm = C_p*(1+0.859*mxratio) ;Specific heat of moist air @ ct. P J/KgK ; Jacobson (1999) Tmean = (Ta + Ts)/2. ;----------------------------------------------------------------------------- ;---------------------- MONIN - OBUKHOV ----------------------------- ;----------------------------------------------------------------------------- IF keyword_set(rinum) THEN GOTO, RiNumonly ;----------- Iterations for turbulent fluxes -------------------------------- itercounter=0 WHILE (maxchg GT 0.001) DO BEGIN maxchg = MAX(ABS ([L-L0, u_star-u_star0, H-H0,EV-EV0]) ) itercounter=itercounter+1 L0 = L u_star0 = u_star H0 = H EV0 = EV ;----- friction velocity, Brutsaert:4.34' u_star = (u_bar*karm)/(ALOG((Z_u-d_0)/z_0)-Psism) IF (U_STAR LT 0) THEN BEGIN PRINT, 'u_star,u_bar,Psism' PRINT, u_star,u_bar,Psism ENDIF ;---------- scalar roughness lengths from Andreas (1987) pag. 177 ;--- for Re<1000 See also van der Broeke (96), p. 53 Re = u_star*z_0/nu ;; roughness Reynolds number ;;Andreas polynomials for ratio ln(z_s/z_0) for different values of Re p.177 b0_andreasT = 0.0 & b0_andreasV = 0.0 b1_andreasT = 0.0 & b1_andreasV = 0.0 b2_andreasT = 0.0 & b2_andreasV = 0.0 IF (Re LE 0.135) THEN BEGIN b0_andreasT = 1.250 & b0_andreasV = 1.610 ENDIF IF (Re GT 0.135 AND Re LE 2.5) THEN BEGIN b0_andreasT = 0.149 & b0_andreasV = 0.351 b1_andreasT = -0.550 & b1_andreasV = -0.628 ENDIF IF (Re GE 2.5 ) THEN BEGIN b0_andreasT = 0.317 & b0_andreasV = 0.396 b1_andreasT = -0.565 & b1_andreasV = -0.512 b2_andreasT = -0.183 & b2_andreasV = -0.180 ENDIF ;IF (Re GT 1000.0) THEN Re= 1000.0 ; this is fiddling: Andreas max Re=1000. ;--- A) ln(z_s/z_0)= b0+b1*ln(Re)+b2*(ln(Re)^2) z_h = z_0*EXP( b0_AndreasT + b1_AndreasT*ALOG(Re) + b2_AndreasT*(ALOG(Re)^2) ) z_q = z_0*EXP( b0_AndreasV + b1_AndreasV*ALOG(Re) + b2_AndreasV*(ALOG(Re)^2) ) ; B) bluff-rough surfaces ---------------------------------- ;z_hb = z_0*7.4*EXP(-2.46*Re^.25) ;Brutsaert: 5.29 ;z_qb = z_0*7.4*EXP(-2.25*Re^.25) ;Brutsaert: 5.28 ; C) rough flow over water -------------------------------- ;z_hc = 0.169*exp(-1.53*u_star^.25)/100. ;Brutsaert: 5.34 ;z_qc = 0.169*exp(-1.40*u_star^.25)/100. ;Brutsaert: 5.33 ; Average B), C) for very wet snow ;z_h = (z_hb + z_hc)/2. ;z_q = (z_qb + z_qc)/2. ;z_h = z_hb ;z_q = z_qb IF KEYWORD_SET(zzz) then BEGIN z_h = z_0 z_q = z_0 ENDIF ;Formulation of Marks and Dozier H = ((theta_a - theta_s)*ah*karm*u_star*rho*C_p)/(ALOG((Z_T-d_0)/z_h) -Psish) EV = ((q_a - q_s)*ae*karm*u_star*rho)/(ALOG( (z_RH-d_0)/z_q )-Psisv) L = ((u_star^3.)*rho)/( karm*earth_g*((H/(Ta*C_p))+0.61*EV) ) ;----- stability functions Brutsaert p. 62-71 ; phisv = phish = phism^2 Brutsaert 4.45 ; zeta <= 0: unstable, zeta > 0: stable Zeta = (Z_U-d_0)/L ; Brutsaert 4.24 (Obukhov 1954) IF(Zeta LE 0) THEN BEGIN ; Unstable conditions ------------ chi = (1-16*Zeta)^0.25 Psisv = 2*ALOG((1+chi^2.)/2.0) ; Brutsaert 4.49 Psism = 2*ALOG((1+chi)/2.0)+ALOG((1+chi^2)/2.0)-2*ATAN(chi)+!PI/2.0 Psish = 2*ALOG((1+chi^2.)/2.0) ; Brutsaert 4.51 ENDIF ELSE BEGIN ; Stable conditions ------------ Psism = -5.2*Zeta ;Brutsaert 4.55 IF zeta GT 1 THEN Psism = -5.2 Psish = Psism Psisv = Psism ENDELSE ;------------------------------------------------------------------- ;----------- in case iterations do not converge -------------------- if (itercounter GT 1000) THEN begin H = 0.0 EV = 0.0 u_star = 0.0 L = 0.0 maxchg = 0.0 IF KEYWORD_SET(verbose) THEN BEGIN text='Ta= '+string(Ta)+'Ts= '+string(Ts)+string(RH)+string(u_bar)+$ ' iteration did not converge after 1000 cicles' message,text,/INFO ENDIF endif ENDWHILE ;print, 'zeta', zeta ;print, 'theta_a - theta_s',(theta_a - theta_s) ;print, 'ALOG((Z_T-d_0)/z_h)',ALOG((Z_T-d_0)/z_h) ;print, 'psism',psism ;print, 'Zeta', zeta ;print, 'L',L ;print, 'u_star', u_star ;-------------------- END ITERATIONS ------------------------------------- Q_EV = EV*L_s ; latent heat in Watts funcs = [z_h,z_q,Zeta,Re] heatdat = [H,Q_EV,u_star,L] RETURN, heatdat ;----------------------------------------------------------------------------- ;----------------------------- Ri Number ----------------------------------- ;----------------------------------------------------------------------------- RiNumonly: thetav_a = theta_a*(1+0.61*q_a) ; Stull 1.5.1b & Brutsaert 3.37a thetav_s = theta_s*(1+0.61*q_s) ; ; q in g/g (or kg/kg) e_star= wVapSat(Ta) p_v_a = e_star*RH/100. p_v_s = wVapSat(Ts,/ICE) Ri=(earth_G*(thetav_a-thetav_s)*Z_T)/(thetav_a*(u_bar^2)) ;** Stull, 5.6.3 IF (Ri GT 0.21) THEN BEGIN ;laminar flow, no turbulence Stull 1988 5.6.2 H = 0.0 Q_EV = 0.0 Cn = 0.0 GOTO, noturbulent ENDIF ;*** coefficients after Martin & Lejeune (1998) Cn = (karm^2.)/( ALOG(Z_U/z_0)*ALOG(Z_T/z_0) ) Ch = Cn IF (Ri GT 0) THEN Ch = Cn*(1-5*Ri)^2 IF (Ri LT 0) THEN BEGIN alpha_par = 0.83*Cn^(-0.62) Ch = Cn*( 1+7*ALOG(1-alpha_par*Ri)/alpha_par ) ;** p_v_a already incorporates RH ENDIF Ce = Ch ; assumes Ch=Ce; check Andreas H = rho*C_p*Ch*u_bar*(Ta-Ts) Q_EV = (L_s*rho/Pz)*(Mv/Md)*Ce*u_bar*(p_v_a - p_v_s) ;** p_v_a already incorporates RH noturbulent: ;jumps here if flow is laminar, not turbulent heatdat = [H,Q_EV,Ri,Cn] RETURN, heatdat END