Context Navigation

pbl_parameters_mod.F90 @ 3225

Last change on this file since 3225 was 3219, checked in by llange, 12 months ago
Mars PCM Small correction for the surface layer scheme: the value of the critical Richardson varys if one used the classic yamada scheme (0.2 in this case) or can be a tunning parameter if one uses the atke scheme Also add the possibily to write the tke in the diagfi.nc when calling pbl_parameter LL
File size: 18.0 KB

Line
1	MODULE pbl_parameters_mod
2
3	!======================================================================================================================!
4	! Subject:
5	!---------
6	! Module used to compute the friction velocity, temperature, and monin_obukhov length from temperature, wind field and thermals outputs.
7	! Also interpolate theta and u;v in the surface layer based on the Monin Obukhov theory
8	! These are diagnostics only and do not influence the code.
9	!----------------------------------------------------------------------------------------------------------------------!
10	! Reference:
11	!-----------
12	! Colaïtis, A., Spiga, A., Hourdin, F., Rio, C., Forget, F., and Millour, E. (2013), A thermal plume model for the Martian convective boundary layer, J. Geophys. Res. Planets, 118, 1468–1487, doi:10.1002/jgre.20104.
13	!
14	!======================================================================================================================!
15
16
17	IMPLICIT NONE
18
19	CONTAINS
20
21	SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0,pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,tke,pts,ph,z_out,n_out, &
22	T_out,u_out,ustar,tstar,vhf,vvv)
23
24	USE comcstfi_h
25	use write_output_mod, only: write_output
26	use turb_mod, only: turb_resolved
27	use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi
28
29	IMPLICIT NONE
30
31	!=======================================================================
32	!
33	! Anlysis of the PBL from input temperature, wind field and thermals outputs: compute the friction velocity, temperature, and monin_obukhov length
34	! and interpolate the potential temperature and winds in the surface layer using the Monin Obukhov theory
35	!
36	! -------
37	!
38	! Author: Arnaud Colaitis 09/01/12; adapted in F90 by Lucas Lange 12/2023
39	! -------
40	!
41	! Arguments:
42	! ----------
43	!
44	! inputs:
45	! ------
46	! ngrid size of the horizontal grid
47	! nlay size of the vertical grid
48	! ps(ngrid) surface pressure (Pa)
49	! pplay(ngrid,nlay) pressure levels
50	! pz0(ngrid) surface roughness length (m)
51	! pg gravity (m s -2)
52	! zzlay(ngrid,nlay) height of mid-layers (m)
53	! zzlev(ngrid,nlay+1) height of layers interfaces (m)
54	! pu(ngrid,nlay) u component of the wind (m/s)
55	! pv(ngrid,nlay) v component of the wind (m/s)
56	! wstar_in(ngrid) free convection velocity in PBL (m/s)
57	! hfmax(ngrid) maximum vertical turbulent heat flux in thermals (W/m^2)
58	! zmax(ngrid) height reached by the thermals (pbl height) (m)
59	! tke(ngrid,nlay+1) turbulent kinetic energy (J/kg)
60	! pts(ngrid) surface temperature (K)
61	! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa (k)
62	! z_out(n_out) heights of interpolation (m)
63	! n_out number of points for interpolation (1)
64	!
65	! outputs:
66	! ------
67	!
68	! T_out(ngrid,n_out) interpolated potential temperature on z_out (K)
69	! u_out(ngrid,n_out) interpolated winds on z_out (m/s)
70	! ustar(ngrid) friction velocity (m/s)
71	! tstar(ngrid) friction temperature (K)
72	!
73	!
74	!=======================================================================
75	!
76	!-----------------------------------------------------------------------
77	! Declarations:
78	! -------------
79
80	#include "callkeys.h"
81
82	! Arguments:
83	! ----------
84
85	INTEGER, INTENT(IN) :: ngrid,nlay,n_out ! size of the horizontal and vertical grid, interpolated altitudes for the surface layer
86	REAL, INTENT(IN) :: pz0(ngrid) ! surface roughness
87	REAL, INTENT(IN) :: ps(ngrid),pplay(ngrid,nlay) ! surface pressure and pressure levels
88	REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay) ! surface gravity, altitude of the interface and mid-layers
89	REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) ! winds
90	REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid) ! free convection velocity , vertical turbulent heat, pbl height
91	REAL, INTENT(IN) :: tke(ngrid,nlay+1) ! Turbulent kinetic energy
92	REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) ! surface temperature, potentiel temperature
93	REAL, INTENT(IN) :: z_out(n_out) ! altitudes of the interpolation in the surface layer
94
95	! Outputs:
96	! --------
97
98	REAL, INTENT(OUT) :: T_out(ngrid,n_out),u_out(ngrid,n_out) ! Temperature and wind of the interpolation in the surface layer
99	REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! Friction velocity and temperature
100
101	! Local:
102	! ------
103
104	INTEGER ig,k,n ! Loop variables
105	INTEGER ii(1) ! Index to compute the pbl mixed layer temperature
106	REAL karman,nu ! Von Karman constant and fluid kinematic viscosity
107	DATA karman,nu/.41,0.001/
108
109	!$OMP THREADPRIVATE(karman,nu)
110
111	SAVE karman,nu
112
113	REAL zout ! altitude to interpolate (local variable during the loop on z_out) (m)
114	REAL rib(ngrid) ! Bulk Richardson number (1)
115	REAL fm(ngrid) ! stability function for momentum (1)
116	REAL fh(ngrid) ! stability function for heat (1)
117	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T (1)
118	! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where
119	! phim = 1+betamzeta or (1-bmzeta)**am
120	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
121	! ah and am are assumed to be -0.25 and -0.5 respectively
122	! lambda is an intermediate variable to simplify the expression
123	REAL betam, betah, alphah, bm, bh, lambda
124
125	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient (1)
126	REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m)
127	REAL ric_colaitis ! critical richardson number (1)
128	REAL reynolds(ngrid) ! Reynolds number (1)
129	REAL prandtl(ngrid) ! Prandtl number (1)
130	REAL pz0tcomp(ngrid) ! Computed z0t during the iterative scheme (m)
131	REAL ite ! Number of iteration (1)
132	REAL residual ! Residual between pz0t and pz0tcomp (m)
133	REAL itemax ! Maximum number of iteration (1)
134	REAL tol_iter ! Tolerance for the residual to ensure convergence (1)
135	REAL pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient (1)
136	REAL zu2(ngrid) ! Near surface winds (m/s)
137	REAL x(ngrid) ! z/zi (1)
138	REAL dvhf(ngrid),dvvv(ngrid) ! dimensionless vertical heat flux and dimensionless vertical velocity variance
139	REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux (W/m^2) and vertical velocity variance (m/s)
140	REAL Teta_out(ngrid,n_out) ! Temporary variable to compute interpolated potential temperature (K)
141	REAL sm ! Stability function computed with the ATKE scheme
142	REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme
143	REAL ric_4interp ! critical richardson number which is either the one from Colaitis et al. (2013) or ATKE (1)
144	! Code:
145	! --------
146
147	!------------------------------------------------------------------------
148	!------------------------------------------------------------------------
149	! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
150	!------------------------------------------------------------------------
151	!------------------------------------------------------------------------
152
153
154
155	! Initialisation :
156
157	ustar(:)=0.
158	tstar(:)=0.
159	reynolds(:)=0.
160	pz0t = 0.
161	pz0tcomp(:) = 0.1*pz0(:)
162	rib(:)=0.
163	pcdv(:)=0.
164	pcdh(:)=0.
165	itemax= 10
166	tol_iter = 0.01
167	f_ri_cd_min = 0.01
168
169	! this formulation assumes alphah=1., implying betah=betam
170	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
171	bm = 16.
172	bh = 16.
173	alphah = 1.
174	betam = 5.
175	betah = 5.
176	lambda = (sqrt(bh/bm))/alphah
177	ric_colaitis = betah/(betam**2)
178
179	if(callatke) then
180	ric_4interp = ric
181	else
182	ric_4interp = ric_colaitis
183	endif
184
185
186	DO ig=1,ngrid
187	ite = 0.
188	residual = abs(pz0tcomp(ig)-pz0t)
189	z1z0 = zzlay(ig,1)/pz0(ig)
190	cdn(ig) = karman/log(z1z0)
191	cdn(ig) = cdn(ig)*cdn(ig)
192	zu2(ig) = pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1) + (log(1.+0.7wstar_in(ig) + 2.3wstar_in(ig)2))2 ! near surface winds, accounting for buyoncy
193	! IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.)
194	DO WHILE((residual .gt. tol_iter*pz0(ig)) .and. (ite .lt. itemax))
195	! Computations of z0T; iterated until z0T converges
196	pz0t = pz0tcomp(ig)
197	z1z0t=zzlay(ig,1)/pz0t
198	chn(ig) = cdn(ig)*log(z1z0)/log(z1z0t)
199	IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then
200	! Richardson number formulation proposed by D.E. England et al. (1995)
201	rib(ig) = (pg/pts(ig))sqrt(zzlay(ig,1)pz0(ig))(((log(zzlay(ig,1)/pz0(ig)))2)/(log(zzlay(ig,1)/pz0t)))(ph(ig,1)-pts(ig))/zu2(ig)
202	ELSE
203	print*,'warning, infinite Richardson at surface'
204	print*,pu(ig,1),pv(ig,1)
205	rib(ig) = ric_colaitis
206	ENDIF
207
208	! Compute the stability functions fm; fh depending on the stability of the surface layer
209
210	IF(callatke) THEN
211	! Computation following parameterizaiton from ATKE
212	IF(rib(ig) .gt. 0) THEN
213	! STABLE BOUNDARY LAYER :
214	sm = MAX(smmin,cn*(1.-rib(ig)/ric))
215	! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
216	prandtl(ig) = pr_neutexp(-pr_slope/pr_neutrib(ig)+rib(ig)/pr_neut) + rib(ig) * pr_slope
217	fm(ig) = max((sm*(3./2.) sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
218	fh(ig) = max((fm(ig) / prandtl(ig)), f_ri_cd_min)
219	ELSE
220	! UNSTABLE BOUNDARY LAYER :
221	sm = 2./rpi * (cinf - cn) * atan(-rib(ig)/ri0) + cn
222	prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig)/ri1) + pr_neut
223	fm(ig) = MAX((sm*(3./2.) sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
224	fh(ig) = MAX((fm(ig) / prandtl(ig)), f_ri_cd_min)
225	ENDIF ! Rib < 0
226	ELSE
227	! Computation following parameterizaiton from from D.E. England et al. (95)
228	IF (rib(ig) .gt. 0.) THEN
229	! STABLE BOUNDARY LAYER :
230	prandtl(ig) = 1.
231	IF(rib(ig) .lt. ric_colaitis) THEN
232	! Assuming alphah=1. and bh=bm for stable conditions :
233	fm(ig) = ((ric_colaitis-rib(ig))/ric_colaitis)**2
234	fh(ig) = ((ric_colaitis-rib(ig))/ric_colaitis)**2
235	ELSE
236	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
237	fm(ig) = 1.
238	fh(ig) = 1.
239	ENDIF
240	ELSE
241	! UNSTABLE BOUNDARY LAYER :
242	fm(ig) = sqrt(1.-lambdabmrib(ig))
243	fh(ig) = (1./alphah)((1.-lambdabhrib(ig))0.5)(1.-lambdabmrib(ig))**0.25
244	prandtl(ig) = alphah((1.-lambdabmrib(ig))0.25)/((1.-lambdabhrib(ig))*0.5)
245	ENDIF ! Rib < 0
246	ENDIF ! atke
247	! Recompute the reynolds and z0t
248	reynolds(ig) = karmansqrt(fm(ig))sqrt(zu2(ig))pz0(ig)/(log(z1z0)nu)
249	pz0tcomp(ig) = pz0(ig)exp(-karman7.3(reynolds(ig)0.25)(prandtl(ig)*0.5)+5karman)
250	residual = abs(pz0t-pz0tcomp(ig))
251	ite = ite+1
252	ENDDO ! of while
253	! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions
254	pcdv(ig) = cdn(ig)*fm(ig)
255	pcdh(ig) = chn(ig)*fh(ig)
256	pz0t = 0. ! for next grid point
257	ENDDO ! of ngrid
258
259	!------------------------------------------------------------------------
260	!------------------------------------------------------------------------
261	! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
262	!------------------------------------------------------------------------
263	!------------------------------------------------------------------------
264
265	! u* theta* computation
266	DO n=1,n_out
267	zout = z_out(n)
268	DO ig=1,ngrid
269	IF (rib(ig) .ge. ric_4interp) THEN !stable case, no turbulence
270	ustar(ig) = 0.
271	tstar(ig) = 0.
272	ELSE
273	ustar(ig) = sqrt(pcdv(ig))sqrt(zu2(ig)) ! By definition, u = sqrt(tau/U^2) = sqrt(Cdv)
274	tstar(ig) = -pcdh(ig)(pts(ig)-ph(ig,1))/sqrt(pcdv(ig)) ! theta = (T1-Tsurf)Cdh(ig)/u
275	ENDIF
276
277	! Interpolation:
278
279	IF(zout .lt. pz0tcomp(ig)) THEN ! If z_out is in the surface layer , theta = tsurf; u = 0
280	u_out(ig,n) = 0.
281	Teta_out(ig,n) = pts(ig)
282	ELSEIF (zout .lt. pz0(ig)) THEN
283	u_out(ig,n) = 0.
284	ELSE
285	IF (rib(ig) .ge. ric_4interp) THEN ! ustar=tstar=0 (and fm=fh=0) because no turbulence
286	u_out(ig,n) = 0
287	Teta_out(ig,n) = pts(ig)
288	ELSE
289	! We use the MO theory: du/dz = u/(kappa z) 1 /sqrt(fm); dtheta/zd = theta/(Pr kappa z) fm/fh
290	u_out(ig,n)= ustar(ig)log(zout/pz0(ig))/(karmansqrt(fm(ig)))
291	Teta_out(ig,n)=pts(ig)+(tstar(ig)sqrt(fm(ig))log(zout/(pz0tcomp(ig)))/(karman*fh(ig)))
292	ENDIF
293	ENDIF
294	ENDDO ! loop on n_out
295
296	! when using convective adjustment without thermals, a vertical potential temperature
297	! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
298	! interpolated at any height in the surface layer is theta at the first level.
299
300	IF ((.not.calltherm).and.(calladj)) THEN
301	Teta_out(:,n)=ph(:,1)
302	u_out(:,n)=(sqrt(cdn(:))sqrt(zu2(ig))/karman)log(zout/pz0(:))
303	ENDIF
304	T_out(:,n) = Teta_out(:,n)(exp((zout/zzlay(:,1))(log(pplay(:,1)/ps))))**rcp ! not sure of what is done here: hydrostatic correction to account for the difference of pressure between zout and z1 ?
305	ENDDO !of n=1,n_out
306
307
308	!------------------------------------------------------------------------
309	!------------------------------------------------------------------------
310	! PART III : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
311	!------------------------------------------------------------------------
312	!------------------------------------------------------------------------
313
314	! We follow Spiga et. al 2010 (QJRMS)
315	! ------------
316	IF (calltherm) THEN
317	DO ig=1, ngrid
318	IF (zmax(ig) .gt. 0.) THEN
319	x(ig) = zout/zmax(ig)
320	ELSE
321	x(ig) = 999.
322	ENDIF
323	ENDDO
324
325	DO ig=1, ngrid
326	! dimensionless vertical heat flux
327	IF (x(ig) .le. 0.3) THEN
328	dvhf(ig) = ((-3.85/log(x(ig)))+0.07log(x(ig))) exp(-4.61*x(ig))
329	ELSEIF (x(ig) .le. 1.) THEN
330	dvhf(ig) = -1.52*x(ig) + 1.24
331	ELSE
332	dvhf(ig) = 0.
333	ENDIF
334	! dimensionless vertical velocity variance
335	IF (x(ig) .le. 1.) THEN
336	dvvv(ig) = 2.05(x(ig)(2./3.))(1.-0.64x(ig))*2
337	ELSE
338	dvvv(ig) = 0.
339	ENDIF
340	ENDDO
341
342	vhf(:) = dvhf(:)*hfmax(:)
343	vvv(:) = dvvv(:)(wstar_in(:))*2
344
345	ENDIF ! of if calltherm
346
347
348	!------------------------------------------------------------------------
349	!------------------------------------------------------------------------
350	! PART IV : Outputs
351	!------------------------------------------------------------------------
352	!------------------------------------------------------------------------
353
354
355	#ifndef MESOSCALE
356	call write_output('tke_pbl','Tke in the pbl after physic','J/kg',tke(:,:))
357	call write_output('rib_pbl','Richardson in pbl parameter','m/s',rib(:))
358	call write_output('cdn_pbl','neutral momentum coef','m/s',cdn(:))
359	call write_output('fm_pbl','momentum stab function','m/s',fm(:))
360	call write_output('uv','wind norm first layer','m/s',sqrt(zu2(:)))
361	call write_output('uvtrue','wind norm first layer','m/s',sqrt(pu(:,1)2+pv(:,1)2))
362	call write_output('chn_pbl','neutral momentum coef','m/s',chn(:))
363	call write_output('fh_pbl','momentum stab function','m/s',fh(:))
364	call write_output('B_pbl','buoyancy','m/',(ph(:,1)-pts(:))/pts(:))
365	call write_output('zot_pbl','buoyancy','ms',pz0tcomp(:))
366	call write_output('zz1','1st layer altitude','m',zzlay(:,1))
367	#endif
368
369	RETURN
370	END SUBROUTINE pbl_parameters
371
372	END MODULE pbl_parameters_mod
373

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

Download in other formats:

Original Format