Context Navigation

coef_diff_turb_mod.F90 @ 5128

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

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

Download in other formats:

Original Format