Context Navigation

lmdz_lscp_old.f90 @ 5353

Last change on this file since 5353 was 5285, checked in by abarral, 7 weeks ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
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: 64.6 KB

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

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

Download in other formats:

Original Format