Context Navigation

fisrtilp.F90 @ 2923

Last change on this file since 2923 was 2923, checked in by fhourdin, 7 years ago

Correction sous drapeaux d'une errance sur la reevaporation de la pluie.
flag_evap_prec = 0 : pas de reevaporation

1 : comme avant
2 : nouveau ou on utilise une surface max des nuages au

dessus comme fraction de la maille dans laquelle la
pluie peut se reevaporer (alors que c'etait la fraction
de la maille juste au dessus avant).

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

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

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

Download in other formats:

Original Format