Context Navigation

fisrtilp.F90 @ 4999

Last change on this file since 4999 was 4669, checked in by Laurent Fairhead, 15 months ago
Merged with trunk revision 4586 corresponding to june 2023 testing
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: 65.4 KB

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

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

Download in other formats:

Original Format