Context Navigation

lmdz_lscp_precip.f90 @ 5440

Last change on this file since 5440 was 5431, checked in by aborella, 2 days ago
Merge with trunk
File size: 79.1 KB

Rev	Line
[5408]	1	MODULE lmdz_lscp_precip
	2	!----------------------------------------------------------------
	3	! Module for the process-oriented treament of precipitation
	4	! that are called in LSCP
	5	! Authors: Atelier Nuage (G. Riviere, L. Raillard, M. Wimmer,
	6	! N. Dutrievoz, E. Vignon, A. Borella, et al.)
	7	! Jan. 2024
	8
	9
	10	IMPLICIT NONE
	11
	12	CONTAINS
	13
	14	!----------------------------------------------------------------
	15	! historical (till CMIP6) version of the pre-cloud formation
	16	! precipitation scheme containing precip evaporation and melting
	17
	18	SUBROUTINE histprecip_precld( &
	19	klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
[5411]	20	zmqc, zneb, znebprecip, znebprecipclr, &
[5408]	21	zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
	22	)
	23
	24	USE lmdz_lscp_ini, ONLY : iflag_evap_prec
	25	USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, ztfondue
	26	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
	27	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
	28
	29	IMPLICIT NONE
	30
	31
	32	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
	33	REAL, INTENT(IN) :: dtime !--time step [s]
	34	LOGICAL, INTENT(IN) :: iftop !--if top of the column
	35
	36
	37	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
	38	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
	39	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
	40
	41	REAL, INTENT(INOUT), DIMENSION(klon) :: zt !--current temperature [K]
	42	REAL, INTENT(IN), DIMENSION(klon) :: ztupnew !--updated temperature of the overlying layer [K]
	43
	44	REAL, INTENT(INOUT), DIMENSION(klon) :: zq !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
	45	REAL, INTENT(INOUT), DIMENSION(klon) :: zmqc !--specific humidity in the precipitation falling from the upper layer [kg/kg]
	46
	47	REAL, INTENT(IN), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-]
[5411]	48	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecip !--fraction of precipitation IN THE LAYER ABOVE [-]
[5408]	49	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
	50
	51	REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
	52	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
	53	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
	54	REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
	55	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
	56	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
	57
	58	REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s]
	59	REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
	60
	61
	62	REAL :: zmair
	63	REAL :: zcpair, zcpeau
	64	REAL, DIMENSION(klon) :: qzero
	65	REAL, DIMENSION(klon) :: zqs, zdqs
	66	REAL, DIMENSION(klon) :: qsl, qsi
	67	REAL, DIMENSION(klon) :: dqsl, dqsi
	68
	69	REAL :: zqev0
	70	REAL :: zqev, zqevt
	71	REAL :: zqevi, zqevti
	72
	73	REAL, DIMENSION(klon) :: zrfln, zifln
	74
	75	REAL :: zmelt
	76	INTEGER :: i
	77
	78	qzero(:) = 0.
	79
	80
	81	! --------------------------------------------------------------------
	82	! P1> Thermalization of precipitation falling from the overlying layer
	83	! --------------------------------------------------------------------
	84	! Computes air temperature variation due to enthalpy transported by
	85	! precipitation. Precipitation is then thermalized with the air in the
	86	! layer.
	87	! The precipitation should remain thermalized throughout the different
	88	! thermodynamical transformations.
	89	! The corresponding water mass should
	90	! be added when calculating the layer's enthalpy change with
	91	! temperature
	92	! ---------------------------------------------------------------------
	93
	94	IF (iftop) THEN
	95
	96	DO i = 1, klon
	97	zmqc(i) = 0.
	98	ENDDO
	99
	100	ELSE
	101
	102	DO i = 1, klon
	103
	104	zmair=(paprsdn(i)-paprsup(i))/RG
	105	! no condensed water so cp=cp(vapor+dry air)
	106	! RVTMP2=rcpv/rcpd-1
	107	zcpair=RCPD(1.0+RVTMP2zq(i))
	108	zcpeau=RCPD*RVTMP2
	109
	110	! zmqc: precipitation mass that has to be thermalized with
	111	! layer's air so that precipitation at the ground has the
	112	! same temperature as the lowermost layer
	113	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
	114	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
	115	zt(i) = ( ztupnew(i)zmqc(i)zcpeau + zcpair*zt(i) ) &
	116	/ (zcpair + zmqc(i)*zcpeau)
	117
	118	ENDDO
	119
	120	ENDIF
	121
	122	! --------------------------------------------------------------------
	123	! P2> Precipitation evaporation/sublimation/melting
	124	! --------------------------------------------------------------------
	125	! A part of the precipitation coming from above is evaporated/sublimated/melted.
	126	! --------------------------------------------------------------------
	127
	128	IF (iflag_evap_prec.GE.1) THEN
	129
	130	! Calculation of saturation specific humidity
	131	! depending on temperature:
	132	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,0,.false.,zqs,zdqs)
	133	! wrt liquid water
	134	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,1,.false.,qsl,dqsl)
	135	! wrt ice
	136	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,2,.false.,qsi,dqsi)
	137
	138	DO i = 1, klon
	139
	140	! if precipitation
	141	IF (zrfl(i)+zifl(i).GT.0.) THEN
	142
	143	! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
	144	! c_iso: likely important to distinguish cs from neb isotope precipitation
	145
	146	IF (iflag_evap_prec.GE.4) THEN
	147	zrfl(i) = zrflclr(i)
	148	zifl(i) = ziflclr(i)
	149	ENDIF
	150
	151	IF (iflag_evap_prec.EQ.1) THEN
	152	znebprecip(i)=zneb(i)
	153	ELSE
	154	znebprecip(i)=MAX(zneb(i),znebprecip(i))
	155	ENDIF
	156
	157	IF (iflag_evap_prec.GT.4) THEN
	158	! Max evaporation not to saturate the clear sky precip fraction
	159	! i.e. the fraction where evaporation occurs
	160	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
	161	ELSEIF (iflag_evap_prec .EQ. 4) THEN
	162	! Max evaporation not to saturate the whole mesh
	163	! Pay attention -> lead to unrealistic and excessive evaporation
	164	zqev0 = MAX(0.0, zqs(i)-zq(i))
	165	ELSE
	166	! Max evap not to saturate the fraction below the cloud
	167	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
	168	ENDIF
	169
	170	! Evaporation of liquid precipitation coming from above
	171	! dP/dz=beta(1-q/qsat)sqrt(P)
	172	! formula from Sundquist 1988, Klemp & Wilhemson 1978
	173	! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
	174
	175	IF (iflag_evap_prec.EQ.3) THEN
	176	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl(i)) &
	177	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
	178	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
	179	ELSE IF (iflag_evap_prec.GE.4) THEN
	180	zqevt = znebprecipclr(i)coef_eva(1.0-zq(i)/qsl(i)) &
	181	*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
	182	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
	183	ELSE
	184	zqevt = 1.coef_eva(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
	185	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
	186	ENDIF
	187
	188	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * RG*dtime/(paprsdn(i)-paprsup(i))
	189
	190	! sublimation of the solid precipitation coming from above
	191	IF (iflag_evap_prec.EQ.3) THEN
	192	zqevti = znebprecip(i)coef_sub(1.0-zq(i)/qsi(i)) &
	193	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
	194	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
	195	ELSE IF (iflag_evap_prec.GE.4) THEN
	196	zqevti = znebprecipclr(i)coef_sub(1.0-zq(i)/qsi(i)) &
	197	*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
	198	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
	199	ELSE
	200	zqevti = 1.coef_sub(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
	201	(paprsdn(i)-paprsup(i))/pplay(i)zt(i)*RD/RG
	202	ENDIF
	203
	204	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) * RG*dtime/(paprsdn(i)-paprsup(i))
	205
	206	! A. JAM
	207	! Evaporation limit: we ensure that the layer's fraction below
	208	! the cloud or the whole mesh (depending on iflag_evap_prec)
	209	! does not reach saturation. In this case, we
	210	! redistribute zqev0 conserving the ratio liquid/ice
	211
	212	IF (zqevt+zqevti.GT.zqev0) THEN
	213	zqev=zqev0*zqevt/(zqevt+zqevti)
	214	zqevi=zqev0*zqevti/(zqevt+zqevti)
	215	ELSE
	216	zqev=zqevt
	217	zqevi=zqevti
	218	ENDIF
	219
	220	!--Diagnostics
	221	dqreva(i) = - zqev / dtime
	222	dqssub(i) = - zqevti / dtime
	223
	224	! New solid and liquid precipitation fluxes after evap and sublimation
	225	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprsdn(i)-paprsup(i))/RG/dtime)
	226	zifln(i) = Max(0.,zifl(i) - zqevi*(paprsdn(i)-paprsup(i))/RG/dtime)
	227
	228
	229	! vapor, temperature, precip fluxes update
	230	! vapor is updated after evaporation/sublimation (it is increased)
	231	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	232	* (RG/(paprsdn(i)-paprsup(i)))*dtime
	233	! zmqc is the total condensed water in the precip flux (it is decreased)
	234	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	235	* (RG/(paprsdn(i)-paprsup(i)))*dtime
	236	! air and precip temperature (i.e., gridbox temperature)
	237	! is updated due to latent heat cooling
	238	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	239	* (RG/(paprsdn(i)-paprsup(i)))*dtime &
	240	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
	241	+ (zifln(i)-zifl(i)) &
	242	* (RG/(paprsdn(i)-paprsup(i)))*dtime &
	243	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
	244
	245	! New values of liquid and solid precipitation
	246	zrfl(i) = zrfln(i)
	247	zifl(i) = zifln(i)
	248
	249	! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
	250	! due to evap + sublim
	251
	252
	253	IF (iflag_evap_prec.GE.4) THEN
	254	zrflclr(i) = zrfl(i)
	255	ziflclr(i) = zifl(i)
	256	IF(zrflclr(i) + ziflclr(i).LE.0) THEN
	257	znebprecipclr(i) = 0.0
	258	ENDIF
	259	zrfl(i) = zrflclr(i) + zrflcld(i)
	260	zifl(i) = ziflclr(i) + ziflcld(i)
	261	ENDIF
	262
	263	! c_iso duplicate for isotopes or loop on isotopes
	264
	265	! Melting:
	266	zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
	267	! precip fraction that is melted
	268	zmelt = MIN(MAX(zmelt,0.),1.)
	269
	270	! update of rainfall and snowfall due to melting
	271	IF (iflag_evap_prec.GE.4) THEN
	272	zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
	273	zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
	274	zrfl(i)=zrflclr(i)+zrflcld(i)
	275	ELSE
	276	zrfl(i)=zrfl(i)+zmelt*zifl(i)
	277	ENDIF
	278
	279
	280	! c_iso: melting of isotopic precipi with zmelt (no fractionation)
	281
	282	! Latent heat of melting because of precipitation melting
	283	! NB: the air + precip temperature is simultaneously updated
	284	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprsdn(i)-paprsup(i)) &
	285	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
	286
	287	IF (iflag_evap_prec.GE.4) THEN
	288	ziflclr(i)=ziflclr(i)*(1.-zmelt)
	289	ziflcld(i)=ziflcld(i)*(1.-zmelt)
	290	zifl(i)=ziflclr(i)+ziflcld(i)
	291	ELSE
	292	zifl(i)=zifl(i)*(1.-zmelt)
	293	ENDIF
	294
	295	ELSE
	296	! if no precip, we reinitialize the cloud fraction used for the precip to 0
	297	znebprecip(i)=0.
	298
	299	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	300
	301	ENDDO ! loop on klon
	302
	303	ENDIF ! (iflag_evap_prec>=1)
	304
	305	END SUBROUTINE histprecip_precld
	306
	307
