Context Navigation

fisrtilp.F90 @ 4766

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

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

Download in other formats:

Original Format