Context Navigation

lmdz_lscp.f90 @ 5445

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

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

Download in other formats:

Original Format