Context Navigation

pbl_parameters.F @ 545

Last change on this file since 545 was 529, checked in by acolaitis, 13 years ago
Physiq: changed called to pbl_parameter and vdifc. Vdifc: added a flag to deactivate mass-variation scheme for LES computations for two reason: it is costy for nothing and does not work in polar LES.
File size: 15.3 KB

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

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

Download in other formats:

Original Format