[5413]	308	SUBROUTINE histprecip_postcld( &
	309	klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, &
	310	zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
	311	rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
	312	zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
	313	zradocond, zradoice, dqrauto, dqsauto &
	314	)
	315
	316	USE lmdz_lscp_ini, ONLY : RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT
	317	USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, ffallv_con, &
	318	cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, ffallv_lsc, &
	319	seuil_neb, rain_int_min, cice_velo, dice_velo
	320	USE lmdz_lscp_ini, ONLY : iflag_evap_prec, iflag_cloudth_vert, iflag_rain_incloud_vol, &
	321	iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
	322	niter_lscp
	323	USE lmdz_lscp_tools, ONLY : fallice_velocity
	324
	325	IMPLICIT NONE
	326
	327
	328	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
	329	REAL, INTENT(IN) :: dtime !--time step [s]
	330	LOGICAL, INTENT(IN) :: iftop !--if top of the column
	331
	332	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
	333	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
	334	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
	335	REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-]
	336	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point
	337	REAL, INTENT(IN), DIMENSION(klon) :: zdqsdT_raw !--derivative of qsat w.r.t. temperature times L/Cp [SI]
	338
	339	REAL, INTENT(INOUT), DIMENSION(klon) :: zt !--current temperature [K]
	340	REAL, INTENT(INOUT), DIMENSION(klon) :: zq !--current water vapor specific humidity [kg/kg]
	341	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliq !--current liquid+ice water specific humidity [kg/kg]
	342	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliql !--current liquid water specific humidity [kg/kg]
	343	REAL, INTENT(INOUT), DIMENSION(klon) :: zoliqi !--current ice water specific humidity [kg/kg]
	344	REAL, INTENT(INOUT), DIMENSION(klon) :: zcond !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg]
	345	REAL, INTENT(IN), DIMENSION(klon) :: zfice !--ice fraction AFTER CONDENSATION [-]
	346	REAL, INTENT(IN), DIMENSION(klon) :: zmqc !--specific humidity in the precipitation falling from the upper layer [kg/kg]
	347
	348	REAL, INTENT(IN), DIMENSION(klon) :: rneb !--cloud fraction [-]
	349	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
	350	REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipcld !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-]
	351	REAL, INTENT(INOUT), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-]
	352	REAL, INTENT(INOUT), DIMENSION(klon) :: tot_zneb !--total cloud cover above the current mesh [-]
	353	REAL, INTENT(INOUT), DIMENSION(klon) :: zrho_up !--air density IN THE LAYER ABOVE [kg/m3]
	354	REAL, INTENT(INOUT), DIMENSION(klon) :: zvelo_up !--ice fallspeed velocity IN THE LAYER ABOVE [m/s]
	355
	356	REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean [kg/s/m2]
	357	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky [kg/s/m2]
	358	REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air [kg/s/m2]
	359	REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean [kg/s/m2]
	360	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky [kg/s/m2]
	361	REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
	362
	363	REAL, INTENT(OUT), DIMENSION(klon) :: zradocond !--condensed water used in the radiation scheme [kg/kg]
	364	REAL, INTENT(OUT), DIMENSION(klon) :: zradoice !--condensed ice water used in the radiation scheme [kg/kg]
	365	REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
	366	REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
	367
	368
	369	! Local variables for precip fraction update
	370	REAL :: smallestreal
	371	REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb
	372	REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld
	373	REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
	374	REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
	375
	376	! Local variables for autoconversion
	377	REAL :: zct, zcl, zexpo, ffallv
	378	REAL :: zchau, zfroi
	379	REAL :: zrain, zsnow, zprecip
	380	REAL :: effective_zneb
	381	REAL, DIMENSION(klon) :: zrho, zvelo
	382	REAL, DIMENSION(klon) :: zdz, iwc
	383
	384	! Local variables for Bergeron process
	385	REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
	386	REAL, DIMENSION(klon) :: zqpreci, zqprecl
	387
	388	! Local variables for calculation of quantity of condensates seen by the radiative scheme
	389	REAL, DIMENSION(klon) :: zradoliq
	390	REAL, DIMENSION(klon) :: ziflprev
	391
	392	! Misc
	393	INTEGER :: i, n
	394
	395	! Initialisation
	396	smallestreal=1.e-9
	397	zqprecl(:) = 0.
	398	zqpreci(:) = 0.
	399	ziflprev(:) = 0.
	400
	401
	402	IF (iflag_evap_prec .GE. 4) THEN
	403
	404	!Partitionning between precipitation coming from clouds and that coming from CS
	405
	406	!0) Calculate tot_zneb, total cloud fraction above the cloud
	407	!assuming a maximum-random overlap (voir Jakob and Klein, 2000)
	408
	409	DO i=1, klon
	410	tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i),zneb(i))) &
	411	/(1.-min(zneb(i),1.-smallestreal))
	412	d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
	413	tot_zneb(i) = tot_znebn(i)
	414
	415
	416	!1) Cloudy to clear air
	417	d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i),znebprecipcld(i))
	418	IF (znebprecipcld(i) .GT. 0.) THEN
	419	d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
	420	d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
	421	ELSE
	422	d_zrfl_cld_clr(i) = 0.
	423	d_zifl_cld_clr(i) = 0.
	424	ENDIF
	425
	426	!2) Clear to cloudy air
	427	d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i) - d_tot_zneb(i) - zneb(i)))
	428	IF (znebprecipclr(i) .GT. 0) THEN
	429	d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
	430	d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
	431	ELSE
	432	d_zrfl_clr_cld(i) = 0.
	433	d_zifl_clr_cld(i) = 0.
	434	ENDIF
	435
	436	!Update variables
	437	znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
	438	znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
	439	zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
	440	ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
	441	zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
	442	ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
	443
	444	! c_iso do the same thing for isotopes precip
	445	ENDDO
	446	ENDIF
	447
	448
	449	! Autoconversion
	450	! -------------------------------------------------------------------------------
	451
	452	! Initialisation
	453	DO i = 1, klon
	454	zrho(i) = pplay(i) / zt(i) / RD
	455	zdz(i) = (paprsdn(i)-paprsup(i)) / (zrho(i)*RG)
	456	iwc(i) = 0.
	457	zneb(i) = MAX(rneb(i),seuil_neb)
	458
	459	! variables for calculation of quantity of condensates seen by the radiative scheme
	460	! NB. only zradocond and zradoice are outputed, but zradoliq is used if ok_radocond_snow
	461	! is TRUE
	462	zradocond(i) = zoliq(i)/REAL(niter_lscp+1)
	463	zradoliq(i) = zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
	464	zradoice(i) = zoliq(i)*zfice(i)/REAL(niter_lscp+1)
	465	ENDDO
	466
	467
	468	DO n = 1, niter_lscp
	469
	470	DO i=1,klon
	471	IF (rneb(i).GT.0.0) THEN
	472	iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
	473	ENDIF
	474	ENDDO
	475
	476	CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo)
	477
	478	DO i = 1, klon
	479
	480	IF (rneb(i).GT.0.0) THEN
	481
	482	! Initialization of zrain, zsnow and zprecip:
	483	zrain=0.
	484	zsnow=0.
	485	zprecip=0.
	486	! c_iso same init for isotopes. Externalisation?
	487
	488	IF (zneb(i).EQ.seuil_neb) THEN
	489	zprecip = 0.0
	490	zsnow = 0.0
	491	zrain= 0.0
	492	ELSE
	493
	494	IF (ptconv(i)) THEN ! if convective point
	495	zcl=cld_lc_con
	496	zct=1./cld_tau_con
	497	zexpo=cld_expo_con
	498	ffallv=ffallv_con
	499	ELSE
	500	zcl=cld_lc_lsc
	501	zct=1./cld_tau_lsc
	502	zexpo=cld_expo_lsc
	503	ffallv=ffallv_lsc
	504	ENDIF
	505
	506
	507	! if vertical heterogeneity is taken into account, we use
	508	! the "true" volume fraction instead of a modified
	509	! surface fraction (which is larger and artificially
	510	! reduces the in-cloud water).
	511	IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
	512	effective_zneb=ctot_vol(i)
	513	ELSE
	514	effective_zneb=zneb(i)
	515	ENDIF
	516
	517
	518	! Liquid water quantity to remove: zchau (Sundqvist, 1978)
	519	! dqliq/dt=-qliq/tau(1-exp(-qcin/clw)*2)
	520	!.........................................................
	521	IF (iflag_autoconversion .EQ. 2) THEN
	522	! two-steps resolution with niter_lscp=1 sufficient
	523	! we first treat the second term (with exponential) in an explicit way
	524	! and then treat the first term (-q/tau) in an exact way
	525	zchau=zoliql(i)(1.-exp(-dtime/REAL(niter_lscp)zct &
	526	(1.-exp(-(zoliql(i)/effective_zneb/zcl)*zexpo))))
	527	ELSE
	528	! old explicit resolution with subtimesteps
	529	zchau = zctdtime/REAL(niter_lscp)zoliql(i) &
	530	(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)*zexpo))
	531	ENDIF
	532
	533
	534	! Ice water quantity to remove (Zender & Kiehl, 1997)
	535	! dqice/dt=1/rhod(rhowice*qice)/dz
	536	!....................................
	537	IF (iflag_autoconversion .EQ. 2) THEN
	538	! exact resolution, niter_lscp=1 is sufficient but works only
	539	! with iflag_vice=0
	540	IF (zoliqi(i) .GT. 0.) THEN
	541	zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
	542	+dice_velodtime/REAL(niter_lscp)cice_velo/zdz(i)ffallv)*(-1./dice_velo))
	543	ELSE
	544	zfroi=0.
	545	ENDIF
	546	ELSE
	547	! old explicit resolution with subtimesteps
	548	zfroi = dtime/REAL(niter_lscp)/zdz(i)zoliqi(i)zvelo(i)
	549	ENDIF
	550
	551	zrain = MIN(MAX(zchau,0.0),zoliql(i))
	552	zsnow = MIN(MAX(zfroi,0.0),zoliqi(i))
	553	zprecip = MAX(zrain + zsnow,0.0)
	554
	555	ENDIF
	556
	557
	558	IF (iflag_autoconversion .GE. 1) THEN
	559	! debugged version with phase conservation through the autoconversion process
	560	zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0)
	561	zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0)
	562	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
	563	ELSE
	564	! bugged version with phase resetting
	565	zoliql(i) = MAX(zoliq(i)(1.-zfice(i))-1.zrain , 0.0)
	566	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zsnow , 0.0)
	567	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
	568	ENDIF
	569
	570	! c_iso: call isotope_conversion (for readibility) or duplicate
	571
	572	! variables for calculation of quantity of condensates seen by the radiative scheme
	573	zradocond(i) = zradocond(i) + zoliq(i)/REAL(niter_lscp+1)
	574	zradoliq(i) = zradoliq(i) + zoliql(i)/REAL(niter_lscp+1)
	575	zradoice(i) = zradoice(i) + zoliqi(i)/REAL(niter_lscp+1)
	576
	577	!--Diagnostics
	578	dqrauto(i) = dqrauto(i) + zrain / dtime
	579	dqsauto(i) = dqsauto(i) + zsnow / dtime
	580
	581	ENDIF ! rneb >0
	582
	583	ENDDO ! i = 1,klon
	584
	585	ENDDO ! n = 1,niter
	586
	587	! Precipitation flux calculation
	588
	589	DO i = 1, klon
	590
	591	IF (iflag_evap_prec.GE.4) THEN
	592	ziflprev(i)=ziflcld(i)
	593	ELSE
	594	ziflprev(i)=zifl(i)*zneb(i)
	595	ENDIF
	596
	597	IF (rneb(i) .GT. 0.0) THEN
	598
	599	! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
	600	! If T<0C, liquid precip are converted into ice, which leads to an increase in
	601	! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
	602	! taken into account through a linearization of the equations and by approximating
	603	! the liquid precipitation process with a "threshold" process. We assume that
	604	! condensates are not modified during this operation. Liquid precipitation is
	605	! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
	606	! Water vapor increases as well
	607	! Note that compared to fisrtilp, we always assume iflag_bergeron=2
	608
	609	IF (ok_bug_phase_lscp) THEN
	610	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
	611	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
	612	ELSE
	613	zqpreci(i)=zcond(i)*zfice(i)-zoliqi(i)
	614	zqprecl(i)=zcond(i)*(1.-zfice(i))-zoliql(i)
	615	ENDIF
	616	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
	617	coef1 = rneb(i)RLSTT/zcpzdqsdT_raw(i)
	618	! Computation of DT if all the liquid precip freezes
	619	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
	620	! T should not exceed the freezing point
	621	! that is Delta > RTT-zt(i)
	622	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
	623	zt(i) = zt(i) + DeltaT
	624	! water vaporization due to temp. increase
	625	Deltaq = rneb(i)zdqsdT_raw(i)DeltaT
	626	! we add this vaporized water to the total vapor and we remove it from the precip
	627	zq(i) = zq(i) + Deltaq
	628	! The three "max" lines herebelow protect from rounding errors
	629	zcond(i) = max( zcond(i)- Deltaq, 0. )
	630	! liquid precipitation converted to ice precip
	631	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
	632	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
	633	! iced water budget
	634	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
	635
	636	! c_iso : mv here condensation of isotopes + redispatchage en precip
	637
	638	IF (iflag_evap_prec.GE.4) THEN
	639	zrflcld(i) = zrflcld(i)+zqprecl(i) &
	640	(paprsdn(i)-paprsup(i))/(RGdtime)
	641	ziflcld(i) = ziflcld(i)+ zqpreci(i) &
	642	(paprsdn(i)-paprsup(i))/(RGdtime)
	643	znebprecipcld(i) = rneb(i)
	644	zrfl(i) = zrflcld(i) + zrflclr(i)
	645	zifl(i) = ziflcld(i) + ziflclr(i)
	646	ELSE
	647	zrfl(i) = zrfl(i)+ zqprecl(i) &
	648	(paprsdn(i)-paprsup(i))/(RGdtime)
	649	zifl(i) = zifl(i)+ zqpreci(i) &
	650	(paprsdn(i)-paprsup(i))/(RGdtime)
	651	ENDIF
	652	! c_iso : same for isotopes
	653
	654	ENDIF ! rneb>0
	655
	656	ENDDO
	657
	658	! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
	659	! if iflag_evap_prec>=4
	660	IF (iflag_evap_prec.GE.4) THEN
	661
	662	DO i=1,klon
	663
	664	IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
	665	znebprecipclr(i) = min(znebprecipclr(i),max( &
	666	zrflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), &
	667	ziflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
	668	ELSE
	669	znebprecipclr(i)=0.0
	670	ENDIF
	671
	672	IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
	673	znebprecipcld(i) = min(znebprecipcld(i), max( &
	674	zrflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), &
	675	ziflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
	676	ELSE
	677	znebprecipcld(i)=0.0
	678	ENDIF
	679	ENDDO
	680
	681	ENDIF
	682
	683	! we recalculate zradoice to account for contributions from the precipitation flux
	684	! if ok_radocond_snow is true
	685	IF ( ok_radocond_snow ) THEN
	686	IF ( .NOT. iftop ) THEN
	687	! for the solid phase (crystals + snowflakes)
	688	! we recalculate zradoice by summing
	689	! the ice content calculated in the mesh
	690	! + the contribution of the non-evaporated snowfall
	691	! from the overlying layer
	692	DO i=1,klon
	693	IF (ziflprev(i) .NE. 0.0) THEN
	694	zradoice(i)=zoliqi(i)+zqpreci(i)+ziflprev(i)/zrho_up(i)/zvelo_up(i)
	695	ELSE
	696	zradoice(i)=zoliqi(i)+zqpreci(i)
	697	ENDIF
	698	zradocond(i)=zradoliq(i)+zradoice(i)
	699	ENDDO
	700	ENDIF
	701	! keep in memory air density and ice fallspeed velocity
	702	zrho_up(:) = zrho(:)
	703	zvelo_up(:) = zvelo(:)
	704	ENDIF
	705
	706	END SUBROUTINE histprecip_postcld
	707
	708
[5408]	709	!----------------------------------------------------------------
	710	! Computes the processes-oriented precipitation formulations for
	711	! evaporation and sublimation
	712	!
	713	SUBROUTINE poprecip_precld( &
	714	klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
	715	qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
[5431]	716	cldfra, rvc_seri, qliq, qice, &
[5408]	717	rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
	718	)
	719
	720	USE lmdz_lscp_ini, ONLY : prt_level, lunout
	721	USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
	722	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
[5431]	723	USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
	724	USE lmdz_lscp_ini, ONLY : eps, temp_nowater
[5408]	725	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
	726
	727	IMPLICIT NONE
	728
	729
	730	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
	731	REAL, INTENT(IN) :: dtime !--time step [s]
	732	LOGICAL, INTENT(IN) :: iftop !--if top of the column
	733
	734
	735	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
	736	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
	737	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
	738
	739	REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K]
	740	REAL, INTENT(IN), DIMENSION(klon) :: tempupnew !--updated temperature of the overlying layer [K]
	741
	742	REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
	743	REAL, INTENT(INOUT), DIMENSION(klon) :: qprecip !--specific humidity in the precipitation falling from the upper layer [kg/kg]
	744
	745	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
	746	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
	747
	748	REAL, INTENT(IN), DIMENSION(klon) :: qvapclrup !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
	749	REAL, INTENT(IN), DIMENSION(klon) :: qtotupnew !--total specific humidity IN THE LAYER ABOVE [kg/kg]
	750
[5431]	751	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
	752	REAL, INTENT(IN), DIMENSION(klon) :: rvc_seri !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg]
	753	REAL, INTENT(IN), DIMENSION(klon) :: qliq !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
	754	REAL, INTENT(IN), DIMENSION(klon) :: qice !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
	755
	756
