Context Navigation

lmdz_lscp_old.f90 @ 5300

Last change on this file since 5300 was 5285, checked in by abarral, 5 weeks ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
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: 64.6 KB

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

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

Download in other formats:

Original Format