Context Navigation

fisrtilp.F @ 1068

Last change on this file since 1068 was 814, checked in by Laurent Fairhead, 17 years ago
Rajout de la physique utilisant les thermiques FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.1 KB

Line
1	!
2	! $Header$
3	!
4	c
5	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
6	s d_t, d_q, d_ql, rneb, radliq, rain, snow,
7	s pfrac_impa, pfrac_nucl, pfrac_1nucl,
8	s frac_impa, frac_nucl,
9	s prfl, psfl, rhcl)
10
11	c
12	USE dimphy
13	IMPLICIT none
14	c======================================================================
15	c Auteur(s): Z.X. Li (LMD/CNRS)
16	c Date: le 20 mars 1995
17	c Objet: condensation et precipitation stratiforme.
18	c schema de nuage
19	c======================================================================
20	c======================================================================
21	cym#include "dimensions.h"
22	cym#include "dimphy.h"
23	#include "YOMCST.h"
24	#include "tracstoke.h"
25	#include "fisrtilp.h"
26	c
27	c Arguments:
28	c
29	REAL dtime ! intervalle du temps (s)
30	REAL paprs(klon,klev+1) ! pression a inter-couche
31	REAL pplay(klon,klev) ! pression au milieu de couche
32	REAL t(klon,klev) ! temperature (K)
33	REAL q(klon,klev) ! humidite specifique (kg/kg)
34	REAL d_t(klon,klev) ! incrementation de la temperature (K)
35	REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
36	REAL d_ql(klon,klev) ! incrementation de l'eau liquide
37	REAL rneb(klon,klev) ! fraction nuageuse
38	REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
39	REAL rhcl(klon,klev) ! humidite relative en ciel clair
40	REAL rain(klon) ! pluies (mm/s)
41	REAL snow(klon) ! neige (mm/s)
42	REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
43	REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
44	cAA
45	c Coeffients de fraction lessivee : pour OFF-LINE
46	c
47	REAL pfrac_nucl(klon,klev)
48	REAL pfrac_1nucl(klon,klev)
49	REAL pfrac_impa(klon,klev)
50	c
51	c Fraction d'aerosols lessivee par impaction et par nucleation
52	c POur ON-LINE
53	c
54	REAL frac_impa(klon,klev)
55	REAL frac_nucl(klon,klev)
56	real zct(klon),zcl(klon)
57	cAA
58	c
59	c Options du programme:
60	c
61	REAL seuil_neb ! un nuage existe vraiment au-dela
62	PARAMETER (seuil_neb=0.001)
63
64	INTEGER ninter ! sous-intervals pour la precipitation
65	PARAMETER (ninter=5)
66	LOGICAL evap_prec ! evaporation de la pluie
67	PARAMETER (evap_prec=.TRUE.)
68	REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
69	logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
70
71	real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
72	real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
73	real erf
74	c
75	LOGICAL cpartiel ! condensation partielle
76	PARAMETER (cpartiel=.TRUE.)
77	REAL t_coup
78	PARAMETER (t_coup=234.0)
79	c
80	c Variables locales:
81	c
82	INTEGER i, k, n, kk
83	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
84	REAL zrfl(klon), zrfln(klon), zqev, zqevt
85	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
86	REAL ztglace, zt(klon)
87	INTEGER nexpo ! exponentiel pour glace/eau
88	REAL zdz(klon),zrho(klon),ztot(klon), zrhol(klon)
89	REAL zchau(klon),zfroi(klon),zfice(klon),zneb(klon)
90	c
91	LOGICAL appel1er
92	SAVE appel1er
93	c$OMP THREADPRIVATE(appel1er)
94	c
95	c---------------------------------------------------------------
96	c
97	cAA Variables traceurs:
98	cAA Provisoire !!! Parametres alpha du lessivage
99	cAA A priori on a 4 scavenging # possibles
100	c
101	REAL a_tr_sca(4)
102	save a_tr_sca
103	c$OMP THREADPRIVATE(a_tr_sca)
104	c
105	c Variables intermediaires
106	c
107	REAL zalpha_tr
108	REAL zfrac_lessi
109	REAL zprec_cond(klon)
110	cAA
111	REAL zmair, zcpair, zcpeau
112	C Pour la conversion eau-neige
113	REAL zlh_solid(klon), zm_solid
114	cIM
115	cym INTEGER klevm1
116	c---------------------------------------------------------------
117	c
118	c Fonctions en ligne:
119	c
120	REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
121	REAL zzz
122	#include "YOETHF.h"
123	#include "FCTTRE.h"
124	fallvc (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_con
125	fallvs (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_lsc
126	c
127	DATA appel1er /.TRUE./
128	cym
129	zdelq=0.0
130
131	IF (appel1er) THEN
132	c
133	PRINT*, 'fisrtilp, ninter:', ninter
134	PRINT*, 'fisrtilp, evap_prec:', evap_prec
135	PRINT*, 'fisrtilp, cpartiel:', cpartiel
136	IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
137	PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
138	PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
139	c CALL abort
140	ENDIF
141	appel1er = .FALSE.
142	c
143	cAA initialiation provisoire
144	a_tr_sca(1) = -0.5
145	a_tr_sca(2) = -0.5
146	a_tr_sca(3) = -0.5
147	a_tr_sca(4) = -0.5
148	c
149	cAA Initialisation a 1 des coefs des fractions lessivees
150	c
151	DO k = 1, klev
152	DO i = 1, klon
153	pfrac_nucl(i,k)=1.
154	pfrac_1nucl(i,k)=1.
155	pfrac_impa(i,k)=1.
156	ENDDO
157	ENDDO
158
159	ENDIF ! test sur appel1er
160	c
161	cMAf Initialisation a 0 de zoliq
162	DO i = 1, klon
163	zoliq(i)=0.
164	ENDDO
165	c Determiner les nuages froids par leur temperature
166	c nexpo regle la raideur de la transition eau liquide / eau glace.
167	c
168	ztglace = RTT - 15.0
169	nexpo = 6
170	ccc nexpo = 1
171	c
172	c Initialiser les sorties:
173	c
174	DO k = 1, klev+1
175	DO i = 1, klon
176	prfl(i,k) = 0.0
177	psfl(i,k) = 0.0
178	ENDDO
179	ENDDO
180
181	DO k = 1, klev
182	DO i = 1, klon
183	d_t(i,k) = 0.0
184	d_q(i,k) = 0.0
185	d_ql(i,k) = 0.0
186	rneb(i,k) = 0.0
187	radliq(i,k) = 0.0
188	frac_nucl(i,k) = 1.
189	frac_impa(i,k) = 1.
190	ENDDO
191	ENDDO
192	DO i = 1, klon
193	rain(i) = 0.0
194	snow(i) = 0.0
195	ENDDO
196	c
197	c Initialiser le flux de precipitation a zero
198	c
199	DO i = 1, klon
200	zrfl(i) = 0.0
201	zneb(i) = seuil_neb
202	ENDDO
203	c
204	c
205	cAA Pour plus de securite
206
207	zalpha_tr = 0.
208	zfrac_lessi = 0.
209
210	cAA----------------------------------------------------------
211	c
212	c Boucle verticale (du haut vers le bas)
213	c
214	cIM : klevm1
215	cym klevm1=klev-1
216	DO 9999 k = klev, 1, -1
217	c
218	cAA----------------------------------------------------------
219	c
220	DO i = 1, klon
221	zt(i)=t(i,k)
222	zq(i)=q(i,k)
223	ENDDO
224	c
225	c Calculer la varition de temp. de l'air du a la chaleur sensible
226	C transporter par la pluie.
227	C Il resterait a rajouter cet effet de la chaleur sensible sur les
228	C flux de surface, du a la diff. de temp. entre le 1er niveau et la
229	C surface.
230	C
231	DO i = 1, klon
232	cIM
233	IF(k.LE.klevm1) THEN
234	zmair=(paprs(i,k)-paprs(i,k+1))/RG
235	zcpair=RCPD(1.0+RVTMP2zq(i))
236	zcpeau=RCPD*RVTMP2
237	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zrfl(i)dtime*zcpeau
238	$ + zmairzcpairzt(i) )
239	$ / (zmairzcpair + zrfl(i)dtime*zcpeau)
240	CC WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
241	ENDIF
242	ENDDO
243	c
244	c
245	c Calculer l'evaporation de la precipitation
246	c
247
248
249	IF (evap_prec) THEN
250	DO i = 1, klon
251	IF (zrfl(i) .GT.0.) THEN
252	IF (thermcep) THEN
253	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
254	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
255	zqs(i)=MIN(0.5,zqs(i))
256	zcor=1./(1.-RETV*zqs(i))
257	zqs(i)=zqs(i)*zcor
258	ELSE
259	IF (zt(i) .LT. t_coup) THEN
260	zqs(i) = qsats(zt(i)) / pplay(i,k)
261	ELSE
262	zqs(i) = qsatl(zt(i)) / pplay(i,k)
263	ENDIF
264	ENDIF
265	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
266	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
267	. * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
268	zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
269	. * RG*dtime/(paprs(i,k)-paprs(i,k+1))
270	zqev = MIN (zqev, zqevt)
271	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
272	. /RG/dtime
273
274	c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous
275	c la glace venant de la couche du dessus est simplement dans la couche
276	c du dessous.
277
278	IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
279
280	zq(i) = zq(i) - (zrfln(i)-zrfl(i))
281	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
282	zt(i) = zt(i) + (zrfln(i)-zrfl(i))
283	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
284	. * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
285	zrfl(i) = zrfln(i)
286	ENDIF
287	ENDDO
288	ENDIF
289	c
290	c Calculer Qs et L/Cp*dQs/dT:
291	c
292	IF (thermcep) THEN
293	DO i = 1, klon
294	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
295	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
296	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
297	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
298	zqs(i) = MIN(0.5,zqs(i))
299	zcor = 1./(1.-RETV*zqs(i))
300	zqs(i) = zqs(i)*zcor
301	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
302	ENDDO
303	ELSE
304	DO i = 1, klon
305	IF (zt(i).LT.t_coup) THEN
306	zqs(i) = qsats(zt(i))/pplay(i,k)
307	zdqs(i) = dqsats(zt(i),zqs(i))
308	ELSE
309	zqs(i) = qsatl(zt(i))/pplay(i,k)
310	zdqs(i) = dqsatl(zt(i),zqs(i))
311	ENDIF
312	ENDDO
313	ENDIF
314	c
315	c Determiner la condensation partielle et calculer la quantite
316	c de l'eau condensee:
317	c
318	IF (cpartiel) THEN
319
320	c print*,'Dans partiel k=',k
321	c
322	c Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
323	c nuageuse a partir des PDF de Sandrine Bony.
324	c rneb : fraction nuageuse
325	c zqn : eau totale dans le nuage
326	c zcond : eau condensee moyenne dans la maille.
327	c on prend en compte le r�hauffement qui diminue la partie condensee
328	c
329	c Version avec les raqts
330
331	if (iflag_pdf.eq.0) then
332
333	do i=1,klon
334	zdelq = min(ratqs(i,k),0.99) * zq(i)
335	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
336	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
337	enddo
338
339	else
340	c
341	c Version avec les nouvelles PDFs.
342	do i=1,klon
343	if(zq(i).lt.1.e-15) then
344	print*,'ZQ(',i,',',k,')=',zq(i)
345	zq(i)=1.e-15
346	endif
347	enddo
348	do i=1,klon
349	zpdf_sig(i)=ratqs(i,k)*zq(i)
350	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
351	zpdf_delta(i)=log(zq(i)/zqs(i))
352	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
353	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
354	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
355	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
356	zpdf_e1(i)=1.-erf(zpdf_e1(i))
357	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
358	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
359	zpdf_e2(i)=1.-erf(zpdf_e2(i))
360	if (zpdf_e1(i).lt.1.e-10) then
361	rneb(i,k)=0.
362	zqn(i)=zqs(i)
363	else
364	rneb(i,k)=0.5*zpdf_e1(i)
365	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
366	endif
367
368	enddo
369
370	endif ! iflag_pdf
371
372	do i=1,klon
373	IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
374	IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
375	rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
376	c zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
377	c On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
378	c la convection.
379	c ATTENTION !!! Il va falloir verifier tout ca.
380	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
381	c print*,'ZDQS ',zdqs(i)
382	c--Olivier
383	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
384	IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
385	IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
386	c--fin
387	ENDDO
388	ELSE
389	DO i = 1, klon
390	IF (zq(i).GT.zqs(i)) THEN
391	rneb(i,k) = 1.0
392	ELSE
393	rneb(i,k) = 0.0
394	ENDIF
395	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
396	ENDDO
397	ENDIF
398	c
399	DO i = 1, klon
400	zq(i) = zq(i) - zcond(i)
401	c zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
402	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
403	ENDDO
404	c
405	c Partager l'eau condensee en precipitation et eau liquide nuageuse
406	c
407	DO i = 1, klon
408	IF (rneb(i,k).GT.0.0) THEN
409	zoliq(i) = zcond(i)
410	zrho(i) = pplay(i,k) / zt(i) / RD
411	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
412	zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
413	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
414	zfice(i) = zfice(i)**nexpo
415	zneb(i) = MAX(rneb(i,k), seuil_neb)
416	radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
417	ENDIF
418	ENDDO
419	c
420	DO n = 1, ninter
421	DO i = 1, klon
422	IF (rneb(i,k).GT.0.0) THEN
423	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
424
425	if (ptconv(i,k)) then
426	zcl(i)=cld_lc_con
427	zct(i)=1./cld_tau_con
428	else
429	zcl(i)=cld_lc_lsc
430	zct(i)=1./cld_tau_lsc
431	endif
432	c quantit�d'eau ��minier.
433	zchau(i) = zct(i)dtime/FLOAT(ninter) zoliq(i)
434	. (1.0-EXP(-(zoliq(i)/zneb(i)/zcl(i))2)) (1.-zfice(i))
435	c meme chose pour la glace.
436	if (ptconv(i,k)) then
437	zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
438	. fallvc(zrhol(i)) zfice(i)
439	else
440	zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
441	. fallvs(zrhol(i)) zfice(i)
442	endif
443	ztot(i) = zchau(i) + zfroi(i)
444	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
445	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
446	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
447	radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
448	ENDIF
449	ENDDO
450	ENDDO
451	c
452	DO i = 1, klon
453	IF (rneb(i,k).GT.0.0) THEN
454	d_ql(i,k) = zoliq(i)
455	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
456	. * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
457	ENDIF
458	IF (zt(i).LT.RTT) THEN
459	psfl(i,k)=zrfl(i)
460	ELSE
461	prfl(i,k)=zrfl(i)
462	ENDIF
463	ENDDO
464	c
465	c Calculer les tendances de q et de t:
466	c
467	DO i = 1, klon
468	d_q(i,k) = zq(i) - q(i,k)
469	d_t(i,k) = zt(i) - t(i,k)
470	ENDDO
471	c
472	cAA--------------- Calcul du lessivage stratiforme -------------
473
474	DO i = 1,klon
475	c
476	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
477	. * (paprs(i,k)-paprs(i,k+1))/RG
478	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
479	cAA lessivage nucleation LMD5 dans la couche elle-meme
480	if (t(i,k) .GE. ztglace) THEN
481	zalpha_tr = a_tr_sca(3)
482	else
483	zalpha_tr = a_tr_sca(4)
484	endif
485	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
486	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
487	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
488	c
489	c nucleation avec un facteur -1 au lieu de -0.5
490	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
491	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
492	ENDIF
493	c
494	ENDDO ! boucle sur i
495	c
496	cAA Lessivage par impaction dans les couches en-dessous
497	DO kk = k-1, 1, -1
498	DO i = 1, klon
499	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
500	if (t(i,kk) .GE. ztglace) THEN
501	zalpha_tr = a_tr_sca(1)
502	else
503	zalpha_tr = a_tr_sca(2)
504	endif
505	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
506	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
507	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
508	ENDIF
509	ENDDO
510	ENDDO
511	c
512	cAA----------------------------------------------------------
513	c FIN DE BOUCLE SUR K
514	9999 CONTINUE
515	c
516	cAA-----------------------------------------------------------
517	c
518	c Pluie ou neige au sol selon la temperature de la 1ere couche
519	c
520	DO i = 1, klon
521	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
522	snow(i) = zrfl(i)
523	zlh_solid(i) = RLSTT-RLVTT
524	ELSE
525	rain(i) = zrfl(i)
526	zlh_solid(i) = 0.
527	ENDIF
528	ENDDO
529	C
530	C For energy conservation : when snow is present, the solification
531	c latent heat is considered.
532	DO k = 1, klev
533	DO i = 1, klon
534	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
535	zmair=(paprs(i,k)-paprs(i,k+1))/RG
536	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
537	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair)
538	END DO
539	END DO
540	c
541	RETURN
542	END

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

Download in other formats:

Original Format