[5408]	757	REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
	758	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
[5431]	759	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
[5408]	760	REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
	761	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
[5431]	762	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
[5408]	763
	764	REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s]
	765	REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
	766
	767
	768	!--Integer for interating over klon
	769	INTEGER :: i
	770	!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
	771	REAL, DIMENSION(klon) :: dhum_to_dflux
[5431]	772	!--Air density [kg/m3] and layer thickness [m]
[5408]	773	REAL, DIMENSION(klon) :: rho, dz
[5431]	774	!--Temporary precip fractions and fluxes for the evaporation
	775	REAL, DIMENSION(klon) :: precipfracclr_tmp, precipfraccld_tmp
	776	REAL, DIMENSION(klon) :: rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp
[5408]	777
	778	!--Saturation values
	779	REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
[5431]	780	!--Vapor in the clear sky and cloud
	781	REAL :: qvapclr, qvapcld
[5408]	782	!--Fluxes tendencies because of evaporation and sublimation
[5431]	783	REAL :: dprecip_evasub_max, dprecip_evasub_tot
	784	REAL :: drainclreva, draincldeva, dsnowclrsub, dsnowcldsub
[5408]	785	!--Specific humidity tendencies because of evaporation and sublimation
	786	REAL :: dqrevap, dqssubl
	787	!--Specific heat constant
	788	REAL :: cpair, cpw
	789
	790	!--Initialisation
	791	qzero(:) = 0.
	792	dqreva(:) = 0.
	793	dqssub(:) = 0.
	794	dqrevap = 0.
	795	dqssubl = 0.
	796
	797	!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
	798	dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime
	799	rho(:) = pplay(:) / temp(:) / RD
	800	dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:)
	801
	802	!--Calculation of saturation specific humidity
	803	!--depending on temperature:
	804	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
	805	!--wrt liquid water
	806	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
	807	!--wrt ice
	808	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
	809
	810
	811
	812	!--First step consists in "thermalizing" the layer:
	813	!--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature
	814	!--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the
	815	!--layer have.
	816
	817	IF (iftop) THEN
	818
	819	DO i = 1, klon
	820	qprecip(i) = 0.
	821	ENDDO
	822
	823	ELSE
	824
	825	DO i = 1, klon
	826	!--No condensed water so cp=cp(vapor+dry air)
	827	!-- RVTMP2=rcpv/rcpd-1
	828	cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
	829	cpw = RCPD * RVTMP2
	830	!--qprecip has to be thermalized with
	831	!--layer's air so that precipitation at the ground has the
	832	!--same temperature as the lowermost layer
	833	!--we convert the flux into a specific quantity qprecip
	834	qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i)
	835	!-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer
	836	temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) &
	837	/ ( cpair + qprecip(i) * cpw )
	838	ENDDO
	839
	840	ENDIF
	841
	842
