Context Navigation

coef_diff_turb_mod.F90 @ 4627

Last change on this file since 4627 was 4230, checked in by fhourdin, 2 years ago
New code for the wakes. Ready for new adventures. Jean-Yves, Lamine, Frederic
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.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(i),i=1,knon)
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	ELSE
185	! No TKE for Standard Physics
186	yq2=0.
187	ENDIF !(iflag_pbl.ge.3)
188
189	END SUBROUTINE coef_diff_turb
190	!
191	!****************************************************************************************
192	!
193	SUBROUTINE coefkz(nsrf, knon, paprs, pplay, &
194	ksta, ksta_ter, &
195	ts, &
196	u,v,t,q, &
197	qsurf, &
198	pcfm, pcfh)
199
200	USE dimphy
201	USE indice_sol_mod
202	USE print_control_mod, ONLY: prt_level, lunout
203
204	!======================================================================
205	! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
206	! (une version strictement identique a l'ancien modele)
207	! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
208	! coefficients d'echange turbulent dans l'atmosphere.
209	! Arguments:
210	! nsrf-----input-I- indicateur de la nature du sol
211	! knon-----input-I- nombre de points a traiter
212	! paprs----input-R- pregssion a chaque intercouche (en Pa)
213	! pplay----input-R- pression au milieu de chaque couche (en Pa)
214	! ts-------input-R- temperature du sol (en Kelvin)
215	! u--------input-R- vitesse u
216	! v--------input-R- vitesse v
217	! t--------input-R- temperature (K)
218	! q--------input-R- vapeur d'eau (kg/kg)
219	!
220	! pcfm-----output-R- coefficients a calculer (vitesse)
221	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
222	!======================================================================
223	INCLUDE "YOETHF.h"
224	INCLUDE "YOMCST.h"
225	INCLUDE "FCTTRE.h"
226	INCLUDE "compbl.h"
227	!
228	! Arguments:
229	!
230	INTEGER, INTENT(IN) :: knon, nsrf
231	REAL, INTENT(IN) :: ksta, ksta_ter
232	REAL, DIMENSION(klon), INTENT(IN) :: ts
233	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
234	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
235	REAL, DIMENSION(klon,klev), INTENT(IN) :: u, v, t, q
236	REAL, DIMENSION(klon), INTENT(IN) :: qsurf
237
238	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pcfm, pcfh
239
240	!
241	! Local variables:
242	!
243	INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite
244	!
245	! Quelques constantes et options:
246	!
247	REAL, PARAMETER :: cepdu2=0.1**2
248	REAL, PARAMETER :: CKAP=0.4
249	REAL, PARAMETER :: cb=5.0
250	REAL, PARAMETER :: cc=5.0
251	REAL, PARAMETER :: cd=5.0
252	REAL, PARAMETER :: clam=160.0
253	REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau
254	LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide
255	REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
256	REAL, PARAMETER :: prandtl=0.4
257	REAL kstable ! diffusion minimale (situation stable)
258	! GKtest
259	! PARAMETER (kstable=1.0e-10)
260	!IM: 261103 REAL kstable_ter, kstable_sinon
261	!IM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
262	!IM: 261103 PARAMETER (kstable_ter = 1.0e-8)
263	!IM: 261103 PARAMETER (kstable_ter = 1.0e-10)
264	!IM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
265	! fin GKtest
266	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
267	INTEGER isommet ! le sommet de la couche limite
268	LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
269	LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere
270
271	!
272	! Variables locales:
273	INTEGER i, k !IM 120704
274	REAL zgeop(klon,klev)
275	REAL zmgeom(klon)
276	REAL zri(klon)
277	REAL zl2(klon)
278	REAL zdphi, zdu2, ztvd, ztvu, zcdn
279	REAL zscf
280	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
281	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
282	REAL, PARAMETER :: t_coup=273.15
283	LOGICAL, PARAMETER :: check=.FALSE.
284	!
285	! contre-gradient pour la chaleur sensible: Kelvin/metre
286	REAL gamt(2:klev)
287
288	LOGICAL, SAVE :: appel1er=.TRUE.
289	!$OMP THREADPRIVATE(appel1er)
290	!
291	! Fonctions thermodynamiques et fonctions d'instabilite
292	REAL fsta, fins, x
293
294	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
295	fins(x) = SQRT(1.0-18.0*x)
296
297	isommet=klev
298
299	IF (appel1er) THEN
300	IF (prt_level > 9) THEN
301	WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
302	WRITE(lunout,*)'coefkz, richum:', richum
303	IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
304	WRITE(lunout,*)'coefkz, isommet:', isommet
305	WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
306	appel1er = .FALSE.
307	ENDIF
308	ENDIF
309	!
310	! Initialiser les sorties
311	!
312	DO k = 1, klev
313	DO i = 1, knon
314	pcfm(i,k) = 0.0
315	pcfh(i,k) = 0.0
316	ENDDO
317	ENDDO
318	DO i = 1, knon
319	itop(i) = 0
320	ENDDO
321
322	!
323	! Prescrire la valeur de contre-gradient
324	!
325	IF (iflag_pbl.EQ.1) THEN
326	DO k = 3, klev
327	gamt(k) = -1.0E-03
328	ENDDO
329	gamt(2) = -2.5E-03
330	ELSE
331	DO k = 2, klev
332	gamt(k) = 0.0
333	ENDDO
334	ENDIF
335	!IM cf JLD/ GKtest
336	IF ( nsrf .NE. is_oce ) THEN
337	!IM 261103 kstable = kstable_ter
338	kstable = ksta_ter
339	ELSE
340	!IM 261103 kstable = kstable_sinon
341	kstable = ksta
342	ENDIF
343	!IM cf JLD/ GKtest fin
344
345	!
346	! Calculer les geopotentiels de chaque couche
347	!
348	DO i = 1, knon
349	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
350	* (paprs(i,1)-pplay(i,1))
351	ENDDO
352	DO k = 2, klev
353	DO i = 1, knon
354	zgeop(i,k) = zgeop(i,k-1) &
355	+ RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) &
356	* (pplay(i,k-1)-pplay(i,k))
357	ENDDO
358	ENDDO
359
360	!
361	! Calculer les coefficients turbulents dans l'atmosphere
362	!
363	DO i = 1, knon
364	itop(i) = isommet
365	ENDDO
366
367
368	DO k = 2, isommet
369	DO i = 1, knon
370	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 &
371	+(v(i,k)-v(i,k-1))**2)
372	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
373	zdphi =zmgeom(i) / 2.0
374	zt = (t(i,k)+t(i,k-1)) * 0.5
375	zq = (q(i,k)+q(i,k-1)) * 0.5
376
377	!
378	! Calculer Qs et dQs/dT:
379	!
380	IF (thermcep) THEN
381	zdelta = MAX(0.,SIGN(1.,RTT-zt))
382	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta) &
383	+ R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
384	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
385	zqs = MIN(0.5,zqs)
386	zcor = 1./(1.-RETV*zqs)
387	zqs = zqs*zcor
388	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
389	ELSE
390	IF (zt .LT. t_coup) THEN
391	zqs = qsats(zt) / pplay(i,k)
392	zdqs = dqsats(zt,zqs)
393	ELSE
394	zqs = qsatl(zt) / pplay(i,k)
395	zdqs = dqsatl(zt,zqs)
396	ENDIF
397	ENDIF
398	!
399	! calculer la fraction nuageuse (processus humide):
400	!
401	if (zq /= 0.) then
402	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
403	else
404	zfr = 0.
405	end if
406	zfr = MAX(0.0,MIN(1.0,zfr))
407	IF (.NOT.richum) zfr = 0.0
408	!
409	! calculer le nombre de Richardson:
410	!
411	IF (tvirtu) THEN
412	ztvd =( t(i,k) &
413	+ zdphi/RCPD/(1.+RVTMP2*zq) &
414	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
415	)(1.+RETVq(i,k))
416	ztvu =( t(i,k-1) &
417	- zdphi/RCPD/(1.+RVTMP2*zq) &
418	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
419	)(1.+RETVq(i,k-1))
420	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
421	zri(i) = zri(i) &
422	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
423	(paprs(i,k)/101325.0)*RKAPPA &
424	/(zdu20.5(ztvd+ztvu))
425
426	ELSE ! calcul de Ridchardson compatible LMD5
427
428	zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
429	-RD0.5(t(i,k)+t(i,k-1))/paprs(i,k) &
430	*(pplay(i,k)-pplay(i,k-1)) &
431	)zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
432	zri(i) = zri(i) + &
433	zmgeom(i)zmgeom(i)gamt(k)/RG &
434	(paprs(i,k)/101325.0)*RKAPPA &
435	/(zdu20.5(t(i,k-1)+t(i,k)))
436	ENDIF
437	!
438	! finalement, les coefficients d'echange sont obtenus:
439	!
440	zcdn=SQRT(zdu2) / zmgeom(i) * RG
441
442	IF (opt_ec) THEN
443	z2geomf=zgeop(i,k-1)+zgeop(i,k)
444	zalm2=(0.5ckap/RGz2geomf &
445	/(1.+0.5ckap/rg/clamz2geomf))**2
446	zalh2=(0.5ckap/rgz2geomf &
447	/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
448	IF (zri(i).LT.0.0) THEN ! situation instable
449	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3 &
450	/ (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
451	zscf = SQRT(-zri(i)*zscf)
452	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
453	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
454	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
455	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
456	ELSE ! situation stable
457	zscf=SQRT(1.+cd*zri(i))
458	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
459	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
460	ENDIF
461	ELSE
462	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
463	/(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
464	pcfm(i,k)=SQRT(MAX(zcdnzcdn(ric-zri(i))/ric, kstable))
465	pcfm(i,k)= zl2(i)* pcfm(i,k)
466	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
467	ENDIF
468	ENDDO
469	ENDDO
470
471	!
472	! Au-dela du sommet, pas de diffusion turbulente:
473	!
474	DO i = 1, knon
475	IF (itop(i)+1 .LE. klev) THEN
476	DO k = itop(i)+1, klev
477	pcfh(i,k) = 0.0
478	pcfm(i,k) = 0.0
479	ENDDO
480	ENDIF
481	ENDDO
482
483	END SUBROUTINE coefkz
484	!
485	!****************************************************************************************
486	!
487	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
488	pcfm, pcfh)
489
490	USE dimphy
491	USE indice_sol_mod
492
493	!======================================================================
494	! J'introduit un peu de diffusion sauf dans les endroits
495	! ou une forte inversion est presente
496	! On peut dire qu'il represente la convection peu profonde
497	!
498	! Arguments:
499	! nsrf-----input-I- indicateur de la nature du sol
500	! knon-----input-I- nombre de points a traiter
501	! paprs----input-R- pression a chaque intercouche (en Pa)
502	! pplay----input-R- pression au milieu de chaque couche (en Pa)
503	! t--------input-R- temperature (K)
504	!
505	! pcfm-----output-R- coefficients a calculer (vitesse)
506	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
507	!======================================================================
508	!
509	! Arguments:
510	!
511	INTEGER, INTENT(IN) :: knon, nsrf
512	REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
513	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
514	REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
515
516	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
517	!
518	! Quelques constantes et options:
519	!
520	REAL, PARAMETER :: prandtl=0.4
521	REAL, PARAMETER :: kstable=0.002
522	! REAL, PARAMETER :: kstable=0.001
523	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
524	REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
525	! PARAMETER (seuil=-0.04)
526	! PARAMETER (seuil=-0.06)
527	! PARAMETER (seuil=-0.09)
528
529	!
530	! Variables locales:
531	!
532	INTEGER i, k, invb(knon)
533	REAL zl2(knon)
534	REAL zdthmin(knon), zdthdp
535
536	INCLUDE "YOMCST.h"
537	!
538	! Initialiser les sorties
539	!
540	DO k = 1, klev
541	DO i = 1, knon
542	pcfm(i,k) = 0.0
543	pcfh(i,k) = 0.0
544	ENDDO
545	ENDDO
546
547	!
548	! Chercher la zone d'inversion forte
549	!
550	DO i = 1, knon
551	invb(i) = klev
552	zdthmin(i)=0.0
553	ENDDO
554	DO k = 2, klev/2-1
555	DO i = 1, knon
556	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
557	- RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
558	zdthdp = zdthdp * 100.0
559	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
560	zdthdp.LT.zdthmin(i) ) THEN
561	zdthmin(i) = zdthdp
562	invb(i) = k
563	ENDIF
564	ENDDO
565	ENDDO
566
567	!
568	! Introduire une diffusion:
569	!
570	IF ( nsrf.EQ.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).EQ.klev) .OR. (zdthmin(i).GT.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