Context Navigation

fisrtilp.F90 @ 5361

Last change on this file since 5361 was 2969, checked in by fhourdin, 7 years ago
Petit bug dans iflag_t_glace>=1 vu uniquement en debug en 1D quand exposant < 1. C'est bien de l'avoir quand même dans les sources au cas où.
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: 57.2 KB

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

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

Download in other formats:

Original Format