Context Navigation

lmdz_lscp.f90 @ 5601

Last change on this file since 5601 was 5601, checked in by aborella, 3 months ago

Multiple changes:

tracers which were ratios are now absolute quantities. This is needed because when the ratio

is not defined, some aberrations may occur

added a new tracer for total specific humidity in contrails
rework of the mixing process for cirrus clouds (and contrails)
changed the numerical integration of ice crystals' sublimation
subroutines do not take real inputs anymore (at least klon tables)
added more radiative diagnostics for contrails

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

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

Download in other formats:

Original Format