Context Navigation

fisrtilp.F90 @ 2873

Last change on this file since 2873 was 2839, checked in by Laurent Fairhead, 7 years ago
Merged trunk changes r2785:2838 into testing branch
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: 54.6 KB

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

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

Download in other formats:

Original Format