Context Navigation

yamada4.F @ 869

Last change on this file since 869 was 789, checked in by Laurent Fairhead, 17 years ago
Mauvais indice pour les boucles FH 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)
[766]	7	use dimphy
[541]	8	IMPLICIT NONE
	9	c.......................................................................
[766]	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
[766]	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/
[766]	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
[766]	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/
[766]	94	c$OMP THREADPRIVATE(ric,rifc,b1,kap)
[541]	95	real frif,falpha,fsm
	96	real fl,zzz,zl0,zq2,zn2
	97
[766]	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
[766]	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====================================================================
[766]	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
[789]	436	do ig=1,ngrid
[541]	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