Context Navigation

pbl_parameters.F @ 2228

Last change on this file since 2228 was 1393, checked in by aslmd, 10 years ago
bug fix pbl_parameter in cases of specific Ri
File size: 12.3 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)
[1226]	4	USE comcstfi_h
[496]	5	IMPLICIT NONE
	6	!=======================================================================
	7	!
	8	! Anlysis of the PBL from input temperature, wind field and thermals outputs.
	9	!
	10	! -------
	11	!
	12	! Author: Arnaud Colaitis 09/01/12
	13	! -------
	14	!
	15	! Arguments:
	16	! ----------
	17	!
	18	! inputs:
	19	! ------
	20	! ngrid size of the horizontal grid
	21	! nlay size of the vertical grid
	22	! pz0(ngrid) surface roughness length
	23	! pg gravity (m s -2)
[529]	24	! zzlay(ngrid,nlay) height of mid-layers
	25	! zzlev(ngrid,nlay+1) height of layers interfaces
[496]	26	! pu(ngrid,nlay) u component of the wind
	27	! pv(ngrid,nlay) v component of the wind
[662]	28	! wstar_in(ngrid) free convection velocity in PBL
[496]	29	! hfmax(ngrid) maximum vertical turbulent heat flux in thermals
	30	! zmax(ngrid) height reached by the thermals (pbl height)
	31	! pts(ngrid) surface temperature
	32	! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa
[566]	33	! z_out(n_out) heights of interpolation
	34	! n_out number of points for interpolation
[496]	35	!
	36	! outputs:
	37	! ------
	38	!
[566]	39	! Teta_out(ngrid,n_out) interpolated teta
	40	! u_out(ngrid,n_out) interpolated u
[496]	41	! ustar(ngrid) friction velocity
	42	! tstar(ngrid) friction temperature
	43	! L_mo(ngrid) monin_obukhov length
	44	!
	45	!
	46	!=======================================================================
	47	!
	48	!-----------------------------------------------------------------------
	49	! Declarations:
	50	! -------------
	51
	52	#include "callkeys.h"
	53
	54	! Arguments:
	55	! ----------
	56
[566]	57	INTEGER, INTENT(IN) :: ngrid,nlay,n_out
[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)
[566]	63	REAL, INTENT(IN) :: z_out(n_out)
[529]	64
	65	! Outputs:
	66	! --------
	67
[605]	68	REAL, INTENT(OUT) :: T_out(ngrid,n_out),u_out(ngrid,n_out)
	69	REAL Teta_out(ngrid,n_out)
[529]	70	REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid)
[496]	71	REAL, INTENT(OUT) :: L_mo(ngrid)
	72
	73	! Local:
	74	! ------
	75
[566]	76	INTEGER ig,k,n
[496]	77	REAL karman,nu
	78	DATA karman,nu/.41,0.001/
	79	SAVE karman,nu
	80
	81	! Local(2):
	82	! ---------
	83
	84	REAL zout
	85	REAL rib(ngrid) ! Bulk Richardson number
	86	REAL fm(ngrid) ! stability function for momentum
	87	REAL fh(ngrid) ! stability function for heat
	88	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
	89	! phim = 1+betamzeta or (1-bmzeta)**am
	90	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
	91	REAL betam, betah, alphah, bm, bh, lambda
	92	! ah and am are assumed to be -0.25 and -0.5 respectively
	93	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient
	94	REAL pz0t ! initial thermal roughness length. (local)
	95	REAL ric ! critical richardson number
	96	REAL reynolds(ngrid) ! reynolds number for UBL
	97	REAL prandtl(ngrid) ! prandtl number for UBL
	98	REAL pz0tcomp(ngrid) ! computed z0t
	99	REAL ite
	100	REAL residual,zcd0,z1
	101	REAL pcdv(ngrid),pcdh(ngrid)
	102	REAL zu2(ngrid) ! Large-scale wind at first layer
	103	REAL pbl_teta(ngrid) ! mixed-layer potential temperature
	104	INTEGER pbl_height_index(ngrid) ! index of nearest vertical grid point for zmax
	105	REAL dteta(ngrid,nlay),x(ngrid) ! potential temperature gradient and z/zi
	106	REAL dvhf(ngrid),dvvv(ngrid) ! dimensionless vertical heat flux and
	107	! dimensionless vertical velocity variance
	108	REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux and vertical velocity variance
	109	INTEGER ii(1)
	110	! temporary
	111	INTEGER dimout
	112
	113	!------------------------------------------------------------------------
	114	!------------------------------------------------------------------------
	115	! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
	116	!------------------------------------------------------------------------
	117	!------------------------------------------------------------------------
	118
