Context Navigation

fisrtilp.F90 @ 2957

Last change on this file since 2957 was 2956, checked in by jbmadeleine, 8 years ago
Corrected rneblsvol in fisrtilp to make it equal to rnebls in the part controlled by the large-scale lognormal PDF. Nb: rneblsvol stays equal to 0 if iflag_cloudth_vert<3.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 56.7 KB

Line
1	!
2	! $Id: fisrtilp.F90 2956 2017-07-17 22:01:04Z jbmadeleine $
3	!
4	!
5	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6	d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, &
7	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
8	frac_impa, frac_nucl, beta, &
9	prfl, psfl, rhcl, zqta, fraca, &
10	ztv, zpspsk, ztla, zthl, iflag_cld_th, &
11	iflag_ice_thermo)
12
13	!
14	USE dimphy
15	USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
16	USE print_control_mod, ONLY: prt_level, lunout
17	USE cloudth_mod
18	USE ioipsl_getin_p_mod, ONLY : getin_p
19	USE phys_local_var_mod, ONLY: ql_seri,qs_seri
20	USE phys_local_var_mod, ONLY: rneblsvol
21	! flag to include modifications to ensure energy conservation (if flag >0)
22	USE add_phys_tend_mod, only : fl_cor_ebil
23	IMPLICIT none
24	!======================================================================
25	! Auteur(s): Z.X. Li (LMD/CNRS)
26	! Date: le 20 mars 1995
27	! Objet: condensation et precipitation stratiforme.
28	! schema de nuage
29	! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
30	! et ilp (il pleut, L. Li)
31	! Principales parties:
32	! P0> Thermalisation des precipitations venant de la couche du dessus
33	! P1> Evaporation de la precipitation (qui vient du niveau k+1)
34	! P2> Formation du nuage (en k)
35	! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
36	! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
37	! les valeurs de T et Q initiales
38	! P2.A.2> Prise en compte du couplage entre eau condensee et T.
39	! P2.A.3> Calcul des valeures finales associees a la formation des nuages
40	! P2.B> Nuage "tout ou rien"
41	! P2.C> Prise en compte de la Chaleur latente apres formation nuage
42	! P3> Formation de la precipitation (en k)
43	!======================================================================
44	! JLD:
45	! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
46	! * fl_cor_ebil doit etre > 0 ;
47	! fl_cor_ebil= 0 pour reproduire anciens bugs
48	!======================================================================
49	include "YOMCST.h"
50	include "fisrtilp.h"
51	include "nuage.h" ! JBM (3/14)
52
53	!
54	! Principaux inputs:
55	!
56	REAL, INTENT(IN) :: dtime ! intervalle du temps (s)
57	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche
58	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression au milieu de couche
59	REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
60	REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! humidite specifique (kg/kg)
61	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! points ou le schema de conv. prof. est actif
62	INTEGER, INTENT(IN) :: iflag_cld_th
63	INTEGER, INTENT(IN) :: iflag_ice_thermo
64	!
65	! Inputs lies aux thermiques
66	!
67	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv
68	REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta, fraca
69	REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk, ztla
70	REAL, DIMENSION(klon,klev), INTENT(IN) :: zthl
71	!
72	! Input/output
73	REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! determine la largeur de distribution de vapeur
74	!
75	! Principaux outputs:
76	!
77	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! incrementation de la temperature (K)
78	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! incrementation de la vapeur d'eau
79	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! incrementation de l'eau liquide
80	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! incrementation de l'eau glace
81	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! fraction nuageuse
82	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radliq ! eau liquide utilisee dans rayonnements
83	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! humidite relative en ciel clair
84	REAL, DIMENSION(klon), INTENT(OUT) :: rain
85	REAL, DIMENSION(klon), INTENT(OUT) :: snow
86	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl
87	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl
88
89	!AA
90	! Coeffients de fraction lessivee : pour OFF-LINE
91	!
92	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_nucl
93	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_1nucl
94	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_impa
95	!
96	! Fraction d'aerosols lessivee par impaction et par nucleation
97	! POur ON-LINE
98	!
99	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa
100	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl
101	!AA
102	! --------------------------------------------------------------------------------
103	!
104	! Options du programme:
105	!
106	REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
107	!$OMP THREADPRIVATE(seuil_neb)
108
109
110	INTEGER ninter ! sous-intervals pour la precipitation
111	PARAMETER (ninter=5)
112	INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
113	!$OMP THREADPRIVATE(iflag_evap_prec)
114	!
115	LOGICAL cpartiel ! condensation partielle
116	PARAMETER (cpartiel=.TRUE.)
117	REAL t_coup
118	PARAMETER (t_coup=234.0)
119	REAL DDT0
120	PARAMETER (DDT0=.01)
121	REAL ztfondue
122	PARAMETER (ztfondue=278.15)
123	! --------------------------------------------------------------------------------
124	!
125	! Variables locales:
126	!
127	INTEGER i, k, n, kk
128	INTEGER,save::itap=0
129	!$OMP THREADPRIVATE(itap)
130
131	REAL qsl, qsi
132	real zct ,zcl
133	INTEGER ncoreczq
134	REAL ctot(klon,klev)
135	REAL ctot_vol(klon,klev)
136	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
137	REAL zdqsdT_raw(klon)
138	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
139
140	logical lognormale(klon)
141	logical ice_thermo
142	LOGICAL convergence(klon)
143	INTEGER n_i(klon), iter
144	REAL cste
145
146	real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
147	real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
148	real erf
149	REAL qcloud(klon)
150
151	REAL zrfl(klon), zrfln(klon), zqev, zqevt
152	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
153	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
154	REAL zoliqp(klon), zoliqi(klon)
155	REAL zt(klon)
156	! JBM (3/14) nexpo is replaced by exposant_glace
157	! REAL nexpo ! exponentiel pour glace/eau
158	! INTEGER, PARAMETER :: nexpo=6
159	INTEGER exposant_glace_old
160	REAL t_glace_min_old
161	REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
162	REAL zchau ,zfroi ,zfice(klon),zneb(klon),znebprecip(klon)
163	REAL zmelt, zpluie, zice
164	REAL dzfice(klon)
165	REAL zsolid
166	!!!!
167	! Variables pour Bergeron
168	REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
169	REAL zqpreci(klon), zqprecl(klon)
170	! Variable pour conservation enegie des precipitations
171	REAL zmqc(klon)
172	!
173	LOGICAL appel1er
174	SAVE appel1er
175	!$OMP THREADPRIVATE(appel1er)
176	!
177	! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
178	! iflag_oldbug_fisrtilp=1 ajoute le BUG
179	INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug
180	!$OMP THREADPRIVATE(iflag_oldbug_fisrtilp)
181	!---------------------------------------------------------------
182	!
183	!AA Variables traceurs:
184	!AA Provisoire !!! Parametres alpha du lessivage
185	!AA A priori on a 4 scavenging # possibles
186	!
187	REAL a_tr_sca(4)
188	save a_tr_sca
189	!$OMP THREADPRIVATE(a_tr_sca)
190	!
191	! Variables intermediaires
192	!
193	REAL zalpha_tr
194	REAL zfrac_lessi
195	REAL zprec_cond(klon)
196	!AA
197	! RomP >>> 15 nov 2012
198	REAL beta(klon,klev) ! taux de conversion de l'eau cond
199	! RomP <<<
200	REAL zmair(klon), zcpair, zcpeau
201	! Pour la conversion eau-neige
202	REAL zlh_solid(klon), zm_solid
203	!---------------------------------------------------------------
204	!
205	! Fonctions en ligne:
206	!
207	REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
208	! (Heymsfield & Donner, 1990)
209	REAL zzz
210
211	include "YOETHF.h"
212	include "FCTTRE.h"
213	fallvc (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_con
214	fallvs (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_lsc
215	!
216	DATA appel1er /.TRUE./
217	!ym
218	!CR: pour iflag_ice_thermo=2, on active que la convection
219	! ice_thermo = iflag_ice_thermo .GE. 1
220
221	itap=itap+1
222	znebprecip(:)=0.
223
224	ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
225	zdelq=0.0
226	ctot_vol(1:klon,1:klev)=0.0
227	rneblsvol(1:klon,1:klev)=0.0
228
229	if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
230	IF (appel1er) THEN
231	CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
232	CALL getin_p('iflag_evap_prec',iflag_evap_prec)
233	CALL getin_p('seuil_neb',seuil_neb)
234	write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp
235	!
236	WRITE(lunout,*) 'fisrtilp, ninter:', ninter
237	WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
238	WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
239
240	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
241	WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
242	WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
243	! CALL abort
244	ENDIF
245	appel1er = .FALSE.
246	!
247	!AA initialiation provisoire
248	a_tr_sca(1) = -0.5
249	a_tr_sca(2) = -0.5
250	a_tr_sca(3) = -0.5
251	a_tr_sca(4) = -0.5
252	!
253	!AA Initialisation a 1 des coefs des fractions lessivees
254	!
255	!cdir collapse
256	DO k = 1, klev
257	DO i = 1, klon
258	pfrac_nucl(i,k)=1.
259	pfrac_1nucl(i,k)=1.
260	pfrac_impa(i,k)=1.
261	beta(i,k)=0. !RomP initialisation
262	ENDDO
263	ENDDO
264
265	ENDIF ! test sur appel1er
266	!
267	!MAf Initialisation a 0 de zoliq
268	! DO i = 1, klon
269	! zoliq(i)=0.
270	! ENDDO
271	! Determiner les nuages froids par leur temperature
272	! nexpo regle la raideur de la transition eau liquide / eau glace.
273	!
274	!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
275	IF (iflag_t_glace.EQ.0) THEN
276	! ztglace = RTT - 15.0
277	t_glace_min_old = RTT - 15.0
278	!AJ<
279	IF (ice_thermo) THEN
280	! nexpo = 2
281	exposant_glace_old = 2
282	ELSE
283	! nexpo = 6
284	exposant_glace_old = 6
285	ENDIF
286
287	ENDIF
288
289	!! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
290	!! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
291	!>AJ
292	!cc nexpo = 1
293	!
294	! Initialiser les sorties:
295	!
296	!cdir collapse
297	DO k = 1, klev+1
298	DO i = 1, klon
299	prfl(i,k) = 0.0
300	psfl(i,k) = 0.0
301	ENDDO
302	ENDDO
303
304	!cdir collapse
305	DO k = 1, klev
306	DO i = 1, klon
307	d_t(i,k) = 0.0
308	d_q(i,k) = 0.0
309	d_ql(i,k) = 0.0
310	d_qi(i,k) = 0.0
311	rneb(i,k) = 0.0
312	radliq(i,k) = 0.0
313	frac_nucl(i,k) = 1.
314	frac_impa(i,k) = 1.
315	ENDDO
316	ENDDO
317	DO i = 1, klon
318	rain(i) = 0.0
319	snow(i) = 0.0
320	zoliq(i)=0.
321	! ENDDO
322	!
323	! Initialiser le flux de precipitation a zero
324	!
325	! DO i = 1, klon
326	zrfl(i) = 0.0
327	zifl(i) = 0.0
328	zneb(i) = seuil_neb
329	ENDDO
330	!
331	!
332	!AA Pour plus de securite
333
334	zalpha_tr = 0.
335	zfrac_lessi = 0.
336
337	!AA==================================================================
338	!
339	ncoreczq=0
340	! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
341	!
342	DO k = klev, 1, -1
343	!
344	!AA===============================================================
345	!
346	! Initialisation temperature et vapeur
347	DO i = 1, klon
348	zt(i)=t(i,k)
349	zq(i)=q(i,k)
350	ENDDO
351	!
352	! ----------------------------------------------------------------
353	! P0> Thermalisation des precipitations venant de la couche du dessus
354	! ----------------------------------------------------------------
355	! Calculer la varition de temp. de l'air du a la chaleur sensible
356	! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
357	! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
358	! des differentes transformations thermodynamiques. Cette masse d'eau doit
359	! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
360	! de l'enthalpie de la couche avec la temperature
361	! Variables calculees ou modifiees:
362	! - zt: temperature de la cocuhe
363	! - zmqc: masse de precip qui doit etre thermalisee
364	!
365	IF(k.LE.klevm1) THEN
366	DO i = 1, klon
367	!IM
368	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
369	! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
370	zcpair=RCPD(1.0+RVTMP2zq(i))
371	zcpeau=RCPD*RVTMP2
372	if (fl_cor_ebil .GT. 0) then
373	! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
374	! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
375	! derniere couche
376	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
377	! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
378	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zmqc(i)zcpeau + zcpair*zt(i) ) &
379	/ (zcpair + zmqc(i)*zcpeau)
380	else ! si on maintient les anciennes erreurs
381	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zrfl(i)dtime*zcpeau &
382	+ zmair(i)zcpairzt(i) ) &
383	/ (zmair(i)zcpair + zrfl(i)dtime*zcpeau)
384	end if
385	ENDDO
386	ELSE ! IF(k.LE.klevm1)
387	DO i = 1, klon
388	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
389	zmqc(i) = 0.
390	ENDDO
391	ENDIF ! end IF(k.LE.klevm1)
392	!
393	! ----------------------------------------------------------------
394	! P1> Calcul de l'evaporation de la precipitation
395	! ----------------------------------------------------------------
396	! On evapore une partie des precipitations venant de la maille du dessus.
397	! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
398	! ce que la fraction de cette couche qui est sous le nuage soit saturee.
399	! Variables calculees ou modifiees:
400	! - zrfl et zifl: flux de precip liquide et glace
401	! - zq, zt: humidite et temperature de la cocuhe
402	! - zmqc: masse de precip qui doit etre thermalisee
403	!
404	IF (iflag_evap_prec>=1) THEN
405	DO i = 1, klon
406	! S'il y a des precipitations
407	IF (zrfl(i)+zifl(i).GT.0.) THEN
408	! Calcul du qsat
409	IF (thermcep) THEN
410	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
411	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
412	zqs(i)=MIN(0.5,zqs(i))
413	zcor=1./(1.-RETV*zqs(i))
414	zqs(i)=zqs(i)*zcor
415	ELSE
416	IF (zt(i) .LT. t_coup) THEN
417	zqs(i) = qsats(zt(i)) / pplay(i,k)
418	ELSE
419	zqs(i) = qsatl(zt(i)) / pplay(i,k)
420	ENDIF
421	ENDIF
422	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
423	ENDDO
424	!AJ<
425
426	IF (.NOT. ice_thermo) THEN
427	DO i = 1, klon
428	! S'il y a des precipitations
429	IF (zrfl(i)+zifl(i).GT.0.) THEN
430	! Evap max pour ne pas saturer la fraction sous le nuage
431	! Evap max jusqu'à atteindre la saturation dans la partie
432	! de la maille qui est sous le nuage de la couche du dessus
433	!!! On ne tient compte de cette fraction que sous une seule
434	!!! couche sous le nuage
435	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
436	! Ajout de la prise en compte des precip a thermiser
437	! avec petite reecriture
438	if (fl_cor_ebil .GT. 0) then ! nouveau
439	! Calcul de l'evaporation du flux de precip herite
440	! d'au-dessus
441	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
442	* zmair(i)/pplay(i,k)zt(i)RD
443	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
444
445	! Seuil pour ne pas saturer la fraction sous le nuage
446	zqev = MIN (zqev, zqevt)
447	! Nouveau flux de precip
448	zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
449	! Aucun flux liquide pour T < t_coup, on reevapore tout.
450	IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
451	zrfln(i)=0.
452	zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
453	END IF
454	! Nouvelle vapeur
455	zq(i) = zq(i) + zqev
456	zmqc(i) = zmqc(i)-zqev
457	! Nouvelle temperature (chaleur latente)
458	zt(i) = zt(i) - zqev &
459	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
460	!!JLD debut de partie a supprimer a terme
461	else ! if (fl_cor_ebil .GT. 0)
462	! Calcul de l'evaporation du flux de precip herite
463	! d'au-dessus
464	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
465	* (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
466	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
467	* RG*dtime/(paprs(i,k)-paprs(i,k+1))
468	! Seuil pour ne pas saturer la fraction sous le nuage
469	zqev = MIN (zqev, zqevt)
470	! Nouveau flux de precip
471	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
472	/RG/dtime
473	! Aucun flux liquide pour T < t_coup
474	IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
475	! Nouvelle vapeur
476	zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
477	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
478	! Nouvelle temperature (chaleur latente)
479	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
480	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
481	* RLVTT/RCPD/(1.0+RVTMP2*zq(i))
482	end if ! if (fl_cor_ebil .GT. 0)
483	!!JLD fin de partie a supprimer a terme
484	zrfl(i) = zrfln(i)
485	zifl(i) = 0.
486	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
487	ENDDO
488	!
489	ELSE ! (.NOT. ice_thermo)
490	! ================================
491	! Avec thermodynamique de la glace
492	! ================================
493	DO i = 1, klon
494	!AJ<
495	! S'il y a des precipitations
496	IF (zrfl(i)+zifl(i).GT.0.) THEN
497
498	IF (iflag_evap_prec==1) THEN
499	znebprecip(i)=zneb(i)
500	ELSE
501	znebprecip(i)=MAX(zneb(i),znebprecip(i))
502	ENDIF
503
504	! Evap max pour ne pas saturer la fraction sous le nuage
505	zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
506
507	!JAM
508	! On differencie qsat pour l'eau et la glace
509	! Si zdelta=1. --> glace
510	! Si zdelta=0. --> eau liquide
511
512	! Calcul du qsat par rapport a l'eau liquide
513	qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
514	qsl= MIN(0.5,qsl)
515	zcor= 1./(1.-RETV*qsl)
516	qsl= qsl*zcor
517
518	! Calcul de l'evaporation du flux de precip venant du dessus
519	! Formulation en racine du flux de precip
520	! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
521	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
522	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
523	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
524	RGdtime/(paprs(i,k)-paprs(i,k+1))
525
526	! Calcul du qsat par rapport a la glace
527	qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
528	qsi= MIN(0.5,qsi)
529	zcor= 1./(1.-RETV*qsi)
530	qsi= qsi*zcor
531
532	! Calcul de la sublimation du flux de precip solide herite
533	! d'au-dessus
534	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
535	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
536	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
537	RGdtime/(paprs(i,k)-paprs(i,k+1))
538
539	!JAM
540	! Limitation de l'evaporation. On s'assure qu'on ne sature pas
541	! la fraction de la couche sous le nuage sinon on repartit zqev0
542	! en conservant la proportion liquide / glace
543
544	IF (zqevt+zqevti.GT.zqev0) THEN
545	zqev=zqev0*zqevt/(zqevt+zqevti)
546	zqevi=zqev0*zqevti/(zqevt+zqevti)
547	ELSE
548	!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
549	! liquides et solides meme si on ne sature pas la couche.
550	! A mon avis, le test est inutile, et il faudrait tout remplacer par:
551	! zqev=zqevt
552	! zqevi=zqevti
553	IF (zqevt+zqevti.GT.0.) THEN
554	zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
555	zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
556	ELSE
557	zqev=0.
558	zqevi=0.
559	ENDIF
560	ENDIF
561	! Nouveaux flux de precip liquide et solide
562	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
563	/RG/dtime)
564	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
565	/RG/dtime)
566
567	! Mise a jour de la vapeur, temperature et flux de precip
568	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
569	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
570	if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
571	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
572	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
573	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
574	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
575	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
576	+ (zifln(i)-zifl(i)) &
577	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
578	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
579	else ! sans correction thermalisation des precips
580	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
581	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
582	* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
583	+ (zifln(i)-zifl(i)) &
584	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
585	* RLSTT/RCPD/(1.0+RVTMP2*zq(i))
586	end if
587	! Nouvelles vaeleurs des precips liquides et solides
588	zrfl(i) = zrfln(i)
589	zifl(i) = zifln(i)
590	! print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
591
592	!CR ATTENTION: deplacement de la fonte de la glace
593	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
594	!!! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 !!!!!!!!! jyg
595	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
596	zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! jyg
597	! fraction de la precip solide qui est fondue
598	zmelt = MIN(MAX(zmelt,0.),1.)
599	! Fusion de la glace
600	zrfl(i)=zrfl(i)+zmelt*zifl(i)
601	if (fl_cor_ebil .LE. 0) then
602	! the following line should not be here. Indeed, if zifl is modified
603	! now, zifl(i)*zmelt is no more the amount of ice that has melt
604	! and therefore the change in temperature computed below is wrong
605	zifl(i)=zifl(i)*(1.-zmelt)
606	end if
607	! Chaleur latente de fusion
608	if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
609	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
610	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
611	else ! sans correction thermalisation des precips
612	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
613	RLMLT/RCPD/(1.0+RVTMP2zq(i))
614	end if
615	if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
616	zifl(i)=zifl(i)*(1.-zmelt)
617	end if
618
619	ELSE
620	! Si on n'a plus de pluies, on reinitialise a 0 la farcion
621	! sous nuageuse utilisee pour la pluie.
622	znebprecip(i)=0.
623	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
624	ENDDO
625
626	ENDIF ! (.NOT. ice_thermo)
627
628	! ----------------------------------------------------------------
629	! Fin evaporation de la precipitation
630	! ----------------------------------------------------------------
631	ENDIF ! (iflag_evap_prec>=1)
632	!
633	! Calculer Qs et L/Cp*dQs/dT:
634	!
635	IF (thermcep) THEN
636	DO i = 1, klon
637	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
638	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
639	if (fl_cor_ebil .GT. 0) then ! nouveau
640	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
641	else
642	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
643	endif
644	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
645	zqs(i) = MIN(0.5,zqs(i))
646	zcor = 1./(1.-RETV*zqs(i))
647	zqs(i) = zqs(i)*zcor
648	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
649	zdqsdT_raw(i) = zdqs(i)* &
650	& RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTTzdelta)
651	ENDDO
652	ELSE
653	DO i = 1, klon
654	IF (zt(i).LT.t_coup) THEN
655	zqs(i) = qsats(zt(i))/pplay(i,k)
656	zdqs(i) = dqsats(zt(i),zqs(i))
657	ELSE
658	zqs(i) = qsatl(zt(i))/pplay(i,k)
659	zdqs(i) = dqsatl(zt(i),zqs(i))
660	ENDIF
661	ENDDO
662	ENDIF
663	!
664	! Determiner la condensation partielle et calculer la quantite
665	! de l'eau condensee:
666	!
667	!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
668	! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
669	! write(,) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
670	! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
671	! stop
672	! endif
673
674	! ----------------------------------------------------------------
675	! P2> Formation du nuage
676	! ----------------------------------------------------------------
677	! Variables calculees:
678	! rneb : fraction nuageuse
679	! zcond : eau condensee moyenne dans la maille.
680	! rhcl: humidite relative ciel-clair
681	! zt : temperature de la maille
682	! ----------------------------------------------------------------
683	!
684	IF (cpartiel) THEN
685	! -------------------------
686	! P2.A> Nuage fractionnaire
687	! -------------------------
688	!
689	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
690	! nuageuse a partir des PDF de Sandrine Bony.
691	! rneb : fraction nuageuse
692	! zqn : eau totale dans le nuage
693	! zcond : eau condensee moyenne dans la maille.
694	! on prend en compte le réchauffement qui diminue la partie
695	! condensee
696	!
697	! Version avec les raqts
698
699	! ----------------------------------------------------------------
700	! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
701	! ----------------------------------------------------------------
702	if (iflag_pdf.eq.0) then
703
704	! version creneau de (Li, 1998)
705	do i=1,klon
706	zdelq = min(ratqs(i,k),0.99) * zq(i)
707	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
708	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
709	enddo
710
711	else ! if (iflag_pdf.eq.0)
712	! ----------------------------------------------------------------
713	! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
714	! les valeurs de T et Q initiales
715	! ----------------------------------------------------------------
716	do i=1,klon
717	if(zq(i).lt.1.e-15) then
718	ncoreczq=ncoreczq+1
719	zq(i)=1.e-15
720	endif
721	enddo
722
723	if (iflag_cld_th>=5) then
724
725	if (iflag_cloudth_vert<=2) then
726	call cloudth(klon,klev,k,ztv, &
727	zq,zqta,fraca, &
728	qcloud,ctot,zpspsk,paprs,ztla,zthl, &
729	ratqs,zqs,t)
730	elseif (iflag_cloudth_vert>=3) then
731	call cloudth_v3(klon,klev,k,ztv, &
732	zq,zqta,fraca, &
733	qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
734	ratqs,zqs,t)
735	endif
736	do i=1,klon
737	rneb(i,k)=ctot(i,k)
738	rneblsvol(i,k)=ctot_vol(i,k)
739	zqn(i)=qcloud(i)
740	enddo
741
742	endif
743
744	if (iflag_cld_th <= 4) then
745	lognormale = .true.
746	elseif (iflag_cld_th >= 6) then
747	! lognormale en l'absence des thermiques
748	lognormale = fraca(:,k) < 1e-10
749	else
750	! Dans le cas iflag_cld_th=5, on prend systématiquement la
751	! bi-gaussienne
752	lognormale = .false.
753	end if
754
755	!CR: variation de qsat avec T en presence de glace ou non
756	!initialisations
757	do i=1,klon
758	DT(i) = 0.
759	n_i(i)=0
760	Tbef(i)=zt(i)
761	qlbef(i)=0.
762	enddo
763
764	! ----------------------------------------------------------------
765	! P2.A.2> Prise en compte du couplage entre eau condensee et T.
766	! Calcul des grandeurs nuageuses en tenant compte de l'effet de
767	! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
768	! qui en dependent. Ce changement de temperature est provisoire, et
769	! la valeur definitive sera calcule plus tard.
770	! Variables calculees:
771	! rneb : nebulosite
772	! zcond: eau condensee en moyenne dans la maille
773	! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
774	! pas clair, il n'y a probablement pas de prise en compte de l'effet de
775	! T sur qsat
776	! ----------------------------------------------------------------
777
778	!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
779	if (iflag_fisrtilp_qsat.ge.0) then
780	! Iteration pour condensation avec variation de qsat(T)
781	! -----------------------------------------------------
782	do iter=1,iflag_fisrtilp_qsat+1
783
784	do i=1,klon
785	! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
786	! !! convergence = .true. tant que l'on n'a pas converge !!
787	! ------------------------------
788	convergence(i)=abs(DT(i)).gt.DDT0
789	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
790	! si on n'a pas converge
791	!
792	! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
793	! ---------------------------------------------------------------
794	! Variables calculees:
795	! rneb : nebulosite
796	! zqn : eau condensee, dans le nuage (in cloud water content)
797	! zcond: eau condensee en moyenne dans la maille
798	! rhcl: humidite relative ciel-clair
799	!
800	Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
801	if (.not.ice_thermo) then
802	zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
803	else
804	if (iflag_t_glace.eq.0) then
805	zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
806	else if (iflag_t_glace.ge.1) then
807	if (iflag_oldbug_fisrtilp.EQ.0) then
808	zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
809	else
810	!avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
811	zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
812	endif
813	endif
814	endif
815	! Calcul de rneb, qzn et zcond pour les PDF lognormales
816	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
817	if (fl_cor_ebil .GT. 0) then
818	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
819	else
820	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
821	end if
822	zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
823	zqs(i) = MIN(0.5,zqs(i))
824	zcor = 1./(1.-RETV*zqs(i))
825	zqs(i) = zqs(i)*zcor
826	zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
827	zpdf_sig(i)=ratqs(i,k)*zq(i)
828	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
829	zpdf_delta(i)=log(zq(i)/zqs(i))
830	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
831	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
832	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
833	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
834	zpdf_e1(i)=1.-erf(zpdf_e1(i))
835	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
836	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
837	zpdf_e2(i)=1.-erf(zpdf_e2(i))
838
839	if (zpdf_e1(i).lt.1.e-10) then
840	rneb(i,k)=0.
841	zqn(i)=zqs(i)
842	else
843	rneb(i,k)=0.5*zpdf_e1(i)
844	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
845	endif
846
847	! If vertical heterogeneity, change fraction by volume as well
848	if (iflag_cloudth_vert>=3) then
849	ctot_vol(i,k)=rneb(i,k)
850	rneblsvol(i,k)=ctot_vol(i,k)
851	endif
852
853	endif !convergence
854
855	enddo ! boucle en i
856
857	! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
858	! due a la condensation.
859	! ---------------------------------------------------------------
860	! Variables calculees:
861	! DT : variation de temperature due a la condensation
862
863	if (.not. ice_thermo) then
864	! --------------------------
865	do i=1,klon
866	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
867
868	qlbef(i)=max(0.,zqn(i)-zqs(i))
869	if (fl_cor_ebil .GT. 0) then
870	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
871	else
872	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
873	end if
874	denom=1.+rneb(i,k)*zdqs(i)
875	DT(i)=num/denom
876	n_i(i)=n_i(i)+1
877	endif
878	enddo
879
880	else ! if (.not. ice_thermo)
881	! --------------------------
882	if (iflag_t_glace.ge.1) then
883	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
884	endif
885
886	do i=1,klon
887	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
888
889	if (iflag_t_glace.eq.0) then
890	zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
891	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
892	zfice(i) = zfice(i)**exposant_glace_old
893	dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
894	& / (t_glace_min_old - RTT)
895	endif
896
897	if (iflag_t_glace.ge.1) then
898	dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
899	& / (t_glace_min - t_glace_max)
900	endif
901
902	if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
903	dzfice(i)=0.
904	endif
905
906	if (zfice(i).lt.1) then
907	cste=RLVTT
908	else
909	cste=RLSTT
910	endif
911
912	qlbef(i)=max(0.,zqn(i)-zqs(i))
913	if (fl_cor_ebil .GT. 0) then
914	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
915	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
916	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
917	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
918	& qlbef(i)dzfice(i)
919	else
920	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
921	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
922	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
923	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))rneb(i,k)qlbef(i)dzfice(i)
924	end if
925	DT(i)=num/denom
926	n_i(i)=n_i(i)+1
927
928	endif ! fin convergence
929	enddo ! fin boucle i
930
931	endif !ice_thermo
932
933	enddo ! iter=1,iflag_fisrtilp_qsat+1
934	! Fin d'iteration pour condensation avec variation de qsat(T)
935	! -----------------------------------------------------------
936	endif ! if (iflag_fisrtilp_qsat.ge.0)
937	! ----------------------------------------------------------------
938	! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
939	! ----------------------------------------------------------------
940
941	endif ! iflag_pdf
942
943	! if (iflag_fisrtilp_qsat.eq.-1) then
944	!------------------------------------------
945	!CR: ATTENTION option fausse mais a existe:
946	! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
947	! activer les lignes suivantes:
948	IF (1.eq.0) THEN
949	DO i=1,klon
950	IF (rneb(i,k) .LE. 0.0) THEN
951	zqn(i) = 0.0
952	rneb(i,k) = 0.0
953	zcond(i) = 0.0
954	rhcl(i,k)=zq(i)/zqs(i)
955	ELSE IF (rneb(i,k) .GE. 1.0) THEN
956	zqn(i) = zq(i)
957	rneb(i,k) = 1.0
958	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
959	rhcl(i,k)=1.0
960	ELSE
961	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
962	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
963	ENDIF
964	ENDDO
965	ENDIF
966	!------------------------------------------
967
968	! ELSE
969	! ----------------------------------------------------------------
970	! P2.A.3> Calcul des valeures finales associees a la formation des nuages
971	! Variables calculees:
972	! rneb : nebulosite
973	! zcond: eau condensee en moyenne dans la maille
974	! zq : eau vapeur dans la maille
975	! zt : temperature de la maille
976	! rhcl: humidite relative ciel-clair
977	! ----------------------------------------------------------------
978	!
979	! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
980	! Calcule de l'eau condensee moyenne dans la maille (zcond),
981	! et de l'humidite relative ciel-clair (rhcl)
982	DO i=1,klon
983	IF (rneb(i,k) .LE. 0.0) THEN
984	zqn(i) = 0.0
985	rneb(i,k) = 0.0
986	zcond(i) = 0.0
987	rhcl(i,k)=zq(i)/zqs(i)
988	ELSE IF (rneb(i,k) .GE. 1.0) THEN
989	zqn(i) = zq(i)
990	rneb(i,k) = 1.0
991	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
992	rhcl(i,k)=1.0
993	ELSE
994	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
995	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
996	ENDIF
997	ENDDO
998	! If vertical heterogeneity, change fraction by volume as well
999	if (iflag_cloudth_vert>=3) then
1000	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1001	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1002	endif
1003
1004	! ENDIF
1005
1006	ELSE ! de IF (cpartiel)
1007	! -------------------------
1008	! P2.B> Nuage "tout ou rien"
1009	! -------------------------
1010	! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
1011	DO i = 1, klon
1012	IF (zq(i).GT.zqs(i)) THEN
1013	rneb(i,k) = 1.0
1014	ELSE
1015	rneb(i,k) = 0.0
1016	ENDIF
1017	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
1018	ENDDO
1019	ENDIF ! de IF (cpartiel)
1020	!
1021	! Mise a jour vapeur d'eau
1022	! -------------------------
1023	DO i = 1, klon
1024	zq(i) = zq(i) - zcond(i)
1025	! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
1026	ENDDO
1027	!AJ<
1028	! ------------------------------------
1029	! P2.C> Prise en compte de la Chaleur latente apres formation nuage
1030	! -------------------------------------
1031	! Variable calcule:
1032	! zt : temperature de la maille
1033	!
1034	IF (.NOT. ice_thermo) THEN
1035	if (iflag_fisrtilp_qsat.lt.1) then
1036	DO i = 1, klon
1037	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
1038	ENDDO
1039	else if (iflag_fisrtilp_qsat.gt.0) then
1040	DO i= 1, klon
1041	if (fl_cor_ebil .GT. 0) then
1042	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1043	else
1044	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1045	end if
1046	ENDDO
1047	endif
1048	ELSE
1049	if (iflag_t_glace.ge.1) then
1050	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1051	endif
1052	if (iflag_fisrtilp_qsat.lt.1) then
1053	DO i = 1, klon
1054	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1055	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1056	! t_glace_max, exposant_glace)
1057	if (iflag_t_glace.eq.0) then
1058	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1059	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1060	zfice(i) = zfice(i)**exposant_glace_old
1061	endif
1062	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
1063	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
1064	ENDDO
1065	else
1066	DO i=1, klon
1067	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1068	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1069	! t_glace_max, exposant_glace)
1070	if (iflag_t_glace.eq.0) then
1071	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1072	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1073	zfice(i) = zfice(i)**exposant_glace_old
1074	endif
1075	if (fl_cor_ebil .GT. 0) then
1076	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1077	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1078	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1079	else
1080	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
1081	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1082	end if
1083	ENDDO
1084	endif
1085	! print*,zt(i),zrfl(i),zifl(i),'temp1'
1086	ENDIF
1087	!>AJ
1088
1089	! ----------------------------------------------------------------
1090	! P3> Formation des precipitations
1091	! ----------------------------------------------------------------
1092	!
1093	! Partager l'eau condensee en precipitation et eau liquide nuageuse
1094	!
1095
1096	! Initialisation de zoliq (eau condensee moyenne dans la maille)
1097	DO i = 1, klon
1098	IF (rneb(i,k).GT.0.0) THEN
1099	zoliq(i) = zcond(i)
1100	zrho(i) = pplay(i,k) / zt(i) / RD
1101	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
1102	ENDIF
1103	ENDDO
1104	!AJ<
1105	IF (.NOT. ice_thermo) THEN
1106	IF (iflag_t_glace.EQ.0) THEN
1107	DO i = 1, klon
1108	IF (rneb(i,k).GT.0.0) THEN
1109	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
1110	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1111	zfice(i) = zfice(i)**exposant_glace_old
1112	! zfice(i) = zfice(i)**nexpo
1113	!! zfice(i)=0.
1114	ENDIF
1115	ENDDO
1116	ELSE ! of IF (iflag_t_glace.EQ.0)
1117	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1118	! DO i = 1, klon
1119	! IF (rneb(i,k).GT.0.0) THEN
1120	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1121	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1122	! t_glace_max, exposant_glace)
1123	! ENDIF
1124	! ENDDO
1125	ENDIF
1126	ENDIF
1127
1128	! Calcul de radliq (eau condensee pour le rayonnement)
1129	! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
1130	! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
1131	! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
1132	! transmise au rayonnement;
1133	! ----------------------------------------------------------------
1134	DO i = 1, klon
1135	IF (rneb(i,k).GT.0.0) THEN
1136	zneb(i) = MAX(rneb(i,k), seuil_neb)
1137	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
1138	! print*,zt(i),'fractionglace'
1139	!>AJ
1140	radliq(i,k) = zoliq(i)/REAL(ninter+1)
1141	ENDIF
1142	ENDDO
1143	!
1144	DO n = 1, ninter
1145	DO i = 1, klon
1146	IF (rneb(i,k).GT.0.0) THEN
1147	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
1148	! Initialization of zpluie and zice:
1149	zpluie=0
1150	zice=0
1151	IF (zneb(i).EQ.seuil_neb) THEN
1152	ztot = 0.0
1153	ELSE
1154	! quantite d'eau a eliminer: zchau (Sundqvist, 1978)
1155	! meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
1156	if (ptconv(i,k)) then
1157	zcl =cld_lc_con
1158	zct =1./cld_tau_con
1159	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1160	fallvc(zrhol(i)) zfice(i)
1161	else
1162	zcl =cld_lc_lsc
1163	zct =1./cld_tau_lsc
1164	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1165	fallvs(zrhol(i)) zfice(i)
1166	endif
1167
1168	! si l'heterogeneite verticale est active, on utilise
1169	! la fraction volumique "vraie" plutot que la fraction
1170	! surfacique modifiee, qui est plus grande et reduit
1171	! sinon l'eau in-cloud de facon artificielle
1172	if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
1173	zchau = zct dtime/REAL(ninter) zoliq(i) &
1174	(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl )2)) (1.-zfice(i))
1175	else
1176	zchau = zct dtime/REAL(ninter) zoliq(i) &
1177	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
1178	endif
1179	!AJ<
1180	IF (.NOT. ice_thermo) THEN
1181	ztot = zchau + zfroi
1182	ELSE
1183	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
1184	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
1185	ztot = zpluie + zice
1186	ENDIF
1187	!>AJ
1188	ztot = MAX(ztot ,0.0)
1189	ENDIF
1190	ztot = MIN(ztot,zoliq(i))
1191	!AJ<
1192	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
1193	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
1194	!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
1195	! temporaires et ne doivent pas etre calcule (alors qu'elles le sont
1196	! si iflag_bergeron <> 2
1197	! A SUPPRIMER A TERME
1198	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
1199	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
1200	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
1201	!>AJ
1202	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
1203	ENDIF
1204	ENDDO ! i = 1,klon
1205	ENDDO ! n = 1,ninter
1206
1207	! ----------------------------------------------------------------
1208	!
1209	IF (.NOT. ice_thermo) THEN
1210	DO i = 1, klon
1211	IF (rneb(i,k).GT.0.0) THEN
1212	d_ql(i,k) = zoliq(i)
1213	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
1214	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1215	ENDIF
1216	ENDDO
1217	ELSE
1218	!
1219	!CR&JYG<
1220	! On prend en compte l'effet Bergeron dans les flux de precipitation :
1221	! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
1222	! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
1223	! et les precipitations est grossierement pris en compte en linearisant les equations
1224	! et en approximant le processus de precipitation liquide par un processus a seuil.
1225	! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
1226	! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
1227	! Le condensat precipitant solide est augmente.
1228	! La vapeur d'eau est augmentee.
1229	!
1230	IF ((iflag_bergeron .EQ. 2)) THEN
1231	DO i = 1, klon
1232	IF (rneb(i,k) .GT. 0.0) THEN
1233	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1234	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1235	if (fl_cor_ebil .GT. 0) then
1236	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
1237	coef1 = rneb(i,k)RLSTT/zcpzdqsdT_raw(i)
1238	! Calcul de DT si toute les precips liquides congelent
1239	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
1240	! On ne veut pas que T devienne superieur a la temp. de congelation.
1241	! donc que Delta > RTT-zt(i
1242	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1243	zt(i) = zt(i) + DeltaT
1244	! Eau vaporisee du fait de l'augmentation de T
1245	Deltaq = rneb(i,k)zdqsdT_raw(i)DeltaT
1246	! on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
1247	zq(i) = zq(i) + Deltaq
1248	! Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
1249	zcond(i) = max( zcond(i)- Deltaq, 0. )
1250	! precip liquide qui congele ou qui s'evapore
1251	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
1252	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1253	! bilan eau glacee
1254	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1255	else ! if (fl_cor_ebil .GT. 0)
1256	! ancien calcul
1257	zcp=RCPD(1.0+RVTMP2(zq(i)+zcond(i)))
1258	coef1 = RLMLT*zdqs(i)/RLVTT
1259	DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
1260	zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
1261	zqprecl(i) = max( zqprecl(i) - zcp/RLMLT(1.+coef1)DeltaT, 0. )
1262	zcond(i) = max( zcond(i) - zcp/RLVTTzdqs(i)DeltaT, 0. )
1263	zq(i) = zq(i) + zcp/RLVTTzdqs(i)DeltaT
1264	zt(i) = zt(i) + DeltaT
1265	end if ! if (fl_cor_ebil .GT. 0)
1266	ENDIF ! rneb(i,k) .GT. 0.0
1267	ENDDO
1268	DO i = 1, klon
1269	IF (rneb(i,k).GT.0.0) THEN
1270	d_ql(i,k) = (1-zfice(i))*zoliq(i)
1271	d_qi(i,k) = zfice(i)*zoliq(i)
1272	zrfl(i) = zrfl(i)+ zqprecl(i) &
1273	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1274	zifl(i) = zifl(i)+ zqpreci(i) &
1275	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1276	ENDIF
1277	ENDDO
1278	!!
1279	ELSE ! iflag_bergeron
1280	!>CR&JYG
1281	!!
1282	DO i = 1, klon
1283	IF (rneb(i,k).GT.0.0) THEN
1284	!CR on prend en compte la phase glace
1285	!JLD inutile car on ne passe jamais ici si .not.ice_thermo
1286	! if (.not.ice_thermo) then
1287	! d_ql(i,k) = zoliq(i)
1288	! d_qi(i,k) = 0.
1289	! else
1290	d_ql(i,k) = (1-zfice(i))*zoliq(i)
1291	d_qi(i,k) = zfice(i)*zoliq(i)
1292	! endif
1293	!AJ<
1294	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1295	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1296	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1297	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1298	! zrfl(i) = zrfl(i)+ zpluie &
1299	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
1300	! zifl(i) = zifl(i)+ zice &
1301	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
1302
1303	!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
1304	IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
1305	zsolid = zrfl(i)
1306	zifl(i) = zifl(i)+zrfl(i)
1307	zrfl(i) = 0.
1308	if (fl_cor_ebil .GT. 0) then
1309	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
1310	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
1311	else
1312	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
1313	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))
1314	end if
1315	ENDIF ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
1316	!RC
1317
1318	ENDIF ! rneb(i,k).GT.0.0
1319	ENDDO
1320
1321	ENDIF ! iflag_bergeron .EQ. 2
1322	ENDIF ! .NOT. ice_thermo
1323
1324	!CR: la fonte est faite au debut
1325	! IF (ice_thermo) THEN
1326	! DO i = 1, klon
1327	! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1328	! zmelt = MIN(MAX(zmelt,0.),1.)
1329	! zrfl(i)=zrfl(i)+zmelt*zifl(i)
1330	! zifl(i)=zifl(i)*(1.-zmelt)
1331	! print*,zt(i),'octavio1'
1332	! zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1333	! RLMLT/RCPD/(1.0+RVTMP2zq(i))
1334	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
1335	! ENDDO
1336	! ENDIF
1337
1338
1339	IF (.NOT. ice_thermo) THEN
1340	DO i = 1, klon
1341	IF (zt(i).LT.RTT) THEN
1342	psfl(i,k)=zrfl(i)
1343	ELSE
1344	prfl(i,k)=zrfl(i)
1345	ENDIF
1346	ENDDO
1347	ELSE
1348	! JAM*************************************************
1349	! Revoir partie ci-dessous: a quoi servent psfl et prfl?
1350	! *****************************************************
1351
1352	DO i = 1, klon
1353	! IF (zt(i).LT.RTT) THEN
1354	psfl(i,k)=zifl(i)
1355	! ELSE
1356	prfl(i,k)=zrfl(i)
1357	! ENDIF
1358	!>AJ
1359	ENDDO
1360	ENDIF
1361	! ----------------------------------------------------------------
1362	! Fin de formation des precipitations
1363	! ----------------------------------------------------------------
1364	!
1365	! Calculer les tendances de q et de t:
1366	!
1367	DO i = 1, klon
1368	d_q(i,k) = zq(i) - q(i,k)
1369	d_t(i,k) = zt(i) - t(i,k)
1370	ENDDO
1371	!
1372	!AA--------------- Calcul du lessivage stratiforme -------------
1373
1374	DO i = 1,klon
1375	!
1376	if(zcond(i).gt.zoliq(i)+1.e-10) then
1377	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1378	else
1379	beta(i,k) = 0.
1380	endif
1381	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1382	* (paprs(i,k)-paprs(i,k+1))/RG
1383	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1384	!AA lessivage nucleation LMD5 dans la couche elle-meme
1385	IF (iflag_t_glace.EQ.0) THEN
1386	if (t(i,k) .GE. t_glace_min_old) THEN
1387	zalpha_tr = a_tr_sca(3)
1388	else
1389	zalpha_tr = a_tr_sca(4)
1390	endif
1391	ELSE ! of IF (iflag_t_glace.EQ.0)
1392	if (t(i,k) .GE. t_glace_min) THEN
1393	zalpha_tr = a_tr_sca(3)
1394	else
1395	zalpha_tr = a_tr_sca(4)
1396	endif
1397	ENDIF
1398	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1399	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
1400	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1401	!
1402	! nucleation avec un facteur -1 au lieu de -0.5
1403	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1404	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
1405	ENDIF
1406	!
1407	ENDDO ! boucle sur i
1408	!
1409	!AA Lessivage par impaction dans les couches en-dessous
1410	DO kk = k-1, 1, -1
1411	DO i = 1, klon
1412	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1413	IF (iflag_t_glace.EQ.0) THEN
1414	if (t(i,kk) .GE. t_glace_min_old) THEN
1415	zalpha_tr = a_tr_sca(1)
1416	else
1417	zalpha_tr = a_tr_sca(2)
1418	endif
1419	ELSE ! of IF (iflag_t_glace.EQ.0)
1420	if (t(i,kk) .GE. t_glace_min) THEN
1421	zalpha_tr = a_tr_sca(1)
1422	else
1423	zalpha_tr = a_tr_sca(2)
1424	endif
1425	ENDIF
1426	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1427	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
1428	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1429	ENDIF
1430	ENDDO
1431	ENDDO
1432	!
1433	!AA===============================================================
1434	! FIN DE LA BOUCLE VERTICALE
1435	end DO
1436	!
1437	!AA==================================================================
1438	!
1439	! Pluie ou neige au sol selon la temperature de la 1ere couche
1440	!
1441	!CR: si la thermo de la glace est active, on calcule zifl directement
1442	IF (.NOT.ice_thermo) THEN
1443	DO i = 1, klon
1444	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1445	!AJ<
1446	! snow(i) = zrfl(i)
1447	snow(i) = zrfl(i)+zifl(i)
1448	!>AJ
1449	zlh_solid(i) = RLSTT-RLVTT
1450	ELSE
1451	rain(i) = zrfl(i)
1452	zlh_solid(i) = 0.
1453	ENDIF
1454	ENDDO
1455
1456	ELSE
1457	DO i = 1, klon
1458	snow(i) = zifl(i)
1459	rain(i) = zrfl(i)
1460	ENDDO
1461
1462	ENDIF
1463	!
1464	! For energy conservation : when snow is present, the solification
1465	! latent heat is considered.
1466	!CR: si thermo de la glace, neige deja prise en compte
1467	IF (.not.ice_thermo) THEN
1468	DO k = 1, klev
1469	DO i = 1, klon
1470	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
1471	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
1472	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1473	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair(i))
1474	END DO
1475	END DO
1476	ENDIF
1477	!
1478
1479	if (ncoreczq>0) then
1480	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1481	endif
1482
1483	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format