Context Navigation

lmdz_lscp_old.F90 @ 4664

Last change on this file since 4664 was 4664, checked in by fhourdin, 10 months ago

standardisatio des noms pour lscp et fisrtilp

fisrtilp passe dans le module lmdz_lscp_old.F90
Prepartation de la replaysation de fisrtilp (deja fait pour lscp)

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

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

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

Download in other formats:

Original Format