Context Navigation

fisrtilp.F90 @ 2368

Last change on this file since 2368 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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: 35.8 KB

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

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

Download in other formats:

Original Format