Context Navigation

pbl_parameters.F @ 508

Last change on this file since 508 was 496, checked in by acolaitis, 13 years ago

M 489 libf/phymars/thermcell_main_mars.F90
^{--------------------> New parameters for thermals, following latest version of the relevant article}

M 489 libf/phymars/physiq.F
^{--------------------> Minor modifications}

D 489 libf/phymars/surflayer_interpol.F
A 0 libf/phymars/pbl_parameters.F
^{--------------------> Replaced surflayer_interpol.F by a cleaner and richer version called pbl_parameters.F}

This routine is the base of what will later be implemented in the MCD and as a tool to
compute surface properties from output fields (so that it will ultimately disappear from libf
and might become an utility)

M 489 libf/phymars/vdif_cd.F
^{--------------------> Minor modification to the Richardson computation}

M 489 libf/phymars/vdifc.F
^{--------------------> Minor modification to the aerodynamic conductances computation, replacing wmax by hfmax, which}

is a better proxy for convective activity. Might be replaced by w_star later.

File size: 15.0 KB

Line
1	SUBROUTINE pbl_parameters(ngrid,nlay,pz0,
2	& pg,pz,pu,pv,wmax,hfmax,zmax,pts,ph,z_out,
3	& Teta_out,u_out,ustar,tstar,wstar,L_mo,vhf,vvv)
4	IMPLICIT NONE
5	!=======================================================================
6	!
7	! Anlysis of the PBL from input temperature, wind field and thermals outputs.
8	!
9	! -------
10	!
11	! Author: Arnaud Colaitis 09/01/12
12	! -------
13	!
14	! Arguments:
15	! ----------
16	!
17	! inputs:
18	! ------
19	! ngrid size of the horizontal grid
20	! nlay size of the vertical grid
21	! pz0(ngrid) surface roughness length
22	! pg gravity (m s -2)
23	! pz(ngrid,nlay) height of layers
24	! pu(ngrid,nlay) u component of the wind
25	! pv(ngrid,nlay) v component of the wind
26	! wmax(ngrid) maximum vertical velocity in thermals (might not be needed
27	! if the computation of w* works)
28	! hfmax(ngrid) maximum vertical turbulent heat flux in thermals
29	! zmax(ngrid) height reached by the thermals (pbl height)
30	! pts(ngrid) surface temperature
31	! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa
32	! z_out height of interpolation
33	!
34	! outputs:
35	! ------
36	!
37	! Teta_out(ngrid) interpolated teta
38	! u_out(ngrid) interpolated u
39	! ustar(ngrid) friction velocity
40	! tstar(ngrid) friction temperature
41	! wstar(ngrid) free convection velocity
42	! L_mo(ngrid) monin_obukhov length
43	!
44	!
45	!=======================================================================
46	!
47	!-----------------------------------------------------------------------
48	! Declarations:
49	! -------------
50
51	#include "comcstfi.h"
52	#include "callkeys.h"
53
54	! Arguments:
55	! ----------
56
57	INTEGER, INTENT(IN) :: ngrid,nlay
58	REAL, INTENT(IN) :: pz0(ngrid)
59	REAL, INTENT(IN) :: pg,pz(ngrid,nlay)
60	REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
61	REAL, INTENT(IN) :: wmax(ngrid),hfmax(ngrid),zmax(ngrid)
62	REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
63	REAL, INTENT(IN) :: z_out
64	REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)
65	REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid),wstar(ngrid)
66	REAL, INTENT(OUT) :: L_mo(ngrid)
67
68	! Local:
69	! ------
70
71	INTEGER ig,k
72	REAL karman,nu
73	DATA karman,nu/.41,0.001/
74	SAVE karman,nu
75
76	! Local(2):
77	! ---------
78
79	REAL zout
80	REAL rib(ngrid) ! Bulk Richardson number
81	REAL fm(ngrid) ! stability function for momentum
82	REAL fh(ngrid) ! stability function for heat
83	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
84	! phim = 1+betamzeta or (1-bmzeta)**am
85	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
86	REAL betam, betah, alphah, bm, bh, lambda
87	! ah and am are assumed to be -0.25 and -0.5 respectively
88	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient
89	REAL pz0t ! initial thermal roughness length. (local)
90	REAL ric ! critical richardson number
91	REAL reynolds(ngrid) ! reynolds number for UBL
92	REAL prandtl(ngrid) ! prandtl number for UBL
93	REAL pz0tcomp(ngrid) ! computed z0t
94	REAL ite
95	REAL residual,zcd0,z1
96	REAL pcdv(ngrid),pcdh(ngrid)
97	REAL zu2(ngrid) ! Large-scale wind at first layer
98	REAL pbl_teta(ngrid) ! mixed-layer potential temperature
99	INTEGER pbl_height_index(ngrid) ! index of nearest vertical grid point for zmax
100	REAL dteta(ngrid,nlay),x(ngrid) ! potential temperature gradient and z/zi
101	REAL dvhf(ngrid),dvvv(ngrid) ! dimensionless vertical heat flux and
102	! dimensionless vertical velocity variance
103	REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux and vertical velocity variance
104	INTEGER ii(1)
105	! temporary
106	INTEGER dimout
107
108	!------------------------------------------------------------------------
109	!------------------------------------------------------------------------
110	! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
111	!------------------------------------------------------------------------
112	!------------------------------------------------------------------------
113
114	c Initialisation :
115
116	L_mo(:)=0.
117	ustar(:)=0.
118	tstar(:)=0.
119	zout=z_out
120	reynolds(:)=0.
121	pz0t = 0.
122	pz0tcomp(:) = 0.1*pz0(:)
123	rib(:)=0.
124	pcdv(:)=0.
125	pcdh(:)=0.
126
127	! this formulation assumes alphah=1., implying betah=betam
128	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
129
130	bm=16. !UBL
131	bh=16. !UBL
132	alphah=1.
133	betam=5. !SBL
134	betah=5. !SBL
135	lambda=(sqrt(bh/bm))/alphah
136	ric=betah/(betam**2)
137
138	DO ig=1,ngrid
139	ite=0.
140	residual=abs(pz0tcomp(ig)-pz0t)
141	zu2(ig)=MAX(pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1),wmax(ig)**2)
142
143	DO WHILE((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.))
144
145	pz0t=pz0tcomp(ig)
146	IF (zu2(ig) .ne. 0.) THEN
147	! Richardson number formulation proposed by D.E. England et al. (1995)
148	rib(ig) = (pg/ph(ig,1))
149	& sqrt(pz(ig,1)pz0(ig))
150	& (((log(pz(ig,1)/pz0(ig)))*2)/(log(pz(ig,1)/pz0t)))
151	& *(ph(ig,1)-pts(ig))/zu2(ig)
152	ELSE
153	print*,'warning, infinite Richardson at surface'
154	print*,pu(ig,1),pv(ig,1)
155	rib(ig) = ric
156	ENDIF
157
158	z1z0=pz(ig,1)/pz0(ig)
159	z1z0t=pz(ig,1)/pz0t
160
161	cdn(ig)=karman/log(z1z0)
162	cdn(ig)=cdn(ig)*cdn(ig)
163	chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
164
165	! STABLE BOUNDARY LAYER :
166	IF (rib(ig) .gt. 0.) THEN
167	! From D.E. England et al. (95)
168	prandtl(ig)=1.
169	if(rib(ig) .lt. ric) then
170	! Assuming alphah=1. and bh=bm for stable conditions :
171	fm(ig)=((ric-rib(ig))/ric)**2
172	fh(ig)=((ric-rib(ig))/ric)**2
173	else
174	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
175	fm(ig)=0.
176	fh(ig)=0.
177	endif
178	! UNSTABLE BOUNDARY LAYER :
179	ELSE
180	! From D.E. England et al. (95)
181	fm(ig)=sqrt(1.-lambdabmrib(ig))
182	fh(ig)=(1./alphah)((1.-lambdabhrib(ig))0.5)
183	& (1.-lambdabmrib(ig))**0.25
184	prandtl(ig)=alphah((1.-lambdabmrib(ig))*0.25)/
185	& ((1.-lambdabhrib(ig))**0.5)
186	ENDIF
187
188	reynolds(ig)=karman*sqrt(fm(ig))
189	& *sqrt(zu2(ig))
190	& pz0(ig)/(log(z1z0)nu)
191	pz0tcomp(ig)=pz0(ig)exp(-karman7.3*
192	& (reynolds(ig)*0.25)(prandtl(ig)**0.5))
193	residual = abs(pz0t-pz0tcomp(ig))
194	ite = ite+1
195
196	ENDDO ! of while
197	pz0t=0.
198
199	! Drag computation:
200
201	pcdv(ig)=cdn(ig)*fm(ig)
202	pcdh(ig)=chn(ig)*fh(ig)
203
204	ENDDO ! of ngrid
205
206	!------------------------------------------------------------------------
207	!------------------------------------------------------------------------
208	! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
209	!------------------------------------------------------------------------
210	!------------------------------------------------------------------------
211
212
213	! Large-scale wind at first layer (without gustiness) and
214	! u* theta* computation
215
216	DO ig=1,ngrid
217	IF (rib(ig) .ge. ric) THEN
218	ustar(ig)=0.
219	tstar(ig)=0.
220	ELSE
221	ustar(ig)=sqrt(pcdv(ig))
222	& sqrt(pu(ig,1)pu(ig,1)+pv(ig,1)*pv(ig,1))
223	tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
224	& /sqrt(pcdv(ig))
225	ENDIF
226	ENDDO
227
228	! Monin Obukhov length:
229
230	DO ig=1,ngrid
231	IF (rib(ig) .gt. ric) THEN
232	L_mo(ig)=0.
233	ELSE
234	L_mo(ig)=pts(ig)(ustar(ig)2)/(pgkarman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL
235	ENDIF
236	ENDDO
237
238	! Interpolation:
239
240	DO ig=1,ngrid
241	IF(zout .lt. pz0tcomp(ig)) THEN
242	u_out(ig)=0.
243	Teta_out(ig)=pts(ig)
244	ELSEIF (L_mo(ig) .gt. 0.) THEN
245	u_out(ig)=(ustar(ig)/karman)*log(zout/pz0(ig)) +
246	& 5.(ustar(ig)/(karmanL_mo(ig)))*(zout-pz0(ig))
247	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))
248	& *log(zout/pz0tcomp(ig)) +
249	& 5.(tstar(ig)/(prandtl(ig)karman*L_mo(ig)))
250	& *(zout-pz0tcomp(ig))
251	ELSEIF (L_mo(ig) .lt. 0.) THEN
252
253	IF(L_mo(ig) .gt. -1000.) THEN
254
255	u_out(ig)=(ustar(ig)/karman)*(
256	& 2.atan((1.-16.zout/L_mo(ig))**0.25)
257	& -2.atan((1.-16.pz0(ig)/L_mo(ig))**0.25)
258	& -2.log(1.+(1.-16.zout/L_mo(ig))**0.25)
259	& +2.log(1.+(1.-16.pz0(ig)/L_mo(ig))**0.25)
260	& - log(1.+sqrt(1.-16.*zout/L_mo(ig)))
261	& + log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
262	& + log(zout/pz0(ig))
263	& )
264
265	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
266	& 2.log(1.+sqrt(1.-16.pz0tcomp(ig)/L_mo(ig)))
267	& -2.log(1.+sqrt(1.-16.zout/L_mo(ig)))
268	& + log(zout/pz0tcomp(ig))
269	& )
270
271	ELSE
272
273	! We have to treat the case where L is very large and negative,
274	! corresponding to a very nearly stable atmosphere (but not quite !)
275	! i.e. teta* <0 and almost zero. We choose the cutoff
276	! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
277	! between the two expression for z/L = -1/1000 is -1.54324*10^-9
278	! (we do that to avoid using r*4 precision, otherwise, we get -inf values)
279
280	u_out(ig)=(ustar(ig)/karman)*(
281	& (4./L_mo(ig))*(zout-pz0(ig))
282	& + (20./(L_mo(ig))*2)(zout2-pz0(ig)2)
283	& + (160./(L_mo(ig))*3)(zout3-pz0(ig)3)
284	& + log(zout/pz0(ig))
285	& )
286
287	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
288	& (8./L_mo(ig))*(zout-pz0tcomp(ig))
289	& + (48./(L_mo(ig))*2)(zout2-pz0tcomp(ig)2)
290	& + (1280./(3.(L_mo(ig))3))(zout3-pz0tcomp(ig)3)
291	& + log(zout/pz0tcomp(ig))
292	& )
293
294	ENDIF
295	ELSE
296	u_out(ig)=0.
297	Teta_out(ig)=pts(ig)
298	ENDIF
299	IF(zout .lt. pz0(ig)) THEN
300	u_out(ig)=0.
301	ENDIF
302	ENDDO
303
304	! when using convective adjustment without thermals, a vertical potential temperature
305	! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
306	! interpolated at any height in the surface layer is theta at the first level.
307
308	IF ((.not.calltherm).and.(calladj)) THEN
309	Teta_out(:)=ph(:,1)
310	ENDIF
311
312	!------------------------------------------------------------------------
313	!------------------------------------------------------------------------
314	! PART III : WSTAR COMPUTATION
315	!------------------------------------------------------------------------
316	!------------------------------------------------------------------------
317
318	! Detection of the mixed-layer potential temperature
319	! ------------
320
321	! Nearest index for the pbl height
322
323	pbl_height_index(:)=1
324
325	DO k=1,nlay-1
326	DO ig=1, ngrid
327	IF (abs(zmax(ig)-pz(ig,k)) .lt.
328	& abs(zmax(ig)-pz(ig,pbl_height_index(ig)))) THEN
329	pbl_height_index(ig)=k
330	ENDIF
331	ENDDO
332	ENDDO
333
334	! Potential temperature gradient
335
336	dteta(:,nlay)=0.
337	DO k=1,nlay-1
338	DO ig=1, ngrid
339	dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(pz(ig,k+1)-pz(ig,k))
340	ENDDO
341	ENDDO
342
343	! Computation of the pbl mixed layer temperature
344
345	DO ig=1, ngrid
346	ii=MINLOC(abs(dteta(ig,1:pbl_height_index(ig))))
347	pbl_teta(ig) = ph(ig,ii(1))
348	ENDDO
349
350	! We follow Spiga et. al 2010 (QJRMS)
351	! ------------
352
353	DO ig=1, ngrid
354	IF (zmax(ig) .gt. 0.) THEN
355	wstar(ig)=(pgzmax(ig)hfmax(ig)/pbl_teta(ig))**(1./3.)
356	ELSE
357	wstar(ig)=0.
358	ENDIF
359	ENDDO
360
361	!------------------------------------------------------------------------
362	!------------------------------------------------------------------------
363	! PART IV : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
364	!------------------------------------------------------------------------
365	!------------------------------------------------------------------------
366
367	! We follow Spiga et. al 2010 (QJRMS)
368	! ------------
369
370	DO ig=1, ngrid
371	IF (zmax(ig) .gt. 0.) THEN
372	x(ig) = zout/zmax(ig)
373	ELSE
374	x(ig) = 999.
375	ENDIF
376	ENDDO
377
378	DO ig=1, ngrid
379	! dimensionless vertical heat flux
380	IF (x(ig) .le. 0.3) THEN
381	dvhf(ig) = ((-3.85/log(x(ig)))+0.07*log(x(ig)))
382	& exp(-4.61x(ig))
383	ELSEIF (x(ig) .le. 1.) THEN
384	dvhf(ig) = -1.52*x(ig) + 1.24
385	ELSE
386	dvhf(ig) = 0.
387	ENDIF
388	! dimensionless vertical velocity variance
389	IF (x(ig) .le. 1.) THEN
390	dvvv(ig) = 2.05(x(ig)(2./3.))(1.-0.64x(ig))*2
391	ELSE
392	dvvv(ig) = 0.
393	ENDIF
394	ENDDO
395
396	vhf(:) = dvhf(:)*hfmax(:)
397	vvv(:) = dvvv(:)(wstar(:))*2
398
399	!------------------------------------------------------------------------
400	!------------------------------------------------------------------------
401	! OUTPUTS
402	!------------------------------------------------------------------------
403	!------------------------------------------------------------------------
404
405	IF (ngrid .eq. 1) THEN
406	dimout=0
407	ELSE
408	dimout=2
409	ENDIF
410
411	call WRITEDIAGFI(ngrid,'Teta_out',
412	& 'potential temperature at z_out','K',
413	& dimout,Teta_out)
414	call WRITEDIAGFI(ngrid,'u_out',
415	& 'horizontal velocity norm at z_out','m/s',
416	& dimout,u_out)
417	call WRITEDIAGFI(ngrid,'u_star',
418	& 'friction velocity','m/s',
419	& dimout,ustar)
420	call WRITEDIAGFI(ngrid,'teta_star',
421	& 'friction potential temperature','K',
422	& dimout,tstar)
423	call WRITEDIAGFI(ngrid,'L',
424	& 'Monin Obukhov length','m',
425	& dimout,L_mo)
426	call WRITEDIAGFI(ngrid,'w_star',
427	& 'Free convection velocity','m',
428	& dimout,wstar)
429	call WRITEDIAGFI(ngrid,'z_i',
430	& 'PBL height','m',
431	& dimout,zmax)
432	call WRITEDIAGFI(ngrid,'hf_max',
433	& 'Maximum vertical heat flux','m',
434	& dimout,hfmax)
435	call WRITEDIAGFI(ngrid,'vvv',
436	& 'Vertical velocity variance at zout','m',
437	& dimout,vvv)
438	call WRITEDIAGFI(ngrid,'vhf',
439	& 'Vertical heat flux at zout','m',
440	& dimout,vhf)
441
442	RETURN
443	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format