Context Navigation

coef_diff_turb_mod.f90 @ 5282

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

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

Download in other formats:

Original Format