Context Navigation

coef_diff_turb_mod.F90 @ 1067

Last change on this file since 1067 was 1067, checked in by Laurent Fairhead, 15 years ago

Modifications lie au premier niveau du modele pour la diffusion turbulent

du vent.

Preparation pour un couplage des courrant oceaniques.

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 18.7 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	!
20	! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
21	! atmosphere
22	! NB! No values are calculated between surface and the first model layer.
23	! ycoefm(:,1) and ycoefh(:,1) are not valid !!!
24	!
25	!
26	! Input arguments
27	!****************************************************************************************
28	REAL, INTENT(IN) :: dtime
29	INTEGER, INTENT(IN) :: nsrf, knon
30	INTEGER, DIMENSION(klon), INTENT(IN) :: ni
31	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs
32	REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay
33	REAL, DIMENSION(klon,klev), INTENT(IN) :: yu, yv
34	REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt
35	REAL, DIMENSION(klon), INTENT(IN) :: yts, yrugos, yqsurf
36	REAL, DIMENSION(klon), INTENT(IN) :: ycdragm
37
38	! InOutput arguments
39	!****************************************************************************************
40	REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2
41
42	! Output arguments
43	!****************************************************************************************
44	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefh
45	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm
46
47	! Other local variables
48	!****************************************************************************************
49	INTEGER :: k, i, j
50	REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta
51	REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq
52	REAL, DIMENSION(klon) :: yustar
53
54	! Include
55	!****************************************************************************************
56	INCLUDE "clesphys.h"
57	INCLUDE "indicesol.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.11) 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
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
190	!======================================================================
191	! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
192	! (une version strictement identique a l'ancien modele)
193	! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
194	! coefficients d'echange turbulent dans l'atmosphere.
195	! Arguments:
196	! nsrf-----input-I- indicateur de la nature du sol
197	! knon-----input-I- nombre de points a traiter
198	! paprs----input-R- pregssion a chaque intercouche (en Pa)
199	! pplay----input-R- pression au milieu de chaque couche (en Pa)
200	! ts-------input-R- temperature du sol (en Kelvin)
201	! rugos----input-R- longeur de rugosite (en m)
202	! u--------input-R- vitesse u
203	! v--------input-R- vitesse v
204	! t--------input-R- temperature (K)
205	! q--------input-R- vapeur d'eau (kg/kg)
206	!
207	! pcfm-----output-R- coefficients a calculer (vitesse)
208	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
209	!======================================================================
210	INCLUDE "YOETHF.h"
211	INCLUDE "FCTTRE.h"
212	INCLUDE "iniprint.h"
213	INCLUDE "indicesol.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	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
392	zfr = MAX(0.0,MIN(1.0,zfr))
393	IF (.NOT.richum) zfr = 0.0
394	!
395	! calculer le nombre de Richardson:
396	!
397	IF (tvirtu) THEN
398	ztvd =( t(i,k) &
399	+ zdphi/RCPD/(1.+RVTMP2*zq) &
400	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
401	)(1.+RETVq(i,k))
402	ztvu =( t(i,k-1) &
403	- zdphi/RCPD/(1.+RVTMP2*zq) &
404	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
405	)(1.+RETVq(i,k-1))
406	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
407	zri(i) = zri(i) &
408	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
409	(paprs(i,k)/101325.0)*RKAPPA &
410	/(zdu20.5(ztvd+ztvu))
411
412	ELSE ! calcul de Ridchardson compatible LMD5
413
414	zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
415	-RD0.5(t(i,k)+t(i,k-1))/paprs(i,k) &
416	*(pplay(i,k)-pplay(i,k-1)) &
417	)zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
418	zri(i) = zri(i) + &
419	zmgeom(i)zmgeom(i)gamt(k)/RG &
420	(paprs(i,k)/101325.0)*RKAPPA &
421	/(zdu20.5(t(i,k-1)+t(i,k)))
422	ENDIF
423	!
424	! finalement, les coefficients d'echange sont obtenus:
425	!
426	zcdn=SQRT(zdu2) / zmgeom(i) * RG
427
428	IF (opt_ec) THEN
429	z2geomf=zgeop(i,k-1)+zgeop(i,k)
430	zalm2=(0.5ckap/RGz2geomf &
431	/(1.+0.5ckap/rg/clamz2geomf))**2
432	zalh2=(0.5ckap/rgz2geomf &
433	/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
434	IF (zri(i).LT.0.0) THEN ! situation instable
435	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3 &
436	/ (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
437	zscf = SQRT(-zri(i)*zscf)
438	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
439	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
440	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
441	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
442	ELSE ! situation stable
443	zscf=SQRT(1.+cd*zri(i))
444	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
445	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
446	ENDIF
447	ELSE
448	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
449	/(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
450	pcfm(i,k)=SQRT(MAX(zcdnzcdn(ric-zri(i))/ric, kstable))
451	pcfm(i,k)= zl2(i)* pcfm(i,k)
452	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
453	ENDIF
454	ENDDO
455	ENDDO
456
457	!
458	! Au-dela du sommet, pas de diffusion turbulente:
459	!
460	DO i = 1, knon
461	IF (itop(i)+1 .LE. klev) THEN
462	DO k = itop(i)+1, klev
463	pcfh(i,k) = 0.0
464	pcfm(i,k) = 0.0
465	ENDDO
466	ENDIF
467	ENDDO
468
469	END SUBROUTINE coefkz
470	!
471	!****************************************************************************************
472	!
473	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
474	pcfm, pcfh)
475
476	USE dimphy
477
478	!======================================================================
479	! J'introduit un peu de diffusion sauf dans les endroits
480	! ou une forte inversion est presente
481	! On peut dire qu'il represente la convection peu profonde
482	!
483	! Arguments:
484	! nsrf-----input-I- indicateur de la nature du sol
485	! knon-----input-I- nombre de points a traiter
486	! paprs----input-R- pression a chaque intercouche (en Pa)
487	! pplay----input-R- pression au milieu de chaque couche (en Pa)
488	! t--------input-R- temperature (K)
489	!
490	! pcfm-----output-R- coefficients a calculer (vitesse)
491	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
492	!======================================================================
493	!
494	! Arguments:
495	!
496	INTEGER, INTENT(IN) :: knon, nsrf
497	REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
498	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
499	REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
500
501	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
502	!
503	! Quelques constantes et options:
504	!
505	REAL, PARAMETER :: prandtl=0.4
506	REAL, PARAMETER :: kstable=0.002
507	! REAL, PARAMETER :: kstable=0.001
508	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
509	REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
510	! PARAMETER (seuil=-0.04)
511	! PARAMETER (seuil=-0.06)
512	! PARAMETER (seuil=-0.09)
513
514	!
515	! Variables locales:
516	!
517	INTEGER i, k, invb(knon)
518	REAL zl2(knon)
519	REAL zdthmin(knon), zdthdp
520
521	INCLUDE "indicesol.h"
522	INCLUDE "YOMCST.h"
523	!
524	! Initialiser les sorties
525	!
526	DO k = 1, klev
527	DO i = 1, knon
528	pcfm(i,k) = 0.0
529	pcfh(i,k) = 0.0
530	ENDDO
531	ENDDO
532
533	!
534	! Chercher la zone d'inversion forte
535	!
536	DO i = 1, knon
537	invb(i) = klev
538	zdthmin(i)=0.0
539	ENDDO
540	DO k = 2, klev/2-1
541	DO i = 1, knon
542	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
543	- RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
544	zdthdp = zdthdp * 100.0
545	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
546	zdthdp.LT.zdthmin(i) ) THEN
547	zdthmin(i) = zdthdp
548	invb(i) = k
549	ENDIF
550	ENDDO
551	ENDDO
552
553	!
554	! Introduire une diffusion:
555	!
556	IF ( nsrf.EQ.is_oce ) THEN
557	DO k = 2, klev
558	DO i = 1, knon
559	!IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
560	!IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
561	!IM cf JLD/ GKtest TERkz2
562	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
563	! fin GKtest
564
565
566	! s'il n'y a pas d'inversion ou si l'inversion est trop faible
567	! IF ( (nsrf.EQ.is_oce) .AND. &
568	IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN
569	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
570	/(paprs(i,2)-paprs(i,klev+1)) ))**2
571	pcfm(i,k)= zl2(i)* kstable
572	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
573	ENDIF
574	ENDDO
575	ENDDO
576	ENDIF
577
578	END SUBROUTINE coefkz2
579	!
580	!****************************************************************************************
581	!
582	END MODULE coef_diff_turb_mod

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

Download in other formats:

Original Format