Context Navigation

coef_diff_turb_mod.F90 @ 1907

Last change on this file since 1907 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
Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 18.8 KB

Line
1	!
2	MODULE coef_diff_turb_mod
3	!
4	! This module contains some procedures for calculation of the coefficients of the
5	! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion
6	! at surface(cdrag)
7	!
8	IMPLICIT NONE
9
10	CONTAINS
11	!
12	!****************************************************************************************
13	!
14	SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
15	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
16	ycoefm, ycoefh ,yq2)
17
18	USE dimphy
19	USE indice_sol_mod
20	!
21	! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
22	! atmosphere
23	! NB! No values are calculated between surface and the first model layer.
24	! ycoefm(:,1) and ycoefh(:,1) are not valid !!!
25	!
26	!
27	! Input arguments
28	!****************************************************************************************
29	REAL, INTENT(IN) :: dtime
30	INTEGER, INTENT(IN) :: nsrf, knon
31	INTEGER, DIMENSION(klon), INTENT(IN) :: ni
32	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs
33	REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay
34	REAL, DIMENSION(klon,klev), INTENT(IN) :: yu, yv
35	REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt
36	REAL, DIMENSION(klon), INTENT(IN) :: yts, yrugos, yqsurf
37	REAL, DIMENSION(klon), INTENT(IN) :: ycdragm
38
39	! InOutput arguments
40	!****************************************************************************************
41	REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2
42
43	! Output arguments
44	!****************************************************************************************
45	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefh
46	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm
47
48	! Other local variables
49	!****************************************************************************************
50	INTEGER :: k, i, j
51	REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta
52	REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq
53	REAL, DIMENSION(klon) :: yustar
54
55	! Include
56	!****************************************************************************************
57	INCLUDE "clesphys.h"
58	INCLUDE "iniprint.h"
59	INCLUDE "compbl.h"
60	INCLUDE "YOETHF.h"
61	INCLUDE "YOMCST.h"
62
63
64	!****************************************************************************************
65	! Calcul de coefficients de diffusion turbulent de l'atmosphere :
66	! ycoefm(:,2:klev), ycoefh(:,2:klev)
67	!
68	!****************************************************************************************
69
70	CALL coefkz(nsrf, knon, ypaprs, ypplay, &
71	ksta, ksta_ter, &
72	yts, yrugos, yu, yv, yt, yq, &
73	yqsurf, &
74	ycoefm, ycoefh)
75
76	!****************************************************************************************
77	! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere :
78	! ycoefm(:,2:klev), ycoefh(:,2:klev)
79	!
80	!****************************************************************************************
81
82	IF (iflag_pbl.EQ.1) THEN
83	CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, &
84	ycoefm0, ycoefh0)
85
86	DO k = 2, klev
87	DO i = 1, knon
88	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
89	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
90	ENDDO
91	ENDDO
92	ENDIF
93
94
95	!****************************************************************************************
96	! Calcul d'une diffusion minimale pour les conditions tres stables
97	!
98	!****************************************************************************************
99	IF (ok_kzmin) THEN
100	CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycdragm, &
101	ycoefm0,ycoefh0)
102
103	DO k = 2, klev
104	DO i = 1, knon
105	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
106	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
107	ENDDO
108	ENDDO
109
110	ENDIF
111
112
113	!****************************************************************************************
114	! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
115	!
116	!****************************************************************************************
117
118	IF (iflag_pbl.GE.3) THEN
119
120	yzlay(1:knon,1)= &
121	RDyt(1:knon,1)/(0.5(ypaprs(1:knon,1)+ypplay(1:knon,1))) &
122	*(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
123	DO k=2,klev
124	DO i = 1, knon
125	yzlay(i,k)= &
126	yzlay(i,k-1)+RD0.5(yt(i,k-1)+yt(i,k)) &
127	/ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG
128	END DO
129	END DO
130
131	DO k=1,klev
132	DO i = 1, knon
133	yteta(i,k)= &
134	yt(i,k)(ypaprs(i,1)/ypplay(i,k))*RKAPPA &
135	(1.+0.61yq(i,k))
136	END DO
137	END DO
138
139	yzlev(1:knon,1)=0.
140	yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
141	DO k=2,klev
142	DO i = 1, knon
143	yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1))
144	END DO
145	END DO
146
147	!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
148	!!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un
149	!!$! bug sur les coefficients de surface :
150	!!$! ycdragh(1:knon) = ycoefm(1:knon,1)
151	!!$! ycdragm(1:knon) = ycoefh(1:knon,1)
152	!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
153	CALL ustarhb(knon,yu,yv,ycdragm, yustar)
154
155	IF (prt_level > 9) THEN
156	WRITE(lunout,*) 'USTAR = ',yustar
157	ENDIF
158
159	! iflag_pbl peut etre utilise comme longuer de melange
160	IF (iflag_pbl.GE.31) THEN
161	CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
162	yzlev,yzlay,yu,yv,yteta, &
163	ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &
164	iflag_pbl)
165	ELSE IF (iflag_pbl<20) THEN
166	CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, &
167	yzlev,yzlay,yu,yv,yteta, &
168	ycdragm,yq2,ykmm,ykmn,ykmq,yustar, &
169	iflag_pbl)
170	ENDIF
171
172	ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
173	ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
174
175	ENDIF !(iflag_pbl.ge.3)
176
177	END SUBROUTINE coef_diff_turb
178	!
179	!****************************************************************************************
180	!
181	SUBROUTINE coefkz(nsrf, knon, paprs, pplay, &
182	ksta, ksta_ter, &
183	ts, rugos, &
184	u,v,t,q, &
185	qsurf, &
186	pcfm, pcfh)
187
188	USE dimphy
189	USE indice_sol_mod
190
191	!======================================================================
192	! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
193	! (une version strictement identique a l'ancien modele)
194	! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
195	! coefficients d'echange turbulent dans l'atmosphere.
196	! Arguments:
197	! nsrf-----input-I- indicateur de la nature du sol
198	! knon-----input-I- nombre de points a traiter
199	! paprs----input-R- pregssion a chaque intercouche (en Pa)
200	! pplay----input-R- pression au milieu de chaque couche (en Pa)
201	! ts-------input-R- temperature du sol (en Kelvin)
202	! rugos----input-R- longeur de rugosite (en m)
203	! u--------input-R- vitesse u
204	! v--------input-R- vitesse v
205	! t--------input-R- temperature (K)
206	! q--------input-R- vapeur d'eau (kg/kg)
207	!
208	! pcfm-----output-R- coefficients a calculer (vitesse)
209	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
210	!======================================================================
211	INCLUDE "YOETHF.h"
212	INCLUDE "FCTTRE.h"
213	INCLUDE "iniprint.h"
214	INCLUDE "compbl.h"
215	INCLUDE "YOMCST.h"
216	!
217	! Arguments:
218	!
219	INTEGER, INTENT(IN) :: knon, nsrf
220	REAL, INTENT(IN) :: ksta, ksta_ter
221	REAL, DIMENSION(klon), INTENT(IN) :: ts
222	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
223	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
224	REAL, DIMENSION(klon,klev), INTENT(IN) :: u, v, t, q
225	REAL, DIMENSION(klon), INTENT(IN) :: rugos
226	REAL, DIMENSION(klon), INTENT(IN) :: qsurf
227
228	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pcfm, pcfh
229
230	!
231	! Local variables:
232	!
233	INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite
234	!
235	! Quelques constantes et options:
236	!
237	REAL, PARAMETER :: cepdu2=0.1**2
238	REAL, PARAMETER :: CKAP=0.4
239	REAL, PARAMETER :: cb=5.0
240	REAL, PARAMETER :: cc=5.0
241	REAL, PARAMETER :: cd=5.0
242	REAL, PARAMETER :: clam=160.0
243	REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau
244	LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide
245	REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
246	REAL, PARAMETER :: prandtl=0.4
247	REAL kstable ! diffusion minimale (situation stable)
248	! GKtest
249	! PARAMETER (kstable=1.0e-10)
250	!IM: 261103 REAL kstable_ter, kstable_sinon
251	!IM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
252	!IM: 261103 PARAMETER (kstable_ter = 1.0e-8)
253	!IM: 261103 PARAMETER (kstable_ter = 1.0e-10)
254	!IM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
255	! fin GKtest
256	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
257	INTEGER isommet ! le sommet de la couche limite
258	LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
259	LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere
260
261	!
262	! Variables locales:
263	INTEGER i, k !IM 120704
264	REAL zgeop(klon,klev)
265	REAL zmgeom(klon)
266	REAL zri(klon)
267	REAL zl2(klon)
268	REAL zdphi, zdu2, ztvd, ztvu, zcdn
269	REAL zscf
270	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
271	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
272	REAL, PARAMETER :: t_coup=273.15
273	LOGICAL, PARAMETER :: check=.FALSE.
274	!
275	! contre-gradient pour la chaleur sensible: Kelvin/metre
276	REAL gamt(2:klev)
277
278	LOGICAL, SAVE :: appel1er=.TRUE.
279	!$OMP THREADPRIVATE(appel1er)
280	!
281	! Fonctions thermodynamiques et fonctions d'instabilite
282	REAL fsta, fins, x
283
284	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
285	fins(x) = SQRT(1.0-18.0*x)
286
287	isommet=klev
288
289	IF (appel1er) THEN
290	IF (prt_level > 9) THEN
291	WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
292	WRITE(lunout,*)'coefkz, richum:', richum
293	IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
294	WRITE(lunout,*)'coefkz, isommet:', isommet
295	WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
296	appel1er = .FALSE.
297	ENDIF
298	ENDIF
299	!
300	! Initialiser les sorties
301	!
302	DO k = 1, klev
303	DO i = 1, knon
304	pcfm(i,k) = 0.0
305	pcfh(i,k) = 0.0
306	ENDDO
307	ENDDO
308	DO i = 1, knon
309	itop(i) = 0
310	ENDDO
311
312	!
313	! Prescrire la valeur de contre-gradient
314	!
315	IF (iflag_pbl.EQ.1) THEN
316	DO k = 3, klev
317	gamt(k) = -1.0E-03
318	ENDDO
319	gamt(2) = -2.5E-03
320	ELSE
321	DO k = 2, klev
322	gamt(k) = 0.0
323	ENDDO
324	ENDIF
325	!IM cf JLD/ GKtest
326	IF ( nsrf .NE. is_oce ) THEN
327	!IM 261103 kstable = kstable_ter
328	kstable = ksta_ter
329	ELSE
330	!IM 261103 kstable = kstable_sinon
331	kstable = ksta
332	ENDIF
333	!IM cf JLD/ GKtest fin
334
335	!
336	! Calculer les geopotentiels de chaque couche
337	!
338	DO i = 1, knon
339	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
340	* (paprs(i,1)-pplay(i,1))
341	ENDDO
342	DO k = 2, klev
343	DO i = 1, knon
344	zgeop(i,k) = zgeop(i,k-1) &
345	+ RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) &
346	* (pplay(i,k-1)-pplay(i,k))
347	ENDDO
348	ENDDO
349
350	!
351	! Calculer les coefficients turbulents dans l'atmosphere
352	!
353	DO i = 1, knon
354	itop(i) = isommet
355	ENDDO
356
357
358	DO k = 2, isommet
359	DO i = 1, knon
360	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 &
361	+(v(i,k)-v(i,k-1))**2)
362	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
363	zdphi =zmgeom(i) / 2.0
364	zt = (t(i,k)+t(i,k-1)) * 0.5
365	zq = (q(i,k)+q(i,k-1)) * 0.5
366
367	!
368	! Calculer Qs et dQs/dT:
369	!
370	IF (thermcep) THEN
371	zdelta = MAX(0.,SIGN(1.,RTT-zt))
372	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta) &
373	+ R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
374	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
375	zqs = MIN(0.5,zqs)
376	zcor = 1./(1.-RETV*zqs)
377	zqs = zqs*zcor
378	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
379	ELSE
380	IF (zt .LT. t_coup) THEN
381	zqs = qsats(zt) / pplay(i,k)
382	zdqs = dqsats(zt,zqs)
383	ELSE
384	zqs = qsatl(zt) / pplay(i,k)
385	zdqs = dqsatl(zt,zqs)
386	ENDIF
387	ENDIF
388	!
389	! calculer la fraction nuageuse (processus humide):
390	!
391	if (zq /= 0.) then
392	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
393	else
394	zfr = 0.
395	end if
396	zfr = MAX(0.0,MIN(1.0,zfr))
397	IF (.NOT.richum) zfr = 0.0
398	!
399	! calculer le nombre de Richardson:
400	!
401	IF (tvirtu) THEN
402	ztvd =( t(i,k) &
403	+ zdphi/RCPD/(1.+RVTMP2*zq) &
404	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
405	)(1.+RETVq(i,k))
406	ztvu =( t(i,k-1) &
407	- zdphi/RCPD/(1.+RVTMP2*zq) &
408	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
409	)(1.+RETVq(i,k-1))
410	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
411	zri(i) = zri(i) &
412	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
413	(paprs(i,k)/101325.0)*RKAPPA &
414	/(zdu20.5(ztvd+ztvu))
415
416	ELSE ! calcul de Ridchardson compatible LMD5
417
418	zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
419	-RD0.5(t(i,k)+t(i,k-1))/paprs(i,k) &
420	*(pplay(i,k)-pplay(i,k-1)) &
421	)zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
422	zri(i) = zri(i) + &
423	zmgeom(i)zmgeom(i)gamt(k)/RG &
424	(paprs(i,k)/101325.0)*RKAPPA &
425	/(zdu20.5(t(i,k-1)+t(i,k)))
426	ENDIF
427	!
428	! finalement, les coefficients d'echange sont obtenus:
429	!
430	zcdn=SQRT(zdu2) / zmgeom(i) * RG
431
432	IF (opt_ec) THEN
433	z2geomf=zgeop(i,k-1)+zgeop(i,k)
434	zalm2=(0.5ckap/RGz2geomf &
435	/(1.+0.5ckap/rg/clamz2geomf))**2
436	zalh2=(0.5ckap/rgz2geomf &
437	/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
438	IF (zri(i).LT.0.0) THEN ! situation instable
439	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3 &
440	/ (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
441	zscf = SQRT(-zri(i)*zscf)
442	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
443	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
444	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
445	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
446	ELSE ! situation stable
447	zscf=SQRT(1.+cd*zri(i))
448	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
449	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
450	ENDIF
451	ELSE
452	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
453	/(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
454	pcfm(i,k)=SQRT(MAX(zcdnzcdn(ric-zri(i))/ric, kstable))
455	pcfm(i,k)= zl2(i)* pcfm(i,k)
456	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
457	ENDIF
458	ENDDO
459	ENDDO
460
461	!
462	! Au-dela du sommet, pas de diffusion turbulente:
463	!
464	DO i = 1, knon
465	IF (itop(i)+1 .LE. klev) THEN
466	DO k = itop(i)+1, klev
467	pcfh(i,k) = 0.0
468	pcfm(i,k) = 0.0
469	ENDDO
470	ENDIF
471	ENDDO
472
473	END SUBROUTINE coefkz
474	!
475	!****************************************************************************************
476	!
477	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
478	pcfm, pcfh)
479
480	USE dimphy
481	USE indice_sol_mod
482
483	!======================================================================
484	! J'introduit un peu de diffusion sauf dans les endroits
485	! ou une forte inversion est presente
486	! On peut dire qu'il represente la convection peu profonde
487	!
488	! Arguments:
489	! nsrf-----input-I- indicateur de la nature du sol
490	! knon-----input-I- nombre de points a traiter
491	! paprs----input-R- pression a chaque intercouche (en Pa)
492	! pplay----input-R- pression au milieu de chaque couche (en Pa)
493	! t--------input-R- temperature (K)
494	!
495	! pcfm-----output-R- coefficients a calculer (vitesse)
496	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
497	!======================================================================
498	!
499	! Arguments:
500	!
501	INTEGER, INTENT(IN) :: knon, nsrf
502	REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
503	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
504	REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
505
506	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
507	!
508	! Quelques constantes et options:
509	!
510	REAL, PARAMETER :: prandtl=0.4
511	REAL, PARAMETER :: kstable=0.002
512	! REAL, PARAMETER :: kstable=0.001
513	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
514	REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
515	! PARAMETER (seuil=-0.04)
516	! PARAMETER (seuil=-0.06)
517	! PARAMETER (seuil=-0.09)
518
519	!
520	! Variables locales:
521	!
522	INTEGER i, k, invb(knon)
523	REAL zl2(knon)
524	REAL zdthmin(knon), zdthdp
525
526	INCLUDE "YOMCST.h"
527	!
528	! Initialiser les sorties
529	!
530	DO k = 1, klev
531	DO i = 1, knon
532	pcfm(i,k) = 0.0
533	pcfh(i,k) = 0.0
534	ENDDO
535	ENDDO
536
537	!
538	! Chercher la zone d'inversion forte
539	!
540	DO i = 1, knon
541	invb(i) = klev
542	zdthmin(i)=0.0
543	ENDDO
544	DO k = 2, klev/2-1
545	DO i = 1, knon
546	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
547	- RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
548	zdthdp = zdthdp * 100.0
549	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
550	zdthdp.LT.zdthmin(i) ) THEN
551	zdthmin(i) = zdthdp
552	invb(i) = k
553	ENDIF
554	ENDDO
555	ENDDO
556
557	!
558	! Introduire une diffusion:
559	!
560	IF ( nsrf.EQ.is_oce ) THEN
561	DO k = 2, klev
562	DO i = 1, knon
563	!IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
564	!IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
565	!IM cf JLD/ GKtest TERkz2
566	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
567	! fin GKtest
568
569
570	! s'il n'y a pas d'inversion ou si l'inversion est trop faible
571	! IF ( (nsrf.EQ.is_oce) .AND. &
572	IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN
573	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
574	/(paprs(i,2)-paprs(i,klev+1)) ))**2
575	pcfm(i,k)= zl2(i)* kstable
576	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
577	ENDIF
578	ENDDO
579	ENDDO
580	ENDIF
581
582	END SUBROUTINE coefkz2
583	!
584	!****************************************************************************************
585	!
586	END MODULE coef_diff_turb_mod

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

Download in other formats:

Original Format