Context Navigation

lmdz_lscp.F90 @ 4803

Last change on this file since 4803 was 4803, checked in by evignon, 17 months ago
implementation sous flag des premiers changements concernant le traitement des precipitations grande echelle dans le cadre de l'atelier nuages Audran, Lea, Niels, Gwendal et Etienne
File size: 60.0 KB

Rev	Line
[4664]	1	MODULE lmdz_lscp
[3999]	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4380]	8	SUBROUTINE lscp(klon,klev,dtime,missing_val, &
[4686]	9	paprs,pplay,temp,qt,ptconv,ratqs, &
[4380]	10	d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, &
[4530]	11	pfraclr,pfracld, &
	12	radocond, radicefrac, rain, snow, &
[3999]	13	frac_impa, frac_nucl, beta, &
	14	prfl, psfl, rhcl, zqta, fraca, &
	15	ztv, zpspsk, ztla, zthl, iflag_cld_th, &
[4380]	16	iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, &
[4639]	17	distcltop,temp_cltop, &
[4380]	18	qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
[4651]	19	Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
[4803]	20	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
	21	dqreva,dqssub,dqrauto,dqrcol,dqrmelt,dqrfreez,dqsauto, &
	22	dqsagg,dqsrim,dqsmelt,dqsfreez)
[3999]	23
	24	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	25	! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
	26	! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
	27	!--------------------------------------------------------------------------------
	28	! Date: 01/2021
	29	!--------------------------------------------------------------------------------
	30	! Aim: Large Scale Clouds and Precipitation (LSCP)
[4226]	31	!
[3999]	32	! This code is a new version of the fisrtilp.F90 routine, which is itself a
[4072]	33	! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
[3999]	34	! and 'ilp' (il pleut, L. Li)
	35	!
[4380]	36	! Compared to the original fisrtilp code, lscp
[3999]	37	! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
	38	! -> consider always precipitation thermalisation (fl_cor_ebil>0)
	39	! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
	40	! -> option "oldbug" by JYG has been removed
	41	! -> iflag_t_glace >0 always
	42	! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
	43	! -> rectangular distribution from L. Li no longer available
	44	! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
	45	!--------------------------------------------------------------------------------
	46	! References:
	47	!
[4412]	48	! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
[3999]	49	! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
	50	! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
[4412]	51	! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
[3999]	52	! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
[4412]	53	! - Touzze-Peifert Ludo, PhD thesis, p117-124
[3999]	54	! -------------------------------------------------------------------------------
	55	! Code structure:
	56	!
[4226]	57	! P0> Thermalization of the precipitation coming from the overlying layer
[3999]	58	! P1> Evaporation of the precipitation (falling from the k+1 level)
	59	! P2> Cloud formation (at the k level)
[4412]	60	! P2.A.1> With the PDFs, calculation of cloud properties using the inital
[3999]	61	! values of T and Q
	62	! P2.A.2> Coupling between condensed water and temperature
	63	! P2.A.3> Calculation of final quantities associated with cloud formation
[4412]	64	! P2.B> Release of Latent heat after cloud formation
[3999]	65	! P3> Autoconversion to precipitation (k-level)
	66	! P4> Wet scavenging
	67	!------------------------------------------------------------------------------
	68	! Some preliminary comments (JBM) :
	69	!
	70	! The cloud water that the radiation scheme sees is not the same that the cloud
	71	! water used in the physics and the dynamics
	72	!
[4412]	73	! During the autoconversion to precipitation (P3 step), radocond (cloud water used
[3999]	74	! by the radiation scheme) is calculated as an average of the water that remains
	75	! in the cloud during the precipitation and not the water remaining at the end
	76	! of the time step. The latter is used in the rest of the physics and advected
	77	! by the dynamics.
	78	!
	79	! In summary:
	80	!
	81	! Radiation:
	82	! xflwc(newmicro)+xfiwc(newmicro) =
[4412]	83	! radocond=lwcon(nc)+iwcon(nc)
[3999]	84	!
	85	! Notetheless, be aware of:
	86	!
[4412]	87	! radocond .NE. ocond(nc)
[3999]	88	! i.e.:
	89	! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
[4412]	90	! but oliq+(ocond-oliq) .EQ. ocond
[3999]	91	! (which is not trivial)
	92	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4226]	93	!
[4654]	94
	95	! USE de modules contenant des fonctions.
[4651]	96	USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
[4664]	97	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat
	98	USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
[4380]	99	USE ice_sursat_mod, ONLY : ice_sursat
[4803]	100	USE lmdz_lscp_poprecip, ONLY : poprecip_evapsub, poprecip_postcld
[3999]	101
[4664]	102	! Use du module lmdz_lscp_ini contenant les constantes
[4666]	103	USE lmdz_lscp_ini, ONLY : prt_level, lunout
[4664]	104	USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
	105	USE lmdz_lscp_ini, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
	106	USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
	107	USE lmdz_lscp_ini, ONLY : coef_eva, coef_eva_i,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
	108	USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
	109	USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc
	110	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
[4803]	111	USE lmdz_lscp_ini, ONLY : ok_poprecip
[4380]	112
[3999]	113	IMPLICIT NONE
	114
	115	!===============================================================================
[4380]	116	! VARIABLES DECLARATION
[3999]	117	!===============================================================================
	118
	119	! INPUT VARIABLES:
	120	!-----------------
	121
[4380]	122	INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels
[3999]	123	REAL, INTENT(IN) :: dtime ! time step [s]
[4380]	124	REAL, INTENT(IN) :: missing_val ! missing value for output
[4059]	125
[3999]	126	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
	127	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
[4654]	128	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K)
[4686]	129	REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg]
[3999]	130	INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds
	131	INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
[4226]	132	! CR: if iflag_ice_thermo=2, only convection is active
[4072]	133	LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated
[4059]	134
[3999]	135	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
	136
	137	!Inputs associated with thermal plumes
[4226]	138
[3999]	139	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! virtual potential temperature [K]
	140	REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg]
	141	REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-]
	142	REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp)
	143	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K]
	144
	145	! INPUT/OUTPUT variables
	146	!------------------------
[4562]	147
	148	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K]
[4686]	149	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs ! function of pressure that sets the large-scale
[4059]	150
[4380]	151
[4059]	152	! Input sursaturation en glace
	153	REAL, DIMENSION(klon,klev), INTENT(INOUT):: rneb_seri ! fraction nuageuse en memoire
[3999]	154
	155	! OUTPUT variables
	156	!-----------------
	157
	158	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K]
	159	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg]
	160	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg]
	161	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg]
	162	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-]
[4380]	163	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-]
[4530]	164	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-]
	165	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-]
[4412]	166	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg]
[3999]	167	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme
	168	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-]
[4412]	169	REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2]
	170	REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2]
[4380]	171	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg]
	172	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg]
[3999]	173	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
	174	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
[4562]	175	REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m]
[4639]	176	REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K]
[4686]	177	REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water
	178
