Context Navigation

lmdz_lscp.f90 @ 5602

Last change on this file since 5602 was 5602, checked in by aborella, 2 months ago
Removed deep convection clouds from prognostic cloud properties
File size: 57.1 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, &
[5396]	9	paprs, pplay, omega, temp, qt, ql_seri, qi_seri, &
[5602]	10	ptconv, cldfracv, qcondcv, ratqs, sigma_qtherm, &
[5204]	11	d_t, d_q, d_ql, d_qi, rneb, rneblsvol, &
[5007]	12	pfraclr, pfracld, &
	13	cldfraliq, sigma2_icefracturb,mean_icefracturb, &
[4530]	14	radocond, radicefrac, rain, snow, &
[3999]	15	frac_impa, frac_nucl, beta, &
[5007]	16	prfl, psfl, rhcl, qta, fraca, &
	17	tv, pspsk, tla, thl, iflag_cld_th, &
[5204]	18	iflag_ice_thermo, distcltop, temp_cltop, &
	19	tke, tke_dissip, &
[5575]	20	cell_area, stratomask, &
[5601]	21	cf_seri, qvc_seri, u_seri, v_seri, &
[5204]	22	qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
[5456]	23	dcf_sub, dcf_con, dcf_mix, &
[5204]	24	dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, &
	25	dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati, &
[5601]	26	cfa_seri, qta_seri, flight_dist, flight_h2o, &
	27	qice_cont, Tcritcont, &
	28	qcritcont, potcontfraP, potcontfraNP, &
	29	cloudth_sth, &
	30	cloudth_senv, cloudth_sigmath, cloudth_sigmaenv, &
[4830]	31	qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, &
	32	dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
[5589]	33	dqsmelt, dqsfreez, dqised, dcfsed, dqvcsed)
[3999]	34
	35	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	36	! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
	37	! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
	38	!--------------------------------------------------------------------------------
	39	! Date: 01/2021
	40	!--------------------------------------------------------------------------------
	41	! Aim: Large Scale Clouds and Precipitation (LSCP)
[4226]	42	!
[3999]	43	! This code is a new version of the fisrtilp.F90 routine, which is itself a
[4072]	44	! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
[3999]	45	! and 'ilp' (il pleut, L. Li)
	46	!
[4380]	47	! Compared to the original fisrtilp code, lscp
[3999]	48	! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
	49	! -> consider always precipitation thermalisation (fl_cor_ebil>0)
	50	! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
	51	! -> option "oldbug" by JYG has been removed
	52	! -> iflag_t_glace >0 always
	53	! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
	54	! -> rectangular distribution from L. Li no longer available
	55	! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
	56	!--------------------------------------------------------------------------------
	57	! References:
	58	!
[4412]	59	! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
[3999]	60	! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
	61	! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
[4412]	62	! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
[3999]	63	! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
[4412]	64	! - Touzze-Peifert Ludo, PhD thesis, p117-124
[3999]	65	! -------------------------------------------------------------------------------
	66	! Code structure:
	67	!
[4226]	68	! P0> Thermalization of the precipitation coming from the overlying layer
[3999]	69	! P1> Evaporation of the precipitation (falling from the k+1 level)
	70	! P2> Cloud formation (at the k level)
[4412]	71	! P2.A.1> With the PDFs, calculation of cloud properties using the inital
[3999]	72	! values of T and Q
	73	! P2.A.2> Coupling between condensed water and temperature
	74	! P2.A.3> Calculation of final quantities associated with cloud formation
[4412]	75	! P2.B> Release of Latent heat after cloud formation
[3999]	76	! P3> Autoconversion to precipitation (k-level)
	77	! P4> Wet scavenging
	78	!------------------------------------------------------------------------------
	79	! Some preliminary comments (JBM) :
	80	!
	81	! The cloud water that the radiation scheme sees is not the same that the cloud
	82	! water used in the physics and the dynamics
	83	!
[4412]	84	! During the autoconversion to precipitation (P3 step), radocond (cloud water used
[3999]	85	! by the radiation scheme) is calculated as an average of the water that remains
	86	! in the cloud during the precipitation and not the water remaining at the end
	87	! of the time step. The latter is used in the rest of the physics and advected
	88	! by the dynamics.
	89	!
	90	! In summary:
	91	!
	92	! Radiation:
	93	! xflwc(newmicro)+xfiwc(newmicro) =
[4412]	94	! radocond=lwcon(nc)+iwcon(nc)
[3999]	95	!
	96	! Notetheless, be aware of:
	97	!
[4412]	98	! radocond .NE. ocond(nc)
[3999]	99	! i.e.:
	100	! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
[4412]	101	! but oliq+(ocond-oliq) .EQ. ocond
[3999]	102	! (which is not trivial)
	103	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4226]	104	!
[4654]	105
	106	! USE de modules contenant des fonctions.
[4651]	107	USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
[5007]	108	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
[5413]	109	USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
[5241]	110	USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
[5413]	111	USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
	112	USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
[3999]	113
[4664]	114	! Use du module lmdz_lscp_ini contenant les constantes
[5204]	115	USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
[5413]	116	USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
	117	USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
	118	USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
	119	USE lmdz_lscp_ini, ONLY : t_glace_min, min_frac_th_cld
	120	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
[5412]	121	USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
[5211]	122	USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
[5452]	123	USE lmdz_lscp_ini, ONLY : ok_plane_contrail
[5204]	124
[5575]	125	! Temporary call for Lamquin et al (2012) diagnostics
	126	USE phys_local_var_mod, ONLY : issrfra100to150, issrfra150to200, issrfra200to250
	127	USE phys_local_var_mod, ONLY : issrfra250to300, issrfra300to400, issrfra400to500
[5601]	128	USE phys_local_var_mod, ONLY : dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub
	129	USE phys_local_var_mod, ONLY : dcfa_cir, dqta_cir, dcfa_mix, dqia_mix, dqta_mix
[5573]	130
[3999]	131	IMPLICIT NONE
	132
	133	!===============================================================================
[4380]	134	! VARIABLES DECLARATION
[3999]	135	!===============================================================================
	136
	137	! INPUT VARIABLES:
	138	!-----------------
	139
[4380]	140	INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels
[3999]	141	REAL, INTENT(IN) :: dtime ! time step [s]
[4380]	142	REAL, INTENT(IN) :: missing_val ! missing value for output
[4059]	143
[3999]	144	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
	145	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
[4654]	146	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K)
[5383]	147	REAL, DIMENSION(klon,klev), INTENT(IN) :: omega ! vertical velocity [Pa/s]
[4686]	148	REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg]
[5396]	149	REAL, DIMENSION(klon,klev), INTENT(IN) :: ql_seri ! liquid specific content [kg/kg]
	150	REAL, DIMENSION(klon,klev), INTENT(IN) :: qi_seri ! ice specific content [kg/kg]
