Context Navigation

surflayer_interpol.F @ 289

Last change on this file since 289 was 272, checked in by acolaitis, 14 years ago
Modifications du surflayer_interpol.F : needed local duplicate for z_out, and included subgrid gustiness in friction quantities estimation.
File size: 10.9 KB

Rev	Line
[267]	1	SUBROUTINE surflayer_interpol(ngrid,nlay,pz0,
	2	& pg,pz,pu,pv,wmax,pts,ph,z_out,Teta_out,u_out,ustar,tstar)
	3	IMPLICIT NONE
	4	!=======================================================================
	5	!
	6	! Subject: interpolation of temperature and velocity norm in the surface layer
[268]	7	! by recomputation of surface layer quantities (Richardson, Prandtl, Reynolds, u, teta)
[267]	8	! -------
	9	!
	10	! Author: Arnaud Colaitis 05/08/11
	11	! -------
	12	!
	13	! Arguments:
	14	! ----------
	15	!
	16	! inputs:
	17	! ------
	18	! ngrid size of the horizontal grid
	19	! pg gravity (m s -2)
	20	! pz(ngrid,nlay) height of layers
	21	! pu(ngrid,nlay) u component of the wind
	22	! pv(ngrid,nlay) v component of the wind
	23	! pts(ngrid) surface temperature
	24	! ph(ngrid) potential temperature T*(p/ps)^kappa
	25	!
	26	!
	27	!=======================================================================
	28	!
	29	!-----------------------------------------------------------------------
	30	! Declarations:
	31	! -------------
	32
	33	#include "comcstfi.h"
	34
	35	! Arguments:
	36	! ----------
	37
	38	INTEGER, INTENT(IN) :: ngrid,nlay
	39	REAL, INTENT(IN) :: pz0(ngrid)
	40	REAL, INTENT(IN) :: pg,pz(ngrid,nlay)
	41	REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
	42	REAL, INTENT(IN) :: wmax(ngrid)
	43	REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
[272]	44	REAL, INTENT(IN) :: z_out ! output height (in m above surface)
[268]	45	REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)! interpolated fields at z_out : potential temperature and norm(uv)
[267]	46	REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! u* and teta*
	47
	48	! Local:
	49	! ------
	50
	51	INTEGER ig
	52
	53	REAL karman,nu
	54	DATA karman,nu/.41,0.001/
	55	SAVE karman,nu
	56
	57	! Local(2):
	58	! ---------
[272]	59	REAL zout
[267]	60
	61	REAL rib(ngrid) ! Bulk Richardson number
	62	REAL fm(ngrid) ! stability function for momentum
	63	REAL fh(ngrid) ! stability function for heat
	64	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
	65
	66	! phim = 1+betamzeta or (1-bmzeta)**am
	67	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
	68	REAL betam, betah, alphah, bm, bh, lambda
	69	! ah and am are assumed to be -0.25 and -0.5 respectively
	70
	71	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient
	72	REAL pz0t ! initial thermal roughness length. (local)
	73	REAL ric ! critical richardson number
	74	REAL reynolds(ngrid) ! reynolds number for UBL
	75	REAL prandtl(ngrid) ! prandtl number for UBL
	76	REAL pz0tcomp(ngrid) ! computed z0t
	77	REAL ite
	78	REAL residual
	79	REAL pcdv(ngrid),pcdh(ngrid)
	80	! For output :
	81
[268]	82	REAL zu2(ngrid) ! Large-scale wind at first layer
[267]	83	REAL L_mo(ngrid) ! Monin-Obukhov length
	84	!-----------------------------------------------------------------------
	85	! couche de surface:
	86	! ------------------
[272]	87	zout=z_out
[267]	88	tstar(:)=0.
	89	ustar(:)=0.
	90	reynolds(:)=0.
	91
	92	! New formulation (AC) :
	93
	94	! phim = 1+betamzeta or (1-bmzeta)**am
	95	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
	96	! am=-0.25, ah=-0.5
	97
	98	pz0t = 0. ! for the sake of simplicity
	99	pz0tcomp(:) = 0.1*pz0(:)
	100	rib(:)=0.
	101	pcdv(:)=0.
	102	pcdh(:)=0.
	103
	104	! this formulation assumes alphah=1., implying betah=betam
	105	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
	106	bm=16. !UBL
	107	bh=16. !UBL
	108	alphah=1.
	109	betam=5. !SBL
	110	betah=5. !SBL
	111	lambda=(sqrt(bh/bm))/alphah
	112	ric=betah/(betam**2)
	113
	114	DO ig=1,ngrid
	115
	116	ite=0.
	117	residual=abs(pz0tcomp(ig)-pz0t)
	118
	119	do while((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.))
	120
	121	pz0t=pz0tcomp(ig)
	122
	123	if ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then
	124
	125	! Classical Richardson number formulation
	126
	127	! rib(ig) = (pg/ph(ig,1))pz(ig,1)(ph(ig,1)-pts(ig))
	128	! & /(pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1))
	129
	130	! Richardson number formulation proposed by D.E. England et al. (1995)
	131
