Context Navigation

yamada4.F @ 1287

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

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

Download in other formats:

Original Format