Context Navigation

yamada4.F90 @ 2561

Last change on this file since 2561 was 2561, checked in by fhourdin, 8 years ago
Parametrisation d'une longueur de melange verticale minimum associee aux circulations meso-echelle introduites par le relief sous maille. D'apres Etienne Vignon et Frédéric Hourdin
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: 29.1 KB

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

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

Download in other formats:

Original Format