Context Navigation

lmdz_lscp_old.F90 @ 5226

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

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

Download in other formats:

Original Format