[566]	119	DO n=1,n_out
	120
[496]	121	c Initialisation :
	122
	123	L_mo(:)=0.
	124	ustar(:)=0.
	125	tstar(:)=0.
[566]	126	zout=z_out(n)
[496]	127	reynolds(:)=0.
	128	pz0t = 0.
	129	pz0tcomp(:) = 0.1*pz0(:)
	130	rib(:)=0.
	131	pcdv(:)=0.
	132	pcdh(:)=0.
	133
	134	! this formulation assumes alphah=1., implying betah=betam
	135	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
	136
	137	bm=16. !UBL
	138	bh=16. !UBL
	139	alphah=1.
	140	betam=5. !SBL
	141	betah=5. !SBL
	142	lambda=(sqrt(bh/bm))/alphah
	143	ric=betah/(betam**2)
	144	DO ig=1,ngrid
	145	ite=0.
	146	residual=abs(pz0tcomp(ig)-pz0t)
[1377]	147
[648]	148	zu2(ig)=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
[636]	149	& + (log(1.+0.7wstar_in(ig) + 2.3wstar_in(ig)2))2
[496]	150
	151	DO WHILE((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.))
	152
	153	pz0t=pz0tcomp(ig)
	154	IF (zu2(ig) .ne. 0.) THEN
	155	! Richardson number formulation proposed by D.E. England et al. (1995)
[605]	156	rib(ig) = (pg/pts(ig))
[636]	157	& sqrt(zzlay(ig,1)pz0(ig))
	158	& (((log(zzlay(ig,1)/pz0(ig)))*2)/(log(zzlay(ig,1)/pz0t)))
[1377]	159	& (ph(ig,1)-pts(ig))/(pu(ig,1)pu(ig,1)+pv(ig,1)*pv(ig,1))
[496]	160	ELSE
	161	print*,'warning, infinite Richardson at surface'
	162	print*,pu(ig,1),pv(ig,1)
	163	rib(ig) = ric
	164	ENDIF
	165
[636]	166	z1z0=zzlay(ig,1)/pz0(ig)
	167	z1z0t=zzlay(ig,1)/pz0t
[496]	168
	169	cdn(ig)=karman/log(z1z0)
	170	cdn(ig)=cdn(ig)*cdn(ig)
	171	chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
	172
	173	! STABLE BOUNDARY LAYER :
	174	IF (rib(ig) .gt. 0.) THEN
	175	! From D.E. England et al. (95)
	176	prandtl(ig)=1.
	177	if(rib(ig) .lt. ric) then
	178	! Assuming alphah=1. and bh=bm for stable conditions :
	179	fm(ig)=((ric-rib(ig))/ric)**2
	180	fh(ig)=((ric-rib(ig))/ric)**2
	181	else
	182	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
	183	fm(ig)=0.
	184	fh(ig)=0.
	185	endif
	186	! UNSTABLE BOUNDARY LAYER :
	187	ELSE
	188	! From D.E. England et al. (95)
	189	fm(ig)=sqrt(1.-lambdabmrib(ig))
	190	fh(ig)=(1./alphah)((1.-lambdabhrib(ig))0.5)
	191	& (1.-lambdabmrib(ig))**0.25
	192	prandtl(ig)=alphah((1.-lambdabmrib(ig))*0.25)/
	193	& ((1.-lambdabhrib(ig))**0.5)
	194	ENDIF
	195
	196	reynolds(ig)=karman*sqrt(fm(ig))
	197	& *sqrt(zu2(ig))
	198	& pz0(ig)/(log(z1z0)nu)
	199	pz0tcomp(ig)=pz0(ig)exp(-karman7.3*
	200	& (reynolds(ig)*0.25)(prandtl(ig)**0.5))
	201	residual = abs(pz0t-pz0tcomp(ig))
	202	ite = ite+1
