Context Navigation

lmdz_lscp_old.f90 @ 5480

Last change on this file since 5480 was 5480, checked in by yann meurdesoif, 14 hours ago

For GPU porting :

add wrapper for source to source tools
Separate initialization phase of lscp_old (firstilp), SAVE variable in compute subroutine cannot be managed properly.

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

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

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

Download in other formats:

Original Format