Context Navigation

yamada4.F @ 1566

Last change on this file since 1566 was 634, checked in by Laurent Fairhead, 20 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.8 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)
[634]	7	use dimphy
[541]	8	IMPLICIT NONE
	9	c.......................................................................
[634]	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
[634]	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/
	74
	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
[634]	85	real l(klon,klev+1)
	86	real,allocatable,save :: l0(:)
	87
[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/
	94
	95	real frif,falpha,fsm
	96	real fl,zzz,zl0,zq2,zn2
	97
[634]	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	logical,save :: firstcall=.true.
	105
[541]	106	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
	107	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
	108	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
	109	fl(zzz,zl0,zq2,zn2)=
	110	s max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	111	s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
	112
[634]	113
	114	nlay=klev
	115	nlev=klev+1
	116
	117	if (firstcall) then
	118	allocate(rino(klon,klev+1),smyam(klon,klev),styam(klon,klev))
	119	allocate(lyam(klon,klev),knyam(klon,klev))
	120	allocate(w2yam(klon,klev),t2yam(klon,klev))
	121	allocate(l0(klon))
	122	firstcall=.false.
	123	endif
	124
	125
[541]	126	if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
	127	stop'probleme de coherence dans appel a MY'
	128	endif
	129
	130	ipas=ipas+1
	131	if (0.eq.1.and.first) then
	132	do ig=1,1000
	133	ri=(ig-800.)/500.
	134	if (ri.lt.ric) then
	135	zrif=frif(ri)
	136	else
	137	zrif=rifc
	138	endif
	139	if(zrif.lt.0.16) then
	140	zalpha=falpha(zrif)
	141	zsm=fsm(zrif)
	142	else
	143	zalpha=1.12
	144	zsm=0.085
	145	endif
	146	c print*,ri,rif,zalpha,zsm
	147	enddo
	148	endif
	149
	150	c.......................................................................
	151	c les increments verticaux
	152	c.......................................................................
	153	c
	154	c!!!!! allerte !!!!!c
	155	c!!!!! zlev n'est pas declare a nlev !!!!!c
	156	c!!!!! ---->
	157	DO ig=1,ngrid
	158	zlev(ig,nlev)=zlay(ig,nlay)
	159	& +( zlay(ig,nlay) - zlev(ig,nlev-1) )
	160	ENDDO
	161	c!!!!! <----
	162	c!!!!! allerte !!!!!c
	163	c
	164	DO k=1,nlay
	165	DO ig=1,ngrid
	166	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
	167	ENDDO
	168	ENDDO
	169	DO ig=1,ngrid
	170	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
	171	ENDDO
	172	DO k=2,nlay
	173	DO ig=1,ngrid
	174	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
	175	ENDDO
	176	ENDDO
	177	DO ig=1,ngrid
	178	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
	179	ENDDO
	180	c
	181	c.......................................................................
	182
	183	do k=2,klev
	184	do ig=1,ngrid
	185	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
	186	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2)
	187	s /(dz(ig,k)*dz(ig,k))
	188	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
	189	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
	190	c n2(ig,k)=0.
	191	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
	192	if (ri.lt.ric) then
	193	rif(ig,k)=frif(ri)
	194	else
	195	rif(ig,k)=rifc
	196	endif
	197	if(rif(ig,k).lt.0.16) then
	198	alpha(ig,k)=falpha(rif(ig,k))
	199	sm(ig,k)=fsm(rif(ig,k))
	200	else
	201	alpha(ig,k)=1.12
	202	sm(ig,k)=0.085
	203	endif
	204	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
	205	c print,'RIF L=',k,rif(ig,k),rialpha(ig,k)
	206
	207
	208	enddo
	209	enddo
	210
	211
	212	c====================================================================
	213	c Au premier appel, on determine l et q2 de facon iterative.
	214	c iterration pour determiner la longueur de melange
	215
	216
	217	if (first.or.iflag_pbl.eq.6) then
	218	do ig=1,ngrid
	219	l0(ig)=10.
	220	enddo
	221	do k=2,klev-1
	222	do ig=1,ngrid
	223	l(ig,k)=l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	224	enddo
	225	enddo
	226
	227	do iter=1,10
	228	do ig=1,ngrid
	229	sq(ig)=1.e-10
	230	sqz(ig)=1.e-10
	231	enddo
	232	do k=2,klev-1
	233	do ig=1,ngrid
	234	q2(ig,k)=l(ig,k)*2zz(ig,k)
	235	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	236	zq=sqrt(q2(ig,k))
	237	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	238	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	239	enddo
	240	enddo
	241	do ig=1,ngrid
	242	l0(ig)=0.2*sqz(ig)/sq(ig)
	243	c l0(ig)=30.
	244	enddo
	245	c print*,'ITER=',iter,' L0=',l0
	246
	247	enddo
	248
	249	c print*,'Fin de l initialisation de q2 et l0'
	250
	251	endif ! first
	252
	253	c====================================================================
	254	c Calcul de la longueur de melange.
	255	c====================================================================
	256
	257	c Mise a jour de l0
	258	do ig=1,ngrid
	259	sq(ig)=1.e-10
	260	sqz(ig)=1.e-10
	261	enddo
	262	do k=2,klev-1
	263	do ig=1,ngrid
	264	zq=sqrt(q2(ig,k))
	265	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	266	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	267	enddo
	268	enddo
	269	do ig=1,ngrid
	270	l0(ig)=0.2*sqz(ig)/sq(ig)
	271	c l0(ig)=30.
	272	enddo
	273	c print*,'ITER=',iter,' L0=',l0
	274	c calcul de l(z)
	275	do k=2,klev
	276	do ig=1,ngrid
	277	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	278	if(first) then
	279	q2(ig,k)=l(ig,k)*2zz(ig,k)
	280	endif
	281	enddo
	282	enddo
	283
	284	c====================================================================
	285	c Yamada 2.0
	286	c====================================================================
	287	if (iflag_pbl.eq.6) then
	288
	289	do k=2,klev
	290	do ig=1,ngrid
	291	q2(ig,k)=l(ig,k)*2zz(ig,k)
	292	enddo
	293	enddo
	294
	295
	296	else if (iflag_pbl.eq.7) then
	297	c====================================================================
	298	c Yamada 2.Fournier
	299	c====================================================================
	300
	301	c Calcul de l, km, au pas precedent
	302	do k=2,klev
	303	do ig=1,ngrid
	304	c print*,'SMML=',sm(ig,k),l(ig,k)
	305	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	306	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	307	mpre(ig,k)=sqrt(m2(ig,k))
	308	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	309	enddo
	310	enddo
	311
	312	do k=2,klev-1
	313	do ig=1,ngrid
	314	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
	315	mcstat=sqrt(m2cstat)
	316
	317	c print*,'M2 L=',k,mpre(ig,k),mcstat
	318	c
	319	c -----{puis on ecrit la valeur de q qui annule l'equation de m
	320	c supposee en q3}
	321	c
	322	IF (k.eq.2) THEN
	323	kmcstat=1.E+0 / mcstat
	324	& ( unsdz(ig,k)kmpre(ig,k+1)
	325	& *mpre(ig,k+1)
	326	& +unsdz(ig,k-1)
	327	& *cd(ig)
	328	& ( sqrt(u(ig,3)2+v(ig,3)*2)
	329	& -mcstat/unsdzdec(ig,k)
	330	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
	331	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	332	ELSE
	333	kmcstat=1.E+0 / mcstat
	334	& ( unsdz(ig,k)kmpre(ig,k+1)
	335	& *mpre(ig,k+1)
	336	& +unsdz(ig,k-1)*kmpre(ig,k-1)
	337	& *mpre(ig,k-1) )
	338	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	339	ENDIF
	340	c print*,'T2 L=',k,tmp2
	341	tmp2=kmcstat
	342	& /( sm(ig,k)/q2(ig,k) )
	343	& /l(ig,k)
	344	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
	345	c print*,'Q2 L=',k,q2(ig,k)
	346	c
	347	enddo
	348	enddo
	349
	350	else if (iflag_pbl.ge.8) then
	351	c====================================================================
	352	c Yamada 2.5 a la Didi
	353	c====================================================================
	354
	355
	356	c Calcul de l, km, au pas precedent
	357	do k=2,klev
	358	do ig=1,ngrid
	359	c print*,'SMML=',sm(ig,k),l(ig,k)
	360	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	361	if (delta(ig,k).lt.1.e-20) then
	362	c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
	363	delta(ig,k)=1.e-20
	364	endif
	365	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	366	aa0=
	367	s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
	368	aa1=
	369	s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
	370	c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
	371	aa(ig,k)=aa1dt/(delta(ig,k)l(ig,k))
	372	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	373	qpre=sqrt(q2(ig,k))
	374	if (iflag_pbl.eq.8 ) then
	375	if (aa(ig,k).gt.0.) then
	376	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
	377	else
	378	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	379	endif
	380	else ! iflag_pbl=9
	381	if (aa(ig,k)*qpre.gt.0.9) then
	382	q2(ig,k)=(qpre10.)*2
	383	else
	384	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	385	endif
	386	endif
	387	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
	388	c print,'Q2 L=',k,q2(ig,k),qpreqpre
	389	enddo
	390	enddo
	391
	392	endif ! Fin du cas 8
	393
	394	c print*,'OK8'
	395
	396	c====================================================================
	397	c Calcul des coefficients de mélange
	398	c====================================================================
	399	do k=2,klev
	400	c print*,'k=',k
	401	do ig=1,ngrid
	402	cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
	403	zq=sqrt(q2(ig,k))
	404	km(ig,k)=l(ig,k)zqsm(ig,k)
	405	kn(ig,k)=km(ig,k)*alpha(ig,k)
	406	kq(ig,k)=l(ig,k)zq0.2
	407	c print*,'KML=',km(ig,k),kn(ig,k)
	408	enddo
	409	enddo
	410
	411	c if (iflag_pbl.ge.7..and.0.eq.1) then
	412	c q2(:,1)=q2(:,2)
	413	c call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
	414	c endif
	415
	416	c Traitement des cas noctrunes avec l'introduction d'une longueur
	417	c minilale.
	418
	419	c====================================================================
	420	c Traitement particulier pour les cas tres stables.
	421	c D'apres Holtslag Boville.
	422
	423	print*,'YAMADA4 0'
	424
	425	do ig=1,ngrid
	426	coriol(ig)=1.e-4
	427	pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
	428	enddo
	429
	430	print*,'pblhmin ',pblhmin
	431	CTest a remettre 21 11 02
	432	c test abd 13 05 02 if(0.eq.1) then
	433	if(1.eq.1) then
	434	do k=2,klev
	435	do ig=1,klon
	436	if (teta(ig,2).gt.teta(ig,1)) then
	437	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	438	kmin=kapzlev(ig,k)qmin
	439	else
	440	kmin=-1. ! kmin n'est utilise que pour les SL stables.
	441	endif
	442	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
	443	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
	444	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
	445	kn(ig,k)=kmin
	446	km(ig,k)=kmin
	447	kq(ig,k)=kmin
	448	c la longueur de melange est suposee etre l= kap z
	449	c K=l q Sm d'ou q2=(K/l Sm)**2
	450	q2(ig,k)=(qmin/sm(ig,k))**2
	451	endif
	452	enddo
	453	enddo
	454	endif
	455
	456	print*,'YAMADA4 1'
	457	c Diagnostique pour stokage
	458
	459	rino=rif
	460	smyam(:,1:klev)=sm(:,1:klev)
	461	styam=sm(:,1:klev)*alpha(:,1:klev)
	462	lyam(1:klon,1:klev)=l(:,1:klev)
	463	knyam(1:klon,1:klev)=kn(:,1:klev)
	464
	465	c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
	466
	467	if(1.eq.0)then
	468	w2yam=q2(:,1:klev)*0.24
	469	s +lyam(:,1:klev)5.17kn(:,1:klev)*n2(:,1:klev)
	470	s /sqrt(q2(:,1:klev))
	471
	472	t2yam=9.1kn(:,1:klev)dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
	473	s *lyam(:,1:klev)
	474	endif
	475
	476	c print*,'OKFIN'
	477	first=.false.
	478	return
	479	end

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

Download in other formats:

Original Format