Context Navigation

lmdz_cloudth.F90 @ 5449

Last change on this file since 5449 was 5226, checked in by abarral, 4 months ago
Merge r5208 r5209
File size: 92.7 KB

Rev	Line
[4651]	1	MODULE lmdz_cloudth
[2686]	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[5144]	7	SUBROUTINE cloudth(ngrid, klev, ind2, &
	8	ztv, po, zqta, fraca, &
	9	qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, &
	10	ratqs, zqs, t, &
	11	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	12
[5144]	13	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs
	14	USE lmdz_yoethf
[5153]	15
[5144]	16	USE lmdz_yomcst
[2686]	17
[5144]	18	IMPLICIT NONE
[5153]	19	INCLUDE "FCTTRE.h"
[4535]	20
[2686]	21
[5144]	22	!===========================================================================
	23	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	24	! Date : 25 Mai 2010
	25	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	26	!===========================================================================
[2686]	27
[5144]	28	INTEGER itap, ind1, ind2
	29	INTEGER ngrid, klev, klon, l, ig
	30	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
[2686]	31
[5144]	32	REAL ztv(ngrid, klev)
	33	REAL po(ngrid)
	34	REAL zqenv(ngrid)
	35	REAL zqta(ngrid, klev)
[2686]	36
[5144]	37	REAL fraca(ngrid, klev + 1)
	38	REAL zpspsk(ngrid, klev)
	39	REAL paprs(ngrid, klev + 1)
	40	REAL pplay(ngrid, klev)
	41	REAL ztla(ngrid, klev)
	42	REAL zthl(ngrid, klev)
[2686]	43
[5144]	44	REAL zqsatth(ngrid, klev)
	45	REAL zqsatenv(ngrid, klev)
[2686]	46
[5144]	47	REAL sigma1(ngrid, klev)
	48	REAL sigma2(ngrid, klev)
	49	REAL qlth(ngrid, klev)
	50	REAL qlenv(ngrid, klev)
	51	REAL qltot(ngrid, klev)
	52	REAL cth(ngrid, klev)
	53	REAL cenv(ngrid, klev)
	54	REAL ctot(ngrid, klev)
	55	REAL rneb(ngrid, klev)
	56	REAL t(ngrid, klev)
	57	REAL qsatmmussig1, qsatmmussig2, sqrt2pi, pi
	58	REAL rdd, cppd, Lv
	59	REAL alth, alenv, ath, aenv
	60	REAL sth, senv, sigma1s, sigma2s, xth, xenv
	61	REAL Tbef, zdelta, qsatbef, zcor
	62	REAL qlbef
	63	REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur
[2686]	64
[5144]	65	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	66	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	67	REAL zqs(ngrid), qcloud(ngrid)
	68	REAL erf
[2686]	69
	70
	71
	72
[5144]	73	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	74	! Gestion de deux versions de cloudth
	75	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	76
	77	IF (iflag_cloudth_vert>=1) THEN
	78	CALL cloudth_vert(ngrid, klev, ind2, &
	79	ztv, po, zqta, fraca, &
	80	qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, &
	81	ratqs, zqs, t)
[2686]	82	RETURN
[5144]	83	ENDIF
	84	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2686]	85
	86
[5144]	87	!-------------------------------------------------------------------------------
	88	! Initialisation des variables r?elles
	89	!-------------------------------------------------------------------------------
	90	sigma1(:, :) = 0.
	91	sigma2(:, :) = 0.
	92	qlth(:, :) = 0.
	93	qlenv(:, :) = 0.
	94	qltot(:, :) = 0.
	95	rneb(:, :) = 0.
	96	qcloud(:) = 0.
	97	cth(:, :) = 0.
	98	cenv(:, :) = 0.
	99	ctot(:, :) = 0.
	100	qsatmmussig1 = 0.
	101	qsatmmussig2 = 0.
	102	rdd = 287.04
	103	cppd = 1005.7
	104	pi = 3.14159
	105	Lv = 2.5e6
	106	sqrt2pi = sqrt(2. * pi)
[2686]	107
	108
	109
[5144]	110	!-------------------------------------------------------------------------------
	111	! Calcul de la fraction du thermique et des ?cart-types des distributions
	112	!-------------------------------------------------------------------------------
[5158]	113	DO ind1 = 1, ngrid
[2686]	114
[5144]	115	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
	116	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2))
[2686]	117
	118
[5144]	119	! zqenv(ind1)=po(ind1)
	120	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	121	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	122	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	123	qsatbef = MIN(0.5, qsatbef)
	124	zcor = 1. / (1. - retv * qsatbef)
	125	qsatbef = qsatbef * zcor
	126	zqsatenv(ind1, ind2) = qsatbef
[2686]	127
[5144]	128	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	129	aenv = 1. / (1. + (alenv * Lv / cppd))
	130	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	131
[5144]	132	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	133	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	134	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	135	qsatbef = MIN(0.5, qsatbef)
	136	zcor = 1. / (1. - retv * qsatbef)
	137	qsatbef = qsatbef * zcor
	138	zqsatth(ind1, ind2) = qsatbef
[2686]	139
[5144]	140	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2)
	141	ath = 1. / (1. + (alth * Lv / cppd))
	142	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2))
[2686]	143
	144
	145
[5144]	146	!------------------------------------------------------------------------------
	147	! Calcul des ?cart-types pour s
	148	!------------------------------------------------------------------------------
[2686]	149
[5144]	150	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	151	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
	152	! if (paprs(ind1,ind2).gt.90000) THEN
	153	! ratqs(ind1,ind2)=0.002
	154	! else
	155	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	156	! endif
	157	sigma1s = (1.1*0.5) (fraca(ind1, ind2)*0.6) / (1 - fraca(ind1, ind2)) ((sth - senv)2)0.5 + 0.002 * po(ind1)
	158	sigma2s = 0.11 * ((sth - senv)2)0.5 / (fraca(ind1, ind2) + 0.01)*0.4 + 0.002 zqta(ind1, ind2)
	159	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	160	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[2686]	161
[5144]	162	!------------------------------------------------------------------------------
	163	! Calcul de l'eau condens?e et de la couverture nuageuse
	164	!------------------------------------------------------------------------------
	165	sqrt2pi = sqrt(2. * pi)
	166	xth = sth / (sqrt(2.) * sigma2s)
	167	xenv = senv / (sqrt(2.) * sigma1s)
	168	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	169	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	170	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	171
[5144]	172	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt(2.) * cth(ind1, ind2))
	173	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
	174	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	175
[5144]	176	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	177	IF (ctot(ind1, ind2)<1.e-10) THEN
	178	ctot(ind1, ind2) = 0.
	179	qcloud(ind1) = zqsatenv(ind1, ind2)
[2686]	180
[5144]	181	else
[2686]	182
[5144]	183	ctot(ind1, ind2) = ctot(ind1, ind2)
	184	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
[2686]	185
[5144]	186	endif
[2686]	187
	188
[5144]	189	! PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
[2686]	190
	191	else ! gaussienne environnement seule
	192
[5144]	193	zqenv(ind1) = po(ind1)
	194	Tbef = t(ind1, ind2)
	195	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	196	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	197	qsatbef = MIN(0.5, qsatbef)
	198	zcor = 1. / (1. - retv * qsatbef)
	199	qsatbef = qsatbef * zcor
	200	zqsatenv(ind1, ind2) = qsatbef
[2686]	201
	202
[5144]	203	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	204	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	205	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	206	aenv = 1. / (1. + (alenv * Lv / cppd))
	207	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	208
[5144]	209	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
[2686]	210
[5144]	211	sqrt2pi = sqrt(2. * pi)
	212	xenv = senv / (sqrt(2.) * sigma1s)
	213	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	214	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
[2686]	215
[5144]	216	IF (ctot(ind1, ind2)<1.e-3) THEN
	217	ctot(ind1, ind2) = 0.
	218	qcloud(ind1) = zqsatenv(ind1, ind2)
[2686]	219
[5144]	220	else
[2686]	221
[5144]	222	ctot(ind1, ind2) = ctot(ind1, ind2)
	223	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
[2686]	224
[5144]	225	endif
[2686]	226
[5144]	227	endif
	228	enddo
[2686]	229
[5144]	230	RETURN
	231	! end
	232	END SUBROUTINE cloudth
[2686]	233
	234
[5144]	235	!===========================================================================
	236	SUBROUTINE cloudth_vert(ngrid, klev, ind2, &
	237	ztv, po, zqta, fraca, &
	238	qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, &
	239	ratqs, zqs, t)
[5143]	240
[5144]	241	!===========================================================================
	242	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	243	! Date : 25 Mai 2010
	244	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	245	!===========================================================================
[2686]	246
[5144]	247	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, vert_alpha
	248	USE lmdz_yoethf
[5153]	249
[5144]	250	USE lmdz_yomcst
[2686]	251
[5144]	252	IMPLICIT NONE
[5153]	253	INCLUDE "FCTTRE.h"
[2686]	254
[5144]	255	INTEGER itap, ind1, ind2
	256	INTEGER ngrid, klev, klon, l, ig
[2686]	257
[5144]	258	REAL ztv(ngrid, klev)
	259	REAL po(ngrid)
	260	REAL zqenv(ngrid)
	261	REAL zqta(ngrid, klev)
[2686]	262
[5144]	263	REAL fraca(ngrid, klev + 1)
	264	REAL zpspsk(ngrid, klev)
	265	REAL paprs(ngrid, klev + 1)
	266	REAL pplay(ngrid, klev)
	267	REAL ztla(ngrid, klev)
	268	REAL zthl(ngrid, klev)
[2686]	269
[5144]	270	REAL zqsatth(ngrid, klev)
	271	REAL zqsatenv(ngrid, klev)
[2686]	272
[5144]	273	REAL sigma1(ngrid, klev)
	274	REAL sigma2(ngrid, klev)
	275	REAL qlth(ngrid, klev)
	276	REAL qlenv(ngrid, klev)
	277	REAL qltot(ngrid, klev)
	278	REAL cth(ngrid, klev)
	279	REAL cenv(ngrid, klev)
	280	REAL ctot(ngrid, klev)
	281	REAL rneb(ngrid, klev)
	282	REAL t(ngrid, klev)
	283	REAL qsatmmussig1, qsatmmussig2, sqrt2pi, pi
	284	REAL rdd, cppd, Lv, sqrt2, sqrtpi
	285	REAL alth, alenv, ath, aenv
	286	REAL sth, senv, sigma1s, sigma2s, xth, xenv
	287	REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv
	288	REAL IntJ, IntI1, IntI2, IntI3, coeffqlenv, coeffqlth
	289	REAL Tbef, zdelta, qsatbef, zcor
	290	REAL qlbef
	291	REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur
	292	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	293	! (J Jouhaud, JL Dufresne, JB Madeleine)
[2686]	294
[5144]	295	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	296	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	297	REAL zqs(ngrid), qcloud(ngrid)
	298	REAL erf
[2686]	299
[5144]	300	!------------------------------------------------------------------------------
	301	! Initialisation des variables r?elles
	302	!------------------------------------------------------------------------------
	303	sigma1(:, :) = 0.
	304	sigma2(:, :) = 0.
	305	qlth(:, :) = 0.
	306	qlenv(:, :) = 0.
	307	qltot(:, :) = 0.
	308	rneb(:, :) = 0.
	309	qcloud(:) = 0.
	310	cth(:, :) = 0.
	311	cenv(:, :) = 0.
	312	ctot(:, :) = 0.
	313	qsatmmussig1 = 0.
	314	qsatmmussig2 = 0.
	315	rdd = 287.04
	316	cppd = 1005.7
	317	pi = 3.14159
	318	Lv = 2.5e6
	319	sqrt2pi = sqrt(2. * pi)
	320	sqrt2 = sqrt(2.)
	321	sqrtpi = sqrt(pi)
[2686]	322
[5144]	323	!-------------------------------------------------------------------------------
	324	! Calcul de la fraction du thermique et des ?cart-types des distributions
	325	!-------------------------------------------------------------------------------
[5158]	326	DO ind1 = 1, ngrid
[2686]	327
[5144]	328	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
	329	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2))
[2686]	330
	331
[5144]	332	! zqenv(ind1)=po(ind1)
	333	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	334	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	335	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	336	qsatbef = MIN(0.5, qsatbef)
	337	zcor = 1. / (1. - retv * qsatbef)
	338	qsatbef = qsatbef * zcor
	339	zqsatenv(ind1, ind2) = qsatbef