[268]	132	! IF((pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1) .lt. 1.)
	133	! & .and. (wmax(ig) .gt. 0.)) THEN
	134	zu2(ig)=
	135	! & (MAX(pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1),wmax(ig)**2))
	136	& ( pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1) + wmax(ig)**2)
	137	! & pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1)
	138	! ELSE
	139	! zu2(ig)=pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1)
	140	! ENDIF
	141
[267]	142	rib(ig) = (pg/ph(ig,1))
[268]	143	! & pz(ig,1)pz0(ig)/sqrt(pz(ig,1)*pz0t)
	144	& sqrt(pz(ig,1)pz0(ig))
[267]	145	& (((log(pz(ig,1)/pz0(ig)))*2)/(log(pz(ig,1)/pz0t)))
[268]	146	& *(ph(ig,1)-pts(ig))/zu2(ig)
[267]	147
[268]	148	! & /(MAX(pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1),wmax(ig)**2))
	149	! & /( pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1) + wmax(ig)**2)
	150
[267]	151	else
	152	print*,'warning, infinite Richardson at surface'
	153	print*,pu(ig,1),pv(ig,1)
	154	rib(ig) = ric ! traiter ce cas !
	155	endif
	156
	157	z1z0=pz(ig,1)/pz0(ig)
	158	z1z0t=pz(ig,1)/pz0t
	159
	160	cdn(ig)=karman/log(z1z0)
	161	cdn(ig)=cdn(ig)*cdn(ig)
	162	chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
	163
	164	! Stable case :
	165	if (rib(ig) .gt. 0.) then
	166	! From D.E. England et al. (95)
	167	prandtl(ig)=1.
	168	if(rib(ig) .lt. ric) then
	169	! Assuming alphah=1. and bh=bm for stable conditions :
	170	fm(ig)=((ric-rib(ig))/ric)**2
	171	fh(ig)=((ric-rib(ig))/ric)**2
	172	else
	173	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
	174	fm(ig)=0.
	175	fh(ig)=0.
	176	endif
	177	! Unstable case :
	178	else
	179	! From D.E. England et al. (95)
	180	fm(ig)=sqrt(1.-lambdabmrib(ig))
	181	fh(ig)=(1./alphah)((1.-lambdabhrib(ig))0.5)
	182	& (1.-lambdabmrib(ig))**0.25
	183	prandtl(ig)=alphah((1.-lambdabmrib(ig))*0.25)/
	184	& ((1.-lambdabhrib(ig))**0.5)
	185	endif
	186
[272]	187	reynolds(ig)=karman*sqrt(fm(ig))
	188	& *sqrt(zu2(ig))
	189	! & sqrt(pu(ig,1)2 + pv(ig,1)*2)
[267]	190	& pz0(ig)/(log(z1z0)nu)
	191	pz0tcomp(ig)=pz0(ig)exp(-karman7.3*
	192	& (reynolds(ig)*0.25)(prandtl(ig)**0.5))
	193
	194
	195	residual = abs(pz0t-pz0tcomp(ig))
	196	ite = ite+1
	197	! if(ite .eq. 10) then
	198	! print*, 'iteration max reached'
	199	! endif
	200	! print*, "iteration nnumber, residual",ite,residual
	201
	202	enddo ! of while
	203
	204	pz0t=0.
	205
	206	! Drag computation :
	207
	208	pcdv(ig)=cdn(ig)*fm(ig)
	209	pcdh(ig)=chn(ig)*fh(ig)
	210
	211	ENDDO
	212
[268]	213	! Large-scale wind at first layer (without gustiness) and
[267]	214	! u* theta* computation
	215	DO ig=1,ngrid
	216
	217	if (rib(ig) .gt. ric) then
	218	ustar(ig)=0.
	219	tstar(ig)=0.
	220	else
[268]	221
	222	! ustar(ig)=karmansqrt(fm(ig)zu2(ig))/(log(pz(ig,1)/pz0(ig)))
	223	! tstar(ig)=karmanfh(ig)(ph(ig,1)-pts(ig))
	224	! & /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig)))
	225
	226	!simpler definition of u* and teta*, equivalent to the formula above :
	227
	228	ustar(ig)=sqrt(pcdv(ig))*sqrt(zu2(ig))
	229	tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))/sqrt(pcdv(ig))
	230
	231	if (tstar(ig) .lt. -50) then
	232	print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig))
	233	& ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig))
	234	endif
