; $Id: HeatFluxarray.pro,v 5.4 ; ; Copyright (c) 2001 Javier G Corripio. ; jgc@geo.ed.ac.uk ;+ ; NAME: ; HEATFLUXARRAY ; ; PURPOSE: ; Computes heat fluxes (sensible and turbulent) to the atmosphere ; computes a 2D array at a time ; ; CATEGORY: ; Atmospheric Sciences ; ; CALLING SEQUENCE: ; Result = HeatFluxarray(Ta(K),Ts(K),RH(%),u_bar,Z_u,Z_T,Z_RH,z,z_0) ; ; 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 heatfluxarray, Ta,Ts,RH,u_bar,Z_u,Z_T,Z_RH,z,z_0 ON_ERROR, 2 IF (n_params() EQ 0) THEN BEGIN message, 'USAGE: result = heatfluxarray(Ta(K),Ts(K),RH(%),u_bar'+$ ',Z_u,Z_T,Z_RH,z,z_0])',/INFO message, 'RETURNS: [H,Q_EV,u_star,L]',/INFO RETURN,!values.F_NaN ENDIF ; Ta 2D array ; Ts, RH, u_bar, z_0 may be 2D array ;-------- 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.4 ;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 --------------------------------------------------- ah = 1. ;ratio of eddy diffusivity and viscosity for heat ae = 1. ;ratio of eddy diffusivity and viscosity for water vapor 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 ;-------- Initial values for iterations sz=size(Ta,/dimensions) cols = sz[0] rows = sz[1] u_star = FLTARR(cols,rows) + 1. u_star0 = FLTARR(cols,rows) ; friction velocity_0 EV = FLTARR(cols,rows) - 20. EV0 = FLTARR(cols,rows) ; Evaporation_0 H = FLTARR(cols,rows) + 30. H0 = FLTARR(cols,rows) ; sensible heat_0 L = FLTARR(cols,rows) + 1e37 ; Obukhov stability length initial v. large value L0 = FLTARR(cols,rows) Psism = FLTARR(cols,rows) ; stability function for momentum Psish = FLTARR(cols,rows) ; stability function for heat Psisv = FLTARR(cols,rows) ; stability function for water vapor Zetasm = FLTARR(cols,rows) Zetash = FLTARR(cols,rows) Zetasv = FLTARR(cols,rows) chi = FLTARR(cols,rows) z_h = FLTARR(cols,rows) z_q = FLTARR(cols,rows) maxchg = 1000. numcells = cols*rows ; ------------------- 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) q_s = rh2sh(100,Ts,Pz,/ice) ;specific humidity of air over snow saturated 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 ----------------------------- ;----------------------------------------------------------------------------- ;----------- Iterations for turbulent fluxes -------------------------------- itercounter=0 WHILE (maxchg GT 1.5) DO BEGIN maxchg = MAX(ABS ([H-H0,EV-EV0]) ) itercounter=itercounter+1 change = WHERE(ABS(L-L0) GT 0.01 OR ABS(u_star-u_star0) GT 0.01 OR $ ABS(H-H0) GT 0.01 OR ABS(EV-EV0) GT 0.01 , countnchg) ;print,'iteration No', itercounter ;print,'time ',systime(/utc) ;print, 'maximum change',maxchg,' change % =',(size(change,/dim)/float(numcells))*100 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) ;---------- 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 ;--- A) ln(z_s/z_0)= b0+b1*ln(Re)+b2*(ln(Re)^2) ReLE0_135 = WHERE (Re LE 0.135, count1) ReLE2_5 = WHERE (Re GT 0.135 AND Re LE 2.5, count2) ReGT2_5 = WHERE (Re GT 2.5, count3) b0_andreasT = 1.250 & b0_andreasV = 1.610 IF (count1 NE 0) THEN BEGIN z_h[ReLE0_135] = z_0[ReLE0_135]*EXP( b0_AndreasT + $ b1_AndreasT*ALOG(Re[ReLE0_135]) + $ b2_AndreasT*(ALOG(Re[ReLE0_135])^2) ) z_q[ReLE0_135] = z_0[ReLE0_135]*EXP( b0_AndreasV + $ b1_AndreasV*ALOG(Re[ReLE0_135]) + $ b2_AndreasV*(ALOG(Re[ReLE0_135])^2) ) ENDIF b0_andreasT = 0.149 & b0_andreasV = 0.351 b1_andreasT = -0.550 & b1_andreasV = -0.628 IF (count2 NE 0) THEN BEGIN z_h[ReLE2_5] = z_0[ReLE2_5]*EXP( b0_AndreasT + $ b1_AndreasT*ALOG(Re[ReLE2_5]) + $ b2_AndreasT*(ALOG(Re[ReLE2_5])^2) ) z_q[ReLE2_5] = z_0[ReLE2_5]*EXP( b0_AndreasV + $ b1_AndreasV*ALOG(Re[ReLE2_5]) + $ b2_AndreasV*(ALOG(Re[ReLE2_5])^2) ) ENDIF b0_andreasT = 0.317 & b0_andreasV = 0.396 b1_andreasT = -0.565 & b1_andreasV = -0.512 b2_andreasT = -0.183 & b2_andreasV = -0.180 IF (count3 NE 0) THEN BEGIN z_h[ReGT2_5] = z_0[ReGT2_5]*EXP( b0_AndreasT + $ b1_AndreasT*ALOG(Re[ReGT2_5]) + $ b2_AndreasT*(ALOG(Re[ReGT2_5])^2) ) z_q[ReGT2_5] = z_0[ReGT2_5]*EXP( b0_AndreasV + $ b1_AndreasV*ALOG(Re[ReGT2_5]) + $ b2_AndreasV*(ALOG(Re[ReGT2_5])^2) ) ENDIF ; 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 ;!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ;------ heat, Brutsaert:4.35' IF countnchg NE 0 THEN $ H[change] = ((theta_a[change] - theta_s[change])*ah*karm*u_star[change]*rho*C_p)/$ (ALOG(((Z_T-d_0[change])/z_h[change]))-Psish[change]) ;------ water vapour, Brutsaert:4.33' IF countnchg NE 0 THEN $ EV[change] = ((q_a[change] - q_s[change])*ae*karm*u_star[change]*rho)/$ (ALOG( ((z_RH-d_0[change])/z_q[change]) )-Psisv[change]) ;------ Obukhov Length, Brutsaert:4.25 IF countnchg NE 0 THEN $ L[change] = (-(u_star[change]^3)*rho)/( karm*earth_g*((H[change]/$ (Tmean[change]*C_pm))+0.61*EV[change]) ) ;----- 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) zetalt0 = WHERE(Zeta LE 0, countz1) zetagt0 = WHERE(Zeta GT 0, countz2) IF (countz1 NE 0) THEN BEGIN ; Unstable conditions ------------ chi[zetalt0] = (1-16*Zeta[zetalt0])^0.25 Psisv[zetalt0] = 2*ALOG((1+chi[zetalt0]^2)/2.0) ; Brutsaert 4.49 Psism[zetalt0] = 2*ALOG((1+chi[zetalt0])/2.0)+ $ ALOG((1+chi[zetalt0]^2)/2.0)-2*ATAN(chi[zetalt0])+!PI/2.0 Psish[zetalt0] = 2*ALOG((1+chi[zetalt0]^2)/2.0) ; Brutsaert 4.51 ENDIF IF (countz2 NE 0) THEN BEGIN ; Stable conditions ------------ beta = 5. ; Marks and Dozier (1992) Psisv[zetagt0] = -beta*Zeta[zetagt0] ; Brutsaert 4.55 Psism[zetagt0] = -beta*Zeta[zetagt0] ; Brutsaert 4.56 Psish[zetagt0] = -beta*Zeta[zetagt0] ; Brutsaert 4.57 ENDIF ;------------------------------------------------------------------- ;----------- in case iterations do not converge -------------------- if (itercounter GT 50) THEN begin H[*,*] = 0.0 EV[*,*] = 0.0 u_star[*,*] = 0.0 L[*,*] = 0.0 maxchg = 0.0 endif ENDWHILE ;-------------------- END ITERATIONS ------------------------------------- Q_EV = EV*L_s ; latent heat in Watts funcs = [z_h,z_q,Re] ;heatdat = [H,Q_EV,u_star,L] heatdat = fltarr(2,cols,rows) heatdat[0,*,*]=H heatdat[1,*,*]=Q_EV RETURN, heatdat END