Context Navigation

yamada_c.F90 @ 5434

Last change on this file since 5434 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 22.9 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE yamada_c(ngrid,timestep,plev,play &
5	& ,pu,pv,pt,d_u,d_v,d_t,cd,q2,km,kn,kq,d_t_diss,ustar &
6	& ,iflag_pbl,okiophys)
7	use dimphy
8	IMPLICIT NONE
9	#include "iniprint.h"
10	#include "YOMCST.h"
11	!.......................................................................
12	!ym#include "dimensions.h"
13	!ym#include "dimphy.h"
14	!.......................................................................
15	!
16	! timestep : pas de temps
17	! g : g
18	! zlev : altitude a chaque niveau (interface inferieure de la couche
19	! de meme indice)
20	! zlay : altitude au centre de chaque couche
21	! u,v : vitesse au centre de chaque couche
22	! (en entree : la valeur au debut du pas de temps)
23	! teta : temperature potentielle au centre de chaque couche
24	! (en entree : la valeur au debut du pas de temps)
25	! cd : cdrag
26	! (en entree : la valeur au debut du pas de temps)
27	! q2 : $q^2$ au bas de chaque couche
28	! (en entree : la valeur au debut du pas de temps)
29	! (en sortie : la valeur a la fin du pas de temps)
30	! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
31	! couche)
32	! (en sortie : la valeur a la fin du pas de temps)
33	! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
34	! (en sortie : la valeur a la fin du pas de temps)
35	!
36	! iflag_pbl doit valoir entre 6 et 9
37	! l=6, on prend systematiquement une longueur d'equilibre
38	! iflag_pbl=6 : MY 2.0
39	! iflag_pbl=7 : MY 2.0.Fournier
40	! iflag_pbl=8/9 : MY 2.5
41	! iflag_pbl=8 with special obsolete treatments for convergence
42	! with Cmpi5 NPv3.1 simulations
43	! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact
44	! iflag_pbl=12 = 11 with vertical diffusion off q2
45	!
46	! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
47	! Correction for very stable PBLs (iflag_pbl=10 and 11)
48	! iflag_pbl=8 converges numerically with NPv3.1
49	! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
50	! -> the model can run with longer time-steps.
51	!.......................................................................
52
53	REAL, DIMENSION(klon,klev) :: d_u,d_v,d_t
54	REAL, DIMENSION(klon,klev) :: pu,pv,pt
55	REAL, DIMENSION(klon,klev) :: d_t_diss
56	INTEGER okiophys
57
58	REAL timestep
59	real plev(klon,klev+1)
60	real play(klon,klev)
61	real ustar(klon)
62	real kmin,qmin,pblhmin(klon),coriol(klon)
63	REAL zlev(klon,klev+1)
64	REAL zlay(klon,klev)
65	REAL zu(klon,klev)
66	REAL zv(klon,klev)
67	REAL zt(klon,klev)
68	REAL teta(klon,klev)
69	REAL cd(klon)
70	REAL q2(klon,klev+1),qpre
71	REAL unsdz(klon,klev)
72	REAL unsdzdec(klon,klev+1)
73
74	REAL km(klon,klev+1)
75	REAL kmpre(klon,klev+1),tmp2
76	REAL mpre(klon,klev+1)
77	REAL kn(klon,klev+1)
78	REAL kq(klon,klev+1)
79	real ff(klon,klev+1),delta(klon,klev+1)
80	real aa(klon,klev+1),aa0,aa1
81	integer iflag_pbl,ngrid
82	integer nlay,nlev
83
84	logical first
85	integer ipas
86	save first,ipas
87	!FH/IM data first,ipas/.true.,0/
88	data first,ipas/.false.,0/
89	!$OMP THREADPRIVATE( first,ipas)
90
91	integer ig,k
92
93
94	real ri,zrif,zalpha,zsm,zsn
95	real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
96
97	real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
98	REAL, DIMENSION(klon,klev+1) :: km2,kn2,sqrtq
99	real dtetadz(klon,klev+1)
100	real m2cstat,mcstat,kmcstat
101	real l(klon,klev+1)
102	real leff(klon,klev+1)
103	real,allocatable,save :: l0(:)
104	!$OMP THREADPRIVATE(l0)
105	real sq(klon),sqz(klon),zz(klon,klev+1)
106	integer iter
107
108	real ric,rifc,b1,kap
109	save ric,rifc,b1,kap
110	data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
111	!$OMP THREADPRIVATE(ric,rifc,b1,kap)
112	real frif,falpha,fsm
113	real fl,zzz,zl0,zq2,zn2
114
115	real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
116	real lyam(klon,klev),knyam(klon,klev)
117	real w2yam(klon,klev),t2yam(klon,klev)
118	logical,save :: firstcall=.true.
119
120	REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt
121	REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt
122	REAL, DIMENSION(klon,klev) :: exner,masse
123	REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg
124
125	!$OMP THREADPRIVATE(firstcall)
126	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
127	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
128	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
129	fl(zzz,zl0,zq2,zn2)= &
130	& max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig)) &
131	& ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
132
133
134	if (firstcall) then
135	allocate(l0(klon))
136	#ifdef IOPHYS
137	call iophys_ini
138	#endif
139	firstcall=.false.
140	endif
141
142
143	#ifdef IOPHYS
144	if (okiophys==1) then
145	call iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev))
146	call iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev))
147	endif
148	#endif
149
150	nlay=klev
151	nlev=klev+1
152
153	!-------------------------------------------------------------------------
154	! Computation of conservative source terms from the turbulent tendencies
155	!-------------------------------------------------------------------------
156
157
158	zu(:,:)=pu(:,:)+0.5*d_u(:,:)
159	zv(:,:)=pv(:,:)+0.5*d_v(:,:)
160	zt(:,:)=pt(:,:)+0.5*d_t(:,:)
161	do k=1,klev
162	exner(:,k)=(play(:,k)/plev(:,1))**RKAPPA
163	masse(:,k)=(plev(:,k)-plev(:,k+1))/RG
164	enddo
165	teta(:,:)=zt(:,:)/exner(:,:)
166
167	! Atmospheric mass at layer interfaces, where the TKE is computed
168	masseb(:,:)=0.
169	do k=1,klev
170	masseb(:,k)=masseb(:,k)+masse(:,k)
171	masseb(:,k+1)=masseb(:,k+1)+masse(:,k)
172	enddo
173	masseb(:,:)=0.5*masseb(:,:)
174
175
176
177	zlev(:,1)=0.
178	zlay(:,1)=RCPDteta(:,1)(1.-exner(:,1))
179	do k=1,klev-1
180	zlay(:,k+1)=zlay(:,k)+0.5RCPD(teta(:,k)+teta(:,k+1))*(exner(:,k)-exner(:,k+1))/RG
181	zlev(:,k)=0.5*(zlay(:,k)+zlay(:,k+1)) ! PASBO
182	enddo
183
184	fluxu(:,klev+1)=0.
185	fluxv(:,klev+1)=0.
186	fluxt(:,klev+1)=0.
187
188	do k=klev,1,-1
189	fluxu(:,k)=fluxu(:,k+1)+masse(:,k)*d_u(:,k)
190	fluxv(:,k)=fluxv(:,k+1)+masse(:,k)*d_v(:,k)
191	fluxt(:,k)=fluxt(:,k+1)+masse(:,k)*d_t(:,k)/exner(:,k) ! Flux de theta
192	enddo
193
194	dddu(:,1)=2zu(:,1)fluxu(:,1)
195	dddv(:,1)=2zv(:,1)fluxv(:,1)
196	dddt(:,1)=(exner(:,1)-1.)*fluxt(:,1)
197
198	do k=2,klev
199	dddu(:,k)=(zu(:,k)-zu(:,k-1))*fluxu(:,k)
200	dddv(:,k)=(zv(:,k)-zv(:,k-1))*fluxv(:,k)
201	dddt(:,k)=(exner(:,k)-exner(:,k-1))*fluxt(:,k)
202	enddo
203	dddu(:,klev+1)=0.
204	dddv(:,klev+1)=0.
205	dddt(:,klev+1)=0.
206
207	#ifdef IOPHYS
208	if (okiophys==1) then
209	call iophys_ecrit('zlay',klev,'Geop','m',zlay)
210	call iophys_ecrit('teta',klev,'teta','K',teta)
211	call iophys_ecrit('temp',klev,'temp','K',zt)
212	call iophys_ecrit('pt',klev,'temp','K',pt)
213	call iophys_ecrit('d_u',klev,'d_u','m/s2',d_u)
214	call iophys_ecrit('d_v',klev,'d_v','m/s2',d_v)
215	call iophys_ecrit('d_t',klev,'d_t','K/s',d_t)
216	call iophys_ecrit('exner',klev,'exner','',exner)
217	call iophys_ecrit('masse',klev,'masse','',masse)
218	call iophys_ecrit('masseb',klev,'masseb','',masseb)
219	call iophys_ecrit('Cm2',klev,'m2 conserv','m/s',(dddu(:,1:klev)+dddv(:,1:klev))/(masseb(:,1:klev)*timestep))
220	call iophys_ecrit('Cn2',klev,'m2 conserv','m/s',(rcpd*dddt(:,1:klev)/masseb(:,1:klev))/timestep)
221	call iophys_ecrit('rifc',klev,'rif conservative','',rcpd*dddt(:,1:klev)/min(dddu(:,1:klev)+dddv(:,1:klev),-1.e-20))
222	endif
223	#endif
224
225
226
227	ipas=ipas+1
228
229
230	!.......................................................................
231	! les increments verticaux
232	!.......................................................................
233	!
234	!!!!!! allerte !!!!!c
235	!!!!!! zlev n'est pas declare a nlev !!!!!c
236	!!!!!! ---->
237	DO ig=1,ngrid
238	zlev(ig,nlev)=zlay(ig,nlay) &
239	& +( zlay(ig,nlay) - zlev(ig,nlev-1) )
240	ENDDO
241	!!!!!! <----
242	!!!!!! allerte !!!!!c
243	!
244	DO k=1,nlay
245	DO ig=1,ngrid
246	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
247	ENDDO
248	ENDDO
249	DO ig=1,ngrid
250	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
251	ENDDO
252	DO k=2,nlay
253	DO ig=1,ngrid
254	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
255	ENDDO
256	ENDDO
257	DO ig=1,ngrid
258	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
259	ENDDO
260	!
261	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
262	! Computing M^2, N^2, Richardson numbers, stability functions
263	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
264
265
266	do k=2,klev
267	do ig=1,ngrid
268	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
269	m2(ig,k)=((zu(ig,k)-zu(ig,k-1))2+(zv(ig,k)-zv(ig,k-1))2)/(dz(ig,k)*dz(ig,k))
270	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
271	n2(ig,k)=RG2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
272	! n2(ig,k)=0.
273	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
274	if (ri.lt.ric) then
275	rif(ig,k)=frif(ri)
276	else
277	rif(ig,k)=rifc
278	endif
279	if(rif(ig,k).lt.0.16) then
280	alpha(ig,k)=falpha(rif(ig,k))
281	sm(ig,k)=fsm(rif(ig,k))
282	else
283	alpha(ig,k)=1.12
284	sm(ig,k)=0.085
285	endif
286	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
287	enddo
288	enddo
289
290
291
292	!====================================================================
293	! Computing the mixing length
294	!====================================================================
295
296	! Mise a jour de l0
297	if (iflag_pbl==8.or.iflag_pbl==10) then
298
299	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
300	! Iterative computation of l0
301	! This version is kept for iflag_pbl only for convergence
302	! with NPv3.1 Cmip5 simulations
303	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304
305	do ig=1,ngrid
306	sq(ig)=1.e-10
307	sqz(ig)=1.e-10
308	enddo
309	do k=2,klev-1
310	do ig=1,ngrid
311	zq=sqrt(q2(ig,k))
312	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
313	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
314	enddo
315	enddo
316	do ig=1,ngrid
317	l0(ig)=0.2*sqz(ig)/sq(ig)
318	enddo
319	do k=2,klev
320	do ig=1,ngrid
321	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
322	enddo
323	enddo
324	! print*,'L0 cas 8 ou 10 ',l0
325
326	else
327
328	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329	! In all other case, the assymptotic mixing length l0 is imposed (100m)
330	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
331
332	l0(:)=150.
333	do k=2,klev
334	do ig=1,ngrid
335	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
336	enddo
337	enddo
338	! print*,'L0 cas autres ',l0
339
340	endif
341
342
343	#ifdef IOPHYS
344	if (okiophys==1) then
345	call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev))
346	call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev))
347	call iophys_ecrit('Km2',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
348	call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev))
349	endif
350	#endif
351
352
353	IF (iflag_pbl<20) then
354	! For diagnostics only
355	RETURN
356
357	ELSE
358
359	! print*,'OK1'
360
361	! Evolution of TKE under source terms K M2 and K N2
362	leff(:,:)=max(l(:,:),1.)
363	IF (iflag_pbl==29) THEN
364	km2(:,:)=km(:,:)*m2(:,:)
365	kn2(:,:)=kn2(:,:)*rif(:,:)
366	ELSEIF (iflag_pbl==25) THEN
367	DO k=1,klev
368	km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) &
369	& /(masse(:,k)*timestep)
370	kn2(:,k)=rcpd0.5(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep)
371	leff(:,k)=0.5*(leff(:,k)+leff(:,k+1))
372	ENDDO
373	km2(:,klev+1)=0. ; kn2(:,klev+1)=0.
374	ELSE
375	km2(:,:)=-(dddu(:,:)+dddv(:,:))/(masseb(:,:)*timestep)
376	kn2(:,:)=rcpddddt(:,:)/(masseb(:,:)timestep)
377	ENDIF
378	q2neg(:,:)=q2(:,:)+timestep*(km2(:,:)-kn2(:,:))
379	q2(:,:)=min(max(q2neg(:,:),1.e-10),1.e4)
380
381	! Dissipation of TKE
382	q2old(:,:)=q2(:,:)
383	q2(:,:)=1./(1./sqrt(q2(:,:))+timestep/(2leff(:,:)b1))
384	q2(:,:)=q2(:,:)*q2(:,:)
385	IF (iflag_pbl<=24) THEN
386	DO k=1,klev
387	d_t_diss(:,k)=(masseb(:,k)(q2neg(:,k)-q2(:,k))+masseb(:,k+1)(q2neg(:,k+1)-q2(:,k+1)))/(2.rcpdmasse(:,k))
388	ENDDO
389	ELSE IF (iflag_pbl<=27) THEN
390	DO k=1,klev
391	d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd
392	ENDDO
393	ENDIF
394	! print*,'iflag_pbl ',d_t_diss
395
396
397	! Compuation of stability functions
398	IF (iflag_pbl/=29) THEN
399	DO k=1,klev
400	DO ig=1,ngrid
401	IF (ABS(km2(ig,k))<=1.e-20) THEN
402	rif(ig,k)=0.
403	ELSE
404	rif(ig,k)=min(kn2(ig,k)/km2(ig,k),rifc)
405	ENDIF
406	IF (rif(ig,k).lt.0.16) THEN
407	alpha(ig,k)=falpha(rif(ig,k))
408	sm(ig,k)=fsm(rif(ig,k))
409	else
410	alpha(ig,k)=1.12
411	sm(ig,k)=0.085
412	endif
413	ENDDO
414	ENDDO
415	ENDIF
416
417	! Computation of turbulent diffusivities
418	IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN
419	DO k=2,klev
420	sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1)))
421	ENDDO
422	ELSE
423	DO k=2,klev
424	sqrtq(:,k)=sqrt(q2(:,k))
425	ENDDO
426	ENDIF
427	DO k=2,klev
428	DO ig=1,ngrid
429	km(ig,k)=leff(ig,k)sqrtq(ig,k)sm(ig,k)
430	kn(ig,k)=km(ig,k)*alpha(ig,k)
431	kq(ig,k)=leff(ig,k)zq0.2
432	! print*,q2(ig,k),zq,km(ig,k)
433	ENDDO
434	ENDDO
435
436
437
438	#ifdef IOPHYS
439	if (okiophys==1) then
440	call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev))
441	call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev))
442	call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
443	call iophys_ecrit('q2neg',klev,'KTE non bornee','m2/s',q2neg(:,1:klev))
444	call iophys_ecrit('alpha',klev,'alpha','',alpha(:,1:klev))
445	call iophys_ecrit('sm',klev,'sm','',sm(:,1:klev))
446	call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
447	call iophys_ecrit('kmf',klev,'Kz final','m2/s',km(:,1:klev))
448	call iophys_ecrit('knf',klev,'Kz final','m2/s',kn(:,1:klev))
449	call iophys_ecrit('kqf',klev,'Kz final','m2/s',kq(:,1:klev))
450	endif
451	#endif
452
453	ENDIF
454
455
456	! print*,'OK2'
457	RETURN
458	!====================================================================
459	! Yamada 2.0
460	!====================================================================
461	if (iflag_pbl.eq.6) then
462
463	do k=2,klev
464	q2(:,k)=l(:,k)*2zz(:,k)
465	enddo
466
467
468	else if (iflag_pbl.eq.7) then
469	!====================================================================
470	! Yamada 2.Fournier
471	!====================================================================
472
473	! Calcul de l, km, au pas precedent
474	do k=2,klev
475	do ig=1,ngrid
476	! print*,'SMML=',sm(ig,k),l(ig,k)
477	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
478	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
479	mpre(ig,k)=sqrt(m2(ig,k))
480	! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
481	enddo
482	enddo
483
484	do k=2,klev-1
485	do ig=1,ngrid
486	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
487	mcstat=sqrt(m2cstat)
488
489	! print*,'M2 L=',k,mpre(ig,k),mcstat
490	!
491	! -----{puis on ecrit la valeur de q qui annule l'equation de m
492	! supposee en q3}
493	!
494	IF (k.eq.2) THEN
495	kmcstat=1.E+0 / mcstat &
496	& ( unsdz(ig,k)kmpre(ig,k+1) &
497	& *mpre(ig,k+1) &
498	& +unsdz(ig,k-1) &
499	& *cd(ig) &
500	& ( sqrt(zu(ig,3)2+zv(ig,3)*2) &
501	& -mcstat/unsdzdec(ig,k) &
502	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2) &
503	& /( unsdz(ig,k)+unsdz(ig,k-1) )
504	ELSE
505	kmcstat=1.E+0 / mcstat &
506	& ( unsdz(ig,k)kmpre(ig,k+1) &
507	& *mpre(ig,k+1) &
508	& +unsdz(ig,k-1)*kmpre(ig,k-1) &
509	& *mpre(ig,k-1) ) &
510	& /( unsdz(ig,k)+unsdz(ig,k-1) )
511	ENDIF
512	! print*,'T2 L=',k,tmp2
513	tmp2=kmcstat &
514	& /( sm(ig,k)/q2(ig,k) ) &
515	& /l(ig,k)
516	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
517	! print*,'Q2 L=',k,q2(ig,k)
518	!
519	enddo
520	enddo
521
522	else if (iflag_pbl==8.or.iflag_pbl==9) then
523	!====================================================================
524	! Yamada 2.5 a la Didi
525	!====================================================================
526
527
528	! Calcul de l, km, au pas precedent
529	do k=2,klev
530	do ig=1,ngrid
531	! print*,'SMML=',sm(ig,k),l(ig,k)
532	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
533	if (delta(ig,k).lt.1.e-20) then
534	! print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
535	delta(ig,k)=1.e-20
536	endif
537	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
538	aa0=(m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
539	aa1=(m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
540	! abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
541	aa(ig,k)=aa1timestep/(delta(ig,k)l(ig,k))
542	! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
543	qpre=sqrt(q2(ig,k))
544	! if (iflag_pbl.eq.8 ) then
545	if (aa(ig,k).gt.0.) then
546	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
547	else
548	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
549	endif
550	! else ! iflag_pbl=9
551	! if (aa(ig,k)*qpre.gt.0.9) then
552	! q2(ig,k)=(qpre10.)*2
553	! else
554	! q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
555	! endif
556	! endif
557	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
558	! print,'Q2 L=',k,q2(ig,k),qpreqpre
559	enddo
560	enddo
561
562	else if (iflag_pbl>=10) then
563
564	! print*,'Schema mixte D'
565	! print*,'Longueur ',l(:,:)
566	do k=2,klev-1
567	l(:,k)=max(l(:,k),1.)
568	km(:,k)=l(:,k)sqrt(q2(:,k))sm(:,k)
569	q2(:,k)=q2(:,k)+timestepkm(:,k)m2(:,k)*(1.-rif(:,k))
570	q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
571	q2(:,k)=1./(1./sqrt(q2(:,k))+timestep/(2l(:,k)b1))
572	q2(:,k)=q2(:,k)*q2(:,k)
573	enddo
574
575
576	else
577	stop'Cas nom prevu dans yamada4'
578
579	endif ! Fin du cas 8
580
581	! print*,'OK8'
582
583	!====================================================================
584	! Calcul des coefficients de m�ange
585	!====================================================================
586	do k=2,klev
587	! print*,'k=',k
588	do ig=1,ngrid
589	!abde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
590	zq=sqrt(q2(ig,k))
591	km(ig,k)=l(ig,k)zqsm(ig,k)
592	kn(ig,k)=km(ig,k)*alpha(ig,k)
593	kq(ig,k)=l(ig,k)zq0.2
594	! print*,'KML=',km(ig,k),kn(ig,k)
595	enddo
596	enddo
597
598	! Transport diffusif vertical de la TKE.
599	if (iflag_pbl.ge.12) then
600	! print*,'YAMADA VDIF'
601	q2(:,1)=q2(:,2)
602	call vdif_q2(timestep,RG,RD,ngrid,plev,zt,kq,q2)
603	endif
604
605	! Traitement des cas noctrunes avec l'introduction d'une longueur
606	! minilale.
607
608	!====================================================================
609	! Traitement particulier pour les cas tres stables.
610	! D'apres Holtslag Boville.
611
612	if (prt_level>1) THEN
613	print*,'YAMADA4 0'
614	endif !(prt_level>1) THEN
615	do ig=1,ngrid
616	coriol(ig)=1.e-4
617	pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
618	enddo
619
620	! print*,'pblhmin ',pblhmin
621	!Test a remettre 21 11 02
622	! test abd 13 05 02 if(0.eq.1) then
623	if(1==1) then
624	if(iflag_pbl==8.or.iflag_pbl==10) then
625
626	do k=2,klev
627	do ig=1,ngrid
628	if (teta(ig,2).gt.teta(ig,1)) then
629	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
630	kmin=kapzlev(ig,k)qmin
631	else
632	kmin=-1. ! kmin n'est utilise que pour les SL stables.
633	endif
634	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
635	! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
636	! s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
637	kn(ig,k)=kmin
638	km(ig,k)=kmin
639	kq(ig,k)=kmin
640	! la longueur de melange est suposee etre l= kap z
641	! K=l q Sm d'ou q2=(K/l Sm)**2
642	q2(ig,k)=(qmin/sm(ig,k))**2
643	endif
644	enddo
645	enddo
646
647	else
648
649	do k=2,klev
650	do ig=1,ngrid
651	if (teta(ig,2).gt.teta(ig,1)) then
652	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
653	kmin=kapzlev(ig,k)qmin
654	else
655	kmin=-1. ! kmin n'est utilise que pour les SL stables.
656	endif
657	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
658	! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
659	! s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
660	kn(ig,k)=kmin
661	km(ig,k)=kmin
662	kq(ig,k)=kmin
663	! la longueur de melange est suposee etre l= kap z
664	! K=l q Sm d'ou q2=(K/l Sm)**2
665	sm(ig,k)=1.
666	alpha(ig,k)=1.
667	q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
668	zq=sqrt(q2(ig,k))
669	km(ig,k)=l(ig,k)zqsm(ig,k)
670	kn(ig,k)=km(ig,k)*alpha(ig,k)
671	kq(ig,k)=l(ig,k)zq0.2
672	endif
673	enddo
674	enddo
675	endif
676
677	endif
678
679	if (prt_level>1) THEN
680	print*,'YAMADA4 1'
681	endif !(prt_level>1) THEN
682	! Diagnostique pour stokage
683
684	if(1.eq.0)then
685	rino=rif
686	smyam(1:ngrid,1)=0.
687	styam(1:ngrid,1)=0.
688	lyam(1:ngrid,1)=0.
689	knyam(1:ngrid,1)=0.
690	w2yam(1:ngrid,1)=0.
691	t2yam(1:ngrid,1)=0.
692
693	smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
694	styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
695	lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
696	knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
697
698	! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
699
700	w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24 &
701	& +lyam(1:ngrid,2:klev)5.17kn(1:ngrid,2:klev) &
702	& *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
703
704	t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev) &
705	& dtetadz(1:ngrid,2:klev)*2 &
706	& /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
707	endif
708
709	! print*,'OKFIN'
710	first=.false.
711	return
712	end

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

Download in other formats:

Original Format