[3999]	151	INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds
	152	INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
[4226]	153	! CR: if iflag_ice_thermo=2, only convection is active
[3999]	154	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
[5602]	155	REAL, DIMENSION(klon,klev), INTENT(IN) :: cldfracv ! cloud fraction from deep convection [-]
	156	REAL, DIMENSION(klon,klev), INTENT(IN) :: qcondcv ! in-cloud condensed specific humidity from deep convection [kg/kg]
[3999]	157
	158	!Inputs associated with thermal plumes
[4226]	159
[5007]	160	REAL, DIMENSION(klon,klev), INTENT(IN) :: tv ! virtual potential temperature [K]
	161	REAL, DIMENSION(klon,klev), INTENT(IN) :: qta ! specific humidity within thermals [kg/kg]
	162	REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-]
	163	REAL, DIMENSION(klon,klev), INTENT(IN) :: pspsk ! exner potential (p/100000)**(R/cp)
	164	REAL, DIMENSION(klon,klev), INTENT(IN) :: tla ! liquid temperature within thermals [K]
[5413]	165	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke !--turbulent kinetic energy [m2/s2]
	166	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip !--TKE dissipation [m2/s3]
[3999]	167
	168	! INPUT/OUTPUT variables
	169	!------------------------
[4562]	170
[5208]	171	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: thl ! liquid potential temperature [K]
	172	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs,sigma_qtherm ! function of pressure that sets the large-scale
[4059]	173
[5208]	174
[5204]	175	! INPUT/OUTPUT condensation and ice supersaturation
	176	!--------------------------------------------------
	177	REAL, DIMENSION(klon,klev), INTENT(INOUT):: cf_seri ! cloud fraction [-]
[5601]	178	REAL, DIMENSION(klon,klev), INTENT(INOUT):: qvc_seri ! cloudy water vapor [kg/kg]
[5204]	179	REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri ! eastward wind [m/s]
	180	REAL, DIMENSION(klon,klev), INTENT(IN) :: v_seri ! northward wind [m/s]
	181	REAL, DIMENSION(klon), INTENT(IN) :: cell_area ! area of each cell [m2]
[5575]	182	REAL, DIMENSION(klon,klev), INTENT(IN) :: stratomask ! fraction of stratosphere (0 or 1)
[4380]	183
[5204]	184	! INPUT/OUTPUT aviation
	185	!--------------------------------------------------
[5601]	186	REAL, DIMENSION(klon,klev), INTENT(INOUT):: cfa_seri ! linear contrails fraction [-]
	187	REAL, DIMENSION(klon,klev), INTENT(INOUT):: qta_seri ! linear contrails total specific humidity [kg/kg]
[5575]	188	REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_dist ! aviation distance flown within the mesh [m/s/mesh]
	189	REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_h2o ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
[3999]	190
	191	! OUTPUT variables
	192	!-----------------
	193
[5204]	194	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K]
	195	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg]
	196	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg]
	197	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg]
	198	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-]
	199	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-]
	200	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-]
	201	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-]
[5007]	202	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliq ! liquid fraction of cloud [-]
	203	REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturb ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
	204	REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturb ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
[5204]	205	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg]
	206	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme
	207	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-]
	208	REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2]
	209	REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2]
	210	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
	211	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
	212	REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m]
	213	REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K]
	214	REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water
[4686]	215
[5007]	216	! fraction of aerosol scavenging through impaction and nucleation (for on-line)
[3999]	217
[5007]	218	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-]
	219	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
[3999]	220
[5204]	221	! for condensation and ice supersaturation
[3999]	222
[5204]	223	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsub !--specific total water content in sub-saturated clear sky region [kg/kg]
	224	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qissr !--specific total water content in supersat region [kg/kg]
	225	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld !--specific total water content in cloudy region [kg/kg]
	226	REAL, DIMENSION(klon,klev), INTENT(OUT) :: subfra !--mesh fraction of subsaturated clear sky [-]
	227	REAL, DIMENSION(klon,klev), INTENT(OUT) :: issrfra !--mesh fraction of ISSR [-]
	228	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_cond !--coefficient governing the ice nucleation RHi threshold [-]
	229	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_sub !--cloud fraction tendency because of sublimation [s-1]
	230	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_con !--cloud fraction tendency because of condensation [s-1]
	231	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_mix !--cloud fraction tendency because of cloud mixing [s-1]
	232	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_adj !--specific ice content tendency because of temperature adjustment [kg/kg/s]
	233	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_sub !--specific ice content tendency because of sublimation [kg/kg/s]
	234	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_con !--specific ice content tendency because of condensation [kg/kg/s]
	235	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_mix !--specific ice content tendency because of cloud mixing [kg/kg/s]
	236	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_adj !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
	237	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_sub !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
	238	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_con !--specific cloud water vapor tendency because of condensation [kg/kg/s]
	239	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_mix !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
	240	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl !--saturation specific humidity wrt liquid [kg/kg]
	241	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsati !--saturation specific humidity wrt ice [kg/kg]
	242
	243	! for contrails and aviation
	244
[5601]	245	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qice_cont !--condensed water in linear contrails used in the radiation scheme [kg/kg]
[5456]	246	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcritcont !--critical temperature for contrail formation [K]
	247	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcritcont !--critical specific humidity for contrail formation [kg/kg]
	248	REAL, DIMENSION(klon,klev), INTENT(OUT) :: potcontfraP !--potential persistent contrail fraction [-]
	249	REAL, DIMENSION(klon,klev), INTENT(OUT) :: potcontfraNP !--potential non-persistent contrail fraction [-]
[5204]	250
	251
[4803]	252	! for POPRECIP
[4686]	253
[4830]	254	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg]
	255	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg]
[4819]	256	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqreva !--rain tendendy due to evaporation [kg/kg/s]
	257	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
	258	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
	259	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
	260	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
	261	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
	262	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsrim !--snow tendency due to riming [kg/kg/s]
	263	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsmelt !--snow tendency due to melting [kg/kg/s]
	264	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrmelt !--rain tendency due to melting [kg/kg/s]
	265	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
	266	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
[5589]	267	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqised !--ice water content tendency due to sedmentation of ice crystals [kg/kg/s]
	268	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcfsed !--cloud fraction tendency due to sedimentation of ice crystals [kg/kg/s]
	269	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvcsed !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
[4686]	270
[5204]	271	! for thermals
[4803]	272
[5204]	273	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth !--mean saturation deficit in thermals
	274	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv !--mean saturation deficit in environment
	275	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath !--std of saturation deficit in thermals
	276	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv !--std of saturation deficit in environment
	277
	278
