Context Navigation

yamada4.F @ 537

Last change on this file since 537 was 102, checked in by slebonnois, 14 years ago
SL : corrections et modifications dans phytitan correspondant a celles faites apres compilation Venus. Titan pas encore compile.
File size: 15.7 KB

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

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

Download in other formats:

Original Format