Context Navigation

yamada4.F90 @ 3513

Last change on this file since 3513 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

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
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 36.0 KB

Line
1	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
3	SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
4	cd, tke, km, kn, kq, ustar, iflag_pbl, drgpro)
5
6	USE dimphy
7	USE ioipsl_getin_p_mod, ONLY : getin_p
8
9	IMPLICIT NONE
10	include "iniprint.h"
11	! .......................................................................
12	! ym#include "dimensions.h"
13	! ym#include "dimphy.h"
14	! ************************************************************************************************
15	!
16	! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5
17	!
18	! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model
19	! Yamada T.
20	! J. Atmos. Sci, 40, 91-106, 1983
21	!
22	!************************************************************************************************
23	! Input :
24	!'======
25	! ni: indice horizontal sur la grille de base, non restreinte
26	! nsrf: type de surface
27	! ngrid: nombre de mailles concern??es sur l'horizontal
28	! dt : pas de temps
29	! g : g
30	! rconst: constante de l'air sec
31	! zlev : altitude a chaque niveau (interface inferieure de la couche
32	! de meme indice)
33	! zlay : altitude au centre de chaque couche
34	! u,v : vitesse au centre de chaque couche
35	! (en entree : la valeur au debut du pas de temps)
36	! teta : temperature potentielle virtuelle au centre de chaque couche
37	! (en entree : la valeur au debut du pas de temps)
38	! cd : cdrag pour la quantit?? de mouvement
39	! (en entree : la valeur au debut du pas de temps)
40	! ustar: vitesse de friction calcul??e par une formule diagnostique
41	! iflag_pbl: flag pour choisir des options du sch??ma de turbulence
42	!
43	! iflag_pbl doit valoir entre 6 et 9
44	! l=6, on prend systematiquement une longueur d'equilibre
45	! iflag_pbl=6 : MY 2.0
46	! iflag_pbl=7 : MY 2.0.Fournier
47	! iflag_pbl=8/9 : MY 2.5
48	! iflag_pbl=8 with special obsolete treatments for convergence
49	! with Cmpi5 NPv3.1 simulations
50	! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact
51	! iflag_pbl=12 = 11 with vertical diffusion off q2
52	!
53	! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
54	! Correction for very stable PBLs (iflag_pbl=10 and 11)
55	! iflag_pbl=8 converges numerically with NPv3.1
56	! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
57	! -> the model can run with longer time-steps.
58	! 2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
59	! On met tke (=q2/2) en entr??e plut??t que q2
60	! On corrige l'update de la tke
61	!
62	! Inpout/Output :
63	!==============
64	! tke : tke au bas de chaque couche
65	! (en entree : la valeur au debut du pas de temps)
66	! (en sortie : la valeur a la fin du pas de temps)
67
68	! Outputs:
69	!==========
70	! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
71	! couche)
72	! (en sortie : la valeur a la fin du pas de temps)
73	! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
74	! (en sortie : la valeur a la fin du pas de temps)
75	!
76	!.......................................................................
77
78	!=======================================================================
79	! Declarations:
80	!=======================================================================
81
82
83	! Inputs/Outputs
84	!----------------
85	REAL dt, g, rconst
86	REAL plev(klon, klev+1), temp(klon, klev)
87	REAL ustar(klon)
88	REAL kmin, qmin, pblhmin(klon), coriol(klon)
89	REAL zlev(klon, klev+1)
90	REAL zlay(klon, klev)
91	REAL u(klon, klev)
92	REAL v(klon, klev)
93	REAL teta(klon, klev)
94	REAL cd(klon)
95	REAL tke(klon, klev+1)
96	REAL unsdz(klon, klev)
97	REAL unsdzdec(klon, klev+1)
98	REAL kn(klon, klev+1)
99	REAL km(klon, klev+1)
100	INTEGER iflag_pbl, ngrid, nsrf
101	INTEGER ni(klon)
102
103	!FC
104	REAL drgpro(klon,klev)
105	REAL winds(klon,klev)
106
107	! Local
108	!-------
109
110	INCLUDE "clesphys.h"
111
112	REAL q2(klon, klev+1)
113	REAL kmpre(klon, klev+1), tmp2, qpre
114	REAL mpre(klon, klev+1)
115	REAL kq(klon, klev+1)
116	REAL ff(klon, klev+1), delta(klon, klev+1)
117	REAL aa(klon, klev+1), aa0, aa1
118	INTEGER nlay, nlev
119
120	LOGICAL,SAVE :: hboville=.TRUE.
121	REAL,SAVE :: viscom,viscoh
122	!$OMP THREADPRIVATE( hboville,viscom,viscoh)
123	INTEGER ig, k
124	REAL ri, zrif, zalpha, zsm, zsn
125	REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
126	REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
127	REAL dtetadz(klon, klev+1)
128	REAL m2cstat, mcstat, kmcstat
129	REAL l(klon, klev+1)
130	REAL zz(klon, klev+1)
131	INTEGER iter
132	REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev)
133	REAL :: disseff
134
135	REAL,SAVE :: ric0,ric,rifc, b1, kap
136	!$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
137	DATA b1, kap/16.6, 0.4/
138	REAL,SAVE :: seuilsm, seuilalpha
139	!$OMP THREADPRIVATE(seuilsm, seuilalpha)
140	REAL,SAVE :: lmixmin
141	!$OMP THREADPRIVATE(lmixmin)
142	LOGICAL, SAVE :: new_yamada4
143	INTEGER, SAVE :: yamada4_num
144	!$OMP THREADPRIVATE(new_yamada4,yamada4_num)
145	REAL, SAVE :: yun,ydeux
146	!$OMP THREADPRIVATE(yun,ydeux)
147
148	REAL frif, falpha, fsm
149	REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
150	lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
151	LOGICAL, SAVE :: firstcall = .TRUE.
152	!$OMP THREADPRIVATE(firstcall)
153
154
155
156	! Fonctions utilis??es
157	!--------------------
158
159	frif(ri) = 0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
160	falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
161	fsm(ri) = 1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
162
163
164	IF (firstcall) THEN
165	! Seuil dans le code de turbulence
166	new_yamada4=.false.
167	CALL getin_p('new_yamada4',new_yamada4)
168
169	IF (new_yamada4) THEN
170	! Corrections et reglages issus du travail de these d'Etienne Vignon.
171	ric=0.143 ! qui donne des valeurs proches des seuils proposes
172	! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
173	! sm=1.1213 (1.12 dans Y83)
174	CALL getin_p('yamada4_ric',ric)
175	ric0=0.19489 ! ric=0.195 originalement, mais produisait sm<0
176	ric=min(ric,ric0) ! Au dela de ric0, sm devient n??gatif
177	rifc=frif(ric)
178	seuilsm=fsm(frif(ric))
179	seuilalpha=falpha(frif(ric))
180	yun=1.
181	ydeux=2.
182	hboville=.FALSE.
183	viscom=1.46E-5
184	viscoh=2.06E-5
185	lmixmin=0.
186	yamada4_num=5
187	ELSE
188	ric=0.195
189	rifc=0.191
190	seuilalpha=1.12
191	seuilsm=0.085
192	yun=2.
193	ydeux=1.
194	hboville=.TRUE.
195	viscom=0.
196	viscoh=0.
197	lmixmin=1.
198	yamada4_num=0
199	ENDIF
200
201	PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
202	firstcall = .FALSE.
203	CALL getin_p('lmixmin',lmixmin)
204	CALL getin_p('yamada4_hboville',hboville)
205	CALL getin_p('yamada4_num',yamada4_num)
206	END IF
207
208
209
210	!===============================================================================
211	! Flags, tests et ??valuations de constantes
212	!===============================================================================
213
214	! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
215
216
217	IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN
218	STOP 'probleme de coherence dans appel a MY'
219	END IF
220
221
222	nlay = klev
223	nlev = klev + 1
224
225
226	!========================================================================
227	! Calcul des increments verticaux
228	!=========================================================================
229
230
231	! Attention: zlev n'est pas declare a nlev
232	DO ig = 1, ngrid
233	zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
234	END DO
235
236
237	DO k = 1, nlay
238	DO ig = 1, ngrid
239	unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k))
240	END DO
241	END DO
242	DO ig = 1, ngrid
243	unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1))
244	END DO
245	DO k = 2, nlay
246	DO ig = 1, ngrid
247	unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1))
248	END DO
249	END DO
250	DO ig = 1, ngrid
251	unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
252	END DO
253
254	!=========================================================================
255	! Richardson number and stability functions
256	!=========================================================================
257
258	! initialize arrays:
259
260	m2(1:ngrid, :) = 0.0
261	sm(1:ngrid, :) = 0.0
262	rif(1:ngrid, :) = 0.0
263
264	!------------------------------------------------------------
265	DO k = 2, klev
266
267	DO ig = 1, ngrid
268	dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
269	m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
270	k-1))*2)/(dz(ig,k)dz(ig,k))
271	dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
272	n2(ig, k) = g2.dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
273	ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
274	IF (ri<ric) THEN
275	rif(ig, k) = frif(ri)
276	ELSE
277	rif(ig, k) = rifc
278	END IF
279	if (new_yamada4) then
280	alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha)
281	sm(ig, k) = max(fsm(rif(ig,k)),seuilsm)
282	else
283	IF (rif(ig,k)<0.16) THEN
284	alpha(ig, k) = falpha(rif(ig,k))
285	sm(ig, k) = fsm(rif(ig,k))
286	ELSE
287	alpha(ig, k) = seuilalpha
288	sm(ig, k) = seuilsm
289	END IF
290
291	end if
292	zz(ig, k) = b1m2(ig, k)(1.-rif(ig,k))*sm(ig, k)
293	END DO
294	END DO
295
296
297
298
299
300	!=======================================================================
301	! DIFFERENT TYPES DE SCHEMA de YAMADA
302	!=======================================================================
303
304	! On commence par calculer q2 a partir de la tke
305
306	IF (new_yamada4) THEN
307	DO k=1,klev+1
308	q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux
309	ENDDO
310	ELSE
311	DO k=1,klev+1
312	q2(1:ngrid,k)=tke(1:ngrid,k)
313	ENDDO
314	ENDIF
315
316	! ====================================================================
317	! Computing the mixing length
318	! ====================================================================
319
320
321	CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
322
323
324	!--------------
325	! Yamada 2.0
326	!--------------
327	IF (iflag_pbl==6) THEN
328
329	DO k = 2, klev
330	q2(1:ngrid, k) = l(1:ngrid, k)*2zz(1:ngrid, k)
331	END DO
332
333
334	!------------------
335	! Yamada 2.Fournier
336	!------------------
337
338	ELSE IF (iflag_pbl==7) THEN
339
340
341	! Calcul de l, km, au pas precedent
342	!....................................
343	DO k = 2, klev
344	DO ig = 1, ngrid
345	delta(ig, k) = q2(ig, k)/(l(ig,k)*2sm(ig,k))
346	kmpre(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
347	mpre(ig, k) = sqrt(m2(ig,k))
348	END DO
349	END DO
350
351	DO k = 2, klev - 1
352	DO ig = 1, ngrid
353	m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12)
354	mcstat = sqrt(m2cstat)
355
356	! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3
357	!.........................................................................
358
359	IF (k==2) THEN
360	kmcstat = 1.E+0/mcstat(unsdz(ig,k)kmpre(ig,k+1)*mpre(ig,k+1)+ &
361	unsdz(ig,k-1)cd(ig)(sqrt(u(ig,3)2+v(ig,3)2)-mcstat/unsdzdec &
362	(ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &
363	-1))
364	ELSE
365	kmcstat = 1.E+0/mcstat(unsdz(ig,k)kmpre(ig,k+1)*mpre(ig,k+1)+ &
366	unsdz(ig,k-1)kmpre(ig,k-1)mpre(ig,k-1))/ &
367	(unsdz(ig,k)+unsdz(ig,k-1))
368	END IF
369
370	tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
371	q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
372
373	END DO
374	END DO
375
376
377	! ------------------------
378	! Yamada 2.5 a la Didi
379	!-------------------------
380
381	ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
382
383	! Calcul de l, km, au pas precedent
384	!....................................
385	DO k = 2, klev
386	DO ig = 1, ngrid
387	delta(ig, k) = q2(ig, k)/(l(ig,k)*2sm(ig,k))
388	IF (delta(ig,k)<1.E-20) THEN
389	delta(ig, k) = 1.E-20
390	END IF
391	km(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
392	aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
393	aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
394	aa(ig, k) = aa1dt/(delta(ig,k)l(ig,k))
395	qpre = sqrt(q2(ig,k))
396	IF (aa(ig,k)>0.) THEN
397	q2(ig, k) = (qpre+aa(ig,k)qpreqpre)**2
398	ELSE
399	q2(ig, k) = (qpre/(1.-aa(ig,k)qpre))*2
400	END IF
401	! else ! iflag_pbl=9
402	! if (aa(ig,k)*qpre.gt.0.9) then
403	! q2(ig,k)=(qpre10.)*2
404	! else
405	! q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
406	! endif
407	! endif
408	q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
409	END DO
410	END DO
411
412	ELSE IF (iflag_pbl>=10) THEN
413
414	IF (yamada4_num>=1) THEN
415
416	DO k = 2, klev - 1
417	DO ig=1,ngrid
418	q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
419	km(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
420	shear(ig,k)=km(ig, k)*m2(ig, k)
421	buoy(ig,k)=km(ig, k)m2(ig, k)(-1.*rif(ig,k))
422	dissip(ig,k)=((sqrt(q2(ig,k)))*3)/(b1l(ig,k))
423	ENDDO
424	ENDDO
425
426	IF (yamada4_num==1) THEN ! Schema du MAR tel quel
427	DO k = 2, klev - 1
428	DO ig=1,ngrid
429	tkeprov=q2(ig,k)/ydeux
430	tkeprov= tkeprov* &
431	& (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &
432	& (tkeprov+dt((-1.)min(0.,buoy(ig,k))+dissip(ig,k)))
433	q2(ig,k)=tkeprov*ydeux
434	ENDDO
435	ENDDO
436	ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation
437	DO k = 2, klev - 1
438	DO ig=1,ngrid
439	tkeprov=q2(ig,k)/ydeux
440	disseff=dissip(ig,k)-min(0.,buoy(ig,k))
441	tkeprov = tkeprov/(1.+dtdisseff/(2.tkeprov))**2
442	tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
443	q2(ig,k)=tkeprov*ydeux
444	! En cas stable, on traite la flotabilite comme la
445	! dissipation, en supposant que buoy/q2^3 est constant.
446	! Puis on prend la solution exacte
447	ENDDO
448	ENDDO
449	ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation
450	DO k = 2, klev - 1
451	DO ig=1,ngrid
452	tkeprov=q2(ig,k)/ydeux
453	disseff=dissip(ig,k)-min(0.,buoy(ig,k))
454	tkeprov=tkeprovexp(-dtdisseff/tkeprov)
455	tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
456	q2(ig,k)=tkeprov*ydeux
457	! En cas stable, on traite la flotabilite comme la
458	! dissipation, en supposant que buoy/q2^3 est constant.
459	! Puis on prend la solution exacte
460	ENDDO
461	ENDDO
462	ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation
463	DO k = 2, klev - 1
464	DO ig=1,ngrid
465	tkeprov=q2(ig,k)/ydeux
466	tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
467	tkeprov= tkeprov* &
468	& tkeprov/ &
469	& (tkeprov+dt((-1.)min(0.,buoy(ig,k))+dissip(ig,k)))
470	q2(ig,k)=tkeprov*ydeux
471	! En cas stable, on traite la flotabilite comme la
472	! dissipation, en supposant que buoy/q2^3 est constant.
473	! Puis on prend la solution exacte
474	ENDDO
475	ENDDO
476	ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
477	DO k = 2, klev - 1
478	DO ig=1,ngrid
479	tkeprov=q2(ig,k)/ydeux
480
481	! if(ifl_pbltree .eq. 0) then
482	! disseff=dissip(ig,k)-min(0.,buoy(ig,k))
483	! tkeexp=exp(-dt*disseff/tkeprov)
484	! tkeprov= shear(ig,k)tkeprov/disseff(1.-tkeexp)+tkeprov*tkeexp
485	! else
486	!FC on ajoute la dissipation due aux arbres
487	disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + drgpro(ig,k)*tkeprov
488	tkeexp=exp(-dt*disseff/tkeprov)
489	! on prend en compte la tke cree par les arbres
490	winds(ig,k)=sqrt(u(ig,k)2+v(ig,k)2)
491	tkeprov= (shear(ig,k)+ &
492	& drgpro(ig,k)(winds(ig,k))3)tkeprov/disseff(1.-tkeexp)+tkeprovtkeexp
493	! endif
494
495	q2(ig,k)=tkeprov*ydeux
496
497	! En cas stable, on traite la flotabilite comme la
498	! dissipation, en supposant que buoy/q2^3 est constant.
499	! Puis on prend la solution exacte
500	ENDDO
501	ENDDO
502	ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation
503	DO k = 2, klev - 1
504	DO ig=1,ngrid
505	tkeprov=q2(ig,k)/ydeux
506	tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
507	disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k))
508	tkeexp=exp(-dt*disseff/tkeprov)
509	tkeprov= tkeprov*tkeexp
510	q2(ig,k)=tkeprov*ydeux
511	! En cas stable, on traite la flotabilite comme la
512	! dissipation, en supposant que buoy/q2^3 est constant.
513	! Puis on prend la solution exacte
514	ENDDO
515	ENDDO
516	ENDIF
517
518	DO k = 2, klev - 1
519	DO ig=1,ngrid
520	q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
521	ENDDO
522	ENDDO
523
524	ELSE
525
526	DO k = 2, klev - 1
527	km(1:ngrid, k) = l(1:ngrid, k)sqrt(q2(1:ngrid,k))sm(1:ngrid, k)
528	q2(1:ngrid, k) = q2(1:ngrid, k) + ydeuxdtkm(1:ngrid, k)m2(1:ngrid, k)(1.-rif(1:ngrid,k))
529	! q2(1:ngrid, k) = q2(1:ngrid, k) + dtkm(1:ngrid, k)m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
530	q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
531	q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yunl(1:ngrid,k)b1))
532	! q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2l(1:ngrid,k)b1))
533	q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
534	END DO
535
536	ENDIF
537
538	ELSE
539	STOP 'Cas nom prevu dans yamada4'
540
541	END IF ! Fin du cas 8
542
543
544	! ====================================================================
545	! Calcul des coefficients de melange
546	! ====================================================================
547
548	DO k = 2, klev
549	DO ig = 1, ngrid
550	zq = sqrt(q2(ig,k))
551	km(ig, k) = l(ig, k)zqsm(ig, k) ! For momentum
552	kn(ig, k) = km(ig, k)*alpha(ig, k) ! For scalars
553	kq(ig, k) = l(ig, k)zq0.2 ! For TKE
554	END DO
555	END DO
556
557
558	!====================================================================
559	! Transport diffusif vertical de la TKE par la TKE
560	!====================================================================
561
562
563	! initialize near-surface and top-layer mixing coefficients
564	!...........................................................
565
566	kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
567	kq(1:ngrid, klev+1) = 0 ! zero at the top
568
569	! Transport diffusif vertical de la TKE.
570	!.......................................
571
572	IF (iflag_pbl>=12) THEN
573	q2(1:ngrid, 1) = q2(1:ngrid, 2)
574	CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
575	END IF
576
577
578	!====================================================================
579	! Traitement particulier pour les cas tres stables, introduction d'une
580	! longueur de m??lange minimale
581	!====================================================================
582	!
583	! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model
584	! Holtslag A.A.M. and Boville B.A.
585	! J. Clim., 6, 1825-1842, 1993
586
587
588	IF (hboville) THEN
589
590
591	IF (prt_level>1) THEN
592	PRINT *, 'YAMADA4 0'
593	END IF
594
595	DO ig = 1, ngrid
596	coriol(ig) = 1.E-4
597	pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
598	END DO
599
600	IF (1==1) THEN
601	IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
602
603	DO k = 2, klev
604	DO ig = 1, ngrid
605	IF (teta(ig,2)>teta(ig,1)) THEN
606	qmin = ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
607	kmin = kapzlev(ig, k)qmin
608	ELSE
609	kmin = -1. ! kmin n'est utilise que pour les SL stables.
610	END IF
611	IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
612
613	kn(ig, k) = kmin
614	km(ig, k) = kmin
615	kq(ig, k) = kmin
616
617	! la longueur de melange est suposee etre l= kap z
618	! K=l q Sm d'ou q2=(K/l Sm)**2
619
620	q2(ig, k) = (qmin/sm(ig,k))**2
621	END IF
622	END DO
623	END DO
624
625	ELSE
626	DO k = 2, klev
627	DO ig = 1, ngrid
628	IF (teta(ig,2)>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	END IF
634	IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
635	kn(ig, k) = kmin
636	km(ig, k) = kmin
637	kq(ig, k) = kmin
638	! la longueur de melange est suposee etre l= kap z
639	! K=l q Sm d'ou q2=(K/l Sm)**2
640	sm(ig, k) = 1.
641	alpha(ig, k) = 1.
642	q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
643	zq = sqrt(q2(ig,k))
644	km(ig, k) = l(ig, k)zqsm(ig, k)
645	kn(ig, k) = km(ig, k)*alpha(ig, k)
646	kq(ig, k) = l(ig, k)zq0.2
647	END IF
648	END DO
649	END DO
650	END IF
651
652	END IF
653
654	END IF ! hboville
655
656	! Ajout d'une viscosite moleculaire
657	km(1:ngrid,2:klev)=km(1:ngrid,2:klev)+viscom
658	kn(1:ngrid,2:klev)=kn(1:ngrid,2:klev)+viscoh
659	kq(1:ngrid,2:klev)=kq(1:ngrid,2:klev)+viscoh
660
661	IF (prt_level>1) THEN
662	PRINT *, 'YAMADA4 1'
663	END IF !(prt_level>1) THEN
664
665
666	!======================================================
667	! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
668	!======================================================
669	!
670	! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer
671	! Abdella K and McFarlane N
672	! J. Atmos. Sci., 54, 1850-1867, 1997
673
674	! Diagnostique pour stokage
675	!..........................
676
677	IF (1==0) THEN
678	rino = rif
679	smyam(1:ngrid, 1) = 0.
680	styam(1:ngrid, 1) = 0.
681	lyam(1:ngrid, 1) = 0.
682	knyam(1:ngrid, 1) = 0.
683	w2yam(1:ngrid, 1) = 0.
684	t2yam(1:ngrid, 1) = 0.
685
686	smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
687	styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
688	lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
689	knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
690
691
692	! Calcul de w'2 et T'2
693	!.......................
694
695	w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
696	lyam(1:ngrid, 2:klev)5.17kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
697	sqrt(q2(1:ngrid,2:klev))
698
699	t2yam(1:ngrid, 2:klev) = 9.1kn(1:ngrid, 2:klev) &
700	dtetadz(1:ngrid, 2:klev)*2/sqrt(q2(1:ngrid,2:klev)) &
701	lyam(1:ngrid, 2:klev)
702	END IF
703
704
705
706	!============================================================================
707	! Mise a jour de la tke
708	!============================================================================
709
710	IF (new_yamada4) THEN
711	DO k=1,klev+1
712	tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux
713	ENDDO
714	ELSE
715	DO k=1,klev+1
716	tke(1:ngrid,k)=q2(1:ngrid,k)
717	ENDDO
718	ENDIF
719
720
721	!============================================================================
722
723	RETURN
724
725
726	END SUBROUTINE yamada4
727
728	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
729
730	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
731	SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
732
733	USE dimphy
734	IMPLICIT NONE
735
736	! vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE
737	! avec un schema implicite en temps avec
738	! inversion d'un syst??me tridiagonal
739	!
740	! Reference: Description of the interface with the surface and
741	! the computation of the turbulet diffusion in LMDZ
742	! Technical note on LMDZ
743	! Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y
744	!
745	!============================================================================
746	! Declarations
747	!============================================================================
748
749	REAL plev(klon, klev+1)
750	REAL temp(klon, klev)
751	REAL timestep
752	REAL gravity, rconst
753	REAL kstar(klon, klev+1), zz
754	REAL kmy(klon, klev+1)
755	REAL q2(klon, klev+1)
756	REAL deltap(klon, klev+1)
757	REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
758	INTEGER ngrid
759
760	INTEGER i, k
761
762
763	!=========================================================================
764	! Calcul
765	!=========================================================================
766
767	DO k = 1, klev
768	DO i = 1, ngrid
769	zz = (plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
770	kstar(i, k) = 0.125(kmy(i,k+1)+kmy(i,k))zz*zz/ &
771	(plev(i,k)-plev(i,k+1))*timestep
772	END DO
773	END DO
774
775	DO k = 2, klev
776	DO i = 1, ngrid
777	deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
778	END DO
779	END DO
780	DO i = 1, ngrid
781	deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
782	deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
783	denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
784	alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
785	beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
786	END DO
787
788	DO k = klev, 2, -1
789	DO i = 1, ngrid
790	denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
791	kstar(i, k-1)
792	alpha(i, k) = (q2(i,k)deltap(i,k)+kstar(i,k)alpha(i,k+1))/denom(i, k)
793	beta(i, k) = kstar(i, k-1)/denom(i, k)
794	END DO
795	END DO
796
797	! Si on recalcule q2(1)
798	!.......................
799	IF (1==0) THEN
800	DO i = 1, ngrid
801	denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
802	q2(i, 1) = (q2(i,1)deltap(i,1)+kstar(i,1)alpha(i,2))/denom(i, 1)
803	END DO
804	END IF
805
806
807	DO k = 2, klev + 1
808	DO i = 1, ngrid
809	q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
810	END DO
811	END DO
812
813	RETURN
814	END SUBROUTINE vdif_q2
815	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
816
817
818
819	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
820	SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
821
822	USE dimphy
823	IMPLICIT NONE
824
825	! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE
826	! avec un schema explicite en temps
827
828
829	!====================================================
830	! Declarations
831	!====================================================
832
833	REAL plev(klon, klev+1)
834	REAL temp(klon, klev)
835	REAL timestep
836	REAL gravity, rconst
837	REAL kstar(klon, klev+1), zz
838	REAL kmy(klon, klev+1)
839	REAL q2(klon, klev+1)
840	REAL deltap(klon, klev+1)
841	REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
842	INTEGER ngrid
843	INTEGER i, k
844
845
846	!==================================================
847	! Calcul
848	!==================================================
849
850	DO k = 1, klev
851	DO i = 1, ngrid
852	zz = (plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
853	kstar(i, k) = 0.125(kmy(i,k+1)+kmy(i,k))zz*zz/ &
854	(plev(i,k)-plev(i,k+1))*timestep
855	END DO
856	END DO
857
858	DO k = 2, klev
859	DO i = 1, ngrid
860	deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
861	END DO
862	END DO
863	DO i = 1, ngrid
864	deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
865	deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
866	END DO
867
868	DO k = klev, 2, -1
869	DO i = 1, ngrid
870	q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
871	k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
872	END DO
873	END DO
874
875	DO i = 1, ngrid
876	q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
877	q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
878	klev)))/deltap(i, klev+1)
879	END DO
880
881	RETURN
882	END SUBROUTINE vdif_q2e
883
884	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
885
886
887	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
888
889	SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix)
890
891
892
893	USE dimphy
894	USE phys_state_var_mod, only: zstd, zsig, zmea
895	USE phys_local_var_mod, only: l_mixmin, l_mix
896
897	! zstd: ecart type de la'altitud e sous-maille
898	! zmea: altitude moyenne sous maille
899	! zsig: pente moyenne de le maille
900
901	USE geometry_mod, only: cell_area
902	! aire_cell: aire de la maille
903
904	IMPLICIT NONE
905	!*************************************************************************
906	! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence
907	! avec la formule de Blackadar
908	! Calcul d'un minimum en fonction de l'orographie sous-maille:
909	! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange
910	! induit par les circulations meso et submeso ??chelles.
911	!
912	! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere
913	! Blackadar A.K.
914	! J. Geophys. Res., 64, No 8, 1962
915	!
916	! * An evaluation of neutral and convective planetary boundary-layer parametrisations relative
917	! to large eddy simulations
918	! Ayotte K et al
919	! Boundary Layer Meteorology, 79, 131-175, 1996
920	!
921	!
922	! * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts
923	! Van de Wiel B.J.H et al
924	! Boundary-Lay Meteorol, 128, 103-166, 2008
925	!
926	!
927	! Histoire:
928	!----------
929	! * premi??re r??daction, Etienne et Frederic, 09/06/2016
930	!
931	! ***********************************************************************
932
933	!==================================================================
934	! Declarations
935	!==================================================================
936
937	! Inputs
938	!-------
939	INTEGER ni(klon) ! indice sur la grille original (non restreinte)
940	INTEGER nsrf ! Type de surface
941	INTEGER ngrid ! Nombre de points concern??s sur l'horizontal
942	INTEGER iflag_pbl ! Choix du sch??ma de turbulence
943	REAL pbl_lmixmin_alpha ! on active ou non le calcul de la longueur de melange minimum
944	REAL lmixmin ! Minimum absolu de la longueur de m??lange
945	REAL zlay(klon, klev) ! altitude du centre de la couche
946	REAL zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche
947	REAL u(klon, klev) ! vitesse du vent zonal
948	REAL v(klon, klev) ! vitesse du vent meridional
949	REAL q2(klon, klev+1) ! energie cin??tique turbulente
950	REAL n2(klon, klev+1) ! frequence de Brunt-Vaisala
951
952	!In/out
953	!-------
954
955	LOGICAL, SAVE :: firstcall = .TRUE.
956	!$OMP THREADPRIVATE(firstcall)
957
958	! Outputs
959	!---------
960
961	REAL lmix(klon, klev+1) ! Longueur de melange
962
963
964	! Local
965	!-------
966
967	INTEGER ig,jg, k
968	REAL h_oro(klon)
969	REAL hlim(klon)
970	REAL, SAVE :: kap=0.4,kapb=0.4
971	!$OMP THREADPRIVATE(kap,kapb)
972	REAL zq
973	REAL sq(klon), sqz(klon)
974	REAL, ALLOCATABLE, SAVE :: l0(:)
975	!$OMP THREADPRIVATE(l0)
976	REAL fl, zzz, zl0, zq2, zn2
977	REAL famorti, zzzz, zh_oro, zhlim
978	REAL l1(klon, klev+1), l2(klon,klev+1)
979	REAL winds(klon, klev)
980	REAL xcell
981	REAL zstdslope(klon)
982	REAL lmax
983	REAL l2strat, l2neutre, extent
984	REAL l2limit(klon)
985	!===============================================================
986	! Fonctions utiles
987	!===============================================================
988
989	! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996
990	!..........................................................................
991
992	fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)kapzlev(ig, &
993	k)/(kapzlev(ig,k)+l0(ig)),0.5sqrt(q2(ig,k))/sqrt( &
994	max(n2(ig,k),1.E-10))), 1.E-5)
995
996	! Fonction d'amortissement de la turbulence au dessus de la montagne
997	! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code
998	!.....................................................................
999
1000	famorti(zzzz, zh_oro, zhlim)=(-1.)ATAN((zzzz-zh_oro)/(zhlim-zh_oro))2./3.1416+1.
1001
1002	IF (ngrid==0) RETURN
1003
1004	IF (firstcall) THEN
1005	ALLOCATE (l0(klon))
1006	firstcall = .FALSE.
1007	END IF
1008
1009
1010	!=====================================================================
1011	! CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1
1012	!=====================================================================
1013
1014
1015	IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
1016
1017
1018	! Iterative computation of l0
1019	! This version is kept for iflag_pbl only for convergence
1020	! with NPv3.1 Cmip5 simulations
1021	!...................................................................
1022
1023	DO ig = 1, ngrid
1024	sq(ig) = 1.E-10
1025	sqz(ig) = 1.E-10
1026	END DO
1027	DO k = 2, klev - 1
1028	DO ig = 1, ngrid
1029	zq = sqrt(q2(ig,k))
1030	sqz(ig) = sqz(ig) + zqzlev(ig, k)(zlay(ig,k)-zlay(ig,k-1))
1031	sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
1032	END DO
1033	END DO
1034	DO ig = 1, ngrid
1035	l0(ig) = 0.2*sqz(ig)/sq(ig)
1036	END DO
1037	DO k = 2, klev
1038	DO ig = 1, ngrid
1039	l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
1040	END DO
1041	END DO
1042
1043	ELSE
1044
1045
1046	! In all other case, the assymptotic mixing length l0 is imposed (150m)
1047	!......................................................................
1048
1049	l0(1:ngrid) = 150.
1050	DO k = 2, klev
1051	DO ig = 1, ngrid
1052	l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
1053	END DO
1054	END DO
1055
1056	END IF
1057
1058	!=================================================================================
1059	! CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2
1060	! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les
1061	! glacier, la glace de mer et les oc??ans)
1062	!=================================================================================
1063
1064	l2(1:ngrid,:)=0.0
1065	l_mixmin(1:ngrid,:,nsrf)=0.
1066	l_mix(1:ngrid,:,nsrf)=0.
1067
1068	IF (nsrf .EQ. 1) THEN
1069
1070	! coefficients
1071	!--------------
1072
1073	extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1.
1074	lmax=150. ! Longueur de m??lange max dans l'absolu
1075
1076	! calculs
1077	!---------
1078
1079	DO ig=1,ngrid
1080
1081	! On calcule la hauteur du relief
1082	!.................................
1083	! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille
1084	! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon)
1085	! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief
1086	! (en gros, une maille de taille xcell avec une pente constante zstdslope)
1087	jg=ni(ig)
1088	! IF (zsig(jg) .EQ. 0.) THEN
1089	! zstdslope(ig)=0.
1090	! ELSE
1091	! xcell=sqrt(cell_area(jg))
1092	! zstdslope(ig)=max((xcellzsig(jg)-zmea(jg))3 /(3.zsig(jg)),0.)
1093	! zstdslope(ig)=sqrt(zstdslope(ig))
1094	! END IF
1095
1096	! h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.) ! Hauteur du relief
1097	h_oro(ig)=zstd(jg)
1098	hlim(ig)=extent*h_oro(ig)
1099	ENDDO
1100
1101	l2limit(1:ngrid)=0.
1102
1103	DO k=2,klev
1104	DO ig=1,ngrid
1105	winds(ig,k)=sqrt(u(ig,k)2+v(ig,k)2)
1106	IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie
1107	l2strat= kapbpbl_lmixmin_alphawinds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008)
1108	l2neutre=kapzlev(ig,k)h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h
1109	l2neutre=MIN(l2neutre,lmax) ! On majore par lmax
1110	l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e)
1111	l2(ig,k)=l2limit(ig)
1112
1113	ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles
1114
1115	! Au dessus des montagnes, on prend la l2limit au sommet des montagnes
1116	! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z)
1117	! et on multiplie l2limit par une fonction qui d??croit entre h et hlim
1118	l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig))
1119	ELSE ! Au dessus de extent*h, on prend l2=l0
1120	l2(ig,k)=0.
1121	END IF
1122	ENDDO
1123	ENDDO
1124	ENDIF ! pbl_lmixmin_alpha
1125
1126	!==================================================================================
1127	! On prend le max entre la longueur de melange de blackadar et celle calcul??e
1128	! en fonction de la topographie
1129	!===================================================================================
1130
1131
1132	DO k=2,klev
1133	DO ig=1,ngrid
1134	lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
1135	ENDDO
1136	ENDDO
1137
1138	! Diagnostics
1139
1140	DO k=2,klev
1141	DO ig=1,ngrid
1142	jg=ni(ig)
1143	l_mix(jg,k,nsrf)=lmix(ig,k)
1144	l_mixmin(jg,k,nsrf)=l2(ig,k)
1145	ENDDO
1146	ENDDO
1147	DO ig=1,ngrid
1148	jg=ni(ig)
1149	l_mix(jg,1,nsrf)=hlim(ig)
1150	ENDDO
1151
1152
1153
1154	END SUBROUTINE mixinglength

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

Download in other formats:

Original Format