Context Navigation

fisrtilp.F90 @ 4500

Last change on this file since 4500 was 4491, checked in by crisi, 18 months ago
Bug corrections in LMDZiso, especially for water tagging
File size: 101.7 KB

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

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

Download in other formats:

Original Format