Context Navigation

pbl_parameters.F @ 3026

Last change on this file since 3026 was 2932, checked in by romain.vande, 21 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

File size: 12.4 KB

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

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

Download in other formats:

Original Format