Context Navigation

lmdz_lscp_old.F90 @ 5213

Last change on this file since 5213 was 5209, checked in by Laurent Fairhead, 7 weeks ago
Once again modifications were not applied to the iso library
File size: 102.5 KB

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

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

Download in other formats:

Original Format