Context Navigation

fisrtilp.F90 @ 2600

Last change on this file since 2600 was 2507, checked in by jbmadeleine, 9 years ago
Added the iflag_t_glace = 2 option that turns off the convertion of cloud water to ice close to the surface when T< t_glace_max (273.15K). This option is now obsolete because the liquid precipitation is now converted to ice when the temperature is negative using the iflag_bergeron option.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 40.7 KB

Line
1	!
2	! $Id: fisrtilp.F90 2507 2016-05-04 14:36:48Z emillour $
3	!
4	!
5	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6	d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, &
7	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
8	frac_impa, frac_nucl, beta, &
9	prfl, psfl, rhcl, zqta, fraca, &
10	ztv, zpspsk, ztla, zthl, iflag_cld_th, &
11	iflag_ice_thermo)
12
13	!
14	USE dimphy
15	USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
16	USE print_control_mod, ONLY: prt_level, lunout
17	IMPLICIT none
18	!======================================================================
19	! Auteur(s): Z.X. Li (LMD/CNRS)
20	! Date: le 20 mars 1995
21	! Objet: condensation et precipitation stratiforme.
22	! schema de nuage
23	! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
24	! et ilp (il pleut, L. Li)
25	! Principales parties:
26	! P1> Evaporation de la precipitation (qui vient du niveau k+1)
27	! P2> Formation du nuage (en k)
28	! P3> Formation de la precipitation (en k)
29	!======================================================================
30	!======================================================================
31	include "YOMCST.h"
32	include "fisrtilp.h"
33	include "nuage.h" ! JBM (3/14)
34
35	!
36	! Principaux inputs:
37	!
38	REAL dtime ! intervalle du temps (s)
39	REAL paprs(klon,klev+1) ! pression a inter-couche
40	REAL pplay(klon,klev) ! pression au milieu de couche
41	REAL t(klon,klev) ! temperature (K)
42	REAL q(klon,klev) ! humidite specifique (kg/kg)
43	!
44	! Principaux outputs:
45	!
46	REAL d_t(klon,klev) ! incrementation de la temperature (K)
47	REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
48	REAL d_ql(klon,klev) ! incrementation de l'eau liquide
49	REAL d_qi(klon,klev) ! incrementation de l'eau glace
50	REAL rneb(klon,klev) ! fraction nuageuse
51	REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
52	REAL rhcl(klon,klev) ! humidite relative en ciel clair
53	REAL rain(klon) ! pluies (mm/s)
54	REAL snow(klon) ! neige (mm/s)
55	REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
56	REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
57	!
58	! Autres arguments
59	!
60	REAL ztv(klon,klev)
61	REAL zqta(klon,klev),fraca(klon,klev)
62	REAL sigma1(klon,klev),sigma2(klon,klev)
63	REAL qltot(klon,klev),ctot(klon,klev)
64	REAL zpspsk(klon,klev),ztla(klon,klev)
65	REAL zthl(klon,klev)
66	REAL ztfondue, qsl, qsi
67
68	logical lognormale(klon)
69	logical ice_thermo
70
71	!AA
72	! Coeffients de fraction lessivee : pour OFF-LINE
73	!
74	REAL pfrac_nucl(klon,klev)
75	REAL pfrac_1nucl(klon,klev)
76	REAL pfrac_impa(klon,klev)
77	!
78	! Fraction d'aerosols lessivee par impaction et par nucleation
79	! POur ON-LINE
80	!
81	REAL frac_impa(klon,klev)
82	REAL frac_nucl(klon,klev)
83	real zct ,zcl
84	!AA
85	!
86	! Options du programme:
87	!
88	REAL seuil_neb ! un nuage existe vraiment au-dela
89	PARAMETER (seuil_neb=0.001)
90
91	INTEGER ninter ! sous-intervals pour la precipitation
92	INTEGER ncoreczq
93	INTEGER iflag_cld_th
94	INTEGER iflag_ice_thermo
95	PARAMETER (ninter=5)
96	LOGICAL evap_prec ! evaporation de la pluie
97	PARAMETER (evap_prec=.TRUE.)
98	REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
99	logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
100
101	real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
102	real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
103	real erf
104	REAL qcloud(klon)
105	!
106	LOGICAL cpartiel ! condensation partielle
107	PARAMETER (cpartiel=.TRUE.)
108	REAL t_coup
109	PARAMETER (t_coup=234.0)
110	!
111	! Variables locales:
112	!
113	INTEGER i, k, n, kk
114	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
115	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
116	LOGICAL convergence(klon)
117	REAL DDT0
118	PARAMETER (DDT0=.01)
119	INTEGER n_i(klon), iter
120	REAL cste
121
122	REAL zrfl(klon), zrfln(klon), zqev, zqevt
123	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
124	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
125	REAL zoliqp(klon), zoliqi(klon)
126	REAL zt(klon)
127	! JBM (3/14) nexpo is replaced by exposant_glace
128	! REAL nexpo ! exponentiel pour glace/eau
129	! INTEGER, PARAMETER :: nexpo=6
130	INTEGER exposant_glace_old
131	REAL t_glace_min_old
132	REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
133	REAL zchau ,zfroi ,zfice(klon),zneb(klon)
134	REAL zmelt, zpluie, zice, zcondold
135	PARAMETER (ztfondue=278.15)
136	REAL dzfice(klon)
137	REAL zsolid
138	!!!!
139	! Variables pour Bergeron
140	REAL zcp, coef1, DeltaT
141	REAL zqpreci(klon), zqprecl(klon)
142	!
143	LOGICAL appel1er
144	SAVE appel1er
145	!$OMP THREADPRIVATE(appel1er)
146	!
147	!---------------------------------------------------------------
148	!
149	!AA Variables traceurs:
150	!AA Provisoire !!! Parametres alpha du lessivage
151	!AA A priori on a 4 scavenging # possibles
152	!
153	REAL a_tr_sca(4)
154	save a_tr_sca
155	!$OMP THREADPRIVATE(a_tr_sca)
156	!
157	! Variables intermediaires
158	!
159	REAL zalpha_tr
160	REAL zfrac_lessi
161	REAL zprec_cond(klon)
162	!AA
163	! RomP >>> 15 nov 2012
164	REAL beta(klon,klev) ! taux de conversion de l'eau cond
165	! RomP <<<
166	REAL zmair, zcpair, zcpeau
167	! Pour la conversion eau-neige
168	REAL zlh_solid(klon), zm_solid
169	!---------------------------------------------------------------
170	!
171	! Fonctions en ligne:
172	!
173	REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
174	! (Heymsfield & Donner, 1990)
175	REAL zzz
176	include "YOETHF.h"
177	include "FCTTRE.h"
178	fallvc (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_con
179	fallvs (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_lsc
180	!
181	DATA appel1er /.TRUE./
182	!ym
183	!CR: pour iflag_ice_thermo=2, on active que la convection
184	! ice_thermo = iflag_ice_thermo .GE. 1
185	ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
186	zdelq=0.0
187
188	if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
189	IF (appel1er) THEN
190	!
191	WRITE(lunout,*) 'fisrtilp, ninter:', ninter
192	WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
193	WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
194	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
195	WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
196	WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
197	! CALL abort
198	ENDIF
199	appel1er = .FALSE.
200	!
201	!AA initialiation provisoire
202	a_tr_sca(1) = -0.5
203	a_tr_sca(2) = -0.5
204	a_tr_sca(3) = -0.5
205	a_tr_sca(4) = -0.5
206	!
207	!AA Initialisation a 1 des coefs des fractions lessivees
208	!
209	!cdir collapse
210	DO k = 1, klev
211	DO i = 1, klon
212	pfrac_nucl(i,k)=1.
213	pfrac_1nucl(i,k)=1.
214	pfrac_impa(i,k)=1.
215	beta(i,k)=0. !RomP initialisation
216	ENDDO
217	ENDDO
218
219	ENDIF ! test sur appel1er
220	!
221	!MAf Initialisation a 0 de zoliq
222	! DO i = 1, klon
223	! zoliq(i)=0.
224	! ENDDO
225	! Determiner les nuages froids par leur temperature
226	! nexpo regle la raideur de la transition eau liquide / eau glace.
227	!
228	!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
229	IF (iflag_t_glace.EQ.0) THEN
230	! ztglace = RTT - 15.0
231	t_glace_min_old = RTT - 15.0
232	!AJ<
233	IF (ice_thermo) THEN
234	! nexpo = 2
235	exposant_glace_old = 2
236	ELSE
237	! nexpo = 6
238	exposant_glace_old = 6
239	ENDIF
240
241	ENDIF
242
243	!! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
244	!! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
245	!>AJ
246	!cc nexpo = 1
247	!
248	! Initialiser les sorties:
249	!
250	!cdir collapse
251	DO k = 1, klev+1
252	DO i = 1, klon
253	prfl(i,k) = 0.0
254	psfl(i,k) = 0.0
255	ENDDO
256	ENDDO
257
258	!cdir collapse
259	DO k = 1, klev
260	DO i = 1, klon
261	d_t(i,k) = 0.0
262	d_q(i,k) = 0.0
263	d_ql(i,k) = 0.0
264	d_qi(i,k) = 0.0
265	rneb(i,k) = 0.0
266	radliq(i,k) = 0.0
267	frac_nucl(i,k) = 1.
268	frac_impa(i,k) = 1.
269	ENDDO
270	ENDDO
271	DO i = 1, klon
272	rain(i) = 0.0
273	snow(i) = 0.0
274	zoliq(i)=0.
275	! ENDDO
276	!
277	! Initialiser le flux de precipitation a zero
278	!
279	! DO i = 1, klon
280	zrfl(i) = 0.0
281	zifl(i) = 0.0
282	zneb(i) = seuil_neb
283	ENDDO
284	!
285	!
286	!AA Pour plus de securite
287
288	zalpha_tr = 0.
289	zfrac_lessi = 0.
290
291	!AA==================================================================
292	!
293	ncoreczq=0
294	! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
295	!
296	DO k = klev, 1, -1
297	!
298	!AA===============================================================
299	!
300	! Initialisation temperature et vapeur
301	DO i = 1, klon
302	zt(i)=t(i,k)
303	zq(i)=q(i,k)
304	ENDDO
305	!
306	! Calculer la varition de temp. de l'air du a la chaleur sensible
307	! transporter par la pluie.
308	! Il resterait a rajouter cet effet de la chaleur sensible sur les
309	! flux de surface, du a la diff. de temp. entre le 1er niveau et la
310	! surface.
311	!
312	IF(k.LE.klevm1) THEN
313	DO i = 1, klon
314	!IM
315	zmair=(paprs(i,k)-paprs(i,k+1))/RG
316	zcpair=RCPD(1.0+RVTMP2zq(i))
317	zcpeau=RCPD*RVTMP2
318	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zrfl(i)dtime*zcpeau &
319	+ zmairzcpairzt(i) ) &
320	/ (zmairzcpair + zrfl(i)dtime*zcpeau)
321	! C WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
322	ENDDO
323	ENDIF
324	! ----------------------------------------------------------------
325	! P1> Debut evaporation de la precipitation
326	! ----------------------------------------------------------------
327	IF (evap_prec) THEN
328	DO i = 1, klon
329	!AJ<
330	!! IF (zrfl(i) .GT.0.) THEN
331	IF (zrfl(i)+zifl(i).GT.0.) THEN
332	!>AJ
333	! Calcul du qsat
334	IF (thermcep) THEN
335	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
336	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
337	zqs(i)=MIN(0.5,zqs(i))
338	zcor=1./(1.-RETV*zqs(i))
339	zqs(i)=zqs(i)*zcor
340	ELSE
341	IF (zt(i) .LT. t_coup) THEN
342	zqs(i) = qsats(zt(i)) / pplay(i,k)
343	ELSE
344	zqs(i) = qsatl(zt(i)) / pplay(i,k)
345	ENDIF
346	ENDIF
347	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
348	ENDDO
349	!AJ<
350	IF (.NOT. ice_thermo) THEN
351	DO i = 1, klon
352	!AJ<
353	!! IF (zrfl(i) .GT.0.) THEN
354	IF (zrfl(i)+zifl(i).GT.0.) THEN
355	!>AJ
356	! Evap max pour ne pas saturer la fraction sous le nuage
357	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
358	! Calcul de l'evaporation du flux de precip herite
359	! d'au-dessus
360	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
361	* (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
362	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
363	* RG*dtime/(paprs(i,k)-paprs(i,k+1))
364	! Seuil pour ne pas saturer la fraction sous le nuage
365	zqev = MIN (zqev, zqevt)
366	! Nouveau flux de precip
367	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
368	/RG/dtime
369	! Aucun flux liquide pour T < t_coup
370	IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
371	! Nouvelle vapeur
372	zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
373	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
374	! Nouvelle temperature (chaleur latente)
375	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
376	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
377	* RLVTT/RCPD/(1.0+RVTMP2*zq(i))
378	zrfl(i) = zrfln(i)
379	zifl(i) = 0.
380	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
381	ENDDO
382	!
383	ELSE ! (.NOT. ice_thermo)
384	!
385	DO i = 1, klon
386	!AJ<
387	!! IF (zrfl(i) .GT.0.) THEN
388	IF (zrfl(i)+zifl(i).GT.0.) THEN
389	!>AJ
390	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
391	! Modification de l'evaporation avec la glace
392	! Differentiation entre precipitation liquide et solide
393	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
394
395	! Evap max pour ne pas saturer la fraction sous le nuage
396	zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
397	! zqev0 = MAX (0.0, zqs(i)-zq(i) )
398
399	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
400	! On differencie qsat pour l'eau et la glace
401	! Si zdelta=1. --> glace
402	! Si zdelta=0. --> eau liquide
403	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
404
405	! Calcul du qsat par rapport a l'eau liquide
406	qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
407	qsl= MIN(0.5,qsl)
408	zcor= 1./(1.-RETV*qsl)
409	qsl= qsl*zcor
410
411	! Calcul de l'evaporation du flux de precip herite
412	! d'au-dessus
413	! Formulation en racine du flux de precip
414	! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
415	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
416	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
417	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
418	RGdtime/(paprs(i,k)-paprs(i,k+1))
419
420
421	! Calcul du qsat par rapport a la glace
422	qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
423	qsi= MIN(0.5,qsi)
424	zcor= 1./(1.-RETV*qsi)
425	qsi= qsi*zcor
426
427	! Calcul de la sublimation du flux de precip solide herite
428	! d'au-dessus
429	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
430	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
431	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
432	RGdtime/(paprs(i,k)-paprs(i,k+1))
433
434
435	!JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
436	! Verification sur l'evaporation
437	! On s'assure qu'on ne sature pas
438	! la fraction sous le nuage sinon on
439	! repartit zqev0 en gardant la proportion
440	! liquide / glace
441	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
442
443	IF (zqevt+zqevti.GT.zqev0) THEN
444	zqev=zqev0*zqevt/(zqevt+zqevti)
445	zqevi=zqev0*zqevti/(zqevt+zqevti)
446
447	ELSE
448	IF (zqevt+zqevti.GT.0.) THEN
449	zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
450	zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
451	ELSE
452	zqev=0.
453	zqevi=0.
454	ENDIF
455	ENDIF
456	! Nouveaux flux de precip liquide et solide
457	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
458	/RG/dtime)
459	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
460	/RG/dtime)
461
462	! Mise a jour de la vapeur, temperature et flux de precip
463	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
464	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
465	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
466	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
467	* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
468	+ (zifln(i)-zifl(i)) &
469	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
470	* RLSTT/RCPD/(1.0+RVTMP2*zq(i))
471
472	zrfl(i) = zrfln(i)
473	zifl(i) = zifln(i)
474
475	!CR ATTENTION: deplacement de la fonte de la glace
476	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
477	!!! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 !!!!!!!!! jyg
478	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
479	zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! jyg
480	zmelt = MIN(MAX(zmelt,0.),1.)
481	! Fusion de la glace
482	zrfl(i)=zrfl(i)+zmelt*zifl(i)
483	zifl(i)=zifl(i)*(1.-zmelt)
484	! print*,zt(i),'octavio1'
485	! Chaleur latente de fusion
486	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
487	RLMLT/RCPD/(1.0+RVTMP2zq(i))
488	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
489	!fin CR
490
491
492
493	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
494	ENDDO
495
496	ENDIF ! (.NOT. ice_thermo)
497
498	! ----------------------------------------------------------------
499	! Fin evaporation de la precipitation
500	! ----------------------------------------------------------------
501	ENDIF ! (evap_prec)
502	!
503	! Calculer Qs et L/Cp*dQs/dT:
504	!
505	IF (thermcep) THEN
506	DO i = 1, klon
507	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
508	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
509	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
510	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
511	zqs(i) = MIN(0.5,zqs(i))
512	zcor = 1./(1.-RETV*zqs(i))
513	zqs(i) = zqs(i)*zcor
514	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
515	ENDDO
516	ELSE
517	DO i = 1, klon
518	IF (zt(i).LT.t_coup) THEN
519	zqs(i) = qsats(zt(i))/pplay(i,k)
520	zdqs(i) = dqsats(zt(i),zqs(i))
521	ELSE
522	zqs(i) = qsatl(zt(i))/pplay(i,k)
523	zdqs(i) = dqsatl(zt(i),zqs(i))
524	ENDIF
525	ENDDO
526	ENDIF
527	!
528	! Determiner la condensation partielle et calculer la quantite
529	! de l'eau condensee:
530	!
531	!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
532	! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
533	! write(,) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
534	! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
535	! stop
536	! endif
537
538	! ----------------------------------------------------------------
539	! P2> Formation du nuage
540	! ----------------------------------------------------------------
541	IF (cpartiel) THEN
542
543	! print*,'Dans partiel k=',k
544	!
545	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
546	! nuageuse a partir des PDF de Sandrine Bony.
547	! rneb : fraction nuageuse
548	! zqn : eau totale dans le nuage
549	! zcond : eau condensee moyenne dans la maille.
550	! on prend en compte le réchauffement qui diminue la partie
551	! condensee
552	!
553	! Version avec les raqts
554
555	if (iflag_pdf.eq.0) then
556
557	! version creneau de (Li, 1998)
558	do i=1,klon
559	zdelq = min(ratqs(i,k),0.99) * zq(i)
560	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
561	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
562	enddo
563
564	else
565	!
566	! Version avec les nouvelles PDFs.
567	do i=1,klon
568	if(zq(i).lt.1.e-15) then
569	ncoreczq=ncoreczq+1
570	zq(i)=1.e-15
571	endif
572	enddo
573
574	if (iflag_cld_th>=5) then
575
576	call cloudth(klon,klev,k,ztv, &
577	zq,zqta,fraca, &
578	qcloud,ctot,zpspsk,paprs,ztla,zthl, &
579	ratqs,zqs,t)
580
581	do i=1,klon
582	rneb(i,k)=ctot(i,k)
583	zqn(i)=qcloud(i)
584	enddo
585
586	endif
587
588	if (iflag_cld_th <= 4) then
589	lognormale = .true.
590	elseif (iflag_cld_th >= 6) then
591	! lognormale en l'absence des thermiques
592	lognormale = fraca(:,k) < 1e-10
593	else
594	! Dans le cas iflag_cld_th=5, on prend systématiquement la
595	! bi-gaussienne
596	lognormale = .false.
597	end if
598
599	!CR: variation de qsat avec T en presence de glace ou non
600	!initialisations
601	do i=1,klon
602	DT(i) = 0.
603	n_i(i)=0
604	Tbef(i)=zt(i)
605	qlbef(i)=0.
606	enddo
607
608
609	!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
610	if (iflag_fisrtilp_qsat.ge.0) then
611	! Iteration pour condensation avec variation de qsat(T)
612	! -----------------------------------------------------
613	do iter=1,iflag_fisrtilp_qsat+1
614
615	do i=1,klon
616	! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
617	convergence(i)=abs(DT(i)).gt.DDT0
618	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
619	Tbef(i)=Tbef(i)+DT(i)
620	if (.not.ice_thermo) then
621	zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
622	else
623	if (iflag_t_glace.eq.0) then
624	zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
625	else if (iflag_t_glace.ge.1) then
626	zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
627	endif
628	endif
629	! Calcul des PDF lognormales
630	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
631	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
632	zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
633	zqs(i) = MIN(0.5,zqs(i))
634	zcor = 1./(1.-RETV*zqs(i))
635	zqs(i) = zqs(i)*zcor
636	zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
637	zpdf_sig(i)=ratqs(i,k)*zq(i)
638	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
639	zpdf_delta(i)=log(zq(i)/zqs(i))
640	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
641	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
642	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
643	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
644	zpdf_e1(i)=1.-erf(zpdf_e1(i))
645	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
646	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
647	zpdf_e2(i)=1.-erf(zpdf_e2(i))
648
649	if (zpdf_e1(i).lt.1.e-10) then
650	rneb(i,k)=0.
651	zqn(i)=zqs(i)
652	else
653	rneb(i,k)=0.5*zpdf_e1(i)
654	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
655	endif
656
657	endif !convergence
658	enddo ! boucle en i
659
660	if (.not. ice_thermo) then
661
662	do i=1,klon
663	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
664
665	qlbef(i)=max(0.,zqn(i)-zqs(i))
666	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
667	denom=1.+rneb(i,k)*zdqs(i)
668	DT(i)=num/denom
669	n_i(i)=n_i(i)+1
670	endif
671	enddo
672
673	else
674	! Iteration pour convergence avec qsat(T)
675	if (iflag_t_glace.ge.1) then
676	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
677	endif
678
679	do i=1,klon
680	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
681
682	if (iflag_t_glace.eq.0) then
683	zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
684	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
685	zfice(i) = zfice(i)**exposant_glace_old
686	dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
687	endif
688
689	if (iflag_t_glace.ge.1) then
690	dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
691	endif
692
693	if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
694	dzfice(i)=0.
695	endif
696
697	if (zfice(i).lt.1) then
698	cste=RLVTT
699	else
700	cste=RLSTT
701	endif
702
703	qlbef(i)=max(0.,zqn(i)-zqs(i))
704	num=-Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
705	denom=1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
706	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))rneb(i,k)qlbef(i)dzfice(i)
707	DT(i)=num/denom
708	n_i(i)=n_i(i)+1
709
710	endif ! fin convergence
711	enddo ! fin boucle i
712
713	endif !ice_thermo
714
715	! endif
716	! enddo
717
718
719	enddo ! iter=1,iflag_fisrtilp_qsat+1
720	! Fin d'iteration pour condensation avec variation de qsat(T)
721	! -----------------------------------------------------------
722	endif
723
724
725	endif ! iflag_pdf
726
727
728	! if (iflag_fisrtilp_qsat.eq.-1) then
729	!------------------------------------------
730	!CR: ATTENTION option fausse mais a existe:
731	! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
732	! activer les lignes suivantes:
733	IF (1.eq.0) THEN
734	DO i=1,klon
735	IF (rneb(i,k) .LE. 0.0) THEN
736	zqn(i) = 0.0
737	rneb(i,k) = 0.0
738	zcond(i) = 0.0
739	rhcl(i,k)=zq(i)/zqs(i)
740	ELSE IF (rneb(i,k) .GE. 1.0) THEN
741	zqn(i) = zq(i)
742	rneb(i,k) = 1.0
743	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
744	rhcl(i,k)=1.0
745	ELSE
746	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
747	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
748	ENDIF
749	ENDDO
750	ENDIF
751	!------------------------------------------
752
753	! ELSE
754
755	! Calcul de l'eau in-cloud (zqn),
756	! moyenne dans la maille (zcond),
757	! fraction nuageuse (rneb) et
758	! humidite relative ciel-clair (rhcl)
759	DO i=1,klon
760	IF (rneb(i,k) .LE. 0.0) THEN
761	zqn(i) = 0.0
762	rneb(i,k) = 0.0
763	zcond(i) = 0.0
764	rhcl(i,k)=zq(i)/zqs(i)
765	ELSE IF (rneb(i,k) .GE. 1.0) THEN
766	zqn(i) = zq(i)
767	rneb(i,k) = 1.0
768	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
769	rhcl(i,k)=1.0
770	ELSE
771	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
772	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
773	ENDIF
774	ENDDO
775
776
777	! ENDIF
778
779	ELSE ! de IF (cpartiel)
780	! Cas "tout ou rien"
781	DO i = 1, klon
782	IF (zq(i).GT.zqs(i)) THEN
783	rneb(i,k) = 1.0
784	ELSE
785	rneb(i,k) = 0.0
786	ENDIF
787	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
788	ENDDO
789	ENDIF
790	! ----------------------------------------------------------------
791	! Fin de formation du nuage
792	! ----------------------------------------------------------------
793	!
794	! Mise a jour vapeur d'eau
795	DO i = 1, klon
796	zq(i) = zq(i) - zcond(i)
797	! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
798	ENDDO
799	!AJ<
800	! Chaleur latente apres formation nuage
801	! -------------------------------------
802	IF (.NOT. ice_thermo) THEN
803	if (iflag_fisrtilp_qsat.lt.1) then
804	DO i = 1, klon
805	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
806	ENDDO
807	else if (iflag_fisrtilp_qsat.gt.0) then
808	DO i= 1, klon
809	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
810	ENDDO
811	endif
812	ELSE
813	if (iflag_t_glace.ge.1) then
814	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
815	endif
816	if (iflag_fisrtilp_qsat.lt.1) then
817	DO i = 1, klon
818	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
819	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
820	! t_glace_max, exposant_glace)
821	if (iflag_t_glace.eq.0) then
822	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
823	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
824	zfice(i) = zfice(i)**exposant_glace_old
825	endif
826	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
827	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
828	ENDDO
829	else
830	DO i=1, klon
831	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
832	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
833	! t_glace_max, exposant_glace)
834	if (iflag_t_glace.eq.0) then
835	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
836	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
837	zfice(i) = zfice(i)**exposant_glace_old
838	endif
839	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
840	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
841	ENDDO
842	endif
843	! print*,zt(i),zrfl(i),zifl(i),'temp1'
844	ENDIF
845	!>AJ
846	! ----------------------------------------------------------------
847	! P3> Formation des precipitations
848	! ----------------------------------------------------------------
849	!
850	! Partager l'eau condensee en precipitation et eau liquide nuageuse
851	!
852
853	! Initialisation de zoliq (eau condensee moyenne dans la maille)
854	DO i = 1, klon
855	IF (rneb(i,k).GT.0.0) THEN
856	zoliq(i) = zcond(i)
857	zrho(i) = pplay(i,k) / zt(i) / RD
858	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
859	ENDIF
860	ENDDO
861	!AJ<
862	IF (.NOT. ice_thermo) THEN
863	IF (iflag_t_glace.EQ.0) THEN
864	DO i = 1, klon
865	IF (rneb(i,k).GT.0.0) THEN
866	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
867	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
868	zfice(i) = zfice(i)**exposant_glace_old
869	! zfice(i) = zfice(i)**nexpo
870	!! zfice(i)=0.
871	ENDIF
872	ENDDO
873	ELSE ! of IF (iflag_t_glace.EQ.0)
874	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
875	! DO i = 1, klon
876	! IF (rneb(i,k).GT.0.0) THEN
877	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
878	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
879	! t_glace_max, exposant_glace)
880	! ENDIF
881	! ENDDO
882	ENDIF
883	ENDIF
884
885	! Calcul de radliq (eau condensee pour le rayonnement)
886	! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
887	! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
888	! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
889	! transmise au rayonnement;
890	! ----------------------------------------------------------------
891	DO i = 1, klon
892	IF (rneb(i,k).GT.0.0) THEN
893	zneb(i) = MAX(rneb(i,k), seuil_neb)
894	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
895	! print*,zt(i),'fractionglace'
896	!>AJ
897	radliq(i,k) = zoliq(i)/REAL(ninter+1)
898	ENDIF
899	ENDDO
900	!
901	DO n = 1, ninter
902	DO i = 1, klon
903	IF (rneb(i,k).GT.0.0) THEN
904	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
905	! Initialization of zpluie and zice:
906	zpluie=0
907	zice=0
908	IF (zneb(i).EQ.seuil_neb) THEN
909	ztot = 0.0
910	ELSE
911	! quantite d'eau a eliminer: zchau (Sundqvist, 1978)
912	! meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
913	if (ptconv(i,k)) then
914	zcl =cld_lc_con
915	zct =1./cld_tau_con
916	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
917	fallvc(zrhol(i)) zfice(i)
918	else
919	zcl =cld_lc_lsc
920	zct =1./cld_tau_lsc
921	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
922	fallvs(zrhol(i)) zfice(i)
923	endif
924	zchau = zct dtime/REAL(ninter) zoliq(i) &
925	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
926	!AJ<
927	IF (.NOT. ice_thermo) THEN
928	ztot = zchau + zfroi
929	ELSE
930	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
931	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
932	ztot = zpluie + zice
933	ENDIF
934	!>AJ
935	ztot = MAX(ztot ,0.0)
936	ENDIF
937	ztot = MIN(ztot,zoliq(i))
938	!AJ<
939	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
940	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
941	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
942	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
943	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
944	!>AJ
945	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
946	ENDIF
947	ENDDO ! i = 1,klon
948	ENDDO ! n = 1,ninter
949	! ----------------------------------------------------------------
950	!
951	IF (.NOT. ice_thermo) THEN
952	DO i = 1, klon
953	IF (rneb(i,k).GT.0.0) THEN
954	d_ql(i,k) = zoliq(i)
955	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
956	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
957	ENDIF
958	ENDDO
959	ELSE
960	!
961	!CR&JYG<
962	! On prend en compte l'effet Bergeron dans les flux de precipitation :
963	! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
964	! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
965	! et les precipitations est grossierement pris en compte en linearisant les equations
966	! et en approximant le processus de precipitation liquide par un processus a seuil.
967	! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
968	! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
969	! Le condensat precipitant solide est augmente.
970	! La vapeur d'eau est augmentee.
971	!
972	IF ((iflag_bergeron .EQ. 2)) THEN
973	DO i = 1, klon
974	IF (rneb(i,k) .GT. 0.0) THEN
975	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
976	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
977	zcp=RCPD(1.0+RVTMP2(zq(i)+zcond(i)))
978	coef1 = RLMLT*zdqs(i)/RLVTT
979	DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
980	zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
981	zqprecl(i) = max( zqprecl(i) - zcp/RLMLT(1.+coef1)DeltaT, 0. )
982	zcond(i) = max( zcond(i) - zcp/RLVTTzdqs(i)DeltaT, 0. )
983	zq(i) = zq(i) + zcp/RLVTTzdqs(i)DeltaT
984	zt(i) = zt(i) + DeltaT
985	ENDIF ! rneb(i,k) .GT. 0.0
986	ENDDO
987	DO i = 1, klon
988	IF (rneb(i,k).GT.0.0) THEN
989	d_ql(i,k) = (1-zfice(i))*zoliq(i)
990	d_qi(i,k) = zfice(i)*zoliq(i)
991	zrfl(i) = zrfl(i)+ zqprecl(i) &
992	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
993	zifl(i) = zifl(i)+ zqpreci(i) &
994	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
995	ENDIF
996	ENDDO
997	!!
998	ELSE ! iflag_bergeron
999	!>CR&JYG
1000	!!
1001	DO i = 1, klon
1002	IF (rneb(i,k).GT.0.0) THEN
1003	!CR on prend en compte la phase glace
1004	if (.not.ice_thermo) then
1005	d_ql(i,k) = zoliq(i)
1006	d_qi(i,k) = 0.
1007	else
1008	d_ql(i,k) = (1-zfice(i))*zoliq(i)
1009	d_qi(i,k) = zfice(i)*zoliq(i)
1010	endif
1011	!AJ<
1012	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1013	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1014	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1015	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1016	! zrfl(i) = zrfl(i)+ zpluie &
1017	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
1018	! zifl(i) = zifl(i)+ zice &
1019	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
1020
1021	!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
1022	IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
1023	zsolid = zrfl(i)
1024	zifl(i) = zifl(i)+zrfl(i)
1025	zrfl(i) = 0.
1026	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
1027	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))
1028	ENDIF ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
1029	!RC
1030
1031	ENDIF ! rneb(i,k).GT.0.0
1032	ENDDO
1033
1034	ENDIF ! iflag_bergeron .EQ. 2
1035	ENDIF ! .NOT. ice_thermo
1036
1037	!CR: la fonte est faite au debut
1038	! IF (ice_thermo) THEN
1039	! DO i = 1, klon
1040	! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1041	! zmelt = MIN(MAX(zmelt,0.),1.)
1042	! zrfl(i)=zrfl(i)+zmelt*zifl(i)
1043	! zifl(i)=zifl(i)*(1.-zmelt)
1044	! print*,zt(i),'octavio1'
1045	! zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1046	! RLMLT/RCPD/(1.0+RVTMP2zq(i))
1047	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
1048	! ENDDO
1049	! ENDIF
1050
1051
1052	IF (.NOT. ice_thermo) THEN
1053	DO i = 1, klon
1054	IF (zt(i).LT.RTT) THEN
1055	psfl(i,k)=zrfl(i)
1056	ELSE
1057	prfl(i,k)=zrfl(i)
1058	ENDIF
1059	ENDDO
1060	ELSE
1061	! JAM*************************************************
1062	! Revoir partie ci-dessous: a quoi servent psfl et prfl?
1063	! *****************************************************
1064
1065	DO i = 1, klon
1066	! IF (zt(i).LT.RTT) THEN
1067	psfl(i,k)=zifl(i)
1068	! ELSE
1069	prfl(i,k)=zrfl(i)
1070	! ENDIF
1071	!>AJ
1072	ENDDO
1073	ENDIF
1074	! ----------------------------------------------------------------
1075	! Fin de formation des precipitations
1076	! ----------------------------------------------------------------
1077	!
1078	!
1079	! Calculer les tendances de q et de t:
1080	!
1081	DO i = 1, klon
1082	d_q(i,k) = zq(i) - q(i,k)
1083	d_t(i,k) = zt(i) - t(i,k)
1084	ENDDO
1085	!
1086	!AA--------------- Calcul du lessivage stratiforme -------------
1087
1088	DO i = 1,klon
1089	!
1090	if(zcond(i).gt.zoliq(i)+1.e-10) then
1091	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1092	else
1093	beta(i,k) = 0.
1094	endif
1095	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1096	* (paprs(i,k)-paprs(i,k+1))/RG
1097	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1098	!AA lessivage nucleation LMD5 dans la couche elle-meme
1099	IF (iflag_t_glace.EQ.0) THEN
1100	if (t(i,k) .GE. t_glace_min_old) THEN
1101	zalpha_tr = a_tr_sca(3)
1102	else
1103	zalpha_tr = a_tr_sca(4)
1104	endif
1105	ELSE ! of IF (iflag_t_glace.EQ.0)
1106	if (t(i,k) .GE. t_glace_min) THEN
1107	zalpha_tr = a_tr_sca(3)
1108	else
1109	zalpha_tr = a_tr_sca(4)
1110	endif
1111	ENDIF
1112	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1113	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
1114	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1115	!
1116	! nucleation avec un facteur -1 au lieu de -0.5
1117	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1118	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
1119	ENDIF
1120	!
1121	ENDDO ! boucle sur i
1122	!
1123	!AA Lessivage par impaction dans les couches en-dessous
1124	DO kk = k-1, 1, -1
1125	DO i = 1, klon
1126	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1127	IF (iflag_t_glace.EQ.0) THEN
1128	if (t(i,kk) .GE. t_glace_min_old) THEN
1129	zalpha_tr = a_tr_sca(1)
1130	else
1131	zalpha_tr = a_tr_sca(2)
1132	endif
1133	ELSE ! of IF (iflag_t_glace.EQ.0)
1134	if (t(i,kk) .GE. t_glace_min) THEN
1135	zalpha_tr = a_tr_sca(1)
1136	else
1137	zalpha_tr = a_tr_sca(2)
1138	endif
1139	ENDIF
1140	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1141	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
1142	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1143	ENDIF
1144	ENDDO
1145	ENDDO
1146	!
1147	!AA===============================================================
1148	! FIN DE LA BOUCLE VERTICALE
1149	end DO
1150	!
1151	!AA==================================================================
1152	!
1153	! Pluie ou neige au sol selon la temperature de la 1ere couche
1154	!
1155	!CR: si la thermo de la glace est active, on calcule zifl directement
1156	IF (.NOT.ice_thermo) THEN
1157	DO i = 1, klon
1158	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1159	!AJ<
1160	! snow(i) = zrfl(i)
1161	snow(i) = zrfl(i)+zifl(i)
1162	!>AJ
1163	zlh_solid(i) = RLSTT-RLVTT
1164	ELSE
1165	rain(i) = zrfl(i)
1166	zlh_solid(i) = 0.
1167	ENDIF
1168	ENDDO
1169
1170	ELSE
1171	DO i = 1, klon
1172	snow(i) = zifl(i)
1173	rain(i) = zrfl(i)
1174	ENDDO
1175
1176	ENDIF
1177	!
1178	! For energy conservation : when snow is present, the solification
1179	! latent heat is considered.
1180	!CR: si thermo de la glace, neige deja prise en compte
1181	IF (.not.ice_thermo) THEN
1182	DO k = 1, klev
1183	DO i = 1, klon
1184	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
1185	zmair=(paprs(i,k)-paprs(i,k+1))/RG
1186	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1187	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair)
1188	END DO
1189	END DO
1190	ENDIF
1191	!
1192
1193	if (ncoreczq>0) then
1194	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1195	endif
1196
1197	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format