!23456789012345678901234567890123456789012345678901234567890123456789012 ! ! compute height-adjusted wind speed and air temperature based on the input ! reference heights, wind speed, air temperature and sea surface temperature. ! ! Written by Dr. Pat Fitzpatrick (fitz@gri.msstate.edu) ! and Yee Lau (lau@gri.msstate.edu) ! ! Affiliation: ! ! GeoSystems Research Institute ! Mississippi State University ! ! ! Input arguments: ! 1. input file with 3 columns as follows: ! i. reference wind speed in m/s ! ii. reference air temperature in K ! iii. sea surface temperature in K ! ! 2. wind speed reference height in meter ! ! 3. air temperature reference height in meter ! ! 4. target adjusted height in meter ! ! 5. output file name ! The output file will have 2 columns as follows: ! i. adjusted wind speed in m/s ! ii. adjusted air temperature in K ! ! all units are metric (meters, m/s, Kelvin) ! ! k = von Karman constant = 0.4 ! uStar = friction wind speed ! thetaStar = potential temperature scale ! Spdh = wind speed at adjusted height Zhm ! Spdr = wind speed at reference height Zrm ! sst = sea surface temperature ! thetaVr = potential virtual temperature at reference height Zr ! tempVr = air temperature at reference height Zrt ! temph = air temperature at adjusted height Zht ! Z0t = thermal roughness length ! Z0m = momentum roughness length ! Zrm = reference level height for wind ! Zrt = reference level height for temperature ! Prt = turbulent Prandt number = 0.95 ! Rib = bulk Richardson number ! g = gravity ! implicit none character*200 inputFile, outputFile character*10 str_Zrm, str_Zrt, str_Zhm, str_Zht, str_Spdr character*10 str_tempVr, str_sst integer iargc, nlines real Gm, Rib, Zrm, Zrt, Gh, Spdh, Spdr, tmp, Zhm, Zht, Z0m, Z0t real thetaStar, uStar, thetaVr, sst, L, temph, tempVr real xleft, xmid, xright real x, Psim, Psit, stabFm, stabFt external f real, PARAMETER :: pi = 3.14159265359 real, PARAMETER :: k = 0.4 real, PARAMETER :: g = 9.8 ! real, PARAMETER :: Z0m = 0.01 ! real, PARAMETER :: Z0t = 0.0001 real, PARAMETER :: vis = 14.0E-6 real, PARAMETER :: alpha = 0.015 real, PARAMETER :: alpha2 = 5.0 real, PARAMETER :: gamma = 16.0 real, PARAMETER :: Prt = 0.95 ! ! note: Jacobsen's book suggests alpha2=6.0, and gamma=19.3...although the stability equations ! also appear different. See Chapter 8, pages 238-249. Jacobsen also notes the coefficients ! Pr, alpha2, and gamma depend on k ! ! Reference: ! ! Jacobsen, M.Z., 2005: Fundamentals of Atmospheric Modeling, 2nd edition. Cambridge ! University Press, 813 pp. ! ! Also see his class notes at: http://web.stanford.edu/group/efmh/jacobson/FAMbook2dEd/index.html ! ! I've also seen gamma=15, alpha2=4.7, and alpha2=7. ! !*************************************************************** ! check the number of input arguments and assign them to variables !*************************************************************** if (iargc() .ne. 6) then print *, 'usage: louis_windTempProfile_InOutFiles' // & ' inputFile windSpdRefHeightInM airTmpRefHeightInM' // & ' windSpdAdjustedHeightInM airTmpAdjustedHeightInM outputFile' print *, 'where inputFile should have 3 columns (no header) =>' // & ' 1. reference wind Speed in m/s ' // & ' 2. reference air temperature in K ' // & ' 3. sea surface temperature in K' print *, 'The outputFile will have 2 columns (no header) =>' // & ' 1. adjusted wind speed in m/s ' // & ' 2. adjusted air temperature in K' stop endif call getarg(1,inputFile) call getarg(2,str_Zrm) call getarg(3,str_Zrt) call getarg(4,str_Zhm) call getarg(5,str_Zht) call getarg(6,outputFile) read(str_Zrm, *) Zrm read(str_Zrt, *) Zrt read(str_Zhm, *) Zhm read(str_Zht, *) Zht open(unit=10, file=inputFile, status='old') open(unit=20, file=outputFile, status='unknown') nlines = 0 do read(10,*,end=50) Spdr, tempVr, sst nlines = nlines + 1 ! ! future change: if specific humidity q is available, change temperature ! to virtual temperature (where q has units g/kg). If its kg/kg, then use ! 0.609 . ! ! tempVr=tempVr(1.+0.000609*q) ! ! ! Potential temperature differs from the actual temperature by 0.1C per 10 meters ! Adjust for adiabatic adjustment to surface ! thetaVr=tempVr+0.01*Zrt !******************************************************** ! Iterate for initial ustar based on Charnock's equation !******************************************************** xleft = 0.000001 xright = 5.0 call bisect(f,xleft,xright,uStar,Zrm,Spdr) ! print *, "bisect uStar=", uStar !******************************************************** ! Compute roughness length Z0m using Charnock's equation !******************************************************** ! Z0m = alpha * uStar * uStar / g Z0m = alpha * uStar * uStar / g + 0.11*vis/uStar !***************************************************************** ! Compute thermal roughness length Z0t from Eq. 10 in Smith (1988) ! ! Smith, S.S., 1988: Coefficients for sea surface wind stress, heat flux, ! and wind profiles as a function of wind speed and temperature, Journal ! of Geophysical Research, 93, No. C12, 15467-15472. ! !***************************************************************** ! Z0t = Z0m Z0t=Zrt/exp((k**2)/(0.00115*log(Zrt/Z0m))) ! Z0t=Zrt/exp((k**2)/(0.001*log(Zrt/Z0m))) ! print *, "Z0m=", Z0m," Z0t=", Z0t !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! note.....may need to iterate for z0m and/or z0t ! in future versions...or recompute below !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !*************************************************************** ! calculate Rib, Gm, Gh, uStar, thetaStar and L ! ! uStar and thetaStar are computed using the method of Louis (1979). ! This method avoids iterating for a solution. ! ! Reference ! ! Louis, J.-F., 1979: A parametric model of vertical eddy fluxes ! in the atmosphere. Boundary-Layer Meteorology, 17, 187-202. ! !*************************************************************** Rib = (g * (thetaVr - sst) * (Zrm - Z0m)**2.0) / & (0.5*(sst+thetaVr) * (Spdr**2) * (Zrt-Z0t)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! not sure if some of these should be Zrt. setting all to Zrm for now. ! same uncertainty with Z0t. We may need a tmpm and tmpt in the future. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (Rib <= 0.0) then tmp = k * k * ((abs(Rib)*Zrm/Z0m)**0.5) / (log(Zrm/Z0m)**2) Gm = 1.0 - ((9.4 * Rib) / (1.0 + 70.0*tmp)) Gh = 1.0 - ((9.4 * Rib) / (1.0 + 50.0*tmp)) else Gm = 1.0 / ((1.0 + 4.7 * Rib)**2.0) Gh = Gm endif uStar = k * abs(Spdr) / log(Zrm/Z0m) * sqrt(Gm) thetaStar = ((k**2 * abs(Spdr) * (thetaVr - sst)) / & (uStar * Prt * log(Zrm/Z0m)**2)) * Gh L = (uStar**2 * 0.5*(sst + thetaVr) ) / (k * g * thetaStar) ! print *, "Rib=", Rib ! print *, "Gm=", Gm ! print *, "Gh=", Gh ! print *, "uStar=", uStar ! print *, "thetaStar=", thetaStar ! print *, "L=", L ! print *, "Zrm/L=", Zrm/L !*************************************************************** ! calculate integrated stability functions for stable & unstable stratifications ! ! For stable stratification (air temp > water temp), use Dyer (1974) ! ! For stable stratification (air temp < water temp), use Paulson (1970) ! ! These were suggested by Smith (1988) ! ! References: ! ! Dyer, A.J., 1974: A review of flux-profile relationships. Boundary Layer ! Meteorology, 7, 363-372. ! ! Dyers, A.J., and B.B. Hicks, 1970: Flux-gradient relationships in the ! constant flux layer, Quarterly Journal of the Royal Meteorological Society, ! 96, 715-721. ! ! Paulson, C.A., 1970: The mathematical representation of wind speed and ! temperature profiles in the unstable atmospheric surface layer, Journal ! Applied Meteorology, 9, 857-861. ! ! Smith, S.S., 1988: Coefficients for sea surface wind stress, heat flux, ! and wind profiles as a function of wind speed and temperature, Journal ! of Geophysical Research, 93, No. C12, 15467-15472. ! !*************************************************************** stabFm = Zhm / L stabFt = Zht / L ! print *, "stabFm=", stabFm ! print *, "stabFt=", stabFt if (stabFm >= (-0.01).and.stabFm <= 0.01) then Psim = 0.0 else if (stabFm > 0.0) then Psim = -alpha2 * stabFm else x = (1.0 - gamma * stabFm)**0.25 Psim = 2.0 * log((1.+x)/2.0) + log((1.+x*x)/2.0) - & 2.0 * atan(x) + (pi/2.0) endif if (stabFt >= (-0.01).and.stabFt <= 0.01) then Psit = 0.0 else if (stabFt > 0.0) then Psit = -alpha2 * stabFt else x = (1.0 - gamma * stabFt)**0.25 Psit = 2.0 * log((1.+x*x)/2.0) endif ! print *, "Psit=", Psit ! print *, "Psim=", Psim !************************************************************* ! Recompute Z0m and Z0t based on new ustar !************************************************************* Z0m = alpha * uStar * uStar / g + 0.11*vis/uStar Z0t=Zrt/exp((k**2)/(0.00115*log(Zrt/Z0m))) ! print *, "Recomputed Z0m=", Z0m," Z0t=", Z0t !************************************************************* ! Compute wind profile at Zhm !************************************************************* Spdh=(uStar/k)*(log(Zhm/Z0m)-Psim) ! print *, "Adjusted speed=", Spdh, ! & ' m/s at adjusted height=',Zhm ,' m' ! print *,'Spdh/Spdr=',Spdh/Spdr !************************************************************* ! Compute potential temperature profile at Zht, and adjust to ! temperature by subtracting 0.01*Zht ! ! Will need to change from virtual temperature to temperature in the ! future if q is available ! !************************************************************* temph=sst+(Prt*thetaStar/k)*(log(Zht/Z0t)-Psit)-0.01*Zht ! print *, 'Adjusted temperature=', temph, ! & ' at adjusted height=',Zht ,' m' ! print *,'temph/tempVr=',temph/tempVr write(20, *) Spdh, temph end do 50 close(10) close(20) end