[3999]	279	! LOCAL VARIABLES:
	280	!----------------
[5601]	281	REAL, DIMENSION(klon) :: qliq_in, qice_in
[4654]	282	REAL, DIMENSION(klon,klev) :: ctot
	283	REAL, DIMENSION(klon,klev) :: ctot_vol
[5450]	284	REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
[5413]	285	REAL :: zdelta
[4654]	286	REAL, DIMENSION(klon) :: zdqsdT_raw
[5007]	287	REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T
	288	REAL, DIMENSION(klon) :: Tbef,qlbef,DT ! temperature, humidity and temp. variation during lognormal iteration
[4654]	289	REAL :: num,denom
	290	REAL :: cste
[4910]	291	REAL, DIMENSION(klon) :: zfice_th
[4654]	292	REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
[5413]	293	REAL, DIMENSION(klon) :: zrfl, zifl
[5450]	294	REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn, zqnl
[4654]	295	REAL, DIMENSION(klon) :: zoliql, zoliqi
	296	REAL, DIMENSION(klon) :: zt
[5450]	297	REAL, DIMENSION(klon) :: zfice,zneb, znebl
[4654]	298	REAL, DIMENSION(klon) :: dzfice
[5007]	299	REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
[4654]	300	REAL, DIMENSION(klon) :: qtot, qzero
[3999]	301	! Variables precipitation energy conservation
[4654]	302	REAL, DIMENSION(klon) :: zmqc
	303	REAL :: zalpha_tr
	304	REAL :: zfrac_lessi
	305	REAL, DIMENSION(klon) :: zprec_cond
	306	REAL, DIMENSION(klon) :: zlh_solid
[4803]	307	REAL, DIMENSION(klon) :: ztupnew
[4913]	308	REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
[5589]	309	REAL, DIMENSION(klon) :: zqice_sedim ! for sedimentation of ice crystals
[4654]	310	REAL, DIMENSION(klon) :: zrflclr, zrflcld
	311	REAL, DIMENSION(klon) :: ziflclr, ziflcld
[5411]	312	REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
[5413]	313	REAL, DIMENSION(klon) :: tot_zneb
[4654]	314	REAL :: qlmpc, qimpc, rnebmpc
[5007]	315	REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
	316	REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl ! for icefrac_lscp_turb
[5413]	317
	318	! for quantity of condensates seen by radiation
	319	REAL, DIMENSION(klon) :: zradocond, zradoice
	320	REAL, DIMENSION(klon) :: zrho_up, zvelo_up
[5204]	321
	322	! for condensation and ice supersaturation
[5450]	323	REAL, DIMENSION(klon) :: qvc, qvcl, shear
[5204]	324	REAL :: delta_z
[5452]	325	! for contrails
[5601]	326	REAL, DIMENSION(klon) :: contfra, qcont
[5452]	327	!--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
[5204]	328	! Constants used for calculating ratios that are advected (using a parent-child
	329	! formalism). This is not done in the dynamical core because at this moment,
	330	! only isotopes can use this parent-child formalism. Note that the two constants
	331	! are the same as the one use in the dynamical core, being also defined in
	332	! dyn3d_common/infotrac.F90
	333	REAL :: min_qParent, min_ratio
[5575]	334	!--for Lamquin et al 2012 diagnostics
	335	REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP
	336	REAL, DIMENSION(klon) :: issrfra250to300UP, issrfra300to400UP, issrfra400to500UP
[3999]	337
[5413]	338	INTEGER i, k, kk, iter
[4654]	339	INTEGER, DIMENSION(klon) :: n_i
[4226]	340	INTEGER ncoreczq
[4654]	341	INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
[4803]	342	LOGICAL iftop
[3999]	343
[4654]	344	LOGICAL, DIMENSION(klon) :: lognormale
	345	LOGICAL, DIMENSION(klon) :: keepgoing
[3999]	346
	347	CHARACTER (len = 20) :: modname = 'lscp'
	348	CHARACTER (len = 80) :: abort_message
	349
	350
	351	!===============================================================================
	352	! INITIALISATION
	353	!===============================================================================
	354
[4226]	355	! Few initial checks
[3999]	356
[4654]	357
[3999]	358	IF (iflag_fisrtilp_qsat .LT. 0) THEN
[4226]	359	abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
	360	CALL abort_physic(modname,abort_message,1)
[3999]	361	ENDIF
	362
	363	! Few initialisations
	364
	365	ctot_vol(1:klon,1:klev)=0.0
	366	rneblsvol(1:klon,1:klev)=0.0
[5411]	367	znebprecip(:)=0.0
[4226]	368	znebprecipclr(:)=0.0
	369	znebprecipcld(:)=0.0
[3999]	370	mpc_bl_points(:,:)=0
	371
	372	IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
	373
	374	! AA for 'safety' reasons
	375	zalpha_tr = 0.
	376	zfrac_lessi = 0.
[4380]	377	beta(:,:)= 0.
[3999]	378
[4380]	379	! Initialisation of variables:
	380
[3999]	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
[4530]	388	pfraclr(:,:)=0.0
	389	pfracld(:,:)=0.0
[5007]	390	cldfraliq(:,:)=0.
	391	sigma2_icefracturb(:,:)=0.
	392	mean_icefracturb(:,:)=0.
[4412]	393	radocond(:,:) = 0.0
[3999]	394	radicefrac(:,:) = 0.0
[4226]	395	frac_nucl(:,:) = 1.0
	396	frac_impa(:,:) = 1.0
[3999]	397	rain(:) = 0.0
	398	snow(:) = 0.0
[5450]	399	zfice(:)=1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
[4114]	400	dzfice(:)=0.0
[5007]	401	zfice_turb(:)=0.0
	402	dzfice_turb(:)=0.0
[3999]	403	zrfl(:) = 0.0
	404	zifl(:) = 0.0
	405	zneb(:) = seuil_neb
[4226]	406	zrflclr(:) = 0.0
	407	ziflclr(:) = 0.0
	408	zrflcld(:) = 0.0
	409	ziflcld(:) = 0.0
	410	tot_zneb(:) = 0.0
	411	qzero(:) = 0.0
[5007]	412	zdistcltop(:)=0.0
	413	ztemp_cltop(:) = 0.0
[4803]	414	ztupnew(:)=0.0
[5204]	415
[4654]	416	distcltop(:,:)=0.
	417	temp_cltop(:,:)=0.
