Context Navigation

yamada4.F @ 746

Last change on this file since 746 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.9 KB

Rev	Line
[541]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
	5	s ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
	6	s ,iflag_pbl)
[704]	7	use dimphy
[541]	8	IMPLICIT NONE
	9	c.......................................................................
[704]	10	cym#include "dimensions.h"
	11	cym#include "dimphy.h"
[541]	12	c.......................................................................
	13	c
	14	c dt : pas de temps
	15	c g : g
	16	c zlev : altitude a chaque niveau (interface inferieure de la couche
	17	c de meme indice)
	18	c zlay : altitude au centre de chaque couche
	19	c u,v : vitesse au centre de chaque couche
	20	c (en entree : la valeur au debut du pas de temps)
	21	c teta : temperature potentielle au centre de chaque couche
	22	c (en entree : la valeur au debut du pas de temps)
	23	c cd : cdrag
	24	c (en entree : la valeur au debut du pas de temps)
	25	c q2 : $q^2$ au bas de chaque couche
	26	c (en entree : la valeur au debut du pas de temps)
	27	c (en sortie : la valeur a la fin du pas de temps)
	28	c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
	29	c couche)
	30	c (en sortie : la valeur a la fin du pas de temps)
	31	c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
	32	c (en sortie : la valeur a la fin du pas de temps)
	33	c
	34	c iflag_pbl doit valoir entre 6 et 9
	35	c l=6, on prend systematiquement une longueur d'equilibre
	36	c iflag_pbl=6 : MY 2.0
	37	c iflag_pbl=7 : MY 2.0.Fournier
	38	c iflag_pbl=8 : MY 2.5
	39	c iflag_pbl=9 : un test ?
	40
	41	c.......................................................................
	42	REAL dt,g,rconst
	43	real plev(klon,klev+1),temp(klon,klev)
	44	real ustar(klon)
	45	real kmin,qmin,pblhmin(klon),coriol(klon)
	46	REAL zlev(klon,klev+1)
	47	REAL zlay(klon,klev)
	48	REAL u(klon,klev)
	49	REAL v(klon,klev)
	50	REAL teta(klon,klev)
	51	REAL cd(klon)
	52	REAL q2(klon,klev+1),qpre
	53	REAL unsdz(klon,klev)
	54	REAL unsdzdec(klon,klev+1)
	55
	56	REAL km(klon,klev+1)
	57	REAL kmpre(klon,klev+1),tmp2
	58	REAL mpre(klon,klev+1)
	59	REAL kn(klon,klev+1)
	60	REAL kq(klon,klev+1)
	61	real ff(klon,klev+1),delta(klon,klev+1)
	62	real aa(klon,klev+1),aa0,aa1
	63	integer iflag_pbl,ngrid
	64
	65
	66	integer nlay,nlev
[704]	67	cym PARAMETER (nlay=klev)
	68	cym PARAMETER (nlev=klev+1)
[541]	69
	70	logical first
	71	integer ipas
	72	save first,ipas
	73	data first,ipas/.true.,0/
[704]	74	c$OMP THREADPRIVATE( first,ipas)
[541]	75
	76	integer ig,k
	77
	78
	79	real ri,zrif,zalpha,zsm,zsn
	80	real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
	81
	82	real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
	83	real dtetadz(klon,klev+1)
	84	real m2cstat,mcstat,kmcstat
[704]	85	real l(klon,klev+1)
	86	real,allocatable,save :: l0(:)
	87	c$OMP THREADPRIVATE(l0)
[541]	88	real sq(klon),sqz(klon),zz(klon,klev+1)
	89	integer iter
	90
	91	real ric,rifc,b1,kap
	92	save ric,rifc,b1,kap
	93	data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
[704]	94	c$OMP THREADPRIVATE(ric,rifc,b1,kap)
[541]	95	real frif,falpha,fsm
	96	real fl,zzz,zl0,zq2,zn2
	97
