Context Navigation

lmdz_lscp_old.F90 @ 5119

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

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

Download in other formats:

Original Format