[6] | 1 | ;+ |
---|
| 2 | ; |
---|
[12] | 3 | ; ========== |
---|
| 4 | ; cor30a.pro |
---|
| 5 | ; ========== |
---|
[6] | 6 | ; |
---|
[12] | 7 | ; .. function:: cor30a(u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave) |
---|
[6] | 8 | ; |
---|
[81] | 9 | ; DESCRIPTION |
---|
| 10 | ; =========== |
---|
[6] | 11 | ; |
---|
[5] | 12 | ; COARE v3 algorithm to compute fluxes |
---|
| 13 | ; |
---|
[6] | 14 | ; version with shortened iteration modified Rt and Rq |
---|
[5] | 15 | ; |
---|
[6] | 16 | ; uses wave information wave period in s and wave ht in m |
---|
| 17 | ; no wave, standard coare 2.6 charnock: jwave=0 |
---|
[5] | 18 | ; |
---|
[6] | 19 | ; Oost et al. zo=50/2/pi L (u*/c)^4.5 if jwave=1 |
---|
[5] | 20 | ; |
---|
[6] | 21 | ; taylor and yelland zo=1200 h*(L/h)^4.5 jwave=2 |
---|
[5] | 22 | ; |
---|
[6] | 23 | ; :param u: wind speed (m/s) at height zu (m) |
---|
| 24 | ; :param us: surface current speed in the wind direction (m/s) |
---|
| 25 | ; :param ts: bulk water temperature (C) if jcool=1, interface water T if jcool=0 |
---|
| 26 | ; :param t: bulk air temperature (C), height zt |
---|
| 27 | ; :param Qs: bulk water spec hum (g/kg) if jcool=1, ... |
---|
| 28 | ; :param Q: bulk air spec hum (g/kg), height zq |
---|
| 29 | ; :param Rs: downward solar flux (W/m^2) (modified because of cool skin) |
---|
| 30 | ; :param Rl: downard IR flux (W/m^2) (modified because of cool skin) |
---|
| 31 | ; :param rain: rain rate (mm/hr) |
---|
| 32 | ; :param zi: PBL depth (m) |
---|
| 33 | ; :param P: Atmos surface pressure (mb) |
---|
| 34 | ; :param zu: wind speed measurement height (m) |
---|
| 35 | ; :param zt: air T measurement height (m) |
---|
| 36 | ; :param zq: air q measurement height (m) |
---|
| 37 | ; :param lat: latitude (deg, N=+) |
---|
| 38 | ; :param jcool: implement cool calculation skin switch, 0=no, 1=yes |
---|
| 39 | ; :param jwave: implement wave dependent roughness model |
---|
| 40 | ; :param twave: wave period (s) |
---|
| 41 | ; :param hwave: wave height (m) |
---|
| 42 | ; |
---|
| 43 | ; EXAMPLES |
---|
| 44 | ; ======== |
---|
| 45 | ; |
---|
| 46 | ; :: |
---|
| 47 | ; |
---|
| 48 | ; u us ts ta qs qa Qsw IRd r pbl Ps zu zt zq lat |
---|
| 49 | ; |
---|
| 50 | ; 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.) |
---|
| 51 | ; |
---|
| 52 | ; Result with these sample values with Matlab code:: |
---|
| 53 | ; |
---|
| 54 | ; 8.64830 101.640 0.0352910 2.17780e-05 0.000115000 0.000115000 |
---|
| 55 | ; -29.5800 0.175430 -0.0423670 -0.000205610 0.250950 0.000351300 |
---|
| 56 | ; 0.000969740 0.00000 8.11390e-05 0.000997340 0.00121410 0.00121400 |
---|
| 57 | ; 0.000941350 0.00107910 0.00107910 0.780060 |
---|
| 58 | ; |
---|
| 59 | ; Result obtained with this idl routine:: |
---|
| 60 | ; |
---|
| 61 | ; 8.64829 101.640 0.0352913 2.17780e-05 0.000115000 0.000115000 |
---|
| 62 | ; -29.5802 0.175432 -0.0423667 -0.000205610 0.250948 0.000351304 |
---|
| 63 | ; 0.000969737 0.00000 8.11394e-05 0.000997343 0.00121407 0.00121400 |
---|
| 64 | ; 0.000941351 0.00107908 0.00107908 0.780056 |
---|
| 65 | ; |
---|
[5] | 66 | ; Maximum error on any parameter: .002 % validated! |
---|
| 67 | ; |
---|
[12] | 68 | ; SEE ALSO |
---|
| 69 | ; ======== |
---|
| 70 | ; |
---|
[175] | 71 | ; used by :func:`tropflux` |
---|
[12] | 72 | ; |
---|
[90] | 73 | ; call :func:`psiu`, :func:`psit`. |
---|
| 74 | ; |
---|
[6] | 75 | ; TODO |
---|
| 76 | ; ==== |
---|
| 77 | ; |
---|
| 78 | ; coding rules |
---|
| 79 | ; |
---|
[81] | 80 | ; check for the new module provide by pk 20110811 by mail to fp |
---|
| 81 | ; "You can use this coare algorithm in your depositary. It is same as the one you |
---|
| 82 | ; have, but with some functions activated, which will be useful in the next stage." |
---|
| 83 | ; |
---|
[6] | 84 | ; EVOLUTIONS |
---|
| 85 | ; ========== |
---|
| 86 | ; |
---|
[81] | 87 | ; $Id$ |
---|
| 88 | ; |
---|
| 89 | ; $URL$ |
---|
| 90 | ; |
---|
[97] | 91 | ; - fplod 20110830T085907Z aedon.locean-ipsl.upmc.fr (Darwin) |
---|
| 92 | ; |
---|
| 93 | ; * add ug |
---|
| 94 | ; |
---|
[90] | 95 | ; - fplod 20110822T085404Z aedon.locean-ipsl.upmc.fr (Darwin) |
---|
| 96 | ; |
---|
| 97 | ; * add ref to psiu and psit called her and provided by pk today |
---|
| 98 | ; |
---|
[6] | 99 | ; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin) |
---|
| 100 | ; |
---|
| 101 | ; * minimal header |
---|
| 102 | ; |
---|
| 103 | ; - pbk 2008 |
---|
| 104 | ; |
---|
| 105 | ; * creation |
---|
| 106 | ; |
---|
| 107 | ;- |
---|
| 108 | function cor30a, u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave |
---|
[5] | 109 | |
---|
| 110 | Qs=Qs/1000. |
---|
| 111 | Q=Q/1000. |
---|
| 112 | |
---|
| 113 | ;*********** set constants ************* |
---|
| 114 | pi=!pi |
---|
[6] | 115 | Beta=1.2 |
---|
| 116 | von=0.4 |
---|
| 117 | fdg=1.00 |
---|
| 118 | tdk=273.16 |
---|
| 119 | ;grav=grv(lat) |
---|
[5] | 120 | grav=9.8 |
---|
| 121 | ;************* air constants ************ |
---|
[6] | 122 | Rgas=287.1 |
---|
| 123 | LLe=(2.501-.00237*ts)*1e6 |
---|
| 124 | cpa=1004.67 |
---|
| 125 | cpv=cpa*(1+0.84*Q) |
---|
| 126 | rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*Q)) |
---|
| 127 | visa=1.326e-5*(1+6.542e-3*t+8.301e-6*t*t-4.84e-9*t*t*t) |
---|
[5] | 128 | ;************ cool skin constants ******* |
---|
[6] | 129 | Al=2.1e-5*(ts+3.2)^0.79 |
---|
| 130 | be=0.026 |
---|
| 131 | cpw=4000 |
---|
| 132 | rhow=1022 |
---|
| 133 | visw=1e-6 |
---|
| 134 | tcw=0.6 |
---|
| 135 | bigc=16*grav*cpw*(rhow*visw)^3/(tcw*tcw*rhoa*rhoa) |
---|
| 136 | wetc=0.622*LLe*Qs/(Rgas*(ts+tdk)^2) |
---|
| 137 | |
---|
[5] | 138 | ;*************** wave parameters ********* |
---|
[6] | 139 | lwave=grav/2/pi*twave^2 |
---|
| 140 | cwave=grav/2/pi*twave |
---|
| 141 | |
---|
[5] | 142 | ;************** compute aux stuff ******* |
---|
[6] | 143 | Rns=Rs*.945 |
---|
| 144 | Rnl=0.97*(5.67e-8*(ts-0.3*jcool+tdk)^4-Rl) |
---|
[5] | 145 | |
---|
| 146 | |
---|
| 147 | |
---|
| 148 | ;*************** Begin bulk loop ******* |
---|
[6] | 149 | |
---|
[5] | 150 | ;*************** first guess ************ |
---|
[6] | 151 | du=u-us |
---|
| 152 | dt=ts-t-.0098*zt |
---|
| 153 | dq=Qs-Q |
---|
| 154 | ta=t+tdk |
---|
| 155 | ug=0.5 |
---|
| 156 | dter=0.3 |
---|
| 157 | dqer=wetc*dter |
---|
| 158 | ut=sqrt(du*du+ug*ug) |
---|
| 159 | u10=ut*alog(10./1e-4)/alog(zu/1e-4) |
---|
| 160 | usr=.035*u10 |
---|
| 161 | zo10=0.011*usr*usr/grav+0.11*visa/usr |
---|
| 162 | Cd10=(von/alog(10./zo10))^2 |
---|
| 163 | Ch10=0.00115 |
---|
| 164 | Ct10=Ch10/sqrt(Cd10) |
---|
| 165 | zot10=10/exp(von/Ct10) |
---|
| 166 | Cd=(von/alog(zu/zo10))^2 |
---|
| 167 | Ct=von/alog(zt/zot10) |
---|
| 168 | CC=von*Ct/Cd |
---|
| 169 | Ribcu=-zu/zi/.004/Beta^3 |
---|
| 170 | Ribu=-grav*zu/ta*((dt-dter*jcool)+.61*ta*dq)/ut^2 |
---|
| 171 | nits=3 |
---|
[97] | 172 | ;if (Ribu le 0.) then begin |
---|
| 173 | ; zetu=CC*Ribu/(1+Ribu/Ribcu) |
---|
| 174 | ;endif else begin |
---|
| 175 | ; zetu=CC*Ribu*(1+27./9*Ribu/CC) |
---|
| 176 | ;endelse |
---|
[5] | 177 | sw=(Ribu le 0.) |
---|
| 178 | zetu=sw*(CC*Ribu/(1+Ribu/Ribcu))+(1-sw)*(CC*Ribu*(1+27./9*Ribu/CC)) |
---|
[97] | 179 | ; |
---|
[6] | 180 | L10=zu/zetu |
---|
[97] | 181 | ;if (zetu gt 50 ) then nits=1 |
---|
[6] | 182 | usr=ut*von/(alog(zu/zo10)-psiu(zu/L10)) |
---|
| 183 | tsr=-(dt-dter*jcool)*von*fdg/(alog(zt/zot10)-psit(zt/L10)) |
---|
| 184 | qsr=-(dq-wetc*dter*jcool)*von*fdg/(alog(zq/zot10)-psit(zq/L10)) |
---|
[5] | 185 | |
---|
[6] | 186 | tkt=.001 |
---|
| 187 | |
---|
[97] | 188 | ;charn=0.011 |
---|
| 189 | ;if (ut gt 10.) then charn=0.011+(ut-10)/(18.-10)*(0.018-0.011) |
---|
| 190 | ;if (ut gt 18.) then charn=0.018 |
---|
[5] | 191 | charn=(((0.011+(ut-10)/(18.-10)*(0.018-0.011)) > .011) < .018) |
---|
[97] | 192 | ; |
---|
[6] | 193 | |
---|
[5] | 194 | ;*************** bulk loop ************ |
---|
| 195 | for i=1,nits do begin |
---|
[6] | 196 | zet=von*grav*zu/ta*(tsr*(1+0.61*Q)+.61*ta*qsr)/(usr*usr)/(1+0.61*Q) |
---|
[5] | 197 | case jwave of |
---|
[6] | 198 | 0: zo=charn*usr*usr/grav+0.11*visa/usr |
---|
[5] | 199 | 1: zo=50./2/pi*lwave*(usr/cwave)^4.5+0.11*visa/usr ;Oost et al |
---|
| 200 | 2: zo=1200*hwave*(hwave/lwave)^4.5+0.11*visa/usr ;Taylor and Yelland |
---|
| 201 | endcase |
---|
[6] | 202 | rr=zo*usr/visa |
---|
| 203 | L=zu/zet |
---|
[97] | 204 | ;zoq=min([1.15e-4,5.5e-5/rr^.6]) |
---|
[5] | 205 | zoq=(5.5e-5/rr^.6 < 1.15e-4) |
---|
[97] | 206 | ; |
---|
[6] | 207 | zot=zoq |
---|
| 208 | usr=ut*von/(alog(zu/zo)-psiu(zu/L)) |
---|
| 209 | tsr=-(dt-dter*jcool)*von*fdg/(alog(zt/zot)-psit(zt/L)) |
---|
| 210 | qsr=-(dq-wetc*dter*jcool)*von*fdg/(alog(zq/zoq)-psit(zq/L)) |
---|
| 211 | Bf=-grav/ta*usr*(tsr+.61*ta*qsr) |
---|
[97] | 212 | ;if (Bf gt 0) then begin |
---|
| 213 | ; ug=Beta*(Bf*zi)^.333 |
---|
| 214 | ;endif else begin |
---|
| 215 | ; ug=.2 |
---|
| 216 | ;endelse |
---|
[5] | 217 | sw=(Bf gt 0) |
---|
| 218 | ug=sw*(Beta*(Bf*zi)^.333)+(1-sw)*.2 |
---|
[97] | 219 | ; |
---|
[6] | 220 | ut=sqrt(du*du+ug*ug) |
---|
| 221 | Rnl=0.97*(5.67e-8*(ts-dter*jcool+tdk)^4-Rl) |
---|
| 222 | hsb=-rhoa*cpa*usr*tsr |
---|
| 223 | hlb=-rhoa*LLe*usr*qsr |
---|
| 224 | qout=Rnl+hsb+hlb |
---|
[5] | 225 | dels=Rns*(.065+11*tkt-6.6e-5/tkt*(1-exp(-tkt/8.0e-4))) ; Eq.16 Shortwave |
---|
[6] | 226 | qcol=qout-dels |
---|
[5] | 227 | alq=Al*qcol+be*hlb*cpw/LLe ; Eq. 7 Buoy flux water |
---|
| 228 | |
---|
[97] | 229 | ; if (alq gt 0) then begin |
---|
| 230 | ; xlamx=6./(1+(bigc*alq/usr^4)^.75)^.333 ; Eq 13 Saunders |
---|
| 231 | ; tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr) ;Eq.11 Sub. thk |
---|
| 232 | ; endif else begin |
---|
| 233 | ; xlamx=6.0 |
---|
| 234 | ; tkt=min([.01,xlamx*visw/(sqrt(rhoa/rhow)*usr)]) ;Eq.11 Sub. thk |
---|
| 235 | ; tkt=(xlamx*visw/(sqrt(rhoa/rhow)*usr) < .01) |
---|
| 236 | ; |
---|
| 237 | ; endelse |
---|
[5] | 238 | sw=(alq gt 0) |
---|
| 239 | xlamx=sw*(6./(1+(bigc*alq/usr^4)^.75)^.333)+(1-sw)*6.0 |
---|
| 240 | tkt=sw*(xlamx*visw/(sqrt(rhoa/rhow)*usr))+(1-sw)*(xlamx*visw/(sqrt(rhoa/rhow)*usr) < .01) |
---|
[97] | 241 | ; |
---|
[6] | 242 | |
---|
[5] | 243 | dter=qcol*tkt/tcw ; Eq.12 Cool skin |
---|
[6] | 244 | dqer=wetc*dter |
---|
| 245 | |
---|
[5] | 246 | endfor ;bulk iter loop |
---|
| 247 | |
---|
| 248 | tau=rhoa*usr*usr*du/ut ;stress |
---|
| 249 | hsb=-rhoa*cpa*usr*tsr ;sens |
---|
| 250 | hlb=-rhoa*LLe*usr*qsr ;lat |
---|
| 251 | ; net solar Rns |
---|
| 252 | ; net lw Rnl |
---|
| 253 | |
---|
| 254 | ;**************** rain heat flux ******** |
---|
| 255 | dwat=2.11e-5*((t+tdk)/tdk)^1.94 ;! water vapour diffusivity |
---|
| 256 | dtmp=(1.+3.309e-3*t-1.44e-6*t*t)*0.02411/(rhoa*cpa) ;!heat diffusivity |
---|
| 257 | alfac= 1/(1+(wetc*LLe*dwat)/(cpa*dtmp)) ;! wet bulb factor |
---|
| 258 | RF= rain*alfac*cpw*((ts-t-dter*jcool)+(Qs-Q-dqer*jcool)*LLe/cpa)/3600. |
---|
| 259 | |
---|
| 260 | |
---|
| 261 | ;y=[[Rns],[-1.*Rnl],[-1.*hlb],[-1.*hsb],[-1.*RF],[tau]] |
---|
| 262 | |
---|
[97] | 263 | ;**************** Webb et al. correection ************ |
---|
[5] | 264 | ;wbar=1.61*hlb/LLe/(1+1.61*Q)/rhoa+hsb/rhoa/cpa/ta ;formulation in hlb already includes webb |
---|
[6] | 265 | ;hl_webb=rhoa*wbar*Q*LLe |
---|
[97] | 266 | ;************** compute transfer coeffs relative to ut @meas. ht ********** |
---|
| 267 | ;Cd=tau/rhoa/ut/max([.1,du]) |
---|
[6] | 268 | ;Cd=tau/rhoa/ut/(du > .1) |
---|
[97] | 269 | ; |
---|
[6] | 270 | Ch=-usr*tsr/ut/(dt-dter*jcool) |
---|
| 271 | Ce=-usr*qsr/(dq-dqer*jcool)/ut |
---|
[97] | 272 | ;************ 10-m neutral coeff realtive to ut ******** |
---|
[6] | 273 | ;Cdn_10=von*von/alog(10./zo)/alog(10./zo) |
---|
| 274 | ;Chn_10=von*von*fdg/alog(10./zo)/alog(10./zot) |
---|
| 275 | ;Cen_10=von*von*fdg/alog(10./zo)/alog(10./zoq) |
---|
[5] | 276 | |
---|
[97] | 277 | y=[[Rns],[-1.*Rnl],[-1.*hlb],[-1.*hsb],[-1.*RF],[tau],[Ch],[Ce],[ug]] |
---|
| 278 | ;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 ] |
---|
[5] | 279 | ; 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
---|
| 280 | ; hsb= sensible heat flux (w/m^2) |
---|
| 281 | ; hlb= latent heat flux (w/m^2) |
---|
| 282 | ; RF= rain heat flux(w/m^2) |
---|
[97] | 283 | ; wbar= webb mean w (m/s) |
---|
[5] | 284 | ; tau= stress (nt/m^2) |
---|
| 285 | ; zo= velocity roughness length (m) |
---|
| 286 | ; zot temperature roughness length (m) |
---|
| 287 | ; zoq= moisture roughness length (m) |
---|
| 288 | ; L= Monin_Obukhov stability length |
---|
| 289 | ; usr= turbulent friction velocity (m/s), including gustiness |
---|
| 290 | ; tsr temperature scaling parameter (K) |
---|
| 291 | ; qsr humidity scaling parameter (g/g) |
---|
| 292 | ; dter= cool skin temperature depression (K) |
---|
| 293 | ; dqer= cool skin humidity depression (g/g) |
---|
| 294 | ; tkt= cool skin thickness (m) |
---|
| 295 | ; Cd= velocity drag coefficient at zu, referenced to u |
---|
| 296 | ; Ch= heat transfer coefficient at zt |
---|
| 297 | ; Ce= moisture transfer coefficient at zq |
---|
| 298 | ; Cdn_10= 10-m velocity drag coeeficient, including gustiness |
---|
| 299 | ; Chn_10= 10-m heat transfer coeeficient, including gustiness |
---|
| 300 | ; Cen_10= 10-m humidity transfer coeeficient, including gustiness |
---|
| 301 | |
---|
| 302 | return, y |
---|
| 303 | |
---|
| 304 | end |
---|