Context Navigation

surflayer_interpol.F @ 324

Last change on this file since 324 was 319, checked in by aslmd, 14 years ago
LMDZ.MARS : minor changes to thermals [output L_mo + default settings]. transparent to users.
File size: 10.8 KB

Rev	Line
[267]	1	SUBROUTINE surflayer_interpol(ngrid,nlay,pz0,
[319]	2	& pg,pz,pu,pv,wmax,pts,ph,z_out,Teta_out,u_out,ustar,tstar,L_mo)
[267]	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*
[319]	47	REAL, INTENT(OUT) :: L_mo(ngrid) ! Monin-Obukhov length
[267]	48
	49	! Local:
	50	! ------
	51
	52	INTEGER ig
	53
	54	REAL karman,nu
	55	DATA karman,nu/.41,0.001/
	56	SAVE karman,nu
	57
	58	! Local(2):
	59	! ---------
[272]	60	REAL zout
[267]	61
	62	REAL rib(ngrid) ! Bulk Richardson number
	63	REAL fm(ngrid) ! stability function for momentum
	64	REAL fh(ngrid) ! stability function for heat
	65	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
	66
	67	! phim = 1+betamzeta or (1-bmzeta)**am
	68	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
	69	REAL betam, betah, alphah, bm, bh, lambda
	70	! ah and am are assumed to be -0.25 and -0.5 respectively
	71
	72	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient
	73	REAL pz0t ! initial thermal roughness length. (local)
	74	REAL ric ! critical richardson number
	75	REAL reynolds(ngrid) ! reynolds number for UBL
	76	REAL prandtl(ngrid) ! prandtl number for UBL
	77	REAL pz0tcomp(ngrid) ! computed z0t
	78	REAL ite
	79	REAL residual
	80	REAL pcdv(ngrid),pcdh(ngrid)
	81	! For output :
	82
[268]	83	REAL zu2(ngrid) ! Large-scale wind at first layer
[267]	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))
[292]	136	& ( pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1))
[268]	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)))
[292]	146	& (ph(ig,1)-pts(ig))/(zu2(ig) + (0.5wmax(ig))**2)
[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))
[292]	188	& sqrt(zu2(ig)+(0.5wmax(ig))**2)
[272]	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
[292]	228	ustar(ig)=sqrt(pcdv(ig))*sqrt(zu2(ig))
	229	tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
	230	& /sqrt(pcdv(ig))
[268]	231
	232	if (tstar(ig) .lt. -50) then
	233	print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig))
	234	& ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig))
	235	endif
[267]	236	endif
	237	ENDDO
	238
	239	! Monin Obukhov length :
	240
	241	DO ig=1,ngrid
	242	if (rib(ig) .gt. ric) then
	243	L_mo(ig)=0.
	244	else
	245	L_mo(ig)=pts(ig)(ustar(ig)2)/(pgkarman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL
	246	endif
	247	ENDDO
	248
[268]	249	DO ig=1,ngrid
[272]	250	IF(zout .ge. pz(ig,1)) THEN
	251	zout=1.
[267]	252	print*, 'Warning, z_out should be less than the first
[272]	253	& vertical grid-point'
[267]	254	print*, 'z1 =',pz(ig,1)
	255	print*, 'z_out=',z_out
	256	print*, 'z_out has been changed to 1m
[272]	257	& and computation is resumed'
[267]	258	ENDIF
	259
[272]	260	IF(zout .lt. pz0tcomp(ig)) THEN
	261	zout=pz0tcomp(ig)
	262	print*, 'Warning, z_out should be more than the thermal
	263	& roughness length'
	264	print*, 'z0 =',pz0tcomp(ig)
[267]	265	print*, 'z_out=',z_out
[272]	266	print*, 'z_out has been changed to z0t
	267	& and computation is resumed'
[267]	268	ENDIF
[268]	269	ENDDO
[267]	270
[272]	271	print*, 'interpolation of u and teta at z_out=',zout
[267]	272
	273	DO ig=1,ngrid
	274	IF (L_mo(ig) .gt. 0.) THEN
[272]	275	u_out(ig)=(ustar(ig)/karman)*log(zout/pz0(ig)) +
	276	& 5.(ustar(ig)/(karmanL_mo(ig)))*(zout-pz0(ig))
[267]	277	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))
[272]	278	& *log(zout/pz0tcomp(ig)) +
[267]	279	& 5.(tstar(ig)/(prandtl(ig)karman*L_mo(ig)))
[272]	280	& *(zout-pz0tcomp(ig))
[267]	281	ELSEIF (L_mo(ig) .lt. 0.) THEN
[268]	282
	283	IF(L_mo(ig) .gt. -1000.) THEN
	284
	285	u_out(ig)=(ustar(ig)/karman)*(
[272]	286	& 2.atan((1.-16.zout/L_mo(ig))**0.25)
[268]	287	& -2.atan((1.-16.pz0(ig)/L_mo(ig))**0.25)
[272]	288	& -2.log(1.+(1.-16.zout/L_mo(ig))**0.25)
[268]	289	& +2.log(1.+(1.-16.pz0(ig)/L_mo(ig))**0.25)
[272]	290	& - log(1.+sqrt(1.-16.*zout/L_mo(ig)))
[268]	291	& + log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
[272]	292	& + log(zout/pz0(ig))
[267]	293	& )
[268]	294
	295	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
	296	& 2.log(1.+sqrt(1.-16.pz0tcomp(ig)/L_mo(ig)))
[272]	297	& -2.log(1.+sqrt(1.-16.zout/L_mo(ig)))
	298	& + log(zout/pz0tcomp(ig))
[267]	299	& )
[268]	300
	301	ELSE
	302
	303	! We have to treat the case where L is very large and negative,
	304	! corresponding to a very nearly stable atmosphere (but not quite !)
	305	! i.e. teta* <0 and almost zero. We choose the cutoff
	306	! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
	307	! between the two expression for z/L = -1/1000 is -1.54324*10^-9
	308	! (we do that to avoid using r*4 precision, otherwise, we get -inf values)
	309
	310	u_out(ig)=(ustar(ig)/karman)*(
[272]	311	& (4./L_mo(ig))*(zout-pz0(ig))
	312	& + (20./(L_mo(ig))*2)(zout2-pz0(ig)2)
	313	& + (160./(L_mo(ig))*3)(zout3-pz0(ig)3)
	314	& + log(zout/pz0(ig))
[268]	315	& )
	316
	317	Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
[272]	318	& (8./L_mo(ig))*(zout-pz0tcomp(ig))
	319	& + (48./(L_mo(ig))*2)(zout2-pz0tcomp(ig)2)
	320	& + (1280./(3.(L_mo(ig))3))(zout3-pz0tcomp(ig)3)
	321	& + log(zout/pz0tcomp(ig))
[268]	322	& )
	323
	324	ENDIF
[267]	325	ELSE
	326	u_out(ig)=0.
[268]	327	Teta_out(ig)=pts(ig)
[267]	328	ENDIF
	329	ENDDO
	330
[268]	331	! Usefull diagnostics for the interpolation routine :
	332
[272]	333	! call WRITEDIAGFI(ngrid,'z0T',
	334	! & 'thermal roughness length','m',
	335	! & 2,pz0tcomp)
	336	! call WRITEDIAGFI(ngrid,'z0',
	337	! & 'roughness length','m',
	338	! & 2,pz0)
[267]	339
	340
	341	RETURN
	342	END

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

Download in other formats:

Original Format