Context Navigation

pbl_parameters.F @ 2997

Last change on this file since 2997 was 2932, checked in by romain.vande, 21 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

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

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

Download in other formats:

Original Format