[5007]	418
[5204]	419	!--Ice supersaturation
	420	gamma_cond(:,:) = 1.
	421	qissr(:,:) = 0.
	422	issrfra(:,:) = 0.
	423	dcf_sub(:,:) = 0.
	424	dcf_con(:,:) = 0.
	425	dcf_mix(:,:) = 0.
	426	dqi_adj(:,:) = 0.
	427	dqi_sub(:,:) = 0.
	428	dqi_con(:,:) = 0.
	429	dqi_mix(:,:) = 0.
	430	dqvc_adj(:,:) = 0.
	431	dqvc_sub(:,:) = 0.
	432	dqvc_con(:,:) = 0.
	433	dqvc_mix(:,:) = 0.
	434	qvc(:) = 0.
	435	shear(:) = 0.
	436	min_qParent = 1.e-30
	437	min_ratio = 1.e-16
	438
[5575]	439	!--for Lamquin et al (2012) diagnostics
	440	issrfra100to150(:) = 0.
	441	issrfra100to150UP(:) = 0.
	442	issrfra150to200(:) = 0.
	443	issrfra150to200UP(:) = 0.
	444	issrfra200to250(:) = 0.
	445	issrfra200to250UP(:) = 0.
	446	issrfra250to300(:) = 0.
	447	issrfra250to300UP(:) = 0.
	448	issrfra300to400(:) = 0.
	449	issrfra300to400UP(:) = 0.
	450	issrfra400to500(:) = 0.
	451	issrfra400to500UP(:) = 0.
	452
[4803]	453	!-- poprecip
[4830]	454	qraindiag(:,:)= 0.
	455	qsnowdiag(:,:)= 0.
[4819]	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.
[4913]	467	zqupnew(:) = 0.
	468	zqvapclr(:) = 0.
[5589]	469	zqice_sedim(:)= 0.
[3999]	470
[4803]	471
	472
[4392]	473	!c_iso: variable initialisation for iso
	474
[3999]	475	!===============================================================================
	476	! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
	477	!===============================================================================
	478
	479	ncoreczq=0
	480
	481	DO k = klev, 1, -1
	482
[4803]	483	IF (k.LE.klev-1) THEN
	484	iftop=.false.
	485	ELSE
	486	iftop=.true.
	487	ENDIF
	488
[3999]	489	! Initialisation temperature and specific humidity
[5007]	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
[3999]	495	DO i = 1, klon
[4654]	496	zt(i)=temp(i,k)
[4686]	497	zq(i)=qt(i,k)
[5601]	498	qliq_in(i) = ql_seri(i,k)
	499	qice_in(i) = qi_seri(i,k)
[4803]	500	IF (.not. iftop) THEN
[4913]	501	ztupnew(i) = temp(i,k+1) + d_t(i,k+1)
	502	zqupnew(i) = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
	503	!--zqs(i) is the saturation specific humidity in the layer above
	504	zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
[4803]	505	ENDIF
[4392]	506	!c_iso init of iso
[3999]	507	ENDDO
[5408]	508
	509	! --------------------------------------------------------------------
	510	! P1> Precipitation processes, before cloud formation:
	511	! Thermalization of precipitation falling from the overlying layer AND
	512	! Precipitation evaporation/sublimation/melting
	513	!---------------------------------------------------------------------
[4803]	514
	515	!================================================================
	516	! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
[5431]	517	IF ( ok_poprecip ) THEN
[4803]	518
[5408]	519	CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	520	zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
[5589]	521	zqvapclr, zqupnew, zqice_sedim, &
[5601]	522	cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
[5408]	523	zrfl, zrflclr, zrflcld, &
	524	zifl, ziflclr, ziflcld, &
[5579]	525	dqreva(:,k), dqssub(:,k), &
[5589]	526	dcfsed(:,k), dqvcsed(:,k) &
[5408]	527	)
[4803]	528
	529	!================================================================
	530	ELSE
	531
[5408]	532	CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
[5411]	533	zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
[5408]	534	zrfl, zrflclr, zrflcld, &
	535	zifl, ziflclr, ziflcld, &
	536	dqreva(:,k), dqssub(:,k) &
	537	)
[3999]	538
[4803]	539	ENDIF ! (ok_poprecip)
[3999]	540
	541	! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
	542	!-------------------------------------------------------
[4226]	543
[4072]	544	qtot(:)=zq(:)+zmqc(:)
[5383]	545	CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
[4072]	546	DO i = 1, klon
[3999]	547	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
[4059]	548	zdqsdT_raw(i) = zdqs(i)RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTT*zdelta)
[3999]	549	IF (zq(i) .LT. 1.e-15) THEN
	550	ncoreczq=ncoreczq+1
	551	zq(i)=1.e-15
	552	ENDIF
[4392]	553	! c_iso: do something similar for isotopes
[3999]	554
[4072]	555	ENDDO
[3999]	556
	557	! --------------------------------------------------------------------
[5408]	558	! P2> Cloud formation
[3999]	559	!---------------------------------------------------------------------
	560	!
	561	! Unlike fisrtilp, we always assume a 'fractional cloud' approach
	562	! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
	563	! is prescribed and depends on large scale variables and boundary layer
	564	! properties)
	565	! The decrease in condensed part due tu latent heating is taken into
	566	! account
	567	! -------------------------------------------------------------------
	568
[5408]	569	! P2.1> With the PDFs (log-normal, bigaussian)
[4424]	570	! cloud properties calculation with the initial values of t and q
[3999]	571	! ----------------------------------------------------------------
	572
	573	! initialise gammasat and qincloud_mpc
	574	gammasat(:)=1.
	575	qincloud_mpc(:)=0.
	576
	577	IF (iflag_cld_th.GE.5) THEN
[4910]	578	! Cloud cover and content in meshes affected by shallow convection,
	579	! are retrieved from a bi-gaussian distribution of the saturation deficit
	580	! following Jam et al. 2013
[3999]	581
[4910]	582	IF (iflag_cloudth_vert.LE.2) THEN
	583	! Old version of Arnaud Jam