[2686]	340
[5144]	341	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	342	aenv = 1. / (1. + (alenv * Lv / cppd))
	343	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	344
[5144]	345	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	346	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	347	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	348	qsatbef = MIN(0.5, qsatbef)
	349	zcor = 1. / (1. - retv * qsatbef)
	350	qsatbef = qsatbef * zcor
	351	zqsatth(ind1, ind2) = qsatbef
[2686]	352
[5144]	353	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2)
	354	ath = 1. / (1. + (alth * Lv / cppd))
	355	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2))
[2686]	356
	357
	358
[5144]	359	!------------------------------------------------------------------------------
	360	! Calcul des ?cart-types pour s
	361	!------------------------------------------------------------------------------
[2686]	362
[5144]	363	sigma1s = (0.92*0.5) (fraca(ind1, ind2)*0.5) / (1 - fraca(ind1, ind2)) ((sth - senv)2)0.5 + ratqs(ind1, ind2) * po(ind1)
	364	sigma2s = 0.09 * ((sth - senv)2)0.5 / (fraca(ind1, ind2) + 0.02)*0.5 + 0.002 zqta(ind1, ind2)
	365	! if (paprs(ind1,ind2).gt.90000) THEN
	366	! ratqs(ind1,ind2)=0.002
	367	! else
	368	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	369	! endif
	370	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	371	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	372	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	373	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[2686]	374
[5144]	375	!------------------------------------------------------------------------------
	376	! Calcul de l'eau condens?e et de la couverture nuageuse
	377	!------------------------------------------------------------------------------
	378	sqrt2pi = sqrt(2. * pi)
	379	xth = sth / (sqrt(2.) * sigma2s)
	380	xenv = senv / (sqrt(2.) * sigma1s)
	381	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	382	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	383	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	384
[5144]	385	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt(2.) * cth(ind1, ind2))
	386	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
	387	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	388
[5144]	389	IF (iflag_cloudth_vert == 1) THEN
	390	!-------------------------------------------------------------------------------
	391	! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
	392	!-------------------------------------------------------------------------------
	393	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	394	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	395	deltasenv = aenv * ratqs(ind1, ind2) * zqsatenv(ind1, ind2)
	396	deltasth = ath * ratqs(ind1, ind2) * zqsatth(ind1, ind2)
	397	! deltasenv=aenv0.01po(ind1)
	398	! deltasth=ath0.01zqta(ind1,ind2)
	399	xenv1 = (senv - deltasenv) / (sqrt(2.) * sigma1s)
	400	xenv2 = (senv + deltasenv) / (sqrt(2.) * sigma1s)
	401	xth1 = (sth - deltasth) / (sqrt(2.) * sigma2s)
	402	xth2 = (sth + deltasth) / (sqrt(2.) * sigma2s)
	403	coeffqlenv = (sigma1s)*2 / (2 sqrtpi * deltasenv)
	404	coeffqlth = (sigma2s)*2 / (2 sqrtpi * deltasth)
[2686]	405
[5144]	406	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth2))
	407	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv2))
	408	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	409
[5144]	410	IntJ = sigma1s * (exp(-1. * xenv1*2) / sqrt2pi) + 0.5 senv * (1 + erf(xenv1))
	411	IntI1 = coeffqlenv * 0.5 * (0.5 * sqrtpi * (erf(xenv2) - erf(xenv1)) + xenv1 * exp(-1. * xenv1*2) - xenv2 exp(-1. * xenv2**2))
	412	IntI2 = coeffqlenv * xenv2 * (exp(-1. * xenv2*2) - exp(-1. xenv1**2))
	413	IntI3 = coeffqlenv * 0.5 * sqrtpi * xenv2*2 (erf(xenv2) - erf(xenv1))
[2686]	414
[5144]	415	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	416	! qlenv(ind1,ind2)=IntJ
	417	! PRINT*, qlenv(ind1,ind2),'VERIF EAU'
[2686]	418
[5144]	419	IntJ = sigma2s * (exp(-1. * xth1*2) / sqrt2pi) + 0.5 sth * (1 + erf(xth1))
	420	! IntI1=coeffqlth((0.5xth1-xth2)exp(-1.xth1*2)+0.5xth2exp(-1.xth2**2))
	421	! IntI2=coeffqlth0.5sqrtpi(0.5+xth22)(erf(xth2)-erf(xth1))
	422	IntI1 = coeffqlth * 0.5 * (0.5 * sqrtpi * (erf(xth2) - erf(xth1)) + xth1 * exp(-1. * xth1*2) - xth2 exp(-1. * xth2**2))
	423	IntI2 = coeffqlth * xth2 * (exp(-1. * xth2*2) - exp(-1. xth1**2))
	424	IntI3 = coeffqlth * 0.5 * sqrtpi * xth2*2 (erf(xth2) - erf(xth1))
	425	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	426	! qlth(ind1,ind2)=IntJ
	427	! PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	428	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	429
[5144]	430	ELSE IF (iflag_cloudth_vert == 2) THEN
[2686]	431
[5144]	432	!-------------------------------------------------------------------------------
	433	! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
	434	!-------------------------------------------------------------------------------
	435	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	436	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	437	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	438	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	439	deltasenv = aenv * vert_alpha * sigma1s
	440	deltasth = ath * vert_alpha * sigma2s
[2686]	441
[5144]	442	xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s)
	443	xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s)
	444	xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s)
	445	xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s)
	446	! coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	447	! coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
[2686]	448
[5144]	449	cth(ind1, ind2) = 0.5 * (1. - 1. * erf(xth1))
	450	cenv(ind1, ind2) = 0.5 * (1. - 1. * erf(xenv1))
	451	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	452
[5144]	453	IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp(-1. * xenv2**2)
	454	IntI1 = (((senv + deltasenv)2 + (sigma1s)2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1))
	455	IntI2 = (sigma1s*2 / (4 deltasenv * sqrtpi)) * (xenv1 * exp(-1. * xenv1*2) - xenv2 exp(-1. * xenv2**2))
	456	IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp(-1. * xenv1*2) - exp(-1. xenv2**2))
[2686]	457
[5144]	458	! IntI1=0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv12)-xenv2exp(-1.xenv2*2))
	459	! IntI2=xenv2(exp(-1.xenv2*2)-exp(-1.xenv1**2))
	460	! IntI3=0.5sqrtpixenv2*2(erf(xenv2)-erf(xenv1))
[2686]	461
[5144]	462	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	463	! qlenv(ind1,ind2)=IntJ
	464	! PRINT*, qlenv(ind1,ind2),'VERIF EAU'
[2686]	465
[5144]	466	IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp(-1. * xth2**2)
	467	IntI1 = (((sth + deltasth)2 + (sigma2s)2) / (8 * deltasth)) * (erf(xth2) - erf(xth1))
	468	IntI2 = (sigma2s*2 / (4 deltasth * sqrtpi)) * (xth1 * exp(-1. * xth1*2) - xth2 exp(-1. * xth2**2))
	469	IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp(-1. * xth1*2) - exp(-1. xth2**2))
[2686]	470
[5144]	471	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	472	! qlth(ind1,ind2)=IntJ
	473	! PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	474	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
	475
	476	ENDIF ! of if (iflag_cloudth_vert==1 or 2)
	477
	478	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	479
	480	IF (cenv(ind1, ind2)<1.e-10.OR.cth(ind1, ind2)<1.e-10) THEN
	481	ctot(ind1, ind2) = 0.
	482	qcloud(ind1) = zqsatenv(ind1, ind2)
	483
	484	else
	485
	486	ctot(ind1, ind2) = ctot(ind1, ind2)
	487	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
	488	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	489	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	490
	491	endif
	492
	493
	494
	495	! PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	496
[2686]	497	else ! gaussienne environnement seule
	498
[5144]	499	zqenv(ind1) = po(ind1)
	500	Tbef = t(ind1, ind2)
	501	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	502	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	503	qsatbef = MIN(0.5, qsatbef)
	504	zcor = 1. / (1. - retv * qsatbef)
	505	qsatbef = qsatbef * zcor
	506	zqsatenv(ind1, ind2) = qsatbef
[2686]	507
	508
[5144]	509	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	510	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	511	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	512	aenv = 1. / (1. + (alenv * Lv / cppd))
	513	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	514
[5144]	515	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
[2686]	516
[5144]	517	sqrt2pi = sqrt(2. * pi)
	518	xenv = senv / (sqrt(2.) * sigma1s)
	519	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	520	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
[2686]	521
[5144]	522	IF (ctot(ind1, ind2)<1.e-3) THEN
	523	ctot(ind1, ind2) = 0.
	524	qcloud(ind1) = zqsatenv(ind1, ind2)
[3493]	525
[5144]	526	else
[3493]	527
[5144]	528	ctot(ind1, ind2) = ctot(ind1, ind2)
	529	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
[3493]	530
[5144]	531	endif
[2686]	532
[5144]	533	endif
	534	enddo
[2686]	535
[5144]	536	RETURN
	537	! end
	538	END SUBROUTINE cloudth_vert
[2686]	539
	540
[5144]	541	SUBROUTINE cloudth_v3(ngrid, klev, ind2, &
	542	ztv, po, zqta, fraca, &
	543	qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
[5226]	544	ratqs, sigma_qtherm,zqs, t, &
[5144]	545	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	546
[5144]	547	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert
	548	USE lmdz_yoethf
[5153]	549
[5144]	550	USE lmdz_yomcst
[2686]	551
[5144]	552	IMPLICIT NONE
[5153]	553	INCLUDE "FCTTRE.h"
[2686]	554
[4674]	555
[5144]	556	!===========================================================================
	557	! Author : Arnaud Octavio Jam (LMD/CNRS)
	558	! Date : 25 Mai 2010
	559	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	560	!===========================================================================
[4674]	561
[5144]	562	INTEGER, INTENT(IN) :: ind2
	563	INTEGER, INTENT(IN) :: ngrid, klev
[2686]	564
[5144]	565	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ztv
	566	REAL, DIMENSION(ngrid), INTENT(IN) :: po
	567	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: zqta
	568	REAL, DIMENSION(ngrid, klev + 1), INTENT(IN) :: fraca
	569	REAL, DIMENSION(ngrid), INTENT(OUT) :: qcloud
	570	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: ctot
	571	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: ctot_vol
	572	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: zpspsk
	573	REAL, DIMENSION(ngrid, klev + 1), INTENT(IN) :: paprs
	574	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: pplay
	575	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ztla
	576	REAL, DIMENSION(ngrid, klev), INTENT(INOUT) :: zthl
[5226]	577	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ratqs, sigma_qtherm
[5144]	578	REAL, DIMENSION(ngrid), INTENT(IN) :: zqs
	579	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: t
	580	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
[2686]	581
[5144]	582	REAL zqenv(ngrid)
	583	REAL zqsatth(ngrid, klev)
	584	REAL zqsatenv(ngrid, klev)
[4674]	585
[5144]	586	REAL sigma1(ngrid, klev)
	587	REAL sigma2(ngrid, klev)
	588	REAL qlth(ngrid, klev)
	589	REAL qlenv(ngrid, klev)
	590	REAL qltot(ngrid, klev)
	591	REAL cth(ngrid, klev)
	592	REAL cenv(ngrid, klev)
	593	REAL cth_vol(ngrid, klev)
	594	REAL cenv_vol(ngrid, klev)
	595	REAL rneb(ngrid, klev)
	596	REAL qsatmmussig1, qsatmmussig2, sqrt2pi, sqrt2, sqrtpi, pi
	597	REAL rdd, cppd, Lv
	598	REAL alth, alenv, ath, aenv
	599	REAL sth, senv, sigma1s, sigma2s, xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2
	600	REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks
	601	REAL Tbef, zdelta, qsatbef, zcor
	602	REAL qlbef
	603	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	604	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	605	REAL erf
	606
	607	INTEGER :: ind1, l, ig
	608
	609	IF (iflag_cloudth_vert>=1) THEN
	610	CALL cloudth_vert_v3(ngrid, klev, ind2, &
	611	ztv, po, zqta, fraca, &
	612	qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
[5226]	613	ratqs, sigma_qtherm, zqs, t, &
[5144]	614	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	615	RETURN
[5144]	616	ENDIF
	617	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2686]	618
	619
[5144]	620	!-------------------------------------------------------------------------------
	621	! Initialisation des variables r?elles
	622	!-------------------------------------------------------------------------------
	623	sigma1(:, :) = 0.
	624	sigma2(:, :) = 0.
	625	qlth(:, :) = 0.
	626	qlenv(:, :) = 0.
	627	qltot(:, :) = 0.
	628	rneb(:, :) = 0.
	629	qcloud(:) = 0.
	630	cth(:, :) = 0.
	631	cenv(:, :) = 0.
	632	ctot(:, :) = 0.
	633	cth_vol(:, :) = 0.
	634	cenv_vol(:, :) = 0.
	635	ctot_vol(:, :) = 0.
	636	qsatmmussig1 = 0.
	637	qsatmmussig2 = 0.
	638	rdd = 287.04
	639	cppd = 1005.7
	640	pi = 3.14159
	641	Lv = 2.5e6
	642	sqrt2pi = sqrt(2. * pi)
	643	sqrt2 = sqrt(2.)
	644	sqrtpi = sqrt(pi)
[2686]	645
	646
[5144]	647	!-------------------------------------------------------------------------------
	648	! Cloud fraction in the thermals and standard deviation of the PDFs
	649	!-------------------------------------------------------------------------------
[5158]	650	DO ind1 = 1, ngrid
[2686]	651
[5144]	652	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
	653	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2))
