Context Navigation

fisrtilp.F90 @ 4657

Last change on this file since 4657 was 4651, checked in by fhourdin, 16 months ago

Pre-replayisation pour cloudth

Cration de lmdz_cloudth_ini
cloudth_mod.F90 -> lmdz_cloudth.F90
Supression des constantes dans lmdz_cloudth
En fait, cloudth lui même n'est pas vraiment replayisable car appeler ligne par ligne.

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

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

Download in other formats:

Original Format