Context Navigation

fisrtilp.F90 @ 4046

Last change on this file since 4046 was 3992, checked in by musat, 3 years ago
Correction de bug IM
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.0 KB

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

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

Download in other formats:

Original Format