Context Navigation

fisrtilp.F90 @ 3809

Last change on this file since 3809 was 3809, checked in by ymipsl, 10 years ago
Add LMDZ in aquaplanet configuration YM
File size: 35.9 KB

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

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

Download in other formats:

Original Format