Context Navigation

lmdz_lscp.f90 @ 5631

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

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

Download in other formats:

Original Format