[3999]	179	! fraction of aerosol scavenging through impaction and nucleation (for on-line)
	180
	181	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-]
	182	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
	183
[4380]	184	! for supersaturation and contrails parameterisation
[3999]	185
[4380]	186	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg]
	187	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg]
	188	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg]
	189	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg]
	190	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-]
	191	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-]
	192	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-]
[4425]	193	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K]
	194	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg]
	195	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q)
	196	REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails
	197	REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails
[4686]	198	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
	199	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
	200	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
	201	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
[3999]	202
[4803]	203	! for POPRECIP
[4686]	204
[4803]	205	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqreva !-- rain tendendy due to evaporation [kg/kg/s]
	206	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqssub !-- snow tendency due to sublimation [kg/kg/s]
	207	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrcol !-- rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
	208	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsagg !-- snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
	209	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrauto !-- rain tendency due to autoconversion of cloud liquid [kg/kg/s]
	210	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsauto !-- snow tendency due to autoconversion of cloud ice [kg/kg/s]
	211	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsrim !-- snow tendency due to riming [kg/kg/s]
	212	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsmelt !-- snow tendency due to melting [kg/kg/s]
	213	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrmelt !-- rain tendency due to melting [kg/kg/s]
	214	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !-- snow tendency due to freezing [kg/kg/s]
	215	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !-- rain tendency due to freezing [kg/kg/s]
[4686]	216
[4803]	217
[3999]	218	! LOCAL VARIABLES:
	219	!----------------
	220
[4654]	221	REAL,DIMENSION(klon) :: qsl, qsi
	222	REAL :: zct, zcl,zexpo
	223	REAL, DIMENSION(klon,klev) :: ctot
	224	REAL, DIMENSION(klon,klev) :: ctot_vol
	225	REAL, DIMENSION(klon) :: zqs, zdqs
	226	REAL :: zdelta, zcor, zcvm5
	227	REAL, DIMENSION(klon) :: zdqsdT_raw
	228	REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T
	229	REAL, DIMENSION(klon) :: Tbef,qlbef,DT
	230	REAL :: num,denom
	231	REAL :: cste
	232	REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta
	233	REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2
	234	REAL :: erf
	235	REAL, DIMENSION(klon,klev) :: icefrac_mpc
	236	REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
	237	REAL, DIMENSION(klon) :: zrfl, zrfln
	238	REAL :: zqev, zqevt
	239	REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
	240	REAL :: zqev0,zqevi, zqevti
	241	REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
	242	REAL, DIMENSION(klon) :: zoliql, zoliqi
	243	REAL, DIMENSION(klon) :: zt
	244	REAL, DIMENSION(klon,klev) :: zrho
	245	REAL, DIMENSION(klon) :: zdz,iwc
	246	REAL :: zchau,zfroi
	247	REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
	248	REAL :: zmelt,zrain,zsnow,zprecip
	249	REAL, DIMENSION(klon) :: dzfice
	250	REAL :: zsolid
	251	REAL, DIMENSION(klon) :: qtot, qzero
	252	REAL, DIMENSION(klon) :: dqsl,dqsi
	253	REAL :: smallestreal
[3999]	254	! Variables for Bergeron process
[4654]	255	REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
	256	REAL, DIMENSION(klon) :: zqpreci, zqprecl
[3999]	257	! Variables precipitation energy conservation
[4654]	258	REAL, DIMENSION(klon) :: zmqc
	259	REAL :: zalpha_tr
	260	REAL :: zfrac_lessi
	261	REAL, DIMENSION(klon) :: zprec_cond
[4803]	262	REAL :: zmair
[4654]	263	REAL :: zcpair, zcpeau
	264	REAL, DIMENSION(klon) :: zlh_solid
[4803]	265	REAL, DIMENSION(klon) :: ztupnew
[4654]	266	REAL :: zm_solid ! for liquid -> solid conversion
	267	REAL, DIMENSION(klon) :: zrflclr, zrflcld
	268	REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
	269	REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
	270	REAL, DIMENSION(klon) :: ziflclr, ziflcld
	271	REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
	272	REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
	273	REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
	274	REAL, DIMENSION(klon,klev) :: velo
	275	REAL :: vr, ffallv
	276	REAL :: qlmpc, qimpc, rnebmpc
	277	REAL, DIMENSION(klon,klev) :: radocondi, radocondl
	278	REAL :: effective_zneb
	279	REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D
[3999]	280
[4562]	281
[4654]	282	INTEGER i, k, n, kk, iter
	283	INTEGER, DIMENSION(klon) :: n_i
[4226]	284	INTEGER ncoreczq
[4654]	285	INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
[4803]	286	LOGICAL iftop
[3999]	287
[4654]	288	LOGICAL, DIMENSION(klon) :: lognormale
	289	LOGICAL, DIMENSION(klon) :: keepgoing
[3999]	290
	291	CHARACTER (len = 20) :: modname = 'lscp'
	292	CHARACTER (len = 80) :: abort_message
	293
	294
	295	!===============================================================================
	296	! INITIALISATION
	297	!===============================================================================
	298
[4226]	299	! Few initial checks
[3999]	300
[4654]	301
[3999]	302	IF (iflag_fisrtilp_qsat .LT. 0) THEN
[4226]	303	abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
	304	CALL abort_physic(modname,abort_message,1)
[3999]	305	ENDIF
	306
	307	! Few initialisations
	308
[4226]	309	znebprecip(:)=0.0
[3999]	310	ctot_vol(1:klon,1:klev)=0.0
	311	rneblsvol(1:klon,1:klev)=0.0
[4226]	312	smallestreal=1.e-9
	313	znebprecipclr(:)=0.0
	314	znebprecipcld(:)=0.0
[3999]	315	mpc_bl_points(:,:)=0
	316
	317	IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
	318
	319	! AA for 'safety' reasons
	320	zalpha_tr = 0.
	321	zfrac_lessi = 0.
[4380]	322	beta(:,:)= 0.
[3999]	323
[4380]	324	! Initialisation of variables:
	325
[3999]	326	prfl(:,:) = 0.0
	327	psfl(:,:) = 0.0
	328	d_t(:,:) = 0.0
	329	d_q(:,:) = 0.0
	330	d_ql(:,:) = 0.0
	331	d_qi(:,:) = 0.0
	332	rneb(:,:) = 0.0
[4530]	333	pfraclr(:,:)=0.0
	334	pfracld(:,:)=0.0
[4412]	335	radocond(:,:) = 0.0
[3999]	336	radicefrac(:,:) = 0.0
[4226]	337	frac_nucl(:,:) = 1.0
	338	frac_impa(:,:) = 1.0
[3999]	339	rain(:) = 0.0
	340	snow(:) = 0.0
[4114]	341	zoliq(:)=0.0
	342	zfice(:)=0.0
	343	dzfice(:)=0.0
	344	zqprecl(:)=0.0
	345	zqpreci(:)=0.0
[3999]	346	zrfl(:) = 0.0
	347	zifl(:) = 0.0