[605]	203
[496]	204	ENDDO ! of while
	205	pz0t=0.
	206
	207	! Drag computation:
	208
	209	pcdv(ig)=cdn(ig)*fm(ig)
	210	pcdh(ig)=chn(ig)*fh(ig)
	211
	212	ENDDO ! of ngrid
	213
	214	!------------------------------------------------------------------------
	215	!------------------------------------------------------------------------
	216	! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
	217	!------------------------------------------------------------------------
	218	!------------------------------------------------------------------------
	219
[1393]	220	! u* theta* computation
[496]	221
	222	DO ig=1,ngrid
	223	IF (rib(ig) .ge. ric) THEN
	224	ustar(ig)=0.
	225	tstar(ig)=0.
	226	ELSE
	227	ustar(ig)=sqrt(pcdv(ig))
	228	& sqrt(pu(ig,1)pu(ig,1)+pv(ig,1)*pv(ig,1))
	229	tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
	230	& /sqrt(pcdv(ig))
	231	ENDIF
	232	ENDDO
	233
[1393]	234	! Interpolation:
	235
[496]	236	DO ig=1,ngrid
	237	IF(zout .lt. pz0tcomp(ig)) THEN
[566]	238	u_out(ig,n)=0.
	239	Teta_out(ig,n)=pts(ig)
[1393]	240
	241	ELSE
	242	IF (rib(ig) .ge. ric) THEN ! ustar=tstar=0 (and fm=fh=0)
	243	u_out(ig,n)=0
	244	Teta_out(ig,n)=pts(ig)
	245	ELSE
[1377]	246	u_out(ig,n)= ustar(ig)*log(zout/pz0(ig))/
	247	&(karman*sqrt(fm(ig)))
[636]	248
[1377]	249	Teta_out(ig,n)=pts(ig)+(tstar(ig)sqrt(fm(ig))log(zout/
	250	& (pz0tcomp(ig)))/
	251	&(karman*fh(ig)))
[1393]	252	ENDIF
[1377]	253	ENDIF
[636]	254
[1377]	255	IF (zout .lt. pz0(ig)) THEN
[566]	256	u_out(ig,n)=0.
[1377]	257	ENDIF
[636]	258
[496]	259	ENDDO
	260
	261	! when using convective adjustment without thermals, a vertical potential temperature
	262	! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
	263	! interpolated at any height in the surface layer is theta at the first level.
	264
	265	IF ((.not.calltherm).and.(calladj)) THEN
[605]	266	Teta_out(:,n)=ph(:,1)
	267	u_out(:,n)=(sqrt(cdn(:))sqrt(pu(:,1)pu(:,1)+pv(:,1)*pv(:,1))
	268	& /karman)*log(zout/pz0(:))
[496]	269	ENDIF
[566]	270	T_out(:,n) = Teta_out(:,n)*(exp(
[529]	271	& (zout/zzlay(:,1))*(log(pplay(:,1)/ps))
	272	& )
	273	& )**rcp
	274
[566]	275	ENDDO !of n=1,n_out
[529]	276
[1393]	277
[496]	278	!------------------------------------------------------------------------
	279	!------------------------------------------------------------------------
	280	! PART III : WSTAR COMPUTATION
	281	!------------------------------------------------------------------------
	282	!------------------------------------------------------------------------
	283
	284	! Detection of the mixed-layer potential temperature
	285	! ------------
	286
	287	! Nearest index for the pbl height
	288
[566]	289	IF (calltherm) THEN
	290
[496]	291	pbl_height_index(:)=1
	292
	293	DO k=1,nlay-1
	294	DO ig=1, ngrid
[529]	295	IF (abs(zmax(ig)-zzlay(ig,k)) .lt.
	296	& abs(zmax(ig)-zzlay(ig,pbl_height_index(ig)))) THEN
