Context Navigation

fisrtilp.F90 @ 5435

Last change on this file since 5435 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
File size: 102.0 KB

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

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

Download in other formats:

Original Format