[4114]	348	ziflprev(:)=0.0
[3999]	349	zneb(:) = seuil_neb
[4226]	350	zrflclr(:) = 0.0
	351	ziflclr(:) = 0.0
	352	zrflcld(:) = 0.0
	353	ziflcld(:) = 0.0
	354	tot_zneb(:) = 0.0
	355	tot_znebn(:) = 0.0
	356	d_tot_zneb(:) = 0.0
	357	qzero(:) = 0.0
[4562]	358	distcltop1D(:)=0.0
[4639]	359	temp_cltop1D(:) = 0.0
[4803]	360	ztupnew(:)=0.0
[4380]	361	!--ice supersaturation
[4226]	362	gamma_ss(:,:) = 1.0
	363	qss(:,:) = 0.0
	364	rnebss(:,:) = 0.0
[4059]	365	Tcontr(:,:) = missing_val
	366	qcontr(:,:) = missing_val
	367	qcontr2(:,:) = missing_val
	368	fcontrN(:,:) = 0.0
	369	fcontrP(:,:) = 0.0
[4654]	370	distcltop(:,:)=0.
	371	temp_cltop(:,:)=0.
[4803]	372	!-- poprecip
	373	dqreva(:,:)=0.0
	374	dqrauto(:,:)=0.0
	375	dqrmelt(:,:)=0.0
	376	dqrfreez(:,:)=0.0
	377	dqrcol(:,:)=0.0
	378	dqssub(:,:)=0.0
	379	dqsauto(:,:)=0.0
	380	dqsrim(:,:)=0.0
	381	dqsagg(:,:)=0.0
	382	dqsfreez(:,:)=0.0
	383	dqsmelt(:,:)=0.0
[3999]	384
[4803]	385
	386
	387
[4392]	388	!c_iso: variable initialisation for iso
	389
	390
[3999]	391	!===============================================================================
	392	! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
	393	!===============================================================================
	394
	395	ncoreczq=0
	396
	397	DO k = klev, 1, -1
	398
[4803]	399	IF (k.LE.klev-1) THEN
	400	iftop=.false.
	401	ELSE
	402	iftop=.true.
	403	ENDIF
	404
[3999]	405	! Initialisation temperature and specific humidity
	406	DO i = 1, klon
[4654]	407	zt(i)=temp(i,k)
[4686]	408	zq(i)=qt(i,k)
[4803]	409	IF (.not. iftop) THEN
	410	ztupnew(i)=temp(i,k+1)+d_t(i,k+1)
	411	ENDIF
[4392]	412	!c_iso init of iso
[3999]	413	ENDDO
[4803]	414
	415	!================================================================
	416	! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
	417	IF (ok_poprecip) THEN
	418
	419	CALL poprecip_evapsub(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	420	zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
	421	zrfl, zrflclr, zrflcld, &
	422	zifl, ziflclr, ziflcld, &
	423	dqreva(:,k),dqssub(:,k) &
	424	)
	425
	426	!================================================================
	427	ELSE
	428
[3999]	429	! --------------------------------------------------------------------
[4803]	430	! P1> Thermalization of precipitation falling from the overlying layer
[3999]	431	! --------------------------------------------------------------------
[4412]	432	! Computes air temperature variation due to enthalpy transported by
[3999]	433	! precipitation. Precipitation is then thermalized with the air in the
	434	! layer.
	435	! The precipitation should remain thermalized throughout the different
[4412]	436	! thermodynamical transformations.
	437	! The corresponding water mass should
[3999]	438	! be added when calculating the layer's enthalpy change with
	439	! temperature
[4412]	440	! See lmdzpedia page todoan
	441	! todoan: check consistency with ice phase
	442	! todoan: understand why several steps
[3999]	443	! ---------------------------------------------------------------------
	444
[4803]	445	IF (iftop) THEN
[3999]	446
	447	DO i = 1, klon
[4803]	448	zmqc(i) = 0.
	449	ENDDO
	450
	451	ELSE
	452
	453	DO i = 1, klon
[3999]	454
[4803]	455	zmair=(paprs(i,k)-paprs(i,k+1))/RG
[3999]	456	! no condensed water so cp=cp(vapor+dry air)
	457	! RVTMP2=rcpv/rcpd-1
	458	zcpair=RCPD(1.0+RVTMP2zq(i))
[4226]	459	zcpeau=RCPD*RVTMP2
	460
[3999]	461	! zmqc: precipitation mass that has to be thermalized with
	462	! layer's air so that precipitation at the ground has the
	463	! same temperature as the lowermost layer
[4803]	464	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
[3999]	465	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
[4803]	466	zt(i) = ( ztupnew(i)zmqc(i)zcpeau + zcpair*zt(i) ) &
[3999]	467	/ (zcpair + zmqc(i)*zcpeau)
[4226]	468
[3999]	469	ENDDO
	470
	471	ENDIF
	472
	473	! --------------------------------------------------------------------
[4803]	474	! P2> Precipitation evaporation/sublimation/melting
[3999]	475	! --------------------------------------------------------------------
[4803]	476	! A part of the precipitation coming from above is evaporated/sublimated/melted.
[3999]	477	! --------------------------------------------------------------------
[4803]	478
[4397]	479	IF (iflag_evap_prec.GE.1) THEN
[3999]	480
[4072]	481	! Calculation of saturation specific humidity
	482	! depending on temperature:
[4380]	483	CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
[4072]	484	! wrt liquid water
[4380]	485	CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
[4072]	486	! wrt ice
[4380]	487	CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
[3999]	488
	489	DO i = 1, klon
[4226]	490
[3999]	491	! if precipitation
	492	IF (zrfl(i)+zifl(i).GT.0.) THEN
[4226]	493
[4563]	494	! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
[4392]	495	! c_iso: likely important to distinguish cs from neb isotope precipitation
	496
[4563]	497	IF (iflag_evap_prec.GE.4) THEN
[4226]	498	zrfl(i) = zrflclr(i)
	499	zifl(i) = ziflclr(i)
	500	ENDIF
[3999]	501
	502	IF (iflag_evap_prec.EQ.1) THEN
	503	znebprecip(i)=zneb(i)
	504	ELSE
	505	znebprecip(i)=MAX(zneb(i),znebprecip(i))
	506	ENDIF
	507
[4563]	508	IF (iflag_evap_prec.GT.4) THEN
	509	! Max evaporation not to saturate the clear sky precip fraction
	510	! i.e. the fraction where evaporation occurs
	511	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
	512	ELSEIF (iflag_evap_prec .EQ. 4) THEN
[4226]	513	! Max evaporation not to saturate the whole mesh
[4563]	514	! Pay attention -> lead to unrealistic and excessive evaporation
[4226]	515	zqev0 = MAX(0.0, zqs(i)-zq(i))
	516	ELSE
	517	! Max evap not to saturate the fraction below the cloud
	518	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
	519	ENDIF
