Context Navigation

fisrtilp.F @ 380

Last change on this file since 380 was 380, checked in by lmdzadmin, 22 years ago
PB dans l'include de fisrtilp.h/inc LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 14.7 KB

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

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

Download in other formats:

Original Format