Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

fisrtilp.F90 @ 5441

Last change on this file since 5441 was 3331, checked in by acozic, 7 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 99.0 KB

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

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

Download in other formats:

Original Format