Context Navigation

lmdz_lscp_old.f90 @ 5833

Last change on this file since 5833 was 5765, checked in by rkazeroni, 5 months ago

For GPU porting:

Add one Fortran comment for each ported param to specify possible names for horizontal variables used in the param
Add Fortran comment to exclude a file and its routines from GPU porting

Note:

The "gpum horizontal" comment is placed in the param main. It also applies to dependencies of the main.
The "gpum nogpu" excludes all modules and routines of the file from GPU porting. This is used for routines not targeted for GPU computing.

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

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

Download in other formats:

Original Format