Context Navigation

fisrtilp_tr.F @ 1582

Last change on this file since 1582 was 1575, checked in by jghattas, 13 years ago

Added suffix _mpi_rank to lmdz.out text file. Each processus now write into seperate file but only if lunout/=6 (lunout is set in run.def).

Change some print* into write(lunout,*)

phytrac.F90 : always include ini_histrac and write_histrac. The file histrac.nc is written if ecrit_tra> 0 (set in physiq.def). Change default value of ecrit_tra into 0.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 12.9 KB

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

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

Download in other formats:

Original Format