;+ ; ; ========== ; cor30a.pro ; ========== ; ; .. function:: cor30a(u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave) ; ; DESCRIPTION ; =========== ; ; COARE v3 algorithm to compute fluxes ; ; version with shortened iteration modified Rt and Rq ; ; uses wave information wave period in s and wave ht in m ; no wave, standard coare 2.6 charnock: jwave=0 ; ; Oost et al. zo=50/2/pi L (u*/c)^4.5 if jwave=1 ; ; taylor and yelland zo=1200 h*(L/h)^4.5 jwave=2 ; ; :param u: wind speed (m/s) at height zu (m) ; :param us: surface current speed in the wind direction (m/s) ; :param ts: bulk water temperature (C) if jcool=1, interface water T if jcool=0 ; :param t: bulk air temperature (C), height zt ; :param Qs: bulk water spec hum (g/kg) if jcool=1, ... ; :param Q: bulk air spec hum (g/kg), height zq ; :param Rs: downward solar flux (W/m^2) (modified because of cool skin) ; :param Rl: downard IR flux (W/m^2) (modified because of cool skin) ; :param rain: rain rate (mm/hr) ; :param zi: PBL depth (m) ; :param P: Atmos surface pressure (mb) ; :param zu: wind speed measurement height (m) ; :param zt: air T measurement height (m) ; :param zq: air q measurement height (m) ; :param lat: latitude (deg, N=+) ; :param jcool: implement cool calculation skin switch, 0=no, 1=yes ; :param jwave: implement wave dependent roughness model ; :param twave: wave period (s) ; :param hwave: wave height (m) ; ; EXAMPLES ; ======== ; ; :: ; ; u us ts ta qs qa Qsw IRd r pbl Ps zu zt zq lat ; ; IDL> x=cor30a(5.5,0,28.7,27.2,24.2,18.5,141.,419.,0.,600.,1010.,15.,15.,15.,0.,1,1,5.,1.) ; ; Result with these sample values with Matlab code:: ; ; 8.64830 101.640 0.0352910 2.17780e-05 0.000115000 0.000115000 ; -29.5800 0.175430 -0.0423670 -0.000205610 0.250950 0.000351300 ; 0.000969740 0.00000 8.11390e-05 0.000997340 0.00121410 0.00121400 ; 0.000941350 0.00107910 0.00107910 0.780060 ; ; Result obtained with this idl routine:: ; ; 8.64829 101.640 0.0352913 2.17780e-05 0.000115000 0.000115000 ; -29.5802 0.175432 -0.0423667 -0.000205610 0.250948 0.000351304 ; 0.000969737 0.00000 8.11394e-05 0.000997343 0.00121407 0.00121400 ; 0.000941351 0.00107908 0.00107908 0.780056 ; ; Maximum error on any parameter: .002 % validated! ; ; SEE ALSO ; ======== ; ; used by :ref:`TropFlux_19890101_20091231.pro` ; ; TODO ; ==== ; ; coding rules ; ; check for the new module provide by pk 20110811 by mail to fp ; "You can use this coare algorithm in your depositary. It is same as the one you ; have, but with some functions activated, which will be useful in the next stage." ; ; EVOLUTIONS ; ========== ; ; $Id$ ; ; $URL$ ; ; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin) ; ; * minimal header ; ; - pbk 2008 ; ; * creation ; ;- function cor30a, u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave Qs=Qs/1000. Q=Q/1000. ;*********** set constants ************* pi=!pi Beta=1.2 von=0.4 fdg=1.00 tdk=273.16 ;grav=grv(lat) grav=9.8 ;************* air constants ************ Rgas=287.1 LLe=(2.501-.00237*ts)*1e6 cpa=1004.67 cpv=cpa*(1+0.84*Q) rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*Q)) visa=1.326e-5*(1+6.542e-3*t+8.301e-6*t*t-4.84e-9*t*t*t) ;************ cool skin constants ******* Al=2.1e-5*(ts+3.2)^0.79 be=0.026 cpw=4000 rhow=1022 visw=1e-6 tcw=0.6 bigc=16*grav*cpw*(rhow*visw)^3/(tcw*tcw*rhoa*rhoa) wetc=0.622*LLe*Qs/(Rgas*(ts+tdk)^2) ;*************** wave parameters ********* lwave=grav/2/pi*twave^2 cwave=grav/2/pi*twave ;************** compute aux stuff ******* Rns=Rs*.945 Rnl=0.97*(5.67e-8*(ts-0.3*jcool+tdk)^4-Rl) ;*************** Begin bulk loop ******* ;*************** first guess ************ du=u-us dt=ts-t-.0098*zt dq=Qs-Q ta=t+tdk ug=0.5 dter=0.3 dqer=wetc*dter ut=sqrt(du*du+ug*ug) u10=ut*alog(10./1e-4)/alog(zu/1e-4) usr=.035*u10 zo10=0.011*usr*usr/grav+0.11*visa/usr Cd10=(von/alog(10./zo10))^2 Ch10=0.00115 Ct10=Ch10/sqrt(Cd10) zot10=10/exp(von/Ct10) Cd=(von/alog(zu/zo10))^2 Ct=von/alog(zt/zot10) CC=von*Ct/Cd Ribcu=-zu/zi/.004/Beta^3 Ribu=-grav*zu/ta*((dt-dter*jcool)+.61*ta*dq)/ut^2 nits=3 ;;if (Ribu le 0.) then begin ;; zetu=CC*Ribu/(1+Ribu/Ribcu) ;;endif else begin ;; zetu=CC*Ribu*(1+27./9*Ribu/CC) ;;endelse sw=(Ribu le 0.) zetu=sw*(CC*Ribu/(1+Ribu/Ribcu))+(1-sw)*(CC*Ribu*(1+27./9*Ribu/CC)) ;; L10=zu/zetu ;;if (zetu gt 50 ) then nits=1 usr=ut*von/(alog(zu/zo10)-psiu(zu/L10)) tsr=-(dt-dter*jcool)*von*fdg/(alog(zt/zot10)-psit(zt/L10)) qsr=-(dq-wetc*dter*jcool)*von*fdg/(alog(zq/zot10)-psit(zq/L10)) tkt=.001 ;;charn=0.011 ;;if (ut gt 10.) then charn=0.011+(ut-10)/(18.-10)*(0.018-0.011) ;;if (ut gt 18.) then charn=0.018 charn=(((0.011+(ut-10)/(18.-10)*(0.018-0.011)) > .011) < .018) ;; ;*************** bulk loop ************ for i=1,nits do begin zet=von*grav*zu/ta*(tsr*(1+0.61*Q)+.61*ta*qsr)/(usr*usr)/(1+0.61*Q) case jwave of 0: zo=charn*usr*usr/grav+0.11*visa/usr 1: zo=50./2/pi*lwave*(usr/cwave)^4.5+0.11*visa/usr ;Oost et al 2: zo=1200*hwave*(hwave/lwave)^4.5+0.11*visa/usr ;Taylor and Yelland endcase rr=zo*usr/visa L=zu/zet ;;zoq=min([1.15e-4,5.5e-5/rr^.6]) zoq=(5.5e-5/rr^.6 < 1.15e-4) ;; zot=zoq usr=ut*von/(alog(zu/zo)-psiu(zu/L)) tsr=-(dt-dter*jcool)*von*fdg/(alog(zt/zot)-psit(zt/L)) qsr=-(dq-wetc*dter*jcool)*von*fdg/(alog(zq/zoq)-psit(zq/L)) Bf=-grav/ta*usr*(tsr+.61*ta*qsr) ;;if (Bf gt 0) then begin ;; ug=Beta*(Bf*zi)^.333 ;;endif else begin ;; ug=.2 ;;endelse sw=(Bf gt 0) ug=sw*(Beta*(Bf*zi)^.333)+(1-sw)*.2 ;; ut=sqrt(du*du+ug*ug) Rnl=0.97*(5.67e-8*(ts-dter*jcool+tdk)^4-Rl) hsb=-rhoa*cpa*usr*tsr hlb=-rhoa*LLe*usr*qsr qout=Rnl+hsb+hlb dels=Rns*(.065+11*tkt-6.6e-5/tkt*(1-exp(-tkt/8.0e-4))) ; Eq.16 Shortwave qcol=qout-dels alq=Al*qcol+be*hlb*cpw/LLe ; Eq. 7 Buoy flux water ;; if (alq gt 0) then begin ;; xlamx=6./(1+(bigc*alq/usr^4)^.75)^.333 ; Eq 13 Saunders ;; tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr) ;Eq.11 Sub. thk ;; endif else begin ;; xlamx=6.0 ;;;; tkt=min([.01,xlamx*visw/(sqrt(rhoa/rhow)*usr)]) ;Eq.11 Sub. thk ;; tkt=(xlamx*visw/(sqrt(rhoa/rhow)*usr) < .01) ;;;; ;; endelse sw=(alq gt 0) xlamx=sw*(6./(1+(bigc*alq/usr^4)^.75)^.333)+(1-sw)*6.0 tkt=sw*(xlamx*visw/(sqrt(rhoa/rhow)*usr))+(1-sw)*(xlamx*visw/(sqrt(rhoa/rhow)*usr) < .01) ;; dter=qcol*tkt/tcw ; Eq.12 Cool skin dqer=wetc*dter endfor ;bulk iter loop tau=rhoa*usr*usr*du/ut ;stress hsb=-rhoa*cpa*usr*tsr ;sens hlb=-rhoa*LLe*usr*qsr ;lat ; net solar Rns ; net lw Rnl ;**************** rain heat flux ******** dwat=2.11e-5*((t+tdk)/tdk)^1.94 ;! water vapour diffusivity dtmp=(1.+3.309e-3*t-1.44e-6*t*t)*0.02411/(rhoa*cpa) ;!heat diffusivity alfac= 1/(1+(wetc*LLe*dwat)/(cpa*dtmp)) ;! wet bulb factor RF= rain*alfac*cpw*((ts-t-dter*jcool)+(Qs-Q-dqer*jcool)*LLe/cpa)/3600. ;y=[[Rns],[-1.*Rnl],[-1.*hlb],[-1.*hsb],[-1.*RF],[tau]] ;;**************** Webb et al. correection ************ ;wbar=1.61*hlb/LLe/(1+1.61*Q)/rhoa+hsb/rhoa/cpa/ta ;formulation in hlb already includes webb ;hl_webb=rhoa*wbar*Q*LLe ;;************** compute transfer coeffs relative to ut @meas. ht ********** ;;Cd=tau/rhoa/ut/max([.1,du]) ;Cd=tau/rhoa/ut/(du > .1) ;;; Ch=-usr*tsr/ut/(dt-dter*jcool) Ce=-usr*qsr/(dq-dqer*jcool)/ut ;;************ 10-m neutral coeff realtive to ut ******** ;Cdn_10=von*von/alog(10./zo)/alog(10./zo) ;Chn_10=von*von*fdg/alog(10./zo)/alog(10./zot) ;Cen_10=von*von*fdg/alog(10./zo)/alog(10./zoq) y=[[Rns],[-1.*Rnl],[-1.*hlb],[-1.*hsb],[-1.*RF],[tau],[Ch],[Ce]] ;;y=[hsb,hlb,tau,zo,zot,zoq,L,usr,tsr,qsr,dter,dqer,tkt,RF,wbar,Cd,Ch,Ce,Cdn_10,Chn_10,Cen_10,ug ] ; 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 ; hsb= sensible heat flux (w/m^2) ; hlb= latent heat flux (w/m^2) ; RF= rain heat flux(w/m^2) ; wbar= webb mean w (m/s) ; tau= stress (nt/m^2) ; zo= velocity roughness length (m) ; zot temperature roughness length (m) ; zoq= moisture roughness length (m) ; L= Monin_Obukhov stability length ; usr= turbulent friction velocity (m/s), including gustiness ; tsr temperature scaling parameter (K) ; qsr humidity scaling parameter (g/g) ; dter= cool skin temperature depression (K) ; dqer= cool skin humidity depression (g/g) ; tkt= cool skin thickness (m) ; Cd= velocity drag coefficient at zu, referenced to u ; Ch= heat transfer coefficient at zt ; Ce= moisture transfer coefficient at zq ; Cdn_10= 10-m velocity drag coeeficient, including gustiness ; Chn_10= 10-m heat transfer coeeficient, including gustiness ; Cen_10= 10-m humidity transfer coeeficient, including gustiness return, y end