Context Navigation

coef_diff_turb_mod.F90 @ 4654

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

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

Download in other formats:

Original Format