Context Navigation

lmdz_lscp.f90 @ 5405

Last change on this file since 5405 was 5396, checked in by evignon, 2 months ago
ajout de ql_seri et qi_seri dans lmdz_lscp
File size: 71.7 KB

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

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

Download in other formats:

Original Format