[3999]	520
	521	! Evaporation of liquid precipitation coming from above
	522	! dP/dz=beta(1-q/qsat)sqrt(P)
	523	! formula from Sundquist 1988, Klemp & Wilhemson 1978
[4563]	524	! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
[3999]	525
	526	IF (iflag_evap_prec.EQ.3) THEN
[4072]	527	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl(i)) &
[3999]	528	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
	529	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
[4563]	530	ELSE IF (iflag_evap_prec.GE.4) THEN
[4072]	531	zqevt = znebprecipclr(i)coef_eva(1.0-zq(i)/qsl(i)) &
[3999]	532	*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
	533	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	534	ELSE
[4072]	535	zqevt = 1.coef_eva(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
[3999]	536	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	537	ENDIF
	538
[4226]	539	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
	540	RGdtime/(paprs(i,k)-paprs(i,k+1))
[3999]	541
	542	! sublimation of the solid precipitation coming from above
	543	IF (iflag_evap_prec.EQ.3) THEN
[4072]	544	zqevti = znebprecip(i)coef_eva_i(1.0-zq(i)/qsi(i)) &
[3999]	545	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
	546	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
[4563]	547	ELSE IF (iflag_evap_prec.GE.4) THEN
[4072]	548	zqevti = znebprecipclr(i)coef_eva_i(1.0-zq(i)/qsi(i)) &
[3999]	549	*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
	550	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	551	ELSE
[4072]	552	zqevti = 1.coef_eva_i(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
[3999]	553	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	554	ENDIF
[4226]	555
[3999]	556	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
[4226]	557	RGdtime/(paprs(i,k)-paprs(i,k+1))
[3999]	558
	559	! A. JAM
	560	! Evaporation limit: we ensure that the layer's fraction below
	561	! the cloud or the whole mesh (depending on iflag_evap_prec)
	562	! does not reach saturation. In this case, we
	563	! redistribute zqev0 conserving the ratio liquid/ice
[4226]	564
[3999]	565	IF (zqevt+zqevti.GT.zqev0) THEN
	566	zqev=zqev0*zqevt/(zqevt+zqevti)
	567	zqevi=zqev0*zqevti/(zqevt+zqevti)
	568	ELSE
	569	zqev=zqevt
	570	zqevi=zqevti
	571	ENDIF
	572
	573
[4392]	574	! New solid and liquid precipitation fluxes after evap and sublimation
[3999]	575	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
	576	/RG/dtime)
	577	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
	578	/RG/dtime)
	579
[4392]	580
[3999]	581	! vapor, temperature, precip fluxes update
[4803]	582	! vapor is updated after evaporation/sublimation (it is increased)
[3999]	583	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	584	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[4803]	585	! zmqc is the total condensed water in the precip flux (it is decreased)
[3999]	586	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	587	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[4803]	588	! air and precip temperature (i.e., gridbox temperature)
	589	! is updated due to latent heat cooling
[3999]	590	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	591	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	592	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
	593	+ (zifln(i)-zifl(i)) &
	594	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	595	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
	596
	597	! New values of liquid and solid precipitation
	598	zrfl(i) = zrfln(i)
	599	zifl(i) = zifln(i)
	600
[4392]	601	! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
	602	! due to evap + sublim
	603
	604
[4563]	605	IF (iflag_evap_prec.GE.4) THEN
[4226]	606	zrflclr(i) = zrfl(i)
	607	ziflclr(i) = zifl(i)
	608	IF(zrflclr(i) + ziflclr(i).LE.0) THEN
	609	znebprecipclr(i) = 0.0
	610	ENDIF
	611	zrfl(i) = zrflclr(i) + zrflcld(i)
	612	zifl(i) = ziflclr(i) + ziflcld(i)
	613	ENDIF
[3999]	614
[4392]	615	! c_iso duplicate for isotopes or loop on isotopes
[3999]	616
[4392]	617	! Melting:
[4412]	618	zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
[3999]	619	! precip fraction that is melted
	620	zmelt = MIN(MAX(zmelt,0.),1.)
[4392]	621
[4412]	622	! update of rainfall and snowfall due to melting
[4563]	623	IF (iflag_evap_prec.GE.4) THEN
[4226]	624	zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
	625	zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
	626	zrfl(i)=zrflclr(i)+zrflcld(i)
[4412]	627
	628	ziflclr(i)=ziflclr(i)*(1.-zmelt)
	629	ziflcld(i)=ziflcld(i)*(1.-zmelt)
	630	zifl(i)=ziflclr(i)+ziflcld(i)
	631
[4226]	632	ELSE
	633	zrfl(i)=zrfl(i)+zmelt*zifl(i)
[4412]	634
	635	zifl(i)=zifl(i)*(1.-zmelt)
[4226]	636	ENDIF
[4412]	637
	638
[4392]	639	! c_iso: melting of isotopic precipi with zmelt (no fractionation)
[3999]	640
[4803]	641	! Latent heat of melting because of precipitation melting
	642	! NB: the air + precip temperature is simultaneously updated
[3999]	643	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
	644	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
	645
	646
	647	ELSE
	648	! if no precip, we reinitialize the cloud fraction used for the precip to 0
	649	znebprecip(i)=0.
	650
	651	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	652
[4803]	653	ENDDO ! loop on klon
[4226]	654
[3999]	655	ENDIF ! (iflag_evap_prec>=1)
	656
[4803]	657	ENDIF ! (ok_poprecip)
	658
[3999]	659	! --------------------------------------------------------------------
	660	! End precip evaporation
	661	! --------------------------------------------------------------------
	662
	663	! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
	664	!-------------------------------------------------------
[4226]	665
[4072]	666	qtot(:)=zq(:)+zmqc(:)
[4380]	667	CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
[4072]	668	DO i = 1, klon
[3999]	669	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
[4059]	670	zdqsdT_raw(i) = zdqs(i)RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTT*zdelta)
[3999]	671	IF (zq(i) .LT. 1.e-15) THEN
	672	ncoreczq=ncoreczq+1
	673	zq(i)=1.e-15
	674	ENDIF
[4392]	675	! c_iso: do something similar for isotopes
[3999]	676
[4072]	677	ENDDO
[3999]	678
	679	! --------------------------------------------------------------------
	680	! P2> Cloud formation
	681	!---------------------------------------------------------------------
	682	!
	683	! Unlike fisrtilp, we always assume a 'fractional cloud' approach
	684	! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
	685	! is prescribed and depends on large scale variables and boundary layer
	686	! properties)
	687	! The decrease in condensed part due tu latent heating is taken into
	688	! account
	689	! -------------------------------------------------------------------
	690
	691	! P2.1> With the PDFs (log-normal, bigaussian)
