Context Navigation

coef_diff_turb_mod.F90 @ 3603

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

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

Download in other formats:

Original Format