Context Navigation

lmdz_lscp.F90 @ 5434

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

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

Download in other formats:

Original Format