[4424]	692	! cloud properties calculation with the initial values of t and q
[3999]	693	! ----------------------------------------------------------------
	694
	695	! initialise gammasat and qincloud_mpc
	696	gammasat(:)=1.
	697	qincloud_mpc(:)=0.
	698
	699	IF (iflag_cld_th.GE.5) THEN
	700
	701	IF (iflag_mpc_bl .LT. 1) THEN
	702
	703	IF (iflag_cloudth_vert.LE.2) THEN
	704
	705	CALL cloudth(klon,klev,k,ztv, &
	706	zq,zqta,fraca, &
	707	qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[4654]	708	ratqs,zqs,temp, &
[4651]	709	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	710
[4651]	711
[3999]	712	ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
	713
	714	CALL cloudth_v3(klon,klev,k,ztv, &
	715	zq,zqta,fraca, &
	716	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[4654]	717	ratqs,zqs,temp, &
[4651]	718	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	719
[4226]	720	!Jean Jouhaud's version, Decembre 2018
	721	!-------------------------------------
[3999]	722
	723	ELSEIF (iflag_cloudth_vert.EQ.6) THEN
	724
	725	CALL cloudth_v6(klon,klev,k,ztv, &
	726	zq,zqta,fraca, &
	727	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[4654]	728	ratqs,zqs,temp, &
[4651]	729	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	730
	731	ENDIF
	732
	733	ELSE
	734	! cloudth_v3 + cold microphysical considerations
	735	! + stationary mixed-phase cloud model
	736
[4380]	737	CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
[4654]	738	temp,ztv,zq,zqta,fraca, zpspsk, &
[3999]	739	paprs,pplay,ztla,zthl,ratqs,zqs,psfl, &
[4651]	740	qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol, &
	741	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	742
	743	ENDIF ! iflag_mpc_bl
	744
	745	DO i=1,klon
	746	rneb(i,k)=ctot(i,k)
	747	rneblsvol(i,k)=ctot_vol(i,k)
	748	zqn(i)=qcloud(i)
	749	ENDDO
	750
	751	ENDIF
	752
	753	IF (iflag_cld_th .LE. 4) THEN
	754
	755	! lognormal
	756	lognormale = .TRUE.
	757
	758	ELSEIF (iflag_cld_th .GE. 6) THEN
	759
	760	! lognormal distribution when no thermals
	761	lognormale = fraca(:,k) < 1e-10
	762
	763	ELSE
[4226]	764	! When iflag_cld_th=5, we always assume
[3999]	765	! bi-gaussian distribution
	766	lognormale = .FALSE.
	767
	768	ENDIF
	769
	770	DT(:) = 0.
	771	n_i(:)=0
	772	Tbef(:)=zt(:)
	773	qlbef(:)=0.
	774
	775	! Treatment of non-boundary layer clouds (lognormale)
	776	! condensation with qsat(T) variation (adaptation)
[4424]	777	! Iterative resolution to converge towards qsat
	778	! with update of temperature, ice fraction and qsat at
	779	! each iteration
[4226]	780
[4424]	781	! todoan -> sensitivity to iflag_fisrtilp_qsat
[3999]	782	DO iter=1,iflag_fisrtilp_qsat+1
[4226]	783
[3999]	784	DO i=1,klon
	785
[4424]	786	! keepgoing = .true. while convergence is not satisfied
	787	keepgoing(i)=ABS(DT(i)).GT.DDT0
[4226]	788
[4424]	789	IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
[4226]	790
[3999]	791	! if not convergence:
	792
	793	! P2.2.1> cloud fraction and condensed water mass calculation
	794	! Calculated variables:
	795	! rneb : cloud fraction
	796	! zqn : total water within the cloud
	797	! zcond: mean condensed water within the mesh
	798	! rhcl: clear-sky relative humidity
	799	!---------------------------------------------------------------
	800
[4424]	801	! new temperature that only serves in the iteration process:
[3999]	802	Tbef(i)=Tbef(i)+DT(i)
	803
	804	! Rneb, qzn and zcond for lognormal PDFs
[4072]	805	qtot(i)=zq(i)+zmqc(i)
[4226]	806
[4072]	807	ENDIF
[3999]	808
[4072]	809	ENDDO
	810
[4380]	811	! Calculation of saturation specific humidity and ice fraction
[4226]	812	CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
	813	CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
[4072]	814	! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
	815	zdqs(:) = gammasat(:)zdqs(:)+zqs(:)dgammasatdt(:)
[4562]	816	! cloud phase determination
	817	IF (iflag_t_glace.GE.4) THEN
	818	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
[4654]	819	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
[4562]	820	ENDIF
[4639]	821	CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:))
[4072]	822
[4424]	823	DO i=1,klon !todoan : check if loop in i is needed
[4072]	824
[4424]	825	IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
[4072]	826
[3999]	827	zpdf_sig(i)=ratqs(i,k)*zq(i)
	828	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
	829	zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
	830	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
	831	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
	832	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
	833	zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
	834	zpdf_e1(i)=1.-erf(zpdf_e1(i))
	835	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
	836	zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
	837	zpdf_e2(i)=1.-erf(zpdf_e2(i))
[4226]	838
[4062]	839	IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN
[4059]	840
	841	IF (zpdf_e1(i).LT.1.e-10) THEN
	842	rneb(i,k)=0.
	843	zqn(i)=gammasat(i)*zqs(i)
	844	ELSE
	845	rneb(i,k)=0.5*zpdf_e1(i)
	846	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
	847	ENDIF
	848
[4072]	849	rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time)
[4059]	850	fcontrN(i,k)=0.0 !--idem
	851	fcontrP(i,k)=0.0 !--idem
	852	qss(i,k)=0.0 !--idem
[4226]	853
[4424]	854	ELSE ! in case of ice supersaturation by Audran
[4226]	855
[4059]	856	!------------------------------------
[4072]	857	! ICE SUPERSATURATION
[4059]	858	!------------------------------------
[3999]	859
[4654]	860	CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
[4226]	861	gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), &
	862	rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), &
[4059]	863	Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
	864
	865	ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
	866
[3999]	867	! If vertical heterogeneity, change fraction by volume as well
	868	IF (iflag_cloudth_vert.GE.3) THEN
	869	ctot_vol(i,k)=rneb(i,k)
	870	rneblsvol(i,k)=ctot_vol(i,k)
	871	ENDIF
	872
	873
[4072]	874	! P2.2.2> Approximative calculation of temperature variation DT
	875	! due to condensation.
	876	! Calculated variables:
	877	! dT : temperature change due to condensation
	878	!---------------------------------------------------------------
[3999]	879
	880
	881	IF (zfice(i).LT.1) THEN
	882	cste=RLVTT
	883	ELSE
	884	cste=RLSTT
	885	ENDIF
	886
	887	qlbef(i)=max(0.,zqn(i)-zqs(i))
	888	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
[4226]	889	+zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
[3999]	890	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
	891	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
[4226]	892	qlbef(i)dzfice(i)
[4424]	893	! here we update a provisory temperature variable that only serves in the iteration
	894	! process
[3999]	895	DT(i)=num/denom
	896	n_i(i)=n_i(i)+1
	897
[4424]	898	ENDIF ! end keepgoing
[3999]	899
	900	ENDDO ! end loop on i
	901
	902	ENDDO ! iter=1,iflag_fisrtilp_qsat+1
	903
	904	! P2.3> Final quantities calculation
	905	! Calculated variables:
	906	! rneb : cloud fraction
	907	! zcond: mean condensed water in the mesh
	908	! zqn : mean water vapor in the mesh
