Context Navigation

lmdz_lscp_old.f90 @ 5500

Last change on this file since 5500 was 5480, checked in by yann meurdesoif, 7 days ago

For GPU porting :

add wrapper for source to source tools
Separate initialization phase of lscp_old (firstilp), SAVE variable in compute subroutine cannot be managed properly.

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.3 KB

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

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

Download in other formats:

Original Format