[2686]	654
[5144]	655	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	656	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	657	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	658	qsatbef = MIN(0.5, qsatbef)
	659	zcor = 1. / (1. - retv * qsatbef)
	660	qsatbef = qsatbef * zcor
	661	zqsatenv(ind1, ind2) = qsatbef
[2686]	662
[5144]	663	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84
	664	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84
	665	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84
[2686]	666
[5144]	667	!po = qt de l'environnement ET des thermique
	668	!zqenv = qt environnement
	669	!zqsatenv = qsat environnement
	670	!zthl = Tl environnement
[2686]	671
[5144]	672	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	673	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	674	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	675	qsatbef = MIN(0.5, qsatbef)
	676	zcor = 1. / (1. - retv * qsatbef)
	677	qsatbef = qsatbef * zcor
	678	zqsatth(ind1, ind2) = qsatbef
[2686]	679
[5144]	680	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84
	681	ath = 1. / (1. + (alth * Lv / cppd)) !al, p84
	682	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84
[2686]	683
[5144]	684	!zqta = qt thermals
	685	!zqsatth = qsat thermals
	686	!ztla = Tl thermals
[2686]	687
[5144]	688	!------------------------------------------------------------------------------
	689	! s standard deviations
	690	!------------------------------------------------------------------------------
[2686]	691
[5144]	692	! tests
	693	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	694	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+ratqs(ind1,ind2)*po(ind1)
	695	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
	696	! final option
	697	sigma1s = (1.1*0.5) (fraca(ind1, ind2)*0.6) / (1 - fraca(ind1, ind2)) ((sth - senv)2)0.5 + ratqs(ind1, ind2) * po(ind1)
	698	sigma2s = 0.11 * ((sth - senv)2)0.5 / (fraca(ind1, ind2) + 0.01)*0.4 + 0.002 zqta(ind1, ind2)
[2686]	699
[5144]	700	!------------------------------------------------------------------------------
	701	! Condensed water and cloud cover
	702	!------------------------------------------------------------------------------
	703	xth = sth / (sqrt2 * sigma2s)
	704	xenv = senv / (sqrt2 * sigma1s)
	705	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	706	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	707	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
	708	ctot_vol(ind1, ind2) = ctot(ind1, ind2)
[2686]	709
[5144]	710	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt2 * cth(ind1, ind2))
	711	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * cenv(ind1, ind2))
	712	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	713
[5144]	714	IF (ctot(ind1, ind2)<1.e-10) THEN
	715	ctot(ind1, ind2) = 0.
	716	qcloud(ind1) = zqsatenv(ind1, ind2)
	717	else
	718	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
	719	endif
	720
[2686]	721	else ! Environnement only, follow the if l.110
	722
[5144]	723	zqenv(ind1) = po(ind1)
	724	Tbef = t(ind1, ind2)
	725	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	726	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	727	qsatbef = MIN(0.5, qsatbef)
	728	zcor = 1. / (1. - retv * qsatbef)
	729	qsatbef = qsatbef * zcor
	730	zqsatenv(ind1, ind2) = qsatbef
[2686]	731
[5144]	732	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	733	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	734	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	735	aenv = 1. / (1. + (alenv * Lv / cppd))
	736	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	737
[5144]	738	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
[2686]	739
[5144]	740	xenv = senv / (sqrt2 * sigma1s)
	741	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	742	ctot_vol(ind1, ind2) = ctot(ind1, ind2)
	743	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * cenv(ind1, ind2))
[2686]	744
[5144]	745	IF (ctot(ind1, ind2)<1.e-3) THEN
	746	ctot(ind1, ind2) = 0.
	747	qcloud(ind1) = zqsatenv(ind1, ind2)
	748	else
	749	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
	750	endif
	751
[2686]	752	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
[5144]	753	enddo ! from the loop on ngrid l.108
	754	RETURN
	755	! end
	756	END SUBROUTINE cloudth_v3
[2686]	757
	758
[5144]	759	!===========================================================================
	760	SUBROUTINE cloudth_vert_v3(ngrid, klev, ind2, &
	761	ztv, po, zqta, fraca, &
	762	qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
[5226]	763	ratqs, sigma_qtherm, zqs, t, &
[5144]	764	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	765
[5144]	766	!===========================================================================
	767	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	768	! Date : 25 Mai 2010
	769	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	770	!===========================================================================
[2686]	771
[5144]	772	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs
	773	USE lmdz_cloudth_ini, ONLY: vert_alpha, vert_alpha_th, sigma1s_factor, sigma1s_power, sigma2s_factor, sigma2s_power, cloudth_ratqsmin, iflag_cloudth_vert_noratqs
	774	USE lmdz_yoethf
[5153]	775
[5144]	776	USE lmdz_yomcst
[2686]	777
[5144]	778	IMPLICIT NONE
[5153]	779	INCLUDE "FCTTRE.h"
[2686]	780
[5144]	781	INTEGER itap, ind1, ind2
	782	INTEGER ngrid, klev, klon, l, ig
	783	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
[2686]	784
[5144]	785	REAL ztv(ngrid, klev)
	786	REAL po(ngrid)
	787	REAL zqenv(ngrid)
	788	REAL zqta(ngrid, klev)
[4651]	789
[5144]	790	REAL fraca(ngrid, klev + 1)
	791	REAL zpspsk(ngrid, klev)
	792	REAL paprs(ngrid, klev + 1)
	793	REAL pplay(ngrid, klev)
	794	REAL ztla(ngrid, klev)
	795	REAL zthl(ngrid, klev)
[2686]	796
[5144]	797	REAL zqsatth(ngrid, klev)
	798	REAL zqsatenv(ngrid, klev)
[2957]	799
[5144]	800	REAL sigma1(ngrid, klev)
	801	REAL sigma2(ngrid, klev)
	802	REAL qlth(ngrid, klev)
	803	REAL qlenv(ngrid, klev)
	804	REAL qltot(ngrid, klev)
	805	REAL cth(ngrid, klev)
	806	REAL cenv(ngrid, klev)
	807	REAL ctot(ngrid, klev)
	808	REAL cth_vol(ngrid, klev)
	809	REAL cenv_vol(ngrid, klev)
	810	REAL ctot_vol(ngrid, klev)
	811	REAL rneb(ngrid, klev)
	812	REAL t(ngrid, klev)
	813	REAL qsatmmussig1, qsatmmussig2, sqrtpi, sqrt2, sqrt2pi, pi
	814	REAL rdd, cppd, Lv
	815	REAL alth, alenv, ath, aenv
	816	REAL sth, senv, sigma1s, sigma2s, sigma1s_fraca, sigma1s_ratqs
	817	REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks
	818	REAL xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2
	819	REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv
	820	REAL IntJ, IntI1, IntI2, IntI3, IntJ_CF, IntI1_CF, IntI3_CF, coeffqlenv, coeffqlth
	821	REAL Tbef, zdelta, qsatbef, zcor
	822	REAL qlbef
[5226]	823	REAL ratqs(ngrid, klev), sigma_qtherm(ngrid,klev) ! determine la largeur de distribution de vapeur
[5144]	824	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	825	! (J Jouhaud, JL Dufresne, JB Madeleine)
[2686]	826
[5144]	827	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	828	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	829	REAL zqs(ngrid), qcloud(ngrid)
	830	REAL erf
[3493]	831
[5144]	832	REAL rhodz(ngrid, klev)
	833	REAL zrho(ngrid, klev)
	834	REAL dz(ngrid, klev)
[3493]	835
[5144]	836	DO ind1 = 1, ngrid
	837	!Layer calculation
	838	rhodz(ind1, ind2) = (paprs(ind1, ind2) - paprs(ind1, ind2 + 1)) / rg !kg/m2
	839	zrho(ind1, ind2) = pplay(ind1, ind2) / t(ind1, ind2) / rd !kg/m3
	840	dz(ind1, ind2) = rhodz(ind1, ind2) / zrho(ind1, ind2) !m : epaisseur de la couche en metre
	841	END DO
[3999]	842
[5144]	843	!------------------------------------------------------------------------------
	844	! Initialize
	845	!------------------------------------------------------------------------------
[2686]	846
[5144]	847	sigma1(:, :) = 0.
	848	sigma2(:, :) = 0.
	849	qlth(:, :) = 0.
	850	qlenv(:, :) = 0.
	851	qltot(:, :) = 0.
	852	rneb(:, :) = 0.
	853	qcloud(:) = 0.
	854	cth(:, :) = 0.
	855	cenv(:, :) = 0.
	856	ctot(:, :) = 0.
	857	cth_vol(:, :) = 0.
	858	cenv_vol(:, :) = 0.
	859	ctot_vol(:, :) = 0.
	860	qsatmmussig1 = 0.
	861	qsatmmussig2 = 0.
	862	rdd = 287.04
	863	cppd = 1005.7
	864	pi = 3.14159
	865	Lv = 2.5e6
	866	sqrt2pi = sqrt(2. * pi)
	867	sqrt2 = sqrt(2.)
	868	sqrtpi = sqrt(pi)
[2957]	869
[2686]	870
	871
[5144]	872	!-------------------------------------------------------------------------------
	873	! Calcul de la fraction du thermique et des ecart-types des distributions
	874	!-------------------------------------------------------------------------------
[5158]	875	DO ind1 = 1, ngrid
[2686]	876
[5144]	877	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) then !Thermal and environnement
[2686]	878
[5144]	879	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2)) !qt = aqtth + (1-a)qtenv
[2686]	880
[5144]	881	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	882	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	883	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	884	qsatbef = MIN(0.5, qsatbef)
	885	zcor = 1. / (1. - retv * qsatbef)
	886	qsatbef = qsatbef * zcor
	887	zqsatenv(ind1, ind2) = qsatbef
[2686]	888
[5144]	889	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84
	890	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84
	891	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84
[2686]	892
[5144]	893	!zqenv = qt environnement
	894	!zqsatenv = qsat environnement
	895	!zthl = Tl environnement
[2686]	896
[5144]	897	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	898	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	899	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	900	qsatbef = MIN(0.5, qsatbef)
	901	zcor = 1. / (1. - retv * qsatbef)
	902	qsatbef = qsatbef * zcor
	903	zqsatth(ind1, ind2) = qsatbef
[2686]	904
[5144]	905	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84
	906	ath = 1. / (1. + (alth * Lv / cppd)) !al, p84
	907	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84
[2686]	908
	909
[5144]	910	!zqta = qt thermals
	911	!zqsatth = qsat thermals
	912	!ztla = Tl thermals
	913	!------------------------------------------------------------------------------
	914	! s standard deviation
	915	!------------------------------------------------------------------------------