[4562]	909	! zfice: ice fraction in clouds
[3999]	910	! zt : temperature
	911	! rhcl : clear-sky relative humidity
	912	! ----------------------------------------------------------------
	913
	914
[4562]	915	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
	916	IF (iflag_t_glace.GE.4) THEN
[4654]	917	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
[4562]	918	distcltop(:,k)=distcltop1D(:)
[4639]	919	temp_cltop(:,k)=temp_cltop1D(:)
[4562]	920	ENDIF
	921
[3999]	922	! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
[4803]	923	CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice, dzfice)
[3999]	924
[4562]	925
	926	! Water vapor update, Phase determination and subsequent latent heat exchange
[3999]	927	DO i=1, klon
	928
	929	IF (mpc_bl_points(i,k) .GT. 0) THEN
	930	zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
	931	! following line is very strange and probably wrong
	932	rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
	933	! water vapor update and partition function if thermals
	934	zq(i) = zq(i) - zcond(i)
	935	zfice(i)=icefrac_mpc(i,k)
	936
	937	ELSE
	938
	939	! Checks on rneb, rhcl and zqn
	940	IF (rneb(i,k) .LE. 0.0) THEN
	941	zqn(i) = 0.0
	942	rneb(i,k) = 0.0
	943	zcond(i) = 0.0
	944	rhcl(i,k)=zq(i)/zqs(i)
	945	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	946	zqn(i) = zq(i)
	947	rneb(i,k) = 1.0
	948	zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))
	949	rhcl(i,k)=1.0
	950	ELSE
	951	zcond(i) = MAX(0.0,zqn(i)-gammasat(i)zqs(i))rneb(i,k)
	952	! following line is very strange and probably wrong:
	953	rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
	954	ENDIF
	955
[4226]	956	! water vapor update
	957	zq(i) = zq(i) - zcond(i)
[3999]	958
	959	ENDIF
	960
[4392]	961	! c_iso : routine that computes in-cloud supersaturation
	962	! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
	963
[3999]	964	! temperature update due to phase change
	965	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
	966	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
	967	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
	968	ENDDO
	969
	970	! If vertical heterogeneity, change volume fraction
	971	IF (iflag_cloudth_vert .GE. 3) THEN
	972	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
	973	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
	974	ENDIF
	975
	976	! ----------------------------------------------------------------
[4114]	977	! P3> Precipitation formation
[3999]	978	! ----------------------------------------------------------------
[4226]	979
[4803]	980	!================================================================
	981	IF (ok_poprecip) THEN
	982
	983	DO i = 1, klon
	984	zoliql(i) = zcond(i) * ( 1. - zfice(i) )
	985	zoliqi(i) = zcond(i) * zfice(i)
	986	ENDDO
	987
	988	CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	989	ctot_vol(:,k), ptconv(:,k), &
	990	zt, zq, zoliql, zoliqi, zfice, &
	991	rneb(:,k), znebprecipclr, znebprecipcld, &
	992	zrfl, zrflclr, zrflcld, &
	993	zifl, ziflclr, ziflcld, &
	994	dqrauto(:,k),dqrcol(:,k),dqrmelt(:,k),dqrfreez(:,k), &
	995	dqsauto(:,k),dqsagg(:,k),dqsrim(:,k),dqsmelt(:,k),dqsfreez(:,k) &
	996	)
	997
	998	DO i = 1, klon
	999	zoliq(i) = zoliql(i) + zoliqi(i)
	1000	IF ( zoliq(i) .GT. 0. ) THEN
	1001	zfice(i) = zoliqi(i)/zoliq(i)
	1002	ELSE
	1003	zfice(i) = 0.0
	1004	ENDIF
	1005	! when poprecip activated, radiation does not see any precipitation content
	1006	radocond(i,k) = zoliq(i)
	1007	radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
	1008	radocondi(i,k)= radocond(i,k)*zfice(i)
	1009	ENDDO
	1010
	1011	!================================================================
	1012	ELSE
	1013
[3999]	1014	! LTP:
[4563]	1015	IF (iflag_evap_prec .GE. 4) THEN
[3999]	1016
	1017	!Partitionning between precipitation coming from clouds and that coming from CS
	1018
	1019	!0) Calculate tot_zneb, total cloud fraction above the cloud
	1020	!assuming a maximum-random overlap (voir Jakob and Klein, 2000)
	1021
	1022	DO i=1, klon
[4424]	1023	tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
	1024	/(1.-min(zneb(i),1.-smallestreal))
[3999]	1025	d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
	1026	tot_zneb(i) = tot_znebn(i)
	1027
	1028
	1029	!1) Cloudy to clear air
	1030	d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
[4424]	1031	IF (znebprecipcld(i) .GT. 0.) THEN
[3999]	1032	d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
	1033	d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
	1034	ELSE
	1035	d_zrfl_cld_clr(i) = 0.
	1036	d_zifl_cld_clr(i) = 0.
	1037	ENDIF
	1038
	1039	!2) Clear to cloudy air
[4226]	1040	d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
[3999]	1041	IF (znebprecipclr(i) .GT. 0) THEN
	1042	d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
	1043	d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
	1044	ELSE
	1045	d_zrfl_clr_cld(i) = 0.
	1046	d_zifl_clr_cld(i) = 0.
	1047	ENDIF
	1048
	1049	!Update variables
[4226]	1050	znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
[3999]	1051	znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
	1052	zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
	1053	ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
	1054	zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
	1055	ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
	1056
[4392]	1057	! c_iso do the same thing for isotopes precip
[3999]	1058	ENDDO
	1059	ENDIF
	1060
[4559]	1061
	1062	! Autoconversion
	1063	! -------------------------------------------------------------------------------
	1064
	1065
[4412]	1066	! Initialisation of zoliq and radocond variables
[3999]	1067
	1068	DO i = 1, klon
	1069	zoliq(i) = zcond(i)
[4072]	1070	zoliqi(i)= zoliq(i)*zfice(i)
	1071	zoliql(i)= zoliq(i)*(1.-zfice(i))
[4392]	1072	! c_iso : initialisation of zoliq* also for isotopes
[4114]	1073	zrho(i,k) = pplay(i,k) / zt(i) / RD
	1074	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
[4072]	1075	iwc(i) = 0.
	1076	zneb(i) = MAX(rneb(i,k),seuil_neb)
[4559]	1077	radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
	1078	radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
	1079	radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
[3999]	1080	ENDDO
	1081
[4559]	1082
	1083	DO n = 1, niter_lscp
[3999]	1084
	1085	DO i=1,klon
	1086	IF (rneb(i,k).GT.0.0) THEN
[4114]	1087	iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
[3999]	1088	ENDIF
	1089	ENDDO
	1090
