Context Navigation

fisrtilp.F90 @ 2684

Last change on this file since 2684 was 2664, checked in by fhourdin, 8 years ago
Gros bug petite modif "min" -> "max"
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: 40.8 KB

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

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

Download in other formats:

Original Format