Context Navigation

fisrtilp.F90 @ 3503

Last change on this file since 3503 was 3493, checked in by idelkadi, 6 years ago
Amelioration de la representation des nuages bas (Jean Jouhaud) : Calcul des nouveaux ecarts types sigmas de la bi-Gaussienne en prenant en compte l'epaisseur des mailles, CF_vol avec ces sigmas, et CF_surf avec la methode de Brooks
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: 57.6 KB

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

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

Download in other formats:

Original Format