[704]	98	cym real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
	99	cym s ,lyam(klon,klev),knyam(klon,klev)
	100	cym s ,w2yam(klon,klev),t2yam(klon,klev)
	101	real,allocatable,save,dimension(:,:) :: rino,smyam,styam,lyam,
	102	s knyam,w2yam,t2yam
	103	cym common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam
	104	c$OMP THREADPRIVATE(rino,smyam,styam,lyam,knyam,w2yam,t2yam)
	105	logical,save :: firstcall=.true.
	106	c$OMP THREADPRIVATE(firstcall)
[541]	107	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
	108	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
	109	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
	110	fl(zzz,zl0,zq2,zn2)=
	111	s max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	112	s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
	113
[704]	114
	115	nlay=klev
	116	nlev=klev+1
	117
	118	if (firstcall) then
	119	allocate(rino(klon,klev+1),smyam(klon,klev),styam(klon,klev))
	120	allocate(lyam(klon,klev),knyam(klon,klev))
	121	allocate(w2yam(klon,klev),t2yam(klon,klev))
	122	allocate(l0(klon))
	123	firstcall=.false.
	124	endif
	125
	126
[541]	127	if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
	128	stop'probleme de coherence dans appel a MY'
	129	endif
	130
	131	ipas=ipas+1
	132	if (0.eq.1.and.first) then
	133	do ig=1,1000
	134	ri=(ig-800.)/500.
	135	if (ri.lt.ric) then
	136	zrif=frif(ri)
	137	else
	138	zrif=rifc
	139	endif
	140	if(zrif.lt.0.16) then
	141	zalpha=falpha(zrif)
	142	zsm=fsm(zrif)
	143	else
	144	zalpha=1.12
	145	zsm=0.085
	146	endif
	147	c print*,ri,rif,zalpha,zsm
	148	enddo
	149	endif
	150
	151	c.......................................................................
	152	c les increments verticaux
	153	c.......................................................................
	154	c
	155	c!!!!! allerte !!!!!c
	156	c!!!!! zlev n'est pas declare a nlev !!!!!c
	157	c!!!!! ---->
	158	DO ig=1,ngrid
	159	zlev(ig,nlev)=zlay(ig,nlay)
	160	& +( zlay(ig,nlay) - zlev(ig,nlev-1) )
	161	ENDDO
	162	c!!!!! <----
	163	c!!!!! allerte !!!!!c
	164	c
	165	DO k=1,nlay
	166	DO ig=1,ngrid
	167	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
	168	ENDDO
	169	ENDDO
	170	DO ig=1,ngrid
	171	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
	172	ENDDO
	173	DO k=2,nlay
	174	DO ig=1,ngrid
	175	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
	176	ENDDO
	177	ENDDO
	178	DO ig=1,ngrid
	179	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
	180	ENDDO
	181	c
	182	c.......................................................................
	183
	184	do k=2,klev
	185	do ig=1,ngrid
	186	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
	187	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2)
	188	s /(dz(ig,k)*dz(ig,k))
	189	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
	190	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
	191	c n2(ig,k)=0.
	192	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
	193	if (ri.lt.ric) then
	194	rif(ig,k)=frif(ri)
	195	else
	196	rif(ig,k)=rifc
	197	endif
	198	if(rif(ig,k).lt.0.16) then
	199	alpha(ig,k)=falpha(rif(ig,k))
	200	sm(ig,k)=fsm(rif(ig,k))
	201	else
	202	alpha(ig,k)=1.12
	203	sm(ig,k)=0.085
	204	endif
	205	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
	206	c print,'RIF L=',k,rif(ig,k),rialpha(ig,k)
	207
	208
	209	enddo
	210	enddo
	211
	212
	213	c====================================================================
	214	c Au premier appel, on determine l et q2 de facon iterative.
	215	c iterration pour determiner la longueur de melange
	216
	217
	218	if (first.or.iflag_pbl.eq.6) then
	219	do ig=1,ngrid
	220	l0(ig)=10.
	221	enddo
	222	do k=2,klev-1
	223	do ig=1,ngrid
	224	l(ig,k)=l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	225	enddo
	226	enddo
	227
	228	do iter=1,10
	229	do ig=1,ngrid
	230	sq(ig)=1.e-10
	231	sqz(ig)=1.e-10
	232	enddo
	233	do k=2,klev-1
	234	do ig=1,ngrid
	235	q2(ig,k)=l(ig,k)*2zz(ig,k)
	236	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	237	zq=sqrt(q2(ig,k))
	238	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	239	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	240	enddo
	241	enddo
	242	do ig=1,ngrid
	243	l0(ig)=0.2*sqz(ig)/sq(ig)
	244	c l0(ig)=30.
	245	enddo
	246	c print*,'ITER=',iter,' L0=',l0
	247
	248	enddo
	249
	250	c print*,'Fin de l initialisation de q2 et l0'
	251
	252	endif ! first
	253
	254	c====================================================================
	255	c Calcul de la longueur de melange.
	256	c====================================================================
	257
	258	c Mise a jour de l0
	259	do ig=1,ngrid
	260	sq(ig)=1.e-10
	261	sqz(ig)=1.e-10
	262	enddo
	263	do k=2,klev-1
	264	do ig=1,ngrid
	265	zq=sqrt(q2(ig,k))
	266	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	267	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	268	enddo
	269	enddo
	270	do ig=1,ngrid
	271	l0(ig)=0.2*sqz(ig)/sq(ig)
	272	c l0(ig)=30.
	273	enddo
	274	c print*,'ITER=',iter,' L0=',l0
	275	c calcul de l(z)
	276	do k=2,klev
	277	do ig=1,ngrid
	278	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	279	if(first) then
	280	q2(ig,k)=l(ig,k)*2zz(ig,k)
	281	endif
	282	enddo
	283	enddo
	284
	285	c====================================================================
	286	c Yamada 2.0
	287	c====================================================================
	288	if (iflag_pbl.eq.6) then
	289
	290	do k=2,klev
	291	do ig=1,ngrid
	292	q2(ig,k)=l(ig,k)*2zz(ig,k)
	293	enddo
	294	enddo
	295
	296
	297	else if (iflag_pbl.eq.7) then
	298	c====================================================================
	299	c Yamada 2.Fournier
	300	c====================================================================
	301
	302	c Calcul de l, km, au pas precedent
	303	do k=2,klev
	304	do ig=1,ngrid
	305	c print*,'SMML=',sm(ig,k),l(ig,k)
	306	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	307	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	308	mpre(ig,k)=sqrt(m2(ig,k))
	309	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	310	enddo
	311	enddo
	312
	313	do k=2,klev-1
	314	do ig=1,ngrid
	315	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
	316	mcstat=sqrt(m2cstat)
	317
	318	c print*,'M2 L=',k,mpre(ig,k),mcstat
	319	c
	320	c -----{puis on ecrit la valeur de q qui annule l'equation de m
	321	c supposee en q3}
	322	c
	323	IF (k.eq.2) THEN
	324	kmcstat=1.E+0 / mcstat
	325	& ( unsdz(ig,k)kmpre(ig,k+1)
	326	& *mpre(ig,k+1)
	327	& +unsdz(ig,k-1)
	328	& *cd(ig)
	329	& ( sqrt(u(ig,3)2+v(ig,3)*2)
	330	& -mcstat/unsdzdec(ig,k)
	331	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
	332	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	333	ELSE
	334	kmcstat=1.E+0 / mcstat
	335	& ( unsdz(ig,k)kmpre(ig,k+1)
	336	& *mpre(ig,k+1)
	337	& +unsdz(ig,k-1)*kmpre(ig,k-1)
	338	& *mpre(ig,k-1) )
	339	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	340	ENDIF
	341	c print*,'T2 L=',k,tmp2
	342	tmp2=kmcstat
	343	& /( sm(ig,k)/q2(ig,k) )
	344	& /l(ig,k)
	345	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
	346	c print*,'Q2 L=',k,q2(ig,k)
	347	c
	348	enddo
	349	enddo
	350
	351	else if (iflag_pbl.ge.8) then
	352	c====================================================================
	353	c Yamada 2.5 a la Didi
	354	c====================================================================
	355
	356
	357	c Calcul de l, km, au pas precedent
	358	do k=2,klev
	359	do ig=1,ngrid
	360	c print*,'SMML=',sm(ig,k),l(ig,k)
	361	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	362	if (delta(ig,k).lt.1.e-20) then
	363	c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
	364	delta(ig,k)=1.e-20
	365	endif
	366	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	367	aa0=
	368	s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
	369	aa1=
	370	s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
	371	c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
	372	aa(ig,k)=aa1dt/(delta(ig,k)l(ig,k))
	373	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	374	qpre=sqrt(q2(ig,k))
	375	if (iflag_pbl.eq.8 ) then
	376	if (aa(ig,k).gt.0.) then
	377	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
	378	else
	379	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	380	endif
	381	else ! iflag_pbl=9
	382	if (aa(ig,k)*qpre.gt.0.9) then
	383	q2(ig,k)=(qpre10.)*2
	384	else
	385	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	386	endif
	387	endif
	388	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
	389	c print,'Q2 L=',k,q2(ig,k),qpreqpre
	390	enddo
	391	enddo
	392
	393	endif ! Fin du cas 8
	394
	395	c print*,'OK8'
	396
	397	c====================================================================
[704]	398	c Calcul des coefficients de m�ange
[541]	399	c====================================================================
	400	do k=2,klev
	401	c print*,'k=',k
	402	do ig=1,ngrid
	403	cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
	404	zq=sqrt(q2(ig,k))
	405	km(ig,k)=l(ig,k)zqsm(ig,k)
	406	kn(ig,k)=km(ig,k)*alpha(ig,k)
	407	kq(ig,k)=l(ig,k)zq0.2
	408	c print*,'KML=',km(ig,k),kn(ig,k)
	409	enddo
	410	enddo
	411
	412	c if (iflag_pbl.ge.7..and.0.eq.1) then
	413	c q2(:,1)=q2(:,2)
	414	c call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
	415	c endif
	416
	417	c Traitement des cas noctrunes avec l'introduction d'une longueur
	418	c minilale.
	419
	420	c====================================================================
	421	c Traitement particulier pour les cas tres stables.
	422	c D'apres Holtslag Boville.
	423
	424	print*,'YAMADA4 0'
	425
	426	do ig=1,ngrid
	427	coriol(ig)=1.e-4
	428	pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
	429	enddo
	430
	431	print*,'pblhmin ',pblhmin
	432	CTest a remettre 21 11 02
	433	c test abd 13 05 02 if(0.eq.1) then
	434	if(1.eq.1) then
	435	do k=2,klev
	436	do ig=1,klon
	437	if (teta(ig,2).gt.teta(ig,1)) then
	438	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	439	kmin=kapzlev(ig,k)qmin
	440	else
	441	kmin=-1. ! kmin n'est utilise que pour les SL stables.
	442	endif
	443	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
	444	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
	445	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
	446	kn(ig,k)=kmin
	447	km(ig,k)=kmin
	448	kq(ig,k)=kmin
	449	c la longueur de melange est suposee etre l= kap z
	450	c K=l q Sm d'ou q2=(K/l Sm)**2
	451	q2(ig,k)=(qmin/sm(ig,k))**2
	452	endif
	453	enddo
	454	enddo
	455	endif
	456
	457	print*,'YAMADA4 1'
	458	c Diagnostique pour stokage
	459
	460	rino=rif
	461	smyam(:,1:klev)=sm(:,1:klev)
	462	styam=sm(:,1:klev)*alpha(:,1:klev)
	463	lyam(1:klon,1:klev)=l(:,1:klev)
	464	knyam(1:klon,1:klev)=kn(:,1:klev)
	465
	466	c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
	467
	468	if(1.eq.0)then
	469	w2yam=q2(:,1:klev)*0.24
	470	s +lyam(:,1:klev)5.17kn(:,1:klev)*n2(:,1:klev)
	471	s /sqrt(q2(:,1:klev))
	472
	473	t2yam=9.1kn(:,1:klev)dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
	474	s *lyam(:,1:klev)
	475	endif
	476
	477	c print*,'OKFIN'
	478	first=.false.
	479	return
	480	end

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

Download in other formats:

Original Format