[3999]	584
[4910]	585	CALL cloudth(klon,klev,k,tv, &
	586	zq,qta,fraca, &
	587	qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
[4654]	588	ratqs,zqs,temp, &
[4651]	589	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	590
[4651]	591
[3999]	592	ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
[4910]	593	! Default version of Arnaud Jam
	594
	595	CALL cloudth_v3(klon,klev,k,tv, &
	596	zq,qta,fraca, &
	597	qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
[5208]	598	ratqs,sigma_qtherm,zqs,temp, &
[4651]	599	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	600
	601
	602	ELSEIF (iflag_cloudth_vert.EQ.6) THEN
[4910]	603	! Jean Jouhaud's version, with specific separation between surface and volume
	604	! cloud fraction Decembre 2018
[3999]	605
[4910]	606	CALL cloudth_v6(klon,klev,k,tv, &
	607	zq,qta,fraca, &
	608	qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
[4654]	609	ratqs,zqs,temp, &
[4651]	610	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	611
[4910]	612	ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
	613	! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
[5007]	614	! for boundary-layer mixed phase clouds
[4910]	615	CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
	616	pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
	617	ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
	618	cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
	619
[3999]	620	ENDIF
	621
	622
	623	DO i=1,klon
	624	rneb(i,k)=ctot(i,k)
	625	rneblsvol(i,k)=ctot_vol(i,k)
	626	zqn(i)=qcloud(i)
[5204]	627	!--AB grid-mean vapor in the cloud - we assume saturation adjustment
	628	qvc(i) = rneb(i,k) * zqs(i)
[3999]	629	ENDDO
	630
	631	ENDIF
	632
	633	IF (iflag_cld_th .LE. 4) THEN
	634
	635	! lognormal
[5007]	636	lognormale(:) = .TRUE.
[3999]	637
	638	ELSEIF (iflag_cld_th .GE. 6) THEN
	639
	640	! lognormal distribution when no thermals
[5007]	641	lognormale(:) = fraca(:,k) < min_frac_th_cld
[3999]	642
	643	ELSE
[4226]	644	! When iflag_cld_th=5, we always assume
[3999]	645	! bi-gaussian distribution
[5007]	646	lognormale(:) = .FALSE.
[3999]	647
	648	ENDIF
	649
	650	DT(:) = 0.
	651	n_i(:)=0
	652	Tbef(:)=zt(:)
	653	qlbef(:)=0.
	654
	655	! Treatment of non-boundary layer clouds (lognormale)
[5383]	656	! We iterate here to take into account the change in qsat(T) and ice fraction
	657	! during the condensation process
	658	! the increment in temperature is calculated using the first principle of
	659	! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
	660	! and a clear sky fraction)
	661	! note that we assume that the vapor in the cloud is at saturation for this calculation
[4226]	662
[3999]	663	DO iter=1,iflag_fisrtilp_qsat+1
[4226]	664
[5204]	665	keepgoing(:) = .FALSE.
	666
[3999]	667	DO i=1,klon
	668
[4424]	669	! keepgoing = .true. while convergence is not satisfied
[4226]	670
[5204]	671	IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
[4226]	672
[3999]	673	! if not convergence:
[5204]	674	! we calculate a new iteration
	675	keepgoing(i) = .TRUE.
[3999]	676
	677	! P2.2.1> cloud fraction and condensed water mass calculation
	678	! Calculated variables:
	679	! rneb : cloud fraction
	680	! zqn : total water within the cloud
	681	! zcond: mean condensed water within the mesh
	682	! rhcl: clear-sky relative humidity
	683	!---------------------------------------------------------------
	684
[4424]	685	! new temperature that only serves in the iteration process:
[3999]	686	Tbef(i)=Tbef(i)+DT(i)
	687
	688	! Rneb, qzn and zcond for lognormal PDFs
[4072]	689	qtot(i)=zq(i)+zmqc(i)
[4226]	690
[4072]	691	ENDIF
[3999]	692
[4072]	693	ENDDO
	694
[4380]	695	! Calculation of saturation specific humidity and ice fraction
[5383]	696	CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
[5450]	697
	698	IF (iflag_icefrac .GE. 3) THEN
	699	! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
	700	! and cloud condensed water content. idea following Dietlitcher et al. 2018, GMD
	701	! we make this option works only for the physically-based and tke-depdenent param (iflag_icefrac>=1)
	702	DO i=1,klon
	703	CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,1,.false.,zqsl,zdqsl)
	704	CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,2,.false.,zqsi,zdqsi)
	705	zqs(i)=zfice(i)zqsi(i)+(1.-zfice(i))zqsl(i)
	706	zdqs(i)=zfice(i)zdqsi(i)+zqsi(i)dzfice(i)+(1.-zfice(i))zdqsl(i)-zqsl(i)dzfice(i)
	707	ENDDO
	708	ENDIF
	709
[5383]	710	CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
[4072]	711	! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
	712	zdqs(:) = gammasat(:)zdqs(:)+zqs(:)dgammasatdt(:)
	713