[2686]	916
[5144]	917	sigma1s_fraca = (sigma1s_factor*0.5) (fraca(ind1, ind2)**sigma1s_power) / &
	918	(1 - fraca(ind1, ind2)) * ((sth - senv)2)0.5
	919	! sigma1s_fraca = (1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5
	920	IF (cloudth_ratqsmin>0.) THEN
	921	sigma1s_ratqs = cloudth_ratqsmin * po(ind1)
	922	ELSE
	923	sigma1s_ratqs = ratqs(ind1, ind2) * po(ind1)
	924	ENDIF
	925	sigma1s = sigma1s_fraca + sigma1s_ratqs
[5226]	926	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
	927	IF (iflag_ratqs==10.or.iflag_ratqs==11) THEN
[5144]	928	sigma1s = ratqs(ind1, ind2) * po(ind1) * aenv
[5226]	929	IF (iflag_ratqs==10.and.sigma_qtherm(ind1, ind2) /= 0) then
	930	sigma2s = sigma_qtherm(ind1, ind2)*ath
	931	ENDIF
	932	ENDIF
[5144]	933	! tests
	934	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	935	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
	936	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
	937	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
	938	! if (paprs(ind1,ind2).gt.90000) THEN
	939	! ratqs(ind1,ind2)=0.002
	940	! else
	941	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	942	! endif
	943	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	944	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	945	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	946	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[2686]	947
[5144]	948	IF (iflag_cloudth_vert == 1) THEN
	949	!-------------------------------------------------------------------------------
	950	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
	951	!-------------------------------------------------------------------------------
[2686]	952
[5144]	953	deltasenv = aenv * ratqs(ind1, ind2) * zqsatenv(ind1, ind2)
	954	deltasth = ath * ratqs(ind1, ind2) * zqsatth(ind1, ind2)
[2686]	955
[5144]	956	xenv1 = (senv - deltasenv) / (sqrt(2.) * sigma1s)
	957	xenv2 = (senv + deltasenv) / (sqrt(2.) * sigma1s)
	958	xth1 = (sth - deltasth) / (sqrt(2.) * sigma2s)
	959	xth2 = (sth + deltasth) / (sqrt(2.) * sigma2s)
	960	coeffqlenv = (sigma1s)*2 / (2 sqrtpi * deltasenv)
	961	coeffqlth = (sigma2s)*2 / (2 sqrtpi * deltasth)
[2686]	962
[5144]	963	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth2))
	964	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv2))
	965	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	966
[5144]	967	! Environment
	968	IntJ = sigma1s * (exp(-1. * xenv1*2) / sqrt2pi) + 0.5 senv * (1 + erf(xenv1))
	969	IntI1 = coeffqlenv * 0.5 * (0.5 * sqrtpi * (erf(xenv2) - erf(xenv1)) + xenv1 * exp(-1. * xenv1*2) - xenv2 exp(-1. * xenv2**2))
	970	IntI2 = coeffqlenv * xenv2 * (exp(-1. * xenv2*2) - exp(-1. xenv1**2))
	971	IntI3 = coeffqlenv * 0.5 * sqrtpi * xenv2*2 (erf(xenv2) - erf(xenv1))
[2686]	972
[5144]	973	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
[2945]	974
[5144]	975	! Thermal
	976	IntJ = sigma2s * (exp(-1. * xth1*2) / sqrt2pi) + 0.5 sth * (1 + erf(xth1))
	977	IntI1 = coeffqlth * 0.5 * (0.5 * sqrtpi * (erf(xth2) - erf(xth1)) + xth1 * exp(-1. * xth1*2) - xth2 exp(-1. * xth2**2))
	978	IntI2 = coeffqlth * xth2 * (exp(-1. * xth2*2) - exp(-1. xth1**2))
	979	IntI3 = coeffqlth * 0.5 * sqrtpi * xth2*2 (erf(xth2) - erf(xth1))
	980	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	981	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	982
[5144]	983	ELSE IF (iflag_cloudth_vert >= 3) THEN
	984	IF (iflag_cloudth_vert < 5) THEN
	985	!-------------------------------------------------------------------------------
	986	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
	987	!-------------------------------------------------------------------------------
	988	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	989	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	990	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	991	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	992	IF (iflag_cloudth_vert == 3) THEN
	993	deltasenv = aenv * vert_alpha * sigma1s
	994	deltasth = ath * vert_alpha_th * sigma2s
	995	ELSE IF (iflag_cloudth_vert == 4) THEN
	996	IF (iflag_cloudth_vert_noratqs == 1) THEN
	997	deltasenv = vert_alpha * max(sigma1s_fraca, 1e-10)
	998	deltasth = vert_alpha_th * sigma2s
	999	ELSE
	1000	deltasenv = vert_alpha * sigma1s
	1001	deltasth = vert_alpha_th * sigma2s
	1002	ENDIF
	1003	ENDIF
[2686]	1004
[5144]	1005	xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s)
	1006	xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s)
	1007	exp_xenv1 = exp(-1. * xenv1**2)
	1008	exp_xenv2 = exp(-1. * xenv2**2)
	1009	xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s)
	1010	xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s)
	1011	exp_xth1 = exp(-1. * xth1**2)
	1012	exp_xth2 = exp(-1. * xth2**2)
[2686]	1013
[5144]	1014	!CF_surfacique
	1015	cth(ind1, ind2) = 0.5 * (1. - 1. * erf(xth1))
	1016	cenv(ind1, ind2) = 0.5 * (1. - 1. * erf(xenv1))
	1017	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[3493]	1018
	1019
[5144]	1020	!CF_volumique & eau condense
	1021	!environnement
	1022	IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp_xenv2
	1023	IntJ_CF = 0.5 * (1. - 1. * erf(xenv2))
	1024	IF (deltasenv < 1.e-10) THEN
	1025	qlenv(ind1, ind2) = IntJ
	1026	cenv_vol(ind1, ind2) = IntJ_CF
	1027	else
	1028	IntI1 = (((senv + deltasenv)2 + (sigma1s)2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1))
	1029	IntI2 = (sigma1s*2 / (4 deltasenv * sqrtpi)) * (xenv1 * exp_xenv1 - xenv2 * exp_xenv2)
	1030	IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp_xenv1 - exp_xenv2)
	1031	IntI1_CF = ((senv + deltasenv) * (erf(xenv2) - erf(xenv1))) / (4 * deltasenv)
	1032	IntI3_CF = (sqrt2 * sigma1s * (exp_xenv1 - exp_xenv2)) / (4 * sqrtpi * deltasenv)
	1033	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	1034	cenv_vol(ind1, ind2) = IntJ_CF + IntI1_CF + IntI3_CF
	1035	endif
[3493]	1036
[5144]	1037	!thermique
	1038	IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp_xth2
	1039	IntJ_CF = 0.5 * (1. - 1. * erf(xth2))
	1040	IF (deltasth < 1.e-10) THEN
	1041	qlth(ind1, ind2) = IntJ
	1042	cth_vol(ind1, ind2) = IntJ_CF
	1043	else
	1044	IntI1 = (((sth + deltasth)2 + (sigma2s)2) / (8 * deltasth)) * (erf(xth2) - erf(xth1))
	1045	IntI2 = (sigma2s*2 / (4 deltasth * sqrtpi)) * (xth1 * exp_xth1 - xth2 * exp_xth2)
	1046	IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp_xth1 - exp_xth2)
	1047	IntI1_CF = ((sth + deltasth) * (erf(xth2) - erf(xth1))) / (4 * deltasth)
	1048	IntI3_CF = (sqrt2 * sigma2s * (exp_xth1 - exp_xth2)) / (4 * sqrtpi * deltasth)
	1049	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	1050	cth_vol(ind1, ind2) = IntJ_CF + IntI1_CF + IntI3_CF
	1051	endif
[3493]	1052
[5144]	1053	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
	1054	ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2)
[3493]	1055
[5144]	1056	ELSE IF (iflag_cloudth_vert == 5) THEN
	1057	sigma1s = (0.71794 + 0.000498239 * dz(ind1, ind2)) * (fraca(ind1, ind2)**0.5) &
	1058	/ (1 - fraca(ind1, ind2)) * (((sth - senv)2)0.5) &
	1059	+ ratqs(ind1, ind2) * po(ind1) !Environment
	1060	sigma2s = (0.03218 + 0.000092655 * dz(ind1, ind2)) / ((fraca(ind1, ind2) + 0.02)*0.5) (((sth - senv)2)0.5) + 0.002 * zqta(ind1, ind2) !Thermals
	1061	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	1062	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	1063	xth = sth / (sqrt(2.) * sigma2s)
	1064	xenv = senv / (sqrt(2.) * sigma1s)
[3493]	1065
[5144]	1066	!Volumique
	1067	cth_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	1068	cenv_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	1069	ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2)
[5160]	1070	!PRINT *,'jeanjean_CV=',ctot_vol(ind1,ind2)
[3493]	1071
[5144]	1072	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt(2.) * cth_vol(ind1, ind2))
	1073	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv_vol(ind1, ind2))
	1074	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[3493]	1075
[5144]	1076	!Surfacique
	1077	!Neggers
	1078	!beta=0.0044
	1079	!inverse_rho=1.+beta*dz(ind1,ind2)
[5160]	1080	!PRINT *,'jeanjean : beta=',beta
[5144]	1081	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
	1082	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
	1083	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
[3493]	1084
[5144]	1085	!Brooks
	1086	a_Brooks = 0.6694
	1087	b_Brooks = 0.1882
	1088	A_Maj_Brooks = 0.1635 !-- sans shear
	1089	!A_Maj_Brooks=0.17 !-- ARM LES
	1090	!A_Maj_Brooks=0.18 !-- RICO LES
	1091	!A_Maj_Brooks=0.19 !-- BOMEX LES
	1092	Dx_Brooks = 200000.
	1093	f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)*(a_Brooks)) (Dx_Brooks**(-b_Brooks))
[5160]	1094	!PRINT *,'jeanjean_f=',f_Brooks
[3493]	1095
[5144]	1096	cth(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(cth_vol(ind1, ind2), 1.))) - 1.))
	1097	cenv(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(cenv_vol(ind1, ind2), 1.))) - 1.))
	1098	ctot(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.))
[5160]	1099	!PRINT *,'JJ_ctot_1',ctot(ind1,ind2)
[2686]	1100
[5144]	1101	ENDIF ! of if (iflag_cloudth_vert<5)
	1102	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
[2686]	1103
[5144]	1104	! if (ctot(ind1,ind2).lt.1.e-10) THEN
	1105	IF (cenv(ind1, ind2)<1.e-10.OR.cth(ind1, ind2)<1.e-10) THEN
	1106	ctot(ind1, ind2) = 0.
	1107	ctot_vol(ind1, ind2) = 0.
	1108	qcloud(ind1) = zqsatenv(ind1, ind2)
[2686]	1109
[5144]	1110	else
[2686]	1111
[5144]	1112	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
	1113	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	1114	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	1115
	1116	endif
	1117
[3493]	1118	else ! gaussienne environnement seule
[3999]	1119
[5144]	1120	zqenv(ind1) = po(ind1)
	1121	Tbef = t(ind1, ind2)
	1122	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1123	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1124	qsatbef = MIN(0.5, qsatbef)
	1125	zcor = 1. / (1. - retv * qsatbef)
	1126	qsatbef = qsatbef * zcor
	1127	zqsatenv(ind1, ind2) = qsatbef
[2686]	1128
	1129
[5144]	1130	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	1131	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	1132	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	1133	aenv = 1. / (1. + (alenv * Lv / cppd))
	1134	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
	1135	sth = 0.
[2686]	1136
[5144]	1137	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
	1138	sigma2s = 0.
[2686]	1139
[5144]	1140	sqrt2pi = sqrt(2. * pi)
	1141	xenv = senv / (sqrt(2.) * sigma1s)
	1142	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	1143	ctot_vol(ind1, ind2) = ctot(ind1, ind2)
	1144	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
[2686]	1145
[5144]	1146	IF (ctot(ind1, ind2)<1.e-3) THEN
	1147	ctot(ind1, ind2) = 0.
	1148	qcloud(ind1) = zqsatenv(ind1, ind2)
[3493]	1149
[5144]	1150	else
[3493]	1151
[5144]	1152	! ctot(ind1,ind2)=ctot(ind1,ind2)
	1153	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
[3493]	1154
[5144]	1155	endif
	1156
[2686]	1157	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
[2958]	1158	! Outputs used to check the PDFs
[5144]	1159	cloudth_senv(ind1, ind2) = senv
	1160	cloudth_sth(ind1, ind2) = sth
	1161	cloudth_sigmaenv(ind1, ind2) = sigma1s
	1162	cloudth_sigmath(ind1, ind2) = sigma2s
[2958]	1163
[5144]	1164	enddo ! from the loop on ngrid l.333
	1165	RETURN
	1166	! end
	1167	END SUBROUTINE cloudth_vert_v3