[5431]	843	!--Initialise the precipitation fractions and fluxes
	844	DO i = 1, klon
	845	precipfracclr_tmp(i) = precipfracclr(i)
	846	precipfraccld_tmp(i) = precipfraccld(i)
	847	rainclr_tmp(i) = rainclr(i)
	848	snowclr_tmp(i) = snowclr(i)
	849	raincld_tmp(i) = raincld(i)
	850	snowcld_tmp(i) = snowcld(i)
	851	ENDDO
[5408]	852
[5431]	853	!--If we relax the LTP assumption that the cloud is the same than in the
	854	!--gridbox above, we can change the tmp quantities
	855	IF ( ok_ice_supersat ) THEN
	856	!--Update the precipitation fraction using advected cloud fraction
	857	CALL poprecip_fracupdate( &
	858	klon, cldfra, precipfracclr_tmp, precipfraccld_tmp, &
	859	rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp)
	860
	861	! here we could put an ELSE statement and do one iteration of the condensation scheme
	862	! (but it would only diagnose cloud from lognormal - link with thermals?)
	863
	864	DO i = 1, klon
	865	!--We cannot have a total fraction greater than the one in the gridbox above
	866	!--because new cloud formation has not yet occured.
	867	!--If this happens, we reduce the precip frac in the cloud, because it can
	868	!--only mean that the cloud fraction at this level is greater than the total
	869	!--precipitation fraction of the level above, and the precip frac in clear sky
	870	!--is zero.
	871	!--NB. this only works because we assume a max-random cloud and precip overlap
	872	precipfraccld_tmp(i) = MIN( precipfraccld_tmp(i), precipfracclr(i) + precipfraccld(i) )
	873	ENDDO
	874
	875	ENDIF
	876	!--At this stage, we guarantee that
	877	!-- precipfracclr + precipfraccld = precipfracclr_tmp + precipfraccld_tmp
	878
	879
[5408]	880	DO i = 1, klon
	881
	882	!--If there is precipitation from the layer above
[5431]	883	IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
[5408]	884
[5431]	885	!--Init
	886	drainclreva = 0.
	887	draincldeva = 0.
	888	dsnowclrsub = 0.
	889	dsnowcldsub = 0.
	890
	891	!--Init the properties of the air where reevaporation occurs
	892	IF ( ok_ice_supersat ) THEN
	893	!--We can diagnose the water vapor in the cloud using the advected
	894	!--water vapor in the cloud and cloud fraction
	895	IF ( ok_unadjusted_clouds .AND. ( cldfra(i) .GT. eps ) ) THEN
	896	qvapcld = rvc_seri(i) * qvap(i) / cldfra(i)
	897	ENDIF
	898	!--We can diagnose completely the water vapor in clear sky, because all
	899	!--the needed variables (ice, liq, vapor in cld, cloud fraction) are
	900	!--advected
	901	!--Note that qvap(i) is the total water in the gridbox
	902	IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
	903	qvapclr = ( qvap(i) - qice(i) - qliq(i) - rvc_seri(i) * qvap(i) ) / ( 1. - cldfra(i) )
	904	ENDIF
	905	ELSEIF ( ok_corr_vap_evasub ) THEN
	906	!--Corrected version from Ludo - we use the same water ratio between
[5408]	907	!--the clear and the cloudy sky as in the layer above. This
	908	!--extends the assumption that the cloud fraction is the same
	909	!--as the layer above. This is assumed only for the evap / subl
	910	!--process
	911	!--Note that qvap(i) is the total water in the gridbox, and
	912	!--precipfraccld(i) is the cloud fraction in the layer above
[5431]	913	IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
	914	qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
	915	ENDIF
[5408]	916	ELSE
	917	!--Legacy version from Ludo - we use the total specific humidity
	918	!--for the evap / subl process
	919	qvapclr = qvap(i)
	920	ENDIF
	921
[5431]	922	!---------------------------
	923	!-- EVAP / SUBL IN CLEAR SKY
	924	!---------------------------
[5408]	925
[5431]	926	IF ( precipfracclr_tmp(i) .GT. eps ) THEN
[5408]	927
[5431]	928	!--Evaporation of liquid precipitation coming from above
	929	!--in the clear sky only
	930	!--dprecip/dz = -beta(1-qvap/qsat)(precip**expo_eva)
	931	!--formula from Sundqvist 1988, Klemp & Wilhemson 1978
	932	!--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
	933	!--which does not need a barrier on rainclr, because included in the formula
	934	drainclreva = precipfracclr_tmp(i) * MAX(0., &
	935	- coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
	936	+ ( rainclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_eva ) &
	937	)**( 1. / ( 1. - expo_eva ) ) - rainclr_tmp(i)
	938
	939	!--Evaporation is limited by 0
	940	!--NB. with ok_ice_supersat activated, this barrier should be useless
	941	drainclreva = MIN(0., drainclreva)
[5408]	942
	943
[5431]	944	!--Sublimation of the solid precipitation coming from above
	945	!--(same formula as for liquid precip)
	946	!--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
	947	!--which does not need a barrier on snowclr, because included in the formula
	948	dsnowclrsub = precipfracclr_tmp(i) * MAX(0., &
	949	- coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) &
	950	+ ( snowclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_sub ) &
	951	)**( 1. / ( 1. - expo_sub ) ) - snowclr_tmp(i)
