C C ***************************************************************** C * u* = u*n hypohtesis * C ***************************************************************** C * * C * Bulk flux algorithm based on the u*=u*n hypothesis when * C * converting to neutral parameters (constant flux hypothesis) * C * * C * The parametrization of the bulk coefficients is from * C * Dupuis H. et al. 1997, JGR 102,C9 21115-21129. * C * * C * The convection parameterization of Fairall et al 1996 is also * C * implemented and may or not be enabled in the code * C * * C * Input parameters : * C * zam (m) Common height of sensors * C * Ua (m/s) Wind speed at height zam * C * Qa (g/kg) Specific humidity at zam * C * Ta (Celsius) Air temperature at zam * C * Pa (hPa) Air pressure at zam * C * Zi (m) Mixed layer height * C * ts (Celsius) Sea surface temperature * C * pcv if set to 1, convection parametrization * C * is enabled * C * * C * Output variables : * C * xle (W/m2) Latent heat flux * C * h (W/m2) Sensible heat flux * C * tau (N/m2) Momentum flux * C * zsurl=Zam/L L is the Monin Obukhov Length * C ****************Denis Bourras 1997/11/24************************* C ****************Last revision 1999/04/20************************* C C C************* ATTENTION C The parametrizations of bulk coefficients maybe easily changed in the code. C However, this program is not intended to work with parametrizations such as C Large and Pond 1982 in which the hypothesis to convert to neutral parameters C is a constant wind speed. The hypothesis in this code is u*=u*n, not U=Un C C More details maybe found in C Bourras, D., 1999, Estimation of the latent heat flux over oceans from C satellite data, Thesis, 203 pp., Paris VI University, France (in rench). C Bourras, D., 2000, Calculation of turbulent fluxes C Technical report CETP/RI/2/2000 (in french) C available from bourras@pacific.jpl.nasa.gov. C Grachev et al. 1998, Boundary Layer meteorology, 86, 257-278. C C ************ ATTENTION 2 C At this time nothing prevents from using simultaneously Dupuis et al. and C and Fairall et al parametrizations, but it might be redundent since the C Dupuis coefficients increase at low wind speeds to account for convection. C Up to you to decide (pcv=1 enables Fairall et al. parametrization). C C *********** References C Fairall et al. 1996, JGR, 101,C2, 3747-3764. C Large and Pond, 1982, Journal of Physical Oceanography, 11, 324-336. C C C*********************************************************************** subroutine fbustarus(ua99,qa99,ts99,ta99,pa99, X Zam99,Zi99,xle,h,tau,zsurl,pcv) C REAL ua99,qa99,ts99,ta99,pa99,Zam99,Zi99,xle,h,tau,zsurl integer pcv C ua=ua99 qa=qa99 ts=ts99 ta=ta99 pa=pa99 Zam=zam99 Zi=zi99 rhomoy=1.2 g=9.81 rspeca=287. pi=3.1415925 vkarman=0.4 C print *,' convection : ',pcv C C***** Input parameter conversion C qa=qa/1000. ! in Kgvap/Kg air total ts=ts+273.15 ! in K ta=ta+273.15 pa=pa*100. ! in Pa C C***** Atmospheric pressure at surface level C ps=pa+rhomoy*g*zam ! rho is assumed constant C print*,'ps=',ps C C***** Density at zam. C C* Mixing ratio at zam ra=qa/(1.-qa) ! in Kgvapor/Kgdryair C print*,'ra',ra C C***** virtual temperature at zam tva=ta*(1.+0.608*qa) C print*,'tva=',tva C rhoa=pa/(rspeca*tva) C print*,'rhoa=',rhoa C C***** First estimation of the Monin obukhov length C C**** Potential temperatures at surface and zam C**** and Delta T pot. C tpota=ta*((100.*1000./pa)**0.286) tpots=ts*((100.*1000./ps)**0.286) dtpot=tpota-tpots ! different convention from Large and Pond 1982 if (abs(dtpot).lt.1e-3) then xle=-50. h=-50. tau=-50. zsurl=-50. stop endif C print*,'tpota=',tpota,' tpots=',tpots,' delta T pot',dtpot C C******* Virtual potential temperature at zam C tvpota=tpota*(1+0.608*qa) C C**** saturating water vapor pressures at surface and zam C ewsats=610.78*exp(17.269*(ts-273.15)/(ts-35.86)) ewsata=610.78*exp(17.269*(ta-273.15)/(ta-35.86)) C print*,'ewsats=',ewsats,' ewsata=',ewsata C C**** Mixing ratios at saturation C rsats=0.622*ewsats/(ps-ewsats) rsata=0.622*ewsata/(pa-ewsata) C print*,'rsats=',rsats,' rsata=',rsata C C**** Specific humidities at saturation and delta Q C qsats=rsats/(1+rsats) qsata=rsata/(1+rsata) C print*,'qsats=',qsats,' qsata=',qsata dq=qa-0.98*qsats ! different convention as in Large and Pond 1982 C print*,'dq=',dq if (dq.gt.0) then C print*,'saturation !' dq=0. endif C C***** and the M-O ratio C t0=ta*(1+1.7e-6*ta*qa) if (dtpot.lt.0) then C unstable case zsurl=-100*zam/((ua**2)*t0)*(-dtpot-1.7e-6*(t0**2)*dq) C free convection wstar=0.3 wg=1.2*wstar uabase=ua if (pcv.eq.1) ua=sqrt((uabase**2)+(wg**2)) else zsurl=-70*zam/((ua**2)*t0)*(-dtpot-2.5e-6*(t0**2)*dq) endif C print*,'z/L=',zsurl C C** FIRST GUESS C u10=ua CD10=0.11e-3 CH10=0.11e-3 CE10=0.11e-3 C ustar=sqrt(CD10)*u10 ! as first approximations tstar=CH10*u10*dtpot/ustar qstar=CE10*u10*dq/ustar C C************************ C* Main loop * C************************ C boucle=1 C print*,'iteration number ',boucle C 10 boucle=boucle+1 ustarsauve=ustar ! The stop test tstarsauve=tstar qstarsauve=qstar C if (zsurl.gt.0.2) then C print*,'Z/L outside range (>0.2)' endif if (zsurl.lt.-1.) then C print*,'Strong convection' endif C C***** Structure functions psi (z/L) C C**** Psim stands for tau, PsiT(=psiQ) for H and LE C**** psim(10/L) if (zsurl.gt.0.) then psim=-7.*zsurl psit=psim psim10surl=-7.*(10./(zam/zsurl)) psit10surl=psim10surl endif if (zsurl.le.0.) then x=(1.-16.*zsurl)**0.25 psim=2.*alog((1.+x)/2.)+alog((1.+(x**2))/2.) X -2.*atan(x)+PI/2. psit=2.*alog((1.+(x**2))/2.) x10surl=(1.-16.* (10./(zam/zsurl)) )**0.25 psim10surl=2.*alog((1.+x10surl)/2.)+ X alog((1.+(x10surl**2))/2.)-2.*atan(x10surl)+PI/2. psit10surl=2.*alog((1.+(x10surl**2))/2.) C******* convective part for psim(z/L) et psit(z/L) frac=1./(1.+(zsurl)**2) ! Fairall et al. 96 C y=(1.-16.*zsurl)**(1./3.) psic=1.5*alog((y**2+y+1.)/3.)-sqrt(3.)* X atan((2.*y+1.)/sqrt(3.))+pi/sqrt(3.) C print*,'frac=',frac,' psic=',psic if (pcv.eq.1) then psim=frac*psim+(1.-frac)*psic psit=frac*psit+(1.-frac)*psic endif C******* convective part for psim(10/L) frac=1./(1.+(10./(zam/zsurl))**2) ! Fairall et al. 96 y=(1.-16.* (10./(zam/zsurl)) )**(1./3.) psic=1.5*alog((y**2+y+1.)/3.) X -sqrt(3.)*atan((2.*y+1.)/sqrt(3.))+pi/sqrt(3.) if (pcv.eq.1) then psim10surl=frac*psim10surl+(1.-frac)*psic psit10surl=frac*psit10surl+(1.-frac)*psic endif endif C print*,'psim=',psim,' psit=',psit,' psim(10/L)=',psim10surl C C****** variables at 10 meters for evaluation of z/L C dtpot10=dtpot+tstar/vkarman*( -alog(10./zam)-psit10surl+psit ) tpot10=dtpot10+ts ! potential temperature at 10 m in K dq10=dq+qstar/vkarman*(alog(10./zam)-psit10surl+psit ) q10=0.98*qsats+dq10 !specific humidity a 10 metres tvpot10=tpot10*(1.+0.608*q10) ! pot virt temp C C************ EQUIVALENT NEUTRAL Wind, hum, and temp at 10 METRES C u10n=ua+ustar*(-alog(zam/10.)+psim)/vkarman C !equivalent neutral wind at 10m u10=ua+ustar*(-alog(zam/10.)+psim-psim10surl)/vkarman C ! to compute CD10 C ! (for intercomparisons) dq10n=dq+qstar*(-alog(zam/10.)+psit)/vkarman dtpot10n=dtpot+tstar*(-alog(zam/10.)+psit)/vkarman C if(u10n.lt.0) print*,'U10neq is below 0, C X convection too strong, iter:',boucle C C*********** H. Dupuis parametrisation C if (u10n.le.5.85) then ! sofia + semaphore low winds cdn10 = 0.0117 * (u10n **(-2)) + 0.000668 ! Yelland et Taylor strong winds else cdn10 = 0.00007 *u10n + 0.0006 endif if (u10n.le.5.2) then ! sofia + semaphore cen10= 0.00279 /u10n +0.000658 else cen10= 0.00279 /(5.2) + 0.000658 endif if (zsurl.le.0.) then if (u10n.le.5.2) then ! sofia + semaphore chn10= 0.00279 /u10n +0.000658 else chn10= 0.00279 /(5.2) + 0.000658 endif else chn10=7e-4 endif C C C***** EVALUATION OF USTAR, QSTAR, et TSTAR C ustar=sqrt(CDN10)*u10n CD10=ustar*ustar/u10/u10 ! for intercomparisons tstar=CHN10*u10n*dtpot10n/ustar qstar=CEN10*u10n*dq10n/ustar cd1=qstar*sqrt(cd10)/dq10 C print*,'ustar=',ustar,' tstar=',tstar,' qstar=',qstar C C****** New calculation of z/L for this iteration C******* This time, L is estimated at 10 meters C zsurlsauve=zsurl ! for convergence test zsurl= zam*vkarman*g*( tstar*(1.+0.608*q10)+0.608 X *tpot10*qstar )/(tvpot10*ustar**2) C print*,'z/L=',zsurl C C***** Free convection hypothesis C if (dtpot.lt.0.) then buoy=-g/tvpota*ustar* X (tstar*(1+0.608*qa)+0.608*tpota*qstar ) mbuoy=1 if (buoy.lt.0.) mbuoy=-1 wstar=((zi*abs(buoy)))**(1./3.) wg=1.2*wstar if (pcv.eq.1) ua=sqrt((uabase**2)+(wg**2)) C print*,'buoy=',buoy,' wstar=',wstar,' wg=',wg,' ua=',ua endif C if (boucle.gt.50) goto 20 if (zsurl.gt.0.25) boucle=100. if (abs(ustar-ustarsauve).ge.1.e-4) goto 10 if (abs(tstar-tstarsauve).ge.1.e-6) goto 10 if (abs(qstar-qstarsauve).ge.1.e-7) goto 10 if (abs(zsurlsauve-zsurl).ge.1.e-3) goto 10 C 20 if (boucle.gt.50) then C print*,'No solution found' C print*,xle,h,tau,zsurl xle=-999. h=-999. tau=-999. zsurl=-999. goto 30 endif C**** two following formulae From P. Taylor algorithm cpa=1004.67*(1+0.00084*qa) cle=(2.501-0.00237*(ts-273.16))*1000000 C C********** Final expression of the fluxes H=-rhoa*cpa*ustar*tstar xle=-rhoa*cle*ustar*qstar tau=rhoa*ustar*ustar C 30 print*,'TAU=',tau,' LE=',xle,' H=',h,' z/L=',zsurl 30 if (zsurl.lt.1000) goto 32 31 xle=-999. h=-999. tau=-999. zsurl=-999. C 32 return END REAL xle,h,tau,zsurl,u10,u10n REAL ua,qa,ts,ta,pa,Zam,Zi integer pcv xle=0. h=0. tau=0. zsurl=0. ua=15. qa=7. ts=24. ta=20. pa=1013.15 Zam=17. Zi=2000. pcv=0 call fbustarus(ua,qa,ts,ta,pa,Zam,Zi,xle,h,tau,zsurl,pcv) C C C print*,'---------------' if (pcv.eq.1) then print *, 'Convection parametrisation' else print *, 'No convection parametrisation' endif print*,'UA ',ua,' (m/s-1)' print*,'QA ',qa,' (g/kg-1)' print*,'TS ',ts,' (C)' print*,'TA ',ta,' (C)' print*,'PA ',pa,' (mb)' print*,'za ',zam,' (m)' print*,'zi ',zi,' (m)' print*,'---------------' print*,'LE ',xle,' (W/m2)' print*,'HS ',h,' (W/m2)' print*,'Tau ',tau,' (N/m2)' print*,'z/L ',zsurl end C