[3493]	1168
[5144]	1169	SUBROUTINE cloudth_v6(ngrid, klev, ind2, &
	1170	ztv, po, zqta, fraca, &
	1171	qcloud, ctot_surf, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
	1172	ratqs, zqs, T, &
	1173	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[3493]	1174
[5144]	1175	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert
	1176	USE lmdz_yoethf
[5153]	1177
[5144]	1178	USE lmdz_yomcst
[3493]	1179
[5144]	1180	IMPLICIT NONE
[5153]	1181	INCLUDE "FCTTRE.h"
[3493]	1182
[5144]	1183	!Domain variables
	1184	INTEGER ngrid !indice Max lat-lon
	1185	INTEGER klev !indice Max alt
	1186	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
	1187	INTEGER ind1 !indice in [1:ngrid]
	1188	INTEGER ind2 !indice in [1:klev]
	1189	!thermal plume fraction
	1190	REAL fraca(ngrid, klev + 1) !thermal plumes fraction in the gridbox
	1191	!temperatures
	1192	REAL T(ngrid, klev) !temperature
	1193	REAL zpspsk(ngrid, klev) !factor (p/p0)**kappa (used for potential variables)
	1194	REAL ztv(ngrid, klev) !potential temperature (voir thermcell_env.F90)
	1195	REAL ztla(ngrid, klev) !liquid temperature in the thermals (Tl_th)
	1196	REAL zthl(ngrid, klev) !liquid temperature in the environment (Tl_env)
	1197	!pressure
	1198	REAL paprs(ngrid, klev + 1) !pressure at the interface of levels
	1199	REAL pplay(ngrid, klev) !pressure at the middle of the level
	1200	!humidity
	1201	REAL ratqs(ngrid, klev) !width of the total water subgrid-scale distribution
	1202	REAL po(ngrid) !total water (qt)
	1203	REAL zqenv(ngrid) !total water in the environment (qt_env)
	1204	REAL zqta(ngrid, klev) !total water in the thermals (qt_th)
	1205	REAL zqsatth(ngrid, klev) !water saturation level in the thermals (q_sat_th)
	1206	REAL zqsatenv(ngrid, klev) !water saturation level in the environment (q_sat_env)
	1207	REAL qlth(ngrid, klev) !condensed water in the thermals
	1208	REAL qlenv(ngrid, klev) !condensed water in the environment
	1209	REAL qltot(ngrid, klev) !condensed water in the gridbox
	1210	!cloud fractions
	1211	REAL cth_vol(ngrid, klev) !cloud fraction by volume in the thermals
	1212	REAL cenv_vol(ngrid, klev) !cloud fraction by volume in the environment
	1213	REAL ctot_vol(ngrid, klev) !cloud fraction by volume in the gridbox
	1214	REAL cth_surf(ngrid, klev) !cloud fraction by surface in the thermals
	1215	REAL cenv_surf(ngrid, klev) !cloud fraction by surface in the environment
	1216	REAL ctot_surf(ngrid, klev) !cloud fraction by surface in the gridbox
	1217	!PDF of saturation deficit variables
	1218	REAL rdd, cppd, Lv
	1219	REAL Tbef, zdelta, qsatbef, zcor
	1220	REAL alth, alenv, ath, aenv
	1221	REAL sth, senv !saturation deficits in the thermals and environment
	1222	REAL sigma_env, sigma_th !standard deviations of the biGaussian PDF
	1223	!cloud fraction variables
	1224	REAL xth, xenv
	1225	REAL inverse_rho, beta !Neggers et al. (2011) method
	1226	REAL a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks !Brooks et al. (2005) method
	1227	!Incloud total water variables
	1228	REAL zqs(ngrid) !q_sat
	1229	REAL qcloud(ngrid) !eau totale dans le nuage
	1230	!Some arithmetic variables
	1231	REAL erf, pi, sqrt2, sqrt2pi
	1232	!Depth of the layer
	1233	REAL dz(ngrid, klev) !epaisseur de la couche en metre
	1234	REAL rhodz(ngrid, klev)
	1235	REAL zrho(ngrid, klev)
	1236	DO ind1 = 1, ngrid
	1237	rhodz(ind1, ind2) = (paprs(ind1, ind2) - paprs(ind1, ind2 + 1)) / rg ![kg/m2]
	1238	zrho(ind1, ind2) = pplay(ind1, ind2) / T(ind1, ind2) / rd ![kg/m3]
	1239	dz(ind1, ind2) = rhodz(ind1, ind2) / zrho(ind1, ind2) ![m]
	1240	END DO
[3493]	1241
[5144]	1242	!------------------------------------------------------------------------------
	1243	! Initialization
	1244	!------------------------------------------------------------------------------
	1245	qlth(:, :) = 0.
	1246	qlenv(:, :) = 0.
	1247	qltot(:, :) = 0.
	1248	cth_vol(:, :) = 0.
	1249	cenv_vol(:, :) = 0.
	1250	ctot_vol(:, :) = 0.
	1251	cth_surf(:, :) = 0.
	1252	cenv_surf(:, :) = 0.
	1253	ctot_surf(:, :) = 0.
	1254	qcloud(:) = 0.
	1255	rdd = 287.04
	1256	cppd = 1005.7
	1257	pi = 3.14159
	1258	Lv = 2.5e6
	1259	sqrt2 = sqrt(2.)
	1260	sqrt2pi = sqrt(2. * pi)
[3493]	1261
[5144]	1262	DO ind1 = 1, ngrid
	1263	!-------------------------------------------------------------------------------
	1264	!Both thermal and environment in the gridbox
	1265	!-------------------------------------------------------------------------------
	1266	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
[3493]	1267	!--------------------------------------------
	1268	!calcul de qsat_env
	1269	!--------------------------------------------
[5144]	1270	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	1271	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1272	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1273	qsatbef = MIN(0.5, qsatbef)
	1274	zcor = 1. / (1. - retv * qsatbef)
	1275	qsatbef = qsatbef * zcor
	1276	zqsatenv(ind1, ind2) = qsatbef
[3493]	1277	!--------------------------------------------
	1278	!calcul de s_env
	1279	!--------------------------------------------
[5144]	1280	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 these Arnaud Jam
	1281	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 these Arnaud Jam
	1282	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 these Arnaud Jam
[3493]	1283	!--------------------------------------------
	1284	!calcul de qsat_th
	1285	!--------------------------------------------
[5144]	1286	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	1287	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1288	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1289	qsatbef = MIN(0.5, qsatbef)
	1290	zcor = 1. / (1. - retv * qsatbef)
	1291	qsatbef = qsatbef * zcor
	1292	zqsatth(ind1, ind2) = qsatbef
[3493]	1293	!--------------------------------------------
[5144]	1294	!calcul de s_th
[3493]	1295	!--------------------------------------------
[5144]	1296	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84 these Arnaud Jam
	1297	ath = 1. / (1. + (alth * Lv / cppd)) !al, p84 these Arnaud Jam
	1298	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84 these Arnaud Jam
[3493]	1299	!--------------------------------------------
	1300	!calcul standard deviations bi-Gaussian PDF
	1301	!--------------------------------------------
[5144]	1302	sigma_th = (0.03218 + 0.000092655 * dz(ind1, ind2)) / ((fraca(ind1, ind2) + 0.01)*0.5) (((sth - senv)2)0.5) + 0.002 * zqta(ind1, ind2)
	1303	sigma_env = (0.71794 + 0.000498239 * dz(ind1, ind2)) * (fraca(ind1, ind2)**0.5) &
	1304	/ (1 - fraca(ind1, ind2)) * (((sth - senv)2)0.5) &
	1305	+ ratqs(ind1, ind2) * po(ind1)
	1306	xth = sth / (sqrt2 * sigma_th)
	1307	xenv = senv / (sqrt2 * sigma_env)
[3493]	1308	!--------------------------------------------
	1309	!Cloud fraction by volume CF_vol
	1310	!--------------------------------------------
[5144]	1311	cth_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	1312	cenv_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	1313	ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2)
	1314	!--------------------------------------------
[3493]	1315	!Condensed water qc
	1316	!--------------------------------------------
[5144]	1317	qlth(ind1, ind2) = sigma_th * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt2 * cth_vol(ind1, ind2))
	1318	qlenv(ind1, ind2) = sigma_env * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * cenv_vol(ind1, ind2))
	1319	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[3493]	1320	!--------------------------------------------
	1321	!Cloud fraction by surface CF_surf
	1322	!--------------------------------------------
	1323	!Method Neggers et al. (2011) : ok for cumulus clouds only
[5144]	1324	!beta=0.0044 (Jouhaud et al.2018)
	1325	!inverse_rho=1.+beta*dz(ind1,ind2)
	1326	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
[3493]	1327	!Method Brooks et al. (2005) : ok for all types of clouds
[5144]	1328	a_Brooks = 0.6694
	1329	b_Brooks = 0.1882
	1330	A_Maj_Brooks = 0.1635 !-- sans dependence au cisaillement de vent
	1331	Dx_Brooks = 200000. !-- si l'on considere des mailles de 200km de cote
	1332	f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)*(a_Brooks)) (Dx_Brooks**(-b_Brooks))
	1333	ctot_surf(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.))
[3493]	1334	!--------------------------------------------
	1335	!Incloud Condensed water qcloud
	1336	!--------------------------------------------
[5144]	1337	IF (ctot_surf(ind1, ind2) < 1.e-10) THEN
	1338	ctot_vol(ind1, ind2) = 0.
	1339	ctot_surf(ind1, ind2) = 0.
	1340	qcloud(ind1) = zqsatenv(ind1, ind2)
	1341	else
	1342	qcloud(ind1) = qltot(ind1, ind2) / ctot_vol(ind1, ind2) + zqs(ind1)
	1343	endif
[3493]	1344
	1345
	1346
[5144]	1347	!-------------------------------------------------------------------------------
	1348	!Environment only in the gridbox
	1349	!-------------------------------------------------------------------------------
[3493]	1350	ELSE
	1351	!--------------------------------------------
	1352	!calcul de qsat_env
	1353	!--------------------------------------------
[5144]	1354	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	1355	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1356	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1357	qsatbef = MIN(0.5, qsatbef)
	1358	zcor = 1. / (1. - retv * qsatbef)
	1359	qsatbef = qsatbef * zcor
	1360	zqsatenv(ind1, ind2) = qsatbef
[3493]	1361	!--------------------------------------------
	1362	!calcul de s_env
	1363	!--------------------------------------------
[5144]	1364	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 these Arnaud Jam
	1365	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 these Arnaud Jam
	1366	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 these Arnaud Jam
[3493]	1367	!--------------------------------------------
	1368	!calcul standard deviations Gaussian PDF
	1369	!--------------------------------------------
[5144]	1370	zqenv(ind1) = po(ind1)
	1371	sigma_env = ratqs(ind1, ind2) * zqenv(ind1)
	1372	xenv = senv / (sqrt2 * sigma_env)
[3493]	1373	!--------------------------------------------
	1374	!Cloud fraction by volume CF_vol
	1375	!--------------------------------------------
[5144]	1376	ctot_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
[3493]	1377	!--------------------------------------------
	1378	!Condensed water qc
	1379	!--------------------------------------------
[5144]	1380	qltot(ind1, ind2) = sigma_env * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * ctot_vol(ind1, ind2))
[3493]	1381	!--------------------------------------------
	1382	!Cloud fraction by surface CF_surf
	1383	!--------------------------------------------
	1384	!Method Neggers et al. (2011) : ok for cumulus clouds only
[5144]	1385	!beta=0.0044 (Jouhaud et al.2018)
	1386	!inverse_rho=1.+beta*dz(ind1,ind2)
	1387	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
[3493]	1388	!Method Brooks et al. (2005) : ok for all types of clouds
[5144]	1389	a_Brooks = 0.6694
	1390	b_Brooks = 0.1882
	1391	A_Maj_Brooks = 0.1635 !-- sans dependence au shear
	1392	Dx_Brooks = 200000.
	1393	f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)*(a_Brooks)) (Dx_Brooks**(-b_Brooks))
	1394	ctot_surf(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.))
[3493]	1395	!--------------------------------------------
	1396	!Incloud Condensed water qcloud
	1397	!--------------------------------------------
[5144]	1398	IF (ctot_surf(ind1, ind2) < 1.e-8) THEN
	1399	ctot_vol(ind1, ind2) = 0.
	1400	ctot_surf(ind1, ind2) = 0.
	1401	qcloud(ind1) = zqsatenv(ind1, ind2)
	1402	else
	1403	qcloud(ind1) = qltot(ind1, ind2) / ctot_vol(ind1, ind2) + zqsatenv(ind1, ind2)
	1404	endif
[3493]	1405
	1406	END IF ! From the separation (thermal/envrionnement) et (environnement only)
	1407
	1408	! Outputs used to check the PDFs
[5144]	1409	cloudth_senv(ind1, ind2) = senv
	1410	cloudth_sth(ind1, ind2) = sth
	1411	cloudth_sigmaenv(ind1, ind2) = sigma_env
	1412	cloudth_sigmath(ind1, ind2) = sigma_th
