Context Navigation

fisrtilp.F90 @ 4674

Last change on this file since 4674 was 2886, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2865:2885 into testing branch
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: 54.8 KB

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

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

Download in other formats:

Original Format