Context Navigation

coef_diff_turb_mod.F90 @ 1691

Last change on this file since 1691 was 1604, checked in by lguez, 13 years ago

Removed two "fi" with no corresponding "if" in "makdim".

In procedure "coefkz", when calculating cloud fraction "zfr", "zq" can
be zero at high altitude. Avoid division by zero, set "zfr" to 0 when
"zq" is 0.

In "physiq", allocatable arrays "tabijgcm", "longcm", "latgcm",
"igcm", "jgcm" cannot be arguments of "phys_output_open" if they have
not been allocated. Allocate them with zero size when
'npCFMIP_param.data' cannot be opened.

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	!
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	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
482	!======================================================================
483	! J'introduit un peu de diffusion sauf dans les endroits
484	! ou une forte inversion est presente
485	! On peut dire qu'il represente la convection peu profonde
486	!
487	! Arguments:
488	! nsrf-----input-I- indicateur de la nature du sol
489	! knon-----input-I- nombre de points a traiter
490	! paprs----input-R- pression a chaque intercouche (en Pa)
491	! pplay----input-R- pression au milieu de chaque couche (en Pa)
492	! t--------input-R- temperature (K)
493	!
494	! pcfm-----output-R- coefficients a calculer (vitesse)
495	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
496	!======================================================================
497	!
498	! Arguments:
499	!
500	INTEGER, INTENT(IN) :: knon, nsrf
501	REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
502	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
503	REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
504
505	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
506	!
507	! Quelques constantes et options:
508	!
509	REAL, PARAMETER :: prandtl=0.4
510	REAL, PARAMETER :: kstable=0.002
511	! REAL, PARAMETER :: kstable=0.001
512	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
513	REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
514	! PARAMETER (seuil=-0.04)
515	! PARAMETER (seuil=-0.06)
516	! PARAMETER (seuil=-0.09)
517
518	!
519	! Variables locales:
520	!
521	INTEGER i, k, invb(knon)
522	REAL zl2(knon)
523	REAL zdthmin(knon), zdthdp
524
525	INCLUDE "indicesol.h"
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