[5408]	952
[5431]	953	!--If ice supersaturation is activated, we allow vapor to condense on falling snow
	954	!--i.e., having a positive dsnowclrsub
	955	IF ( .NOT. ok_ice_supersat ) THEN
	956	!--Sublimation is limited by 0
	957	dsnowclrsub = MIN(0., dsnowclrsub)
	958	ENDIF
[5408]	959
	960
[5431]	961	!--Evaporation limit: we ensure that the layer's fraction below
	962	!--the clear sky does not reach saturation. In this case, we
	963	!--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
	964	!--Max evaporation is computed not to saturate the clear sky precip fraction
	965	!--(i.e., the fraction where evaporation occurs)
	966	!--It is expressed as a max flux dprecip_evasub_max
	967
	968	IF ( .NOT. ok_ice_supersat ) THEN
	969	dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr_tmp(i)) &
	970	* dhum_to_dflux(i)
	971	ELSE
	972	dprecip_evasub_max = ( qvapclr - qsat(i) ) * precipfracclr_tmp(i) &
	973	* dhum_to_dflux(i)
	974	ENDIF
	975	dprecip_evasub_tot = drainclreva + dsnowclrsub
	976
	977	!--Barriers
	978	!--If activates if the total is LOWER than the max because
	979	!--everything is negative
	980	IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
	981	( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
	982	(( dprecip_evasub_max .GT. 0. ) .AND. &
	983	( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
	984	drainclreva = dprecip_evasub_max * drainclreva / dprecip_evasub_tot
	985	dsnowclrsub = dprecip_evasub_max * dsnowclrsub / dprecip_evasub_tot
	986	ENDIF
	987
	988	ENDIF ! precipfracclr_tmp .GT. eps
	989
	990
	991	!---------------------------
	992	!-- EVAP / SUBL IN THE CLOUD
	993	!---------------------------
	994
	995	IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
	996	!--Evaporation of liquid precipitation coming from above
	997	!--in the cloud only
	998	!--dprecip/dz = -beta(1-qvap/qsat)(precip**expo_eva)
	999	!--formula from Sundqvist 1988, Klemp & Wilhemson 1978
	1000	!--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
	1001	!--which does not need a barrier on raincld, because included in the formula
	1002	!draincldeva = precipfraccld_tmp(i) * MAX(0., &
	1003	! - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
	1004	! + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
	1005	! )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
	1006
	1007	!--Evaporation is limited by 0
	1008	!--NB. with ok_ice_supersat activated, this barrier should be useless
	1009	!draincldeva = MIN(0., draincldeva)
	1010	draincldeva = 0.
	1011
	1012
	1013	!--Sublimation of the solid precipitation coming from above
	1014	!--(same formula as for liquid precip)
	1015	!--Exact explicit formulation (snowcld is resolved exactly, qvap explicitly)
	1016	!--which does not need a barrier on snowcld, because included in the formula
	1017	dsnowcldsub = precipfraccld_tmp(i) * MAX(0., &
	1018	- coef_sub * ( 1. - expo_sub ) * (1. - qvapcld / qsati(i)) * dz(i) &
	1019	+ ( snowcld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_sub ) &
	1020	)**( 1. / ( 1. - expo_sub ) ) - snowcld_tmp(i)
	1021
	1022	!--There is no barrier because we want deposition to occur if it is possible
	1023
	1024
	1025	!--Evaporation limit: we ensure that the layer's fraction below
	1026	!--the clear sky does not reach saturation. In this case, we
	1027	!--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
	1028	!--Max evaporation is computed not to saturate the clear sky precip fraction
	1029	!--(i.e., the fraction where evaporation occurs)
	1030	!--It is expressed as a max flux dprecip_evasub_max
	1031
	1032	dprecip_evasub_max = ( qvapcld - qsat(i) ) * precipfraccld_tmp(i) * dhum_to_dflux(i)
	1033	dprecip_evasub_tot = draincldeva + dsnowcldsub
	1034
	1035	!--Barriers
	1036	!--If activates if the total is LOWER than the max because
	1037	!--everything is negative
	1038	IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
	1039	( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
	1040	(( dprecip_evasub_max .GT. 0. ) .AND. &
	1041	( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
	1042	draincldeva = dprecip_evasub_max * draincldeva / dprecip_evasub_tot
	1043	dsnowcldsub = dprecip_evasub_max * dsnowcldsub / dprecip_evasub_tot
	1044	ENDIF
	1045
	1046	ENDIF ! ok_unadjusted_clouds
	1047
	1048
[5408]	1049	!--New solid and liquid precipitation fluxes after evap and sublimation
[5431]	1050	dqrevap = ( drainclreva + draincldeva ) / dhum_to_dflux(i)
	1051	dqssubl = ( dsnowclrsub + dsnowcldsub ) / dhum_to_dflux(i)
[5408]	1052
	1053	!--Vapor is updated after evaporation/sublimation (it is increased)
	1054	qvap(i) = qvap(i) - dqrevap - dqssubl
	1055	!--qprecip is the total condensed water in the precip flux (it is decreased)
	1056	qprecip(i) = qprecip(i) + dqrevap + dqssubl
	1057	!--Air and precip temperature (i.e., gridbox temperature)
	1058	!--is updated due to latent heat cooling
	1059	temp(i) = temp(i) &
	1060	+ dqrevap * RLVTT / RCPD &
	1061	/ ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
	1062	+ dqssubl * RLSTT / RCPD &
	1063	/ ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
	1064
	1065	!--Add tendencies
[5431]	1066	!--The MAX is needed because in some cases, the flux can be slightly
	1067	!--negative (numerical precision)
	1068	rainclr_tmp(i) = MAX(0., rainclr_tmp(i) + drainclreva)
	1069	snowclr_tmp(i) = MAX(0., snowclr_tmp(i) + dsnowclrsub)
	1070	raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
	1071	snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)
[5408]	1072
[5431]	1073
	1074	!--We reattribute fluxes according to initial fractions, using proportionnality
	1075	!--Note that trough this process, we conserve total precip fluxes
	1076	!-- rainclr_tmp + raincld_tmp = rainclr + raincld (same for snow)
	1077	!--If the cloud has increased, a part of the precip in clear sky falls in clear sky,
	1078	!--the rest falls in the cloud
	1079	IF ( ( precipfraccld(i) .GT. eps ) .AND. ( precipfraccld_tmp(i) .GT. precipfraccld(i) ) ) THEN
	1080	!--All the precip from cloud falls in the cloud
	1081	raincld(i) = MAX(0., raincld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
	1082	rainclr(i) = MAX(0., rainclr_tmp(i) + raincld_tmp(i) - raincld(i))
	1083	snowcld(i) = MAX(0., snowcld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
	1084	snowclr(i) = MAX(0., snowclr_tmp(i) + snowcld_tmp(i) - snowcld(i))
	1085
	1086	!--If the cloud has narrowed, a part of the precip in cloud falls in clear sky,
	1087	!--the rest falls in the cloud
	1088	ELSEIF ( ( precipfracclr(i) .GT. eps ) .AND. ( precipfracclr_tmp(i) .GT. precipfracclr(i) ) ) THEN
	1089	!--All the precip from clear sky falls in clear sky
	1090	rainclr(i) = MAX(0., rainclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
	1091	raincld(i) = MAX(0., raincld_tmp(i) + rainclr_tmp(i) - rainclr(i))
	1092	snowclr(i) = MAX(0., snowclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
	1093	snowcld(i) = MAX(0., snowcld_tmp(i) + snowclr_tmp(i) - snowclr(i))
	1094
	1095	ELSE
	1096	!--If the cloud stays the same or if there is no cloud above and
	1097	!--in the current layer, there is no reattribution
	1098	rainclr(i) = rainclr_tmp(i)
	1099	raincld(i) = raincld_tmp(i)
	1100	snowclr(i) = snowclr_tmp(i)
	1101	snowcld(i) = snowcld_tmp(i)
	1102	ENDIF
	1103
	1104	!--If there is no more precip fluxes, the precipitation fraction is set to 0
[5408]	1105	IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
[5431]	1106	IF ( ( raincld(i) + snowcld(i) ) .LE. 0. ) precipfraccld(i) = 0.
[5408]	1107
	1108	!--Calculation of the total fluxes
	1109	rain(i) = rainclr(i) + raincld(i)
	1110	snow(i) = snowclr(i) + snowcld(i)
	1111
	1112	ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
	1113	!--If no precip, we reinitialize the cloud fraction used for the precip to 0
	1114	precipfraccld(i) = 0.
	1115	precipfracclr(i) = 0.
	1116
	1117	ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
	1118
	1119	!--Diagnostic tendencies
	1120	dqssub(i) = dqssubl / dtime
	1121	dqreva(i) = dqrevap / dtime
	1122
	1123	ENDDO ! loop on klon
	1124
	1125	END SUBROUTINE poprecip_precld
	1126
	1127
	1128	!----------------------------------------------------------------
[5431]	1129	! Computes the precipitation fraction update
	1130	! The goal of this routine is to reattribute precipitation fractions
	1131	! and fluxes to clear or cloudy air, depending on the variation of
	1132	! the cloud fraction on the vertical dimension. We assume a
	1133	! maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
	1134	! and LTP thesis, 2021)
	1135	! NB. in fact, we assume a maximum-random overlap of the total precip. frac
	1136	!
	1137	SUBROUTINE poprecip_fracupdate( &
	1138	klon, cldfra, precipfracclr, precipfraccld, &
	1139	rainclr, raincld, snowclr, snowcld)
	1140
	1141	USE lmdz_lscp_ini, ONLY : eps
	1142
	1143	IMPLICIT NONE
	1144
	1145	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
	1146
	1147	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-]
	1148	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
	1149	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
	1150	!--NB. at the end of the routine, becomes the fraction of precip
	1151	!--in the current layer
	1152
	1153	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
	1154	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
	1155	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
	1156	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
	1157
	1158	!--Local variables
	1159	INTEGER :: i
	1160	REAL :: dcldfra
	1161	REAL :: precipfractot
	1162	REAL :: dprecipfracclr, dprecipfraccld
	1163	REAL :: drainclr, dsnowclr
	1164	REAL :: draincld, dsnowcld
	1165
	1166
	1167	DO i = 1, klon
	1168
	1169	!--Initialisation
	1170	precipfractot = precipfracclr(i) + precipfraccld(i)
	1171
	1172	!--Instead of using the cloud cover which was use in LTP thesis, we use the
	1173	!--total precip. fraction to compute the maximum-random overlap. This is
	1174	!--because all the information of the cloud cover is embedded into
	1175	!--precipfractot, and this allows for taking into account the potential
	1176	!--reduction of the precipitation fraction because either the flux is too
	1177	!--small (see barrier at the end of poprecip_postcld) or the flux is completely
	1178	!--evaporated (see barrier at the end of poprecip_precld)
	1179	!--NB. precipfraccld(i) is here the cloud fraction of the layer above
	1180	!precipfractot = 1. - ( 1. - precipfractot ) * &
	1181	! ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
	1182	! / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
	1183
	1184
	1185	IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
	1186	precipfractot = 1.
	1187	ELSE
	1188	precipfractot = 1. - ( 1. - precipfractot ) * &
	1189	( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
	1190	/ ( 1. - precipfraccld(i) )
	1191	ENDIF
	1192
	1193	!--precipfraccld(i) is here the cloud fraction of the layer above
	1194	dcldfra = cldfra(i) - precipfraccld(i)
	1195	!--Tendency of the clear-sky precipitation fraction. We add a MAX on the
	1196	!--calculation of the current CS precip. frac.
	1197	!dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
	1198	!--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
	1199	!--if precipfractot < cldfra)
	1200	dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
	1201	!--Tendency of the cloudy precipitation fraction. We add a MAX on the
	1202	!--calculation of the current CS precip. frac.
	1203	!dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
	1204	!--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
	1205	!--if cldfra < 0)
	1206	dprecipfraccld = dcldfra
	1207
	1208
	1209	!--If the cloud extends
	1210	IF ( dprecipfraccld .GT. 0. ) THEN
	1211	!--If there is no CS precip, nothing happens.
	1212	!--If there is, we reattribute some of the CS precip flux
	1213	!--to the cloud precip flux, proportionnally to the
	1214	!--decrease of the CS precip fraction
	1215	IF ( precipfracclr(i) .LT. eps ) THEN
	1216	drainclr = 0.
	1217	dsnowclr = 0.
	1218	ELSE
	1219	drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i)
	1220	dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i)
	1221	ENDIF
	1222	!--If the cloud narrows
	1223	ELSEIF ( dprecipfraccld .LT. 0. ) THEN
	1224	!--We reattribute some of the cloudy precip flux
	1225	!--to the CS precip flux, proportionnally to the
	1226	!--decrease of the cloud precip fraction
	1227	IF ( precipfraccld(i) .LT. eps ) THEN
	1228	draincld = 0.
	1229	dsnowcld = 0.
	1230	ELSE
	1231	draincld = dprecipfraccld / precipfraccld(i) * raincld(i)
	1232	dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i)
	1233	ENDIF
	1234	drainclr = - draincld
	1235	dsnowclr = - dsnowcld
	1236	!--If the cloud stays the same or if there is no cloud above and
	1237	!--in the current layer, nothing happens
	1238	ELSE
	1239	drainclr = 0.
	1240	dsnowclr = 0.
	1241	ENDIF
	1242
	1243	!--We add the tendencies
	1244	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
	1245	precipfraccld(i) = precipfraccld(i) + dprecipfraccld
	1246	precipfracclr(i) = precipfracclr(i) + dprecipfracclr
	1247	rainclr(i) = MAX(0., rainclr(i) + drainclr)
	1248	snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
	1249	raincld(i) = MAX(0., raincld(i) - drainclr)
	1250	snowcld(i) = MAX(0., snowcld(i) - dsnowclr)
	1251
	1252	ENDDO
	1253
	1254	END SUBROUTINE poprecip_fracupdate
	1255
	1256
	1257	!----------------------------------------------------------------
[5408]	1258	! Computes the processes-oriented precipitation formulations for
	1259	! - autoconversion (auto) via a deposition process
	1260	! - aggregation (agg)
	1261	! - riming (rim)
	1262	! - collection (col)
	1263	! - melting (melt)
	1264	! - freezing (freez)
	1265	!
	1266	SUBROUTINE poprecip_postcld( &
	1267	klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
	1268	temp, qvap, qliq, qice, icefrac, cldfra, &
	1269	precipfracclr, precipfraccld, &
	1270	rain, rainclr, raincld, snow, snowclr, snowcld, &
	1271	qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
	1272	dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
	1273
	1274	USE lmdz_lscp_ini, ONLY : prt_level, lunout
	1275	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
	1276	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
	1277
	1278	USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb, &
	1279	cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
	1280	thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
	1281	rho_rain, r_rain, r_snow, rho_ice, &
	1282	tau_auto_snow_min, tau_auto_snow_max, &
	1283	thresh_precip_frac, eps, &
	1284	gamma_melt, alpha_freez, beta_freez, temp_nowater, &
	1285	iflag_cloudth_vert, iflag_rain_incloud_vol, &
	1286	cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez, &
	1287	rain_fallspeed_clr, rain_fallspeed_cld, &
	1288	snow_fallspeed_clr, snow_fallspeed_cld
	1289
	1290
	1291	IMPLICIT NONE
	1292
	1293	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-]
	1294	REAL, INTENT(IN) :: dtime !--time step [s]
	1295
	1296	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
	1297	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
	1298	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
	1299
	1300	REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-]
	1301	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point
	1302
	1303	REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K]
	1304	REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity [kg/kg]
	1305	REAL, INTENT(INOUT), DIMENSION(klon) :: qliq !--current liquid water specific humidity [kg/kg]
	1306	REAL, INTENT(INOUT), DIMENSION(klon) :: qice !--current ice water specific humidity [kg/kg]
	1307	REAL, INTENT(IN), DIMENSION(klon) :: icefrac !--ice fraction [-]
	1308	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-]
	1309
	1310	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
	1311	REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
	1312	!--NB. at the end of the routine, becomes the fraction of precip
	1313	!--in the current layer
	1314
	1315	REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
	1316	REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
	1317	REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
	1318	REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
	1319	REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
	1320	REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
	1321
	1322	REAL, INTENT(OUT), DIMENSION(klon) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg]
	1323	REAL, INTENT(OUT), DIMENSION(klon) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg]
	1324	REAL, INTENT(OUT), DIMENSION(klon) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
	1325	REAL, INTENT(OUT), DIMENSION(klon) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
	1326	REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
	1327	REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
	1328	REAL, INTENT(OUT), DIMENSION(klon) :: dqsrim !--snow tendency due to riming [kg/kg/s]
	1329	REAL, INTENT(OUT), DIMENSION(klon) :: dqsmelt !--snow tendency due to melting [kg/kg/s]
	1330	REAL, INTENT(OUT), DIMENSION(klon) :: dqrmelt !--rain tendency due to melting [kg/kg/s]
	1331	REAL, INTENT(OUT), DIMENSION(klon) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
	1332	REAL, INTENT(OUT), DIMENSION(klon) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
	1333
	1334
	1335
	1336	!--Local variables
	1337
	1338	INTEGER :: i
	1339	!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
	1340	REAL, DIMENSION(klon) :: dhum_to_dflux
	1341	REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
	1342
	1343	!--Collection, aggregation and riming
	1344	REAL :: eff_cldfra
	1345	REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
	1346	REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
	1347	REAL :: rho_snow
	1348	REAL :: dqlcol !--loss of liquid cloud content due to collection by rain [kg/kg/s]
	1349	REAL :: dqiagg !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
	1350	REAL :: dqlrim !--loss of liquid cloud content due to riming on snow [kg/kg/s]
	1351
	1352	!--Autoconversion
	1353	REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
	1354	REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
	1355	REAL :: dqlauto !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
	1356	REAL :: dqiauto !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
	1357
	1358	!--Melting
	1359	REAL :: dqsmelt_max, air_thermal_conduct
	1360	REAL :: nb_snowflake_clr, nb_snowflake_cld
	1361	REAL :: capa_snowflake, temp_wetbulb
	1362	REAL :: rho, r_ice
	1363	REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
	1364	REAL, DIMENSION(klon) :: qzero, qsat, dqsat
	1365
	1366	!--Freezing
	1367	REAL :: dqrfreez_max
	1368	REAL :: tau_freez
	1369	REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
	1370	REAL :: coef_freez
	1371	REAL :: dqifreez !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
	1372	REAL :: Eff_rain_ice
	1373
	1374
	1375	!--Initialisation of variables
	1376
	1377
	1378	qzero(:) = 0.
	1379
	1380	dqrcol(:) = 0.
	1381	dqsagg(:) = 0.
	1382	dqsauto(:) = 0.
	1383	dqrauto(:) = 0.
	1384	dqsrim(:) = 0.
	1385	dqrmelt(:) = 0.
	1386	dqsmelt(:) = 0.
	1387	dqrfreez(:) = 0.
	1388	dqsfreez(:) = 0.
	1389
	1390
