Context Navigation

fisrtilp.F90 @ 2160

Last change on this file since 2160 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
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: 36.0 KB

Line
1	!
2	! $Id: fisrtilp.F90 2160 2014-11-28 15:36:29Z 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_cldcon, &
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_cldcon
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_cldcon>=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_cldcon <= 4) then
562	lognormale = .true.
563	elseif (iflag_cldcon >= 6) then
564	! lognormale en l'absence des thermiques
565	lognormale = fraca(:,k) < 1e-10
566	else
567	! Dans le cas iflag_cldcon=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	if (lognormale(i)) then
786	zt(i)=Tbef(i)
787	else
788	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
789	endif
790	ENDDO
791	endif
792	ELSE
793	if (iflag_t_glace.eq.1) then
794	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
795	endif
796	if (iflag_fisrtilp_qsat.lt.1) then
797	DO i = 1, klon
798	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
799	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
800	! t_glace_max, exposant_glace)
801	if (iflag_t_glace.eq.0) then
802	zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
803	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
804	zfice(i) = zfice(i)**exposant_glace_old
805	endif
806	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
807	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
808	ENDDO
809	else
810	DO i=1, klon
811	if (lognormale(i)) then
812	zt(i)=Tbef(i)
813	else
814	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
815	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
816	! t_glace_max, exposant_glace)
817	if (iflag_t_glace.eq.0) then
818	zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
819	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
820	zfice(i) = zfice(i)**exposant_glace_old
821	endif
822	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
823	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
824	endif
825	ENDDO
826	endif
827	! print*,zt(i),zrfl(i),zifl(i),'temp1'
828	ENDIF
829	!>AJ
830	!
831	! Partager l'eau condensee en precipitation et eau liquide nuageuse
832	!
833	DO i = 1, klon
834	IF (rneb(i,k).GT.0.0) THEN
835	zoliq(i) = zcond(i)
836	zrho(i) = pplay(i,k) / zt(i) / RD
837	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
838	ENDIF
839	ENDDO
840	!AJ<
841	IF (.NOT. ice_thermo) THEN
842	IF (iflag_t_glace.EQ.0) THEN
843	DO i = 1, klon
844	IF (rneb(i,k).GT.0.0) THEN
845	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
846	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
847	zfice(i) = zfice(i)**exposant_glace_old
848	! zfice(i) = zfice(i)**nexpo
849	!! zfice(i)=0.
850	ENDIF
851	ENDDO
852	ELSE ! of IF (iflag_t_glace.EQ.0)
853	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
854	! DO i = 1, klon
855	! IF (rneb(i,k).GT.0.0) THEN
856	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
857	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
858	! t_glace_max, exposant_glace)
859	! ENDIF
860	! ENDDO
861	ENDIF
862	ENDIF
863	DO i = 1, klon
864	IF (rneb(i,k).GT.0.0) THEN
865	zneb(i) = MAX(rneb(i,k), seuil_neb)
866	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
867	! print*,zt(i),'fractionglace'
868	!>AJ
869	radliq(i,k) = zoliq(i)/REAL(ninter+1)
870	ENDIF
871	ENDDO
872	!
873	DO n = 1, ninter
874	DO i = 1, klon
875	IF (rneb(i,k).GT.0.0) THEN
876	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
877	! Initialization of zpluie and zice:
878	zpluie=0
879	zice=0
880	IF (zneb(i).EQ.seuil_neb) THEN
881	ztot = 0.0
882	ELSE
883	! quantite d'eau a eliminer: zchau
884	! meme chose pour la glace: zfroi
885	if (ptconv(i,k)) then
886	zcl =cld_lc_con
887	zct =1./cld_tau_con
888	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
889	fallvc(zrhol(i)) zfice(i)
890	else
891	zcl =cld_lc_lsc
892	zct =1./cld_tau_lsc
893	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
894	fallvs(zrhol(i)) zfice(i)
895	endif
896	zchau = zct dtime/REAL(ninter) zoliq(i) &
897	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
898	!AJ<
899	IF (.NOT. ice_thermo) THEN
900	ztot = zchau + zfroi
901	ELSE
902	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
903	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
904	ztot = zpluie + zice
905	ENDIF
906	!>AJ
907	ztot = MAX(ztot ,0.0)
908	ENDIF
909	ztot = MIN(ztot,zoliq(i))
910	!AJ<
911	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
912	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
913	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
914	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
915	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
916	!>AJ
917	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
918	ENDIF
919	ENDDO
920	ENDDO
921	!
922	IF (.NOT. ice_thermo) THEN
923	DO i = 1, klon
924	IF (rneb(i,k).GT.0.0) THEN
925	d_ql(i,k) = zoliq(i)
926	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
927	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
928	ENDIF
929	ENDDO
930	ELSE
931	DO i = 1, klon
932	IF (rneb(i,k).GT.0.0) THEN
933	!CR on prend en compte la phase glace
934	if (.not.ice_thermo) then
935	d_ql(i,k) = zoliq(i)
936	d_qi(i,k) = 0.
937	else
938	d_ql(i,k) = (1-zfice(i))*zoliq(i)
939	d_qi(i,k) = zfice(i)*zoliq(i)
940	endif
941	!AJ<
942	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
943	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
944	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
945	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
946	! zrfl(i) = zrfl(i)+ zpluie &
947	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
948	! zifl(i) = zifl(i)+ zice &
949	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
950
951	ENDIF
952	ENDDO
953	ENDIF
954
955	!CR: la fonte est faite au debut
956	! IF (ice_thermo) THEN
957	! DO i = 1, klon
958	! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
959	! zmelt = MIN(MAX(zmelt,0.),1.)
960	! zrfl(i)=zrfl(i)+zmelt*zifl(i)
961	! zifl(i)=zifl(i)*(1.-zmelt)
962	! print*,zt(i),'octavio1'
963	! zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
964	! RLMLT/RCPD/(1.0+RVTMP2zq(i))
965	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
966	! ENDDO
967	! ENDIF
968
969
970	IF (.NOT. ice_thermo) THEN
971	DO i = 1, klon
972	IF (zt(i).LT.RTT) THEN
973	psfl(i,k)=zrfl(i)
974	ELSE
975	prfl(i,k)=zrfl(i)
976	ENDIF
977	ENDDO
978	ELSE
979	! JAM*************************************************
980	! Revoir partie ci-dessous: à quoi servent psfl et prfl?
981	! *****************************************************
982
983	DO i = 1, klon
984	! IF (zt(i).LT.RTT) THEN
985	psfl(i,k)=zifl(i)
986	! ELSE
987	prfl(i,k)=zrfl(i)
988	! ENDIF
989	!>AJ
990	ENDDO
991	ENDIF
992	!
993	!
994	! Calculer les tendances de q et de t:
995	!
996	DO i = 1, klon
997	d_q(i,k) = zq(i) - q(i,k)
998	d_t(i,k) = zt(i) - t(i,k)
999	ENDDO
1000	!
1001	!AA--------------- Calcul du lessivage stratiforme -------------
1002
1003	DO i = 1,klon
1004	!
1005	if(zcond(i).gt.zoliq(i)+1.e-10) then
1006	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1007	else
1008	beta(i,k) = 0.
1009	endif
1010	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1011	* (paprs(i,k)-paprs(i,k+1))/RG
1012	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1013	!AA lessivage nucleation LMD5 dans la couche elle-meme
1014	IF (iflag_t_glace.EQ.0) THEN
1015	if (t(i,k) .GE. t_glace_min_old) THEN
1016	zalpha_tr = a_tr_sca(3)
1017	else
1018	zalpha_tr = a_tr_sca(4)
1019	endif
1020	ELSE ! of IF (iflag_t_glace.EQ.0)
1021	if (t(i,k) .GE. t_glace_min) THEN
1022	zalpha_tr = a_tr_sca(3)
1023	else
1024	zalpha_tr = a_tr_sca(4)
1025	endif
1026	ENDIF
1027	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1028	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
1029	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1030	!
1031	! nucleation avec un facteur -1 au lieu de -0.5
1032	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1033	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
1034	ENDIF
1035	!
1036	ENDDO ! boucle sur i
1037	!
1038	!AA Lessivage par impaction dans les couches en-dessous
1039	DO kk = k-1, 1, -1
1040	DO i = 1, klon
1041	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1042	IF (iflag_t_glace.EQ.0) THEN
1043	if (t(i,kk) .GE. t_glace_min_old) THEN
1044	zalpha_tr = a_tr_sca(1)
1045	else
1046	zalpha_tr = a_tr_sca(2)
1047	endif
1048	ELSE ! of IF (iflag_t_glace.EQ.0)
1049	if (t(i,kk) .GE. t_glace_min) THEN
1050	zalpha_tr = a_tr_sca(1)
1051	else
1052	zalpha_tr = a_tr_sca(2)
1053	endif
1054	ENDIF
1055	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1056	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
1057	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1058	ENDIF
1059	ENDDO
1060	ENDDO
1061	!
1062	!AA----------------------------------------------------------
1063	! FIN DE BOUCLE SUR K
1064	end DO
1065	!
1066	!AA-----------------------------------------------------------
1067	!
1068	! Pluie ou neige au sol selon la temperature de la 1ere couche
1069	!
1070	!CR: si la thermo de la glace est active, on calcule zifl directement
1071	IF (.NOT.ice_thermo) THEN
1072	DO i = 1, klon
1073	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1074	!AJ<
1075	! snow(i) = zrfl(i)
1076	snow(i) = zrfl(i)+zifl(i)
1077	!>AJ
1078	zlh_solid(i) = RLSTT-RLVTT
1079	ELSE
1080	rain(i) = zrfl(i)
1081	zlh_solid(i) = 0.
1082	ENDIF
1083	ENDDO
1084
1085	ELSE
1086	DO i = 1, klon
1087	snow(i) = zifl(i)
1088	rain(i) = zrfl(i)
1089	ENDDO
1090
1091	ENDIF
1092	!
1093	! For energy conservation : when snow is present, the solification
1094	! latent heat is considered.
1095	!CR: si thermo de la glace, neige deja prise en compte
1096	IF (.not.ice_thermo) THEN
1097	DO k = 1, klev
1098	DO i = 1, klon
1099	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
1100	zmair=(paprs(i,k)-paprs(i,k+1))/RG
1101	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1102	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair)
1103	END DO
1104	END DO
1105	ENDIF
1106	!
1107
1108	if (ncoreczq>0) then
1109	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1110	endif
1111
1112	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format