Context Navigation

fisrtilp.F90 @ 2865

Last change on this file since 2865 was 2814, checked in by jyg, 8 years ago
Cosmetic changes to fisrtilp
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 54.6 KB

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

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

Download in other formats:

Original Format