Context Navigation

pbl_parameters.F @ 508

Last change on this file since 508 was 496, checked in by acolaitis, 13 years ago

M 489 libf/phymars/thermcell_main_mars.F90
^{--------------------> New parameters for thermals, following latest version of the relevant article}

M 489 libf/phymars/physiq.F
^{--------------------> Minor modifications}

D 489 libf/phymars/surflayer_interpol.F
A 0 libf/phymars/pbl_parameters.F
^{--------------------> Replaced surflayer_interpol.F by a cleaner and richer version called pbl_parameters.F}

This routine is the base of what will later be implemented in the MCD and as a tool to
compute surface properties from output fields (so that it will ultimately disappear from libf
and might become an utility)

M 489 libf/phymars/vdif_cd.F
^{--------------------> Minor modification to the Richardson computation}

M 489 libf/phymars/vdifc.F
^{--------------------> Minor modification to the aerodynamic conductances computation, replacing wmax by hfmax, which}

is a better proxy for convective activity. Might be replaced by w_star later.

File size: 15.0 KB

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

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

Download in other formats:

Original Format