Context Navigation

coef_diff_turb_mod.F90 @ 1023

Last change on this file since 1023 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 20.6 KB

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

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

Download in other formats:

Original Format