[5450]	714	! Cloud condensation based on subgrid pdf
[5204]	715	!--AB Activates a condensation scheme that allows for
	716	!--ice supersaturation and contrails evolution from aviation
	717	IF (ok_ice_supersat) THEN
	718
	719	!--Calculate the shear value (input for condensation and ice supersat)
	720	DO i = 1, klon
	721	!--Cell thickness [m]
	722	delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
	723	IF ( iftop ) THEN
	724	! top
	725	shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
	726	+ ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
	727	ELSEIF ( k .EQ. 1 ) THEN
	728	! surface
	729	shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
	730	+ ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
	731	ELSE
	732	! other layers
	733	shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
	734	- (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
	735	+ ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
	736	- (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
	737	ENDIF
	738	ENDDO
	739
	740	!---------------------------------------------
	741	!-- CONDENSATION AND ICE SUPERSATURATION --
	742	!---------------------------------------------
	743
	744	CALL condensation_ice_supersat( &
[5602]	745	klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), &
	746	cldfracv(:,k), qcondcv(:,k), &
[5601]	747	cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
[5575]	748	shear, tke_dissip(:,k), cell_area, stratomask(:,k), &
[5406]	749	Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
	750	rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
[5204]	751	dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
	752	dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
	753	dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
[5601]	754	cfa_seri(:,k), qta_seri(:,k), flight_dist(:,k), flight_h2o(:,k), &
	755	contfra, qcont, Tcritcont(:,k), qcritcont(:,k), &
	756	potcontfraP(:,k), potcontfraNP(:,k), &
	757	dcfa_ini(:,k), dqia_ini(:,k), dqta_ini(:,k), &
	758	dcfa_sub(:,k), dqia_sub(:,k), dqta_sub(:,k), &
	759	dcfa_cir(:,k), dqta_cir(:,k), &
	760	dcfa_mix(:,k), dqia_mix(:,k), dqta_mix(:,k))
[5204]	761
	762
[5211]	763	ELSE
	764	!--generalised lognormal condensation scheme (Bony and Emanuel 2001)
	765
[5241]	766	CALL condensation_lognormal( &
	767	klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
	768	keepgoing, rneb(:,k), zqn, qvc)
[5204]	769
	770
	771	ENDIF ! .NOT. ok_ice_supersat
[4226]	772
[5450]	773	! cloud phase determination
	774	IF (iflag_icefrac .GE. 2) THEN
	775	! phase partitioning depending on temperature. activates here in the iteration process if iflag_icefrac > 2
[5579]	776	CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
[5450]	777	rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
	778	ELSE
	779	! phase partitioning depending on temperature and eventually distance to cloud top
	780	IF (iflag_t_glace.GE.4) THEN
	781	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
	782	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
	783	ENDIF
	784	CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
	785	ENDIF
	786
	787
[5204]	788	DO i=1,klon
	789	IF (keepgoing(i)) THEN
[3999]	790
	791	! If vertical heterogeneity, change fraction by volume as well
	792	IF (iflag_cloudth_vert.GE.3) THEN
	793	ctot_vol(i,k)=rneb(i,k)
	794	rneblsvol(i,k)=ctot_vol(i,k)
	795	ENDIF
	796
	797
[4072]	798	! P2.2.2> Approximative calculation of temperature variation DT
	799	! due to condensation.
	800	! Calculated variables:
	801	! dT : temperature change due to condensation
	802	!---------------------------------------------------------------
[3999]	803
	804
	805	IF (zfice(i).LT.1) THEN
	806	cste=RLVTT
	807	ELSE
	808	cste=RLSTT
	809	ENDIF
[5007]	810
[5204]	811	IF ( ok_unadjusted_clouds ) THEN
	812	!--AB We relax the saturation adjustment assumption
	813	!-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
	814	IF ( rneb(i,k) .GT. eps ) THEN
	815	qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
	816	ELSE
	817	qlbef(i) = 0.
	818	ENDIF
	819	ELSE
	820	qlbef(i)=max(0.,zqn(i)-zqs(i))
	821	ENDIF
	822
[3999]	823	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
[4226]	824	+zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
[3999]	825	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
	826	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
[4226]	827	qlbef(i)dzfice(i)
[4424]	828	! here we update a provisory temperature variable that only serves in the iteration
	829	! process
[3999]	830	DT(i)=num/denom
	831	n_i(i)=n_i(i)+1
	832
[4424]	833	ENDIF ! end keepgoing
[3999]	834
	835	ENDDO ! end loop on i
	836
	837	ENDDO ! iter=1,iflag_fisrtilp_qsat+1
	838
[5408]	839	! P2.2> Final quantities calculation
[3999]	840	! Calculated variables:
	841	! rneb : cloud fraction
	842	! zcond: mean condensed water in the mesh
	843	! zqn : mean water vapor in the mesh
[4562]	844	! zfice: ice fraction in clouds
[3999]	845	! zt : temperature
	846	! rhcl : clear-sky relative humidity
	847	! ----------------------------------------------------------------
	848
	849
[5450]	850	! Cloud phase final determination
	851	!--------------------------------
[4562]	852	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
	853	IF (iflag_t_glace.GE.4) THEN
[5007]	854	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
	855	distcltop(:,k)=zdistcltop(:)
	856	temp_cltop(:,k)=ztemp_cltop(:)
	857	ENDIF
[5450]	858	! Partition function depending on temperature for all clouds (shallow convective and stratiform)
[5204]	859	CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
[4562]	860
[5450]	861	! Partition function depending on tke for non shallow-convective clouds, erase previous estimation
[5007]	862	IF (iflag_icefrac .GE. 1) THEN
[5579]	863	CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
[5007]	864	rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
	865	ENDIF
	866
[5450]	867	! Water vapor update, subsequent latent heat exchange for each cloud type
	868	!------------------------------------------------------------------------
[3999]	869	DO i=1, klon
[5007]	870	! Overwrite phase partitioning in boundary layer mixed phase clouds when the
	871	! iflag_cloudth_vert=7 and specific param is activated
[3999]	872	IF (mpc_bl_points(i,k) .GT. 0) THEN
	873	zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
	874	! following line is very strange and probably wrong
	875	rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
	876	! water vapor update and partition function if thermals
	877	zq(i) = zq(i) - zcond(i)
[4910]	878	zfice(i)=zfice_th(i)
[3999]	879	ELSE
	880	! Checks on rneb, rhcl and zqn
	881	IF (rneb(i,k) .LE. 0.0) THEN
	882	zqn(i) = 0.0
	883	rneb(i,k) = 0.0
	884	zcond(i) = 0.0
	885	rhcl(i,k)=zq(i)/zqs(i)
	886	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	887	zqn(i) = zq(i)
	888	rneb(i,k) = 1.0
[5204]	889	IF ( ok_unadjusted_clouds ) THEN
	890	!--AB We relax the saturation adjustment assumption
	891	!-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
	892	zcond(i) = MAX(0., zqn(i) - qvc(i))
	893	ELSE
	894	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
	895	ENDIF
[3999]	896	rhcl(i,k)=1.0
	897	ELSE
[5204]	898	IF ( ok_unadjusted_clouds ) THEN
	899	!--AB We relax the saturation adjustment assumption
	900	!-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
	901	zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
	902	ELSE
	903	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
	904	ENDIF
[3999]	905	! following line is very strange and probably wrong:
	906	rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
[5007]	907	! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
	908	IF (iflag_icefrac .GE. 1) THEN
	909	IF (lognormale(i)) THEN
	910	zcond(i) = zqliq(i) + zqice(i)
[5450]	911	zfice(i) = zfice_turb(i)
[5007]	912	rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
	913	ENDIF
	914	ENDIF
[3999]	915	ENDIF
	916
[4226]	917	! water vapor update
	918	zq(i) = zq(i) - zcond(i)
[3999]	919
	920	ENDIF
[5450]	921
	922
[3999]	923	! temperature update due to phase change
	924	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
	925	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
	926	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
	927	ENDDO
	928
	929	! If vertical heterogeneity, change volume fraction
	930	IF (iflag_cloudth_vert .GE. 3) THEN
	931	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
	932	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
	933	ENDIF
	934
[5204]	935
[3999]	936	! ----------------------------------------------------------------
[5408]	937	! P3> Precipitation processes, after cloud formation
	938	! - precipitation formation, melting/freezing
[3999]	939	! ----------------------------------------------------------------
[4226]	940
[5413]	941	! Initiate the quantity of liquid and solid condensates
	942	! Note that in the following, zcond is the total amount of condensates
	943	! including newly formed precipitations (i.e., condensates formed by the
	944	! condensation process), while zoliq is the total amount of condensates in
	945	! the cloud (i.e., on which precipitation processes have applied)
	946	DO i = 1, klon
	947	zoliq(i) = zcond(i)
	948	zoliql(i) = zcond(i) * ( 1. - zfice(i) )
	949	zoliqi(i) = zcond(i) * zfice(i)
	950	ENDDO
	951
[5579]	952	IF (ok_plane_contrail) THEN
[5602]	953	!!--Contrails do not precipitate. We remove then from the variables temporarily
	954	!DO i = 1, klon
	955	! rneb(i,k) = rneb(i,k) - contfra(i)
	956	! zoliqi(i) = zoliqi(i) - ( qcont(i) - zqs(i) * contfra(i) )
	957	!ENDDO
	958	!--Contrails precipitate as natural clouds. We save the partition of ice
	959	!--between natural clouds and contrails
	960	!--NB. we use qcont as a temporary variable to save this partition
[5593]	961	DO i = 1, klon
[5602]	962	IF ( zoliqi(i) .GT. 0. ) THEN
	963	qcont(i) = ( qcont(i) - zqs(i) * contfra(i) ) / zoliqi(i)
	964	ELSE
	965	qcont(i) = 0.
	966	ENDIF
[5593]	967	ENDDO
[5579]	968	ENDIF
	969
[4803]	970	!================================================================
[5408]	971	! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
[4803]	972	IF (ok_poprecip) THEN
	973
	974	CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	975	ctot_vol(:,k), ptconv(:,k), &
	976	zt, zq, zoliql, zoliqi, zfice, &
[5589]	977	rneb(:,k), zqice_sedim, znebprecipclr, znebprecipcld, &
[4803]	978	zrfl, zrflclr, zrflcld, &
	979	zifl, ziflclr, ziflcld, &
[4830]	980	qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
[4819]	981	dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
	982	dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
[5589]	983	dqsmelt(:,k), dqsfreez(:,k), dqised(:,k) &
[4803]	984	)
[5431]	985	DO i = 1, klon
	986	zoliq(i) = zoliql(i) + zoliqi(i)
	987	ENDDO
