Context Navigation

coef_diff_turb_mod.F90 @ 2518

Last change on this file since 2518 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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

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

Download in other formats:

Original Format