[4226]	1091	CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
[3999]	1092
	1093	DO i = 1, klon
	1094
	1095	IF (rneb(i,k).GT.0.0) THEN
	1096
[4072]	1097	! Initialization of zrain, zsnow and zprecip:
	1098	zrain=0.
	1099	zsnow=0.
	1100	zprecip=0.
[4392]	1101	! c_iso same init for isotopes. Externalisation?
[3999]	1102
	1103	IF (zneb(i).EQ.seuil_neb) THEN
[4072]	1104	zprecip = 0.0
	1105	zsnow = 0.0
	1106	zrain= 0.0
[3999]	1107	ELSE
[4420]	1108
[3999]	1109	IF (ptconv(i,k)) THEN ! if convective point
	1110	zcl=cld_lc_con
	1111	zct=1./cld_tau_con
[4559]	1112	zexpo=cld_expo_con
	1113	ffallv=ffallv_con
[3999]	1114	ELSE
	1115	zcl=cld_lc_lsc
	1116	zct=1./cld_tau_lsc
[4559]	1117	zexpo=cld_expo_lsc
	1118	ffallv=ffallv_lsc
[3999]	1119	ENDIF
	1120
[4559]	1121
[3999]	1122	! if vertical heterogeneity is taken into account, we use
	1123	! the "true" volume fraction instead of a modified
	1124	! surface fraction (which is larger and artificially
	1125	! reduces the in-cloud water).
[4072]	1126
[4559]	1127	! Liquid water quantity to remove: zchau (Sundqvist, 1978)
	1128	! dqliq/dt=-qliq/tau(1-exp(-qcin/clw)*2)
	1129	!.........................................................
[3999]	1130	IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
	1131
[4420]	1132	! if vertical heterogeneity is taken into account, we use
	1133	! the "true" volume fraction instead of a modified
	1134	! surface fraction (which is larger and artificially
	1135	! reduces the in-cloud water).
[4559]	1136	effective_zneb=ctot_vol(i,k)
	1137	ELSE
	1138	effective_zneb=zneb(i)
	1139	ENDIF
[4420]	1140
	1141
[4559]	1142	IF (iflag_autoconversion .EQ. 2) THEN
	1143	! two-steps resolution with niter_lscp=1 sufficient
	1144	! we first treat the second term (with exponential) in an explicit way
	1145	! and then treat the first term (-q/tau) in an exact way
	1146	zchau=zoliql(i)(1.-exp(-dtime/REAL(niter_lscp)zct &
	1147	(1.-exp(-(zoliql(i)/effective_zneb/zcl)*zexpo))))
[4420]	1148	ELSE
[4559]	1149	! old explicit resolution with subtimesteps
	1150	zchau = zctdtime/REAL(niter_lscp)zoliql(i) &
	1151	(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)*zexpo))
[4420]	1152	ENDIF
	1153
	1154
[4559]	1155	! Ice water quantity to remove (Zender & Kiehl, 1997)
	1156	! dqice/dt=1/rhod(rhowice*qice)/dz
	1157	!....................................
	1158	IF (iflag_autoconversion .EQ. 2) THEN
	1159	! exact resolution, niter_lscp=1 is sufficient but works only
	1160	! with iflag_vice=0
	1161	IF (zoliqi(i) .GT. 0.) THEN
	1162	zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
	1163	+dice_velodtime/REAL(niter_lscp)cice_velo/zdz(i)ffallv)*(-1./dice_velo))
	1164	ELSE
	1165	zfroi=0.
	1166	ENDIF
	1167	ELSE
	1168	! old explicit resolution with subtimesteps
	1169	zfroi = dtime/REAL(niter_lscp)/zdz(i)zoliqi(i)velo(i,k)
	1170	ENDIF
	1171
[4397]	1172	zrain = MIN(MAX(zchau,0.0),zoliql(i))
	1173	zsnow = MIN(MAX(zfroi,0.0),zoliqi(i))
[4072]	1174	zprecip = MAX(zrain + zsnow,0.0)
[3999]	1175
	1176	ENDIF
[4226]	1177
[4559]	1178	IF (iflag_autoconversion .GE. 1) THEN
	1179	! debugged version with phase conservation through the autoconversion process
	1180	zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0)
	1181	zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0)
	1182	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
	1183	ELSE
	1184	! bugged version with phase resetting
	1185	zoliql(i) = MAX(zoliq(i)(1.-zfice(i))-1.zrain , 0.0)
	1186	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zsnow , 0.0)
	1187	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
	1188	ENDIF
[3999]	1189
[4392]	1190	! c_iso: call isotope_conversion (for readibility) or duplicate
	1191
[4559]	1192	radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
	1193	radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
	1194	radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
[4226]	1195
[4072]	1196	ENDIF ! rneb >0
[3999]	1197
	1198	ENDDO ! i = 1,klon
	1199
	1200	ENDDO ! n = 1,niter
	1201
[4072]	1202	! Precipitation flux calculation
[3999]	1203
	1204	DO i = 1, klon
[4380]	1205
[4563]	1206	IF (iflag_evap_prec.GE.4) THEN
[4380]	1207	ziflprev(i)=ziflcld(i)
	1208	ELSE
	1209	ziflprev(i)=zifl(i)*zneb(i)
	1210	ENDIF
[3999]	1211
	1212	IF (rneb(i,k) .GT. 0.0) THEN
	1213
[4072]	1214	! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
	1215	! If T<0C, liquid precip are converted into ice, which leads to an increase in
	1216	! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
	1217	! taken into account through a linearization of the equations and by approximating
	1218	! the liquid precipitation process with a "threshold" process. We assume that
	1219	! condensates are not modified during this operation. Liquid precipitation is
	1220	! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
	1221	! Water vapor increases as well
	1222	! Note that compared to fisrtilp, we always assume iflag_bergeron=2
	1223
	1224	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
	1225	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
	1226	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
	1227	coef1 = rneb(i,k)RLSTT/zcpzdqsdT_raw(i)
	1228	! Computation of DT if all the liquid precip freezes
	1229	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
	1230	! T should not exceed the freezing point
	1231	! that is Delta > RTT-zt(i)
	1232	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
	1233	zt(i) = zt(i) + DeltaT
	1234	! water vaporization due to temp. increase
	1235	Deltaq = rneb(i,k)zdqsdT_raw(i)DeltaT
	1236	! we add this vaporized water to the total vapor and we remove it from the precip
	1237	zq(i) = zq(i) + Deltaq
	1238	! The three "max" lines herebelow protect from rounding errors
	1239	zcond(i) = max( zcond(i)- Deltaq, 0. )
	1240	! liquid precipitation converted to ice precip
	1241	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
	1242	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
	1243	! iced water budget
	1244	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
[4226]	1245
[4392]	1246	! c_iso : mv here condensation of isotopes + redispatchage en precip
	1247
