Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

fisrtilp.F90 @ 3485

Last change on this file since 3485 was 3331, checked in by acozic, 6 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 99.0 KB

Line
1	!
2	! $Id: fisrtilp.F90 2507 2016-05-04 14:36:48Z jbmadeleine $
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	#ifdef ISO
13	& ,xt,xtrain,xtsnow &
14	& ,d_xt,d_xtl,d_xti,pxtrainfl,pxtsnowfl &
15	#ifdef DIAGISO
16	& ,fcond_ls,taux_cond_ls,taux_precip_ls & ! juste diagnostique
17	#endif
18	#endif
19	& )
20
21	!
22	USE dimphy
23	USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
24	USE print_control_mod, ONLY: prt_level, lunout
25	#ifdef ISO
26	USE infotrac_phy, ONLY: ntraciso,niso,index_trac,ntraceurs_zone
27	USE isotopes_mod, ONLY: essai_convergence,bidouille_anti_divergence, &
28	& Rdefault,ridicule,pxtice,pxtmelt,iso_eau,iso_HDO, &
29	& iso_O18,ridicule_rain
30	USE isotopes_routines_mod, ONLY: iso_revap_fisrtilp, &
31	condiso_liq_ice_vectall
32	#ifdef ISOVERIF
33	USE isotopes_verif_mod, ONLY: errmax,errmaxrel,deltalim_snow,deltalim, &
34	iso_verif_egalite_choix,iso_verif_aberrant_choix, &
35	iso_verif_positif,iso_verif_positif_choix,iso_verif_noNaN, &
36	iso_verif_aberrant,iso_verif_positif_choix_nostop, &
37	iso_verif_aberrant_encadre,iso_verif_egalite_choix_nostop, &
38	iso_verif_aberrant_nostop,iso_verif_o18_aberrant_nostop, &
39	deltaD,deltaO,iso_verif_egalite
40	#endif
41	#ifdef ISOTRAC
42	use isotrac_mod, only: option_traceurs,option_tmin,seuil_tag_tmin_ls, &
43	& izone_poubelle, option_cond,izone_cond,index_iso, index_zone,izone_oce,izone_ddft
44	use isotrac_routines_mod, only: iso_recolorise_condensation
45	use isotopes_routines_mod, only: condiso_liq_ice_vectall_trac
46	#ifdef ISOVERIF
47	USE isotopes_verif_mod, ONLY: iso_verif_traceur_nostop,iso_verif_traceur
48	use isotrac_routines_mod, ONLY: iso_verif_traceur_pbidouille
49	#endif
50	#endif
51	#endif
52	IMPLICIT none
53	!======================================================================
54	! Auteur(s): Z.X. Li (LMD/CNRS)
55	! Date: le 20 mars 1995
56	! Objet: condensation et precipitation stratiforme.
57	! schema de nuage
58	! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
59	! et ilp (il pleut, L. Li)
60	! Principales parties:
61	! P1> Evaporation de la precipitation (qui vient du niveau k+1)
62	! P2> Formation du nuage (en k)
63	! P3> Formation de la precipitation (en k)
64	!======================================================================
65	!======================================================================
66	include "YOMCST.h"
67	include "fisrtilp.h"
68	include "nuage.h" ! JBM (3/14)
69	#ifdef CAMMODIF
70	include "sursat_marine.h"
71	#endif
72
73	!
74	! Principaux inputs:
75	!
76	REAL dtime ! intervalle du temps (s)
77	REAL paprs(klon,klev+1) ! pression a inter-couche
78	REAL pplay(klon,klev) ! pression au milieu de couche
79	REAL t(klon,klev) ! temperature (K)
80	REAL q(klon,klev) ! humidite specifique (kg/kg)
81	!
82	! Principaux outputs:
83	!
84	REAL d_t(klon,klev) ! incrementation de la temperature (K)
85	REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
86	REAL d_ql(klon,klev) ! incrementation de l'eau liquide
87	REAL d_qi(klon,klev) ! incrementation de l'eau glace
88	REAL rneb(klon,klev) ! fraction nuageuse
89	REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
90	REAL rhcl(klon,klev) ! humidite relative en ciel clair
91	REAL rain(klon) ! pluies (mm/s)
92	REAL snow(klon) ! neige (mm/s)
93	REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
94	REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
95	!
96	! Autres arguments
97	!
98	REAL ztv(klon,klev)
99	REAL zqta(klon,klev),fraca(klon,klev)
100	REAL sigma1(klon,klev),sigma2(klon,klev)
101	REAL qltot(klon,klev),ctot(klon,klev)
102	REAL zpspsk(klon,klev),ztla(klon,klev)
103	REAL zthl(klon,klev)
104	REAL ztfondue, qsl, qsi
105
106	logical lognormale(klon)
107	logical ice_thermo
108
109	#ifdef ISO
110	real double_to_real
111	double precision real_to_double
112	!integer niso
113	real xt(ntraciso,klon,klev)
114	real d_xt(ntraciso,klon,klev)
115	real d_xtl(ntraciso,klon,klev)
116	real d_xti(ntraciso,klon,klev)
117	real xtrain(ntraciso,klon)
118	real xtsnow(ntraciso,klon)
119	real pxtrainfl(ntraciso,klon,klev+1)
120	real pxtsnowfl(ntraciso,klon,klev+1)
121	! xt <-> q
122	! xtrain et xtsonw <-> rain et snow
123	! d_xt et d_xtl <-> d_q et d_ql
124	! pxtrainfl et pxtsnowflux <-> prfl et psfl
125	#ifdef DIAGISO
126	real fcond_ls(klon,klev)
127	real taux_cond_ls(klon,klev)
128	real taux_precip_ls(klon,klev)
129	#endif
130	#endif
131
132
133	!AA
134	! Coeffients de fraction lessivee : pour OFF-LINE
135	!
136	REAL pfrac_nucl(klon,klev)
137	REAL pfrac_1nucl(klon,klev)
138	REAL pfrac_impa(klon,klev)
139	!
140	! Fraction d'aerosols lessivee par impaction et par nucleation
141	! POur ON-LINE
142	!
143	REAL frac_impa(klon,klev)
144	REAL frac_nucl(klon,klev)
145	real zct ,zcl
146	!AA
147	!
148	! Options du programme:
149	!
150	REAL seuil_neb ! un nuage existe vraiment au-dela
151	PARAMETER (seuil_neb=0.001)
152
153	INTEGER ninter ! sous-intervals pour la precipitation
154	INTEGER ncoreczq
155	INTEGER iflag_cld_th
156	INTEGER iflag_ice_thermo
157	PARAMETER (ninter=5)
158	LOGICAL evap_prec ! evaporation de la pluie
159	PARAMETER (evap_prec=.TRUE.)
160	REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
161	logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
162
163	real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
164	real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
165	real erf
166	REAL qcloud(klon)
167	!
168	LOGICAL cpartiel ! condensation partielle
169	PARAMETER (cpartiel=.TRUE.)
170	REAL t_coup
171	PARAMETER (t_coup=234.0)
172
173	#ifdef CAMMODIF
174	! pour la sursat de Marine/Bernard
175	real zpdf_beta(klon)
176	real RH_crit, pi
177	REAL rhcl_p0, rhcl_p1, rhcl_p2, rhcl_p
178	real zpdf_delta_o(klon), zpdf_a_o(klon), zpdf_b_o(klon)
179	real zpdf_e1_o(klon)
180	real rneb_o(klon,klev)
181	real flag_ss(klon,klev),clsk_ss(klon,klev),rneb_ss(klon,klev)
182	real flagb_ss(klon,klev)
183	REAL zqs_star(klon), ys(klon), yss(klon)
184	#endif
185	!
186	! Variables locales:
187	!
188	INTEGER i, k, n, kk
189	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
190	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
191	LOGICAL convergence(klon)
192	REAL DDT0
193	PARAMETER (DDT0=.01)
194	INTEGER n_i(klon), iter
195	REAL cste
196
197	REAL zrfl(klon), zrfln(klon), zqev, zqevt
198	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
199	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
200	REAL zoliqp(klon), zoliqi(klon)
201	REAL zt(klon)
202	! JBM (3/14) nexpo is replaced by exposant_glace
203	! REAL nexpo ! exponentiel pour glace/eau
204	! INTEGER, PARAMETER :: nexpo=6
205	INTEGER exposant_glace_old
206	REAL t_glace_min_old
207	REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
208	REAL zchau ,zfroi ,zfice(klon),zneb(klon)
209	REAL zmelt, zpluie, zice, zcondold
210	PARAMETER (ztfondue=278.15)
211	REAL dzfice(klon)
212	REAL zsolid
213	!!!!
214	! Variables pour Bergeron
215	REAL zcp, coef1, DeltaT
216	REAL zqpreci(klon), zqprecl(klon)
217
218	#ifdef ISO
219	real zxtrfl(ntraciso,klon)
220	real zxtrfln(ntraciso,klon)
221	REAL zxtifl(ntraciso,klon), zxtifln(ntraciso,klon)
222	real zxt(ntraciso,klon) ! equivalent local de xt
223	real zxtn(ntraciso,klon) ! isotpes nuageux
224	real zxtn_reduit(niso)
225	real zxtcond(ntraciso,klon)
226	real zxtoliq(ntraciso,klon)
227	real zq_ancien(klon)
228	REAL zxtpreci(ntraciso,klon), zxtprecl(ntraciso,klon)
229	integer ixt
230	#ifdef ISOVERIF
231	! integer iso_verif_aberrant_nostop
232	! integer iso_verif_egalite_choix_nostop
233	! integer iso_verif_positif_nostop
234	! integer iso_verif_positif_choix_nostop
235	integer trace_cas(klon)
236	! integer iso_verif_traceur_nostop
237	#endif
238	integer iso_verif_nonan_nostop
239	real zxtice(ntraciso,klon),zxtliq(ntraciso,klon)
240	real zrfl_ancien(klon)
241	real zqev_diag(klon)
242	real zxtrfl_ancien(ntraciso,klon)
243	real zxt_ancien(ntraciso,klon)
244	#ifdef ISOTRAC
245	real zxtres(ntraciso) ! humidite spec nuageuse après condensation et precip
246	real zqcs(klon) ! humidite spec en ciel clair
247	real zxtcs(ntraciso,klon)
248	real zcondn(klon) ! humidite spec de l'eau liquide dans le nuage
249	real zxtcondn(ntraciso,klon)
250	integer izone
251	#ifdef ISOVERIF
252	real zq_verif
253	#endif
254	#endif
255	! logical negation ! pour nier conditions logiques
256	#endif
257
258	!
259	LOGICAL appel1er
260	SAVE appel1er
261	!$OMP THREADPRIVATE(appel1er)
262	!
263	!---------------------------------------------------------------
264	!
265	!AA Variables traceurs:
266	!AA Provisoire !!! Parametres alpha du lessivage
267	!AA A priori on a 4 scavenging # possibles
268	!
269	REAL a_tr_sca(4)
270	save a_tr_sca
271	!$OMP THREADPRIVATE(a_tr_sca)
272	!
273	! Variables intermediaires
274	!
275	REAL zalpha_tr
276	REAL zfrac_lessi
277	REAL zprec_cond(klon)
278	!AA
279	! RomP >>> 15 nov 2012
280	REAL beta(klon,klev) ! taux de conversion de l'eau cond
281	! RomP <<<
282	REAL zmair, zcpair, zcpeau
283	! Pour la conversion eau-neige
284	REAL zlh_solid(klon), zm_solid
285	!---------------------------------------------------------------
286	!
287	! Fonctions en ligne:
288	!
289	REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
290	! (Heymsfield & Donner, 1990)
291	REAL zzz
292	include "YOETHF.h"
293	include "FCTTRE.h"
294	fallvc (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_con
295	fallvs (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_lsc
296	!
297	DATA appel1er /.TRUE./
298	!ym
299	!CR: pour iflag_ice_thermo=2, on active que la convection
300	! ice_thermo = iflag_ice_thermo .GE. 1
301	ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
302	zdelq=0.0
303
304	if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
305	IF (appel1er) THEN
306	!
307	WRITE(lunout,*) 'fisrtilp, dtime:', dtime
308	WRITE(lunout,*) 'fisrtilp, ninter:', ninter
309	WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
310	WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
311	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
312	WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
313	WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
314	! CALL abort
315	ENDIF
316	appel1er = .FALSE.
317
318	#ifdef ISO
319	! cam verif
320	#ifdef ISOVERIF
321	write(,) 'ilp 310'
322	do k=1,klev
323	do i=1,klon
324	if (iso_eau.gt.0) then
325	call iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &
326	& 'il pleut 205',errmax,errmaxrel)
327	endif !if (iso_eau.gt.0) then
328	if (iso_HDO.gt.0) then
329	if (q(i,k).gt.ridicule) then
330	call iso_verif_aberrant(xt(iso_HDO,i,k)/q(i,k), &
331	& 'il pleut 210')
332	endif !if (zq(i).gt.ridicule) then
333	endif !if (iso_HDO.gt.0) then
334	enddo !do k=1,klev
335	enddo !do i=1,klon
336	do k=1,klev
337	do i=1,klon
338	call iso_verif_positif(370.0-t(i,k),'il pleut 250')
339	call iso_verif_positif(t(i,k)-100.0,'il pleut 251')
340	enddo
341	enddo
342	write(,) 'ilp 331'
343	#endif
344	! end cam verif
345	#endif
346	!
347	!AA initialiation provisoire
348	a_tr_sca(1) = -0.5
349	a_tr_sca(2) = -0.5
350	a_tr_sca(3) = -0.5
351	a_tr_sca(4) = -0.5
352	!
353	!AA Initialisation a 1 des coefs des fractions lessivees
354	!
355	!cdir collapse
356	DO k = 1, klev
357	DO i = 1, klon
358	pfrac_nucl(i,k)=1.
359	pfrac_1nucl(i,k)=1.
360	pfrac_impa(i,k)=1.
361	beta(i,k)=0. !RomP initialisation
362	#ifdef CAMMODIF
363	flag_ss(i,k) = 0.
364	clsk_ss(i,k) = 0.
365	rneb_ss(i,k) = 0.
366	flagb_ss(i,k) = 0.
367	#endif
368	ENDDO
369	ENDDO
370
371	ENDIF ! test sur appel1er
372	!
373	!MAf Initialisation a 0 de zoliq
374	! DO i = 1, klon
375	! zoliq(i)=0.
376	! ENDDO
377	! Determiner les nuages froids par leur temperature
378	! nexpo regle la raideur de la transition eau liquide / eau glace.
379	!
380	!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
381	IF (iflag_t_glace.EQ.0) THEN
382	! ztglace = RTT - 15.0
383	t_glace_min_old = RTT - 15.0
384	!AJ<
385	IF (ice_thermo) THEN
386	! nexpo = 2
387	exposant_glace_old = 2
388	ELSE
389	! nexpo = 6
390	exposant_glace_old = 6
391	ENDIF
392
393	ENDIF
394
395	!! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
396	!! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
397	!>AJ
398	!cc nexpo = 1
399	!
400	! Initialiser les sorties:
401	!
402	!cdir collapse
403	!#ifdef ISOVERIF
404	! write(,) 'ilp 393'
405	!#endif
406	DO k = 1, klev+1
407	DO i = 1, klon
408	prfl(i,k) = 0.0
409	psfl(i,k) = 0.0
410	#ifdef ISO
411	do ixt=1,ntraciso
412	pxtrainfl(ixt,i,k)=0.0
413	pxtsnowfl(ixt,i,k)=0.0
414	enddo !do ixt=1,niso
415	#endif
416	ENDDO
417	ENDDO
418
419	!cdir collapse
420	DO k = 1, klev
421	DO i = 1, klon
422	d_t(i,k) = 0.0
423	d_q(i,k) = 0.0
424	d_ql(i,k) = 0.0
425	d_qi(i,k) = 0.0
426	#ifdef ISO
427	#ifdef DIAGISO
428	fcond_ls(i,k)=0.0
429	#endif
430	do ixt=1,ntraciso
431	d_xt(ixt,i,k) = 0.0
432	d_xtl(ixt,i,k) = 0.0
433	d_xti(ixt,i,k) = 0.0
434	enddo !do ixt=1,niso
435	#endif
436
437	rneb(i,k) = 0.0
438	radliq(i,k) = 0.0
439	frac_nucl(i,k) = 1.
440	frac_impa(i,k) = 1.
441	ENDDO
442	ENDDO
443	DO i = 1, klon
444	rain(i) = 0.0
445	snow(i) = 0.0
446	zoliq(i)=0.
447	#ifdef ISO
448	do ixt=1,ntraciso
449	xtrain(ixt,i) = 0.0
450	xtsnow(ixt,i) = 0.0
451	zxtoliq(ixt,i)=0.
452	enddo !do ixt=1,niso
453	#endif
454	! ENDDO
455	!
456	! Initialiser le flux de precipitation a zero
457	!
458	! DO i = 1, klon
459	zrfl(i) = 0.0
460	zifl(i) = 0.0
461	#ifdef ISO
462	if (essai_convergence) then
463	#ifdef ISOVERIF
464	zrfln(i)=0.0
465	zifln(i)=0.0
466	#endif
467	else
468	zrfln(i)=0.0
469	zifln(i)=0.0
470	endif
471	do ixt=1,ntraciso
472	zxtrfl(ixt,i)=0.0
473	zxtifl(ixt,i)=0.0
474	zxtifln(ixt,i)=0.0
475	enddo
476	#endif
477	zneb(i) = seuil_neb
478	ENDDO
479	!#ifdef ISOVERIF
480	! write(,) 'ilp 465'
481	!#endif
482	!
483	!
484	!AA Pour plus de securite
485
486	zalpha_tr = 0.
487	zfrac_lessi = 0.
488
489	!AA==================================================================
490	!
491	ncoreczq=0
492	! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
493	!
494	DO k = klev, 1, -1
495	!
496	#ifdef ISO
497	! cam debug!
498	! verif en debut de boucle
499	#ifdef ISOVERIF
500	!write(,) 'ilp 478: boucle en k: k=',k
501	do i=1,klon
502	if (iso_eau.gt.0) then
503	call iso_verif_egalite_choix(zxtrfl(iso_eau,i), &
504	& zrfl(i), &
505	& 'il pleut 310',errmax,errmaxrel)
506	endif !if (iso_eau.gt.0) then
507	if (iso_HDO.gt.0) then
508	if (zrfl(i).gt.ridicule_rain) then
509	call iso_verif_aberrant(zxtrfl(iso_HDO,i) &
510	& /zrfl(i),'il pleut 316')
511	endif !if (zrfl(i).gt.ridicule_rain) then
512	endif !if (iso_HDO.gt.0) then
513	#ifdef ISOTRAC
514	call iso_verif_traceur(zxtrfl(1,i), &
515	& 'il pleut 365')
516	call iso_verif_traceur_pbidouille(zxtrfl(1,i), &
517	& 'il pleut 388')
518	#endif
519	enddo !do i=1,klon
520	#endif
521	#ifdef ISOVERIF
522	!write(,) 'ilp tmp 412'
523	do i=1,klon
524	do ixt=1,ntraciso
525	call iso_verif_noNaN(zxtrfl(ixt,i),'ilp 406')
526	call iso_verif_noNaN(zxtoliq(ixt,i),'ilp 407')
527	enddo !do ixt=1,ntraciso
528	enddo !do i=1,klon
529	#endif
530	! end verif en debut de boucle
531	! end cam debug
532	#endif
533	!AA===============================================================
534	!
535	! Initialisation temperature et vapeur
536	DO i = 1, klon
537	zt(i)=t(i,k)
538	zq(i)=q(i,k)
539	#ifdef ISOVERIF
540	if (q(i,k).lt.0.0) then
541	write(,) 'ilp 400: avant clmain 31janv: i,k,q=', &
542	& i,k,q(i,k)
543	CALL abort_physic('ilp 484', 'on stop', 1)
544	endif
545	#endif
546	#ifdef ISO
547	if (essai_convergence) then
548	else
549	q(i,k)=max(0.0,q(i,k))
550	zq(i)=max(0.0,zq(i))
551	endif
552	#ifdef ISOTRAC
553	#ifdef ISOVERIF
554	call iso_verif_traceur(xt(1,i,k),'ilp 404')
555	! write(,) 'ilp 402 temporaire: xt(1:ntraciso,i,k)=', &
556	! & xt(1:ntraciso,i,k)
557	#endif
558	#endif
559	do ixt=1,ntraciso
560	xt(ixt,i,k)=max(xt(ixt,i,k),0.0)
561	enddo !do ixt=1,niso
562	#ifdef ISOTRAC
563	#ifdef ISOVERIF
564	call iso_verif_traceur_pbidouille(xt(1,i,k),'ilp 407')
565	#endif
566	#endif
567	do ixt=1,ntraciso
568	zxt(ixt,i)=xt(ixt,i,k)
569	enddo !do ixt=1,niso
570
571	! cam verif
572	#ifdef ISOVERIF
573	!write(,) 'ilp tmp 562'
574	call iso_verif_positif(370.0-zt(i),'il pleut 425')
575	call iso_verif_positif(zt(i)-100.0,'il pleut 426')
576	if (iso_eau.gt.0) then
577	call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
578	& 'il pleut 289',errmax,errmaxrel)
579	endif !if (iso_eau.gt.0) then
580	if (iso_HDO.gt.0) then
581	if (zq(i).gt.ridicule) then
582	call iso_verif_aberrant(zxt(iso_HDO,i)/zq(i), &
583	& 'il pleut 337')
584	endif !if (zq(i).gt.ridicule) then
585	endif !if (iso_HDO.gt.0) then
586	#ifdef ISOTRAC
587	call iso_verif_traceur(zxt(1,i),'ilp 431')
588	#endif
589	#endif
590	! end cam verif
591	#endif
592	ENDDO
593	!
594	! Calculer la varition de temp. de l'air du a la chaleur sensible
595	! transporter par la pluie.
596	! Il resterait a rajouter cet effet de la chaleur sensible sur les
597	! flux de surface, du a la diff. de temp. entre le 1er niveau et la
598	! surface.
599	!
600	IF(k.LE.klevm1) THEN
601	DO i = 1, klon
602	!IM
603	zmair=(paprs(i,k)-paprs(i,k+1))/RG
604	zcpair=RCPD(1.0+RVTMP2zq(i))
605	zcpeau=RCPD*RVTMP2
606	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zrfl(i)dtime*zcpeau &
607	+ zmairzcpairzt(i) ) &
608	/ (zmairzcpair + zrfl(i)dtime*zcpeau)
609	! C WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
610	ENDDO
611	ENDIF
612	! ----------------------------------------------------------------
613	! P1> Debut evaporation de la precipitation
614	! ----------------------------------------------------------------
615
616	#ifdef ISO
617	DO i = 1, klon
618	zq_ancien(i)=zq(i)
619	zqev_diag(i)=0.0
620	zrfl_ancien(i)=zrfl(i)
621	do ixt=1,ntraciso
622	zxtrfl_ancien(ixt,i)=zxtrfl(ixt,i)
623	zxt_ancien(ixt,i)=zxt(ixt,i)
624	enddo
625	enddo ! DO i = 1, klon
626	#ifdef ISOVERIF
627	!write(,) 'ilp tmp 616'
628	! i=1
629	! write(,) 'ilp 604: i=',i
630	! write(,) 'zrfl_ancien(i)=',zrfl_ancien(i)
631	! write(,) 'zrfln(i)=',zrfln(i)
632	! write(,) 'zifl(i)=',zifl(i)
633	if (iso_eau.gt.0) then
634	DO i = 1, klon
635	call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
636	& 'il pleut 426a',errmax,errmaxrel)
637	call iso_verif_egalite_choix(zxtrfl_ancien(iso_eau,i), &
638	& zrfl_ancien(i), &
639	& 'il pleut 426b',errmax,errmaxrel)
640	enddo
641	endif
642	#ifdef ISOTRAC
643	do i=1,klon
644	call iso_verif_traceur(zxt_ancien(1,i),'ilp 532')
645	call iso_verif_traceur(zxtrfl_ancien(1,i), &
646	& 'ilp 533')
647	enddo !do i=1,klon
648	#endif
649	#endif
650	#endif
651
652	IF (evap_prec) THEN
653	DO i = 1, klon
654	!AJ<
655	!! IF (zrfl(i) .GT.0.) THEN
656	IF (zrfl(i)+zifl(i).GT.0.) THEN
657	!>AJ
658	! Calcul du qsat
659	IF (thermcep) THEN
660	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
661	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
662	zqs(i)=MIN(0.5,zqs(i))
663	zcor=1./(1.-RETV*zqs(i))
664	zqs(i)=zqs(i)*zcor
665	ELSE
666	IF (zt(i) .LT. t_coup) THEN
667	zqs(i) = qsats(zt(i)) / pplay(i,k)
668	ELSE
669	zqs(i) = qsatl(zt(i)) / pplay(i,k)
670	ENDIF
671	ENDIF
672	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
673	ENDDO
674	!AJ<
675	IF (.NOT. ice_thermo) THEN
676	DO i = 1, klon
677	!AJ<
678	!! IF (zrfl(i) .GT.0.) THEN
679	IF (zrfl(i)+zifl(i).GT.0.) THEN
680	!>AJ
681	! Evap max pour ne pas saturer la fraction sous le nuage
682	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
683	! Calcul de l'evaporation du flux de precip herite
684	! d'au-dessus
685	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
686	* (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
687	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
688	* RG*dtime/(paprs(i,k)-paprs(i,k+1))
689	! Seuil pour ne pas saturer la fraction sous le nuage
690	zqev = MIN (zqev, zqevt)
691	! Nouveau flux de precip
692	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
693	/RG/dtime
694
695	#ifdef ISO
696	! write(,) 'ilp tmp 686'
697	zqev_diag(i)=zqev
698	if (essai_convergence) then
699	else
700	zrfln(i)=max(0.0,zrfln(i))
701	endif
702	#endif
703
704
705	! Aucun flux liquide pour T < t_coup
706	IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
707	! Nouvelle vapeur
708	zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
709	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
710	! Nouvelle temperature (chaleur latente)
711	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
712	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
713	* RLVTT/RCPD/(1.0+RVTMP2*zq(i))
714	zrfl(i) = zrfln(i)
715	zifl(i) = 0.
716	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
717	ENDDO
718	#ifdef ISO
719	! write(,) 'ilp tmp 709'
720	! i=1
721	! write(,) 'ilp 701 tmp: i=',i
722	! write(,) 'zrfl(i)=',zrfl(i)
723	! write(,) 'zrfln(i)=',zrfln(i)
724	! write(,) 'zifl(i)=',zifl(i)
725	! write(,) 'zqev_diag(i)=',zqev_diag(i)
726
727	! * traitement de l'évaporation
728	call iso_revap_fisrtilp(klon,klev,k, &
729	& zrfl_ancien,zrfl,zrfln,zt,zxt_ancien, &
730	& zxtrfl,zxtrfl_ancien,zxtrfln,zxt, &
731	& paprs,dtime, &
732	& zqs,zq_ancien,zqev_diag,zq)
733	! write(,) 'ilp tmp 723'
734	#endif
735
736	!
737	ELSE ! (.NOT. ice_thermo)
738
739	#ifdef ISO
740	CALL abort_physic('ilp 664', 'isotopes pas encore prevus dans ice_thermo', 1)
741	#endif
742	!
743	DO i = 1, klon
744	!AJ<
745	!! IF (zrfl(i) .GT.0.) THEN
746	IF (zrfl(i)+zifl(i).GT.0.) THEN
747	!>AJ
748	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
749	! Modification de l'evaporation avec la glace
750	! Differentiation entre precipitation liquide et solide
751	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
752
753	! Evap max pour ne pas saturer la fraction sous le nuage
754	zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
755	! zqev0 = MAX (0.0, zqs(i)-zq(i) )
756
757	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
758	! On differencie qsat pour l'eau et la glace
759	! Si zdelta=1. --> glace
760	! Si zdelta=0. --> eau liquide
761	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
762
763	! Calcul du qsat par rapport a l'eau liquide
764	qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
765	qsl= MIN(0.5,qsl)
766	zcor= 1./(1.-RETV*qsl)
767	qsl= qsl*zcor
768
769	! Calcul de l'evaporation du flux de precip herite
770	! d'au-dessus
771	! Formulation en racine du flux de precip
772	! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
773	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
774	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
775	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
776	RGdtime/(paprs(i,k)-paprs(i,k+1))
777
778
779	! Calcul du qsat par rapport a la glace
780	qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
781	qsi= MIN(0.5,qsi)
782	zcor= 1./(1.-RETV*qsi)
783	qsi= qsi*zcor
784
785	! Calcul de la sublimation du flux de precip solide herite
786	! d'au-dessus
787	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
788	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
789	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
790	RGdtime/(paprs(i,k)-paprs(i,k+1))
791
792
793	!JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
794	! Verification sur l'evaporation
795	! On s'assure qu'on ne sature pas
796	! la fraction sous le nuage sinon on
797	! repartit zqev0 en gardant la proportion
798	! liquide / glace
799	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
800
801	IF (zqevt+zqevti.GT.zqev0) THEN
802	zqev=zqev0*zqevt/(zqevt+zqevti)
803	zqevi=zqev0*zqevti/(zqevt+zqevti)
804
805	ELSE
806	IF (zqevt+zqevti.GT.0.) THEN
807	zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
808	zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
809	ELSE
810	zqev=0.
811	zqevi=0.
812	ENDIF
813	ENDIF
814	! Nouveaux flux de precip liquide et solide
815	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
816	/RG/dtime)
817	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
818	/RG/dtime)
819
820	! Mise a jour de la vapeur, temperature et flux de precip
821	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
822	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
823	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
824	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
825	* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
826	+ (zifln(i)-zifl(i)) &
827	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
828	* RLSTT/RCPD/(1.0+RVTMP2*zq(i))
829
830	zrfl(i) = zrfln(i)
831	zifl(i) = zifln(i)
832
833	!CR ATTENTION: deplacement de la fonte de la glace
834	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
835	!!! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 !!!!!!!!! jyg
836	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
837	zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! jyg
838	zmelt = MIN(MAX(zmelt,0.),1.)
839	! Fusion de la glace
840	zrfl(i)=zrfl(i)+zmelt*zifl(i)
841	zifl(i)=zifl(i)*(1.-zmelt)
842	! print*,zt(i),'octavio1'
843	! Chaleur latente de fusion
844	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
845	RLMLT/RCPD/(1.0+RVTMP2zq(i))
846	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
847	!fin CR
848
849
850
851	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
852	ENDDO
853
854	ENDIF ! (.NOT. ice_thermo)
855
856	! ----------------------------------------------------------------
857	! Fin evaporation de la precipitation
858	! ----------------------------------------------------------------
859	ENDIF ! (evap_prec)
860	!
861	! Calculer Qs et L/Cp*dQs/dT:
862	#ifdef ISOVERIF
863	! write(,) 'ilp tmp 849'
864	do i=1,klon
865	call iso_verif_positif(370.0-zt(i),'il pleut 938')
866	call iso_verif_positif(zt(i)-100.0,'il pleut 939')
867	do ixt=1,ntraciso
868	call iso_verif_positif_choix(zxt(ixt,i),0.0,'ilp 614')
869	enddo
870	if (zq(i).gt.ridicule) then
871	if (iso_HDO.gt.0) then
872	call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
873	& 'il pleut 856')
874	if (iso_O18.gt.0) then
875	if (iso_verif_O18_aberrant_nostop(zxt(iso_HDO,i)/zq(i), &
876	& zxt(iso_O18,i)/zq(i),'il pleut 859').eq.1) then
877	write(,) 'i,k,zq=',i,k,zq(i)
878	write(,) 'zqev_diag(i)=',zqev_diag(i)
879	write(,) 'zrfl_ancien(i)=',zrfl_ancien(i)
880	write(,) 'zrfln(i)=',zrfln(i)
881	write(,) 'zq_ancien(i)=',zq_ancien(i)
882	write(,) 'deltaD(zrfl_ancien)=',deltaD(zxtrfl_ancien(iso_HDO,i))
883	write(,) 'deltaO(zrfl_ancien)=',deltaO(zxtrfl_ancien(iso_O18,i))
884	write(,) 'deltaD(zq_ancien)=',deltaD(zxt_ancien(iso_HDO,i))
885	write(,) 'deltaO(zq_ancien)=',deltaO(zxt_ancien(iso_O18,i))
886	write(,) 'zt=',zt(i)
887	stop
888	endif !if (iso_verif_O18_aberrant_nostop
889	endif !if (iso_O18.gt.0) then
890	endif !if (iso_HDO.gt.0) then
891	endif !if (zq(i).gt.ridicule) then
892	enddo !do i=1,klon
893	#ifdef ISOTRAC
894	! write(,) 'ilp tmp 633'
895	do i=1,klon
896	call iso_verif_traceur(zxtrfl(1,i),'il pleut 632')
897	enddo
898	! write(,) 'ilp tmp 637,thermcep'
899	#endif
900	#endif
901	#ifdef ISOVERIF
902	! write(,) 'ilp tmp 888'
903	do i=1,klon
904	do ixt=1,ntraciso
905	call iso_verif_noNaN(zxtrfl(ixt,i),'ilp 624')
906	enddo !do ixt=1,ntraciso
907	enddo ! do i=1,klon
908	! write(,) 'ilp tmp 638: k,zoliq,zxtoliq(1,363)=',k,zoliq(363),zxtoliq(1,363)
909	#endif
910	!
911	IF (thermcep) THEN
912	DO i = 1, klon
913	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
914	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
915	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
916	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
917	zqs(i) = MIN(0.5,zqs(i))
918	zcor = 1./(1.-RETV*zqs(i))
919	zqs(i) = zqs(i)*zcor
920	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
921	ENDDO
922	ELSE
923	DO i = 1, klon
924	IF (zt(i).LT.t_coup) THEN
925	zqs(i) = qsats(zt(i))/pplay(i,k)
926	zdqs(i) = dqsats(zt(i),zqs(i))
927	ELSE
928	zqs(i) = qsatl(zt(i))/pplay(i,k)
929	zdqs(i) = dqsatl(zt(i),zqs(i))
930	ENDIF
931	ENDDO
932	ENDIF
933	!
934	! Determiner la condensation partielle et calculer la quantite
935	! de l'eau condensee:
936	!
937	!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
938	! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
939	! write(,) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
940	! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
941	! stop
942	! endif
943
944	! ----------------------------------------------------------------
945	! P2> Formation du nuage
946	! ----------------------------------------------------------------
947	IF (cpartiel) THEN
948
949	! print*,'Dans partiel k=',k
950	!
951	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
952	! nuageuse a partir des PDF de Sandrine Bony.
953	! rneb : fraction nuageuse
954	! zqn : eau totale dans le nuage
955	! zcond : eau condensee moyenne dans la maille.
956	! on prend en compte le réchauffement qui diminue la partie
957	! condensee
958	!
959	! Version avec les raqts
960
961	#ifdef CAMMODIF
962	if (ok_sursat_marine.eq.1) then
963	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
964	! Marine: c'est ici qu'intervient la sursaturation
965	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
966	! Calcul de zqs_star, seuil critique de sursaturation
967
968
969	DO i = 1, klon
970
971	IF (zt(i) .LT. t_coup) THEN
972	CALL thresh_supersat(zt(i),RH_crit)
973	zqs_star(i) = zqs(i)*RH_crit
974	zpdf_beta(i) = RH_crit
975	ELSE
976	RH_crit = 1.
977	zqs_star(i) = zqs(i)
978	zpdf_beta(i) = RH_crit
979	ENDIF
980
981	if (zqs(i) .gt. zqs_star(i)) then
982	write(,) 'ATTENTION zqs > zqs_star!!'
983	write(,) 'zqs zqs_star RH_crit=', &
984	& zqs(i), zqs_star(i), RH_crit
985	STOP
986	endif !if (zqs(i) .gt. zqs_star(i)) then
987
988	ENDDO !DO i = 1, klon
989	endif !if (sursat_marine.eq.1) then
990	#endif
991	!#ifdef CAMMODIF
992
993	#ifdef ISO
994	#ifdef ISOVERIF
995	do i=1,klon
996	if (iso_eau.gt.0) then
997	if (iso_verif_egalite_choix_nostop(zxt(iso_eau,i),zq(i), &
998	& 'il pleut 608',errmax,errmaxrel).eq.1) then
999	write(,) 'i=',i
1000	stop
1001	endif
1002	endif !if (iso_eau.gt.0) then
1003	#ifdef ISOTRAC
1004	call iso_verif_traceur(zxt(1,i),'il pleut 1106')
1005	#endif
1006	enddo !do i=1,klon
1007	#endif
1008	#endif
1009
1010	if (iflag_pdf.eq.0) then
1011
1012	! version creneau de (Li, 1998)
1013	do i=1,klon
1014	zdelq = min(ratqs(i,k),0.99) * zq(i)
1015	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
1016	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
1017	enddo
1018
1019	#ifdef CAMMODIF
1020	if (ok_sursat_marine.eq.1) then
1021	DO i=1,klon
1022	IF (rneb(i,k) .LE. 0.0) THEN
1023	rhcl(i,k)=zq(i)/zqs(i)
1024	ELSE IF (rneb(i,k) .GE. 1.0) THEN
1025	rhcl(i,k)=1.0
1026	ELSE
1027	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1028	ENDIF
1029	ENDDO
1030	endif !if (ok_sursat_marine.eq.1) then
1031	#endif
1032	!#ifdef CAMMODIF
1033	#ifdef ISO
1034	! Hypothèse simplificatrice mais peut-être fausse:
1035	! eau nuageuse a même composition que eau de la maille
1036	! verif
1037	#ifdef ISOVERIF
1038	do i=1,klon
1039	if (iso_eau.gt.0) then
1040	call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
1041	& 'il pleut 618',errmax,errmaxrel)
1042	endif !if (iso_eau.gt.0) then
1043	enddo !do i=1,klon
1044	#endif
1045	! end verif
1046	do i=1,klon
1047	do ixt=1,ntraciso
1048	if (zq(i).gt.0) then
1049	zxtn(ixt,i)=zqn(i)*(zxt(ixt,i)/zq(i))
1050	zxtn(ixt,i)=max(0.0,zxtn(ixt,i))
1051	else
1052	zxtn(ixt,i)=0.
1053	endif
1054	enddo !do ixt=1,niso
1055	enddo !do i=1,klon
1056
1057	! verif
1058	#ifdef ISOVERIF
1059	do i=1,klon
1060	if (iso_eau.gt.0) then
1061	call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
1062	& 'il pleut 625',errmax,errmaxrel)
1063	endif !if (iso_eau.gt.0) then
1064	enddo !do i=1,klon
1065	#endif
1066	! end verif
1067
1068	#endif
1069	! write(,) 'il pleut 641: fin des calculs de zqn et zxtn'
1070	! write(,) 'iflag_pdf=',iflag_pdf
1071	#ifdef ISOTRAC
1072	if (option_tmin.eq.1) then
1073	!write(,) 'option_tmin pas encore prevue ici.'
1074	!write(,) 'seulement prevue dans iflag_pdf=1'
1075	!write(,) 'a programmer'
1076	CALL abort_physic('ilp 970', 'option_tmin pas encore prevue ici', 1)
1077	endif !if (option_tmin.eq.1) then
1078	#endif
1079
1080	else !if (iflag_pdf.eq.0) then
1081
1082	#ifdef ISO
1083	! verif
1084	#ifdef ISOVERIF
1085	do i=1,klon
1086	if (iso_eau.gt.0) then
1087	call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
1088	& 'il pleut 655',errmax,errmaxrel)
1089	endif !if (iso_eau.gt.0) then
1090	#ifdef ISOTRAC
1091	call iso_verif_traceur(zxt(1,i),'il pleut 1164')
1092	#endif
1093	enddo !do i=1,klon
1094	! write(,) 'ilp tmp 753'
1095	#endif
1096	! end verif
1097	#endif
1098	!
1099	! Version avec les nouvelles PDFs.
1100	do i=1,klon
1101	if(zq(i).lt.1.e-15) then
1102	ncoreczq=ncoreczq+1
1103	#ifdef ISO
1104	zq_ancien(i)=zq(i)
1105	#endif
1106	zq(i)=1.e-15
1107	#ifdef ISO
1108	! on conserve le même rapport pour les isotopes
1109	if (zq_ancien(i).gt.0.0) then
1110	do ixt=1,ntraciso
1111	zxt(ixt,i)=zxt(ixt,i)/zq_ancien(i)*zq(i)
1112	zxt(ixt,i)=max(0.0,zxt(ixt,i))
1113	enddo
1114	else !if (zq_ancien(i).gt.0.0) then
1115	do ixt=1,ntraciso
1116	zxt(ixt,i)=0.0
1117	enddo
1118	endif !if (zq_ancien(i).gt.0.0) then
1119
1120	! verif
1121	#ifdef ISOVERIF
1122	if (iso_eau.gt.0) then
1123	call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
1124	& 'il pleut 654',errmax,errmaxrel)
1125	endif !if (iso_eau.gt.0) then
1126	#ifdef ISOTRAC
1127	call iso_verif_traceur(zxt(1,i),'il pleut 1207')
1128	#endif
1129	#endif
1130	! end verif
1131	#endif
1132
1133	endif
1134	enddo
1135
1136	if (iflag_cld_th>=5) then
1137
1138	call cloudth(klon,klev,k,ztv, &
1139	zq,zqta,fraca, &
1140	qcloud,ctot,zpspsk,paprs,ztla,zthl, &
1141	ratqs,zqs,t)
1142
1143	do i=1,klon
1144	rneb(i,k)=ctot(i,k)
1145	zqn(i)=qcloud(i)
1146
1147	#ifdef CAMMODIF
1148	if (ok_sursat_marine.eq.1) then
1149	! on prend le rhcl du cas creneau
1150
1151	IF (rneb(i,k) .LE. 0.0) THEN
1152	rhcl(i,k)=zq(i)/zqs(i)
1153	ELSE IF (rneb(i,k) .GE. 1.0) THEN
1154	rhcl(i,k)=1.0
1155	ELSE
1156	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1157	ENDIF
1158	endif !if (ok_sursat_marine.eq.1) then
1159	! if (ok_sursat_marine.eq.1) then
1160	#endif
1161	!#ifdef CAMMODIF
1162
1163	enddo
1164
1165	endif !if (iflag_cldcon>=5) then
1166
1167	if (iflag_cld_th <= 4) then
1168	lognormale = .true.
1169	elseif (iflag_cld_th >= 6) then
1170	! lognormale en l'absence des thermiques
1171	lognormale = fraca(:,k) < 1e-10
1172	else
1173	! Dans le cas iflag_cld_th=5, on prend systématiquement la
1174	! bi-gaussienne
1175	lognormale = .false.
1176	end if
1177
1178	!CR: variation de qsat avec T en presence de glace ou non
1179	!initialisations
1180	do i=1,klon
1181	DT(i) = 0.
1182	n_i(i)=0
1183	Tbef(i)=zt(i)
1184	qlbef(i)=0.
1185	enddo
1186
1187
1188	!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
1189	if (iflag_fisrtilp_qsat.ge.0) then
1190	! Iteration pour condensation avec variation de qsat(T)
1191	! -----------------------------------------------------
1192	do iter=1,iflag_fisrtilp_qsat+1
1193
1194	do i=1,klon
1195	! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
1196	convergence(i)=abs(DT(i)).gt.DDT0
1197	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
1198	Tbef(i)=Tbef(i)+DT(i)
1199	if (.not.ice_thermo) then
1200	zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
1201	else
1202	if (iflag_t_glace.eq.0) then
1203	zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
1204	else if (iflag_t_glace.ge.1) then
1205	zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
1206	endif
1207	endif
1208	! Calcul des PDF lognormales
1209	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
1210	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
1211	zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
1212	zqs(i) = MIN(0.5,zqs(i))
1213	zcor = 1./(1.-RETV*zqs(i))
1214	zqs(i) = zqs(i)*zcor
1215	zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
1216	zpdf_sig(i)=ratqs(i,k)*zq(i)
1217	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
1218
1219	#ifdef CAMMODIF
1220	if (ok_sursat_marine.eq.1) then
1221	ys(i) = - (log(zqs(i)/zq(i)))/zpdf_k(i) - zpdf_k(i)/2.
1222	yss(i) = ys(i) - log(zpdf_beta(i))/zpdf_k(i)
1223
1224	if (yss(i) .lt. ys(i)) then
1225	write(,) 'attention yss < ys !!'
1226	STOP
1227	endif
1228	endif !if (ok_sursat_marine.eq.1) then
1229	#endif
1230	#ifdef CAMMODIF
1231	if ((ok_sursat_marine.eq.0).or. &
1232	& (zqs(i) .eq. zqs_star(i))) then
1233	#endif
1234
1235	zpdf_delta(i)=log(zq(i)/zqs(i))
1236	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
1237	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
1238	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
1239	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
1240	zpdf_e1(i)=1.-erf(zpdf_e1(i))
1241	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
1242	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
1243	zpdf_e2(i)=1.-erf(zpdf_e2(i))
1244
1245	if (zpdf_e1(i).lt.1.e-10) then
1246	rneb(i,k)=0.
1247	zqn(i)=zqs(i)
1248	else
1249	rneb(i,k)=0.5*zpdf_e1(i)
1250	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
1251	endif
1252
1253	#ifdef CAMMODIF
1254	else ! sursaturation
1255	! on est dans ce cas soit si ok_sursat_marine.eq.1 et si on a
1256	! sursat
1257
1258	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1259	! Cas ou il y a de la sursaturation
1260	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1261
1262	pi = ACOS(-1.)
1263
1264	! Calcul de rneb orig
1265	!cccccccccccccccccccc
1266
1267	zpdf_delta_o(i) = log(zq(i)/zqs(i))
1268	zpdf_a_o(i) = zpdf_delta_o(i)/(zpdf_k(i)*sqrt(2.))
1269	zpdf_b_o(i) = zpdf_k(i)/(2.*sqrt(2.))
1270	zpdf_e1_o(i) = zpdf_a_o(i)-zpdf_b_o(i)
1271	zpdf_e1_o(i) = sign(min(abs(zpdf_e1_o(i)),5.),zpdf_e1_o(i))
1272	zpdf_e1_o(i)=1.-erf(zpdf_e1_o(i))
1273
1274	if (zpdf_e1_o(i) .lt. 1.e-10) then
1275	rneb_o(i,k) = 0.
1276	else
1277	rneb_o(i,k) = 0.5*zpdf_e1_o(i)
1278	endif
1279
1280	! Calcul de rneb
1281	!ccccccccccccccc
1282
1283	if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1284
1285	zpdf_e1(i) = 1.-erf(yss(i)/sqrt(2.))
1286
1287	if (zpdf_e1(i) .lt. 1.e-10) then
1288	rneb(i,k) = exp(-yss(i)*2/2.)(yss(i)-ys(i))/(2.sqrt(2.pi))
1289	else
1290	rneb(i,k) = 0.5*zpdf_e1(i) + &
1291	& exp(-yss(i)*2/2.)(yss(i)-ys(i))/(2.sqrt(2.pi))
1292	endif
1293
1294	if (rneb(i,k) .lt. 0.) then
1295	write(,) 'ys, yss, yss-ys, 0.5*erf(yss(i)/sqrt(2.))',&
1296	& ys(i), yss(i), yss(i)-ys(i), 0.5*erf(yss(i)/sqrt(2.))
1297	endif
1298
1299
1300	else !if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1301
1302	! DL de rneb pour beta proche de 1
1303
1304	zpdf_e1(i) = 1.-erf(ys(i)/sqrt(2.))
1305
1306	if (zpdf_e1(i) .lt. 1.e-10) then
1307	rneb(i,k) = exp(-0.5ys(i)2)(zpdf_beta(i)-1.)/(2.zpdf_k(i)sqrt(2.*pi))- &
1308	& exp(-0.5ys(i)2)(zpdf_beta(i)-1)*2/(4.zpdf_k(i)sqrt(2.pi))+ &
1309	& exp(-0.5ys(i)2)(zpdf_beta(i)-1)*3 &
1310	& (1.+2zpdf_k(i)2-ys(i)2)/(12.zpdf_k(i)*3sqrt(2.*pi))
1311	else
1312	rneb(i,k) = 0.5*zpdf_e1(i) + &
1313	& exp(-0.5ys(i)2)(zpdf_beta(i)-1.)/(2.zpdf_k(i)sqrt(2.*pi))- &
1314	& exp(-0.5ys(i)2)(zpdf_beta(i)-1)*2/(4.zpdf_k(i)sqrt(2.pi))+ &
1315	& exp(-0.5ys(i)2)(zpdf_beta(i)-1)*3 &
1316	& (1.+2zpdf_k(i)2-ys(i)2)/(12.zpdf_k(i)*3sqrt(2.*pi))
1317	endif
1318
1319	if (rneb(i,k) .lt. 0.) then
1320	write(,) 'rneb< 0 ds le cas DL'
1321	endif
1322
1323
1324	endif !if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1325
1326	if ( rneb(i,k) .lt. 0. .or. rneb(i,k) .gt. 1.) then
1327	write(,) 'SURSAT: rneb(i,k)< 0 or rneb(i,k)>1', &
1328	& 'rneb(i,k),ys(i),zpdf_beta(i),zpdf_k(i),pi', &
1329	& rneb(i,k),ys(i),zpdf_beta(i),zpdf_k(i),pi
1330	endif
1331
1332
1333
1334	! Flag presence de sursaturation
1335	!ccccccccccccccccccccccccccccccc
1336
1337	if (abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i)) .and. &
1338	& (rneb_o(i,k) - rneb(i,k) .gt. 1.e-5 .or. &
1339	& (rneb_o(i,k) .gt. 0.98 .and. rneb(i,k) .gt. 0.98))) then
1340
1341	flag_ss(i,k) = 1.
1342
1343	else
1344
1345	flag_ss(i,k) = 0.
1346
1347	if (rneb_o(i,k) - rneb(i,k) .lt. -1.e-5) then
1348	write(,) 'rneb_o < rneb, rneb_o, rneb, yss, ys,&
1349	& 0.5zpdf_e1_o, 0.5zpdf_e1, exp',&
1350	& rneb_o(i,k),rneb(i,k),yss(i), ys(i), &
1351	& 0.5zpdf_e1_o(i),0.5zpdf_e1(i), &
1352	& exp(-yss(i)*2/2.)(yss(i)-ys(i))/(2.sqrt(2.pi))
1353	endif
1354
1355	endif
1356
1357	! Fraction de ciel clair sursature
1358	!ccccccccccccccccccccccccccccccccc
1359
1360	if (abs(zpdf_beta(i) -1.) .gt. epsilon(zpdf_beta(i)) .and. &
1361	& rneb_o(i,k) - rneb(i,k) .gt. 1.e-5) then
1362	clsk_ss(i,k) = rneb_o(i,k) - rneb(i,k)
1363	else
1364	clsk_ss(i,k) = 0.
1365	endif
1366
1367
1368	! Nuage dans un environnement sursature
1369	!ccccccccccccccccccccccccccccccccccccccc
1370
1371	if (rneb(i,k) .gt. 0 .and. &
1372	& abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1373	rneb_ss(i,k) = rneb(i,k)
1374	endif
1375
1376	! Autre detection de sursaturation
1377	!ccccccccccccccccccccccccccccccccc
1378
1379	if (rneb_o(i,k) .gt. 0 .and. &
1380	& abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1381	flagb_ss(i,k) = 1.
1382	else
1383	flagb_ss(i,k) = 0.
1384	endif
1385
1386	! Calcul de zqn
1387	!cccccccccccccc
1388
1389	! Calcul de zqn
1390	!cccccccccccccc
1391
1392	if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1393
1394	zpdf_e1(i) = 0.5zq(i)(1.-erf((yss(i)+zpdf_k(i))/sqrt(2.)))
1395
1396	zqn(i) = zpdf_e1(i) + &
1397	& exp(-zpdf_k(i)2/2.-yss(i)2/2.-zpdf_k(i)(ys(i)+yss(i))) &
1398	& zq(i)* &
1399	& (exp(zpdf_k(i)*yss(i))+ &
1400	& exp(zpdf_k(i)ys(i))(zpdf_k(i)ys(i)-zpdf_k(i)yss(i)-1))/ &
1401	& (zpdf_k(i)*2sqrt(2.pi)(yss(i)-ys(i)))
1402
1403	else !if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1404
1405	! DL de zqn pour beta proche de 1
1406
1407	zpdf_e1(i) = 0.5zq(i)(1.-erf((zpdf_k(i)+ys(i))/sqrt(2.)))
1408
1409	zqn(i) = zpdf_e1(i)+ &
1410	& exp(-0.5(zpdf_k(i)+ys(i))2)zq(i)*(zpdf_beta(i)-1.)/ &
1411	& (2.zpdf_k(i)sqrt(2.*pi))- &
1412	& exp(-0.5(zpdf_k(i)+ys(i))2)zq(i)(zpdf_beta(i)-1.)*2/ &
1413	& (12.zpdf_k(i)sqrt(2.*pi))+ &
1414	& exp(-0.5(zpdf_k(i)+ys(i))2)zq(i)(zpdf_beta(i)-1.)3 &
1415	& (2.+zpdf_k(i)*2-2.ys(i)*2)/(24.zpdf_k(i)*3sqrt(2.*pi))
1416
1417	endif !if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1418
1419	if (rneb(i,k) .gt. 0. .and. rneb(i,k) .lt. 1.) then
1420	zqn(i) = zqn(i)/rneb(i,k) ! a verifier
1421	endif
1422
1423	! Calcul de rhcl
1424	!ccccccccccccccc
1425
1426	! if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1427
1428	! rhcl_p0 = &
1429	! & erf((zpdf_k(i)+yss(i))/sqrt(2.))-erf(zpdf_k(i)/sqrt(2.))
1430
1431	! rhcl_p1 = &
1432	! & exp(-0.5zpdf_k(i)2-0.5yss(i)*2-zpdf_k(i)(ys(i)+yss(i)))
1433
1434	! rhcl_p2 = exp(zpdf_k(i)yss(i))+exp(zpdf_k(i)ys(i))* &
1435	! & (-1.+zpdf_k(i)ys(i)-zpdf_k(i)yss(i))
1436
1437	! rhcl_p = rhcl_p1rhcl_p2sqrt(2./pi)/ &
1438	! & ((ys(i)-yss(i))zpdf_k(i)*2)
1439
1440	! rhcl(i,k) = 0.5zq(i)(rhcl_p0+rhcl_p)
1441
1442	! rhcl(i,k) = rhcl(i,k)/zqs(i)
1443
1444	! else !if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then!
1445
1446	! DL de rhcl pour beta proche de 1
1447
1448	! rhcl(i,k) = &
1449	! & 0.5zq(i) &
1450	! & (erf((zpdf_k(i)+ys(i))/sqrt(2.))-erf(zpdf_k(i)/sqrt(2.)))- &
1451	! & exp(-0.5(zpdf_k(i)+ys(i))2)zq(i)*(zpdf_beta(i)-1.)/ &
1452	! & (2.zpdf_k(i)sqrt(2.*pi))+ &
1453	! & exp(-0.5(zpdf_k(i)+ys(i))2)zq(i)(zpdf_beta(i)-1.)*2/ &
1454	! & (12.zpdf_k(i)sqrt(2.*pi))- &
1455	! & exp(-0.5(zpdf_k(i)+ys(i))2)zq(i)(zpdf_beta(i)-1.)3 &
1456	! & (2.+zpdf_k(i)*2-2.ys(i)**2)/ &
1457	! & (24.zpdf_k(i)3sqrt(2.*pi))
1458
1459	! rhcl(i,k) = rhcl(i,k)/zqs(i)
1460
1461	! endif !if ( abs(zpdf_beta(i)-1.) .gt. epsilon(zpdf_beta(i))) then
1462
1463
1464	! IF (rneb(i,k) .LE. 0.0) THEN
1465	! rhcl(i,k)=zq(i)/zqs(i)
1466	! ELSE IF (rneb(i,k) .GE. 1.0) THEN
1467	! rhcl(i,k)=1.0
1468	! ELSE
1469	! rhcl(i,k)=rhcl(i,k)/(1.-rneb(i,k))
1470	! ENDIF
1471	endif ! sursaturation
1472	#endif
1473	!#ifdef CAMMODIF
1474
1475
1476	endif !convergence
1477	enddo ! boucle en i
1478
1479	if (.not. ice_thermo) then
1480
1481	do i=1,klon
1482	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
1483
1484	qlbef(i)=max(0.,zqn(i)-zqs(i))
1485	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
1486	denom=1.+rneb(i,k)*zdqs(i)
1487	DT(i)=num/denom
1488	n_i(i)=n_i(i)+1
1489	endif
1490	enddo
1491
1492	else
1493	! Iteration pour convergence avec qsat(T)
1494	if (iflag_t_glace.ge.1) then
1495	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1496	endif
1497
1498	do i=1,klon
1499	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
1500
1501	if (iflag_t_glace.eq.0) then
1502	zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1503	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1504	zfice(i) = zfice(i)**exposant_glace_old
1505	dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
1506	endif
1507
1508	if (iflag_t_glace.ge.1) then
1509	dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
1510	endif
1511
1512	if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
1513	dzfice(i)=0.
1514	endif
1515
1516	if (zfice(i).lt.1) then
1517	cste=RLVTT
1518	else
1519	cste=RLSTT
1520	endif
1521
1522	qlbef(i)=max(0.,zqn(i)-zqs(i))
1523	num=-Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
1524	denom=1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
1525	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))rneb(i,k)qlbef(i)dzfice(i)
1526	DT(i)=num/denom
1527	n_i(i)=n_i(i)+1
1528
1529	endif ! fin convergence
1530	enddo ! fin boucle i
1531
1532	endif !ice_thermo
1533
1534	! endif
1535	! enddo
1536
1537
1538	enddo ! iter=1,iflag_fisrtilp_qsat+1
1539	! Fin d'iteration pour condensation avec variation de qsat(T)
1540	! -----------------------------------------------------------
1541	endif
1542
1543
1544	endif ! iflag_pdf
1545
1546
1547	! if (iflag_fisrtilp_qsat.eq.-1) then
1548	!------------------------------------------
1549	!CR: ATTENTION option fausse mais a existe:
1550	! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
1551	! activer les lignes suivantes:
1552	IF (1.eq.0) THEN
1553	DO i=1,klon
1554	IF (rneb(i,k) .LE. 0.0) THEN
1555	zqn(i) = 0.0
1556	rneb(i,k) = 0.0
1557	zcond(i) = 0.0
1558	rhcl(i,k)=zq(i)/zqs(i)
1559	ELSE IF (rneb(i,k) .GE. 1.0) THEN
1560	zqn(i) = zq(i)
1561	rneb(i,k) = 1.0
1562	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
1563	rhcl(i,k)=1.0
1564	ELSE
1565	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
1566	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1567	ENDIF
1568	ENDDO
1569	ENDIF
1570	!------------------------------------------
1571
1572	! ELSE
1573
1574	! Calcul de l'eau in-cloud (zqn),
1575	! moyenne dans la maille (zcond),
1576	! fraction nuageuse (rneb) et
1577	! humidite relative ciel-clair (rhcl)
1578	DO i=1,klon
1579	IF (rneb(i,k) .LE. 0.0) THEN
1580	zqn(i) = 0.0
1581	rneb(i,k) = 0.0
1582	zcond(i) = 0.0
1583	rhcl(i,k)=zq(i)/zqs(i)
1584	ELSE IF (rneb(i,k) .GE. 1.0) THEN
1585	zqn(i) = zq(i)
1586	rneb(i,k) = 1.0
1587	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1588	rhcl(i,k)=1.0
1589	ELSE
1590	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1591	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1592	ENDIF
1593	ENDDO
1594
1595
1596	! ENDIF
1597
1598	ELSE ! de IF (cpartiel)
1599	! Cas "tout ou rien"
1600	DO i = 1, klon
1601	IF (zq(i).GT.zqs(i)) THEN
1602	rneb(i,k) = 1.0
1603	ELSE
1604	rneb(i,k) = 0.0
1605	ENDIF
1606	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
1607	ENDDO
1608	ENDIF
1609
1610	#ifdef ISO
1611	! on calcule la compo iso de l'eau dans le nuage et en ciel
1612	! clair: zxtn/zqn et zxtcs/zqcs
1613
1614	! Hypothèsesimplificatrice mais peut-être fausse:
1615	! eau nuageuse a même composition que eau de la maille
1616
1617	! ce calcul n'eat pas sensé dépendre du type de PDF utilisé
1618	do i=1,klon
1619	#ifdef ISOTRAC
1620	if (option_tmin.eq.1) then
1621	! c'est mieux si zqn>=0.0
1622	if (essai_convergence) then
1623	else
1624	zqn(i)=max(0.0,zqn(i))
1625	endif
1626
1627	if (rneb(i,k).lt.1.0) then
1628	zqcs(i)=(zq(i)-zqn(i)*rneb(i,k))/(1.0-rneb(i,k))
1629	#ifdef ISOVERIF
1630	! write(,) 'rneb,zq,zqn,zqcs=',rneb(i,k),zq(i), &
1631	! & zqn(i),zqcs(i)
1632	call iso_verif_positif(zqcs(i),'ipleut 853')
1633	#endif
1634	else
1635	zqcs(i)=0.0
1636	endif
1637	endif ! if (option_tmin.eq.1) then
1638	#endif
1639	#ifdef ISOVERIF
1640	if (iso_eau.gt.0) then
1641	call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
1642	& 'il pleut 748',errmax,errmaxrel)
1643	endif !if (iso_eau.gt.0) then
1644	do ixt=1,ntraciso
1645	call iso_verif_positif_choix(zxt(ixt,i),0.0,'ilp 857')
1646	enddo
1647	#ifdef ISOTRAC
1648	call iso_verif_traceur_pbidouille(zxt(1,i), &
1649	& 'il pleut 874')
1650	#endif
1651	#endif
1652	if ((bidouille_anti_divergence).and. &
1653	& (iso_eau.gt.0)) then
1654	zxt(iso_eau,i)= max(zq(i),0.0)
1655	! if (zxt(iso_eau,i).gt.0.0) then
1656	! do ixt=1,niso
1657	! zxt(ixt,i)=zxt(ixt,i)*zq(i)/zxt(iso_eau,i)
1658	! enddo
1659	! else !if (zxt(iso_eau,i).gt.0.0) then
1660	! do ixt=1,niso
1661	! zxt(ixt,i)=0.0
1662	! enddo
1663	! endif !if (zxt(iso_eau,i).gt.0.0) then
1664	endif !if ((bidouille_anti_divergence).and.
1665
1666	if (zq(i).gt.1e-15) then
1667	! 2 fev 2010: on modifie seuil de zq
1668	! avant, c'était 0. Vérifier que ça n'a pas de
1669	! conséquences trop graves pour les isotopes.
1670	do ixt=1,ntraciso
1671	zxtn(ixt,i)=zqn(i)*(zxt(ixt,i)/zq(i))
1672	enddo !do ixt=1,niso
1673	#ifdef ISOTRAC
1674	if (option_tmin.eq.1) then
1675	do ixt=1,ntraciso
1676	zxtcs(ixt,i)=(zqcs(i))*(zxt(ixt,i)/zq(i))
1677	enddo !do ixt=1,ntraciso
1678	endif ! if (option_tmin.eq.1) then
1679	#endif
1680	else !if (zq(i).gt.0.0) then
1681	! problème: zq~0 et sans compo, mais zqn >0
1682	! on fait ce qu'on peut dans ce cas pathologique
1683	do ixt=1,ntraciso
1684	zxtn(ixt,i)=0.0
1685	enddo !do ixt=1,niso
1686	#ifdef ISOTRAC
1687	if (option_tmin.eq.1) then
1688	do ixt=1,ntraciso
1689	zxtcs(ixt,i)=0.0
1690	enddo !do ixt=1,niso
1691	endif ! if (option_tmin.eq.1) then
1692	#endif
1693	if ((bidouille_anti_divergence) &
1694	& .and.(iso_eau.gt.0)) then
1695	! write(,) 'ilp tmp 848: warning: bidouille'
1696	zxtn(iso_eau,i)=zqn(i)
1697	#ifdef ISOTRAC
1698	zxtn(index_trac(izone_poubelle,iso_eau),i)=zqn(i)
1699	if (option_tmin.eq.1) then
1700	zxtcs(iso_eau,i)=zqcs(i)
1701	endif ! if (option_tmin.eq.1) then
1702	! write(,) 'zxtn(:,i)=',zxtn(:,i)
1703	#endif
1704	endif !if ((bidouille_anti_divergence)
1705	endif !if (zq(i).gt.0.0) then
1706
1707	#ifdef ISOVERIF
1708	do ixt=1,niso
1709	call iso_verif_positif(zxtn(ixt,i),'il pleut 1062')
1710	enddo !do ixt=1,niso
1711	if (iso_eau.gt.0) then
1712	call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
1713	& 'il pleut 932',errmax,errmaxrel)
1714	endif !if ((iso_eau.gt.0)) then
1715	#ifdef ISOTRAC
1716	! write(,) 'zxt(:,i)=',zxt(:,i)
1717	! write(,) 'zq,zqn=',zq(i) ,zqn(i)
1718	call iso_verif_traceur_pbidouille(zxtn(1,i), &
1719	& 'il pleut 1277')
1720	if (option_tmin.eq.1) then
1721	call iso_verif_positif(zxtcs(ixt,i),'il pleut 895')
1722	do ixt=1,niso
1723	call iso_verif_positif(zxtcs(ixt,i),'il pleut 897')
1724	enddo !do ixt=1,niso
1725	if (iso_eau.gt.0) then
1726	call iso_verif_egalite_choix(zxtcs(iso_eau,i), &
1727	& zqcs(i),'il pleut 733',errmax,errmaxrel)
1728	endif !if ((iso_eau.gt.0)) then
1729	endif ! if (option_tmin.eq.1) then
1730	#endif
1731	#endif
1732	! write(,) 'il pleut 695: calcul de zxtn, i=',i
1733	! write(,) 'zxtn(4,i)=',zxtn(4,i)
1734	! write(,) 'zqn(i)=',zqn(i)
1735	do ixt=1,ntraciso
1736	zxtn(ixt,i)=max(0.0,zxtn(ixt,i))
1737	enddo !do ixt=1,niso
1738	#ifdef ISOTRAC
1739	if (option_tmin.eq.1) then
1740	do ixt=1,ntraciso
1741	zxtcs(ixt,i)=max(0.0,zxtcs(ixt,i))
1742	enddo
1743	endif ! if (option_tmin.eq.1) then
1744	#endif
1745
1746	#ifdef ISOVERIF
1747	if (iso_eau.gt.0) then
1748	call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
1749	& 'il pleut 746',errmax,errmaxrel)
1750	endif !if (iso_eau.gt.0) then
1751	if (iso_HDO.gt.0) then
1752	if (zqn(i).gt.ridicule) then
1753	call iso_verif_aberrant(zxtn(iso_HDO,i)/zqn(i), &
1754	& 'il pleut 977')
1755	endif !if (zqn(i).gt.ridicule) then
1756	endif !if (iso_HDO.gt.0) then
1757	#ifdef ISOTRAC
1758	call iso_verif_traceur_pbidouille &
1759	& (zxtn(1,i),'il pleut 1300')
1760	if (option_tmin.eq.1) then
1761	if (iso_eau.gt.0) then
1762	call iso_verif_egalite_choix(zxtcs(iso_eau,i), &
1763	& zqcs(i), &
1764	& 'il pleut 971',errmax,errmaxrel)
1765	endif !if (iso_eau.gt.0) then
1766	call iso_verif_traceur_pbidouille &
1767	& (zxtcs(1,i),'il pleut 1300')
1768
1769	endif ! if (option_tmin.eq.1) then
1770	#endif
1771	#endif
1772	enddo ! do i=1,klon
1773
1774	! write(,) 'ilp tmp 1311: zoliq,zxtoliq(1,363)=',zoliq(363),zxtoliq(1,363)
1775	#endif
1776
1777
1778	! ----------------------------------------------------------------
1779	! Fin de formation du nuage
1780	! ----------------------------------------------------------------
1781	!
1782	! Mise a jour vapeur d'eau
1783	DO i = 1, klon
1784	zq(i) = zq(i) - zcond(i)
1785	! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
1786	ENDDO
1787	!AJ<
1788	! Chaleur latente apres formation nuage
1789	! -------------------------------------
1790	IF (.NOT. ice_thermo) THEN
1791	if (iflag_fisrtilp_qsat.lt.1) then
1792	DO i = 1, klon
1793	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
1794	ENDDO
1795	else if (iflag_fisrtilp_qsat.gt.0) then
1796	DO i= 1, klon
1797	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1798	ENDDO
1799	endif
1800	ELSE
1801	if (iflag_t_glace.ge.1) then
1802	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1803	endif
1804	if (iflag_fisrtilp_qsat.lt.1) then
1805	DO i = 1, klon
1806	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1807	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1808	! t_glace_max, exposant_glace)
1809	if (iflag_t_glace.eq.0) then
1810	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1811	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1812	zfice(i) = zfice(i)**exposant_glace_old
1813	endif
1814	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
1815	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
1816	ENDDO
1817	else
1818	DO i=1, klon
1819	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1820	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1821	! t_glace_max, exposant_glace)
1822	if (iflag_t_glace.eq.0) then
1823	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1824	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1825	zfice(i) = zfice(i)**exposant_glace_old
1826	endif
1827	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
1828	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1829	ENDDO
1830	endif
1831	! print*,zt(i),zrfl(i),zifl(i),'temp1'
1832	ENDIF
1833
1834	! isotopes lors de la condensation
1835	#ifdef ISO
1836	#ifdef ISOVERIF
1837	! write(,) 'ilp tmp 1823'
1838	do i=1,klon
1839	! verif
1840	if (iso_eau.gt.0) then
1841	call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
1842	& 'il pleut 810',errmax,errmaxrel)
1843	endif !if (iso_eau.gt.0) then
1844	#ifdef ISOTRAC
1845	call iso_verif_traceur(zxtn(1,i),'il pleut 1408')
1846	! write(,) 'fisrtilp 1421: appel condiso_liq_ice_vectall'
1847	#endif
1848	enddo !do i=1,klon
1849	! write(,) 'il pleut 961: avant condensation'
1850	! end verif
1851	#endif
1852
1853	do i=1,klon
1854	! call calcul_zfice(zt(i),zfice(i)) ! inliner pour optimisation:
1855	zfice(i)=1.0-(zt(i)-pxtice)/(pxtmelt-pxtice)
1856	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1857	enddo
1858
1859	! calcul de la composition du condensat
1860	! hypothèse simplificatrice (mais à garder en mémoire!): l'eau
1861	! nuageuse a la même composition que l'eau totale de la maille:
1862	! c'est à dire zxtn=(zxt/zq)*zqn
1863
1864	call condiso_liq_ice_vectall(zxtn,zqn,zcond,zt,zfice, &
1865	& zxtice,zxtliq,klon)
1866	#ifdef ISOTRAC
1867	!#ifdef ISOVERIF
1868	! write(,) 'il pleut 980: condensation pour les traceurs'
1869	!#endif
1870	call condiso_liq_ice_vectall_trac(zxtn,zqn,zcond, &
1871	& zt,zfice,zxtice,zxtliq,klon)
1872	#ifdef ISOVERIF
1873	! write(,) 'ilp tmp 1859'
1874	do i=1,klon
1875	call iso_verif_traceur(zxtn(1,i),'il pleut 1428a')
1876	call iso_verif_traceur(zxtliq(1,i),'il pleut 1428b')
1877	call iso_verif_traceur(zxtice(1,i),'il pleut 1428c')
1878	enddo
1879	#endif
1880	#endif
1881	do i=1,klon
1882	do ixt=1,ntraciso
1883	zxtcond(ixt,i)=zxtice(ixt,i)+zxtliq(ixt,i)
1884	zxtcond(ixt,i)=max(zxtcond(ixt,i),0.0)
1885	zxtcond(ixt,i)=min(zxtcond(ixt,i),zxt(ixt,i))
1886	enddo !do ixt=1,niso
1887	enddo !do i=1,klon
1888
1889	do i=1,klon
1890	#ifdef DIAGISO
1891	! juste diag: fraction condensée
1892	if (zq(i).gt.ridicule) then
1893	fcond_ls(i,k)=zcond(i)/zq(i)
1894	endif
1895	! en juste diag
1896	#endif
1897	enddo !do i=1,klon
1898
1899	!#ifdef ISOVERIF
1900	#ifdef ISOVERIF
1901	do i=1,klon
1902	do ixt=1,ntraciso
1903	call iso_verif_noNAN(zxtcond(ixt,i),'il pleut 638')
1904	enddo
1905	enddo !do i=1,klon
1906	#endif
1907	#ifdef ISOVERIF
1908	if (iso_eau.gt.0) then
1909	do i=1,klon
1910	call iso_verif_egalite_choix(zxtcond(iso_eau,i),zcond(i), &
1911	& 'il pleut 645',errmax,errmaxrel)
1912	enddo ! do i=1,klon
1913	endif
1914	if (iso_HDO.gt.0) then
1915	do i=1,klon
1916	if (zcond(i).gt.ridicule) then
1917	call iso_verif_aberrant(zxtcond(iso_HDO,i) &
1918	& /zcond(i), 'il pleut 637')
1919	endif !if (zcond(i).gt.ridicule) then
1920	IF ((t(i,1)) .LT. RTT) THEN
1921	call iso_verif_aberrant_choix(zxtcond(iso_hdo,i), &
1922	& zcond(i),ridicule_rain,deltalim_snow,'il pleut 1276')
1923	endif
1924	enddo ! do i=1,klon
1925	endif !if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
1926	#ifdef ISOTRAC
1927	do i=1,klon
1928	if (iso_verif_traceur_nostop(zxtcond(1,i),'il pleut 1455') &
1929	& .eq.1) then
1930	write(,) 'zxtn(:,i)=',zxtn(:,i)
1931	write(,) 'zxt(:,i)=',zxt(:,i)
1932	write(,) 'zxtliq(:,i)=',zxtliq(:,i)
1933	write(,) 'zxtice(:,i)=',zxtice(:,i)
1934	stop
1935	endif
1936	enddo
1937	#endif
1938	#endif
1939
1940	! retranchage à zxt
1941	do i=1,klon
1942	do ixt=1,ntraciso
1943	zxt(ixt,i)=zxt(ixt,i)-zxtcond(ixt,i)
1944	enddo !do ixt=1,niso
1945	enddo !do i=1,klon
1946
1947	#ifdef ISOVERIF
1948	do i=1,klon
1949	do ixt=1,niso
1950	call iso_verif_positif(zxt(ixt,i),'il pleut 1088')
1951	enddo !do ixt=1,niso
1952	#ifdef ISOTRAC
1953	! on sépare pour mieux débugger 8 fev 2010
1954	do ixt=1,ntraciso
1955	call iso_verif_positif(zxt(ixt,i),'il pleut 1093')
1956	enddo !do ixt=1,niso
1957	#endif
1958	enddo !do i=1,klon
1959	#endif
1960
1961	#ifdef ISOTRAC
1962	if (option_tmin.ge.1) then
1963	do i=1,klon
1964	! colorier la vapeur résiduelle selon température de
1965	! condensation, et le condensat en un tag spécifique
1966	! Attention: zqn,zxtn ne servent qu'au calcul de zcond,zxtcond.
1967	! C'est en fait zxt qui se fait retranché
1968	! correction le 31 mars 2010: c'est à qn qu'il faut retrancher
1969	! le condensat, car la condensation LS est un processus sous
1970	! maille -> ne tagguer que le résidus de qn.
1971	if ((zcond(i).gt.0.0).and.(zq(i).gt.0.0)) then
1972
1973	if (rneb(i,k).gt.0.0) then
1974	zcondn(i)=zcond(i)/rneb(i,k)
1975	do ixt=1,ntraciso
1976	zxtcondn(ixt,i)=zxtcond(ixt,i)/rneb(i,k)
1977	enddo !do ixt=1,ntraciso
1978	else !if (rneb.eq.0.0)
1979	zcondn(i)=0.0
1980	do ixt=1,ntraciso
1981	zxtcondn(ixt,i)=0.0
1982	enddo !do ixt=1,ntraciso
1983	#ifdef ISOVERIF
1984	call iso_verif_egalite(zcond(i),0.0,'ilp 1225')
1985	#endif
1986	endif !if (rneb.eq.0.0)
1987
1988	#ifdef ISOVERIF
1989	call iso_verif_egalite_choix(zxtcondn(iso_eau,i),zcondn(i), &
1990	& 'il pleut 1244',errmax10,errmaxrel10)
1991	call iso_verif_traceur(zxtcondn(1,i),'cv3_routines 2302')
1992	zq_verif=rneb(i,k)zqn(i)+(1.0-rneb(i,k))zqcs(i)-zcond(i)
1993	call iso_verif_egalite_choix(zq(i),zq_verif, &
1994	& 'il pleut 1235',errmax,errmaxrel)
1995	do ixt=1,niso
1996	zq_verif=rneb(i,k)*zxtn(ixt,i) &
1997	& +(1.0-rneb(i,k))*zxtcs(ixt,i)-zxtcond(ixt,i)
1998	call iso_verif_egalite_choix(zxt(ixt,i),zq_verif, &
1999	& 'il pleut 1235',errmax,errmaxrel)
2000	enddo
2001	#endif
2002	! bidouille
2003	if (bidouille_anti_divergence) then
2004	zxtcondn(iso_eau,i)=zcondn(i)
2005	#ifdef ISOVERIF
2006	call iso_verif_traceur_pbidouille(zxtcondn(1,i), &
2007	& 'il pleut 1261')
2008	#endif
2009	endif
2010
2011	if (option_traceurs.eq.17) then
2012	call iso_recolorise_condensation(zqn(i),zcondn(i), &
2013	& zxtn(1,i),zxtcondn(1,i),zt(i),1.0,zxtres, &
2014	& seuil_tag_tmin_ls)
2015	else !if (option_traceurs.eq.17) then
2016	! cas du tag 18 ou 19
2017	call iso_recolorise_condensation(zqn(i),zcondn(i), &
2018	& zxtn(1,i),zxtcondn(1,i),zqs(i),1.0,zxtres, &
2019	& seuil_tag_tmin_ls)
2020	endif !if (option_traceurs.eq.17) then
2021	do ixt=1+niso,ntraciso
2022	zxt(ixt,i)=rneb(i,k)*zxtres(ixt) &
2023	& +(1.0-rneb(i,k))*zxtcs(ixt,i)
2024	enddo
2025	#ifdef ISOVERIF
2026	call iso_verif_traceur(zxt(1,i),'ilp 1025')
2027	do izone=izone_oce,izone_ddft
2028	if ((izone.eq.izone_oce).or.(izone.eq.izone_ddft)) then
2029	! avec condensation, on ne peut que perdre des traceurs à
2030	! part le tag résuel et le condensat
2031	if (iso_verif_positif_choix_nostop( &
2032	& zxt_ancien(index_trac(izone,iso_eau),i) &
2033	& -zxt(index_trac(izone,iso_eau),i),1e-8,'ilp 1270') &
2034	& .eq.1) then
2035	write(,) 'i,izone,rneb=',i,izone,rneb(i,k)
2036	write(,) 'zxt(:,i)=',zxt(:,i)
2037	write(,) 'zxt_ancien(:,i)=',zxt_ancien(:,i)
2038	write(,) 'zxtcs(:,i)=',zxtcs(:,i)
2039	write(,) 'zxtres(:)=',zxtres(:)
2040	write(,) 'zxtcond(:,i)=',zxtcond(:,i)
2041	write(,) 'zxtn(:,i)=',zxtn(:,i)
2042	write(,) 'zxtcondn(:,i)=',zxtcondn(:,i)
2043	stop
2044	endif !if (iso_verif_positif_nostop(
2045	endif !if ((izone.eq.izone_oce).or.(izone.eq.izone_ddft)) then
2046	enddo !do izone=izone_oce,izone_revap
2047	#endif
2048	endif !if (cond.gt.0.0) then
2049	enddo !do i=1,klon
2050
2051	#ifdef ISOVERIF
2052	! vérifier qu'on recolorise bien le condensat
2053	call iso_verif_positif(float(option_cond-1),'ilp 1310')
2054	#endif
2055	endif ! if (option_tmin.ge.1) then
2056
2057	if (option_cond.eq.1) then
2058	! colorier le condensat en un tag spécifique
2059	do i=1,klon
2060	do ixt=1+niso,ntraciso
2061	if (index_zone(ixt).eq.izone_cond) then
2062	zxtcond(ixt,i)=zxtcond(index_iso(ixt),i)
2063	else
2064	zxtcond(ixt,i)=0.0
2065	endif
2066	enddo !do ixt=1,ntraciso
2067	enddo !do i=1,klon
2068	#ifdef ISOVERIF
2069	do i=1,klon
2070	call iso_verif_egalite_choix(zxtcond(iso_eau,i),zcond(i), &
2071	& 'il pleut 1139',errmax,errmaxrel)
2072	call iso_verif_traceur(zxtcond(1,i),'cv3_routines 2302')
2073	enddo !do i=1,klon
2074	#endif
2075	endif ! if (option_cond.ge.1) then
2076	#endif
2077
2078
2079	#ifdef ISOVERIF
2080	! write(,) 'ilp tmp 2066'
2081	do i=1,klon
2082	do ixt=1,niso
2083	call iso_verif_positif(zxt(ixt,i),'il pleut 778')
2084	enddo !do ixt=1,niso
2085	#ifdef ISOTRAC
2086	! on sépare pour mieux débugger 8 fev 2010
2087	do ixt=1,ntraciso
2088	call iso_verif_positif(zxt(ixt,i),'il pleut 1135')
2089	call iso_verif_traceur(zxt(1,i),'il pleut 1485')
2090	enddo !do ixt=1,niso
2091	#endif
2092	if (iso_eau.gt.0) then
2093	call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
2094	& 'il pleut 644',errmax,errmaxrel)
2095	endif !if (iso_eau.gt.0) then
2096	if (iso_HDO.gt.0) then
2097	if (zq(i).gt.ridicule) then
2098	call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
2099	& 'il pleut 648')
2100	if (iso_O18.gt.0) then
2101	if (iso_verif_O18_aberrant_nostop(zxt(iso_HDO,i)/zq(i), &
2102	& zxt(iso_O18,i)/zq(i),'il pleut 2059').eq.1) then
2103	write(,) 'i,k,zq=',i,k,zq(i)
2104	stop
2105	endif !if (iso_verif_O18_aberrant_nostop
2106	endif !if (iso_O18.gt.0) then
2107	endif !if (zq(i).gt.ridicule) then
2108	endif ! if (iso_HDO.gt.0) then
2109	enddo ! i=1,klon
2110	#ifdef ISOTRAC
2111	do i=1,klon
2112	call iso_verif_traceur(zxt(1,i),'il pleut 1485')
2113	enddo
2114	#endif
2115	#endif
2116	#ifdef ISOVERIF
2117	do i=1,klon
2118	do ixt=1,ntraciso
2119	call iso_verif_noNAN(zxtcond(ixt,i),'il pleut 1560')
2120	call iso_verif_noNAN(zxt(ixt,i),'il pleut 1561')
2121	enddo ! ixt=1,niso
2122	enddo !do i=1,klon
2123	#endif
2124	#endif
2125	!>AJ
2126	! ----------------------------------------------------------------
2127	! P3> Formation des precipitations
2128	! ----------------------------------------------------------------
2129	!
2130	! Partager l'eau condensee en precipitation et eau liquide nuageuse
2131	!
2132
2133	! Initialisation de zoliq (eau condensee moyenne dans la maille)
2134	DO i = 1, klon
2135	IF (rneb(i,k).GT.0.0) THEN
2136	zoliq(i) = zcond(i)
2137	zrho(i) = pplay(i,k) / zt(i) / RD
2138	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
2139	ENDIF
2140	ENDDO
2141	!AJ<
2142	IF (.NOT. ice_thermo) THEN
2143	IF (iflag_t_glace.EQ.0) THEN
2144	DO i = 1, klon
2145	IF (rneb(i,k).GT.0.0) THEN
2146	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
2147	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
2148	zfice(i) = zfice(i)**exposant_glace_old
2149	! zfice(i) = zfice(i)**nexpo
2150	!! zfice(i)=0.
2151	ENDIF
2152	ENDDO
2153	ELSE ! of IF (iflag_t_glace.EQ.0)
2154	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
2155	! DO i = 1, klon
2156	! IF (rneb(i,k).GT.0.0) THEN
2157	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
2158	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
2159	! t_glace_max, exposant_glace)
2160	! ENDIF
2161	! ENDDO
2162	ENDIF
2163	ENDIF
2164
2165	! Calcul de radliq (eau condensee pour le rayonnement)
2166	! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
2167	! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
2168	! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
2169	! transmise au rayonnement;
2170	! ----------------------------------------------------------------
2171	DO i = 1, klon
2172	IF (rneb(i,k).GT.0.0) THEN
2173	zneb(i) = MAX(rneb(i,k), seuil_neb)
2174	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
2175	! print*,zt(i),'fractionglace'
2176	!>AJ
2177	radliq(i,k) = zoliq(i)/REAL(ninter+1)
2178	ENDIF
2179	ENDDO
2180	#ifdef ISO
2181	! initilsation des quantités totales d'eau liquide et solide qui
2182	! tombent
2183	do i=1,klon
2184	if (rneb(i,k).gt.0.0) then
2185	do ixt=1,ntraciso
2186	zxtoliq(ixt,i)=zxtcond(ixt,i)
2187	enddo
2188	endif ! if (rneb(i,k).gt.0.0)
2189	enddo
2190
2191	#ifdef ISOVERIF
2192	! write(,) 'ilp tmp 1626: k,zoliq,zxtoliq(1,363)=',k,zoliq(363),zxtoliq(1,363)
2193	do i=1,klon
2194	do ixt=1,ntraciso
2195	if (iso_verif_noNAN_nostop(zxtoliq(ixt,i), &
2196	& 'il pleut 1628').eq.1) then
2197	write(,) 'i,k,ixt=',i,k,ixt
2198	write(,) 'rneb(i,k)=',rneb(i,k)
2199	write(,) 'zxtcond(ixt,i)=',zxtcond(ixt,i)
2200	stop
2201	endif
2202	enddo ! ixt=1,niso
2203	enddo !do i=1,klon
2204	#endif
2205	#ifdef ISOVERIF
2206	!i=519
2207	!write(,) 'k,i,zoliq(i)=',k,i,zoliq(i)
2208	do i=1,klon
2209	! cam verif
2210	if (iso_eau.gt.0) then
2211	call iso_verif_egalite_choix(zxtoliq(iso_eau,i),zoliq(i), &
2212	& 'il pleut 678',errmax,errmaxrel)
2213	endif !if (iso_eau.gt.0) then
2214	if (zoliq(i).gt.ridicule) then
2215	if (iso_HDO.gt.0) then
2216	call iso_verif_aberrant(zxtoliq(iso_HDO,i)/zoliq(i), &
2217	& 'il pleut 895a')
2218	if (iso_O18.gt.0) then
2219	if (iso_verif_O18_aberrant_nostop(zxtoliq(iso_HDO,i)/zoliq(i), &
2220	& zxtoliq(iso_O18,i)/zoliq(i),'il pleut 895b').eq.1) then
2221	write(,) 'i,k,zoliq,zfice=',i,k,zoliq(i),zfice(i)
2222	stop
2223	endif
2224	endif ! if (iso_HDO.gt.0) then
2225	endif !if (iso_HDO.gt.0) then
2226	endif !if (zoliq(i).gt.ridicule) then
2227	! end cam verif
2228	enddo ! do i=1,klon
2229	#endif
2230	#endif
2231
2232	!
2233	DO n = 1, ninter
2234	DO i = 1, klon
2235	IF (rneb(i,k).GT.0.0) THEN
2236	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
2237	! Initialization of zpluie and zice:
2238	zpluie=0
2239	zice=0
2240	IF (zneb(i).EQ.seuil_neb) THEN
2241	ztot = 0.0
2242	ELSE
2243	! quantite d'eau a eliminer: zchau (Sundqvist, 1978)
2244	! meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
2245	if (ptconv(i,k)) then
2246	zcl =cld_lc_con
2247	zct =1./cld_tau_con
2248	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
2249	fallvc(zrhol(i)) zfice(i)
2250	else
2251	zcl =cld_lc_lsc
2252	zct =1./cld_tau_lsc
2253	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
2254	fallvs(zrhol(i)) zfice(i)
2255	endif
2256	zchau = zct dtime/REAL(ninter) zoliq(i) &
2257	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
2258	!AJ<
2259	IF (.NOT. ice_thermo) THEN
2260	ztot = zchau + zfroi
2261	ELSE
2262	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
2263	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
2264	ztot = zpluie + zice
2265	ENDIF
2266	!>AJ
2267	ztot = MAX(ztot ,0.0)
2268	ENDIF
2269	ztot = MIN(ztot,zoliq(i))
2270	!AJ<
2271	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
2272	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
2273	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
2274	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
2275	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
2276	!>AJ
2277	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
2278	ENDIF
2279	ENDDO ! i = 1,klon
2280	ENDDO ! n = 1,ninter
2281	! ----------------------------------------------------------------
2282	!
2283	#ifdef ISO
2284	do i=1,klon
2285	IF (rneb(i,k).GT.0.0) THEN
2286	if (zcond(i).gt.ridicule**2) then
2287	! le 21 dec: on change .gt.0 en .gt.ridicule**2 pour éviter valeur
2288	! ridiculement petites
2289	do ixt=1,ntraciso
2290	zxtoliq(ixt,i)=zoliq(i)*(zxtcond(ixt,i)/zcond(i))
2291	enddo !do ixt=1,niso
2292	else !if (zcond(i).gt.0.0) then
2293	if (zoliq(i).lt.ridicule) then
2294	do ixt=1,ntraciso
2295	zxtoliq(ixt,i)=0.0
2296	enddo
2297	else !if (zoliq(i).lt.ridicule) then
2298	! modif 28 oct 2008
2299	do ixt=1,niso
2300	zxtoliq(ixt,i)=zoliq(i)*Rdefault(ixt)
2301	enddo !do ixt=1,niso
2302	#ifdef ISOTRAC
2303	do ixt=niso+1,ntraciso
2304	if (index_zone(ixt).eq.izone_poubelle) then
2305	zxtoliq(ixt,i)=zoliq(i)*Rdefault(ixt)
2306	else
2307	zxtoliq(ixt,i)=0.0
2308	endif
2309	enddo
2310	#endif
2311	#ifdef ISOVERIF
2312	write(,) 'il pleut 1412: zcond(i)=',zcond(i)
2313	write(,) 'zoliq(i)=',zoliq(i)
2314	stop
2315	#endif
2316	endif !if (zoliq(i).lt.ridicule) then
2317	endif !if (zcond(i).gt.0.0) then
2318	endif !IF (rneb(i,k).GT.0.0) THEN
2319	enddo ! do i=1,klon
2320
2321	! cam verif
2322	#ifdef ISOVERIF
2323	do i=1,klon
2324	IF (rneb(i,k).GT.0.0) THEN
2325	do ixt=1,niso
2326	call iso_verif_noNAN(zxtoliq(ixt,i),'il pleut 718')
2327	! end cam verif
2328	enddo ! do ixt=1,niso
2329	! cam verif
2330	endif !IF (rneb(i,k).GT.0.0) THEN
2331	enddo !do i=1,klon
2332	#endif
2333	#ifdef ISOVERIF
2334	do i=1,klon
2335	IF (rneb(i,k).GT.0.0) THEN
2336	if (iso_eau.gt.0) then
2337	call iso_verif_egalite_choix(zxtoliq(iso_eau,i),zoliq(i), &
2338	& 'il pleut 749',errmax,errmaxrel)
2339	endif !if (iso_eau.gt.0) then
2340	if (iso_HDO.gt.0) then
2341	if (zoliq(i).gt.ridicule) then
2342	call iso_verif_aberrant(zxtoliq(iso_HDO,i)/zoliq(i), &
2343	& 'il pleut 963')
2344	endif
2345	endif !if (iso_HDO.gt.0) then
2346	! end cam verif
2347	#ifdef ISOTRAC
2348	! write(,) 'ilp 1201 tmp: i=',i
2349	! write(,) 'zoliq(i)=',zoliq(i)
2350	! write(,) 'zxtoliq(i)=',zxtoliq(iso_eau:ntraciso:3,i)
2351	! write(,) 'zcond(i)=',zcond(i)
2352	! write(,) 'zxtcond(i)=',zxtcond(iso_eau:ntraciso:3,i)
2353	call iso_verif_traceur(zxtoliq(1,i),'il pleut 1634')
2354	#endif
2355	endif !IF (rneb(i,k).GT.0.0) THEN
2356	enddo ! do i=1,klon
2357	#endif
2358	#endif
2359
2360	IF (.NOT. ice_thermo) THEN
2361	DO i = 1, klon
2362	IF (rneb(i,k).GT.0.0) THEN
2363	d_ql(i,k) = zoliq(i)
2364	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
2365	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
2366
2367	#ifdef ISO
2368	#ifdef DIAGISO
2369	taux_cond_ls(i,k)=zcond(i)(paprs(i,k)-paprs(i,k+1))/(RGdtime)
2370	taux_precip_ls(i,k)=MAX(zcond(i)-zoliq(i),0.0) &
2371	& (paprs(i,k)-paprs(i,k+1))/(RGdtime)
2372	#endif
2373	do ixt=1,ntraciso
2374	d_xtl(ixt,i,k)=zxtoliq(ixt,i)
2375	enddo ! do ixt=1,niso
2376	! Rq: 28 oct 2008: si zoliq a la même compo que zcond, alors
2377	! on évite pb numériques si on prend directement zrfl avec
2378	! la compo de zcond?
2379	! do ixt=1,niso
2380	! zxtrfl(ixt,i)=zxtrfl(ixt,i) &
2381	! & +max(zxtcond(ixt,i)-zxtoliq(ixt,i),0.0) &
2382	! & * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
2383	! enddo
2384	if (zcond(i).gt.ridicule**2) then
2385	! modif le 19 dec 2012 pour éviter zcond pathologiquement petits
2386	do ixt=1,ntraciso
2387	zxtrfl(ixt,i)=zxtrfl(ixt,i) &
2388	& +(zxtcond(ixt,i)/zcond(i)) &
2389	& * MAX(zcond(i)-zoliq(i),0.0) &
2390	& * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
2391	enddo !do ixt=1,niso
2392	else !if (zcond(i).gt.0.0) then
2393	if (zcond(i)-zoliq(i).gt.0.0) then
2394	do ixt=1,niso
2395	zxtrfl(ixt,i)=zxtrfl(ixt,i) &
2396	& +Rdefault(ixt) &
2397	& * MAX(zcond(i)-zoliq(i),0.0) &
2398	& * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
2399	enddo !do ixt=1,niso
2400	#ifdef ISOVERIF
2401	write(,) 'il pleut 1023: zcond(i)=',zcond(i)
2402	write(,) 'zcond(i),zoliq(i)=',zcond(i),zoliq(i)
2403	stop
2404	#endif
2405	#ifdef ISOTRAC
2406	do ixt=niso+1,ntraciso
2407	if (index_zone(ixt).eq.izone_poubelle) then
2408	zxtrfl(ixt,i)=zxtrfl(ixt,i) &
2409	& +Rdefault(ixt) &
2410	& * MAX(zcond(i)-zoliq(i),0.0) &
2411	& * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
2412	endif
2413	enddo !do ixt=1,niso
2414	#endif
2415	else !if (zcond(i)-zoliq(i).gt.0.0) then
2416	! zrfl non modifié
2417	endif !if (zcond(i)-zoliq(i).gt.0.0) then
2418	endif !if (zcond(i).gt.0.0) then
2419
2420	! cam verif
2421	!#ifdef ISOVERIF
2422	#ifdef ISOVERIF
2423	do ixt=1,ntraciso
2424	if ((iso_verif_noNAN_nostop(zxtrfl(ixt,i), &
2425	& 'il pleut 772').eq.1).or. &
2426	& (iso_verif_noNAN_nostop(zxtoliq(ixt,i), &
2427	& 'il pleut 772b').eq.1)) then
2428	write(,) 'zcond(i)=',zcond(i)
2429	write(,) 'zxtcond(ixt,i)=',zxtcond(ixt,i)
2430	write(,) 'zoliq(i)=',zoliq(i)
2431	write(,) 'Rdefault(ixt)=',Rdefault(ixt)
2432	write(,) 'paprs(i,k),paprs(i,k+1)=',paprs(i,k),paprs(i,k+1)
2433	stop
2434	endif
2435	enddo
2436	#endif
2437	#ifdef ISOVERIF
2438	if (iso_eau.gt.0) then
2439	call iso_verif_egalite_choix(d_xtl(iso_eau,i,k), &
2440	& d_ql(i,k),'il pleut 771',errmax,errmaxrel)
2441	if (iso_verif_egalite_choix_nostop( &
2442	& zxtrfl(iso_eau,i),zrfl(i), &
2443	& 'il pleut 771b',errmax,errmaxrel).eq.1)then
2444	write(,) 'trace_cas(i)=',trace_cas(i)
2445	write(,) 'zxtcond(iso_eau,i)',zxtcond(iso_eau,i)
2446	write(,) 'zcond(i)',zcond(i)
2447	write(,) 'zxtoliq(iso_eau,i)',zxtoliq(iso_eau,i)
2448	write(,) 'zoliq(i)',zoliq(i)
2449	stop
2450	endif
2451	endif !if (iso_eau.gt.0) then
2452	if (iso_HDO.gt.0) then
2453	IF ((t(i,1)) .LT. RTT) THEN
2454	call iso_verif_aberrant_choix(zxtrfl(iso_hdo,i), &
2455	& zrfl(i),ridicule_rain,deltalim_snow,'il pleut 1458')
2456	endif
2457	endif !if (iso_HDO.gt.0) then
2458	#ifdef ISOTRAC
2459	call iso_verif_traceur(zxtrfl(1,i),'il pleut 1729')
2460	#endif
2461	#endif
2462	! bidouille
2463	if ((bidouille_anti_divergence).and.(iso_eau.gt.0)) then
2464	zxtrfl(iso_eau,i)= zrfl(i)
2465	! if (zxtrfl(iso_eau,i).gt.0.0) then
2466	! do ixt=1,niso
2467	! zxtrfl(ixt,i)=zxtrfl(ixt,i)*zrfl(i)/zxtrfl(iso_eau,i)
2468	! enddo
2469	! endif !if (zxtrfl(iso_eau,i).gt.0.0) then
2470	endif ! if (bidouille_anti_divergence) then
2471	! end bidouille
2472	! end cam verif
2473	#endif
2474	ENDIF !IF (rneb(i,k).GT.0.0) THEN
2475	ENDDO
2476	ELSE
2477
2478	#ifdef ISO
2479	CALL abort_physic('ilp 2347', 'isotopes pas encore prevue ici', 1)
2480	#endif
2481	!
2482	!CR&JYG<
2483	! On prend en compte l'effet Bergeron dans les flux de precipitation :
2484	! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
2485	! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
2486	! et les precipitations est grossierement pris en compte en linearisant les equations
2487	! et en approximant le processus de precipitation liquide par un processus a seuil.
2488	! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
2489	! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
2490	! Le condensat precipitant solide est augmente.
2491	! La vapeur d'eau est augmentee.
2492	!
2493	IF ((iflag_bergeron .EQ. 2)) THEN
2494	DO i = 1, klon
2495	IF (rneb(i,k) .GT. 0.0) THEN
2496	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
2497	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
2498	zcp=RCPD(1.0+RVTMP2(zq(i)+zcond(i)))
2499	coef1 = RLMLT*zdqs(i)/RLVTT
2500	DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
2501	zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
2502	zqprecl(i) = max( zqprecl(i) - zcp/RLMLT(1.+coef1)DeltaT, 0. )
2503	zcond(i) = max( zcond(i) - zcp/RLVTTzdqs(i)DeltaT, 0. )
2504	zq(i) = zq(i) + zcp/RLVTTzdqs(i)DeltaT
2505	zt(i) = zt(i) + DeltaT
2506	ENDIF ! rneb(i,k) .GT. 0.0
2507	ENDDO
2508	DO i = 1, klon
2509	IF (rneb(i,k).GT.0.0) THEN
2510	d_ql(i,k) = (1-zfice(i))*zoliq(i)
2511	d_qi(i,k) = zfice(i)*zoliq(i)
2512	zrfl(i) = zrfl(i)+ zqprecl(i) &
2513	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
2514	zifl(i) = zifl(i)+ zqpreci(i) &
2515	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
2516	ENDIF
2517	ENDDO
2518	!!
2519	ELSE ! iflag_bergeron
2520	!>CR&JYG
2521	!!
2522	DO i = 1, klon
2523	IF (rneb(i,k).GT.0.0) THEN
2524	!CR on prend en compte la phase glace
2525	if (.not.ice_thermo) then
2526	d_ql(i,k) = zoliq(i)
2527	d_qi(i,k) = 0.
2528	else
2529	d_ql(i,k) = (1-zfice(i))*zoliq(i)
2530	d_qi(i,k) = zfice(i)*zoliq(i)
2531	endif
2532	!AJ<
2533	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
2534	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
2535	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
2536	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
2537	! zrfl(i) = zrfl(i)+ zpluie &
2538	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
2539	! zifl(i) = zifl(i)+ zice &
2540	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
2541
2542	!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
2543	IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
2544	zsolid = zrfl(i)
2545	zifl(i) = zifl(i)+zrfl(i)
2546	zrfl(i) = 0.
2547	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
2548	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))
2549	ENDIF ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
2550	!RC
2551
2552	ENDIF ! rneb(i,k).GT.0.0
2553	ENDDO
2554
2555	ENDIF ! iflag_bergeron .EQ. 2
2556	ENDIF ! .NOT. ice_thermo
2557
2558	!CR: la fonte est faite au debut
2559	! IF (ice_thermo) THEN
2560	! DO i = 1, klon
2561	! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
2562	! zmelt = MIN(MAX(zmelt,0.),1.)
2563	! zrfl(i)=zrfl(i)+zmelt*zifl(i)
2564	! zifl(i)=zifl(i)*(1.-zmelt)
2565	! print*,zt(i),'octavio1'
2566	! zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
2567	! RLMLT/RCPD/(1.0+RVTMP2zq(i))
2568	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
2569	! ENDDO
2570	! ENDIF
2571
2572
2573	IF (.NOT. ice_thermo) THEN
2574	DO i = 1, klon
2575	IF (zt(i).LT.RTT) THEN
2576	psfl(i,k)=zrfl(i)
2577	#ifdef ISO
2578	do ixt=1,ntraciso
2579	pxtsnowfl(ixt,i,k)=zxtrfl(ixt,i)
2580	enddo
2581	#endif
2582	ELSE
2583	prfl(i,k)=zrfl(i)
2584	#ifdef ISO
2585	do ixt=1,ntraciso
2586	pxtrainfl(ixt,i,k)=zxtrfl(ixt,i)
2587	enddo
2588	#endif
2589	ENDIF
2590	ENDDO
2591	ELSE ! if ice_thermo
2592	! JAM*************************************************
2593	! Revoir partie ci-dessous: a quoi servent psfl et prfl?
2594	! *****************************************************
2595
2596	#ifdef ISO
2597	CALL abort_physic('ilp 2465', 'isos pas prevus ici', 1)
2598	#endif
2599	DO i = 1, klon
2600	! IF (zt(i).LT.RTT) THEN
2601	psfl(i,k)=zifl(i)
2602	! ELSE
2603	prfl(i,k)=zrfl(i)
2604	! ENDIF
2605	!>AJ
2606	ENDDO
2607	ENDIF
2608
2609	#ifdef ISO
2610	#ifdef ISOVERIF
2611	DO i = 1, klon
2612	if (iso_eau.gt.0) then
2613	call iso_verif_egalite_choix(pxtrainfl(iso_eau,i,k),prfl(i,k), &
2614	& 'il pleut 804',errmax,errmaxrel)
2615	call iso_verif_egalite_choix(pxtsnowfl(iso_eau,i,k),psfl(i,k), &
2616	& 'il pleut 805',errmax,errmaxrel)
2617	call iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &
2618	& 'il pleut 1612',errmax,errmaxrel)
2619	if (iso_verif_positif_choix_nostop(pxtrainfl(iso_eau,i,k), &
2620	& ridicule_rain,'ilp 1547').eq.1) then
2621	write(,) 'prfl(i,k)=',prfl(i,k)
2622	stop
2623	endif
2624	endif !if (iso_eau.gt.0) then
2625	if (iso_HDO.gt.0) then
2626	if (prfl(i,k).gt.ridicule_rain) then
2627	call iso_verif_aberrant(pxtrainfl(iso_HDO,i,k)/prfl(i,k), &
2628	& 'il pleut 1020')
2629	endif !if (prfl(i,k).gt.ridicule_rain) then
2630	if (psfl(i,k).gt.ridicule_rain) then
2631	call iso_verif_aberrant(pxtsnowfl(iso_HDO,i,k)/psfl(i,k), &
2632	& 'il pleut 1024')
2633	endif !if (psfl(i,k).gt.ridicule_rain) then
2634	if (zq(i).gt.ridicule) then
2635	call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
2636	& 'il pleut 1204')
2637	endif !if (zq(i).gt.ridicule) then
2638	endif !if (iso_HDO.gt.0) then
2639	#ifdef ISOTRAC
2640	call iso_verif_traceur(zxt(1,i),'il pleut 1792')
2641	call iso_verif_traceur(pxtsnowfl(1,i,k),'il pleut 1793')
2642	call iso_verif_traceur(pxtrainfl(1,i,k),'il pleut 1794')
2643	#endif
2644	ENDDO !DO i = 1, klon
2645	#endif
2646	#endif
2647	! ----------------------------------------------------------------
2648	! Fin de formation des precipitations
2649	! ----------------------------------------------------------------
2650	!
2651	!
2652	! Calculer les tendances de q et de t:
2653	!
2654	DO i = 1, klon
2655	d_q(i,k) = zq(i) - q(i,k)
2656	#ifdef ISO
2657	do ixt=1,ntraciso
2658	d_xt(ixt,i,k)=zxt(ixt,i)-xt(ixt,i,k)
2659	enddo
2660	! cam verif
2661	#ifdef ISOVERIF
2662	if (iso_eau.gt.0) then
2663	if (iso_verif_egalite_choix_nostop(d_xt(iso_eau,i,k), &
2664	& d_q(i,k),'il pleut 838',errmax,errmaxrel).eq.1) then
2665	write(,) 'i,k=',i,k
2666	write(,) 'zxt,zq=',zxt(iso_eau,i),zq(i)
2667	write(,) 'xt,q=',xt(iso_eau,i,k),q(i,k)
2668	stop
2669	endif
2670	endif !if (iso_eau.gt.0) then
2671	if (1.eq.0) then
2672	if (iso_HDO.gt.0) then
2673	if (abs(d_q(i,k)).gt.max(errmaxrel*q(i,k),errmax)) then
2674	if (iso_verif_aberrant_nostop( &
2675	& d_xt(iso_HDO,i,k)/d_q(i,k),'il pleut 1229').eq.1) then
2676	! write(,) 'd_q(i,k)=',d_q(i,k)
2677	! write(,) 'q(i,k)=',q(i,k)
2678	! write(,) 'zq(i)=',zq(i)
2679	! write(,) 'deltaD_prec=',deltaD(xt(iso_HDO,i,k)/q(i,k))
2680	! write(,) 'deltaD=',deltaD(zxt(iso_HDO,i)/zq(i))
2681	! stop
2682	endif ! if (iso_verif_aberrant_nostop(
2683	endif !if (zq(i).gt.ridicule) then
2684	endif !if (iso_HDO.gt.0) then
2685	endif !if (1.eq.0) then
2686	#endif
2687	! end cam verif
2688	#endif
2689	d_t(i,k) = zt(i) - t(i,k)
2690	ENDDO
2691	!
2692	!AA--------------- Calcul du lessivage stratiforme -------------
2693
2694	DO i = 1,klon
2695	!
2696	if(zcond(i).gt.zoliq(i)+1.e-10) then
2697	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
2698	else
2699	beta(i,k) = 0.
2700	endif
2701	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
2702	* (paprs(i,k)-paprs(i,k+1))/RG
2703	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
2704	!AA lessivage nucleation LMD5 dans la couche elle-meme
2705	IF (iflag_t_glace.EQ.0) THEN
2706	if (t(i,k) .GE. t_glace_min_old) THEN
2707	zalpha_tr = a_tr_sca(3)
2708	else
2709	zalpha_tr = a_tr_sca(4)
2710	endif
2711	ELSE ! of IF (iflag_t_glace.EQ.0)
2712	if (t(i,k) .GE. t_glace_min) THEN
2713	zalpha_tr = a_tr_sca(3)
2714	else
2715	zalpha_tr = a_tr_sca(4)
2716	endif
2717	ENDIF
2718	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
2719	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
2720	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
2721	!
2722	! nucleation avec un facteur -1 au lieu de -0.5
2723	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
2724	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
2725	ENDIF
2726	!
2727	ENDDO ! boucle sur i
2728	!
2729	!AA Lessivage par impaction dans les couches en-dessous
2730	DO kk = k-1, 1, -1
2731	DO i = 1, klon
2732	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
2733	IF (iflag_t_glace.EQ.0) THEN
2734	if (t(i,kk) .GE. t_glace_min_old) THEN
2735	zalpha_tr = a_tr_sca(1)
2736	else
2737	zalpha_tr = a_tr_sca(2)
2738	endif
2739	ELSE ! of IF (iflag_t_glace.EQ.0)
2740	if (t(i,kk) .GE. t_glace_min) THEN
2741	zalpha_tr = a_tr_sca(1)
2742	else
2743	zalpha_tr = a_tr_sca(2)
2744	endif
2745	ENDIF
2746	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
2747	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
2748	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
2749	ENDIF
2750	ENDDO
2751	ENDDO !DO kk = k-1, 1, -1
2752
2753	#ifdef ISO
2754	! verif en fin de boucle
2755	#ifdef ISOVERIF
2756	if (iso_eau.gt.0) then
2757	do i=1,klon
2758	call iso_verif_egalite_choix(zxtrfl(iso_eau,i) &
2759	& ,zrfl(i), 'il pleut 1089',errmax,errmaxrel)
2760	call iso_verif_egalite_choix(d_xtl(iso_eau,i,k),d_ql(i,k), &
2761	& 'il pleut 1205',errmax,errmaxrel)
2762	#ifdef ISOTRAC
2763	call iso_verif_traceur(zxtrfl(1,i),'il pleut 1894')
2764	call iso_verif_traceur(d_xtl(1,i,k),'il pleut 1895')
2765	#endif
2766	enddo !do i=1,klon
2767	endif ! if (iso_eau.gt.0) then
2768	#endif
2769	#endif
2770	! end verif en fin de boucle
2771
2772	!
2773	!AA===============================================================
2774	! FIN DE LA BOUCLE VERTICALE
2775	end DO
2776	!
2777	!AA==================================================================
2778	!
2779	! Pluie ou neige au sol selon la temperature de la 1ere couche
2780	!
2781	!CR: si la thermo de la glace est active, on calcule zifl directement
2782	IF (.NOT.ice_thermo) THEN
2783	DO i = 1, klon
2784	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
2785	!AJ<
2786	! snow(i) = zrfl(i)
2787	snow(i) = zrfl(i)+zifl(i)
2788	!>AJ
2789	zlh_solid(i) = RLSTT-RLVTT
2790	#ifdef ISO
2791	do ixt=1,ntraciso
2792	xtsnow(ixt,i) = zxtrfl(ixt,i)
2793	enddo !do ixt=1,ntraciso
2794	#endif
2795	ELSE
2796	rain(i) = zrfl(i)
2797	zlh_solid(i) = 0.
2798	#ifdef ISO
2799	do ixt=1,ntraciso
2800	xtrain(ixt,i) = zxtrfl(ixt,i)
2801	enddo !do ixt=1,ntraciso
2802	#endif
2803	ENDIF
2804	#ifdef ISO
2805	! cam verif
2806	#ifdef ISOVERIF
2807	do ixt=1,ntraciso
2808	call iso_verif_noNAN(xtrain(ixt,i),'il pleut 927')
2809	call iso_verif_noNAN(xtsnow(ixt,i),'il pleut 927')
2810	enddo ! do ixt=1,niso
2811	#endif
2812	#ifdef ISOVERIF
2813	if (iso_eau.gt.0) then
2814	call iso_verif_egalite_choix(xtrain(iso_eau,i),rain(i), &
2815	& 'il pleut 926a',errmax,errmaxrel)
2816	if (snow(i)*dtime.gt.ridicule) then
2817	call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
2818	& 'il pleut 938a',errmax,errmaxrel)
2819	endif !if (snow(i)*dtime.gt.ridicule) then
2820	endif !if (iso_eau.gt.0) then
2821	if (iso_HDO.gt.0) then
2822	call iso_verif_aberrant_choix(xtrain(iso_hdo,i),rain(i), &
2823	& ridicule_rain,deltalim,'il pleut 926b')
2824	call iso_verif_aberrant_choix(xtsnow(iso_hdo,i),snow(i), &
2825	& ridicule_rain,deltalim_snow,'il pleut 938b')
2826	endif !if (iso_eau.gt.0) then
2827	#ifdef ISOTRAC
2828	call iso_verif_traceur(xtrain(1,i),'il pleut 1952')
2829	call iso_verif_traceur(xtsnow(1,i),'il pleut 1953')
2830	#endif
2831	#endif
2832	! end cam verif
2833	#endif
2834
2835	ENDDO
2836
2837	ELSE
2838	#ifdef ISO
2839	CALL abort_physic('ilp 2708', 'isos pas prevus ici', 1)
2840	#endif
2841	DO i = 1, klon
2842	snow(i) = zifl(i)
2843	rain(i) = zrfl(i)
2844	ENDDO
2845
2846	ENDIF
2847	!
2848	! For energy conservation : when snow is present, the solification
2849	! latent heat is considered.
2850	!CR: si thermo de la glace, neige deja prise en compte
2851	IF (.not.ice_thermo) THEN
2852	DO k = 1, klev
2853	DO i = 1, klon
2854	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
2855	zmair=(paprs(i,k)-paprs(i,k+1))/RG
2856	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
2857	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair)
2858	END DO
2859	END DO
2860	ENDIF
2861	!
2862
2863	if (ncoreczq>0) then
2864	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
2865	endif
2866
2867
2868	#ifdef ISO
2869	! verif en sortie de il pleut
2870	#ifdef ISOVERIF
2871	do k=1,klev
2872	do i=1,klon
2873	if (iso_eau.gt.0) then
2874	call iso_verif_egalite_choix(d_xtl(iso_eau,i,k),d_ql(i,k), &
2875	& 'il pleut 1289',errmax,errmaxrel)
2876	endif !if (iso_eau.gt.0) then
2877	#ifdef ISOTRAC
2878	call iso_verif_traceur(d_xtl(1,i,k),'il pleut 1982')
2879	#endif
2880	enddo !do i=1,klon
2881	enddo !do k=1,klev
2882	#endif
2883	! end verif
2884	#endif
2885
2886	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format