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

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	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