[4803]	988
	989	!================================================================
	990	ELSE
	991
[5413]	992	CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	993	ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
	994	zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
	995	rneb(:,k), znebprecipclr, znebprecipcld, &
	996	zneb, tot_zneb, zrho_up, zvelo_up, &
	997	zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
	998	zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
	999	)
[3999]	1000
[4803]	1001	ENDIF ! ok_poprecip
[5579]	1002
	1003	IF (ok_plane_contrail) THEN
[5602]	1004	!!--Contrails are reintroduced in the variables
	1005	!DO i = 1, klon
	1006	! rneb(i,k) = rneb(i,k) + contfra(i)
	1007	! zoliqi(i) = zoliqi(i) + ( qcont(i) - zqs(i) * contfra(i) )
	1008	!ENDDO
	1009	!--Contrails fraction is left unchanged, but contrails water has changed
[5593]	1010	DO i = 1, klon
[5602]	1011	IF ( zoliqi(i) .LE. 0. ) THEN
	1012	contfra(i) = 0.
	1013	qcont(i) = 0.
	1014	ELSE
	1015	qcont(i) = zqs(i) * contfra(i) + zoliqi(i) * qcont(i)
	1016	ENDIF
[5593]	1017	ENDDO
[5579]	1018	ENDIF
[4803]	1019
[5408]	1020	! End of precipitation processes after cloud formation
	1021	! ----------------------------------------------------
[3999]	1022
[5406]	1023	!----------------------------------------------------------------------
[5408]	1024	! P4> Calculation of cloud condensates amount seen by radiative scheme
[5406]	1025	!----------------------------------------------------------------------
[4803]	1026
[5413]	1027	DO i=1,klon
[3999]	1028
[5413]	1029	IF (ok_poprecip) THEN
	1030	IF (ok_radocond_snow) THEN
	1031	radocond(i,k) = zoliq(i)
	1032	zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
	1033	ELSE
	1034	radocond(i,k) = zoliq(i)
	1035	zradoice(i) = zoliqi(i)
	1036	ENDIF
	1037	ELSE
	1038	radocond(i,k) = zradocond(i)
	1039	ENDIF
[4114]	1040
[5413]	1041	! calculate the percentage of ice in "radocond" so cloud+precip seen
	1042	! by the radiation scheme
	1043	IF (radocond(i,k) .GT. 0.) THEN
	1044	radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
	1045	ENDIF
[4114]	1046	ENDDO
	1047
[3999]	1048	! ----------------------------------------------------------------
[5408]	1049	! P5> Wet scavenging
[3999]	1050	! ----------------------------------------------------------------
	1051
	1052	!Scavenging through nucleation in the layer
	1053
	1054	DO i = 1,klon
	1055
	1056	IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
	1057	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
	1058	ELSE
	1059	beta(i,k) = 0.
	1060	ENDIF
	1061
[4059]	1062	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
[3999]	1063
	1064	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1065
[4654]	1066	IF (temp(i,k) .GE. t_glace_min) THEN
[3999]	1067	zalpha_tr = a_tr_sca(3)
	1068	ELSE
	1069	zalpha_tr = a_tr_sca(4)
	1070	ENDIF
	1071
[5450]	1072	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
	1073	frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
[4226]	1074
[3999]	1075	! Nucleation with a factor of -1 instead of -0.5
[5450]	1076	zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
[3999]	1077
	1078	ENDIF
	1079
[4226]	1080	ENDDO
	1081
[3999]	1082	! Scavenging through impaction in the underlying layer
	1083
	1084	DO kk = k-1, 1, -1
	1085
	1086	DO i = 1, klon
	1087
	1088	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1089
[4654]	1090	IF (temp(i,kk) .GE. t_glace_min) THEN
[3999]	1091	zalpha_tr = a_tr_sca(1)
	1092	ELSE
	1093	zalpha_tr = a_tr_sca(2)
	1094	ENDIF
	1095
[5450]	1096	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
	1097	frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
[3999]	1098
	1099	ENDIF
	1100
	1101	ENDDO
	1102
	1103	ENDDO
[5406]	1104
	1105	!------------------------------------------------------------
[5408]	1106	! P6 > write diagnostics and outputs
[5406]	1107	!------------------------------------------------------------
	1108
	1109	!--AB Write diagnostics and tracers for ice supersaturation
[5601]	1110	IF ( ok_plane_contrail ) THEN
	1111	cfa_seri(:,k) = contfra(:)
	1112	qta_seri(:,k) = qcont(:)
	1113	qice_cont(:,k) = qcont(:) - zqs(:) * contfra(:)
	1114	ENDIF
	1115
	1116	IF ( ok_ice_supersat .AND. .NOT. ok_unadjusted_clouds) qvc(:) = zqs(:) * rneb(:,k)
	1117
[5406]	1118	IF ( ok_ice_supersat ) THEN
	1119	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
	1120	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
[4226]	1121
[5406]	1122	DO i = 1, klon
	1123
[5545]	1124	cf_seri(i,k) = rneb(i,k)
	1125
[5406]	1126	IF ( zoliq(i) .LE. 0. ) THEN
	1127	!--If everything was precipitated, the remaining empty cloud is dissipated
	1128	!--and everything is transfered to the subsaturated clear sky region
