Context Navigation

fisrtilp.F90 @ 2699

Last change on this file since 2699 was 2696, checked in by fhourdin, 8 years ago
Correction pour option iflag_cloudth_vert=1
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: 41.1 KB

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

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

Download in other formats:

Original Format