Context Navigation

fisrtilp.F90 @ 5225

Last change on this file since 5225 was 2925, checked in by fcheruy, 7 years ago
Update tree branch to trunk version
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: 55.4 KB

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

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

Download in other formats:

Original Format