!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Program by: D BALA SUBRAHAMANYAM, subrahamanyam@gmail.com ! ! LINUX Version ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! PROGRAM LAST MODIFIED: SEPTEMBER 01, 2006 ! ! PROGRAM FOR ESTIMATION OF AIR-SEA INTERACTION PARAMETERS ! ! PROGRAM IS BASED ON REVISED BULK ALGORITHM SUGGESTED BY ! ! SUBRAHAMANYAM AND RADHIKA, 2002 (J.Atmos.Sol.Terr.Phy.)! ! SUBRAHAMANYAM AND RADHIKA, 2003 (Annal.Geophy.) ! ! ! ! INPUT PARAMETERS ===> ! ! Measurement Height (z, in m) ! ! Wind Speed (ws, in m/s) ! ! Air Temperature (temp, in deg.C) ! ! Sea Surface Temperature (sst, in deg. C) ! ! Relative Humidity (rh, in %) ! ! Pressure (pres, in mb) ! ! ! ! OUTPUT PARAMETERS ===> ! ! All the exchange coeffcients are scaled for 10-m ! ! Drag Coefficient (cd x 10^-3) ! ! Exchange Coefficient for heat and moisture (ch and ce) ! ! (ch x 10^-3 and ce x 10^-3) ! ! Drag Coefficient for neutral stability ! ! (cdn x 10^-3) ! ! Exchange Coefficients for heat and moisture ! ! for neutral stability ! ! (chn x 10^-3 and cen x 10^-3) ! ! Sensible Heat Flux (shf, in Wm^-2) ! ! Latent Heat Flux (xlhf, in Wm^-2) ! ! Momentum Flux (xmf, in Nm^-2) ! ! Friction Velocity (ustarf, in ms^-1) ! ! Stability Parameter (zlf) ! ! Roughness Length - wind (z0) ! ! Roughness Length - temp (z0t) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine dbsrrz (z, sst, temp, rh, ws, pres, z0, z0t, ustarf, & & zlf, cd, ch, ce, cdn, chn, cen, shf, xlhf, xmf, br, ws10) g = 9.8 ! Acceleration due to gravity alpha = 0.011 ! Charnock's formula (After Smith, 1988) rho = 1.225 ! atmospheric density (kg/m3) rhocp = 1231.0 ! rho * Cp rholv = 3013.5 ! rho * Lv xk = 0.4 ! von Karman constant xlv = 2.45 * 1000000. ! Latent heat of vaporization z0i = .0001 ! Initial assumed value u0i = 0. ! Initial assumed wind value z0 = z0i u0 = u0i Call z0tcal (z0, xk, z0t) Call z0tcal (z0, xk, z0q) 100 Continue Call SpHum (temp, rh, pres, sh10, r10) Call SpHum (sst, 98., pres+1.0, sh, r) r10 = r10/1000. ! Mixing Ratio (g/kg) r = r/1000. ! Mixing Ratio (g/kg) pdbt = (temp + 273.16) * ( (1000./pres) ** 0.286 ) psst = (sst + 273.16) vdbt = pdbt * ( 1. + (.61 * r10) ) vsst = psst * ( 1. + (.61 * r) ) ustar = (xk * (ws - u0))/(alog(z/z0)) tstar = (xk * (pdbt - psst))/(alog(z/z0t)) tvstar = (xk * (vdbt - psst))/(alog(z/z0t)) qstar = (xk * (sh10 - sh) * 1000.)/(alog(z/z0q)) If (ws >= 1.0) NI = 5 If (ws < 1.0) NI = 0 Do 210 I = 1, NI xlnum = vdbt * (ustar**2.) xlden = g * xk * tvstar xl = xlnum/xlden zl = z/xl zl10 = 10.0/xlden if (zl < 0.) go to 150 ! For Unstable stratification if (zl == 0.) go to 160 ! For Neutral stratification if (zl > 0.) go to 170 ! For Stable stratification 150 Continue ! For Unstable stratification x = (1. - (16. * zl)) ** .25 x10 = (1. - (16. * zl10)) ** .25 term1 = 2. * alog((1. + x)/2.) term2 = alog ((1. + (x*x))/2.) term3 = -2. * (atan(x)) term4 = 1.5707 psim = term1 + term2 + term3 + term4 psih = 2. * term2 psiq = psih xterm1 = 2. * alog((1. + x10)/2.) xterm2 = alog ((1. + (x10*x10))/2.) xterm3 = -2. * (atan(x10)) xterm4 = 1.5707 psim10 = xterm1 + xterm2 + xterm3 + xterm4 psih10 = 2. * xterm2 psiq10 = psih10 go to 180 160 Continue ! For Neutral stratification psim = 0. psih = 0. psiq = 0. psim10 = 0. psih10 = 0. psiq10 = 0. go to 180 170 Continue ! For Stable stratification psim = -5. * zl psih = psim psiq = psih psim10 = -5. * zl10 psih10 = psim10 psiq10 = psiq10 180 Continue u0 = ustar ustar = ( (ws) * xk ) / ( alog(z/z0) - psim ) tstar = ( (pdbt-psst) * xk ) / ( alog(z/z0t) - psih ) tvstar = ( (vdbt-vsst) * xk ) / ( alog(z/z0t) - psih ) qstar = ( (sh10-sh) * xk * 1000. ) / ( alog(z/z0t) - psiq ) z0c = alpha * (ustar*ustar)/g z0s = (.11 * 14.)/(ustar * 1000000.) z0 = z0c + z0s Call z0tcal (z0, xk, z0t) Call z0tcal (z0, xk, z0q) 210 Continue ws10_1 = ustar/0.4 ws10_2 = alog(10.0/3.0) ws10_3 = ws10_1 * ws10_2 ws10 = ws + ws10_3 ustarp = ustar ustarf = ustar tstarp = tvstar qstarp = qstar zlp = zl zlf = zl10 delta1 = pdbt - psst delta2 = (sh10 - sh) * 1000. cdp = (ustarp/(ws-ustar))**2 chp = (ustarp*tstarp)/(ws * delta1) cep = (ustarp*qstarp)/(ws * delta2) d11 = psst - pdbt d22 = sst - temp xmf = rho * cdp * ((ws) ** 2) shf = rhocp * chp * (ws) * (psst - pdbt) xlhf = rholv * cep * (ws) * (sh - sh10) * 1000. Call cdz (cdp, chp, cep, z, cd10, ch10, ce10, 10.0) cdp = cd10 chp = ch10 cep = ce10 cd = cdp * 1000. ch = chp * 1000. ce = cep * 1000. cdn = ((xk*xk)/(log(10./z0))**2.)*1000. chn = (xk*xk)/( (log(10./z0)) * (log(10./z0t)) ) chn = chn * 1000. cen = (xk*xk)/( (log(10./z0)) * (log(10./z0q)) ) cen = cen * 1000. z0 = z0 * 1000.0 z0t = z0t * 1000.0 br = shf/xlhf Return End Subroutine SpHum (airtmp, relhum, press, spehum, xmrat) ! Subroutine for converting relative humidity to specific humidity eps = 0.622 pow1 = (7.5 * airtmp)/(237.3 + airtmp) xnum1 = 10.**pow1 xnum11 = 6.107 * xnum1 e1 = xnum11*(relhum/100.) spehum = ((eps*e1)/(press - (0.378 * e1))) xmrat = ( (eps*e1)/(press - e1) ) Return End Subroutine z0tcal (zip, vk, zop) vk = 0.4 const1 = 1.15/1000. const2 = alog(10./zip) const3 = const1 * const2 const4 = (vk*vk)/const3 const5 = exp (const4) zop = 10. / const5 Return End !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Subroutine for conversion of air-sea exchange coefficients ...... ! ! from one level 'z1' to 'z2' (Last modified: DBS, August 21, 2006) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine cdz (cd1, ch1, ce1, z1, cd2, ch2, ce2, z2) vk = 0.4 !!! von Korman constant brterm1 = vk/sqrt(cd1) brterm2 = log(z1/z2) brterm = brterm1 - brterm2 cd2 = (vk*vk) / (brterm*brterm) first = vk * sqrt(cd2) brh = (vk/ch1) * sqrt(cd1) bre = (vk/ce1) * sqrt(cd1) brkh = (brh + brterm2) brke = (bre + brterm2) ch2 = first / brkh ce2 = first / brke Return End