Context Navigation

lmdz_lscp_old.F90 @ 5210

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

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

Download in other formats:

Original Format