Context Navigation

lmdz_lscp_old.F90 @ 5218

Last change on this file since 5218 was 5208, checked in by Laurent Fairhead, 8 weeks ago
Louis d'Alencon modifications: Il s'agit des modifs permettant de choisir de calculer la variance dans les thermiques par le nouveau modèle (iflag_ratqs = 10) ou de continuer à calculer cette variance par la param d'Arnaud ce qui fait une version hybride avec variance pronostique dans l'environnement mais variance diagnostique dans les thermiques (iflag_ratqs = 11). Le iflag_ratqs = 4 version standard continue de faire toutes les variances de façon diagnostique.
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: 64.6 KB

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

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

Download in other formats:

Original Format