Context Navigation

coef_diff_turb_mod.f90 @ 5426

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

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

Download in other formats:

Original Format