[4563]	1248	IF (iflag_evap_prec.GE.4) THEN
[4226]	1249	zrflcld(i) = zrflcld(i)+zqprecl(i) &
	1250	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1251	ziflcld(i) = ziflcld(i)+ zqpreci(i) &
	1252	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1253	znebprecipcld(i) = rneb(i,k)
	1254	zrfl(i) = zrflcld(i) + zrflclr(i)
	1255	zifl(i) = ziflcld(i) + ziflclr(i)
	1256	ELSE
[3999]	1257	zrfl(i) = zrfl(i)+ zqprecl(i) &
	1258	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1259	zifl(i) = zifl(i)+ zqpreci(i) &
[4226]	1260	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
[3999]	1261	ENDIF
[4392]	1262	! c_iso : same for isotopes
[3999]	1263
	1264	ENDIF ! rneb>0
	1265
	1266	ENDDO
	1267
[4226]	1268	! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
[4563]	1269	! if iflag_evap_prec>=4
	1270	IF (iflag_evap_prec.GE.4) THEN
[4114]	1271
[4226]	1272	DO i=1,klon
[4114]	1273
[4226]	1274	IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
[3999]	1275	znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
	1276	(MAX(znebprecipclr(i),seuil_neb)rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)rain_int_min)))
	1277	ELSE
[4226]	1278	znebprecipclr(i)=0.0
	1279	ENDIF
	1280
	1281	IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
[3999]	1282	znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
	1283	(MAX(znebprecipcld(i),seuil_neb)rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)rain_int_min)))
[4226]	1284	ELSE
	1285	znebprecipcld(i)=0.0
	1286	ENDIF
[3999]	1287
[4226]	1288	ENDDO
[3999]	1289
[4226]	1290	ENDIF
[3999]	1291
[4803]	1292	ENDIF ! ok_poprecip
	1293
[3999]	1294	! End of precipitation formation
	1295	! --------------------------------
	1296
	1297
[4803]	1298	! Calculation of cloud condensates amount seen by radiative scheme
	1299	!-----------------------------------------------------------------
	1300
[4114]	1301	! Calculation of the concentration of condensates seen by the radiation scheme
[4412]	1302	! for liquid phase, we take radocondl
	1303	! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
[4803]	1304	! we recalculate radocondi to account for contributions from the precipitation flux
	1305	! TODO: for the moment, we deactivate this option when ok_poprecip=.true.
[3999]	1306
[4803]	1307	IF ((ok_radocond_snow) .AND. (k .LT. klev) .AND. (.NOT. ok_poprecip)) THEN
[4114]	1308	! for the solid phase (crystals + snowflakes)
[4412]	1309	! we recalculate radocondi by summing
[4114]	1310	! the ice content calculated in the mesh
	1311	! + the contribution of the non-evaporated snowfall
	1312	! from the overlying layer
	1313	DO i=1,klon
	1314	IF (ziflprev(i) .NE. 0.0) THEN
[4412]	1315	radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
[4114]	1316	ELSE
[4412]	1317	radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
[4114]	1318	ENDIF
[4412]	1319	radocond(i,k)=radocondl(i,k)+radocondi(i,k)
[4114]	1320	ENDDO
	1321	ENDIF
	1322
[4412]	1323	! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
[4114]	1324	DO i=1,klon
[4412]	1325	IF (radocond(i,k) .GT. 0.) THEN
	1326	radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
[4114]	1327	ENDIF
	1328	ENDDO
	1329
[3999]	1330	! ----------------------------------------------------------------
	1331	! P4> Wet scavenging
	1332	! ----------------------------------------------------------------
	1333
	1334	!Scavenging through nucleation in the layer
	1335
	1336	DO i = 1,klon
	1337
	1338	IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
	1339	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
	1340	ELSE
	1341	beta(i,k) = 0.
	1342	ENDIF
	1343
[4059]	1344	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
[3999]	1345
	1346	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1347
[4654]	1348	IF (temp(i,k) .GE. t_glace_min) THEN
[3999]	1349	zalpha_tr = a_tr_sca(3)
	1350	ELSE
	1351	zalpha_tr = a_tr_sca(4)
	1352	ENDIF
	1353
	1354	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1355	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
[4226]	1356
[3999]	1357	! Nucleation with a factor of -1 instead of -0.5
	1358	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
	1359
	1360	ENDIF
	1361
[4226]	1362	ENDDO
	1363
[3999]	1364	! Scavenging through impaction in the underlying layer
	1365
	1366	DO kk = k-1, 1, -1
	1367
	1368	DO i = 1, klon
	1369
	1370	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1371
[4654]	1372	IF (temp(i,kk) .GE. t_glace_min) THEN
[3999]	1373	zalpha_tr = a_tr_sca(1)
	1374	ELSE
	1375	zalpha_tr = a_tr_sca(2)
	1376	ENDIF
	1377
	1378	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1379	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
	1380
	1381	ENDIF
	1382
	1383	ENDDO
	1384
	1385	ENDDO
[4226]	1386
[4072]	1387	!--save some variables for ice supersaturation
[4059]	1388	!
	1389	DO i = 1, klon
[4072]	1390	! for memory
[4059]	1391	rneb_seri(i,k) = rneb(i,k)
[3999]	1392
[4072]	1393	! for diagnostics
[4059]	1394	rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
[3999]	1395
[4059]	1396	qvc(i,k) = zqs(i) * rneb(i,k)
[4072]	1397	qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--added by OB because of pathologic cases with the lognormal
[4059]	1398	qcld(i,k) = qvc(i,k) + zcond(i)
[4226]	1399	ENDDO
[4072]	1400	!q_sat
[4380]	1401	CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
	1402	CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
[4059]	1403
[4803]	1404
	1405
	1406	! Outputs:
	1407	!-------------------------------
	1408	! Precipitation fluxes at layer interfaces
	1409	! + precipitation fractions +
	1410	! temperature and water species tendencies
	1411	DO i = 1, klon
	1412	psfl(i,k)=zifl(i)
	1413	prfl(i,k)=zrfl(i)
	1414	pfraclr(i,k)=znebprecipclr(i)
	1415	pfracld(i,k)=znebprecipcld(i)
	1416	d_ql(i,k) = (1-zfice(i))*zoliq(i)
	1417	d_qi(i,k) = zfice(i)*zoliq(i)
	1418	d_q(i,k) = zq(i) - qt(i,k)
	1419	! c_iso: same for isotopes
	1420	d_t(i,k) = zt(i) - temp(i,k)
	1421	ENDDO
	1422
	1423
[4226]	1424	ENDDO
[4059]	1425
[3999]	1426
	1427	! Rain or snow at the surface (depending on the first layer temperature)
	1428	DO i = 1, klon
	1429	snow(i) = zifl(i)
	1430	rain(i) = zrfl(i)
[4392]	1431	! c_iso final output
[3999]	1432	ENDDO
	1433
	1434	IF (ncoreczq>0) THEN
	1435	WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
	1436	ENDIF
	1437
[4654]	1438
	1439	RETURN
	1440
[4380]	1441	END SUBROUTINE lscp
[3999]	1442	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1443
[4664]	1444	END MODULE lmdz_lscp

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

Download in other formats:

Original Format