Context Navigation

coef_diff_turb_mod.F90 @ 5229

Last change on this file since 5229 was 4881, checked in by evignon, 8 months ago
extraction plus propre de la dissipation de TKE
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.4 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 print_control_mod, 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.EQ.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.GE.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.GE.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 print_control_mod, 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	!
251	! Local variables:
252	!
253	INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite
254	!
255	! Quelques constantes et options:
256	!
257	REAL, PARAMETER :: cepdu2=0.1**2
258	REAL, PARAMETER :: CKAP=0.4
259	REAL, PARAMETER :: cb=5.0
260	REAL, PARAMETER :: cc=5.0
261	REAL, PARAMETER :: cd=5.0
262	REAL, PARAMETER :: clam=160.0
263	REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau
264	LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide
265	REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
266	REAL, PARAMETER :: prandtl=0.4
267	REAL kstable ! diffusion minimale (situation stable)
268	! GKtest
269	! PARAMETER (kstable=1.0e-10)
270	!IM: 261103 REAL kstable_ter, kstable_sinon
271	!IM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
272	!IM: 261103 PARAMETER (kstable_ter = 1.0e-8)
273	!IM: 261103 PARAMETER (kstable_ter = 1.0e-10)
274	!IM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
275	! fin GKtest
276	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
277	INTEGER isommet ! le sommet de la couche limite
278	LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
279	LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere
280
281	!
282	! Variables locales:
283	INTEGER i, k !IM 120704
284	REAL zgeop(klon,klev)
285	REAL zmgeom(klon)
286	REAL zri(klon)
287	REAL zl2(klon)
288	REAL zdphi, zdu2, ztvd, ztvu, zcdn
289	REAL zscf
290	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
291	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
292	REAL, PARAMETER :: t_coup=273.15
293	LOGICAL, PARAMETER :: check=.FALSE.
294	!
295	! contre-gradient pour la chaleur sensible: Kelvin/metre
296	REAL gamt(2:klev)
297
298	LOGICAL, SAVE :: appel1er=.TRUE.
299	!$OMP THREADPRIVATE(appel1er)
300	!
301	! Fonctions thermodynamiques et fonctions d'instabilite
302	REAL fsta, fins, x
303
304	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
305	fins(x) = SQRT(1.0-18.0*x)
306
307	isommet=klev
308
309	IF (appel1er) THEN
310	IF (prt_level > 9) THEN
311	WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
312	WRITE(lunout,*)'coefkz, richum:', richum
313	IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
314	WRITE(lunout,*)'coefkz, isommet:', isommet
315	WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
316	appel1er = .FALSE.
317	ENDIF
318	ENDIF
319	!
320	! Initialiser les sorties
321	!
322	DO k = 1, klev
323	DO i = 1, knon
324	pcfm(i,k) = 0.0
325	pcfh(i,k) = 0.0
326	ENDDO
327	ENDDO
328	DO i = 1, knon
329	itop(i) = 0
330	ENDDO
331
332	!
333	! Prescrire la valeur de contre-gradient
334	!
335	IF (iflag_pbl.EQ.1) THEN
336	DO k = 3, klev
337	gamt(k) = -1.0E-03
338	ENDDO
339	gamt(2) = -2.5E-03
340	ELSE
341	DO k = 2, klev
342	gamt(k) = 0.0
343	ENDDO
344	ENDIF
345	!IM cf JLD/ GKtest
346	IF ( nsrf .NE. is_oce ) THEN
347	!IM 261103 kstable = kstable_ter
348	kstable = ksta_ter
349	ELSE
350	!IM 261103 kstable = kstable_sinon
351	kstable = ksta
352	ENDIF
353	!IM cf JLD/ GKtest fin
354
355	!
356	! Calculer les geopotentiels de chaque couche
357	!
358	DO i = 1, knon
359	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
360	* (paprs(i,1)-pplay(i,1))
361	ENDDO
362	DO k = 2, klev
363	DO i = 1, knon
364	zgeop(i,k) = zgeop(i,k-1) &
365	+ RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) &
366	* (pplay(i,k-1)-pplay(i,k))
367	ENDDO
368	ENDDO
369
370	!
371	! Calculer les coefficients turbulents dans l'atmosphere
372	!
373	DO i = 1, knon
374	itop(i) = isommet
375	ENDDO
376
377
378	DO k = 2, isommet
379	DO i = 1, knon
380	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 &
381	+(v(i,k)-v(i,k-1))**2)
382	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
383	zdphi =zmgeom(i) / 2.0
384	zt = (t(i,k)+t(i,k-1)) * 0.5
385	zq = (q(i,k)+q(i,k-1)) * 0.5
386
387	!
388	! Calculer Qs et dQs/dT:
389	!
390	IF (thermcep) THEN
391	zdelta = MAX(0.,SIGN(1.,RTT-zt))
392	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta) &
393	+ R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
394	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
395	zqs = MIN(0.5,zqs)
396	zcor = 1./(1.-RETV*zqs)
397	zqs = zqs*zcor
398	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
399	ELSE
400	IF (zt .LT. t_coup) THEN
401	zqs = qsats(zt) / pplay(i,k)
402	zdqs = dqsats(zt,zqs)
403	ELSE
404	zqs = qsatl(zt) / pplay(i,k)
405	zdqs = dqsatl(zt,zqs)
406	ENDIF
407	ENDIF
408	!
409	! calculer la fraction nuageuse (processus humide):
410	!
411	if (zq /= 0.) then
412	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
413	else
414	zfr = 0.
415	end if
416	zfr = MAX(0.0,MIN(1.0,zfr))
417	IF (.NOT.richum) zfr = 0.0
418	!
419	! calculer le nombre de Richardson:
420	!
421	IF (tvirtu) THEN
422	ztvd =( t(i,k) &
423	+ zdphi/RCPD/(1.+RVTMP2*zq) &
424	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
425	)(1.+RETVq(i,k))
426	ztvu =( t(i,k-1) &
427	- zdphi/RCPD/(1.+RVTMP2*zq) &
428	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
429	)(1.+RETVq(i,k-1))
430	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
431	zri(i) = zri(i) &
432	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
433	(paprs(i,k)/101325.0)*RKAPPA &
434	/(zdu20.5(ztvd+ztvu))
435
436	ELSE ! calcul de Ridchardson compatible LMD5
437
438	zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
439	-RD0.5(t(i,k)+t(i,k-1))/paprs(i,k) &
440	*(pplay(i,k)-pplay(i,k-1)) &
441	)zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
442	zri(i) = zri(i) + &
443	zmgeom(i)zmgeom(i)gamt(k)/RG &
444	(paprs(i,k)/101325.0)*RKAPPA &
445	/(zdu20.5(t(i,k-1)+t(i,k)))
446	ENDIF
447	!
448	! finalement, les coefficients d'echange sont obtenus:
449	!
450	zcdn=SQRT(zdu2) / zmgeom(i) * RG
451
452	IF (opt_ec) THEN
453	z2geomf=zgeop(i,k-1)+zgeop(i,k)
454	zalm2=(0.5ckap/RGz2geomf &
455	/(1.+0.5ckap/rg/clamz2geomf))**2
456	zalh2=(0.5ckap/rgz2geomf &
457	/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
458	IF (zri(i).LT.0.0) THEN ! situation instable
459	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3 &
460	/ (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
461	zscf = SQRT(-zri(i)*zscf)
462	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
463	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
464	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
465	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
466	ELSE ! situation stable
467	zscf=SQRT(1.+cd*zri(i))
468	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
469	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
470	ENDIF
471	ELSE
472	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
473	/(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
474	pcfm(i,k)=SQRT(MAX(zcdnzcdn(ric-zri(i))/ric, kstable))
475	pcfm(i,k)= zl2(i)* pcfm(i,k)
476	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
477	ENDIF
478	ENDDO
479	ENDDO
480
481	!
482	! Au-dela du sommet, pas de diffusion turbulente:
483	!
484	DO i = 1, knon
485	IF (itop(i)+1 .LE. klev) THEN
486	DO k = itop(i)+1, klev
487	pcfh(i,k) = 0.0
488	pcfm(i,k) = 0.0
489	ENDDO
490	ENDIF
491	ENDDO
492
493	END SUBROUTINE coefkz
494	!
495	!****************************************************************************************
496	!
497	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
498	pcfm, pcfh)
499
500	USE dimphy
501	USE indice_sol_mod
502
503	!======================================================================
504	! J'introduit un peu de diffusion sauf dans les endroits
505	! ou une forte inversion est presente
506	! On peut dire qu'il represente la convection peu profonde
507	!
508	! Arguments:
509	! nsrf-----input-I- indicateur de la nature du sol
510	! knon-----input-I- nombre de points a traiter
511	! paprs----input-R- pression a chaque intercouche (en Pa)
512	! pplay----input-R- pression au milieu de chaque couche (en Pa)
513	! t--------input-R- temperature (K)
514	!
515	! pcfm-----output-R- coefficients a calculer (vitesse)
516	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
517	!======================================================================
518	!
519	! Arguments:
520	!
521	INTEGER, INTENT(IN) :: knon, nsrf
522	REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
523	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
524	REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
525
526	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
527	!
528	! Quelques constantes et options:
529	!
530	REAL, PARAMETER :: prandtl=0.4
531	REAL, PARAMETER :: kstable=0.002
532	! REAL, PARAMETER :: kstable=0.001
533	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
534	REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
535	! PARAMETER (seuil=-0.04)
536	! PARAMETER (seuil=-0.06)
537	! PARAMETER (seuil=-0.09)
538
539	!
540	! Variables locales:
541	!
542	INTEGER i, k, invb(knon)
543	REAL zl2(knon)
544	REAL zdthmin(knon), zdthdp
545
546	INCLUDE "YOMCST.h"
547	!
548	! Initialiser les sorties
549	!
550	DO k = 1, klev
551	DO i = 1, knon
552	pcfm(i,k) = 0.0
553	pcfh(i,k) = 0.0
554	ENDDO
555	ENDDO
556
557	!
558	! Chercher la zone d'inversion forte
559	!
560	DO i = 1, knon
561	invb(i) = klev
562	zdthmin(i)=0.0
563	ENDDO
564	DO k = 2, klev/2-1
565	DO i = 1, knon
566	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
567	- RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
568	zdthdp = zdthdp * 100.0
569	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
570	zdthdp.LT.zdthmin(i) ) THEN
571	zdthmin(i) = zdthdp
572	invb(i) = k
573	ENDIF
574	ENDDO
575	ENDDO
576
577	!
578	! Introduire une diffusion:
579	!
580	IF ( nsrf.EQ.is_oce ) THEN
581	DO k = 2, klev
582	DO i = 1, knon
583	!IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
584	!IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
585	!IM cf JLD/ GKtest TERkz2
586	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
587	! fin GKtest
588
589
590	! s'il n'y a pas d'inversion ou si l'inversion est trop faible
591	! IF ( (nsrf.EQ.is_oce) .AND. &
592	IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN
593	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
594	/(paprs(i,2)-paprs(i,klev+1)) ))**2
595	pcfm(i,k)= zl2(i)* kstable
596	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
597	ENDIF
598	ENDDO
599	ENDDO
600	ENDIF
601
602	END SUBROUTINE coefkz2
603	!
604	!****************************************************************************************
605	!
606	END MODULE coef_diff_turb_mod

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

Download in other formats:

Original Format