Context Navigation

lmdz_lscp.F90 @ 5501

Last change on this file since 5501 was 5227, checked in by abarral, 4 months ago
Merge r5210 5211
File size: 69.6 KB

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

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

Download in other formats:

Original Format