Context Navigation

fisrtilp.F90 @ 5454

Last change on this file since 5454 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
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: 65.1 KB

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

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

Download in other formats:

Original Format