[5431]	1391	!--Update the precipitation fraction following cloud formation
	1392	CALL poprecip_fracupdate( &
	1393	klon, cldfra, precipfracclr, precipfraccld, &
	1394	rainclr, raincld, snowclr, snowcld)
	1395
	1396
[5408]	1397	DO i = 1, klon
	1398
	1399	!--Variables initialisation
	1400	dqlcol = 0.
	1401	dqiagg = 0.
	1402	dqiauto = 0.
	1403	dqlauto = 0.
	1404	dqlrim = 0.
	1405
	1406	!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
	1407	dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
	1408	qtot(i) = qvap(i) + qliq(i) + qice(i) &
	1409	+ ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
	1410
	1411	!--If vertical heterogeneity is taken into account, we use
	1412	!--the "true" volume fraction instead of a modified
	1413	!--surface fraction (which is larger and artificially
	1414	!--reduces the in-cloud water).
	1415	IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
	1416	eff_cldfra = ctot_vol(i)
	1417	ELSE
	1418	eff_cldfra = cldfra(i)
	1419	ENDIF
	1420
	1421
	1422	!--Start precipitation growth processes
	1423
	1424	!--If the cloud is big enough, the precipitation processes activate
	1425	! TODO met on seuil_neb ici ?
	1426	IF ( cldfra(i) .GE. seuil_neb ) THEN
	1427
	1428	!---------------------------------------------------------
	1429	!-- COLLECTION AND AGGREGATION
	1430	!---------------------------------------------------------
	1431	!--Collection: processus through which rain collects small liquid droplets
	1432	!--in suspension, and add it to the rain flux
	1433	!--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
	1434	!--Those processes are treated before autoconversion because we do not
	1435	!--want to collect/aggregate the newly formed fluxes, which already
	1436	!--"saw" the cloud as they come from it
	1437	!--The formulas come from Muench and Lohmann 2020
	1438
	1439	!--gamma_col: tuning coefficient [-]
	1440	!--rho_rain: volumic mass of rain [kg/m3]
	1441	!--r_rain: size of the rain droplets [m]
	1442	!--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
	1443	!--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
	1444	!--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
	1445	!--get in-cloud mean quantities. The two divisions by eff_cldfra are
	1446	!--then simplified.
	1447
	1448	!--The collection efficiency is perfect.
	1449	Eff_rain_liq = 1.
	1450	coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
	1451	IF ( raincld(i) .GT. 0. ) THEN
	1452	!--Exact explicit version, which does not need a barrier because of
	1453	!--the exponential decrease
	1454	dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
	1455
	1456	!--Add tendencies
	1457	qliq(i) = qliq(i) + dqlcol
	1458	raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
	1459
	1460	!--Diagnostic tendencies
	1461	dqrcol(i) = - dqlcol / dtime
	1462	ENDIF
	1463
	1464	!--Same as for aggregation
	1465	!--Eff_snow_liq formula:
	1466	!--it s a product of a collection efficiency and a sticking efficiency
	1467	! Milbrandt and Yau formula that gives very low values:
	1468	! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
	1469	! Lin 1983's formula
	1470	Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) )
	1471	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
	1472	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
	1473	coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
	1474	IF ( snowcld(i) .GT. 0. ) THEN
	1475	!--Exact explicit version, which does not need a barrier because of
	1476	!--the exponential decrease
	1477	dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
	1478
	1479	!--Add tendencies
	1480	qice(i) = qice(i) + dqiagg
	1481	snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
	1482
	1483	!--Diagnostic tendencies
	1484	dqsagg(i) = - dqiagg / dtime
	1485	ENDIF
	1486
	1487
	1488	!---------------------------------------------------------
	1489	!-- AUTOCONVERSION
	1490	!---------------------------------------------------------
	1491	!--Autoconversion converts liquid droplets/ice crystals into
	1492	!--rain drops/snowflakes. It relies on the formulations by
	1493	!--Sundqvist 1978.
	1494
	1495	!--If we are in a convective point, we have different parameters
	1496	!--for the autoconversion
	1497	IF ( ptconv(i) ) THEN
	1498	qthresh_auto_rain = cld_lc_con
	1499	qthresh_auto_snow = cld_lc_con_snow
	1500
	1501	tau_auto_rain = cld_tau_con
	1502	!--tau for snow depends on the ice fraction in mixed-phase clouds
	1503	tau_auto_snow = tau_auto_snow_max &
	1504	+ ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
	1505
	1506	expo_auto_rain = cld_expo_con
	1507	expo_auto_snow = cld_expo_con
	1508	ELSE
	1509	qthresh_auto_rain = cld_lc_lsc
	1510	qthresh_auto_snow = cld_lc_lsc_snow
	1511
	1512	tau_auto_rain = cld_tau_lsc
	1513	!--tau for snow depends on the ice fraction in mixed-phase clouds
	1514	tau_auto_snow = tau_auto_snow_max &
	1515	+ ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
	1516
	1517	expo_auto_rain = cld_expo_lsc
	1518	expo_auto_snow = cld_expo_lsc
	1519	ENDIF
	1520
	1521
	1522	! Liquid water quantity to remove according to (Sundqvist, 1978)
	1523	! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
	1524	!
	1525	!--And same formula for ice
	1526	!
	1527	!--We first treat the second term (with exponential) in an explicit way
	1528	!--and then treat the first term (-q/tau) in an exact way
	1529
	1530	dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
	1531	- ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
	1532
	1533	dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
	1534	- ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
	1535
	1536
	1537	!--Barriers so that we don't create more rain/snow
	1538	!--than there is liquid/ice
	1539	dqlauto = MAX( - qliq(i), dqlauto )
	1540	dqiauto = MAX( - qice(i), dqiauto )
	1541
	1542	!--Add tendencies
	1543	qliq(i) = qliq(i) + dqlauto
	1544	qice(i) = qice(i) + dqiauto
	1545	raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
	1546	snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
	1547
	1548	!--Diagnostic tendencies
	1549	dqsauto(i) = - dqiauto / dtime
	1550	dqrauto(i) = - dqlauto / dtime
	1551
	1552
	1553	!---------------------------------------------------------
	1554	!-- RIMING
	1555	!---------------------------------------------------------
	1556	!--Process which converts liquid droplets in suspension into
	1557	!--snow because of the collision between
	1558	!--those and falling snowflakes.
	1559	!--The formula comes from Muench and Lohmann 2020
	1560	!--NB.: this process needs a temperature adjustment
	1561
	1562	!--Eff_snow_liq formula: following Ferrier 1994,
	1563	!--assuming 1
	1564	Eff_snow_liq = 1.0
	1565	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
	1566	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
	1567	coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
	1568	IF ( snowcld(i) .GT. 0. ) THEN
	1569	!--Exact version, which does not need a barrier because of
	1570	!--the exponential decrease
	1571	dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
	1572
	1573	!--Add tendencies
	1574	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
	1575	qliq(i) = qliq(i) + dqlrim
	1576	snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
	1577
	1578	!--Temperature adjustment with the release of latent
	1579	!--heat because of solid condensation
	1580	temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
	1581	/ ( 1. + RVTMP2 * qtot(i) )
	1582
	1583	!--Diagnostic tendencies
	1584	dqsrim(i) = - dqlrim / dtime
	1585	ENDIF
	1586
	1587	ENDIF ! cldfra .GE. seuil_neb
	1588
	1589	ENDDO ! loop on klon
	1590
	1591
	1592	!--Re-calculation of saturation specific humidity
	1593	!--because riming changed temperature
	1594	CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
	1595
	1596	DO i = 1, klon
	1597
	1598	!---------------------------------------------------------
	1599	!-- MELTING
	1600	!---------------------------------------------------------
	1601	!--Process through which snow melts into rain.
	1602	!--The formula is homemade.
	1603	!--NB.: this process needs a temperature adjustment
	1604
	1605	!--dqsmelt_max : maximum snow melting so that temperature
	1606	!-- stays higher than 273 K [kg/kg]
	1607	!--capa_snowflake : capacitance of a snowflake, equal to
	1608	!-- the radius if the snowflake is a sphere [m]
	1609	!--temp_wetbulb : wet-bulb temperature [K]
	1610	!--snow_fallspeed : snow fall velocity (in clear/cloudy sky) [m/s]
	1611	!--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
	1612	!--gamma_melt : tuning parameter for melting [-]
	1613	!--nb_snowflake : number of snowflakes (in clear/cloudy air) [-]
	1614
	1615	IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
	1616	!--Computed according to
	1617	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
	1618	dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
	1619	* ( 1. + RVTMP2 * qtot(i) ))
	1620
	1621	!--Initialisation
	1622	dqsclrmelt = 0.
	1623	dqscldmelt = 0.
	1624
	1625	!--We assume that the snowflakes are spherical
	1626	capa_snowflake = r_snow
	1627	!--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
	1628	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
	1629	!--rho_snow formula follows Brandes et al. 2007 (JAMC)
	1630	rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
	1631
	1632	!--In clear air
	1633	IF ( ( snowclr(i) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN
	1634	!--Formula for the wet-bulb temperature from ECMWF (IFS)
	1635	!--The vapor used is the vapor in the clear sky
	1636	temp_wetbulb = temp(i) &
	1637	- ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
	1638	* ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
	1639	- 40.637 * ( temp(i) - 275. ) )
	1640	!--Calculated according to
	1641	!-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
	1642	nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
	1643	/ ( 4. / 3. * RPI * r_snow*3. rho_snow )
	1644	dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
	1645	* capa_snowflake / RLMLT * gamma_melt &
	1646	* MAX(0., temp_wetbulb - RTT) * dtime
	1647
	1648	!--Barrier to limit the melting flux to the clr snow flux in the mesh
	1649	dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
	1650	ENDIF
	1651
	1652
	1653	!--In cloudy air
	1654	IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
	1655	!--As the air is saturated, the wet-bulb temperature is equal to the
	1656	!--temperature
	1657	temp_wetbulb = temp(i)
	1658	!--Calculated according to
	1659	!-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
	1660	nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
	1661	/ ( 4. / 3. * RPI * r_snow*3. rho_snow )
	1662	dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
	1663	* capa_snowflake / RLMLT * gamma_melt &
	1664	* MAX(0., temp_wetbulb - RTT) * dtime
	1665
	1666	!--Barrier to limit the melting flux to the cld snow flux in the mesh
	1667	dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
	1668	ENDIF
	1669
	1670
	1671	!--Barrier on temperature. If the total melting flux leads to a
	1672	!--positive temperature, it is limited to keep temperature above 0 degC.
	1673	!--It is activated if the total is LOWER than the max
	1674	!--because everything is negative
	1675	dqstotmelt = dqsclrmelt + dqscldmelt
	1676	IF ( dqstotmelt .LT. dqsmelt_max ) THEN
	1677	!--We redistribute the max melted snow keeping
	1678	!--the clear/cloud partition of the melted snow
	1679	dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
	1680	dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
	1681	dqstotmelt = dqsmelt_max
	1682
	1683	ENDIF
	1684
	1685	!--Add tendencies
	1686	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
	1687	rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
	1688	raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
	1689	snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
	1690	snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
	1691
	1692	!--Temperature adjustment with the release of latent
	1693	!--heat because of melting
	1694	temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
	1695	/ ( 1. + RVTMP2 * qtot(i) )
	1696
	1697	!--Diagnostic tendencies
	1698	dqrmelt(i) = - dqstotmelt / dtime
	1699	dqsmelt(i) = dqstotmelt / dtime
	1700
	1701	ENDIF
	1702
	1703
	1704	!---------------------------------------------------------
	1705	!-- FREEZING
	1706	!---------------------------------------------------------
	1707	!--Process through which rain freezes into snow.
	1708	!-- We parameterize it as a 2 step process:
	1709	!--first: freezing following collision with ice crystals
	1710	!--second: immersion freezing following (inspired by Bigg 1953)
	1711	!--the latter is parameterized as an exponential decrease of the rain
	1712	!--water content with a homemade formulya
	1713	!--This is based on a caracteritic time of freezing, which
	1714	!--exponentially depends on temperature so that it is
	1715	!--equal to 1 for temp_nowater (see below) and is close to
	1716	!--0 for RTT (=273.15 K).
	1717	!--NB.: this process needs a temperature adjustment
	1718	!--dqrfreez_max : maximum rain freezing so that temperature
	1719	!-- stays lower than 273 K [kg/kg]
	1720	!--tau_freez : caracteristic time of freezing [s]
	1721	!--gamma_freez : tuning parameter [s-1]
	1722	!--alpha_freez : tuning parameter for the shape of the exponential curve [-]
	1723	!--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
	1724
	1725	IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
	1726
	1727
	1728	!--1st step: freezing following collision with ice crystals
	1729	!--Sub-process of freezing which quantifies the collision between
	1730	!--ice crystals in suspension and falling rain droplets.
	1731	!--The rain droplets freeze, becoming graupel, and carrying
	1732	!--the ice crystal (which acted as an ice nucleating particle).
	1733	!--The formula is adapted from the riming formula.
	1734	!--it works only in the cloudy part
	1735
	1736	dqifreez = 0.
	1737	dqrtotfreez_step1 = 0.
	1738
	1739	IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
	1740	( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
	1741	dqrclrfreez = 0.
	1742	dqrcldfreez = 0.
	1743
	1744	!--Computed according to
	1745	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
	1746	dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
	1747	* ( 1. + RVTMP2 * qtot(i) ))
	1748
	1749
	1750	!--The collision efficiency is assumed unity
	1751	Eff_rain_ice = 1.
	1752	coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
	1753	!--Exact version, which does not need a barrier because of
	1754	!--the exponential decrease.
	1755	dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
	1756
	1757	!--We add the part of rain water that freezes, limited by a temperature barrier
	1758	!--This quantity is calculated assuming that the number of drop that freeze correspond to the number
	1759	!--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
	1760	!--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
	1761	!--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
	1762	!--The assumption above corresponds to dNi = dNr, i.e.,
	1763	!-- dqr = dqi * (4/3 PI rho_rain * r_rain*3) / (4/3 PI rho_ice r_ice**3)
	1764	!--Dry density [kg/m3]
	1765	rho = pplay(i) / temp(i) / RD
	1766	!--r_ice formula from Sun and Rikus (1999)
	1767	r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
	1768	+ 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )*0.2535 ( temp(i) - RTT + 190. ) ) / 2.
	1769	dqrcldfreez = dqifreez * rho_rain * r_rain*3. / ( rho_ice r_ice**3. )
	1770	dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
	1771	dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
	1772	dqrtotfreez_step1 = dqrcldfreez
	1773
	1774	!--Add tendencies
	1775	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
	1776	qice(i) = qice(i) + dqifreez
	1777	raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
	1778	snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
	1779	temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
	1780	/ ( 1. + RVTMP2 * qtot(i) )
	1781
	1782	ENDIF
	1783
	1784	!-- Second step immersion freezing of rain drops
	1785	!-- with a homemade timeconstant depending on temperature
	1786
	1787	dqrclrfreez = 0.
	1788	dqrcldfreez = 0.
	1789	dqrtotfreez_step2 = 0.
	1790	!--Computed according to
	1791	!--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
	1792
	1793	dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
	1794	* ( 1. + RVTMP2 * qtot(i) ))
	1795
	1796
	1797	tau_freez = 1. / ( beta_freez &
	1798	* EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
	1799
	1800
	1801	!--In clear air
	1802	IF ( rainclr(i) .GT. 0. ) THEN
	1803	!--Exact solution of dqrain/dt = -qrain/tau_freez
	1804	dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
	1805	ENDIF
	1806
	1807	!--In cloudy air
	1808	IF ( raincld(i) .GT. 0. ) THEN
	1809	!--Exact solution of dqrain/dt = -qrain/tau_freez
	1810	dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
	1811	ENDIF
	1812
	1813	!--temperature barrier step 2
	1814	!--It is activated if the total is LOWER than the max
	1815	!--because everything is negative
	1816	dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
	1817
	1818	IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
	1819	!--We redistribute the max freezed rain keeping
	1820	!--the clear/cloud partition of the freezing rain
	1821	dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
	1822	dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
	1823	dqrtotfreez_step2 = dqrfreez_max
	1824	ENDIF
	1825
	1826
	1827	!--Add tendencies
	1828	!--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
	1829	rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
	1830	raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
	1831	snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
	1832	snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
	1833
	1834
	1835	!--Temperature adjustment with the uptake of latent
	1836	!--heat because of freezing
	1837	temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
	1838	/ ( 1. + RVTMP2 * qtot(i) )
	1839
	1840	!--Diagnostic tendencies
	1841	dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2
	1842	dqrfreez(i) = dqrtotfreez / dtime
	1843	dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
	1844
	1845	ENDIF
	1846
	1847
	1848
	1849	!--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
	1850	!--we reduce the precipiration fraction in the clear/cloudy air so that the new
	1851	!--local flux of rain+snow is equal to rain_int_min.
	1852	!--Here, rain+snow is the gridbox-mean flux of precip.
	1853	!--Therefore, (rain+snow)/precipfrac is the local flux of precip.
	1854	!--If the local flux of precip is lower than rain_int_min, i.e.,
	1855	!-- (rain+snow)/precipfrac < rain_int_min , i.e.,
	1856	!-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
	1857	!--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
	1858	!--Note that this is physically different than what is proposed in LTP thesis.
	1859	precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
	1860	precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
	1861
	1862	!--Calculate outputs
	1863	rain(i) = rainclr(i) + raincld(i)
	1864	snow(i) = snowclr(i) + snowcld(i)
	1865
	1866	!--Diagnostics
	1867	!--BEWARE this is indeed a diagnostic: this is an estimation from
	1868	!--the value of the flux at the bottom interface of the mesh and
	1869	!--and assuming an upstream numerical calculation
	1870	!--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
	1871	!--used for computing the total ice water content in the mesh
	1872	!--for radiation only
	1873	qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
	1874	+ raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
	1875	qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
	1876	+ snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
	1877
	1878
	1879	ENDDO ! loop on klon
	1880
	1881	END SUBROUTINE poprecip_postcld
	1882
	1883	END MODULE lmdz_lscp_precip

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

Download in other formats:

Original Format