Context Navigation

fisrtilp.F @ 1203

Last change on this file since 1203 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: 16.1 KB

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

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

Download in other formats:

Original Format