Context Navigation

fisrtilp.F90 @ 4660

Last change on this file since 4660 was 4651, checked in by fhourdin, 15 months ago

Pre-replayisation pour cloudth

Cration de lmdz_cloudth_ini
cloudth_mod.F90 -> lmdz_cloudth.F90
Supression des constantes dans lmdz_cloudth
En fait, cloudth lui même n'est pas vraiment replayisable car appeler ligne par ligne.

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.8 KB

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

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

Download in other formats:

Original Format