Context Navigation

fisrtilp.F90 @ 2791

Last change on this file since 2791 was 2720, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2664:2719 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: 41.7 KB

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

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

Download in other formats:

Original Format