[5545]	1129	!--NB. we do not change rneb, as it is a diagnostic only
	1130	cf_seri(i,k) = 0.
[5406]	1131	qvc(i) = 0.
	1132	ENDIF
	1133
[5601]	1134	qvc_seri(i,k) = qvc(i)
[5406]	1135
	1136	!--Diagnostics
	1137	gamma_cond(i,k) = gammasat(i)
	1138	subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
	1139	qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
[5412]	1140	qcld(i,k) = qvc(i) + zoliq(i)
[5575]	1141
	1142	!--Calculation of the ice supersaturated fraction following Lamquin et al (2012)
	1143	!--methodology: in each layer, we make a maximum random overlap assumption for
	1144	!--ice supersaturation
	1145	IF ( ( paprs(i,k) .GT. 10000. ) .AND. ( paprs(i,k) .LE. 15000. ) ) THEN
	1146	IF ( issrfra100to150UP(i) .GT. ( 1. - eps ) ) THEN
	1147	issrfra100to150(i) = 1.
	1148	ELSE
	1149	issrfra100to150(i) = 1. - ( 1. - issrfra100to150(i) ) * &
	1150	( 1. - MAX( issrfra(i,k), issrfra100to150UP(i) ) ) &
	1151	/ ( 1. - issrfra100to150UP(i) )
	1152	issrfra100to150UP(i) = issrfra(i,k)
	1153	ENDIF
	1154	ELSEIF ( ( paprs(i,k) .GT. 15000. ) .AND. ( paprs(i,k) .LE. 20000. ) ) THEN
	1155	IF ( issrfra150to200UP(i) .GT. ( 1. - eps ) ) THEN
	1156	issrfra150to200(i) = 1.
	1157	ELSE
	1158	issrfra150to200(i) = 1. - ( 1. - issrfra150to200(i) ) * &
	1159	( 1. - MAX( issrfra(i,k), issrfra150to200UP(i) ) ) &
	1160	/ ( 1. - issrfra150to200UP(i) )
	1161	issrfra150to200UP(i) = issrfra(i,k)
	1162	ENDIF
	1163	ELSEIF ( ( paprs(i,k) .GT. 20000. ) .AND. ( paprs(i,k) .LE. 25000. ) ) THEN
	1164	IF ( issrfra200to250UP(i) .GT. ( 1. - eps ) ) THEN
	1165	issrfra200to250(i) = 1.
	1166	ELSE
	1167	issrfra200to250(i) = 1. - ( 1. - issrfra200to250(i) ) * &
	1168	( 1. - MAX( issrfra(i,k), issrfra200to250UP(i) ) ) &
	1169	/ ( 1. - issrfra200to250UP(i) )
	1170	issrfra200to250UP(i) = issrfra(i,k)
	1171	ENDIF
	1172	ELSEIF ( ( paprs(i,k) .GT. 25000. ) .AND. ( paprs(i,k) .LE. 30000. ) ) THEN
	1173	IF ( issrfra250to300UP(i) .GT. ( 1. - eps ) ) THEN
	1174	issrfra250to300(i) = 1.
	1175	ELSE
	1176	issrfra250to300(i) = 1. - ( 1. - issrfra250to300(i) ) * &
	1177	( 1. - MAX( issrfra(i,k), issrfra250to300UP(i) ) ) &
	1178	/ ( 1. - issrfra250to300UP(i) )
	1179	issrfra250to300UP(i) = issrfra(i,k)
	1180	ENDIF
	1181	ELSEIF ( ( paprs(i,k) .GT. 30000. ) .AND. ( paprs(i,k) .LE. 40000. ) ) THEN
	1182	IF ( issrfra300to400UP(i) .GT. ( 1. - eps ) ) THEN
	1183	issrfra300to400(i) = 1.
	1184	ELSE
	1185	issrfra300to400(i) = 1. - ( 1. - issrfra300to400(i) ) * &
	1186	( 1. - MAX( issrfra(i,k), issrfra300to400UP(i) ) ) &
	1187	/ ( 1. - issrfra300to400UP(i) )
	1188	issrfra300to400UP(i) = issrfra(i,k)
	1189	ENDIF
	1190	ELSEIF ( ( paprs(i,k) .GT. 40000. ) .AND. ( paprs(i,k) .LE. 50000. ) ) THEN
	1191	IF ( issrfra400to500UP(i) .GT. ( 1. - eps ) ) THEN
	1192	issrfra400to500(i) = 1.
	1193	ELSE
	1194	issrfra400to500(i) = 1. - ( 1. - issrfra400to500(i) ) * &
	1195	( 1. - MAX( issrfra(i,k), issrfra400to500UP(i) ) ) &
	1196	/ ( 1. - issrfra400to500UP(i) )
	1197	issrfra400to500UP(i) = issrfra(i,k)
	1198	ENDIF
	1199	ENDIF
	1200
[5406]	1201	ENDDO
	1202	ENDIF
	1203
[4803]	1204	! Outputs:
	1205	!-------------------------------
	1206	! Precipitation fluxes at layer interfaces
	1207	! + precipitation fractions +
	1208	! temperature and water species tendencies
	1209	DO i = 1, klon
	1210	psfl(i,k)=zifl(i)
	1211	prfl(i,k)=zrfl(i)
	1212	pfraclr(i,k)=znebprecipclr(i)
	1213	pfracld(i,k)=znebprecipcld(i)
	1214	d_q(i,k) = zq(i) - qt(i,k)
	1215	d_t(i,k) = zt(i) - temp(i,k)
[5412]	1216
	1217	IF (ok_bug_phase_lscp) THEN
	1218	d_ql(i,k) = (1-zfice(i))*zoliq(i)
	1219	d_qi(i,k) = zfice(i)*zoliq(i)
	1220	ELSE
	1221	d_ql(i,k) = zoliql(i)
	1222	d_qi(i,k) = zoliqi(i)
	1223	ENDIF
	1224
[4803]	1225	ENDDO
	1226
	1227
[5413]	1228	ENDDO ! loop on k from top to bottom
[4059]	1229
[3999]	1230
	1231	! Rain or snow at the surface (depending on the first layer temperature)
	1232	DO i = 1, klon
	1233	snow(i) = zifl(i)
	1234	rain(i) = zrfl(i)
[4392]	1235	! c_iso final output
[3999]	1236	ENDDO
	1237
	1238	IF (ncoreczq>0) THEN
	1239	WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
	1240	ENDIF
	1241
[4654]	1242
	1243	RETURN
	1244
[4380]	1245	END SUBROUTINE lscp
[3999]	1246	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1247
[4664]	1248	END MODULE lmdz_lscp

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

Download in other formats:

Original Format