[3493]	1413
[5144]	1414	END DO ! From the loop on ngrid
[3493]	1415
[5144]	1416	END SUBROUTINE cloudth_v6
[5105]	1417
[3999]	1418
[5144]	1419	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1420	SUBROUTINE cloudth_mpc(klon, klev, ind2, mpc_bl_points, &
	1421	& temp, qt, qt_th, frac_th, zpspsk, paprsup, paprsdn, play, thetal_th, &
	1422	& ratqs, qcloud, qincloud, icefrac, ctot, ctot_vol, &
	1423	& cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
	1424	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1425	! Author : Etienne Vignon (LMDZ/CNRS)
	1426	! Date: April 2024
	1427	! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
	1428	! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
	1429	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[3999]	1430
[5144]	1431	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs
	1432	USE lmdz_cloudth_ini, ONLY: C_mpc, Ni, C_cap, Ei, d_top, vert_alpha, vert_alpha_th, sigma1s_factor, sigma1s_power, sigma2s_factor, sigma2s_power, cloudth_ratqsmin, iflag_cloudth_vert_noratqs
	1433	USE lmdz_lscp_tools, only: CALC_QSAT_ECMWF
	1434	USE lmdz_lscp_ini, only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
[3999]	1435
[5144]	1436	IMPLICIT NONE
[3999]	1437
	1438
[5144]	1439	!------------------------------------------------------------------------------
	1440	! Declaration
	1441	!------------------------------------------------------------------------------
[3999]	1442
[5144]	1443	! INPUT/OUTPUT
[3999]	1444
[5144]	1445	INTEGER, INTENT(IN) :: klon, klev, ind2
[3999]	1446
[5144]	1447	REAL, DIMENSION(klon), INTENT(IN) :: temp ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
	1448	REAL, DIMENSION(klon), INTENT(IN) :: qt ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
	1449	REAL, DIMENSION(klon), INTENT(IN) :: qt_th ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
	1450	REAL, DIMENSION(klon), INTENT(IN) :: thetal_th ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
	1451	REAL, DIMENSION(klon), INTENT(IN) :: frac_th ! Fraction of the mesh covered by thermals [0-1]
	1452	REAL, DIMENSION(klon), INTENT(IN) :: zpspsk ! Exner potential
	1453	REAL, DIMENSION(klon), INTENT(IN) :: paprsup ! Pressure at top layer interface [Pa]
	1454	REAL, DIMENSION(klon), INTENT(IN) :: paprsdn ! Pressure at bottom layer interface [Pa]
	1455	REAL, DIMENSION(klon), INTENT(IN) :: play ! Pressure of layers [Pa]
	1456	REAL, DIMENSION(klon), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib.
[3999]	1457
[5144]	1458	INTEGER, DIMENSION(klon, klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds
[3999]	1459
[5144]	1460	REAL, DIMENSION(klon), INTENT(OUT) :: ctot ! Cloud fraction [0-1]
	1461	REAL, DIMENSION(klon), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1]
	1462	REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg]
	1463	REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg]
	1464	REAL, DIMENSION(klon), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1]
	1465	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
	1466	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
	1467	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
	1468	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
[3999]	1469
	1470
[5144]	1471	! LOCAL VARIABLES
[3999]	1472
[5144]	1473	INTEGER itap, ind1, l, ig, iter, k
	1474	INTEGER iflag_topthermals, niter
[4380]	1475
[5144]	1476	REAL qcth(klon)
	1477	REAL qcenv(klon)
	1478	REAL qctot(klon)
	1479	REAL cth(klon)
	1480	REAL cenv(klon)
	1481	REAL cth_vol(klon)
	1482	REAL cenv_vol(klon)
	1483	REAL qt_env(klon), thetal_env(klon)
	1484	REAL icefracenv(klon), icefracth(klon)
	1485	REAL sqrtpi, sqrt2, sqrt2pi
	1486	REAL alth, alenv, ath, aenv
	1487	REAL sth, senv, sigma1s, sigma2s, sigma1s_fraca, sigma1s_ratqs
	1488	REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks
	1489	REAL xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2
	1490	REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv
	1491	REAL IntJ, IntI1, IntI2, IntI3, IntJ_CF, IntI1_CF, IntI3_CF, coeffqlenv, coeffqlth
	1492	REAL zdelta, qsatbef, zcor
	1493	REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
	1494	REAL qlbef
	1495	REAL dqsatenv(klon), dqsatth(klon)
	1496	REAL erf
	1497	REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
	1498	REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)
	1499	REAL rhodz(klon)
	1500	REAL rho(klon)
	1501	REAL dz(klon)
	1502	REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)
	1503	REAL alenvl, aenvl
	1504	REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
	1505	REAL senvi, senvl, qbase, sbase, qliqth, qiceth
	1506	REAL qmax, ttarget, stmp, cout, coutref, fraci
	1507	REAL maxi, mini, pas
[3999]	1508
[5144]	1509	! Modifty the saturation deficit PDF in thermals
	1510	! in the presence of ice crystals
	1511	CHARACTER (len = 80) :: abort_message
	1512	CHARACTER (len = 20) :: routname = 'cloudth_mpc'
[3999]	1513
	1514
[5144]	1515	!------------------------------------------------------------------------------
	1516	! Initialisation
	1517	!------------------------------------------------------------------------------
[3999]	1518
	1519
[5144]	1520	! Few initial checksS
[3999]	1521
[5144]	1522	DO k = 1, klev
[3999]	1523	DO ind1 = 1, klon
[5144]	1524	rhodz(ind1) = (paprsdn(ind1) - paprsup(ind1)) / rg !kg/m2
	1525	zrho(ind1) = play(ind1) / temp(ind1) / rd !kg/m3
	1526	dz(ind1) = rhodz(ind1) / zrho(ind1) !m : epaisseur de la couche en metre
[3999]	1527	END DO
[5144]	1528	END DO
[3999]	1529
[5144]	1530	icefracth(:) = 0.
	1531	icefracenv(:) = 0.
	1532	sqrt2pi = sqrt(2. * rpi)
	1533	sqrt2 = sqrt(2.)
	1534	sqrtpi = sqrt(rpi)
	1535	icefrac(:) = 0.
[3999]	1536
[5144]	1537	!-------------------------------------------------------------------------------
	1538	! Identify grid points with potential mixed-phase conditions
	1539	!-------------------------------------------------------------------------------
[3999]	1540
[5144]	1541	DO ind1 = 1, klon
	1542	IF ((temp(ind1) < RTT) .AND. (temp(ind1) > temp_nowater) &
	1543	.AND. (ind2<=klev - 2) &
	1544	.AND. (frac_th(ind1)>min_frac_th_cld)) THEN
	1545	mpc_bl_points(ind1, ind2) = 1
	1546	ELSE
	1547	mpc_bl_points(ind1, ind2) = 0
	1548	ENDIF
	1549	ENDDO
[3999]	1550
	1551
[5144]	1552	!-------------------------------------------------------------------------------
	1553	! Thermal fraction calculation and standard deviation of the distribution
	1554	!-------------------------------------------------------------------------------
[3999]	1555
[5144]	1556	! initialisations and calculation of temperature, humidity and saturation specific humidity
[3999]	1557
[5144]	1558	DO ind1 = 1, klon
[3999]	1559
[5144]	1560	rhodz(ind1) = (paprsdn(ind1) - paprsup(ind1)) / rg ! kg/m2
	1561	rho(ind1) = play(ind1) / temp(ind1) / rd ! kg/m3
	1562	dz(ind1) = rhodz(ind1) / rho(ind1) ! m : epaisseur de la couche en metre
	1563	Tbefenv(ind1) = temp(ind1)
	1564	thetal_env(ind1) = Tbefenv(ind1) / zpspsk(ind1)
	1565	Tbefth(ind1) = thetal_th(ind1) * zpspsk(ind1)
	1566	rhoth(ind1) = play(ind1) / Tbefth(ind1) / RD
	1567	qt_env(ind1) = (qt(ind1) - frac_th(ind1) * qt_th(ind1)) / (1. - frac_th(ind1)) !qt = aqtth + (1-a)qtenv
[3999]	1568
[5144]	1569	ENDDO
[3999]	1570
[5144]	1571	! Calculation of saturation specific humidity
	1572	CALL CALC_QSAT_ECMWF(klon, Tbefenv, qt_env, play, RTT, 1, .FALSE., qslenv, dqsatenv)
	1573	CALL CALC_QSAT_ECMWF(klon, Tbefenv, qt_env, play, RTT, 2, .FALSE., qsienv, dqsatenv)
	1574	CALL CALC_QSAT_ECMWF(klon, Tbefth, qt_th, play, RTT, 1, .FALSE., qslth, dqsatth)
[4072]	1575
[5144]	1576	DO ind1 = 1, klon
[4072]	1577
[5144]	1578	IF (frac_th(ind1)>min_frac_th_cld) THEN !Thermal and environnement
[4072]	1579
[5144]	1580	! unlike in the other cloudth routine,
	1581	! We consider distributions of the saturation deficit WRT liquid
	1582	! at positive AND negative celcius temperature
	1583	! subsequently, cloud fraction corresponds to the part of the pdf corresponding
	1584	! to superstauration with respect to liquid WHATEVER the temperature
[4072]	1585
[5144]	1586	! Environment:
[4072]	1587
[5144]	1588	alenv = (0.622 * RLVTT * qslenv(ind1)) / (rd * thetal_env(ind1)**2)
	1589	aenv = 1. / (1. + (alenv * RLVTT / rcpd))
	1590	senv = aenv * (qt_env(ind1) - qslenv(ind1))
[4072]	1591
	1592
[5144]	1593	! Thermals:
[3999]	1594
[5144]	1595	alth = (0.622 * RLVTT * qslth(ind1)) / (rd * thetal_th(ind1)**2)
	1596	ath = 1. / (1. + (alth * RLVTT / rcpd))
	1597	sth = ath * (qt_th(ind1) - qslth(ind1))
[3999]	1598
	1599
[5144]	1600	! IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
[3999]	1601
	1602
[5144]	1603	! Standard deviation of the distributions
[3999]	1604
[5144]	1605	sigma1s_fraca = (sigma1s_factor*0.5) (frac_th(ind1)**sigma1s_power) / &
	1606	(1 - frac_th(ind1)) * ((sth - senv)2)0.5
[3999]	1607
[5144]	1608	IF (cloudth_ratqsmin>0.) THEN
	1609	sigma1s_ratqs = cloudth_ratqsmin * qt(ind1)
	1610	ELSE
	1611	sigma1s_ratqs = ratqs(ind1) * qt(ind1)
	1612	ENDIF
[3999]	1613
[5144]	1614	sigma1s = sigma1s_fraca + sigma1s_ratqs
	1615	IF (iflag_ratqs==11) THEN
	1616	sigma1s = ratqs(ind1) * qt(ind1) * aenv
	1617	ENDIF
	1618	sigma2s = (sigma2s_factor * (((sth - senv)2)0.5) / ((frac_th(ind1) + 0.02)*sigma2s_power)) + 0.002 qt_th(ind1)
[3999]	1619
[5144]	1620	deltasenv = aenv * vert_alpha * sigma1s
	1621	deltasth = ath * vert_alpha_th * sigma2s
[3999]	1622
[5144]	1623	xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s)
	1624	xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s)
	1625	exp_xenv1 = exp(-1. * xenv1**2)
	1626	exp_xenv2 = exp(-1. * xenv2**2)
	1627	xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s)
	1628	xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s)
	1629	exp_xth1 = exp(-1. * xth1**2)
	1630	exp_xth2 = exp(-1. * xth2**2)
[3999]	1631
[5144]	1632	!surface CF
[3999]	1633
[5144]	1634	cth(ind1) = 0.5 * (1. - 1. * erf(xth1))
	1635	cenv(ind1) = 0.5 * (1. - 1. * erf(xenv1))
	1636	ctot(ind1) = frac_th(ind1) * cth(ind1) + (1. - 1. * frac_th(ind1)) * cenv(ind1)
[3999]	1637
	1638
[5144]	1639	!volume CF, condensed water and ice fraction
[3999]	1640
[5144]	1641	!environnement
[3999]	1642
[5144]	1643	IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp_xenv2
	1644	IntJ_CF = 0.5 * (1. - 1. * erf(xenv2))
