Context Navigation

pbl_parameters_mod.F90 @ 3162

Last change on this file since 3162 was 3151, checked in by llange, 19 months ago
Mars PCM Switching pbl_parameter.F into a module. Cleaning of unused variables/sections in it. Correction in the computations of ustar, thetastar as it was done in a different way compare to vdif_cd. It is now self-consistent. Maybe pbl_surface in the MCD needs also to be changed ? In discussion with Ehouarn and Aymeric. LL
File size: 15.8 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
28	IMPLICIT NONE
29
30	!=======================================================================
31	!
32	! Anlysis of the PBL from input temperature, wind field and thermals outputs: compute the friction velocity, temperature, and monin_obukhov length
33	! and interpolate the potential temperature and winds in the surface layer using the Monin Obukhov theory
34	!
35	! -------
36	!
37	! Author: Arnaud Colaitis 09/01/12; adapted in F90 by Lucas Lange 12/2023
38	! -------
39	!
40	! Arguments:
41	! ----------
42	!
43	! inputs:
44	! ------
45	! ngrid size of the horizontal grid
46	! nlay size of the vertical grid
47	! ps(ngrid) surface pressure (Pa)
48	! pplay(ngrid,nlay) pressure levels
49	! pz0(ngrid) surface roughness length (m)
50	! pg gravity (m s -2)
51	! zzlay(ngrid,nlay) height of mid-layers (m)
52	! zzlev(ngrid,nlay+1) height of layers interfaces (m)
53	! pu(ngrid,nlay) u component of the wind (m/s)
54	! pv(ngrid,nlay) v component of the wind (m/s)
55	! wstar_in(ngrid) free convection velocity in PBL (m/s)
56	! hfmax(ngrid) maximum vertical turbulent heat flux in thermals (W/m^2)
57	! zmax(ngrid) height reached by the thermals (pbl height) (m)
58	! pts(ngrid) surface temperature (K)
59	! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa (k)
60	! z_out(n_out) heights of interpolation (m)
61	! n_out number of points for interpolation (1)
62	!
63	! outputs:
64	! ------
65	!
66	! T_out(ngrid,n_out) interpolated potential temperature on z_out (K)
67	! u_out(ngrid,n_out) interpolated winds on z_out (m/s)
68	! ustar(ngrid) friction velocity (m/s)
69	! tstar(ngrid) friction temperature (K)
70	!
71	!
72	!=======================================================================
73	!
74	!-----------------------------------------------------------------------
75	! Declarations:
76	! -------------
77
78	#include "callkeys.h"
79
80	! Arguments:
81	! ----------
82
83	INTEGER, INTENT(IN) :: ngrid,nlay,n_out ! size of the horizontal and vertical grid, interpolated altitudes for the surface layer
84	REAL, INTENT(IN) :: pz0(ngrid) ! surface roughness
85	REAL, INTENT(IN) :: ps(ngrid),pplay(ngrid,nlay) ! surface pressure and pressure levels
86	REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay) ! surface gravity, altitude of the interface and mid-layers
87	REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) ! winds
88	REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid) ! free convection velocity , vertical turbulent heat, pbl height
89	REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) ! surface temperature, potentiel temperature
90	REAL, INTENT(IN) :: z_out(n_out) ! altitudes of the interpolation in the surface layer
91
92	! Outputs:
93	! --------
94
95	REAL, INTENT(OUT) :: T_out(ngrid,n_out),u_out(ngrid,n_out) ! Temperature and wind of the interpolation in the surface layer
96	REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! Friction velocity and temperature
97
98	! Local:
99	! ------
100
101	INTEGER ig,k,n ! Loop variables
102	INTEGER ii(1) ! Index to compute the pbl mixed layer temperature
103	REAL karman,nu ! Von Karman constant and fluid kinematic viscosity
104	DATA karman,nu/.41,0.001/
105
106	!$OMP THREADPRIVATE(karman,nu)
107
108	SAVE karman,nu
109
110	REAL zout ! altitude to interpolate (local variable during the loop on z_out) (m)
111	REAL rib(ngrid) ! Bulk Richardson number (1)
112	REAL fm(ngrid) ! stability function for momentum (1)
113	REAL fh(ngrid) ! stability function for heat (1)
114	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T (1)
115	! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where
116	! phim = 1+betamzeta or (1-bmzeta)**am
117	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
118	! ah and am are assumed to be -0.25 and -0.5 respectively
119	! lambda is an intermediate variable to simplify the expression
120	REAL betam, betah, alphah, bm, bh, lambda
121
122	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient (1)
123	REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m)
124	REAL ric ! critical richardson number (1)
125	REAL reynolds(ngrid) ! Reynolds number (1)
126	REAL prandtl(ngrid) ! Prandtl number (1)
127	REAL pz0tcomp(ngrid) ! Computed z0t during the iterative scheme (m)
128	REAL ite ! Number of iteration (1)
129	REAL residual ! Residual between pz0t and pz0tcomp (m)
130	REAL itemax ! Maximum number of iteration (1)
131	REAL tol_iter ! Tolerance for the residual to ensure convergence (1)
132	REAL pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient (1)
133	REAL zu2(ngrid) ! Near surface winds (m/s)
134	REAL x(ngrid) ! z/zi (1)
135	REAL dvhf(ngrid),dvvv(ngrid) ! dimensionless vertical heat flux and dimensionless vertical velocity variance
136	REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux (W/m^2) and vertical velocity variance (m/s)
137	REAL Teta_out(ngrid,n_out) ! Temporary variable to compute interpolated potential temperature (K)
138
139	! Code:
140	! --------
141
142	!------------------------------------------------------------------------
143	!------------------------------------------------------------------------
144	! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
145	!------------------------------------------------------------------------
146	!------------------------------------------------------------------------
147
148
149
150	! Initialisation :
151
152	ustar(:)=0.
153	tstar(:)=0.
154	reynolds(:)=0.
155	pz0t = 0.
156	pz0tcomp(:) = 0.1*pz0(:)
157	rib(:)=0.
158	pcdv(:)=0.
159	pcdh(:)=0.
160	itemax= 10
161	tol_iter = 0.01
162
163	! this formulation assumes alphah=1., implying betah=betam
164	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
165	bm = 16.
166	bh = 16.
167	alphah = 1.
168	betam = 5.
169	betah = 5.
170	lambda = (sqrt(bh/bm))/alphah
171
172	ric = betah/(betam**2)
173
174	DO ig=1,ngrid
175	ite = 0.
176	residual = abs(pz0tcomp(ig)-pz0t)
177	z1z0 = zzlay(ig,1)/pz0(ig)
178	cdn(ig) = karman/log(z1z0)
179	cdn(ig) = cdn(ig)*cdn(ig)
180	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
181	! IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.)
182	DO WHILE((residual .gt. tol_iter*pz0(ig)) .and. (ite .lt. itemax))
183	! Computations of z0T; iterated until z0T converges
184	pz0t = pz0tcomp(ig)
185	z1z0t=zzlay(ig,1)/pz0t
186	chn(ig) = cdn(ig)*log(z1z0)/log(z1z0t)
187	IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then
188	! Richardson number formulation proposed by D.E. England et al. (1995)
189	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)
190	ELSE
191	print*,'warning, infinite Richardson at surface'
192	print*,pu(ig,1),pv(ig,1)
193	rib(ig) = ric
194	ENDIF
195
196	! Compute the stability functions fm; fh depending on the stability of the surface layer, from D.E. England et al. (95)
197	! STABLE BOUNDARY LAYER :
198	IF (rib(ig) .gt. 0.) THEN
199	prandtl(ig) = 1.
200	IF(rib(ig) .lt. ric) THEN
201	! Assuming alphah=1. and bh=bm for stable conditions :
202	fm(ig) = ((ric-rib(ig))/ric)**2
203	fh(ig) = ((ric-rib(ig))/ric)**2
204	ELSE
205	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
206	fm(ig) = 1.
207	fh(ig) = 1.
208	! should be 0 no?
209	ENDIF
210	! UNSTABLE BOUNDARY LAYER :
211	ELSE
212	fm(ig) = sqrt(1.-lambdabmrib(ig))
213	fh(ig) = (1./alphah)((1.-lambdabhrib(ig))0.5)(1.-lambdabmrib(ig))**0.25
214	prandtl(ig) = alphah((1.-lambdabmrib(ig))0.25)/((1.-lambdabhrib(ig))*0.5)
215	ENDIF
216	! Recompute the reynolds and z0t
217	reynolds(ig) = karmansqrt(fm(ig))sqrt(zu2(ig))pz0(ig)/(log(z1z0)nu)
218	pz0tcomp(ig) = pz0(ig)exp(-karman7.3(reynolds(ig)0.25)(prandtl(ig)*0.5)+5karman)
219	residual = abs(pz0t-pz0tcomp(ig))
220	ite = ite+1
221	ENDDO ! of while
222	! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions
223	pcdv(ig) = cdn(ig)*fm(ig)
224	pcdh(ig) = chn(ig)*fh(ig)
225	pz0t = 0. ! for next grid point
226	ENDDO ! of ngrid
227
228	!------------------------------------------------------------------------
229	!------------------------------------------------------------------------
230	! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
231	!------------------------------------------------------------------------
232	!------------------------------------------------------------------------
233
234	! u* theta* computation
235	DO n=1,n_out
236	zout = z_out(n)
237	DO ig=1,ngrid
238	IF (rib(ig) .ge. ric) THEN !stable case, no turbulence
239	ustar(ig) = 0.
240	tstar(ig) = 0.
241	ELSE
242	ustar(ig) = sqrt(pcdv(ig))sqrt(zu2(ig)) ! By definition, u = sqrt(tau/U^2) = sqrt(Cdv)
243	tstar(ig) = -pcdh(ig)(pts(ig)-ph(ig,1))/sqrt(pcdv(ig)) ! theta = (T1-Tsurf)Cdh(ig)/u
244	ENDIF
245
246	! Interpolation:
247
248	IF(zout .lt. pz0tcomp(ig)) THEN ! If z_out is in the surface layer , theta = tsurf; u = 0
249	u_out(ig,n) = 0.
250	Teta_out(ig,n) = pts(ig)
251	ELSEIF (zout .lt. pz0(ig)) THEN
252	u_out(ig,n) = 0.
253	ELSE
254	IF (rib(ig) .ge. ric) THEN ! ustar=tstar=0 (and fm=fh=0) because no turbulence
255	u_out(ig,n) = 0
256	Teta_out(ig,n) = pts(ig)
257	ELSE
258	! We use the MO theory: du/dz = u/(kappa z) 1 /sqrt(fm); dtheta/zd = theta/(Pr kappa z) fm/fh
259	u_out(ig,n)= ustar(ig)log(zout/pz0(ig))/(karmansqrt(fm(ig)))
260	Teta_out(ig,n)=pts(ig)+(tstar(ig)sqrt(fm(ig))log(zout/(pz0tcomp(ig)))/(karman*fh(ig)))
261	ENDIF
262	ENDIF
263	ENDDO ! loop on n_out
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(zu2(ig))/karman)log(zout/pz0(:))
272	ENDIF
273	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 ?
274	ENDDO !of n=1,n_out
275
276
277	!------------------------------------------------------------------------
278	!------------------------------------------------------------------------
279	! PART III : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
280	!------------------------------------------------------------------------
281	!------------------------------------------------------------------------
282
283	! We follow Spiga et. al 2010 (QJRMS)
284	! ------------
285	IF (calltherm) THEN
286	DO ig=1, ngrid
287	IF (zmax(ig) .gt. 0.) THEN
288	x(ig) = zout/zmax(ig)
289	ELSE
290	x(ig) = 999.
291	ENDIF
292	ENDDO
293
294	DO ig=1, ngrid
295	! dimensionless vertical heat flux
296	IF (x(ig) .le. 0.3) THEN
297	dvhf(ig) = ((-3.85/log(x(ig)))+0.07log(x(ig))) exp(-4.61*x(ig))
298	ELSEIF (x(ig) .le. 1.) THEN
299	dvhf(ig) = -1.52*x(ig) + 1.24
300	ELSE
301	dvhf(ig) = 0.
302	ENDIF
303	! dimensionless vertical velocity variance
304	IF (x(ig) .le. 1.) THEN
305	dvvv(ig) = 2.05(x(ig)(2./3.))(1.-0.64x(ig))*2
306	ELSE
307	dvvv(ig) = 0.
308	ENDIF
309	ENDDO
310
311	vhf(:) = dvhf(:)*hfmax(:)
312	vvv(:) = dvvv(:)(wstar_in(:))*2
313
314	ENDIF ! of if calltherm
315
316
317	!------------------------------------------------------------------------
318	!------------------------------------------------------------------------
319	! PART IV : Outputs
320	!------------------------------------------------------------------------
321	!------------------------------------------------------------------------
322
323
324	#ifndef MESOSCALE
325	call write_output('rib_pbl','Richardson in pbl parameter','m/s',rib(:))
326	call write_output('cdn_pbl','neutral momentum coef','m/s',cdn(:))
327	call write_output('fm_pbl','momentum stab function','m/s',fm(:))
328	call write_output('uv','wind norm first layer','m/s',sqrt(zu2(:)))
329	call write_output('uvtrue','wind norm first layer','m/s',sqrt(pu(:,1)2+pv(:,1)2))
330	call write_output('chn_pbl','neutral momentum coef','m/s',chn(:))
331	call write_output('fh_pbl','momentum stab function','m/s',fh(:))
332	call write_output('B_pbl','buoyancy','m/',(ph(:,1)-pts(:))/pts(:))
333	call write_output('zot_pbl','buoyancy','ms',pz0tcomp(:))
334	call write_output('zz1','1st layer altitude','m',zzlay(:,1))
335	#endif
336
337	RETURN
338	END SUBROUTINE pbl_parameters
339
340	END MODULE pbl_parameters_mod
341

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

Download in other formats:

Original Format