Context Navigation

fisrtilp.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

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

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

Download in other formats:

Original Format