Context Navigation

fisrtilp_tr.F @ 1098

Last change on this file since 1098 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: 12.8 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$OMP THREADPRIVATE(appel1er)
105	c
106	c---------------------------------------------------------------
107	c
108	cAA Variables traceurs:
109	cAA Provisoire !!! Parametres alpha du lessivage
110	cAA A priori on a 4 scavenging # possibles
111	c
112	REAL a_tr_sca(4)
113	save a_tr_sca
114	c$OMP THREADPRIVATE(a_tr_sca)
115	c
116	c Variables intermediaires
117	c
118	REAL zalpha_tr
119	REAL zfrac_lessi
120	REAL zprec_cond(klon)
121	cAA
122	c---------------------------------------------------------------
123	c
124	c Fonctions en ligne:
125	c
126	REAL fallv ! vitesse de chute pour crystaux de glace
127	REAL zzz
128	#include "YOETHF.h"
129	#include "FCTTRE.h"
130	fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
131	ccc fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
132	ccc fallv (zzz) = 3.29 * ((zzz)**0.16)
133	c
134	DATA appel1er /.TRUE./
135	c
136	IF (appel1er) THEN
137	c
138	PRINT*, 'fisrtilp, calcrat:', calcrat
139	PRINT*, 'fisrtilp, ninter:', ninter
140	PRINT*, 'fisrtilp, evap_prec:', evap_prec
141	PRINT*, 'fisrtilp, cpartiel:', cpartiel
142	IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
143	PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
144	PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
145	CALL abort
146	ENDIF
147	appel1er = .FALSE.
148	c
149	cAA initialiation provisoire
150	a_tr_sca(1) = -0.5
151	a_tr_sca(2) = -0.5
152	a_tr_sca(3) = -0.5
153	a_tr_sca(4) = -0.5
154	c
155	cAA Initialisation a 1 des coefs des fractions lessivees
156	c
157	DO k = 1, klev
158	DO i = 1, klon
159	pfrac_nucl(i,k)=1.
160	pfrac_1nucl(i,k)=1.
161	pfrac_impa(i,k)=1.
162	ENDDO
163	ENDDO
164
165	ENDIF ! test sur appel1er
166	c
167	cMAf Initialisation a 0 de zoliq
168	DO i = 1, klon
169	zoliq(i)=0.
170	ENDDO
171	c Determiner les nuages froids par leur temperature
172	c
173	ztglace = RTT - 15.0
174	nexpo = 6
175	ccc nexpo = 1
176	c
177	c Initialiser les sorties:
178	c
179	DO k = 1, klev+1
180	DO i = 1, klon
181	prfl(i,k) = 0.0
182	psfl(i,k) = 0.0
183	ENDDO
184	ENDDO
185
186	DO k = 1, klev
187	DO i = 1, klon
188	d_t(i,k) = 0.0
189	d_q(i,k) = 0.0
190	d_ql(i,k) = 0.0
191	rneb(i,k) = 0.0
192	radliq(i,k) = 0.0
193	frac_nucl(i,k) = 1.
194	frac_impa(i,k) = 1.
195	ENDDO
196	ENDDO
197	DO i = 1, klon
198	rain(i) = 0.0
199	snow(i) = 0.0
200	ENDDO
201	c
202	c Initialiser le flux de precipitation a zero
203	c
204	DO i = 1, klon
205	zrfl(i) = 0.0
206	zneb(i) = seuil_neb
207	ENDDO
208	c
209	c
210	cAA Pour plus de securite
211
212	zalpha_tr = 0.
213	zfrac_lessi = 0.
214
215	cAA----------------------------------------------------------
216	c
217	c Boucle verticale (du haut vers le bas)
218	c
219	DO 9999 k = klev, 1, -1
220	c
221	cAA----------------------------------------------------------
222	c
223	DO i = 1, klon
224	zt(i)=t(i,k)
225	zq(i)=q(i,k)
226	ENDDO
227	c
228	c Calculer l'evaporation de la precipitation
229	c
230	IF (evap_prec) THEN
231	DO i = 1, klon
232	IF (zrfl(i) .GT.0.) THEN
233	IF (thermcep) THEN
234	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
235	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
236	zqs(i)=MIN(0.5,zqs(i))
237	zcor=1./(1.-RETV*zqs(i))
238	zqs(i)=zqs(i)*zcor
239	ELSE
240	IF (zt(i) .LT. t_coup) THEN
241	zqs(i) = qsats(zt(i)) / pplay(i,k)
242	ELSE
243	zqs(i) = qsatl(zt(i)) / pplay(i,k)
244	ENDIF
245	ENDIF
246	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
247	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
248	. * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
249	zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
250	. * RG*dtime/(paprs(i,k)-paprs(i,k+1))
251	zqev = MIN (zqev, zqevt)
252	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
253	. /RG/dtime
254	zq(i) = zq(i) - (zrfln(i)-zrfl(i))
255	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
256	zt(i) = zt(i) + (zrfln(i)-zrfl(i))
257	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
258	. * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
259	zrfl(i) = zrfln(i)
260	ENDIF
261	ENDDO
262	ENDIF
263	c
264	c Calculer Qs et L/Cp*dQs/dT:
265	c
266	IF (thermcep) THEN
267	DO i = 1, klon
268	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
269	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
270	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
271	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
272	zqs(i) = MIN(0.5,zqs(i))
273	zcor = 1./(1.-RETV*zqs(i))
274	zqs(i) = zqs(i)*zcor
275	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
276	ENDDO
277	ELSE
278	DO i = 1, klon
279	IF (zt(i).LT.t_coup) THEN
280	zqs(i) = qsats(zt(i))/pplay(i,k)
281	zdqs(i) = dqsats(zt(i),zqs(i))
282	ELSE
283	zqs(i) = qsatl(zt(i))/pplay(i,k)
284	zdqs(i) = dqsatl(zt(i),zqs(i))
285	ENDIF
286	ENDDO
287	ENDIF
288	c
289	c Determiner la condensation partielle et calculer la quantite
290	c de l'eau condensee:
291	c
292	IF (cpartiel) THEN
293	DO i = 1, klon
294	c
295	zdelq = ratqs(i,k) * zq(i)
296	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
297	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
298	IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
299	IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
300	rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
301	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
302
303	c--Olivier
304	RHcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
305	IF (rneb(i,k) .LE. 0.0) RHcl(i,k)=zq(i)/zqs(i)
306	IF (rneb(i,k) .GE. 1.0) RHcl(i,k)=1.0
307	c--fin
308
309	ENDDO
310	ELSE
311	DO i = 1, klon
312	IF (zq(i).GT.zqs(i)) THEN
313	rneb(i,k) = 1.0
314	ELSE
315	rneb(i,k) = 0.0
316	ENDIF
317	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
318	ENDDO
319	ENDIF
320	c
321	DO i = 1, klon
322	zq(i) = zq(i) - zcond(i)
323	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
324	ENDDO
325	c
326	c Partager l'eau condensee en precipitation et eau liquide nuageuse
327	c
328	DO i = 1, klon
329	IF (rneb(i,k).GT.0.0) THEN
330	zoliq(i) = zcond(i)
331	zrho(i) = pplay(i,k) / zt(i) / RD
332	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
333	zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
334	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
335	zfice(i) = zfice(i)**nexpo
336	zneb(i) = MAX(rneb(i,k), seuil_neb)
337	radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
338	ENDIF
339	ENDDO
340	c
341	DO n = 1, ninter
342	DO i = 1, klon
343	IF (rneb(i,k).GT.0.0) THEN
344	zchau(i) = ctdtime/FLOAT(ninter) zoliq(i)
345	. * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) (1.-zfice(i))
346	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
347	zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
348	. fallv(zrhol(i)) zfice(i)
349	ztot(i) = zchau(i) + zfroi(i)
350	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
351	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
352	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
353	radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
354	ENDIF
355	ENDDO
356	ENDDO
357	c
358	DO i = 1, klon
359	IF (rneb(i,k).GT.0.0) THEN
360	d_ql(i,k) = zoliq(i)
361	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
362	. * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
363	ENDIF
364	IF (zt(i).LT.RTT) THEN
365	psfl(i,k)=zrfl(i)
366	ELSE
367	prfl(i,k)=zrfl(i)
368	ENDIF
369	ENDDO
370	c
371	c Calculer les tendances de q et de t:
372	c
373	DO i = 1, klon
374	d_q(i,k) = zq(i) - q(i,k)
375	d_t(i,k) = zt(i) - t(i,k)
376	ENDDO
377	c
378	cAA--------------- Calcul du lessivage stratiforme -------------
379
380	DO i = 1,klon
381	c
382	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
383	. * (paprs(i,k)-paprs(i,k+1))/RG
384	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
385	cAA lessivage nucleation LMD5 dans la couche elle-meme
386	if (t(i,k) .GE. ztglace) THEN
387	zalpha_tr = a_tr_sca(3)
388	else
389	zalpha_tr = a_tr_sca(4)
390	endif
391	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
392	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
393	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
394	c
395	c nucleation avec un facteur -1 au lieu de -0.5
396	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
397	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
398	ENDIF
399	c
400	ENDDO ! boucle sur i
401	c
402	cAA Lessivage par impaction dans les couches en-dessous
403	DO kk = k-1, 1, -1
404	DO i = 1, klon
405	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
406	if (t(i,kk) .GE. ztglace) THEN
407	zalpha_tr = a_tr_sca(1)
408	else
409	zalpha_tr = a_tr_sca(2)
410	endif
411	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
412	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
413	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
414	ENDIF
415	ENDDO
416	ENDDO
417	c
418	cAA----------------------------------------------------------
419	c FIN DE BOUCLE SUR K
420	9999 CONTINUE
421	c
422	cAA-----------------------------------------------------------
423	c
424	c Pluie ou neige au sol selon la temperature de la 1ere couche
425	c
426	DO i = 1, klon
427	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
428	snow(i) = zrfl(i)
429	ELSE
430	rain(i) = zrfl(i)
431	ENDIF
432	ENDDO
433	c
434	RETURN
435	END

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

Download in other formats:

Original Format