Context Navigation

lmdz_lscp_old.F90 @ 5454

Last change on this file since 5454 was 5285, checked in by abarral, 2 months ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 102.5 KB

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

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

Download in other formats:

Original Format