Context Navigation

fisrtilp.F90 @ 5183

Last change on this file since 5183 was 2969, checked in by fhourdin, 7 years ago
Petit bug dans iflag_t_glace>=1 vu uniquement en debug en 1D quand exposant < 1. C'est bien de l'avoir quand même dans les sources au cas où.
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: 57.2 KB

Line
1	!
2	! $Id: fisrtilp.F90 2969 2017-07-25 17:22:52Z dcugnet $
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	IF (iflag_evap_prec==3) THEN
522	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl) &
523	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
524	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
525	ELSE
526	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
527	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
528	ENDIF
529
530
531	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
532	RGdtime/(paprs(i,k)-paprs(i,k+1))
533
534	! Calcul du qsat par rapport a la glace
535	qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
536	qsi= MIN(0.5,qsi)
537	zcor= 1./(1.-RETV*qsi)
538	qsi= qsi*zcor
539
540	! Calcul de la sublimation du flux de precip solide herite
541	! d'au-dessus
542	IF (iflag_evap_prec==3) THEN
543	zqevti = znebprecip(i)coef_eva(1.0-zq(i)/qsi) &
544	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
545	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
546	ELSE
547	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
548	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
549	ENDIF
550	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
551	RGdtime/(paprs(i,k)-paprs(i,k+1))
552
553	!JAM
554	! Limitation de l'evaporation. On s'assure qu'on ne sature pas
555	! la fraction de la couche sous le nuage sinon on repartit zqev0
556	! en conservant la proportion liquide / glace
557
558	IF (zqevt+zqevti.GT.zqev0) THEN
559	zqev=zqev0*zqevt/(zqevt+zqevti)
560	zqevi=zqev0*zqevti/(zqevt+zqevti)
561	ELSE
562	!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
563	! liquides et solides meme si on ne sature pas la couche.
564	! A mon avis, le test est inutile, et il faudrait tout remplacer par:
565	! zqev=zqevt
566	! zqevi=zqevti
567	IF (zqevt+zqevti.GT.0.) THEN
568	zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
569	zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
570	ELSE
571	zqev=0.
572	zqevi=0.
573	ENDIF
574	ENDIF
575	! Nouveaux flux de precip liquide et solide
576	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
577	/RG/dtime)
578	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
579	/RG/dtime)
580
581	! Mise a jour de la vapeur, temperature et flux de precip
582	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
583	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
584	if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
585	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
586	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
587	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
588	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
589	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
590	+ (zifln(i)-zifl(i)) &
591	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
592	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
593	else ! sans correction thermalisation des precips
594	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
595	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
596	* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
597	+ (zifln(i)-zifl(i)) &
598	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
599	* RLSTT/RCPD/(1.0+RVTMP2*zq(i))
600	end if
601	! Nouvelles vaeleurs des precips liquides et solides
602	zrfl(i) = zrfln(i)
603	zifl(i) = zifln(i)
604	! print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
605
606	!CR ATTENTION: deplacement de la fonte de la glace
607	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
608	!!! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 !!!!!!!!! jyg
609	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
610	zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! jyg
611	! fraction de la precip solide qui est fondue
612	zmelt = MIN(MAX(zmelt,0.),1.)
613	! Fusion de la glace
614	zrfl(i)=zrfl(i)+zmelt*zifl(i)
615	if (fl_cor_ebil .LE. 0) then
616	! the following line should not be here. Indeed, if zifl is modified
617	! now, zifl(i)*zmelt is no more the amount of ice that has melt
618	! and therefore the change in temperature computed below is wrong
619	zifl(i)=zifl(i)*(1.-zmelt)
620	end if
621	! Chaleur latente de fusion
622	if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
623	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
624	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
625	else ! sans correction thermalisation des precips
626	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
627	RLMLT/RCPD/(1.0+RVTMP2zq(i))
628	end if
629	if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
630	zifl(i)=zifl(i)*(1.-zmelt)
631	end if
632
633	ELSE
634	! Si on n'a plus de pluies, on reinitialise a 0 la farcion
635	! sous nuageuse utilisee pour la pluie.
636	znebprecip(i)=0.
637	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
638	ENDDO
639
640	ENDIF ! (.NOT. ice_thermo)
641
642	! ----------------------------------------------------------------
643	! Fin evaporation de la precipitation
644	! ----------------------------------------------------------------
645	ENDIF ! (iflag_evap_prec>=1)
646	!
647	! Calculer Qs et L/Cp*dQs/dT:
648	!
649	IF (thermcep) THEN
650	DO i = 1, klon
651	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
652	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
653	if (fl_cor_ebil .GT. 0) then ! nouveau
654	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
655	else
656	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
657	endif
658	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
659	zqs(i) = MIN(0.5,zqs(i))
660	zcor = 1./(1.-RETV*zqs(i))
661	zqs(i) = zqs(i)*zcor
662	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
663	zdqsdT_raw(i) = zdqs(i)* &
664	& RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTTzdelta)
665	ENDDO
666	ELSE
667	DO i = 1, klon
668	IF (zt(i).LT.t_coup) THEN
669	zqs(i) = qsats(zt(i))/pplay(i,k)
670	zdqs(i) = dqsats(zt(i),zqs(i))
671	ELSE
672	zqs(i) = qsatl(zt(i))/pplay(i,k)
673	zdqs(i) = dqsatl(zt(i),zqs(i))
674	ENDIF
675	ENDDO
676	ENDIF
677	!
678	! Determiner la condensation partielle et calculer la quantite
679	! de l'eau condensee:
680	!
681	!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
682	! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
683	! write(,) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
684	! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
685	! stop
686	! endif
687
688	! ----------------------------------------------------------------
689	! P2> Formation du nuage
690	! ----------------------------------------------------------------
691	! Variables calculees:
692	! rneb : fraction nuageuse
693	! zcond : eau condensee moyenne dans la maille.
694	! rhcl: humidite relative ciel-clair
695	! zt : temperature de la maille
696	! ----------------------------------------------------------------
697	!
698	IF (cpartiel) THEN
699	! -------------------------
700	! P2.A> Nuage fractionnaire
701	! -------------------------
702	!
703	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
704	! nuageuse a partir des PDF de Sandrine Bony.
705	! rneb : fraction nuageuse
706	! zqn : eau totale dans le nuage
707	! zcond : eau condensee moyenne dans la maille.
708	! on prend en compte le réchauffement qui diminue la partie
709	! condensee
710	!
711	! Version avec les raqts
712
713	! ----------------------------------------------------------------
714	! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
715	! ----------------------------------------------------------------
716	if (iflag_pdf.eq.0) then
717
718	! version creneau de (Li, 1998)
719	do i=1,klon
720	zdelq = min(ratqs(i,k),0.99) * zq(i)
721	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
722	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
723	enddo
724
725	else ! if (iflag_pdf.eq.0)
726	! ----------------------------------------------------------------
727	! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
728	! les valeurs de T et Q initiales
729	! ----------------------------------------------------------------
730	do i=1,klon
731	if(zq(i).lt.1.e-15) then
732	ncoreczq=ncoreczq+1
733	zq(i)=1.e-15
734	endif
735	enddo
736
737	if (iflag_cld_th>=5) then
738
739	if (iflag_cloudth_vert<=2) then
740	call cloudth(klon,klev,k,ztv, &
741	zq,zqta,fraca, &
742	qcloud,ctot,zpspsk,paprs,ztla,zthl, &
743	ratqs,zqs,t)
744	elseif (iflag_cloudth_vert>=3) then
745	call cloudth_v3(klon,klev,k,ztv, &
746	zq,zqta,fraca, &
747	qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
748	ratqs,zqs,t)
749	endif
750	do i=1,klon
751	rneb(i,k)=ctot(i,k)
752	rneblsvol(i,k)=ctot_vol(i,k)
753	zqn(i)=qcloud(i)
754	enddo
755
756	endif
757
758	if (iflag_cld_th <= 4) then
759	lognormale = .true.
760	elseif (iflag_cld_th >= 6) then
761	! lognormale en l'absence des thermiques
762	lognormale = fraca(:,k) < 1e-10
763	else
764	! Dans le cas iflag_cld_th=5, on prend systématiquement la
765	! bi-gaussienne
766	lognormale = .false.
767	end if
768
769	!CR: variation de qsat avec T en presence de glace ou non
770	!initialisations
771	do i=1,klon
772	DT(i) = 0.
773	n_i(i)=0
774	Tbef(i)=zt(i)
775	qlbef(i)=0.
776	enddo
777
778	! ----------------------------------------------------------------
779	! P2.A.2> Prise en compte du couplage entre eau condensee et T.
780	! Calcul des grandeurs nuageuses en tenant compte de l'effet de
781	! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
782	! qui en dependent. Ce changement de temperature est provisoire, et
783	! la valeur definitive sera calcule plus tard.
784	! Variables calculees:
785	! rneb : nebulosite
786	! zcond: eau condensee en moyenne dans la maille
787	! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
788	! pas clair, il n'y a probablement pas de prise en compte de l'effet de
789	! T sur qsat
790	! ----------------------------------------------------------------
791
792	!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
793	if (iflag_fisrtilp_qsat.ge.0) then
794	! Iteration pour condensation avec variation de qsat(T)
795	! -----------------------------------------------------
796	do iter=1,iflag_fisrtilp_qsat+1
797
798	do i=1,klon
799	! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
800	! !! convergence = .true. tant que l'on n'a pas converge !!
801	! ------------------------------
802	convergence(i)=abs(DT(i)).gt.DDT0
803	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
804	! si on n'a pas converge
805	!
806	! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
807	! ---------------------------------------------------------------
808	! Variables calculees:
809	! rneb : nebulosite
810	! zqn : eau condensee, dans le nuage (in cloud water content)
811	! zcond: eau condensee en moyenne dans la maille
812	! rhcl: humidite relative ciel-clair
813	!
814	Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
815	if (.not.ice_thermo) then
816	zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
817	else
818	if (iflag_t_glace.eq.0) then
819	zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
820	else if (iflag_t_glace.ge.1) then
821	if (iflag_oldbug_fisrtilp.EQ.0) then
822	zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
823	else
824	!avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
825	zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
826	endif
827	endif
828	endif
829	! Calcul de rneb, qzn et zcond pour les PDF lognormales
830	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
831	if (fl_cor_ebil .GT. 0) then
832	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
833	else
834	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
835	end if
836	zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
837	zqs(i) = MIN(0.5,zqs(i))
838	zcor = 1./(1.-RETV*zqs(i))
839	zqs(i) = zqs(i)*zcor
840	zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
841	zpdf_sig(i)=ratqs(i,k)*zq(i)
842	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
843	zpdf_delta(i)=log(zq(i)/zqs(i))
844	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
845	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
846	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
847	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
848	zpdf_e1(i)=1.-erf(zpdf_e1(i))
849	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
850	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
851	zpdf_e2(i)=1.-erf(zpdf_e2(i))
852
853	if (zpdf_e1(i).lt.1.e-10) then
854	rneb(i,k)=0.
855	zqn(i)=zqs(i)
856	else
857	rneb(i,k)=0.5*zpdf_e1(i)
858	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
859	endif
860
861	! If vertical heterogeneity, change fraction by volume as well
862	if (iflag_cloudth_vert>=3) then
863	ctot_vol(i,k)=rneb(i,k)
864	rneblsvol(i,k)=ctot_vol(i,k)
865	endif
866
867	endif !convergence
868
869	enddo ! boucle en i
870
871	! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
872	! due a la condensation.
873	! ---------------------------------------------------------------
874	! Variables calculees:
875	! DT : variation de temperature due a la condensation
876
877	if (.not. ice_thermo) then
878	! --------------------------
879	do i=1,klon
880	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
881
882	qlbef(i)=max(0.,zqn(i)-zqs(i))
883	if (fl_cor_ebil .GT. 0) then
884	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
885	else
886	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
887	end if
888	denom=1.+rneb(i,k)*zdqs(i)
889	DT(i)=num/denom
890	n_i(i)=n_i(i)+1
891	endif
892	enddo
893
894	else ! if (.not. ice_thermo)
895	! --------------------------
896	if (iflag_t_glace.ge.1) then
897	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
898	endif
899
900	do i=1,klon
901	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
902
903	if (iflag_t_glace.eq.0) then
904	zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
905	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
906	zfice(i) = zfice(i)**exposant_glace_old
907	dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
908	& / (t_glace_min_old - RTT)
909	endif
910
911	if (iflag_t_glace.ge.1.and.zfice(i)>0.) then
912	dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
913	& / (t_glace_min - t_glace_max)
914	endif
915
916	if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
917	dzfice(i)=0.
918	endif
919
920	if (zfice(i).lt.1) then
921	cste=RLVTT
922	else
923	cste=RLSTT
924	endif
925
926	qlbef(i)=max(0.,zqn(i)-zqs(i))
927	if (fl_cor_ebil .GT. 0) then
928	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
929	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
930	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
931	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
932	& qlbef(i)dzfice(i)
933	else
934	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
935	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
936	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
937	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))rneb(i,k)qlbef(i)dzfice(i)
938	end if
939	DT(i)=num/denom
940	n_i(i)=n_i(i)+1
941
942	endif ! fin convergence
943	enddo ! fin boucle i
944
945	endif !ice_thermo
946
947	enddo ! iter=1,iflag_fisrtilp_qsat+1
948	! Fin d'iteration pour condensation avec variation de qsat(T)
949	! -----------------------------------------------------------
950	endif ! if (iflag_fisrtilp_qsat.ge.0)
951	! ----------------------------------------------------------------
952	! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
953	! ----------------------------------------------------------------
954
955	endif ! iflag_pdf
956
957	! if (iflag_fisrtilp_qsat.eq.-1) then
958	!------------------------------------------
959	!CR: ATTENTION option fausse mais a existe:
960	! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
961	! activer les lignes suivantes:
962	IF (1.eq.0) THEN
963	DO i=1,klon
964	IF (rneb(i,k) .LE. 0.0) THEN
965	zqn(i) = 0.0
966	rneb(i,k) = 0.0
967	zcond(i) = 0.0
968	rhcl(i,k)=zq(i)/zqs(i)
969	ELSE IF (rneb(i,k) .GE. 1.0) THEN
970	zqn(i) = zq(i)
971	rneb(i,k) = 1.0
972	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
973	rhcl(i,k)=1.0
974	ELSE
975	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
976	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
977	ENDIF
978	ENDDO
979	ENDIF
980	!------------------------------------------
981
982	! ELSE
983	! ----------------------------------------------------------------
984	! P2.A.3> Calcul des valeures finales associees a la formation des nuages
985	! Variables calculees:
986	! rneb : nebulosite
987	! zcond: eau condensee en moyenne dans la maille
988	! zq : eau vapeur dans la maille
989	! zt : temperature de la maille
990	! rhcl: humidite relative ciel-clair
991	! ----------------------------------------------------------------
992	!
993	! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
994	! Calcule de l'eau condensee moyenne dans la maille (zcond),
995	! et de l'humidite relative ciel-clair (rhcl)
996	DO i=1,klon
997	IF (rneb(i,k) .LE. 0.0) THEN
998	zqn(i) = 0.0
999	rneb(i,k) = 0.0
1000	zcond(i) = 0.0
1001	rhcl(i,k)=zq(i)/zqs(i)
1002	ELSE IF (rneb(i,k) .GE. 1.0) THEN
1003	zqn(i) = zq(i)
1004	rneb(i,k) = 1.0
1005	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1006	rhcl(i,k)=1.0
1007	ELSE
1008	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1009	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1010	ENDIF
1011	ENDDO
1012	! If vertical heterogeneity, change fraction by volume as well
1013	if (iflag_cloudth_vert>=3) then
1014	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1015	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1016	endif
1017
1018	! ENDIF
1019
1020	ELSE ! de IF (cpartiel)
1021	! -------------------------
1022	! P2.B> Nuage "tout ou rien"
1023	! -------------------------
1024	! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
1025	DO i = 1, klon
1026	IF (zq(i).GT.zqs(i)) THEN
1027	rneb(i,k) = 1.0
1028	ELSE
1029	rneb(i,k) = 0.0
1030	ENDIF
1031	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
1032	ENDDO
1033	ENDIF ! de IF (cpartiel)
1034	!
1035	! Mise a jour vapeur d'eau
1036	! -------------------------
1037	DO i = 1, klon
1038	zq(i) = zq(i) - zcond(i)
1039	! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
1040	ENDDO
1041	!AJ<
1042	! ------------------------------------
1043	! P2.C> Prise en compte de la Chaleur latente apres formation nuage
1044	! -------------------------------------
1045	! Variable calcule:
1046	! zt : temperature de la maille
1047	!
1048	IF (.NOT. ice_thermo) THEN
1049	if (iflag_fisrtilp_qsat.lt.1) then
1050	DO i = 1, klon
1051	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
1052	ENDDO
1053	else if (iflag_fisrtilp_qsat.gt.0) then
1054	DO i= 1, klon
1055	if (fl_cor_ebil .GT. 0) then
1056	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1057	else
1058	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1059	end if
1060	ENDDO
1061	endif
1062	ELSE
1063	if (iflag_t_glace.ge.1) then
1064	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1065	endif
1066	if (iflag_fisrtilp_qsat.lt.1) then
1067	DO i = 1, klon
1068	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1069	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1070	! t_glace_max, exposant_glace)
1071	if (iflag_t_glace.eq.0) then
1072	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1073	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1074	zfice(i) = zfice(i)**exposant_glace_old
1075	endif
1076	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
1077	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
1078	ENDDO
1079	else
1080	DO i=1, klon
1081	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1082	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1083	! t_glace_max, exposant_glace)
1084	if (iflag_t_glace.eq.0) then
1085	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1086	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1087	zfice(i) = zfice(i)**exposant_glace_old
1088	endif
1089	if (fl_cor_ebil .GT. 0) then
1090	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1091	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1092	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1093	else
1094	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
1095	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1096	end if
1097	ENDDO
1098	endif
1099	! print*,zt(i),zrfl(i),zifl(i),'temp1'
1100	ENDIF
1101	!>AJ
1102
1103	! ----------------------------------------------------------------
1104	! P3> Formation des precipitations
1105	! ----------------------------------------------------------------
1106	!
1107	! Partager l'eau condensee en precipitation et eau liquide nuageuse
1108	!
1109
1110	! Initialisation de zoliq (eau condensee moyenne dans la maille)
1111	DO i = 1, klon
1112	IF (rneb(i,k).GT.0.0) THEN
1113	zoliq(i) = zcond(i)
1114	zrho(i) = pplay(i,k) / zt(i) / RD
1115	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
1116	ENDIF
1117	ENDDO
1118	!AJ<
1119	IF (.NOT. ice_thermo) THEN
1120	IF (iflag_t_glace.EQ.0) THEN
1121	DO i = 1, klon
1122	IF (rneb(i,k).GT.0.0) THEN
1123	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
1124	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1125	zfice(i) = zfice(i)**exposant_glace_old
1126	! zfice(i) = zfice(i)**nexpo
1127	!! zfice(i)=0.
1128	ENDIF
1129	ENDDO
1130	ELSE ! of IF (iflag_t_glace.EQ.0)
1131	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1132	! DO i = 1, klon
1133	! IF (rneb(i,k).GT.0.0) THEN
1134	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1135	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1136	! t_glace_max, exposant_glace)
1137	! ENDIF
1138	! ENDDO
1139	ENDIF
1140	ENDIF
1141
1142	! Calcul de radliq (eau condensee pour le rayonnement)
1143	! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
1144	! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
1145	! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
1146	! transmise au rayonnement;
1147	! ----------------------------------------------------------------
1148	DO i = 1, klon
1149	IF (rneb(i,k).GT.0.0) THEN
1150	zneb(i) = MAX(rneb(i,k), seuil_neb)
1151	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
1152	! print*,zt(i),'fractionglace'
1153	!>AJ
1154	radliq(i,k) = zoliq(i)/REAL(ninter+1)
1155	ENDIF
1156	ENDDO
1157	!
1158	DO n = 1, ninter
1159	DO i = 1, klon
1160	IF (rneb(i,k).GT.0.0) THEN
1161	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
1162	! Initialization of zpluie and zice:
1163	zpluie=0
1164	zice=0
1165	IF (zneb(i).EQ.seuil_neb) THEN
1166	ztot = 0.0
1167	ELSE
1168	! quantite d'eau a eliminer: zchau (Sundqvist, 1978)
1169	! meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
1170	if (ptconv(i,k)) then
1171	zcl =cld_lc_con
1172	zct =1./cld_tau_con
1173	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1174	fallvc(zrhol(i)) zfice(i)
1175	else
1176	zcl =cld_lc_lsc
1177	zct =1./cld_tau_lsc
1178	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1179	fallvs(zrhol(i)) zfice(i)
1180	endif
1181
1182	! si l'heterogeneite verticale est active, on utilise
1183	! la fraction volumique "vraie" plutot que la fraction
1184	! surfacique modifiee, qui est plus grande et reduit
1185	! sinon l'eau in-cloud de facon artificielle
1186	if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
1187	zchau = zct dtime/REAL(ninter) zoliq(i) &
1188	(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl )2)) (1.-zfice(i))
1189	else
1190	zchau = zct dtime/REAL(ninter) zoliq(i) &
1191	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
1192	endif
1193	!AJ<
1194	IF (.NOT. ice_thermo) THEN
1195	ztot = zchau + zfroi
1196	ELSE
1197	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
1198	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
1199	ztot = zpluie + zice
1200	ENDIF
1201	!>AJ
1202	ztot = MAX(ztot ,0.0)
1203	ENDIF
1204	ztot = MIN(ztot,zoliq(i))
1205	!AJ<
1206	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
1207	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
1208	!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
1209	! temporaires et ne doivent pas etre calcule (alors qu'elles le sont
1210	! si iflag_bergeron <> 2
1211	! A SUPPRIMER A TERME
1212	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
1213	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
1214	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
1215	!>AJ
1216	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
1217	ENDIF
1218	ENDDO ! i = 1,klon
1219	ENDDO ! n = 1,ninter
1220
1221	! ----------------------------------------------------------------
1222	!
1223	IF (.NOT. ice_thermo) THEN
1224	DO i = 1, klon
1225	IF (rneb(i,k).GT.0.0) THEN
1226	d_ql(i,k) = zoliq(i)
1227	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
1228	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1229	ENDIF
1230	ENDDO
1231	ELSE
1232	!
1233	!CR&JYG<
1234	! On prend en compte l'effet Bergeron dans les flux de precipitation :
1235	! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
1236	! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
1237	! et les precipitations est grossierement pris en compte en linearisant les equations
1238	! et en approximant le processus de precipitation liquide par un processus a seuil.
1239	! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
1240	! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
1241	! Le condensat precipitant solide est augmente.
1242	! La vapeur d'eau est augmentee.
1243	!
1244	IF ((iflag_bergeron .EQ. 2)) THEN
1245	DO i = 1, klon
1246	IF (rneb(i,k) .GT. 0.0) THEN
1247	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1248	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1249	if (fl_cor_ebil .GT. 0) then
1250	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
1251	coef1 = rneb(i,k)RLSTT/zcpzdqsdT_raw(i)
1252	! Calcul de DT si toute les precips liquides congelent
1253	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
1254	! On ne veut pas que T devienne superieur a la temp. de congelation.
1255	! donc que Delta > RTT-zt(i
1256	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1257	zt(i) = zt(i) + DeltaT
1258	! Eau vaporisee du fait de l'augmentation de T
1259	Deltaq = rneb(i,k)zdqsdT_raw(i)DeltaT
1260	! on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
1261	zq(i) = zq(i) + Deltaq
1262	! Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
1263	zcond(i) = max( zcond(i)- Deltaq, 0. )
1264	! precip liquide qui congele ou qui s'evapore
1265	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
1266	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1267	! bilan eau glacee
1268	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1269	else ! if (fl_cor_ebil .GT. 0)
1270	! ancien calcul
1271	zcp=RCPD(1.0+RVTMP2(zq(i)+zcond(i)))
1272	coef1 = RLMLT*zdqs(i)/RLVTT
1273	DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
1274	zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
1275	zqprecl(i) = max( zqprecl(i) - zcp/RLMLT(1.+coef1)DeltaT, 0. )
1276	zcond(i) = max( zcond(i) - zcp/RLVTTzdqs(i)DeltaT, 0. )
1277	zq(i) = zq(i) + zcp/RLVTTzdqs(i)DeltaT
1278	zt(i) = zt(i) + DeltaT
1279	end if ! if (fl_cor_ebil .GT. 0)
1280	ENDIF ! rneb(i,k) .GT. 0.0
1281	ENDDO
1282	DO i = 1, klon
1283	IF (rneb(i,k).GT.0.0) THEN
1284	d_ql(i,k) = (1-zfice(i))*zoliq(i)
1285	d_qi(i,k) = zfice(i)*zoliq(i)
1286	zrfl(i) = zrfl(i)+ zqprecl(i) &
1287	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1288	zifl(i) = zifl(i)+ zqpreci(i) &
1289	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1290	ENDIF
1291	ENDDO
1292	!!
1293	ELSE ! iflag_bergeron
1294	!>CR&JYG
1295	!!
1296	DO i = 1, klon
1297	IF (rneb(i,k).GT.0.0) THEN
1298	!CR on prend en compte la phase glace
1299	!JLD inutile car on ne passe jamais ici si .not.ice_thermo
1300	! if (.not.ice_thermo) then
1301	! d_ql(i,k) = zoliq(i)
1302	! d_qi(i,k) = 0.
1303	! else
1304	d_ql(i,k) = (1-zfice(i))*zoliq(i)
1305	d_qi(i,k) = zfice(i)*zoliq(i)
1306	! endif
1307	!AJ<
1308	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1309	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1310	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1311	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1312	! zrfl(i) = zrfl(i)+ zpluie &
1313	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
1314	! zifl(i) = zifl(i)+ zice &
1315	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
1316
1317	!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
1318	IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
1319	zsolid = zrfl(i)
1320	zifl(i) = zifl(i)+zrfl(i)
1321	zrfl(i) = 0.
1322	if (fl_cor_ebil .GT. 0) then
1323	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
1324	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
1325	else
1326	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
1327	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))
1328	end if
1329	ENDIF ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
1330	!RC
1331
1332	ENDIF ! rneb(i,k).GT.0.0
1333	ENDDO
1334
1335	ENDIF ! iflag_bergeron .EQ. 2
1336	ENDIF ! .NOT. ice_thermo
1337
1338	!CR: la fonte est faite au debut
1339	! IF (ice_thermo) THEN
1340	! DO i = 1, klon
1341	! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1342	! zmelt = MIN(MAX(zmelt,0.),1.)
1343	! zrfl(i)=zrfl(i)+zmelt*zifl(i)
1344	! zifl(i)=zifl(i)*(1.-zmelt)
1345	! print*,zt(i),'octavio1'
1346	! zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1347	! RLMLT/RCPD/(1.0+RVTMP2zq(i))
1348	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
1349	! ENDDO
1350	! ENDIF
1351
1352
1353	IF (.NOT. ice_thermo) THEN
1354	DO i = 1, klon
1355	IF (zt(i).LT.RTT) THEN
1356	psfl(i,k)=zrfl(i)
1357	ELSE
1358	prfl(i,k)=zrfl(i)
1359	ENDIF
1360	ENDDO
1361	ELSE
1362	! JAM*************************************************
1363	! Revoir partie ci-dessous: a quoi servent psfl et prfl?
1364	! *****************************************************
1365
1366	DO i = 1, klon
1367	! IF (zt(i).LT.RTT) THEN
1368	psfl(i,k)=zifl(i)
1369	! ELSE
1370	prfl(i,k)=zrfl(i)
1371	! ENDIF
1372	!>AJ
1373	ENDDO
1374	ENDIF
1375	! ----------------------------------------------------------------
1376	! Fin de formation des precipitations
1377	! ----------------------------------------------------------------
1378	!
1379	! Calculer les tendances de q et de t:
1380	!
1381	DO i = 1, klon
1382	d_q(i,k) = zq(i) - q(i,k)
1383	d_t(i,k) = zt(i) - t(i,k)
1384	ENDDO
1385	!
1386	!AA--------------- Calcul du lessivage stratiforme -------------
1387
1388	DO i = 1,klon
1389	!
1390	if(zcond(i).gt.zoliq(i)+1.e-10) then
1391	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1392	else
1393	beta(i,k) = 0.
1394	endif
1395	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1396	* (paprs(i,k)-paprs(i,k+1))/RG
1397	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1398	!AA lessivage nucleation LMD5 dans la couche elle-meme
1399	IF (iflag_t_glace.EQ.0) THEN
1400	if (t(i,k) .GE. t_glace_min_old) THEN
1401	zalpha_tr = a_tr_sca(3)
1402	else
1403	zalpha_tr = a_tr_sca(4)
1404	endif
1405	ELSE ! of IF (iflag_t_glace.EQ.0)
1406	if (t(i,k) .GE. t_glace_min) THEN
1407	zalpha_tr = a_tr_sca(3)
1408	else
1409	zalpha_tr = a_tr_sca(4)
1410	endif
1411	ENDIF
1412	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1413	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
1414	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1415	!
1416	! nucleation avec un facteur -1 au lieu de -0.5
1417	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1418	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
1419	ENDIF
1420	!
1421	ENDDO ! boucle sur i
1422	!
1423	!AA Lessivage par impaction dans les couches en-dessous
1424	DO kk = k-1, 1, -1
1425	DO i = 1, klon
1426	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1427	IF (iflag_t_glace.EQ.0) THEN
1428	if (t(i,kk) .GE. t_glace_min_old) THEN
1429	zalpha_tr = a_tr_sca(1)
1430	else
1431	zalpha_tr = a_tr_sca(2)
1432	endif
1433	ELSE ! of IF (iflag_t_glace.EQ.0)
1434	if (t(i,kk) .GE. t_glace_min) THEN
1435	zalpha_tr = a_tr_sca(1)
1436	else
1437	zalpha_tr = a_tr_sca(2)
1438	endif
1439	ENDIF
1440	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1441	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
1442	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1443	ENDIF
1444	ENDDO
1445	ENDDO
1446	!
1447	!AA===============================================================
1448	! FIN DE LA BOUCLE VERTICALE
1449	end DO
1450	!
1451	!AA==================================================================
1452	!
1453	! Pluie ou neige au sol selon la temperature de la 1ere couche
1454	!
1455	!CR: si la thermo de la glace est active, on calcule zifl directement
1456	IF (.NOT.ice_thermo) THEN
1457	DO i = 1, klon
1458	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1459	!AJ<
1460	! snow(i) = zrfl(i)
1461	snow(i) = zrfl(i)+zifl(i)
1462	!>AJ
1463	zlh_solid(i) = RLSTT-RLVTT
1464	ELSE
1465	rain(i) = zrfl(i)
1466	zlh_solid(i) = 0.
1467	ENDIF
1468	ENDDO
1469
1470	ELSE
1471	DO i = 1, klon
1472	snow(i) = zifl(i)
1473	rain(i) = zrfl(i)
1474	ENDDO
1475
1476	ENDIF
1477	!
1478	! For energy conservation : when snow is present, the solification
1479	! latent heat is considered.
1480	!CR: si thermo de la glace, neige deja prise en compte
1481	IF (.not.ice_thermo) THEN
1482	DO k = 1, klev
1483	DO i = 1, klon
1484	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
1485	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
1486	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1487	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair(i))
1488	END DO
1489	END DO
1490	ENDIF
1491	!
1492
1493	if (ncoreczq>0) then
1494	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1495	endif
1496
1497	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format