Context Navigation

pbl_parameters.F @ 627

Last change on this file since 627 was 605, checked in by acolaitis, 13 years ago
Added a number of corrections to pbl_parameter. The routine can now be called either without thermals (it will then use log law for wind and first level temp for t_out) and with thermals and richardson surface layer (it will conduct a Monin Obukhov interpolation). WARNING: The routine default output is now T instead of TETA. (not much change at this altitude though)
File size: 13.9 KB

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

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

Download in other formats:

Original Format