Context Navigation

yamada4.F @ 684

Last change on this file since 684 was 541, checked in by lmdzadmin, 20 years ago

Convergence avec la version d'Olivia Coindreau incluant:

le offline
les thermiques
mellor & yamada dans la couche limite

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format