Context Navigation

lmdz_lscp.F90 @ 4950

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

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

Download in other formats:

Original Format