[3999]	1645
[5144]	1646	IF (deltasenv < 1.e-10) THEN
	1647	qcenv(ind1) = IntJ
	1648	cenv_vol(ind1) = IntJ_CF
	1649	ELSE
	1650	IntI1 = (((senv + deltasenv)2 + (sigma1s)2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1))
	1651	IntI2 = (sigma1s*2 / (4 deltasenv * sqrtpi)) * (xenv1 * exp_xenv1 - xenv2 * exp_xenv2)
	1652	IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp_xenv1 - exp_xenv2)
	1653	IntI1_CF = ((senv + deltasenv) * (erf(xenv2) - erf(xenv1))) / (4 * deltasenv)
	1654	IntI3_CF = (sqrt2 * sigma1s * (exp_xenv1 - exp_xenv2)) / (4 * sqrtpi * deltasenv)
	1655	qcenv(ind1) = IntJ + IntI1 + IntI2 + IntI3
	1656	cenv_vol(ind1) = IntJ_CF + IntI1_CF + IntI3_CF
	1657	IF (Tbefenv(ind1) < temp_nowater) THEN
	1658	! freeze all droplets in cirrus temperature regime
	1659	icefracenv(ind1) = 1.
	1660	ENDIF
	1661	ENDIF
[3999]	1662
	1663
	1664
[5144]	1665	!thermals
[3999]	1666
[5144]	1667	IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp_xth2
	1668	IntJ_CF = 0.5 * (1. - 1. * erf(xth2))
[3999]	1669
[5144]	1670	IF (deltasth < 1.e-10) THEN
	1671	qcth(ind1) = IntJ
	1672	cth_vol(ind1) = IntJ_CF
	1673	ELSE
	1674	IntI1 = (((sth + deltasth)2 + (sigma2s)2) / (8 * deltasth)) * (erf(xth2) - erf(xth1))
	1675	IntI2 = (sigma2s*2 / (4 deltasth * sqrtpi)) * (xth1 * exp_xth1 - xth2 * exp_xth2)
	1676	IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp_xth1 - exp_xth2)
	1677	IntI1_CF = ((sth + deltasth) * (erf(xth2) - erf(xth1))) / (4 * deltasth)
	1678	IntI3_CF = (sqrt2 * sigma2s * (exp_xth1 - exp_xth2)) / (4 * sqrtpi * deltasth)
	1679	qcth(ind1) = IntJ + IntI1 + IntI2 + IntI3
	1680	cth_vol(ind1) = IntJ_CF + IntI1_CF + IntI3_CF
	1681	IF (Tbefth(ind1) < temp_nowater) THEN
	1682	! freeze all droplets in cirrus temperature regime
	1683	icefracth(ind1) = 1.
	1684	ENDIF
[3999]	1685
[5144]	1686	ENDIF
[4910]	1687
[5144]	1688	qctot(ind1) = frac_th(ind1) * qcth(ind1) + (1. - 1. * frac_th(ind1)) * qcenv(ind1)
	1689	ctot_vol(ind1) = frac_th(ind1) * cth_vol(ind1) + (1. - 1. * frac_th(ind1)) * cenv_vol(ind1)
[3999]	1690
[5144]	1691	IF (cenv(ind1)<min_neb_th.AND.cth(ind1)<min_neb_th) THEN
	1692	ctot(ind1) = 0.
	1693	ctot_vol(ind1) = 0.
	1694	qcloud(ind1) = qslenv(ind1)
	1695	qincloud(ind1) = 0.
	1696	icefrac(ind1) = 0.
	1697	ELSE
	1698	qcloud(ind1) = qctot(ind1) / ctot(ind1) + qslenv(ind1)
	1699	qincloud(ind1) = qctot(ind1) / ctot(ind1)
	1700	IF (qctot(ind1) > 0) THEN
	1701	icefrac(ind1) = (frac_th(ind1) * qcth(ind1) * icefracth(ind1) + (1. - 1. * frac_th(ind1)) * qcenv(ind1) * icefracenv(ind1)) / qctot(ind1)
	1702	icefrac(ind1) = max(min(1., icefrac(ind1)), 0.)
	1703	ENDIF
	1704	ENDIF
[3999]	1705
	1706
[5144]	1707	! ELSE ! mpc_bl_points>0
[3999]	1708
[5144]	1709	ELSE ! gaussian for environment only
[3999]	1710
[5144]	1711	alenv = (0.622 * RLVTT * qslenv(ind1)) / (rd * thetal_env(ind1)**2)
	1712	aenv = 1. / (1. + (alenv * RLVTT / rcpd))
	1713	senv = aenv * (qt_env(ind1) - qslenv(ind1))
	1714	sth = 0.
	1715	sigma1s = ratqs(ind1) * qt_env(ind1)
	1716	sigma2s = 0.
[3999]	1717
[5144]	1718	xenv = senv / (sqrt(2.) * sigma1s)
	1719	cenv(ind1) = 0.5 * (1. - 1. * erf(xenv))
	1720	ctot(ind1) = 0.5 * (1. + 1. * erf(xenv))
	1721	ctot_vol(ind1) = ctot(ind1)
	1722	qctot(ind1) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1))
[3999]	1723
[5082]	1724	IF (ctot(ind1)<min_neb_th) THEN
[5144]	1725	ctot(ind1) = 0.
	1726	qcloud(ind1) = qslenv(ind1)
	1727	qincloud(ind1) = 0.
	1728	ELSE
	1729	qcloud(ind1) = qctot(ind1) / ctot(ind1) + qslenv(ind1)
	1730	qincloud(ind1) = MAX(qctot(ind1) / ctot(ind1), 0.)
[3999]	1731	ENDIF
	1732
[5144]	1733	ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
[3999]	1734
[5144]	1735	! Outputs used to check the PDFs
	1736	cloudth_senv(ind1) = senv
	1737	cloudth_sth(ind1) = sth
	1738	cloudth_sigmaenv(ind1) = sigma1s
	1739	cloudth_sigmath(ind1) = sigma2s
[3999]	1740
[4114]	1741	ENDDO !loop on klon
[3999]	1742
[5144]	1743	END SUBROUTINE cloudth_mpc
[4114]	1744
[3999]	1745
[5144]	1746	! ELSE ! mpc_bl_points>0
[3999]	1747
[5144]	1748	! ! Treat boundary layer mixed phase clouds
[5105]	1749
[5144]	1750	! ! thermals
	1751	! !=========
[3999]	1752
[5144]	1753	! ! ice phase
	1754	! !...........
[4910]	1755
[5144]	1756	! qiceth=0.
	1757	! deltazlev_mpc=dz(ind1,:)
	1758	! temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
	1759	! pres_mpc=pplay(ind1,:)
	1760	! fraca_mpc=fraca(ind1,:)
	1761	! snowf_mpc=snowflux(ind1,:)
	1762	! iflag_topthermals=0
	1763	! IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN
	1764	! iflag_topthermals = 1
	1765	! ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
	1766	! .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
	1767	! iflag_topthermals = 2
	1768	! ELSE
	1769	! iflag_topthermals = 0
	1770	! ENDIF
[4910]	1771
[5144]	1772	! CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
	1773	! qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
