Context Navigation

yamada4.F @ 1461

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

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

Download in other formats:

Original Format