[496]	297	pbl_height_index(ig)=k
	298	ENDIF
	299	ENDDO
	300	ENDDO
	301
	302	! Potential temperature gradient
	303
	304	dteta(:,nlay)=0.
	305	DO k=1,nlay-1
	306	DO ig=1, ngrid
[529]	307	dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(zzlay(ig,k+1)-zzlay(ig,k))
[496]	308	ENDDO
	309	ENDDO
	310
	311	! Computation of the pbl mixed layer temperature
	312
	313	DO ig=1, ngrid
	314	ii=MINLOC(abs(dteta(ig,1:pbl_height_index(ig))))
	315	pbl_teta(ig) = ph(ig,ii(1))
	316	ENDDO
	317
[1393]	318
[496]	319	!------------------------------------------------------------------------
	320	!------------------------------------------------------------------------
	321	! PART IV : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
	322	!------------------------------------------------------------------------
	323	!------------------------------------------------------------------------
	324
	325	! We follow Spiga et. al 2010 (QJRMS)
	326	! ------------
	327
	328	DO ig=1, ngrid
	329	IF (zmax(ig) .gt. 0.) THEN
	330	x(ig) = zout/zmax(ig)
	331	ELSE
	332	x(ig) = 999.
	333	ENDIF
	334	ENDDO
	335
	336	DO ig=1, ngrid
	337	! dimensionless vertical heat flux
	338	IF (x(ig) .le. 0.3) THEN
	339	dvhf(ig) = ((-3.85/log(x(ig)))+0.07*log(x(ig)))
	340	& exp(-4.61x(ig))
	341	ELSEIF (x(ig) .le. 1.) THEN
	342	dvhf(ig) = -1.52*x(ig) + 1.24
	343	ELSE
	344	dvhf(ig) = 0.
	345	ENDIF
	346	! dimensionless vertical velocity variance
	347	IF (x(ig) .le. 1.) THEN
	348	dvvv(ig) = 2.05(x(ig)(2./3.))(1.-0.64x(ig))*2
	349	ELSE
	350	dvvv(ig) = 0.
	351	ENDIF
	352	ENDDO
	353
	354	vhf(:) = dvhf(:)*hfmax(:)
[662]	355	vvv(:) = dvvv(:)(wstar_in(:))*2
[496]	356
[566]	357	ENDIF ! of if calltherm
[496]	358
[1212]	359	#ifndef MESOSCALE
[636]	360	call WRITEDIAGFI(ngrid,'rib_pbl',
	361	& 'Richardson in pbl parameter','m/s',2,rib)
	362	call WRITEDIAGFI(ngrid,'cdn_pbl',
	363	& 'neutral momentum coef','m/s',2,cdn)
	364	call WRITEDIAGFI(ngrid,'fm_pbl',
	365	& 'momentum stab function','m/s',2,fm)
	366	call WRITEDIAGFI(ngrid,'uv',
	367	& 'wind norm first layer','m/s',2,sqrt(zu2))
	368	call WRITEDIAGFI(ngrid,'uvtrue',
	369	& 'wind norm first layer','m/s',2,sqrt(pu(:,1)2+pv(:,1)2))
	370	call WRITEDIAGFI(ngrid,'chn_pbl',
	371	& 'neutral momentum coef','m/s',2,chn)
	372	call WRITEDIAGFI(ngrid,'fh_pbl',
	373	& 'momentum stab function','m/s',2,fh)
	374	call WRITEDIAGFI(ngrid,'B_pbl',
[1377]	375	& 'buoyancy','m/',2,(ph(:,1)-pts(:))/pts(:))
[636]	376	call WRITEDIAGFI(ngrid,'zot_pbl',
[1377]	377	& 'buoyancy','ms',2,pz0tcomp)
	378	call WRITEDIAGFI(ngrid,'zz1','buoyancy','m',2,zzlay(:,1))
[1212]	379	#endif
[636]	380
[496]	381	RETURN
	382	END

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

Download in other formats:

Original Format