Context Navigation

yamada4.F @ 4237

Last change on this file since 4237 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

Line
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	use dimphy
8	IMPLICIT NONE
9	c.......................................................................
10	cym#include "dimensions.h"
11	cym#include "dimphy.h"
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
67	cym PARAMETER (nlay=klev)
68	cym PARAMETER (nlev=klev+1)
69
70	logical first
71	integer ipas
72	save first,ipas
73	data first,ipas/.true.,0/
74	c$OMP THREADPRIVATE( first,ipas)
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
85	real l(klon,klev+1)
86	real,allocatable,save :: l0(:)
87	c$OMP THREADPRIVATE(l0)
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	c$OMP THREADPRIVATE(ric,rifc,b1,kap)
95	real frif,falpha,fsm
96	real fl,zzz,zl0,zq2,zn2
97
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)
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
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
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====================================================================
398	c Calcul des coefficients de m�ange
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,ngrid
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