Context Navigation

fisrtilp.F90 @ 5445

Last change on this file since 5445 was 2880, checked in by Laurent Fairhead, 8 years ago
Missing initialisations, made the code unreproducible
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.8 KB

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

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

Download in other formats:

Original Format