Context Navigation

fisrtilp_tr.F @ 2101

Last change on this file since 2101 was 524, checked in by lmdzadmin, 21 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.7 KB

Line
1	!
2	! $Header$
3	!
4	c
5	SUBROUTINE fisrtilp_tr(dtime,paprs,pplay,t,q,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,
10	s RHcl) ! relative humidity in clear sky (needed for aer optical properties; aeropt.F)
11
12	c
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	#include "dimensions.h"
22	#include "dimphy.h"
23	#include "YOMCST.h"
24	#include "tracstoke.h"
25	c
26	c Arguments:
27	c
28	REAL dtime ! intervalle du temps (s)
29	REAL paprs(klon,klev+1) ! pression a inter-couche
30	REAL pplay(klon,klev) ! pression au milieu de couche
31	REAL t(klon,klev) ! temperature (K)
32	REAL q(klon,klev) ! humidite specifique (kg/kg)
33	REAL d_t(klon,klev) ! incrementation de la temperature (K)
34	REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
35	REAL d_ql(klon,klev) ! incrementation de l'eau liquide
36	REAL rneb(klon,klev) ! fraction nuageuse
37	REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
38	REAL rain(klon) ! pluies (mm/s)
39	REAL snow(klon) ! neige (mm/s)
40	REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
41	REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
42
43	Cjq For aerosol opt properties needed (see aeropt.F)
44	REAL RHcl(klon,klev)
45
46	cAA
47	c Coeffients de fraction lessivee : pour OFF-LINE
48	c
49	REAL pfrac_nucl(klon,klev)
50	REAL pfrac_1nucl(klon,klev)
51	REAL pfrac_impa(klon,klev)
52	c
53	c Fraction d'aerosols lessivee par impaction et par nucleation
54	c POur ON-LINE
55	c
56	REAL frac_impa(klon,klev)
57	REAL frac_nucl(klon,klev)
58	cAA
59	c
60	c Options du programme:
61	c
62	REAL seuil_neb ! un nuage existe vraiment au-dela
63	PARAMETER (seuil_neb=0.001)
64	REAL ct ! inverse du temps pour qu'un nuage precipite
65	PARAMETER (ct=1./1800.)
66	REAL cl ! seuil de precipitation
67	PARAMETER (cl=2.6e-4)
68	ccc PARAMETER (cl=2.3e-4)
69	ccc PARAMETER (cl=2.0e-4)
70	INTEGER ninter ! sous-intervals pour la precipitation
71	PARAMETER (ninter=5)
72	LOGICAL evap_prec ! evaporation de la pluie
73	PARAMETER (evap_prec=.TRUE.)
74	REAL coef_eva
75	PARAMETER (coef_eva=2.0E-05)
76	LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
77	REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
78	PARAMETER (calcrat=.TRUE.)
79	REAL zx_min, rat_max
80	PARAMETER (zx_min=1.0, rat_max=0.01)
81	REAL zx_max, rat_min
82	PARAMETER (zx_max=0.1, rat_min=0.3)
83	REAL zx
84	c
85	LOGICAL cpartiel ! condensation partielle
86	PARAMETER (cpartiel=.TRUE.)
87	REAL t_coup
88	PARAMETER (t_coup=234.0)
89	c
90	c Variables locales:
91	c
92	INTEGER i, k, n, kk
93	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
94	REAL zrfl(klon), zrfln(klon), zqev, zqevt
95	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
96	REAL ztglace, zt(klon)
97	INTEGER nexpo ! exponentiel pour glace/eau
98	REAL zdz(klon),zrho(klon),ztot(klon), zrhol(klon)
99	REAL zchau(klon),zfroi(klon),zfice(klon),zneb(klon)
100	c
101	LOGICAL appel1er
102	SAVE appel1er
103	c
104	c---------------------------------------------------------------
105	c
106	cAA Variables traceurs:
107	cAA Provisoire !!! Parametres alpha du lessivage
108	cAA A priori on a 4 scavenging # possibles
109	c
110	REAL a_tr_sca(4)
111	save a_tr_sca
112	c
113	c Variables intermediaires
114	c
115	REAL zalpha_tr
116	REAL zfrac_lessi
117	REAL zprec_cond(klon)
118	cAA
119	c---------------------------------------------------------------
120	c
121	c Fonctions en ligne:
122	c
123	REAL fallv ! vitesse de chute pour crystaux de glace
124	REAL zzz
125	#include "YOETHF.h"
126	#include "FCTTRE.h"
127	fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
128	ccc fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
129	ccc fallv (zzz) = 3.29 * ((zzz)**0.16)
130	c
131	DATA appel1er /.TRUE./
132	c
133	IF (appel1er) THEN
134	c
135	PRINT*, 'fisrtilp, calcrat:', calcrat
136	PRINT*, 'fisrtilp, ninter:', ninter
137	PRINT*, 'fisrtilp, evap_prec:', evap_prec
138	PRINT*, 'fisrtilp, cpartiel:', cpartiel
139	IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
140	PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
141	PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
142	CALL abort
143	ENDIF
144	appel1er = .FALSE.
145	c
146	cAA initialiation provisoire
147	a_tr_sca(1) = -0.5
148	a_tr_sca(2) = -0.5
149	a_tr_sca(3) = -0.5
150	a_tr_sca(4) = -0.5
151	c
152	cAA Initialisation a 1 des coefs des fractions lessivees
153	c
154	DO k = 1, klev
155	DO i = 1, klon
156	pfrac_nucl(i,k)=1.
157	pfrac_1nucl(i,k)=1.
158	pfrac_impa(i,k)=1.
159	ENDDO
160	ENDDO
161
162	ENDIF ! test sur appel1er
163	c
164	cMAf Initialisation a 0 de zoliq
165	DO i = 1, klon
166	zoliq(i)=0.
167	ENDDO
168	c Determiner les nuages froids par leur temperature
169	c
170	ztglace = RTT - 15.0
171	nexpo = 6
172	ccc nexpo = 1
173	c
174	c Initialiser les sorties:
175	c
176	DO k = 1, klev+1
177	DO i = 1, klon
178	prfl(i,k) = 0.0
179	psfl(i,k) = 0.0
180	ENDDO
181	ENDDO
182
183	DO k = 1, klev
184	DO i = 1, klon
185	d_t(i,k) = 0.0
186	d_q(i,k) = 0.0
187	d_ql(i,k) = 0.0
188	rneb(i,k) = 0.0
189	radliq(i,k) = 0.0
190	frac_nucl(i,k) = 1.
191	frac_impa(i,k) = 1.
192	ENDDO
193	ENDDO
194	DO i = 1, klon
195	rain(i) = 0.0
196	snow(i) = 0.0
197	ENDDO
198	c
199	c Initialiser le flux de precipitation a zero
200	c
201	DO i = 1, klon
202	zrfl(i) = 0.0
203	zneb(i) = seuil_neb
204	ENDDO
205	c
206	c
207	cAA Pour plus de securite
208
209	zalpha_tr = 0.
210	zfrac_lessi = 0.
211
212	cAA----------------------------------------------------------
213	c
214	c Boucle verticale (du haut vers le bas)
215	c
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 l'evaporation de la precipitation
226	c
227	IF (evap_prec) THEN
228	DO i = 1, klon
229	IF (zrfl(i) .GT.0.) THEN
230	IF (thermcep) THEN
231	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
232	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
233	zqs(i)=MIN(0.5,zqs(i))
234	zcor=1./(1.-RETV*zqs(i))
235	zqs(i)=zqs(i)*zcor
236	ELSE
237	IF (zt(i) .LT. t_coup) THEN
238	zqs(i) = qsats(zt(i)) / pplay(i,k)
239	ELSE
240	zqs(i) = qsatl(zt(i)) / pplay(i,k)
241	ENDIF
242	ENDIF
243	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
244	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
245	. * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
246	zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
247	. * RG*dtime/(paprs(i,k)-paprs(i,k+1))
248	zqev = MIN (zqev, zqevt)
249	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
250	. /RG/dtime
251	zq(i) = zq(i) - (zrfln(i)-zrfl(i))
252	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
253	zt(i) = zt(i) + (zrfln(i)-zrfl(i))
254	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
255	. * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
256	zrfl(i) = zrfln(i)
257	ENDIF
258	ENDDO
259	ENDIF
260	c
261	c Calculer Qs et L/Cp*dQs/dT:
262	c
263	IF (thermcep) THEN
264	DO i = 1, klon
265	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
266	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
267	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
268	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
269	zqs(i) = MIN(0.5,zqs(i))
270	zcor = 1./(1.-RETV*zqs(i))
271	zqs(i) = zqs(i)*zcor
272	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
273	ENDDO
274	ELSE
275	DO i = 1, klon
276	IF (zt(i).LT.t_coup) THEN
277	zqs(i) = qsats(zt(i))/pplay(i,k)
278	zdqs(i) = dqsats(zt(i),zqs(i))
279	ELSE
280	zqs(i) = qsatl(zt(i))/pplay(i,k)
281	zdqs(i) = dqsatl(zt(i),zqs(i))
282	ENDIF
283	ENDDO
284	ENDIF
285	c
286	c Determiner la condensation partielle et calculer la quantite
287	c de l'eau condensee:
288	c
289	IF (cpartiel) THEN
290	DO i = 1, klon
291	c
292	zdelq = ratqs(i,k) * zq(i)
293	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
294	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
295	IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
296	IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
297	rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
298	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
299
300	c--Olivier
301	RHcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
302	IF (rneb(i,k) .LE. 0.0) RHcl(i,k)=zq(i)/zqs(i)
303	IF (rneb(i,k) .GE. 1.0) RHcl(i,k)=1.0
304	c--fin
305
306	ENDDO
307	ELSE
308	DO i = 1, klon
309	IF (zq(i).GT.zqs(i)) THEN
310	rneb(i,k) = 1.0
311	ELSE
312	rneb(i,k) = 0.0
313	ENDIF
314	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
315	ENDDO
316	ENDIF
317	c
318	DO i = 1, klon
319	zq(i) = zq(i) - zcond(i)
320	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
321	ENDDO
322	c
323	c Partager l'eau condensee en precipitation et eau liquide nuageuse
324	c
325	DO i = 1, klon
326	IF (rneb(i,k).GT.0.0) THEN
327	zoliq(i) = zcond(i)
328	zrho(i) = pplay(i,k) / zt(i) / RD
329	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
330	zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
331	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
332	zfice(i) = zfice(i)**nexpo
333	zneb(i) = MAX(rneb(i,k), seuil_neb)
334	radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
335	ENDIF
336	ENDDO
337	c
338	DO n = 1, ninter
339	DO i = 1, klon
340	IF (rneb(i,k).GT.0.0) THEN
341	zchau(i) = ctdtime/FLOAT(ninter) zoliq(i)
342	. * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) (1.-zfice(i))
343	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
344	zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
345	. fallv(zrhol(i)) zfice(i)
346	ztot(i) = zchau(i) + zfroi(i)
347	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
348	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
349	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
350	radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
351	ENDIF
352	ENDDO
353	ENDDO
354	c
355	DO i = 1, klon
356	IF (rneb(i,k).GT.0.0) THEN
357	d_ql(i,k) = zoliq(i)
358	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
359	. * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
360	ENDIF
361	IF (zt(i).LT.RTT) THEN
362	psfl(i,k)=zrfl(i)
363	ELSE
364	prfl(i,k)=zrfl(i)
365	ENDIF
366	ENDDO
367	c
368	c Calculer les tendances de q et de t:
369	c
370	DO i = 1, klon
371	d_q(i,k) = zq(i) - q(i,k)
372	d_t(i,k) = zt(i) - t(i,k)
373	ENDDO
374	c
375	cAA--------------- Calcul du lessivage stratiforme -------------
376
377	DO i = 1,klon
378	c
379	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
380	. * (paprs(i,k)-paprs(i,k+1))/RG
381	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
382	cAA lessivage nucleation LMD5 dans la couche elle-meme
383	if (t(i,k) .GE. ztglace) THEN
384	zalpha_tr = a_tr_sca(3)
385	else
386	zalpha_tr = a_tr_sca(4)
387	endif
388	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
389	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
390	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
391	c
392	c nucleation avec un facteur -1 au lieu de -0.5
393	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
394	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
395	ENDIF
396	c
397	ENDDO ! boucle sur i
398	c
399	cAA Lessivage par impaction dans les couches en-dessous
400	DO kk = k-1, 1, -1
401	DO i = 1, klon
402	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
403	if (t(i,kk) .GE. ztglace) THEN
404	zalpha_tr = a_tr_sca(1)
405	else
406	zalpha_tr = a_tr_sca(2)
407	endif
408	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
409	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
410	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
411	ENDIF
412	ENDDO
413	ENDDO
414	c
415	cAA----------------------------------------------------------
416	c FIN DE BOUCLE SUR K
417	9999 CONTINUE
418	c
419	cAA-----------------------------------------------------------
420	c
421	c Pluie ou neige au sol selon la temperature de la 1ere couche
422	c
423	DO i = 1, klon
424	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
425	snow(i) = zrfl(i)
426	ELSE
427	rain(i) = zrfl(i)
428	ENDIF
429	ENDDO
430	c
431	RETURN
432	END

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

Download in other formats:

Original Format