Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

lmdz_lscp_old.f90

Last change on this file was 5765, checked in by rkazeroni, 5 weeks ago

For GPU porting:

Add one Fortran comment for each ported param to specify possible names for horizontal variables used in the param
Add Fortran comment to exclude a file and its routines from GPU porting

Note:

The "gpum horizontal" comment is placed in the param main. It also applies to dependencies of the main.
The "gpum nogpu" excludes all modules and routines of the file from GPU porting. This is used for routines not targeted for GPU computing.

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

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

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

Download in other formats:

Original Format