Context Navigation

pbl_parameters_mod.F90 @ 3212

Last change on this file since 3212 was 3167, checked in by llange, 19 months ago
Mars PCM Introducing the scheme from ATKE workshop to solve turbulent diffusion + surface layer parameterization. Works only if callatke = .true. in the callphys.def. By default, it is false and the model runs as usual with the yamada 2.5 closure Tuning of the several parameters for the ATKE in progress LL
File size: 17.4 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,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, 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	! pts(ngrid) surface temperature (K)
60	! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa (k)
61	! z_out(n_out) heights of interpolation (m)
62	! n_out number of points for interpolation (1)
63	!
64	! outputs:
65	! ------
66	!
67	! T_out(ngrid,n_out) interpolated potential temperature on z_out (K)
68	! u_out(ngrid,n_out) interpolated winds on z_out (m/s)
69	! ustar(ngrid) friction velocity (m/s)
70	! tstar(ngrid) friction temperature (K)
71	!
72	!
73	!=======================================================================
74	!
75	!-----------------------------------------------------------------------
76	! Declarations:
77	! -------------
78
79	#include "callkeys.h"
80
81	! Arguments:
82	! ----------
83
84	INTEGER, INTENT(IN) :: ngrid,nlay,n_out ! size of the horizontal and vertical grid, interpolated altitudes for the surface layer
85	REAL, INTENT(IN) :: pz0(ngrid) ! surface roughness
86	REAL, INTENT(IN) :: ps(ngrid),pplay(ngrid,nlay) ! surface pressure and pressure levels
87	REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay) ! surface gravity, altitude of the interface and mid-layers
88	REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) ! winds
89	REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid) ! free convection velocity , vertical turbulent heat, pbl height
90	REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) ! surface temperature, potentiel temperature
91	REAL, INTENT(IN) :: z_out(n_out) ! altitudes of the interpolation in the surface layer
92
93	! Outputs:
94	! --------
95
96	REAL, INTENT(OUT) :: T_out(ngrid,n_out),u_out(ngrid,n_out) ! Temperature and wind of the interpolation in the surface layer
97	REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! Friction velocity and temperature
98
99	! Local:
100	! ------
101
102	INTEGER ig,k,n ! Loop variables
103	INTEGER ii(1) ! Index to compute the pbl mixed layer temperature
104	REAL karman,nu ! Von Karman constant and fluid kinematic viscosity
105	DATA karman,nu/.41,0.001/
106
107	!$OMP THREADPRIVATE(karman,nu)
108
109	SAVE karman,nu
110
111	REAL zout ! altitude to interpolate (local variable during the loop on z_out) (m)
112	REAL rib(ngrid) ! Bulk Richardson number (1)
113	REAL fm(ngrid) ! stability function for momentum (1)
114	REAL fh(ngrid) ! stability function for heat (1)
115	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T (1)
116	! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where
117	! phim = 1+betamzeta or (1-bmzeta)**am
118	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
119	! ah and am are assumed to be -0.25 and -0.5 respectively
120	! lambda is an intermediate variable to simplify the expression
121	REAL betam, betah, alphah, bm, bh, lambda
122
123	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient (1)
124	REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m)
125	REAL ric ! critical richardson number (1)
126	REAL reynolds(ngrid) ! Reynolds number (1)
127	REAL prandtl(ngrid) ! Prandtl number (1)
128	REAL pz0tcomp(ngrid) ! Computed z0t during the iterative scheme (m)
129	REAL ite ! Number of iteration (1)
130	REAL residual ! Residual between pz0t and pz0tcomp (m)
131	REAL itemax ! Maximum number of iteration (1)
132	REAL tol_iter ! Tolerance for the residual to ensure convergence (1)
133	REAL pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient (1)
134	REAL zu2(ngrid) ! Near surface winds (m/s)
135	REAL x(ngrid) ! z/zi (1)
136	REAL dvhf(ngrid),dvvv(ngrid) ! dimensionless vertical heat flux and dimensionless vertical velocity variance
137	REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux (W/m^2) and vertical velocity variance (m/s)
138	REAL Teta_out(ngrid,n_out) ! Temporary variable to compute interpolated potential temperature (K)
139	REAL sm ! Stability function computed with the ATKE scheme
140	REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme
141
142	! Code:
143	! --------
144
145	!------------------------------------------------------------------------
146	!------------------------------------------------------------------------
147	! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
148	!------------------------------------------------------------------------
149	!------------------------------------------------------------------------
150
151
152
153	! Initialisation :
154
155	ustar(:)=0.
156	tstar(:)=0.
157	reynolds(:)=0.
158	pz0t = 0.
159	pz0tcomp(:) = 0.1*pz0(:)
160	rib(:)=0.
161	pcdv(:)=0.
162	pcdh(:)=0.
163	itemax= 10
164	tol_iter = 0.01
165	f_ri_cd_min = 0.01
166
167	! this formulation assumes alphah=1., implying betah=betam
168	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
169	bm = 16.
170	bh = 16.
171	alphah = 1.
172	betam = 5.
173	betah = 5.
174	lambda = (sqrt(bh/bm))/alphah
175
176	ric = betah/(betam**2)
177
178	DO ig=1,ngrid
179	ite = 0.
180	residual = abs(pz0tcomp(ig)-pz0t)
181	z1z0 = zzlay(ig,1)/pz0(ig)
182	cdn(ig) = karman/log(z1z0)
183	cdn(ig) = cdn(ig)*cdn(ig)
184	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
185	! IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.)
186	DO WHILE((residual .gt. tol_iter*pz0(ig)) .and. (ite .lt. itemax))
187	! Computations of z0T; iterated until z0T converges
188	pz0t = pz0tcomp(ig)
189	z1z0t=zzlay(ig,1)/pz0t
190	chn(ig) = cdn(ig)*log(z1z0)/log(z1z0t)
191	IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then
192	! Richardson number formulation proposed by D.E. England et al. (1995)
193	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)
194	ELSE
195	print*,'warning, infinite Richardson at surface'
196	print*,pu(ig,1),pv(ig,1)
197	rib(ig) = ric
198	ENDIF
199
200	! Compute the stability functions fm; fh depending on the stability of the surface layer
201
202	IF(callatke) THEN
203	! Computation following parameterizaiton from ATKE
204	IF(rib(ig) .gt. 0) THEN
205	! STABLE BOUNDARY LAYER :
206	sm = MAX(smmin,cn*(1.-rib(ig)/ric))
207	! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
208	prandtl(ig) = pr_neutexp(-pr_slope/pr_neutrib(ig)+rib(ig)/pr_neut) + rib(ig) * pr_slope
209	fm(ig) = max((sm*(3./2.) sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
210	fh(ig) = max((fm(ig) / prandtl(ig)), f_ri_cd_min)
211	ELSE
212	! UNSTABLE BOUNDARY LAYER :
213	sm = 2./rpi * (cinf - cn) * atan(-rib(ig)/ri0) + cn
214	prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig)/ri1) + pr_neut
215	fm(ig) = MAX((sm*(3./2.) sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
216	fh(ig) = MAX((fm(ig) / prandtl(ig)), f_ri_cd_min)
217	ENDIF ! Rib < 0
218	ELSE
219	! Computation following parameterizaiton from from D.E. England et al. (95)
220	IF (rib(ig) .gt. 0.) THEN
221	! STABLE BOUNDARY LAYER :
222	prandtl(ig) = 1.
223	IF(rib(ig) .lt. ric) THEN
224	! Assuming alphah=1. and bh=bm for stable conditions :
225	fm(ig) = ((ric-rib(ig))/ric)**2
226	fh(ig) = ((ric-rib(ig))/ric)**2
227	ELSE
228	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
229	fm(ig) = 1.
230	fh(ig) = 1.
231	ENDIF
232	ELSE
233	! UNSTABLE BOUNDARY LAYER :
234	fm(ig) = sqrt(1.-lambdabmrib(ig))
235	fh(ig) = (1./alphah)((1.-lambdabhrib(ig))0.5)(1.-lambdabmrib(ig))**0.25
236	prandtl(ig) = alphah((1.-lambdabmrib(ig))0.25)/((1.-lambdabhrib(ig))*0.5)
237	ENDIF ! Rib < 0
238	ENDIF ! atke
239	! Recompute the reynolds and z0t
240	reynolds(ig) = karmansqrt(fm(ig))sqrt(zu2(ig))pz0(ig)/(log(z1z0)nu)
241	pz0tcomp(ig) = pz0(ig)exp(-karman7.3(reynolds(ig)0.25)(prandtl(ig)*0.5)+5karman)
242	residual = abs(pz0t-pz0tcomp(ig))
243	ite = ite+1
244	ENDDO ! of while
245	! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions
246	pcdv(ig) = cdn(ig)*fm(ig)
247	pcdh(ig) = chn(ig)*fh(ig)
248	pz0t = 0. ! for next grid point
249	ENDDO ! of ngrid
250
251	!------------------------------------------------------------------------
252	!------------------------------------------------------------------------
253	! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
254	!------------------------------------------------------------------------
255	!------------------------------------------------------------------------
256
257	! u* theta* computation
258	DO n=1,n_out
259	zout = z_out(n)
260	DO ig=1,ngrid
261	IF (rib(ig) .ge. ric) THEN !stable case, no turbulence
262	ustar(ig) = 0.
263	tstar(ig) = 0.
264	ELSE
265	ustar(ig) = sqrt(pcdv(ig))sqrt(zu2(ig)) ! By definition, u = sqrt(tau/U^2) = sqrt(Cdv)
266	tstar(ig) = -pcdh(ig)(pts(ig)-ph(ig,1))/sqrt(pcdv(ig)) ! theta = (T1-Tsurf)Cdh(ig)/u
267	ENDIF
268
269	! Interpolation:
270
271	IF(zout .lt. pz0tcomp(ig)) THEN ! If z_out is in the surface layer , theta = tsurf; u = 0
272	u_out(ig,n) = 0.
273	Teta_out(ig,n) = pts(ig)
274	ELSEIF (zout .lt. pz0(ig)) THEN
275	u_out(ig,n) = 0.
276	ELSE
277	IF (rib(ig) .ge. ric) THEN ! ustar=tstar=0 (and fm=fh=0) because no turbulence
278	u_out(ig,n) = 0
279	Teta_out(ig,n) = pts(ig)
280	ELSE
281	! We use the MO theory: du/dz = u/(kappa z) 1 /sqrt(fm); dtheta/zd = theta/(Pr kappa z) fm/fh
282	u_out(ig,n)= ustar(ig)log(zout/pz0(ig))/(karmansqrt(fm(ig)))
283	Teta_out(ig,n)=pts(ig)+(tstar(ig)sqrt(fm(ig))log(zout/(pz0tcomp(ig)))/(karman*fh(ig)))
284	ENDIF
285	ENDIF
286	ENDDO ! loop on n_out
287
288	! when using convective adjustment without thermals, a vertical potential temperature
289	! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
290	! interpolated at any height in the surface layer is theta at the first level.
291
292	IF ((.not.calltherm).and.(calladj)) THEN
293	Teta_out(:,n)=ph(:,1)
294	u_out(:,n)=(sqrt(cdn(:))sqrt(zu2(ig))/karman)log(zout/pz0(:))
295	ENDIF
296	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 ?
297	ENDDO !of n=1,n_out
298
299
300	!------------------------------------------------------------------------
301	!------------------------------------------------------------------------
302	! PART III : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
303	!------------------------------------------------------------------------
304	!------------------------------------------------------------------------
305
306	! We follow Spiga et. al 2010 (QJRMS)
307	! ------------
308	IF (calltherm) THEN
309	DO ig=1, ngrid
310	IF (zmax(ig) .gt. 0.) THEN
311	x(ig) = zout/zmax(ig)
312	ELSE
313	x(ig) = 999.
314	ENDIF
315	ENDDO
316
317	DO ig=1, ngrid
318	! dimensionless vertical heat flux
319	IF (x(ig) .le. 0.3) THEN
320	dvhf(ig) = ((-3.85/log(x(ig)))+0.07log(x(ig))) exp(-4.61*x(ig))
321	ELSEIF (x(ig) .le. 1.) THEN
322	dvhf(ig) = -1.52*x(ig) + 1.24
323	ELSE
324	dvhf(ig) = 0.
325	ENDIF
326	! dimensionless vertical velocity variance
327	IF (x(ig) .le. 1.) THEN
328	dvvv(ig) = 2.05(x(ig)(2./3.))(1.-0.64x(ig))*2
329	ELSE
330	dvvv(ig) = 0.
331	ENDIF
332	ENDDO
333
334	vhf(:) = dvhf(:)*hfmax(:)
335	vvv(:) = dvvv(:)(wstar_in(:))*2
336
337	ENDIF ! of if calltherm
338
339
340	!------------------------------------------------------------------------
341	!------------------------------------------------------------------------
342	! PART IV : Outputs
343	!------------------------------------------------------------------------
344	!------------------------------------------------------------------------
345
346
347	#ifndef MESOSCALE
348	call write_output('rib_pbl','Richardson in pbl parameter','m/s',rib(:))
349	call write_output('cdn_pbl','neutral momentum coef','m/s',cdn(:))
350	call write_output('fm_pbl','momentum stab function','m/s',fm(:))
351	call write_output('uv','wind norm first layer','m/s',sqrt(zu2(:)))
352	call write_output('uvtrue','wind norm first layer','m/s',sqrt(pu(:,1)2+pv(:,1)2))
353	call write_output('chn_pbl','neutral momentum coef','m/s',chn(:))
354	call write_output('fh_pbl','momentum stab function','m/s',fh(:))
355	call write_output('B_pbl','buoyancy','m/',(ph(:,1)-pts(:))/pts(:))
356	call write_output('zot_pbl','buoyancy','ms',pz0tcomp(:))
357	call write_output('zz1','1st layer altitude','m',zzlay(:,1))
358	#endif
359
360	RETURN
361	END SUBROUTINE pbl_parameters
362
363	END MODULE pbl_parameters_mod
364

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

Download in other formats:

Original Format