Context Navigation

fisrtilp.F90 @ 2878

Last change on this file since 2878 was 2878, checked in by fhourdin, 7 years ago
Initialisation manquante de zmqc dans fisrtilp.F90. Resoud des plantages et sans doute des problÃmes de reproductibilitÃ©.
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.7 KB

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

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

Download in other formats:

Original Format