Context Navigation

yamada_c.F90 @ 5440

Last change on this file since 5440 was 1761, checked in by Laurent Fairhead, 12 years ago

New version of Mellor et Yamada pronostic TKE

... based on energy transfer from the mean state.

The new version is yamada_c.
It must be called after vertical diffusion rather than just before
since the source terms u_z w'u' + v_z w'theta' and g/theta w'theta'
are diagnosed from the vertical diffusion (as energy loss from the mean
state) rather than computed as K (u_z^2+v_z2) or g/\theta K theta_z
(where _z means vertical derivative).
The call to this version is controled by iflag_pbl.

iflag_pbl = 20 : with TKE computed at interfaces between layers
iflag_pbl = 25 : with TKE computed within the layer
In both cases, the dissipation of turbulence is translated into heat, and
passed to the physics as dtdiss (temperature tendency due to dissipation

of TKE).

The diffusion coefficient being computed after dissipation, it must be
kept for diffusion at the next time step, and thus be stored in the
restartphy file.
This coefficient must be computed and stored for each sub-surface.

A new way of managing subsurface variables is introduced.
For arrays of the form X(:,nbsrf) are extented to X(:,nbsrf+1), where
is_ave=nbsrf+1, is an additional sub-surface which contains the averaged values.

coef_diff_turb_mod.F90 : change of flags.
ener_conserv.F90 : energy conservation must not be applied in those

cases

indicesol.h : definition of is_ave
pbl_surface_mod.F90 : call to yamada_c and changes in the management of

coefh/coefm

physiq.F : Change in the initialisation of pmflxr/s and

modified calls to pbl_surface_mod (introduction
of dtdiss, initialisation
of pbl_tke and coefh in 1D).

phys_local_var_mod.F90 : declaration of d_t_diss
phys_output_mod.F90 : new outputs (bils_tke and bils_diss) and

coefh->coefh(:,:,is_ave)

phys_state_var_mod.F90 : modified declaration for coefh and coefm

(nbsrf -> nbsrf+1)

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