Context Navigation

fisrtilp.F90 @ 4038

Last change on this file since 4038 was 3927, checked in by Laurent Fairhead, 3 years ago
Initial import of the physics wih isotopes from Camille Risi CR
File size: 102.1 KB

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

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

Download in other formats:

Original Format