Context Navigation

yamada4.F @ 922

Last change on this file since 922 was 878, checked in by Laurent Fairhead, 17 years ago
Bascule de la physique du LMD vers la physique avec thermiques LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.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	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	if(1.eq.0)then
461	rino=rif
462	smyam(1:ngrid,1)=0.
463	styam(1:ngrid,1)=0.
464	lyam(1:ngrid,1)=0.
465	knyam(1:ngrid,1)=0.
466	w2yam(1:ngrid,1)=0.
467	t2yam(1:ngrid,1)=0.
468
469	smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
470	styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
471	lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
472	knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
473
474	c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
475
476	w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24
477	s +lyam(1:ngrid,2:klev)5.17kn(1:ngrid,2:klev)
478	s *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
479
480	t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev)
481	s dtetadz(1:ngrid,2:klev)*2
482	s /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
483	endif
484
485	c print*,'OKFIN'
486	first=.false.
487	return
488	end

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

Download in other formats:

Original Format