[267]	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
[268]	248	DO ig=1,ngrid
[272]	249	IF(zout .ge. pz(ig,1)) THEN
	250	zout=1.
[267]	251	print*, 'Warning, z_out should be less than the first
[272]	252	& vertical grid-point'
[267]	253	print*, 'z1 =',pz(ig,1)
	254	print*, 'z_out=',z_out
	255	print*, 'z_out has been changed to 1m
[272]	256	& and computation is resumed'
[267]	257	ENDIF
	258
[272]	259	IF(zout .lt. pz0tcomp(ig)) THEN
	260	zout=pz0tcomp(ig)
	261	print*, 'Warning, z_out should be more than the thermal
	262	& roughness length'
	263	print*, 'z0 =',pz0tcomp(ig)
[267]	264	print*, 'z_out=',z_out
[272]	265	print*, 'z_out has been changed to z0t
	266	& and computation is resumed'
[267]	267	ENDIF
[268]	268	ENDDO
[267]	269
[272]	270	print*, 'interpolation of u and teta at z_out=',zout
[267]	271
	272	DO ig=1,ngrid
	273	IF (L_mo(ig) .gt. 0.) THEN
[272]	274	u_out(ig)=(ustar(ig)/karman)*log(zout/pz0(ig)) +
	275	& 5.(ustar(ig)/(karmanL_mo(ig)))*(zout-pz0(ig))
[267]	276	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))
[272]	277	& *log(zout/pz0tcomp(ig)) +
[267]	278	& 5.(tstar(ig)/(prandtl(ig)karman*L_mo(ig)))
[272]	279	& *(zout-pz0tcomp(ig))
[267]	280	ELSEIF (L_mo(ig) .lt. 0.) THEN
[268]	281
	282	IF(L_mo(ig) .gt. -1000.) THEN
	283
	284	u_out(ig)=(ustar(ig)/karman)*(
[272]	285	& 2.atan((1.-16.zout/L_mo(ig))**0.25)
[268]	286	& -2.atan((1.-16.pz0(ig)/L_mo(ig))**0.25)
[272]	287	& -2.log(1.+(1.-16.zout/L_mo(ig))**0.25)
[268]	288	& +2.log(1.+(1.-16.pz0(ig)/L_mo(ig))**0.25)
[272]	289	& - log(1.+sqrt(1.-16.*zout/L_mo(ig)))
[268]	290	& + log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
[272]	291	& + log(zout/pz0(ig))
[267]	292	& )
[268]	293
	294	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
	295	& 2.log(1.+sqrt(1.-16.pz0tcomp(ig)/L_mo(ig)))
[272]	296	& -2.log(1.+sqrt(1.-16.zout/L_mo(ig)))
	297	& + log(zout/pz0tcomp(ig))
[267]	298	& )
[268]	299
	300	ELSE
	301
	302	! We have to treat the case where L is very large and negative,
	303	! corresponding to a very nearly stable atmosphere (but not quite !)
	304	! i.e. teta* <0 and almost zero. We choose the cutoff
	305	! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
	306	! between the two expression for z/L = -1/1000 is -1.54324*10^-9
	307	! (we do that to avoid using r*4 precision, otherwise, we get -inf values)
	308
	309	u_out(ig)=(ustar(ig)/karman)*(
[272]	310	& (4./L_mo(ig))*(zout-pz0(ig))
	311	& + (20./(L_mo(ig))*2)(zout2-pz0(ig)2)
	312	& + (160./(L_mo(ig))*3)(zout3-pz0(ig)3)
	313	& + log(zout/pz0(ig))
[268]	314	& )
	315
	316	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
[272]	317	& (8./L_mo(ig))*(zout-pz0tcomp(ig))
	318	& + (48./(L_mo(ig))*2)(zout2-pz0tcomp(ig)2)
	319	& + (1280./(3.(L_mo(ig))3))(zout3-pz0tcomp(ig)3)
	320	& + log(zout/pz0tcomp(ig))
[268]	321	& )
	322
	323	ENDIF
[267]	324	ELSE
	325	u_out(ig)=0.
[268]	326	Teta_out(ig)=pts(ig)
[267]	327	ENDIF
	328	ENDDO
	329
[268]	330	! Usefull diagnostics for the interpolation routine :
	331
[272]	332	! call WRITEDIAGFI(ngrid,'L',
	333	! & 'Monin Obukhov length','m',
	334	! & 2,L_mo)
	335	! call WRITEDIAGFI(ngrid,'z0T',
	336	! & 'thermal roughness length','m',
	337	! & 2,pz0tcomp)
	338	! call WRITEDIAGFI(ngrid,'z0',
	339	! & 'roughness length','m',
	340	! & 2,pz0)
[267]	341
	342
	343	RETURN
	344	END

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

Download in other formats:

Original Format