dimension x(4),y(6) data x/5.98,28.82,28.68,20.47/ c obtain constant values for use in bulk algorithm common/para/hu,ht,hq,zi,ps,beta,g,ak,rgas,tok,cpa,sc,pr,z10 call const c call the bulk algorithm call bulk(x(1),x(2),x(3),x(4),y(1),y(2),y(3),y(4),y(5),y(6)) print *,y c 6.419259 5.6866370E-02 90.63017 1.078768 -0.1105995 -3.8373694E-03 stop end c c U. of Arizona (UA) Bulk Aerodynamic Algorithm c 12/22/97 c -------- c Add the `zo/L' term to the equations, because in a global modeling c environment (particularly over land), the absolute value of L could c be smaller than zo (either zom or zoh) (11/04/98) c c Reference: Zeng et al. 1998, Intercomparison of bulk aerodynamic c algorithms for the computation of sea surface fluxes c using the TOGA COARE and TAO data. J. Climate, c 11, 2628-2644. c c For additional information, contact c Prof. Xubin Zeng c Department of Atmospheric Science c PAS Building, #81 c The University of Arizona c Tucson, AZ 85721 c USA c Tel:520-621-4782 c Email:xubin@gogo.atmo.arizona.edu c c input: c u = sqrt(u_x^2 + u_y^2): wind speed in m/s at hu (m) height c ts: surface temperature in (deg C) c t: air temperature in (deg C) at ht (m) height c q: air specific humidity in (g/kg) at hq (m) height c output: c u10: wind speed at 10 meter (m/s) c tau: wind stress (N/m2) c alh: latent heat flux (W/m2) c ash: sensible heat flux (W/m2) c dth: air-surface potential temperature difference c dqh: air-surface specific humidity difference c subroutine bulk(u,ts,t,q,u10,tau,alh,ash,dth,dqh) common/para/hu,ht,hq,zi,ps,beta,g,ak,rgas,tok,cpa,sc,pr,z10 th=(t+tok)*(1000./ps)**(rgas/cpa) !potential T dth=t+0.0098*ht-ts qs=qsat(ts,ps)*0.98 qs=.62197*qs/(ps-0.378*qs) ! in kg/kg dqh=q/1000.-qs thv=th*(1.+0.61*q/1000.) ! virtual potential T dthv=dth*(1.+0.61*q/1000.)+0.61*th*dqh rho=ps*100./(rgas*(ts+tok)*(1.+0.61*qs)) ! density xlv=(2.501-0.00237*ts)*1.e+6 !J/kg c Kinematic viscosity of dry air (m2/s)- Andreas (1989) CRREL Rep. 89-11 visa=1.326e-5*(1+6.542e-3*t+8.301e-6*t*t-4.84e-9*t*t*t) c initial values of u* and convective velocity ustar=0.06 wc=0.5 if(dthv.ge.0.) then um=max(u,0.1) else um=sqrt(u*u+wc*wc) endif c loop to obtain initial and good ustar and zo do i=1,5 zo=0.013*ustar*ustar/g+0.11*visa/ustar ustar=ak*um/alog(hu/zo) enddo rb=g*hu*dthv/(thv*um*um) if(rb.ge.0.) then ! neutral or stable zeta=rb*alog(hu/zo)/(1.-5.*min(rb,0.19)) else !unstable zeta=rb*alog(hu/zo) endif obu=hu/zeta c main iterations (2-10 iterations would be fine) do i=1,10 call rough(zo,zot,zoq,ustar,visa) c wind zeta=hu/obu zetam=1.574 if(zeta.lt.-zetam) then ! zeta < -1 ustar=ak*um/(alog(-zetam*obu/zo)-psi(1,-zetam)+ $ psi(1,zo/obu)+ $ 1.14*((-zeta)**0.333-(zetam)**0.333)) else if (zeta.lt.0.) then ! -1 <= zeta < 0 ustar=ak*um/(alog(hu/zo)-psi(1,zeta)+psi(1,zo/obu)) else if (zeta.le.1.) then ! 0 <= zeta <= 1 ustar=ak*um/(alog(hu/zo) +5.*zeta-5.*zo/obu) else ! 1 < zeta, phi=5+zeta ustar=ak*um/(alog(obu/zo)+5.-5.*zo/obu+ $ (5.*alog(zeta)+zeta-1.)) endif c temperature zeta=ht/obu zetat=0.465 if(zeta.lt.-zetat) then ! zeta < -1 tstar=ak*dth/(alog(-zetat*obu/zot)-psi(2,-zetat)+ $ psi(2,zot/obu)+ $ 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333))) else if (zeta.lt.0.) then ! -1 <= zeta < 0 tstar=ak*dth/(alog(ht/zot)-psi(2,zeta)+psi(2,zot/obu)) else if (zeta.le.1.) then ! 0 <= ztea <= 1 tstar=ak*dth/(alog(ht/zot)+5.*zeta-5.*zot/obu) else ! 1 < zeta, phi=5+zeta tstar=ak*dth/(alog(obu/zot)+5.-5.*zot/obu+ $ (5.*alog(zeta)+zeta-1.)) endif c humidity zeta=hq/obu zetat=0.465 if(zeta.lt.-zetat) then ! zeta < -1 qstar=ak*dqh/(alog(-zetat*obu/zoq)-psi(2,-zetat)+ $ psi(2,zoq/obu)+ $ 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333))) else if (zeta.lt.0.) then ! -1 <= zeta < 0 qstar=ak*dqh/(alog(hq/zoq)-psi(2,zeta)+psi(2,zoq/obu)) else if (zeta.le.1.) then ! 0 <= ztea <= 1 qstar=ak*dqh/(alog(hq/zoq)+5.*zeta-5.*zoq/obu) else ! 1 < zeta, phi=5+zeta qstar=ak*dqh/(alog(obu/zoq)+5.-5.*zoq/obu+ $ (5.*alog(zeta)+zeta-1.)) endif thvstar=tstar*(1.+0.61*q/1000.)+0.61*th*qstar obu=ustar**2*thv/(ak*g*thvstar) if (dthv.ge.0.) then um=max(u,0.1) else wc=beta*(-g*ustar*thvstar*zi/thv)**0.333 um=sqrt(u*u+wc*wc) endif enddo tau=rho*ustar*ustar*u/um alh=-rho*xlv*qstar*ustar ash=-rho*cpa*tstar*ustar c x and y components of tau: c taux=rho*ustar*ustar*u_x/um c tauy=rho*ustar*ustar*u_y/um c 10-meter wind (without w_* part) zeta=z10/obu if(zeta.lt.0.) then u10=u+(ustar/ak)*(alog(z10/hu)-(psi(1,zeta)- $ psi(1,hu/obu))) else u10=u+(ustar/ak)*(alog(z10/hu)+5.*zeta-5.*hu/obu) endif return end c constants subroutine const common/para/hu,ht,hq,zi,ps,beta,g,ak,rgas,tok,cpa,sc,pr,z10 hu=4. ! m ht=3. ! m hq=3. ! m zi=1000. ! m (pbl height) ps=1013. ! mb (sfc pressure) beta=1. ! - (in computing W_*) g=9.8 ! ms-2 ak=0.4 ! - (Von Karman constant) rgas=287.1 ! J(kg)^-1K^-1 gas constant for dry air tok=273.16 ! Celsius to Kelvin cpa=1004.67 ! J/kg/K specific heat of dry air sc=0.60 ! =kinematic viscosity/water vapor molecular diffusivity ! the Schmidt number (from Garrratt, 1992, p.286) pr=0.71 ! =nu/thermal diffusivity (the Prandtl number) z10=10. ! m (reference height) return end c stability function for rb < 0 function psi(k,zeta) chik=(1.-16*zeta)**0.25 if(k.eq.1) then psi=2.*alog((1.+chik)*0.5)+ $ alog((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.) else psi=2.*alog((1.+chik*chik)*0.5) endif return end c Tetens' formula for saturation vp Buck(1981) JAM 20, 1527-1532 c p in mb, t in C, and qsat in mb function qsat(t,p) qsat = (1.0007+3.46e-6*p)*6.1121*exp(17.502*t/(240.97+t)) return end c our formulation for zo,zot,zoq subroutine rough(zo,zot,zoq,ustar,visa) common/para/hu,ht,hq,zi,ps,beta,g,ak,rgas,tok,cpa,sc,pr,z10 zo=0.013*ustar*ustar/g+0.11*visa/ustar re=ustar*zo/visa xq=2.67*re**0.25-2.57 xt= xq zoq=zo/exp(xq) zot=zo/exp(xt) return end