source: trunk/src/cor30a.pro

Last change on this file was 204, checked in by pinsard, 10 years ago

fix thanks to coding rules; typo

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