Context Navigation

fisrtilp_tr.F @ 3919

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

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

Download in other formats:

Original Format