[5099]	1774
[5144]	1775	! ! qmax calculation
	1776	! sigma2s=(sigma2s_factor((MAX((sthl-senvl),0.)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
	1777	! deltasth=athlvert_alpha_thsigma2s
	1778	! xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
	1779	! xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
	1780	! exp_xth1 = exp(-1.xth1*2)
	1781	! exp_xth2 = exp(-1.xth2*2)
	1782	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1783	! IntJ_CF=0.5(1.-1.erf(xth2))
	1784	! IntI1=(((sthl+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1785	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1786	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1787	! IntI1_CF=((sthl+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1788	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
	1789	! qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
[5099]	1790
	1791
[5144]	1792	! ! Liquid phase
	1793	! !................
	1794	! ! We account for the effect of ice crystals in thermals on sthl
	1795	! ! and on the width of the distribution
[5099]	1796
[5144]	1797	! sthlc=sthl1./(1.+C_mpcqith(ind1,ind2)) &
	1798	! + (1.-1./(1.+C_mpcqith(ind1,ind2))) athl*(qsith(ind1,ind2)-qslth(ind1))
[5099]	1799
[5144]	1800	! sigma2sc=(sigma2s_factor((MAX((sthlc-senvl),0.)2)*0.5) &
	1801	! /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
	1802	! +0.002*zqta(ind1,ind2)
	1803	! deltasthc=athlvert_alpha_thsigma2sc
[5099]	1804
	1805
[5144]	1806	! xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
	1807	! xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)
	1808	! exp_xth1 = exp(-1.xth1*2)
	1809	! exp_xth2 = exp(-1.xth2*2)
	1810	! IntJ=0.5sthlc(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
	1811	! IntJ_CF=0.5(1.-1.erf(xth2))
	1812	! IntI1=(((sthlc+deltasthc)2+(sigma2sc)2)/(8deltasthc))(erf(xth2)-erf(xth1))
	1813	! IntI2=(sigma2sc*2/(4deltasthcsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1814	! IntI3=((sqrt2sigma2sc(sthlc+deltasthc))/(4sqrtpideltasthc))*(exp_xth1-exp_xth2)
	1815	! IntI1_CF=((sthlc+deltasthc)(erf(xth2)-erf(xth1)))/(4deltasthc)
	1816	! IntI3_CF=(sqrt2sigma2sc(exp_xth1-exp_xth2))/(4sqrtpideltasthc)
	1817	! qliqth=IntJ+IntI1+IntI2+IntI3
[5099]	1818
[5144]	1819	! qlth(ind1,ind2)=MAX(0.,qliqth)
[5099]	1820
[5144]	1821	! ! Condensed water
[5099]	1822
[5144]	1823	! qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
[5099]	1824
	1825
[5144]	1826	! ! consistency with subgrid distribution
[5099]	1827
[5144]	1828	! IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
	1829	! fraci=qith(ind1,ind2)/qcth(ind1,ind2)
	1830	! qcth(ind1,ind2)=qmax
	1831	! qith(ind1,ind2)=fraci*qmax
	1832	! qlth(ind1,ind2)=(1.-fraci)*qmax
	1833	! ENDIF
[5099]	1834
[5144]	1835	! ! Cloud Fraction
	1836	! !...............
	1837	! ! calculation of qbase which is the value of the water vapor within mixed phase clouds
	1838	! ! such that the total water in cloud = qbase+qliqth+qiceth
	1839	! ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
	1840	! ! sbase and qbase calculation (note that sbase is wrt liq so negative)
	1841	! ! look for an approximate solution with iteration
[5099]	1842
[5144]	1843	! ttarget=qcth(ind1,ind2)
	1844	! mini= athl*(qsith(ind1,ind2)-qslth(ind1))
	1845	! maxi= 0. !athl(3.sqrt(sigma2s))
	1846	! niter=20
	1847	! pas=(maxi-mini)/niter
	1848	! stmp=mini
	1849	! sbase=stmp
	1850	! coutref=1.E6
	1851	! DO iter=1,niter
	1852	! cout=ABS(sigma2s/SQRT(2.RPI)EXP(-((sthl-stmp)/sigma2s)*2)+(sthl-stmp)/SQRT(2.)(1.-erf(-(sthl-stmp)/sigma2s)) &
	1853	! + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
	1854	! IF (cout .LT. coutref) THEN
	1855	! sbase=stmp
	1856	! coutref=cout
	1857	! ELSE
	1858	! stmp=stmp+pas
	1859	! ENDIF
	1860	! ENDDO
	1861	! qbase=MAX(0., sbase/athl+qslth(ind1))
[5099]	1862
[5144]	1863	! ! surface cloud fraction in thermals
	1864	! cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
	1865	! cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
[5099]	1866
	1867
[5144]	1868	! !volume cloud fraction in thermals
	1869	! !to be checked
	1870	! xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
	1871	! xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)
	1872	! exp_xth1 = exp(-1.xth1*2)
	1873	! exp_xth2 = exp(-1.xth2*2)
[5099]	1874
[5144]	1875	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1876	! IntJ_CF=0.5(1.-1.erf(xth2))
[5099]	1877
[5144]	1878	! IF (deltasth .LT. 1.e-10) THEN
	1879	! cth_vol(ind1,ind2)=IntJ_CF
	1880	! ELSE
	1881	! IntI1=(((sthl+deltasth-sbase)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1882	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1883	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1884	! IntI1_CF=((sthl-sbase+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1885	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
	1886	! cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
	1887	! ENDIF
	1888	! cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
[5099]	1889
	1890
	1891
[5144]	1892	! ! Environment
	1893	! !=============
	1894	! ! In the environment/downdrafts, ONLY liquid clouds
	1895	! ! See Shupe et al. 2008, JAS
[5099]	1896
[5144]	1897	! ! standard deviation of the distribution in the environment
	1898	! sth=sthl
	1899	! senv=senvl
	1900	! sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
	1901	! & (1-fraca(ind1,ind2))(MAX((sth-senv),0.)2)*0.5
	1902	! ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
	1903	! ! in the environement
[5099]	1904
[5144]	1905	! sigma1s_ratqs=1E-10
	1906	! IF (cloudth_ratqsmin>0.) THEN
	1907	! sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
	1908	! ENDIF
[5099]	1909
[5144]	1910	! sigma1s = sigma1s_fraca + sigma1s_ratqs
	1911	! IF (iflag_ratqs.EQ.11) THEN
	1912	! sigma1s = ratqs(ind1,ind2)po(ind1)aenv
	1913	! ENDIF
	1914	! IF (iflag_ratqs.EQ.11) THEN
	1915	! sigma1s = ratqs(ind1,ind2)po(ind1)aenvl
	1916	! ENDIF
	1917	! deltasenv=aenvlvert_alphasigma1s
	1918	! xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
	1919	! xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
	1920	! exp_xenv1 = exp(-1.xenv1*2)
	1921	! exp_xenv2 = exp(-1.xenv2*2)
[5099]	1922
[5144]	1923	! !surface CF
	1924	! cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
[5099]	1925
[5144]	1926	! !volume CF and condensed water
	1927	! IntJ=0.5senvl(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
	1928	! IntJ_CF=0.5(1.-1.erf(xenv2))
[5099]	1929
[5144]	1930	! IF (deltasenv .LT. 1.e-10) THEN
	1931	! qcenv(ind1,ind2)=IntJ
	1932	! cenv_vol(ind1,ind2)=IntJ_CF
	1933	! ELSE
	1934	! IntI1=(((senvl+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	1935	! IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
	1936	! IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
	1937	! IntI1_CF=((senvl+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
	1938	! IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
	1939	! qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
	1940	! cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
	1941	! ENDIF
[5099]	1942
[5144]	1943	! qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
	1944	! cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
[5099]	1945
	1946
	1947
[5144]	1948	! ! Thermals + environment
	1949	! ! =======================
	1950	! ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	1951	! qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
	1952	! ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
	1953	! IF (qcth(ind1,ind2) .GT. 0) THEN
	1954	! icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
	1955	! /(fraca(ind1,ind2)*qcth(ind1,ind2) &
	1956	! +(1.-1.fraca(ind1,ind2))qcenv(ind1,ind2))
	1957	! icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
	1958	! ELSE
	1959	! icefrac(ind1,ind2)=0.
	1960	! ENDIF
[5099]	1961
[5144]	1962	! IF (cenv(ind1,ind2).LT.1.e-10.OR.cth(ind1,ind2).LT.1.e-10) THEN
	1963	! ctot(ind1,ind2)=0.
	1964	! ctot_vol(ind1,ind2)=0.
	1965	! qincloud(ind1)=0.
	1966	! qcloud(ind1)=zqsatenv(ind1)
	1967	! ELSE
	1968	! qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
	1969	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
	1970	! qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
	1971	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
	1972	! ENDIF
[5099]	1973
[5144]	1974	! ENDIF ! mpc_bl_points
[5099]	1975
	1976
	1977
[5144]	1978	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1979	SUBROUTINE ICE_MPC_BL_CLOUDS(ind1, ind2, klev, Ni, Ei, C_cap, d_top, iflag_topthermals, temp, pres, qth, qsith, qlth, deltazlev, vith, fraca, qith)
[5099]	1980
[5144]	1981	! parameterization of ice for boundary
	1982	! layer mixed-phase clouds assuming a stationary system
[5099]	1983
[5144]	1984	! Note that vapor deposition on ice crystals and riming of liquid droplets
	1985	! depend on the ice number concentration Ni
	1986	! One could assume that Ni depends on qi, e.g., Ni=beta(qirho)**xi
	1987	! and use values from Hong et al. 2004, MWR for instance
	1988	! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
	1989	! One could also think of a more complex expression of Ni;
	1990	! function of qi, T, the concentration in aerosols or INP ..
	1991	! Here we prefer fixing Ni to a tuning parameter
	1992	! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
	1993	! in Mioche et al. 2017
[4910]	1994
	1995
[5144]	1996	! References:
	1997	!------------
	1998	! This parameterization is thoroughly described in Vignon et al.
[4910]	1999
[5144]	2000	! More specifically
	2001	! for the Water vapor deposition process:
[3999]	2002
[5144]	2003	! Rotstayn, L. D., 1997: A physically based scheme for the treat-
	2004	! ment of stratiform cloudfs and precipitation in large-scale
	2005	! models. I: Description and evaluation of the microphysical
	2006	! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
[5099]	2007
[5144]	2008	! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
	2009	! stratiform cloud microphysics scheme in the NCAR Com-
	2010	! munity Atmosphere Model (CAM3). Part I: Description and
	2011	! numerical tests. J. Climate, 21, 3642–3659
[5099]	2012
[5144]	2013	! for the Riming process:
[5099]	2014
[5144]	2015	! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
	2016	! scale structure and organization of clouds and precipitation in
	2017	! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
	2018	! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
[5099]	2019
[5144]	2020	! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
	2021	! forecasts of winter precipitation using an improved bulk
	2022	! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
	2023	! forecasts of winter precipitation using an improved bulk
	2024	! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
[5099]	2025
[5144]	2026	! For the formation of clouds by thermals:
[5099]	2027
[5144]	2028	! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
	2029	! the Atmospheric Sciences, 65, 407–425.
[5099]	2030
[5144]	2031	! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
	2032	! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
[5099]	2033
	2034
	2035
[5144]	2036	! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
	2037	!=============================================================================
[5099]	2038
[5117]	2039	USE phys_state_var_mod, ONLY: fm_therm, detr_therm, entr_therm
[5144]	2040	USE lmdz_yomcst
[3999]	2041
[5134]	2042	IMPLICIT NONE
[3999]	2043
[5144]	2044	INTEGER, INTENT(IN) :: ind1, ind2, klev ! horizontal and vertical indices and dimensions
[4072]	2045	INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ?
[5144]	2046	REAL, INTENT(IN) :: Ni ! ice number concentration [m-3]
	2047	REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency
	2048	REAL, INTENT(IN) :: C_cap ! ice crystal capacitance
	2049	REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter
	2050	REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals
	2051	REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa]
	2052	REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg]
	2053	REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg]
	2054	REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg]
	2055	REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m]
	2056	REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s]
	2057	REAL, DIMENSION(klev + 1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals
	2058	REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg]
[3999]	2059
[5144]	2060	INTEGER ind2p1, ind2p2
[3999]	2061	REAL rho(klev)
	2062	REAL unsurtaudet, unsurtaustardep, unsurtaurim
[4072]	2063	REAL qi, AA, BB, Ka, Dv, rhoi
	2064	REAL p0, t0, fp1, fp2
[3999]	2065	REAL alpha, flux_term
	2066	REAL det_term, precip_term, rim_term, dep_term
	2067
[5144]	2068	ind2p1 = ind2 + 1
	2069	ind2p2 = ind2 + 2
[3999]	2070
[5144]	2071	rho = pres / temp / RD ! air density kg/m3
[3999]	2072
[5144]	2073	Ka = 2.4e-2 ! thermal conductivity of the air, SI
	2074	p0 = 101325.0 ! ref pressure
	2075	T0 = 273.15 ! ref temp
	2076	rhoi = 500.0 ! cloud ice density following Reisner et al. 1998
	2077	alpha = 700. ! fallvelocity param
[3999]	2078
[5144]	2079	IF (iflag_topthermals > 0) THEN ! uppermost thermals levels
[3999]	2080
[5144]	2081	Dv = 0.0001 * 0.211 * (p0 / pres(ind2)) * ((temp(ind2) / T0)**1.94) ! water vapor diffusivity in air, SI
[3999]	2082
[5144]	2083	! Detrainment term:
[3999]	2084
[5144]	2085	unsurtaudet = detr_therm(ind1, ind2) / rho(ind2) / deltazlev(ind2)
[3999]	2086
[4072]	2087
[5144]	2088	! Deposition term
	2089	AA = RLSTT / Ka / temp(ind2) * (RLSTT / RV / temp(ind2) - 1.)
	2090	BB = 1. / (rho(ind2) * Dv * qsith(ind2))
	2091	unsurtaustardep = C_cap * (Ni*0.66) (qth(ind2) - qsith(ind2)) / qsith(ind2) &
	2092	* 4. * RPI / (AA + BB) * (6. * rho(ind2) / rhoi / RPI / Gamma(4.))**(0.33)
[3999]	2093
[5144]	2094	! Riming term neglected at this level
	2095	!unsurtaurim=rho(ind2)alpha3./rhoi/2.Eiqlth(ind2)((p0/pres(ind2))*0.4)
[3999]	2096
[5144]	2097	qi = fraca(ind2) * unsurtaustardep / MAX((d_top * unsurtaudet), 1E-12)
	2098	qi = MAX(qi, 0.)**(3. / 2.)
[3999]	2099
	2100	ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
	2101
[5144]	2102	Dv = 0.0001 * 0.211 * (p0 / pres(ind2p1)) * ((temp(ind2p1) / T0)**1.94) ! water vapor diffusivity in air, SI
[3999]	2103
[5144]	2104	! Detrainment term:
[3999]	2105
[5144]	2106	unsurtaudet = detr_therm(ind1, ind2p1) / rho(ind2p1) / deltazlev(ind2p1)
	2107	det_term = -unsurtaudet * qith(ind2p1) * rho(ind2p1)
[3999]	2108
	2109
[5144]	2110	! Deposition term
[3999]	2111
[5144]	2112	AA = RLSTT / Ka / temp(ind2p1) * (RLSTT / RV / temp(ind2p1) - 1.)
	2113	BB = 1. / (rho(ind2p1) * Dv * qsith(ind2p1))
	2114	unsurtaustardep = C_cap * (Ni*0.66) (qth(ind2p1) - qsith(ind2p1)) &
	2115	/ qsith(ind2p1) * 4. * RPI / (AA + BB) &
	2116	* (6. * rho(ind2p1) / rhoi / RPI / Gamma(4.))**(0.33)
	2117	dep_term = rho(ind2p1) * fraca(ind2p1) * (qith(ind2p1)*0.33) unsurtaustardep
[3999]	2118
[5144]	2119	! Riming term
[3999]	2120
[5144]	2121	unsurtaurim = rho(ind2p1) * alpha * 3. / rhoi / 2. * Ei * qlth(ind2p1) * ((p0 / pres(ind2p1))**0.4)
	2122	rim_term = rho(ind2p1) * fraca(ind2p1) * qith(ind2p1) * unsurtaurim
[4072]	2123
[5144]	2124	! Precip term
[3999]	2125
[5144]	2126	! We assume that there is no solid precipitation outside thermals
	2127	! so the precipitation flux within thermals is equal to the precipitation flux
	2128	! at mesh-scale divided by thermals fraction
[3999]	2129
[5144]	2130	fp2 = 0.
	2131	fp1 = 0.
	2132	IF (fraca(ind2p1) > 0.) THEN
	2133	fp2 = -qith(ind2p2) * rho(ind2p2) * vith(ind2p2) * fraca(ind2p2)! flux defined positive upward
	2134	fp1 = -qith(ind2p1) * rho(ind2p1) * vith(ind2p1) * fraca(ind2p1)
	2135	ENDIF
[3999]	2136
[5144]	2137	precip_term = -1. / deltazlev(ind2p1) * (fp2 - fp1)
[4072]	2138
[5144]	2139	! Calculation in a top-to-bottom loop
[3999]	2140
[5144]	2141	IF (fm_therm(ind1, ind2p1) > 0.) THEN
	2142	qi = 1. / fm_therm(ind1, ind2p1) * &
	2143	(deltazlev(ind2p1) * (-rim_term - dep_term - det_term - precip_term) + &
	2144	fm_therm(ind1, ind2p2) * (qith(ind2p1)))
	2145	END IF
[5105]	2146
[5144]	2147	ENDIF ! top thermals
[3999]	2148
[5144]	2149	qith(ind2) = MAX(0., qi)
	2150
	2151	END SUBROUTINE ICE_MPC_BL_CLOUDS
	2152
[4651]	2153	END MODULE lmdz_cloudth

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

Download in other formats:

Original Format