Context Navigation

pbl_parameters.F @ 2740

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

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

Download in other formats:

Original Format