Context Navigation

yamada4.F @ 4104

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

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

Download in other formats:

Original Format