Context Navigation

fisrtilp.F90 @ 5454

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

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

Download in other formats:

Original Format